Comparative transcriptome profiling of immune response against Vibrio harveyi infection in Chinese tongue sole

Vibrio harveyi is a major bacterial pathogen that causes fatal vibriosis in Chinese tongue sole (Cynoglossus semilaevis), resulting in massive mortality in the farming industry. However, the molecular mechanisms of C. semilaevis response to V. harveyi infection are poorly understood. Here, we performed transcriptomic analysis of C. semilaevis, comparing resistant and susceptible families in response to V. harveyi challenge (CsRC and CsSC) and control conditions (CsRU and CsSU). RNA libraries were constructed using 12 RNA samples isolated from three biological replicates of the four groups. We performed transcriptome sequencing on an Illumina HiSeq platform, and generating a total of 1,095 million paired-end reads, with the number of clean reads per library ranging from 75.27 M to 99.97 M. Through pairwise comparisons among the four groups, we identified 713 genes exhibiting significant differences at the transcript level. Furthermore, the expression levels were validated by real-time qPCR. Our results provide a valuable resource and new insights into the immune response to V. harveyi infection.


Background & Summary
Knowledge of fish immune systems contributes to understanding the evolution of the immune system, and there is an increasing interest in fish immunology for its unique position in the evolutionary spectrum from lower vertebrates to higher vertebrates 1 . Meanwhile, infectious pathogens, such as bacteria, mould, viruses and protozoans, cause a mass mortality in commercial fish, therefore, it is urgent to study the underlying molecular mechanisms of fish immunology, and to explore novel methods to enhance defences against pathogens in fish 2, 3 .
Previous studies on immune analyses in fish have primarily concentrated on several important genes in model species 4,5 , while the response against bacterial infection in other immune-regulated genes is still unclear. Nevertheless, transcriptomic profiling using next-generation sequencing technologies provides a new approach to studying fish immunology in various marine aquatic species. For example, transcriptomic profiling is conducted to evaluate whole-genome expression patterns in the immune response against bacterial and viral infection to analyze any relevant differences observed. In Epinephelus coioides, transcriptome analysis during Vibrio alginolyticus infection revealed changes in immune gene expression with concomitant induction of innate immune-related complement and hepcidin systems 6 . Transcriptomic analysis of Salmo salar harbouring an infectious salmon anemia virus revealed 3,023 differentially expressed transcripts, with extreme differences in the expression of viral segments between susceptible and resistant groups 7 . Furthermore, transcriptomic profiling sheds lights on potentially novel immune-related transcripts. Transcriptome analysis of C. semilaevis responding to Vibrio anguillarum infection identified multiple differentially expressed annotated and novel genes, which were mostly relevant to the immune response, immune system regulation, and cytokine production 8 . Taken together, these transcriptomic analyses of the response to bacterial and viral infection in teleosts allow us to understand the molecular mechanisms of immune response and to identify novel genes associated with bacterial infection.
C. semilaevis is a valuable marine aquatic species distributed in Northern China 9 . However, vibriosis, which is caused by various bacteria such as Vibrio harveyi, Vibrio anguillarum, Vibrio alginolyticus, Vibrio Parahemolyticus, Vibrio rotiferianus, and Vibrio aestuarianus, has severely disrupted the development of C. semilaevis aquaculture. In C. semilaevis farming, V. harveyi is a major pathogen, causing severe infectious vibriosis with symptoms of putrefied skin, ascites, and tail rot. Although some studies examining C. semilaevis with V. harveyi infection have been reported 10,11 , the underlying molecular mechanisms mounted against V. harveyi infection by the host have not been extensively studied, and the exploitation of genetic resources is required. To address this knowledge gap, we selected two C. semilaevis families based on their significant mortality differences after V. harveyi infection. One family with a high mortality rate (cumulative mortality rate, CMR, >80%) was considered the V. harveyi susceptible family, whereas the other one with a low mortality rate (CMR < 20%) was considered the V. harveyi resistant family. Understanding the different immune molecular mechanisms will be very helpful for enhancing host ability against V. harveyi infection and for breeding V. harveyi resistant strains of C. semilaevis.
Herein, we performed the transcriptome analyses of two phenotypes of C. semilaevis (susceptible and resistant to V. harveyi) under V. harveyi challenge and control conditions. We discribe the detailed procedure of our experimental design including the treatment of fish, tissues collection, library construction and transcriptome  www.nature.com/scientificdata www.nature.com/scientificdata/ sequencing. Quality control was conducted to evaluate the quality of our transcriptome data using FastQC, and a high-quality dataset is presented. Additionally, we performed comparative transcriptomic analyses of four C. semilaevis groups with the aim of screening key genes that cause the differences in disease resistance between resistant and susceptible families and providing an improved understanding of the immune response to V. harveyi infection.

