Sampling, high throughput sequencing of RNA libraries, quality control of raw data
There was no dead fish found in this experiment, all the fish in the tanks were vigorous, healthy, and disease-free. In the hypoxic group, we found that T. rubripes exhibited lying on the bottom of the pool and swimming slowly. In the meanwhile, the feeding rate of fish was also reduced.
After the total RNA of the samples was extracted from gills tissue, a total of 6 RNA (cDNA) libraries were constructed in this study (3 normoxic groups and 3 hypoxic groups). High throughput sequencing results of 12 samples, including normoxic groups and hypoxic groups are summarized in Supplementary Information Table S1-S4. The results showed that the average size of insert fragments ranged from 284.1-294.3bp in these libraries (Supplementary Information Table S1). As listed in Supplementary Information Table S2, a total of 272,683,668 raw data was obtained, ranging from 41,836,262 − 49,513,198 for each group. Following the pre-processing steps of adapter clipping, poly-N, and low reads quality filtering, a total of 263,594,614 clean reads and 39.55 Gb clean bases were generated among the 6 libraries with Q20 (97.47% − 97.79%) and Q30 (93.07% − 93.87%) were calculated. And the GC content of all samples were all exceeded 50%. These results indicated the sequencing data generated in this study were reliable and could be used for further analysis.
Mapping and counting reads
The high-quality reads obtained were compared with the reference genome of the T. rubripes, the results showed that the total mapping rate and unique mapping rate were more than 93.28% and 78.26%, respectively, and the proper mapping rate ranged from 86.11%-88.13% (Supplementary Information Table S3). As shown in Supplementary Information Table S4, we found that the aligned reads in genomes were mainly distributed in exon (87.44% − 88.98%), only 4.16% − 5.08%, and 3.01% − 3.28% reads were localized in intron and intergenic, respectively. Furthermore, Analysis of mapping results using FastQC and MultiQC manifested that the duplicated level of the mapping reads ranged from 23.8% − 28.4% (Supplementary Information Table S5). All these results showed that our sequencing data had a good mapping rate and could be used for further analysis.
Analysis of differentially expressed genes (DEGs)
After mapping these unique clean reads to specific genes and analyzing gene expression by HTSeq, the percent of assigned reads were shown in Supplementary Information Table S6. Differential expression genes (DEGs) analysis of the two groups was performed using the DESeq2 R package, results showed that a total of 683 genes were identified as DEGs among the treatment groups using the criteria of |log2foldchange| >= 1 & p.adj < 0.05 (including 378 up-regulated and 305 down-regulated genes) (Fig. 1a). Among them, only 143 genes had a difference greater than two-fold (including 84 up-regulated genes and 59 down-regulated genes) (Fig. 1b). All the following analysis was based on the analysis with |log2foldchange| >= 1 & p.adj < 0.05, and the detailed DEGs information obtained from the dataset was present in Supplementary Information Table S8. In order to visualize, we used the DESeq2 to extract transformed count data with the method of variance stabilizing transformation. Based on the normalized expression data of DEGs to print heatmap and PCA diagram (Fig. 2, 3), we found that 6 samples could be well clustered into two branches (normoxic group and hypoxic group) according to the treatment condition, indicating that the experimental data was reliable and the sample selection was reasonable in this study.
Annotation and functional analysis on DEGs in gills
Functional classification of DEGs using GO enrichment analysis
To identify the most highly regulated hypoxia-responsive pathways in response to chronic hypoxic stress, a GO enrichment analysis was performed on the up-and down-regulated genes from three aspects: biological processes (BP), molecular functions (MF), and cellular component (CC), respectively. The list of all the GO terms associated with the DEGs was provided in Supplementary Information S9. Results showed that the 143 DEGs in gills of hypoxia-treated fish were classified into 3039 GO terms, including 2499 terms of BP (82.23%), 318 terms of MF (10.46%), and 222 terms of CC (7.31%). In gills, we observed that several GO terms related to metabolic and stress response in the BP, including “drug metabolic process” (GO:0017144), “response to xenobiotic stimulus” (GO:0009410), “estrogen metabolic process” (GO:0008210), “xenobiotic catabolic process” (GO:0042178) and “toxin metabolic process” (GO:0009404) (Fig. 4a) showed the highest enrichment. Moreover, under the classification of MF (Fig. 4b), “heme binding” (GO:0020037), “tetrapyrrole binding” (GO:0046906), and “iron ion binding” (GO:0005506) were the most highly represented terms with high gene numbers. In contrast, there were no GO terms significantly enriched in the category of CC (Fig. 4c). In order to intuitively show our results, we clustered all these enriched terms and interesting found that GO terms with similar functions could be well clustered into one major functional category (Fig. 5), which were grouped into eight categories in the present study, including “arachidonic vitamin catabolic xenobiotic”, “regulation mitochondrial depolarization apoptotic death” and “erythrocyte homeostasis fate commitm” and so on.
Moreover, it was noteworthy that significant changes in the transcription of genes involved in several GO terms related to the metabolic process were significantly enriched, including “drug metabolic process” (GO:0017144), “estrogen metabolic process” (GO:0042178), “xenobiotic catabolic process” (GO:0042178), “xenobiotic metabolic process” (GO:0006805). After clustering all enriched GO terms, all of the above metabolism-related GO terms were enriched to “arachidonic vitamin catabolic xenobiotic”. We further analyzed these DEGs and observed that the mRNA level of cyp1a1 and ugt1a1 were both up-regulated (Fig. 8a, 8b) in those terms. Thus, we speculated that these genes related to metabolism may play an essential role in the response to chronic hypoxic stress.
Functional classification of DEGs using KEGG Enrichment Analysis of DEGs
To identify biological pathways that were affected in the gills of T. rubripes during exposure to hypoxia, 143 DEGs were mapped to the KEGG database (Kanehisa et al. 2017) and the enriched pathways were listed in Supplementary Information Table S10. Meanwhile, the KEGG enrichment analysis output was plotted as a histogram for the top 50 enriched KEGG pathways (Fig. 6). In the gills transcriptome of the T. rubripes under chronic hypoxic stress, immune regulation–related pathway “cytokine-cytokine receptor interaction” (hsa04060, 10 genes) was found to be most significantly enriched using KEGG enrichment analysis, and the expression levels of genes encoding key proteins ccl19, cl19l1 (ccl13) and inflammatory cytokine IL-17F were up-regulated (Fig. 8c-e) under chronic hypoxic stress. Apart from that, “Steroid hormone biosynthesis” and “Viral protein interaction with cytokine and cytokine receptor” were also significantly enriched under chronic hypoxia. This result indicated that hypoxia substantial changes in immune-related defense in the gills of T .rubripes.
Immune-related pathways appear to be a positive response to chronic hypoxic stress using GSEA enrichment analysis
Gene Set Enrichment Analysis (GSEA, http://software.broadinstitute.org/gsea/) is a computational method used to determine whether a set of a priori defined genes shows statistically significant, consistent differences between two samples (Subramanian et al. 2005). To better understand the regulatory mechanism of hypoxia, we further performed GSEA enrichment analysis for both acute and chronic hypoxic stress (Shang et al. 2022a). We found that several immune-related biological processes, apoptosis-related pathway and some process related oxygen-binding were significantly altered in the gills of hypoxic stress through transcriptome technology and functional enrichment analysis.
Based on the GO gene set of T. rubripes, there were 175 GO terms were significantly enriched (q values < 25%) under acute hypoxic stress (Supplementary Information Table S11), including 50 up-regulated terms (28.6%) and 125 down-regulated terms (71.4%). Under chronic hypoxic stress, only 27 terms were significantly enriched (Supplementary Information Table S12), including 16 up-regulated terms (59.3%) and 11 down-regulated terms (40.7%). The heatmap showed that immune-related terms including “positive regulation of interleukin-10 production”, “positive regulation of interleukin-4 production”, “chemokine-mediated signaling pathway”, “CCR chemokine receptor binding”, “G protein-coupled receptor signaling pathway”, “inflammatory response”, “cytokine activity”, “cellular response to interleukin-1”, “epidermal growth factor receptor signaling pathway”, “positive regulation of protein kinase B signaling”, and “small GTPase mediated signal transduction” had different changes at different DO levels. They were all inhibited under the acute hypoxic group and activated under the chronic hypoxic group (Fig. 7a). Among them, “CCR chemokine receptor binding” was significantly enriched under hypoxic stress.
In the KEGG gene set of T. rubripes, there were 69 (27 up-regulated pathways, 42 down-regulated pathways) and 79 (27 up-regulated pathways, 52 down-regulated pathways) pathways were significantly enriched (q values < 25%) under acute and chronic hypoxic stress, respectively (Supplementary Information Table S13, S14). Through the heatmap we also interesting found that some immune-related pathways showed interesting changes under acute and chronic hypoxic stress. For example, “Tight junction”, “Leukocyte transendothelial migration”, “Complement and coagulation cascades”, “Cytosolic DNA-sensing pathway”, “Cytokine-cytokine receptor interaction”, “JAK-STAT signaling pathway”, and “Endocytosis” were inhibited under the acute hypoxic group, while they were all activated under chronic hypoxic group (Fig. 7b). Among them, “Cytokine-cytokine receptor interaction” was significantly enriched in the present study.
Gene set enrichment analysis (GSEA) of other pathways related to hypoxia
Cell death is critical for fundamental physiological processes such as development, immunity, and tissue homeostasis. In the current study, the apoptosis-related pathway “Ferroptosis” was inhibited under acute hypoxia but significantly activated chronic hypoxic stress (Fig. 7b). Meanwhile, the result showed that ferroptosis was enriched with 20 up-regulated genes. Among them, gclc, gclm, and LOC101073513 (frim) were the DEGs (Fig. 8f-h).
Moreover, “oxygen binding” and “oxygen carrier activity” were activated under acute and chronic conditions, and they were both significantly enriched under the chronic hypoxic group (Fig. 7a). In the current study, an oxygen-binding hemoprotein, mb (myoglobin, myg), was enriched in both pathways and showed increased in response to chronic hypoxia (Fig. 8i). Overall, our results may important to explain the function of the gills in T. rubripes under hypoxia.
Validation of RNA-Seq Data with qRT-PCR
Validation of RNA-Seq data was performed using qPCR on randomly selected genes. As shown in Fig. 9, the trends in the expression patterns of these genes detected by qPCR were consistent with the most genes obtained by RNA-Seq data analysis. These data confirm the reliability of the RNA-Seq results used in this study, and the follow-up research analysis could be continued.