Overcrowding as a chronic stressor in teleosts may increase the incidence of physical injuries [1-3] and susceptibility to diseases [4,5] and many suppress the immune response [6,7]. Many studies have demonstrated the influence of high stock density on the immune system [8]. Lymphocytopenia is a classic hematological response to stress and the number of lymphocytes can be decreased under chronic stress via the suppression of lymphocyte proliferation activities and the induction of apoptosis by B cells [9,10]. The stock density level has a direct relationship with the decrease in lymphocytes [11,12]. In addition, overcrowding can inhibit the production of specific antibodies IgM in fish blood by elevating cortisol levels [13]. The negative effect of high stock density on immunoglobulin M (IgM) production was observed when the sub-Antarctic notothenioid fish (Eleginops maclovinus ) was challenged with Piscirickettsia salmonis protein extract [13]. Similarly, elevated cortisol levels inhibited IgM production in Plecoglossus altivelis when kept at a high stock density [14]. By increasing the production of relative oxygen species, a high stock density can also affect the level of oxidative stress and alter cellular functions [15]. However, our knowledge of how crowding stress suppresses the immune response is still limited.

The second-generation sequencing-based whole transcriptome analysis tool RNA-Seq allows the simultaneous and comprehensive characterization of all gene activities in a specific biological process. Accordingly, RNA-Seq has been employed to understand the immune events during a wider perspective on different biological settings in teleost fish [16,17]. RNA-Seq allows us not only to identify the gene and characterize the gene expression levels at the same time but also to elucidate host bacterial interactions [18]. To date, many immune studies in turbot have been performed using RNA-Seq following infection [19-21]. Using this sequencing technique, a large number of immune-related genes in several fish species were identified, including Japanese seabass (Lateolabrax japonicus ), turbot (S. maximus ), mud Loach (Misgurnus anguillicaudatus ) and large yellow croaker (Larimichthys crocea ) [22-25]. However, to our knowledge, no research has systematically characterized the immune strategies in turbot following attenuated Edwardsiella tarda vaccination and crowding stress at the transcriptomic level.

In this study, we aimed to determine the effects of crowding stress on the immune response in vaccinated fish. The juvenile turbot vaccinated with attenuated E. tarda were cultured at two different densities, high density (HD, 20.53 ± 0.05 kg/m²) as experimental groups and low density (LD, 5.25 ± 0.02 kg/m²) as control groups. After 5 weeks, the transcriptomes of the spleen and head kidney of vaccinated turbot from the two different density groups were analyzed. This transcriptomic study has enabled us to gain a better understanding of the effects of crowding stress on the immune response and identify a large number of differentially expressed target genes that bring us closer to the development of strategies designed to effectively combat pathogens. The identification of immune-relevant genes that could be potential markers for disease resistance is a means to establish a successful genetic breeding program. An understanding of the immune-related pathways in defense mechanisms is a relevant factor for enhancing the resistance of cultured fish to diseases.

Fish and experimental conditions

Healthy juvenile turbot weighing 150 ± 10 g were obtained from the farm of Shandong Oriental Ocean Sci-Tech Co. Ltd (Shandong, China), where the study was conducted. The fish were acclimated to the experimental environment for 15 days using a flow-through system. After vaccination with attenuated E. tarda , all of the fish were randomly assigned to two different stock densities: LD with 100 fish per tank (5.25 ± 0.02 kg/m²) and HD with 400 fish per tank (20.53 ± 0.05 kg/m²). Each density was tested with three replicates.

The fish were fed twice daily with commercial turbot feed (Ningbo Tech-Bank Co. Ltd, Zhejiang, China) and the daily ration was 1% of their body weight. During the test period, the water temperature was maintained at 16 ± 0.5°C with salinity at 28.2 ± 3%, dissolved oxygen at 8.0 ± 0.5 mg/l, and pH at 7.5 ± 0.3. The photoperiod was 12-h light and 12-h dark. No mortality was recorded during the experiment.

Sampling

In our previous study, the special antibody IgM in blood in vaccinated turbot approached maximum at the 5th week [26]. Therefore, 15 fish (5 fish per tank) were randomly collected after 5 weeks from each density group immediately after anesthetization with 0.05% tricaine methane sulfonate (MS-222, Sigma, St. Louis, MO, USA). The HD fish weighed 26.02 ± 0.05 kg/m² and the LD weighed 5.35 ± 0.02 kg/m² at sampling. The spleen and head kidney from 15 fish were dissected and pooled (five fish per pool). Samples were flashfrozen in liquid nitrogen and stored at −80°C prior to RNA extraction.