Methods
Ethical approval. The collection and handling of the animals in the study was approved by the Animal Care and Use Committee of Chinese Academy of Fishery Sciences' , and all animals and experiments were conducted in accordance with the guidelines for the care and use of laboratory animals at the Chinese Academy of Fishery Sciences.
Fish rearing and bacterial challenge. The fish (109 ± 24.8 g) used in this experiment were obtained from two C. semilaevis families described above at the Haiyang High-Tech Experimental Base (Shandong, China). Fish were kept in seawater ponds with a continuous supply of seawater at a temperature of 20~23 °C. After 7 days' acclimation, the fish were challenged with Vibrio harveyi (kept by Key Laboratory for Sustainable Utilization of  www.nature.com/scientificdata www.nature.com/scientificdata/ Marine Fisheries Resources). A pre-test was conducted to confirm the concentration of V. harveyi (8*10 4 cfu/ ml). Fish were randomly selected from the two families and challenged with the same concentration of V. harveyi by intraperitoneal injection based on their weights (2 ml/kg). Fish were sampled before injection and 24 h after infection, and the liver, spleen, and kidney tissues were collected from three individual fish in each group and immediately frozen in liquid nitrogen. Tissues were stored at −80 °C until RNA extraction. All fish were anesthetized with a lethal dose of MS-222 (300 ppm) to prevent suffering. The unchallenged and challenged resistant families of C. semilaevis were termed the CsRU and CsRC groups, respectively. The unchallenged and challenged susceptible family of C. semilaevis were termed the CsSU and CsSC groups, respectively. Three samples were used in each group (Table 1).
RNA extraction, library construction, RNA sequencing. Total RNA was extracted with TRIzol reagents (Invitrogen, USA) following the instructions of the manufacturer. Purified RNA was quantified using Qubit ® RNA Assay Kit in a Qubit ® 2.0 Fluorimeter (Life Technologies, CA, USA), and its integrity was evaluated using the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Equal amounts of total RNA from the kidney, spleen, and liver of individual fish were pooled to generate the RNA sample preparation as one biological replicate. Three biological replicates of each group were used to construct cDNA libraries following the Illumina standard operating procedure. Libraries were sequenced on an Illumina HiSeq platform to generate 150 bp paired-end reads.
Quality validation, data cleaning and normalization. We used FastQC 12 to assess the quality of raw reads in fastq format, and all results were merged and visualized using MultiQC 13 . Clean reads were generated from raw reads by removing low quality reads and those containing adapters, poly-N using RNA-QC-Chain 14 with default parameters, then mapped onto the C. semilaevis reference genome (Accession no. GCF_000523025.1) using TopHat software with the parameter of mismatch = 2. We then used Cufflinks with default parameters to construct and identify both known and novel transcripts from TopHat alignment results 15 . Subsequently, we used HTSeq. 16 to count the number of fragments mapped to each gene with the parameters: -m union, -s no, and the expected number of fragments per kilobase of transcript sequence per Millions base pairs (FPKM) were calculated to assess the expression levels.

Downstream analysis.
We used the DESeq package to conduct differential expression analysis 17 and the p-values were adjusted by the Benjamini & Hochberg method for controlling the false discovery rate 18 . Genes with an adjusted p-value < 0.05 were considered differentially expressed genes (DEGs). Furthermore, we calculated the Pearson correlation between samples according to gene expression profiles and the correlation matrix was visualized using ggplot2 19 . Box plots, volcano plots, heat maps and Venn diagrams were drawn using R packages. The analysis workflow is shown in Fig. 1.

Real-time qPCR validation.
In this study, we randomly selected 24 genes for real-time qPCR validation to confirm the results of differential expression analysis. Real-time qRCR was performed with SYBR ® Premix Ex Taq ™ (TaKaRa, Japan) on an ABI 7500 Fast Real-Time PCR system (Applied Biosystems, USA) under the following conditions: denaturation at 95 °C for 30 s, then 40 cycles of 95 °C for 15 s, 60 °C for 20 s, and 72 °C for 10 s. Relative gene expression was analyzed by 2 −ΔΔCt method. β-actin was chosen as the internal control for normalization 20 . We used Prism software to determine statistical significance and draw plots.

Data Records
Raw FASTQ files were deposited into the NCBI Sequence Read Archive (SRA) with accession number SRP186770 (Table 1) Table 3. Statistics analysis of clean reads mapping onto reference genome.

Technical Validation
All RNA samples used for library construction had 260:280 ratios of ≥1.5 and an RNA integrity number (RIN) of ≥8. We constructed 12 RNA libraries of mixed tissues with three biological replicates from four groups (CsSU, CsRU, CsSC, and CsRC) (Fig. 1). We applied FastQC and RNA-QC-Chain to verify that the data was suitable for downstream analysis (Fig. 2, Table 2). www.nature.com/scientificdata www.nature.com/scientificdata/ After clean reads were mapped onto the C. semilaevis reference genome, we calculated the number and percentage of uniquely mapped reads and multiply mapped reads in Table 3. The correlation of gene expression levels between samples is an important index to verify the reliability of an experiment, and the square of the Pearson correlation coefficient (R2) of >0.9 was a prerequisite for differential expression analysis (Fig. 3a). The FPKM boxplot shows the distribution of gene expression levels in Fig. 3b. Additionally, we analyzed the expression profiles among the four groups in the pairwise comparisons. As shown in Fig. 3c, downregulated and upregulated DEGs are highlighted in green and red with a threshold of −log 10 (adjusted p-value) ≥1.3, respectively. Furthermore, a cluster analysis of the DEGs indicated that the expression patterns of those groups differed significantly from each other (Fig. 3d). We identified a total of 713 DEGs in four pairwise comparisons (CsRC vs CsRU, CsRC vs CsSC, CsRUvs CsSU and CsSC vs CsSU) (Fig. 3e). Although the values of the log 2 fold change from the transcriptomic analysis and qPCR analysis were different, the differential expression levels of these selected genes by qPCR were highly consistent with those observed by RNA-Seq (Fig. 3f). The primers for these genes are shown in Table 4.
Taken together, our findings present a high-quality transcriptomic dataset characterizing the C. semilaevis response to V. harveyi infection. Additionally, we screened multiple genes associated with the immune response to V. harveyi infection. The dataset provides a valuable resource for isolating the immune-related genes, for better understanding the biological process of disease resistance, and for exploring reliable ways of host immune defence against V. harveyi.

Code availability
The softwares used for data processing are included in the methods and available in the following list: