Transcriptome sequencing, assembly and GO analysis
To investigate the effects of hypoxia on their molecular pathways and stress responsive genes involved in this processes. Six RNA-seq libraries were constructed using gills from three clams under hypoxia treatment (H1, H2 and H3) and three clams in normoxia control (N1, N2 and N3) to generate a reference transcriptome of R. philippinarum. Library sequencing yielded a total of 297,734,780 raw reads with a mean of 125 bp from six libraries. After filter of low-quality reads, 283,487,670 (95.2%) high quality reads were retained and de novo assembled. Clean reads obtained for H1, H2, H3, N1, N2, and N3 have been submitted to the SRA database in NCBI (accession numbers SRP151767). A total of 290,406 unigenes were obtained with an average length of 935 bp (Minimum length: 201 bp, Maximum length: 46,346 bp, N50 length: 1389 bp) (Table 1). The length distribution of the unigenes is shown in Figure 1. In the distribution of sequences annotated in Nr, 46.5% unigenes have similarities with Crassostrea gigas, followed by Lottia gigantea (10.3%), Hydra vulgaris (7.0%), Aplysia californica (6.5%), Branchiostoma floridae (2.2%), and others (27.4%) (Figure 2).
All the unigenes were enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms (Table 2). Differentially expressed genes involved in peroxisome, mitochondrial envelope, and transferring phosphorus-containing groups were identified in the cellular components (Additional file 1). For the molecular functions, genes related to catalytic activity, binding, oxidoreductase activity, peptidase activity, peroxiredoxin activity and transferase activity highly represented among the molecular functions (Figure 3, Additional file 2). In this study, several hypoxia regulated functions such as regulation of apoptotic process, immune system process, antioxidant activity, immune response, tumor necrosis factor receptor binding, peroxidase activity, and oxidoreductase activity were shown in Table 3. Decrease in defense response and response to stress activity might be related to regulation of response to oxidative stress, cell surface receptor signaling pathway, and reactive oxygen species (Additional file 1). Moreover, enzymes involved in immune response, antioxidant activity, peroxidase activity, oxidoreductase activity, peroxiredoxin activity, cytokine receptor binding and regulation of apoptotic process were enriched under hypoxia conditions (Table 3), indicating that R. philippinarum might acclimated to hypoxia through physiological, immune and defense response to stress. Additionally, a number of genes were involved in important categories, including cellular nitrogen compound biosynthetic process, organonitrogen compound catabolic process, cellular protein metabolic process, and macromolecular complex were also identified (Additional file 1), which may play potential roles in the hypoxia tolerance of R. philippinarum.
KOG analysis and KEGG classification
Eukaryotic Orthologous Groups (KOG) analysis was performed to classify all unigenes into different functional categories (Figure 4). Annotated genes in these categories are probably related to energy production and conversion, signaling transduction, intracellular free amino acid transport and metabolism, cell signaling pathways, and regulation of ionic content and cell volume. A total of 22,862 unigenes were classified into 231 different pathways through KEGG pathway analysis, such as metabolism, cellular processes, and organismal systems. The metabolic pathways with most unigenes divide into 11 groups, including carbohydrate metabolism, energy metabolism, lipid metabolism, amino acid metabolism, glycan biosynthesis and metabolism, and so on (Table 4). All DGEs from the comparison of Hypoxia and control group were mapped to KEGG for identification of the biological pathways in response to hypoxia. It is indicated that immune-related pathways were enriched, such as NF-κB signaling pathway, Toll-like receptor signaling pathway, complement and coagulation cascades, chemokine signaling pathway, and PI3K-Akt signaling pathway (Additional file 3).
Hypoxia responsive genes and molecular pathways in Manila clam
In this study, 148,358,474 and 135,129,196 qualified Illumina read pairs were obtained from the hypoxia challenged group (H1, H2, and H3) and normoxia control group (N1, N2, and N3) of R. philippinarum. RNA-seq analysis revealed that only 75 unigenes were differentially expressed genes (DEGs) under hypoxia challenge (qvalue<0.005 and |log2(foldchange)|>1). Among these, 32 DEGs were down-regulated and 43 DEGs were up-regulated (Figure 5A, Additional file 4). Heatmap of cluster analysis of DEGs from the transcriptomes of hypoxia and normoxia in R. philippinarum was shown in Figure 5B.
A number of hypoxia responsive pathways, including the HIF-1 signaling pathway, NF-kappa B signaling pathway, AMPK signaling pathway, Apoptosis, Jak-STAT signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, Rap1 signaling pathway, Ras signaling pathway, TNF signaling pathway, PPAR signaling pathway, Notch signaling pathway, Wnt signaling pathway, Calcium signaling pathway, cAMP and cGMP/PKG signaling pathway were indentify in this study (Additional file 1). In HIF-1 signaling pathway, the NF-κB, PI3K and p70S6K are upregulated, which could increase the HIF-1α mRNA or HIF-1α (Figure 6). Consequently, the downstream gene of the HIF signaling pathway, enolase (ENO1), is significantly upregulated in the hypoxia challenged clams, which promote anaerobic metabolism and causing reduced oxygen consumption.
Nuclear factor kappa B (NF-κB) activation starts with the phosphorylation, ubiquitination and subsequent proteosomal degradation of IκB (Figure 7). NF-κB complex could enter the nucleus with the IκB degradation, and then it turn on the expression of specific genes with DNA-binding sites for NF-κB. Genes activated by NF-κB result in physiological response, for instance, an immune/inflammatory response, cellular proliferation and survival response. As shown in Figure 7, the NF-κB turns on expression of its own repressor IκBα, and then the IκBα inhibits NF-κB, and thus forms a feedback loop, which regulates the activity of NF-κB [38]. These pathways play potential roles in the signal transduction, intracellular transduction, and sensing of stress signals, endocrine system and immune response of R. philippinarum under hypoxia stress.
The effects of hypoxia on SDH, LDH and AKP enzyme activity
As shown in Figure 8A&B, the succinate dehydrogenase (SDH) and the lactate dehydrogenase (LDH) activity in hepatopancreas of R. philippinarum increased first and reach a peak at 5 d, and then decreased at 8 d. In hypoxia group, the activity of SDH in hepatopancreas at 8 d significantly lower than that in control group, while the activity of LDH in hepatopancreas at 8 d was higher in hypoxia group compared to that in control group, which may be caused by clams switch to anaerobic metabolism at 8 d hypoxia stress. The activity of alkaline phosphatase (AKP) in hepatopancreas was significantly increased to a maximum at 5 d and then decreased at 8 d under hypoxia challenge, while almost no change of the AKP activity in control group was observed (Figure 8C). In the gills, however, no significant change in AKP activity occurs between hypoxia challenge and control group from 0 to 8 d (8 days), and an organ/tissue-specific response was observed, suggesting AKP might be involved in defense response and immune function in hepatopancreas of R. philippinarum.
Time-course of transcripts detected in hypoxia challenged clams
To corroborate the DEGs identified in the RNA-seq analysis, hypoxia responsive genes were selected from the DEGs for quantitative real-time PCR (qPCR) validation from the genes of interest based on their functions. The expressions values of 14 genes were analyzed in two tissues of Manila clam at different times of hypoxia stress. The expressions of 14 hypoxia responsive genes were significantly different between hypoxia challenged group and control group at different time 2, 5, and 8 d under hypoxia stress (Figure 9). Under hypoxia conditions, the Cytochrome P450, HSP70, Peroxisomal membrane protein (PEX), Sterol carrier protein (SCP), Glutathione peroxidase (GPx), Fibropellin-1, and Inhibitor of apoptosis protein (IAP) are expressed at a significant higher level at 2, 5, and 8 d hypoxia exposure than in control groups. While the Cadherin, Calmodulin, Defensin, Ubiquitin-conjugating enzyme E2, Serine/threonine-protein phosphatase, E3 ubiquitin-protein ligase and Ras-related proteins are expressed at a lower level at 2, 5, and 8 d hypoxia exposure compared to the control groups (Figure 9).
As shown in Figure 9, the expression of Fibropellin-1 increased 248.9-fold higher levels in 2 d post hypoxia exposure compared to normoxia control group. HSP70 are significantly up-regulated in R. philippinarum under hypoxia stress with 6.6-fold up-regulation at 8 d compared with control group. The expression value of IAP increased significantly from 7.3 at 2 d to 50.5 at 8 d under hypoxia stress. The expression of GPx and PEX increased 4.66 and 6.12-fold higher levels in 8 d post hypoxia exposure compared to normoxia control group. On the other hand, the expression of Cadherin and Calmodulin was decreased to 0.59 and 0.66-fold lower levels at 2 d compared with control group. The expression value of E3 ubiquitin-protein ligase and Ubiquitin-conjugating enzyme E2 decreased significantly from 0.32 and 0.54 at 2 d to 0.17 and 0.14 at 5 d under hypoxia stress, respectively. The expression value of Defensin was significantly decreased over 10-fold down-regulation at 2 d hypoxia treatment compared with control group. Overall, the qPCR validation results of 7 up-regulated genes and 7 down-regulated genes were consistent with the RNA-seq analysis.