RNA extraction, library construction, and sequencing

Total RNA was extracted from the tissues using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. RNA degradation was monitored on 1% agarose gels. RNA concentration and purity were checked using a Nano-Drop 2000 spectrophotometer (Thermo, NewYork, MA, USA). Further, RNA integrity was assessed with the Agilent 2100 RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA). High-quality samples (RNA integrity number>7.0, RNA concentration>100 ng/μl) were selected for high-throughput sequencing. For each sample, equal amounts of RNA from the three replicates were pooled for RNA-Seq library construction.

After total RNA was extracted, mRNA was enriched by oligo (dT) beads. The purified mRNA was fragmented into short fragments using fragmentation buffer and then used as a template for first-strand cDNA synthesis using random hexamers and reverse transcriptase. The second-strand cDNA was synthesized using DNA polymerase I and RNase H. Then the cDNA fragments were purified with a QiaQuick PCR extraction kit, end-repaired, poly(A) was added, and then ligated to Illumina sequencing adapters. The ligation products were size selected with agarose gel electrophoresis and amplified by PCR. Then the cDNA library was constructed from the cDNA synthesized using NEBNext UltraTM RNA Library Prep Kit for Illumina (New England Biolabs [NEB], Ipswich, MA, USA). The cDNA libraries were sequenced using Illumina HiSeqTM4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

De novo assembly of sequencing reads

