Gill Na+/K+-ATPase activity
The Na+/K+-ATPase activity of gill tissue in each treatment group is illustrated in Fig. 2. In AA, the Na+/K+-ATPase activity was significantly higher than in the CK (P<0.01). At the time of low-temperature stress prolongation, the Na+/K+-ATPase activity of BB and CC was significantly lower than AA (P<0.05). The differences between CC and CK were significant (P<0.05), while the Na+/K+-ATPase activity of BB was similar to CK (P>0.05).
Transcriptome profiles and annotation
A total of 593.19 million raw reads were obtained by sequencing 12 cDNA libraries. After removal of adaptor and primer, ploy-N-containing reads, Ribosome RNA (rRNA), and low-quality reads, about 585.51 (98.70%) million clean reads were obtained (Table 1). The uniquely mapped percentages of these transcripts ranged from 83.26% to 84.62%, while the percentage of multiple mapped transcripts ranged from 2.30% to 2.80%. The average Q20 and average Q30 of AA reached 97.98% and 97.92%, the average GC content of the experimental group was 45.27%, and the average Q20 and Q30 were 95.74% and 89.93%. The error rate was 0.03%, indicating that the clean reads can be used for subsequent biological information analysis.
Differentially expressed genes
During the RNA mapping, at least 83.26% of the clean reads were aligned to the tilapia reference genome. Statistical analysis performed with the DESeq2 package identified the differentially expressed genes between the CK and the cold treatments (AA, BB, CC). A total of 12,448 DEGs were identified between CK and AA, BB, and CC (Fig. 3). Compared to CK, there were 1,784 DEGs, of which 792 were upregulated and 992 were downregulated in AA (Fig 3A). In BB, 4,883 DEGs were identified, of which 1,827 were upregulated and 3,056 were downregulated (Fig 3B). In CC, 5,781 DEGs were identified, of which 1,924 DEGs were upregulated and 3,857 were downregulated (Fig 3C). The response to low temperature was strong in terms of the total number of DEGs and also in the degree of transcriptional changes. A total of 478 DEGs in AA (212 over-expressed and 266 under-expressed), 1,642 DEGs in BB (511 over-expressed and 1,131 under-expressed), and 2,310 DEGs in CC (554 over-expressed and 1,756 under-expressed) showed large FCs (|log2FC| ≥ 2). The global gene expression profiles of tilapia in AA, BB, and CC showed large differences when compared to CK under low-temperature stress conditions. In addition, the number of DEGs increased with treatment time. Through sequence alignment, many functional genes were recognized, such as cell division cycle 42 (cdc42), aquaporin families (aqps), FK506-binding protein 5 (fkbp5), heat shock protein family A (Hsp70) member 4 (hspa4a), and toll-like receptor 5 (tlr5). These genes may increase our understanding of the response of tilapia to low-temperature stress.
Gene Ontology (GO) and KEGG pathway enrichment
GO enrichment and KEGG pathway enrichment analysis were performed using the corresponding annotations of genes to understand the functional relevance of the DEGs. Compared to CK, 49,759 DEGs could be assigned by GO classification. Of these, there were 6,414, 19,656, and 23,689 DEGs in the CK/AA, CK/BB, and CK/CC comparisons, respectively (Fig. 4). Compared to CK, the top five significant annotated GO enrichments in the cold treatment groups were all binding (GO:0005488), cellular process (GO:0009987), single-organism process (GO:0044699), metabolic process (GO:0008152), and catalytic activity (GO:0003824). The genes enriched to these five pathways accounted for 42.33%, 40.39%, and 40.57% of the total genes in each group. The classified genes (groups CK, AA, BB, and CC) produced 599 different pathways, including organismal systems, environmental information processing, human diseases, metabolism, and cellular processes. In the AA/CK comparison, 564 DEGs were assigned to 203 pathways (Fig. 5). In the CK/BB comparison, 1,023 DEGs were assigned to 227 pathways. Compared to CK, 1,234 DEGs were assigned to 216 pathways in CC. Compared to CK, the top three KEGG pathways with the largest number of enriched genes in AA, BB, and CC were all cytokine-cytokine receptor interaction (ko04060), metabolic pathways (ko01100), and cell adhesion molecules (CAMs) (ko04514). According to the KEGG function annotations, a total of 55 significantly enriched pathways (P-value < 0.05) were identified in the CK/AA, CK/BB, and CK/CC comparisons.
Trend analysis of differentially expressed genes
A total of 1,305 DEGs were performed by STEM to study the transcriptome changes of gills during low-temperature stress. A total of 1,305 were classified into eight profiles (Fig. 6), and among them, two significant profiles were identified (P<0.01). These included one upregulated profile (profile covering 445 DEGs) and one downregulated profile (profile covering 761 DEGs).
Validation of RNA-Seq Profiles by RT-PCR
To validate the gene profiles from RNA-seq, the expression of five genes was measured by real-time PCR analyses (Fig. 7). The melting-curve analysis showed a single product for all selected genes. The overall consistency of expression from RNA-seq and RT-PCR for the majority of genes revealed the sequence assembly, and the expression analysis of genes provided transcriptome information.