Ethics statementThis study was carried out in accordance with the recommendations of the Care and Use of Laboratory Animals of the Chinese Academy of Fishery Sciences (Ethics Committee number: YSFRI2019002). The protocol was approved by the Animal Care and Use Committee of the Chinese Academy of Fishery Sciences.Selective breeding and sample preparationThe selective breeding of the vibriosis resistant and susceptible families for C. semilaevis were performed as previously described3. Firstly, the genetic sex of parental fish was identified by a sex-specific AFLP marker20. Full-sib families were then constructed by strip spawning and were reared in separate common tanks under a flow-through system with the same rearing conditions. To trace the lineage, we tagged each family with visible implant elastomers and the pedigree information was precisely recorded3.In each breeding generation, we performed artificial Vibrio challenge tests using the fish with an average size of 10–12 cm. A medial lethal dose (LD50) of V. harveyi challenge was determined as previously described3. After intraperitoneal injection, we recorded the mortality of each family at 12 days post injection, and the families with a survival rate >80% and <30% were considered as Vibrio resistant (VR) and susceptible (VS) families, respectively. To date, the artificial challenge and selection for the resistant families processed continuously for more than 12 years, lasting for 5 generations. Fish from the VR and VS families were used for RNA collection (Fig. 1).Fig. 1The experimental workflow for resistant tongue sole selection and RNA-Seq data analysis workflows.RNA-Seq and comparative transcriptomic analysesTo characterize and compare the gene expression patterns, we collected the liver, spleen and intestine tissues from the resistant and susceptible fish, respectively. Total RNA was extracted from three biological replicates of each group with TRIzol (Invitrogen, USA). Pair-ended (PE) RNA-seq libraries were constructed using the Truseq mRNA stranded RNA-Seq Library Prep Kit (Illumina, USA) according to the standard protocol. Sequencing of totally 18 libraries was conducted with a Illumina HiSeq 2000 sequencing platform, generating raw reads with a read length of 2 × 100 bp and an insert size of 350 bp. The raw data was then quality filtered using RNA-QC-Chain21, filtering out the ambiguous N’s, adaptors, low quality reads with more than 20% of the bases having a quality score < 20. Finally, we obtained 62.26–79.69 million raw reads per sample, amounting to a total of 120.35 Gb clean data (Table 1). We calculated the fragments per kilobase per million mapped sequence reads (FPKM) value for each gene using RSEM (v1.2.12)22, and used NOIseq23. to identify the DEGs, with log2|FoldChange| ≥ 1 with a probability ≥0.9.Table 1 Summary of the RNA-Seq data for liver, spleen, and intestine from resistant and susceptible families.The clean reads were mapped to the reference genome of the tongue sole (NCBI Accession No. GCA_000523025.1) using BWA24 with default parameters. The average mapping rate ranged from 83.95% to 93.18%, averagely at 88.66% (Table 1).KEGG and GO enrichment analysisKEGG and GO enrichment analyses were performed using phyper in R software. KEGG pathways and GO terms with corrected p-values < 0.05 were considered as significantly enriched ones. GO terms are clustered into three main categories: biological processes (BP), molecular functions (MF), and cellular components (CC).In liver, the central metabolic organ, the 1303 DEGs, including 649 up and 654 down-expressed gene, were enriched in “metabolic pathways” (q < 0.05), covering metabolism processes such as “arginine biosynthesis”, “pyruvate metabolism” and “sphingolipid metabolism” (Table 2).Table 2 The enriched KEGG pathways of the DEGs soles in liver, spleen and intestine between the resistant and susceptible tongue soles.In spleen, the DEGs were most significantly abundant in “complement and coagulation cascades” (q < 0.05) (Table 2). Interestingly, the 347 up-expressed genes were mainly involved in immune systems, such as “NOD-like receptor signaling pathway” (q < 0.05), “TNF signaling pathway”, “Th1 and Th2 cell differentiation” and “IL-17 signaling pathway”, while the 703 down-expressed genes were significantly enriched in a variety of metabolic pathways. The significantly enriched GO_MF terms for the down-expressed genes were “aromatic amino acid family metabolic process” and “oxidoreductase activity” (q < 0.05).In intestine, we observed 374 over-expressed and 260 down-expressed genes in the resistant fish. The over-expressed genes were mostly abundant in “pancreatic secretion” pathway (q < 0.05) (Table 2) and GO_MF term of “serine-type endopeptidase activity” (q < 0.05). The down-expressed genes were significantly enriched in “fatty acid metabolism” “phagosome” and “apoptosis”, as well as a GO term of “phospholipid catabolic process” (q < 0.05).Weighted gene correlation network analysis (WGCNA)To gain a holistic and comprehensive gene expression landscape, we used the total 30 RNA-seq datasets, including the newly generated 18 datasets and the previously released 12 datasets, to construct a multiple-organ gene co-expression network.First, we filtered and normalized the genes based on (1) the FPKM value was larger than 1; (2) genes had the FPKM values in more than 16 samples; (3) adjusted the artifacts to make every gene follow an approximate normal distribution. As a result, 20512 genes were used in the normalized expression matrix. Using the mean FPKM values of these genes, we constructed a gene co-expression network and identified the gene co-expression modules with the R package ‘WGCNA’25,26. The optimal soft thresholding power of 30 for adjacency computation was graphically determined when the degree of independence was 0.8. Subsequently, the cutreeDynamic function was used for tree pruning of the gene hierarchical clustering dendrograms, resulting in 15 co-expression modules constructed with corresponding Eigengenes (Fig. 2A). Module-sample associations were estimated using the Spearman correlation coefficient analysis between the module eigengenes and tissues type. As a result, we identified 12 modules that were strongly correlated with different tissue types, which was visualized as a heatmap (Fig. 2B). All these tissues might be involved in the determination or regulation of the disease resistance.Fig. 2Multiple-tissue gene co-expression network for vibriosis resistance in the tongue sole. (A) Hierarchical clustering dendrograms of the co-expressed genes in the modules were identified by WGCNA. A total of 15 modules are colored on the left and top side. The heatmap shows the Person’s association matrix among all genes in the analysis. Yellow color represents low associations and progressively darker red color represents higher associations between genes. (B) Associations between co-expressed gene modules and tissue types. Each row corresponds to a module, and each column to a sample. Each cell contains the corresponding correlation with a p-value. The table is color-coded by correlation with intensity and directions of correlations indicated on the right side of the heatmap.Fig. 3Q-PCR validation of the gene expression levels obtained from RNA-seq. VS: vibriosis susceptible fish, VR: vibriosis resistant fish. Bar height indicates the gene expression level measured by q-PCR (left Y-axis, relative expression measured by ΔΔCt method) and line graph represents RNA-seq measurements (right Y-axis, Fragments Per Kilobase Million, FPKM). X-axis shows gene names with tissues in parentheses.Genes in the salmon module had a significant correlation with the gill samples and were significantly enriched in KEGG pathways of “nitrogen metabolism” and “proximal tubule bicarbonate reclamation” (q < 0.05), as well as GO_MF terms of “ammonium transmembrane transporter activity” (q < 0.05).The purple module had a correlation with spleen. Genes in this module were overrepresented in the digestive system, such as “pancreatic secretion”, “protein digestion and absorption” and “fat digestion and absorption” (q < 0.05), as well as GO_MF terms of a series of “peptidase activity” (p < 0.05).The cyan and green modules, which were strongly correlated with the liver tissue, were enriched in the KEGG pathways of “complement and coagulation cascades” and “phagosome” (q < 0.05).The modules in the red and yellow color were correlated with intestine. Genes in these modules were significantly enriched for “glycosphingolipid biosynthesis”, “endocrine, renin-angiotensin system” and “mineral absorption” pathways (q < 0.05).Genes in the brown, pink, tan, black, magenta and the greenyellow modules were correlated with skin. Among them, the brown, pink, tan and black modules were enriched in a number of cellular community pathways, such as “focal adhesion”, “tight junction”, and “adherens junction” (q < 0.05). The greenyellow module was mainly enriched for “cardiac muscle contraction” and signal transduction pathways, such as “calcium signaling pathway” and “cGMP-PKG signaling pathway” (q < 0.05).Overall, our results suggest that resistance to Vibrio involves a multifaceted response encompassing metabolic adaptations, immune signaling, and tissue-specific gene expression patterns. These findings have implications for understanding the molecular basis of resistance and can inform the development of breeding strategies aimed at enhancing disease resistance in aquaculture species.