To obtain high-quality clean reads for the transcript assembly, raw reads were filtered according to the following rules: removing reads containing more than 10% unknown nucleotides; removing lowquality reads containing more than 50% of low quality bases; removing reads containing adapters. The resulting high-quality sequences were de novo assembled into contigs and transcripts with Trinity software (http://trinityrnaseq.sf.net) [27]. To reduce data redundancy, transcripts with a minimum length of 200 bp were assembled and mapped to the reference transcriptome using the Bowtie2 short reads alignment tool under default parameters [28].

Gene annotation

All assembled unigenes were used as queries in searching the NCBI non-redundant (Nr) protein database and SwissProt protein database using the BLASTX program. The cutoff E-value was set at 1e−5 and only the top gene ID and name were initially assigned to each unigene. The unigenes were functionally annotated by gene ontology (GO) analysis with Blast2GO software (E-value<10−5). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using KEGG Automatic Annotation Server (KASS) with default parameters [29].

Identification of differentially expressed genes

All clean sequencing reads from each of the cDNA libraries were mapped back to the transcriptome assembly using the software Bowtie2 with default settings. The reads aligned to each unigene in the alignment file were counted for each sample. These read counts were normalized as RPKM (reads per kilobase of transcripts per million fragments mapped) values [30] and further analysis to identify differentially expressed genes (DEGs) among the two different groups was conducted using a web tool DESeq. The false discovery rate (FDR) method was introduced to determine the threshold p-value in multiple tests to judge the significance of the difference in gene expression. If the FDR was <0.05 and there was at least a 2-fold difference in RPKM values among samples, the unigene was considered a significant DEG [31]

Experimental validation using qPCR

Twelve DEGs (Table 1) were selected randomly for validation of RNA-Seq data by qPCR using a SYBR Premix Ex Taq kit (Invitrogen) according to the manufacturer's instructions. The same RNA samples were used for both Illumina library synthesis and the qPCR verification assay. The first strand cDNA was obtained from 2 μg of total RNA using a PrimeScript 1^st strand cDNA synthesis kit (Takara, Dalian, China). Melting curve analyses were performed following amplification. The specific primers used for qPCR are listed in Additional Table 1, and β-actin was used as an endogenous control. The thermal profile for SYBR Green qPCR was 95°C for 5 min, followed by 40 cycles of 95°C for 15 s, 60°C for 30 s. The relative gene expression was analyzed using the comparative threshold cycle method established by Livak et al. [32].

Table 1: Primers used for qPCR validation of RNA-Seq.

Ethical standards

All experiments were performed according to the Regulations on Animal Experimentation at Yellow Sea Fisheries Research Institute, Qingdao, China. Following the experiment, all surviving fish were kept in laboratory.

De novo assembly and annotation

To obtain the turbot transcriptome expression profile after vaccination, four cDNA libraries were constructed using spleen and head kidney from the two groups of HD and LD turbot. The results of transcriptome have been submitted to NCBI (SRA accession no. SRP129900). A total of 447 million clean reads were obtained for subsequent analysis after eliminating low-quality sequence and adaptor sequences from the original data sequence by quality analysis. The clean reads showed a high quality with a Q20 of 97.81% and Q30 of 94.83% (Table 2).

Table 2: The quality of clean reads. HS and HK represent the spleen and head kidney from high density respectively; LS and LK represent the spleen and head kidney from low density respectively.

According to the basic characteristics of transcripts, the distribution characteristics of the actual coverage of the reads are as follows: the nearer the 5'- and 3'- ends are, the lower the average sequencing depth is, but the overall homogenization degree is relatively high (Figure 1). Transcript de novo assembly was performed for the clean reads by Trinity. The total number and length of unigenes were 41,136 and 52,423,775 bp, respectively. The maximal length of the unigenes was 23,178 bp with an average length of 1274 bp (N50: 2295), the GC content was 49.13% (Table 3). These results indicate that the sequencing data was high-quality, and the unigenes could be used for subsequent annotation analysis.

Figure 1: The homogeneity distribution of reads.

Table 3: Summary of de novo assembly of transcriptomic profiles of S. maximus.

To obtain comprehensive gene information, 41,136 unigenes were annotated by four databases including Nr, Clusters of orthologous groups for eukaryotic complete genomes (KOG), SwissProt and KEGG (Table 4). Annotation was matched successfully for 22,931 of the unigenes (22,806 in Nr, 20,101 in SwissProt, 16,035 in KOG, and 14,332 in KEGG) and 18,205 unigenes were without annotation. The distribution of 22,931 annotated unigenes from KEGG, KOG, Nr and SwissProt data are shown in Figure 2. In Nr annotation, 22,806 unigenes were matched to multiple species genomes, including Larimichthys crocea (31.58%), Stegastes partitus (25.76%), Oreochromis niloticus (5.68%), Notothenia coriiceps (5.16%), Cynoglossus semilaevis (4.92%), and others (26.88%) (Figure 3).

Table 4: The statistics of annotated unigenes of transcriptomic profiles in four databases.

Figure 2: The venn of four database annotation. The distribution of a total 22931 annotated unigenes in KEGG, KOG, Nr and Swissprot database.

Figure 3: The top 10 species match to unigenes of S. maximus. The x-axis indicates the number of unigenes and the y-axis indicates matched species. From left to right, the matching degree from high to low.

Identification and analysis of DEGs

All of the unigenes were analyzed with DEG-Seq, whose threshold was restricted to q-value <0.005, |log2 (fold-change)|>1. A total of 1155 genes showed statistically significant differential expression in two immunologic organs that were compared in HD and LD turbot. In detail, 397 unigenes were identified as DEGs containing 146 significantly upregulated and 251 significantly downregulated unigenes in spleen. Additionally, 758 unigenes were identified as DEGs containing 100 significantly upregulated and 658 significantly downregulated unigenes in head kidney (Figure 4). In the statistical figures, we observed that the downregulation levels of DEGs were significantly higher than the upregulation levels after attenuated E. tarda administration. Among the 1155 DEGs, a group of 84 DEGs was shared between the two organs, 16 up- and 66 downregulated, and 2 either up- or downregulated depending on the organ.

Figure 4:The statistics of different expression gene in spleen and head kidney. The red represents upregulated expressed genes and the green represents downregulated expressed genes.

GO enrichment and pathway analysis

To further investigate the function of the DEGs, all were classified into three GO categories (biological process, molecular function, and cellular component) by Blast2GO. A total of 36 GO terms including 16 biological process terms, 12 cellular component terms and 8 molecular functions terms were associated with spleen DEGs (Figure 5). A total of 42 GO terms including 19 biological process terms, 13 cellular component terms, and 10 molecular functions terms were associated with head kidney DEGs (Figure 6). Analysis of level 2 GO term distribution showed that cellular process, binding, and cell were the most common annotation terms in the three GO categories. Among these, the immune system process was retained for further pathway analysis, although it was not the most significant.

Figure 5: Gene ontology (GO) enrichment analysis of the differently expressed genes. A total of 36 GO terms including 16 biological process terms (1-16), 12 cellular component terms (17-28) and 8 molecular functions terms in spleen (29-36). Note: 1 Locomotion; 2 immune system process; 3 response to stimulus; 4 biological adhesion; 5 multi-organism process; 6 developmental process; 7 multicellular organismal process; 8 localization; 9 reproduction; 10 reproductive process; 11 biological regulation; 12 metabolic process; 13 single-organism process; 14 signaling; 15 cellular component organization or biogenesis; 16 cellular process; 17 binding; 18 molecular function regulator; 19 nucleic acid binding transcription factor activity; 20 transporter activity; 21 catalytic activity; 22 structural molecule activity; 23 signal transducer activity; 24 molecular transducer activity; 25 extracellular region part; 26 extracellular region; 27 extracellular matrix; 28 cell junction; 29 supramolecular fiber; 30 membrane part; 31 membrane; 32 cell; 33 cell part; 34 macromolecular complex; 35 organelle; 36 organelle part.

Figure 6: Gene ontology (GO) enrichment analysis of the differently expressed genes. A total of 42 GO terms including 19 biological process terms (1-19), 13 cellular component terms (20-32) and 10 molecular functions (33-42) terms in head kidney. Note: 1 biological regulation; 2 immune system process; 3 biological adhesion; 4 multicellular organismal process; 5 developmental process; 6 behavior; 7 growth; 8 locomotion; 9 cellular component organization or biogenesis; 10 reproduction; 11 reproductive process; 12 multi-organism process; 13 rhythmic process; 14 localization; 15 single-organism process; 16 cellular process; 17 metabolic process; 18 signaling; 19 response to stimulus; 20 binding; 21 nucleic acid binding transcription factor activity; 22 molecular transducer activity; 23 transcription factor activity, protein binding; 24 transporter activity; 25 molecular function regulator; 26 catalytic activity; 27 antioxidant activity; 28 signal transducer activity; 29 structural molecule activity; 30 extracellular region; 31 extracellular region part; 32 extracellular matrix component; 33 membrane-enclosed lumen; 34 extracellular matrix; 35 organelle; 36 cell junction; 37 membrane part; 38 cell; 39 cell part; 40 membrane; 41 organelle part; 42 macromolecular complex.

A total of 401 unigenes were homologous to immune-relevant genes based on the KEGG annotation and were classified into 15 functional pathways (Table 5). Furthermore, through immune-related DEG screening, 10 DEGs (Table 6) were assigned to six immune-relevant pathways, including the RIG-I-like receptor signaling pathway, cytosolic DNA-sensing pathway, Toll-like receptor (TLR) signaling pathway, NOD-like receptor signaling pathway, chemokine signaling pathway, and platelet activation. Among these immune-relevant pathways, we focused on significantly changed signaling pathways (P<0.05) such as the TLR signaling pathway, cytosolic DNA-sensing pathway, and platelet activation.

Table 5: KEGG pathway enrichment analyses of the immune related differently expressed genes and annotation genes. The symbol “-” represents the number of differently expressed genes were zero. The number of annotation genes represents expressed genes in the corresponding pathway in this study.

Table 6: The information of immune related DEGs. Bold values indicate significantly changed relative to the low density groups. The positive numbers indicate up-regulation and the negative numbers indicate down-regulation.

Validation of RNA-Seq profiles by qPCR

To validate the RNA-Seq data, 12 DEGs were selected randomly for qPCR. Although, the test unigenes displayed different expression levels (Figure 7), in general, the qPCR results showed a positive correlation with RNA-Seq data, indicating the reliability and accuracy of the RNASeq expression analysis.

Figure 7: Validation of RNA-seq data by using qPCR. X-axis, pairwise comparison groups; Y-axis, fold change in gene expression. CCL20, CISH, IFNA, SOCS3, IL1R2and PIK3R were detected in head kidney; CXCL10, EPOR, CCL19, CXCL13, IL18 and IL1βwere detected in spleen. The qPCR data were given in mean ± SD, n=3.

RNA-Seq technology uses RNA-Seq sequence data to provide a complete picture of the transcriptome under study (global annotation) and simultaneously provide quantitative data related to differential gene expression. In the present study, we have obtained the transcriptome comparative data of S. maximus spleen and head kidney under two different culture densities and have acquired 41,136 unigenes after assembling and filtration. RNA-Seq technology can annotate successfully in the absence of a genome dependent upon the other existent genome resources and expressed sequence tag collections [33]. Although the turbot genome was released and was available in the European Nucleotide Archive, the mapping percent was very low; therefore, we used the BLASTX program to annotate the unigenes in an approach adopted by other studies [20].

In Nr annotation, only 1123 unigenes were matched to Cynoglossus semilaevis (4.92%) genomes. The origin of flatfish has been the focus of long-term controversy [34]. Phylogenetic studies show low interfamily/ suborder statistical support among flatfish, especially between the Psettoidei and Pleuronectoidei, and the low phylogenetic support is considered evidence of a polyphyletic origin [35]. Genetic evidence on the diversification of flatfish have been reported using genome sequencing of the turbot and the tongue sole [36,37], suggesting that different strategies are adapted for demersal life [34].

A total of 1155 DEGs were categorized into many categories based on GO annotation and pathway analyses. These categories included the immune system, which is the key category we focus on. Below we highlight three key constituents of the immune system, the TLR signaling pathway, the cytosolic DNA-sensing pathway, and platelet activation, which were all influenced significantly by crowding stress.

TLR signaling pathway

Recognition of pathogen-associated molecular patterns (PAMPs) is essential for the activation of innate immunity and it is carried out by TLRs. TLRs, specific families of pattern recognition receptors, are responsible for detecting microbial pathogens and generating innate immune responses. Pathogen recognition by TLRs provokes rapid activation of innate immunity by inducing the production of proinflammatory cytokines and upregulation of costimulatory molecules as well as priming the adaptive immune system [38-40]. Efficient immune responses depend upon a close interaction between the innate and adaptive immune systems. TLRs serve as an important link between the innate and adaptive immune responses [41]. Importantly, TLRs participate in the direct regulation of the adaptive immune response, possibly as costimulatory molecules [42]. Dendritic cells stimulated directly or indirectly by TLRs from pathogens mature into a specific form and are able to activate a specific immune response that is appropriate for the elimination of the pathogen [43]. In diseases, TLRs also participate in inflammation and immune responses that are driven by self-, allo- or xeno-antigens [40].

In our study, inflammatory cytokines such as interleukin (IL)-1β and interferon (IFN)- were markedly upregulated in spleen when turbot suffered crowding stress. It is well-known that the protective properties of these inflammatory factors remove the adverse stimulus in the first line of an organism’s defense against pathogen invasion, but a stronger inflammatory response can lead to tissue damage as well [44]. IL-10 is also (anti) inflammatory cytokine, may decrease the production of pro-inflammatory cytokines and keep the immune system in balance. In our study, there was no significant upregulation of IL-10, indicating that the inflammatory cytokines were overexpressed. For similar reasons, markedly upregulated IL-10 inhibited IL-1β expression in Schizothorax prenanti after poly(I:C) challenge [44]. Chemokines have a broad spectrum of effects on innate immunity, such as becoming a promoter to adhere various types of leukocyte to inflammatory loci in early immunity [45,46]. Research has shown that turbot CXC chemokines were dramatically and rapidly induced at 5 h after challenge and the expression of CXC chemokine was not observed until after 12 h [42]. In our study, CXCL9 and CXCL10 were upregulated at the 5th week, which may indicate the overexpression of CXC chemokines. In addition, we found that PIK3R was markedly downregulated in spleen. PI3k inhibition leads to the suppression of IL-10 but induced IL-1β expression [47]. The PI3ks are lipid kinases that are activated by receptor tyrosine and suppress cell apoptosis and promote cell growth and play important roles in the immune response [48]. Downregulated PIK3R signifies that crowding stress is disadvantageous to the immune response.

In head kidney transcriptomic data, no inflammatory cytokines and costimulatory molecules were observed to be upregulated or downregulated. However, PIK3R and MAPK1-3 were dramatically downregulated in the TLR signaling pathway. The downregulation of MAPK1-3 indicates the suppression of extracellular signal-regulated kinases (ERK), which are widely expressed protein kinase intracellular signaling molecules involved in maintaining cell survival and mediating intracellular signal transduction in response to a variety of stimuli [49-51]. Some studies show that upregulated ERK plays a key role in environmental adaptation processes such as protecting cells from stress [52,53]. In our study, downregulated ERK signifies that crowding stress damages spleen cells.

Cytosolic DNA-sensing pathway

Specific families of pattern recognition receptors are responsible for detecting foreign DNA from invading microbes or host cells and generating innate immune responses. Cytosolic DNA sensors detect intracytoplasmic double-stranded DNA and signal through a central adaptor protein, stimulator of IFN genes (STING), which recruits and activates the downstream TBK1/IRF3 signaling axis to induce the production of type IFN [54,55]. Type IFN are critical for innate immune defense against viral infections [56,57]. Furthermore, cytosolic DNA also activates STING-dependent and STINGindependent autophagy to promote host defense against bacterial infection [58,59]. Studies by Rengaraj et al. have shown that the cytosolic DNA-sensing pathway is the most important immune pathway in necrotic enteritis-induced chickens [60]. Similar results were obtained in our study, the DEG enrichment in the cytosolic DNAsensing pathway are significantly up regulated in spleen. IL-1β, IFN-α, and CXCL10 are up regulated, but IL-18 is down regulated in the cytosolic DNA-sensing pathway. Studies have confirmed that IL-18 is important for regulating the Th1 and Th2 immune response [61]. The down regulation of IL-18 will attenuate the ability of the spleen cell immune response. In head kidney, DEGs are not enriched in the cytosolic DNA-sensing pathway indicating that spleen and head kidney play different roles in the immune response.

Platelet activation

It has become obvious that platelets are not only cell fragments that plug the leak in a damaged blood vessel but also key components in the innate immune system, which is supported by the presence of TLRs on platelets [62]. Recently platelets have been reported to express TLR2, 4 and 9, reinforcing their role as primitive immune cells in host defense [63,64]. As platelets are usually the first cells to respond to a wound, they have an important role in the immune response to infection through their activation [65,66]. Platelets acting as components of the innate immune system detect the presence of infectious agents and coordinate the response to the pathogen [67]. When activated, platelets secrete over 300 proteins [68] including the bioactive molecules adenosine diphosphate (ADP) and serotonin to reduce blood loss, cytokines and chemokine recruiting leucocytes to deal with any potential infection, and antimicrobial peptides to kill pathogens. The activation of platelets leads to the secretion of antimicrobial peptides, although many bacteria have become resistant to these peptides [69].

In our study, AC and TLN were significantly downregulated in the platelet activation pathway in head kidney. TLN is a master regulator of platelet integrin activation in vivo. If there is a lack of TLN, the platelet aggregation function will be damaged severely and integrin αIIbβ3 activation will also be inhibited [70,71]. TLN downregulation decreases ligand binding to cellular integrins resulting in restrained platelet and leukocyte function. AC, the second messenger cyclic adenosine monophosphate (cAMP) and cAMP-dependent protein kinase A are important mediators in determining the response of a cell to external stimuli [72,73]. Activated AC increases the intracellular cAMP level, which enhances the host defense [74]. The downregulation of AC will attenuate the response of a cell to external stimuli by decreasing cAMP levels.

Apart from the nine genes mentioned above, one gene DDX3X is significantly downregulated in the RIG-I-like receptor signaling pathway. DDX3X participates in regulating mRNA translation and some signaling pathways [75,76]. DDX3X directly regulates Wnt/β- catenin signaling and INF-induced signal pathways [75,77]. DDX3X is responsible for detecting viral pathogens and generating innate immune responses. Vaccinia virus protein K7 can combine with DDX3X to suppress the innate immune response [78]. Downregulated DDX3X indicates that the innate immune response is inhibited.

Using RNA-Seq-based transcriptome profiling, we assessed, for the first time, the transcriptomic responses of turbot cultured in different densities following vaccination with attenuated E. tarda. A total of 1155 DEGs were annotated in four databases. We identified immunerelated DEGs compared in high and low stock densities. After RNASeq data analysis and qPCR validation, we observed that immunerelated DEGs were primarily enriched in the TLR signaling pathway and cytosolic DNA-sensing pathway in spleen, and platelet activation occurs in head kidney. The inflammatory cytokines, IL-1β, IFN-α, CXCL10, and CXCL9, and signal-regulated cytokines, PIk3R, MAPK1-3, AC, and TLN, are involved in these pathways. Overexpressed inflammatory cytokines and downregulated signalregulated cytokines show that cells suffer damage and the immune response is restrained when turbot is stressed due to crowding. Although further functional studies are needed to characterize the key immune factors governing the turbot immune responses, these results can provide a foundation to evaluate the relationship between immunosuppression and crowding stress.

This work was supported by the National Natural Science Foundation of China (No.31402315), the Central Public-interest Scientific Institution Basal Research Fund, CAFS [grant number 2017HYZD04], the Modern Agriculture Industry System Construction of Special Funds (CARS-50-G10) and Qingdao Shinan District Science and Technology Bureau (2016-3-006-ZH).