DOI: https://doi.org/10.21203/rs.2.17553/v1
Background: Hypoxia is an important environmental stressor in aquatic ecosystems, with increasingly impacts on global biodiversity. Benthic communities are the most sensitive parts of the coastal ecosystem to eutrophication and resulting hypoxia. As a filter-feeding organism living in the seafloor sediment, Ruditapes philippinarum represents an excellent "sentinel" species to assess the quality of marine environment. In order to gain insight into the molecular response and acclimatization mechanisms to hypoxia stress in marine invertebrates, we examined hypoxia-induced changes in immune-related gene expression and gene pathways involved in hypoxia regulation of R. philippinarum.
Results: We investigated the response of the Manila clam R. philippinarum to hypoxia under experimental conditions and focused on the analysis of the differential expression patterns of specific genes associated with hypoxia response by RNA-seq and time course qPCR analysis. A total of 75 genes were captured significantly differentially expressed, and were categorized into antioxidant/oxidative stress response, chaperones/heat shock proteins, immune alteration, and cell proliferation/apoptosis. Fourteen hypoxia responsive genes were validated significantly up/down regulated at different time 0, 2, 5, and 8 d in gills of R. philippinarum in hypoxia challenged group. Functional enrichment analysis revealed the HIF signaling pathway and NF-κB signaling pathway play pivotal roles in hypoxia tolerance and resistance in R. philippinarum.
Conclusion: The HIF signaling pathway and NF-κB signaling pathway play a critical role in hypoxia tolerance and resistance in Manila clam. The immune and defense related genes and pathways obtained here gain a fundamental understanding of the hypoxia stress in marine bivalves and provide important insights into the physiological acclimation, immune response and defense activity under hypoxia challenge. The reduced metabolism is a consequence of counterbalancing investments in immune defense against other physiological processes.
Hypoxia is considered to be one of the most important stressors to aquatic animals [1]. In recent years, hypoxia in marine habitat has draw much more attention due to it is one of the ecological concerns in the world, especially the discharge of nutrient rich water and excessive anthropogenic input of organic matter into the sea [2]. Hypoxic layer often occurs near the sea floor during the summer in coastal areas, and leads to the mass mortality of the benthic animals [3,4]. Notably, excessive nutrient enrichment and human activities may cause water eutrophication and result in hypoxic or anoxic bottom environments [5,6]. Hypoxia may have remarkably impacts on aquatic animals at different levels. It has been reported that the behavioural, physiological and biochemical responses to hypoxia stress in several marine mollusca species [1,7,8].
Marine intertidal animals frequently encounter hypoxia or anoxia during low tide aerial exposure. Most sediment dwelling organisms are well adapted to the hypoxic aquatic realm and marine bivalves live in the intertidal region have to survive in the fluctuating environment at regular intervals. Indeed, intertidal bivalves are not passive to environmental changes and exhibit dynamic changes in their physiology, gene expression and metabolism to meet changes in their environment. When bivalves are challenged by hypoxia, a great number of physiological and molecular reactions are involved to deal with the low oxygen level stress [9]. One response of bivalves to aerial exposure or very low levels of hypoxia is to close their shells and regulate their internal environment. The respiratory activities usually slow down to adjust water flow and oxygen uptake, and reduce oxygen consumption with decreased energy expenditure [10]. In most marine invertebrates, the physiological adaptation to anaerobiosis caused by hypoxia or anoxia include maintaining reserves of energy store such as sugars or aspartate, and producing alternative end products to enhance ATP production [11]. Meanwhile, hypoxia may cause a slowdown of the process of reactive oxygen species (ROS) generation, and thus a decline in antioxidant enzymes activities [12], which could affect the immune and survival of shellfish species. So far, a number of studies have been conducted on the effects of hypoxia or anoxia on metabolic rate and modulation of enzyme activity in marine bivalves [13–15], as well as several research works in some shellfish species on molecular level [8,16,17].
Alkaline phosphatase (AKP) exerts a protective role in physiological function and immune defense such as detoxification of LPS and protection against endotoxin, an equally ubiquitous product of Gram-negative bacteria that may cause lethal complications after an infection with microorganisms [18,19]. Succinate dehydrogenase (SDH) activity represents the physiological conditions of substrate oxidation and energy metabolism, along with mitochondrial function [20]. Lactate dehydrogenase (LDH) plays a critical role in maintaining aerobic metabolism by converting lactate, the major by-product of anaerobic glycolysis, to pyruvate via oxidation in the presence of its coenzyme nicotinamide adenine dinucleotide, (NADH) [21]. Therefore, these enzymes AKP, SDH, and LDH can be used as biochemical and physiological markers (i.e. biomarkers) for the assessment of aquatic animal health conditions. Over the past decade, RNA sequencing (RNA-seq) has become an indispensable tool for transcriptome-wide analysis of differential gene expression and elucidation of the mechanisms underlying suites of ecological processes [22]. The differentially expressed genes with oxygen depletion stress and transcriptional responses to hypoxia were reported in the oyster Crassostrea gigas [8], mussel Mytilus galloprovincialis [16], and abalone Haliotis diversicolor [17].
The Manila clam, Ruditapes philippinarum, is an economically and scientifically important marine bivalve species, which has a wide geographic distribution from Europe to Asia [23]. Manila clam is one of the major aquaculture species in the world and the production of R. philippinarum reached over 4.0 million tons, equivalent to 3.7 million USD, in 2015 [23]. In recent years, however, this species faces greater risk of exposure to hypoxia as eutrophication worsens throughout its coastal habitats. As a filter-feeding organism living in the seafloor sediment, R. philippinarum represents an excellent “sentinel” species to assess the quality of marine environment [24]. In the past decades, hypoxia has become a serious issue in coastal environmental conservation, since some of the clams die because of the hypoxic water [25,26]. Hence, the tolerance of the bivalves for hypoxia and/or anoxic conditions has been investigated by a number of researchers [1,25,27]. More recently, the effects of hypoxia on survival, behavior, metabolism and cellular and tissue damage of R. philippinarum were conducted, suggesting that severe hypoxia significantly affects the physiology of R. philippinarum [28,29]. Until now, however, little effort has been put towards elucidating the molecular mechanisms combined with physiological and biochemical response of clams to the hypoxia resistance and acclimations. The underlying mechanisms of molecular alterations and responses to hypoxia stress in R. philippinarum remains largely unknown.
In this study, RNA-seq approach was adopted to investigate the transcriptome profiles of the gills from R. philippinarum under hypoxia and normoxia conditions. This work aimed to identify the differentially expressed genes and their expression patterns under low oxygen challenge to better understand the transcriptomic regulation in response to hypoxia stress and to investigate hypoxia-induced changes in immune-related gene expression and gene pathways involved in hypoxia regulation of R. philippinarum. Meanwhile, the responses of the energy metabolism and immune related enzymes (SDH, LDH and AKP) were investigated in R. philippinarum under hypoxia challenge at experimental conditions. These results provide new insights into understanding of the integrative regulation mechanism at the physiological, biochemical and the molecular level involved in hypoxia tolerance and resistance in R. philippinarum.
We produced six RNA-seq libraries 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 Manila clam and to investigate the effects of hypoxia on their molecular pathways and stress responsive genes involved in this processes. 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). Trinity software was applied to combine and assemble the clean reads from six libraries into a reference transcriptome. 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 total, 106,689 unigenes were annotated in the Nr and Swiss-Prot databases, with an E value of 1.0e–5 of BLASTX. A total of 61,972 unigenes (21.33%) were annotated in Nr, 68,093 unigenes (23.44%) were annotated in PFAM database, and 44,717 (15.39% of the total) were annotated in Swiss-Prot. All the unigenes were enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms (Table 2). 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).
To determine whether gene transcripts with related functions were significantly enriched in the hypoxia regulation, we analyzed the enrichment of the GO biological processes, molecular functions, and cellular components in R. philippinarum under hypoxia challenge. GO analysis revealed 68,222 unigenes annotated in one or more GO term. Genes related to vital biological processes, such as metabolic process, cellular process, biological regulation, localization, response to stimulus, and cell surface receptor signaling pathway were identified (Figure 3). Differentially expressed genes involved in peroxisome, mitochondrial envelope, and transferring phosphorus-containing groups were identified in the cellular components (Additional file 2). 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 5). In addition, genes involved in signaling receptor activity, calcium ion binding, receptor signaling protein activity, and structural molecule activity were identified (Additional file 2).
In this study, several hypoxia regulated functions such as immune system process, antioxidant activity, immune response, tumor necrosis factor receptor binding, peroxidase activity, regulation of apoptotic process, and oxidoreductase activity were shown in Table 3. 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 Manila clam might acclimated to hypoxia through immune and defense response to stress. Decrease in defense response and response to stress activity was related to regulation of response to oxidative stress, cell surface receptor signaling pathway, and reactive oxygen species (Additional file 2). Additionally, a number of unigenes 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 2), which may play potential roles in the hypoxia tolerance of R. philippinarum.
Eukaryotic Orthologous Groups (KOG) analysis was performed to classify all unigenes into different functional categories (Figure 4). The signal transduction mechanisms contain the largest number of unigenes 3805, the chaperones, protein turnover, and posttranslational modification containing 2502 unigenes. The categories of interest include transcription (1214), RNA processing and modification (1004), intracellular trafficking, secretion, and vesicular transport (1424), and translation, ribosomal structure and biogenesis (1137). 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 pairwise comparisons 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).
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. Gene expression analysis showed that a total of 75 unigenes were differentially expressed genes (DEGs) after exposed to hypoxia stress (qvalue<0.005 and |log2(foldchange)|>1). The DEGs include 43 up-regulate and 32 down-regulate unigenes expressed significantly different between hypoxia stress group and normoxia control group (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. The numbers of specific genes and shared genes between Hypoxia group (H1, H2, and H3) and Normoxia group (N1, N2, and N3) was shown as venn diagram (Figure 5C).
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 3). 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.
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.
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 of 14 genes were analyzed in two tissues at different times of hypoxia stress. As shown in Figure 9, the expressions of 14 hypoxia responsive genes were significantly different between hypoxia group and control group at different time 2, 5, and 8 d under hypoxia stress. These results showed a high consistency with the FPKM analysis of RNA-seq. 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). We found Fibropellin–1 to be expressed at more than 200-fold higher levels in 2 d post hypoxia exposure compared to normoxia control (Figure 9). The qPCR results of 7 upregulated genes and 7 downregulated genes were consistent with the RNA-seq analysis.
Hypoxia activates a variety of complex pathways at both the cellular and organism level, with the ultimate aim of reinstating oxygen homoeostasis. In the past decades, physiological and biochemical responses to hypoxia have been studied in several marine bivalve species [26,28,39,40]. Adenosine–5-triphosphate (ATP), the major energy source for cells in the body, is predominantly supplied by a series of metabolic pathways including glycolysis, the citric acid cycle, and the electron transport chain [41]. Among these, glycolysis occurs (with variations) in nearly all organisms, both aerobic and anaerobic. The wide occurrence of glycolysis indicates that it is one of the most important energy metabolic pathways [42]. It provides not only high-energy compounds like ATP, but also pyruvate, which can be used in the citric acid cycle to generate more ATP, NADH, and FADH2 [43–45]. In this study, the enzymatic assays of LDH indicated that the hepatopancreas of R. philippinarum switch to anaerobic metabolism at 8 d under hypoxia at 15 °C, while the SDH activity is likely to related to mitochondrial aerobic energy production, suggesting that aerobic and anaerobic metabolism were both existed in R. philippinarum. When Manila clams are under hypoxia, their shell valves are frequently closed, likely undergoing low oxygen consumption and reduced metabolisms [46]. The activity of LDH is closely related to cell metabolism, and its activity is often used as an indicator to evaluate the level of anaerobic metabolism [47]. As shown in Figure 8A&B, the level of dissolved oxygen in the water is insufficient with the hypoxia exposure time increased, and the clams may depend on anaerobic respiration to supply the metabolic energy, which is strongly reducing their energy requirement, so the activity of SDH and LDH decreases with prolonged hypoxia challenge. The different changes in the activities of key enzymes involved in glycolysis may indicate a diverse strategy of shellfish under hypoxia [29].
The hypoxia condition is a challenging situation because the energy production is drastically decreased in anaerobic metabolic pathways compared to aerobic energy provision. The decrease in metabolic rate at hypoxic environment is mainly caused by two mechanisms, including the depression of ATP-requiring processes and the inhibition of ATP-generating pathway. The reduction of the metabolic rate under hypoxia challenge provides benefits for facultative anaerobes. The glycogen stores of Mytilus edulis would only be enough for 3 days to provide ATP by anaerobic fermentation. However, mussels are able to sustain their life at anoxic conditions for weeks by reducing energy requirement [48]. The Manila clam exhibited strong tolerance to hypoxia as the 20-day LC50 for dissolved oxygen (DO) was estimated to be 0.57 mg L–1, and the LT50 at 0.5 mg L–1 DO was 422 hours.29 In this study, AKP activity in hepatopancreas was significantly increased at 5 d in response to hypoxia stress, indicating that the hepatopancreas is an important organ that involved in the immune response and hypoxia regulation in R. philippinarum. In the gills, however, no significant change in AKP activity observed between hypoxia challenge and control groups during 8 days, suggesting a tissue-specific response exist in R. philippinarum, which was also reported in C. gigas [1].
Previous studies explored how environmental changes are transduced into coordinated changes in physiological processes by examining gene expression cycles using microarray in intertidal bivalves [49–54]. More recently, transciptomic profiling is widely employed in elucidating the genetic basis of molecular response to different environmental factors in marine bivalve [55,56]. In this study, a transcriptomic approach was firstly employed to investigate the gene expression patterns of the Manila clam under hypoxic challenge. Hypoxia responsive genes associated with stress response, substance and energy metabolism, and immune response were corroborated by qPCR. In this study, all of the 14 DEGs detected by qPCR were consistent with the RNA-seq results. The expression profiles of these genes were characterized at different hypoxia exposure duration (0, 2, 5, and 8 d). These genes reflect different aspects of the stress adaptation and molecular mechanisms of immune response to hypoxia exposure. The KEGG analyses of DEGs show that molecular pathways related to TNF signaling, NOD-like receptor signaling, RIG-I-like receptor signaling, TGF-beta signaling AMPK signaling, NF-kappa B signaling, cAMP signaling pathway, and apoptosis are highly enriched, implicating these processes are important for physiological adaptation and immune response to hypoxia stress in R. philippinarum. Manila clam may have evolved various sophisticated signaling pathways to sense their immediate environment and orchestrate appropriate transcriptional responses that mediate adaptation in benthic environment.
Aerobic organisms have the antioxidant mechanisms to prevent the oxidative damage and protect them against oxidant stress. This biological process involve in a number of cellular molecules and antioxidant enzymes [51]. In the last decade, however, most studies in marine bivalves investigated the activities of the antioxidant enzymes, but only a few studies reported both the biochemical responses and molecular actors governing the antioxidant following hypoxia [57,58]. In this study, we observed a rapid regulation of a set of genes involved in the molecular responses to oxidative stress. The cytochrome P450, peroxisomal membrane protein (PEX), and Glutathione peroxidase (GPx) were both upregulated under hypoxia challenge. Peroxidase, an important class of antioxidation enzymes, is responsible for the defense response to oxidative stress and degradation of ROS [59]. The expression of GPx mRNA increased under hypoxia challenge to protect cells from ROS that can be formed upon reoxygenation, which has been reported in C. gigas [1,8]. Stress-induced immune changes have been reported in many marine invertebrates, including oyster [60,61] and mussel [62,63]. The innate immunity is the mainly immunological defense mechanism in invertebrate metazoan [64]. In the present study, several immune-related genes such as the defensin, IAP and HSP70 were up-regulated under the hypoxia challenge, demonstrating that these genes are involved in defense and immune response. In addition, focal adhesion, NF-kappa B signaling pathway, apoptosis, TNF signaling pathway were enriched by KEGG analysis, suggesting that these pathways play significant roles in the immune response and defense mechanisms against hypoxia stress.
Hypoxia inducible factors (HIFs) are a family of highly conserved transcription factors that act as main regulators of oxygen homeostasis and the adaptive response to hypoxia [65]. NF-κB signaling pathway has been proved to be involved in innate immune response to bacterial infection and hypoxia stress in molluscan shellfish, including R.philippinarum [66], Meretrix meretrix [67], C. gigas [68], and H. diversicolor [17]. It is indicated that genes associated with NF-κB signaling pathway participated in the immunomodulation process to respond to hypoxia stress [17]. Therefore, we focus on the HIF–1 and NF-kappaB signaling pathways in the present study. Nuclear factor kappa B (NF-κB) is a transcription factor regulates the expression of cytokines, effector enzymes and apoptosis inhibitors in response to extracellular signals. NF-κB signaling plays a critical role in immune defense, stress responses and inflammation. NF-κB comprises RelA (p65), RelB, c-Rel, NF-κB1 (p50) and NF-κB2 (p52). NF-κB can stimulate transcription of its target genes in a very quick fashion, as it exists freely in the cytoplasm, albeit inhibited by IκB proteins. In the current study, NF-kB is upregulated, which could increase the HIF–1α mRNA level, as well as the PI3K and p70S6K are upregulated in PI3K-Akt signaling pathway.
Enolase is a hypoxic stress protein, which may contribute to hypoxic tolerance by increasing anaerobic metabolism [69]. As one the downstream genes of the HIF signaling pathway, the enolase (ENO1) is significantly upregulated in the hypoxia challenged clams, which promote anaerobic metabolism and causing reduced oxygen consumption. Concomitantly, hypoxia induced NFκB activity is preceded by temporally sequential IKK activation, IκB phosphorylation, and IκB degradation, indicating that hypoxia activates NFκB through increased IKK activity (Figure 7). The IKK-dependence of hypoxia-induced NFκB activity is consistent with most other stimuli of this pathway [70]. The role of IKK in oncogenesis and inflammation is well established, but there is emerging evidence that IKK may have antiinflammatory functions [71,72]. In this study, we also observed an increase in IKKβ in hypoxia (Figure 7). The observation that IKKβ are hypoxia-sensitive, suggests the possibility of a molecular mechanism where a p50- mediated response could be resolved over a course of hypoxia [70]. Therefore, NF-κB plays important roles in regulating genes responsible for both the immune response and hypoxia stress resistance.
The physiological, biochemical and molecular responses of R. philippinarum under hypoxia stress were firstly investigated in the present work. The results provide comprehensive understanding of the molecular mechanism for hypoxic tolerance and resistance in marine bivalve species. The majority of immune related genes and pathways play vital roles in the immune defense against hypoxia stress. The reduced metabolism might be the consequence of counterbalancing investments in immune defense against other physiological processes. We speculate that it is likely to be an adaptive strategy for Molluscs to survive. This study provides useful information for further study of the critical molecular molecules regulating the hypoxia tolerance and gain new insights on mechanisms of resistance to hypoxia stress in marine bivalves.
Manila clams were collected from Zhuanghe, Dalian, China. After being transported from the field to the laboratory, the clams were cleaned to remove any fouling and were acclimated for one week in aerated 70 L plastic tanks, containing water at 15 ± 0.5 °C with salinity of 32 ppt. Other water quality parameters were measured during the experiment (pH: 8.0 ± 0.2, dissolved oxygen: 7.0 ± 0.5 mg/L). All the clams were fed with Spirulina powder daily for one week before the hypoxia challenge, and water was exchanged once per day to discharge waste products from marine invertebrates. Manila clam is not an endangered or protected species, so no specific permits were required for this study.
For the hypoxia challenge, clams were placed in another closed tank in which O2 concentration was controlled by monitor through nitrogen injection. A feedback was provided to the monitor by an O2 sensing probe (YSI ProPlus, USA) placed in the tank and nitrogen flow was adjusted in order to maintain O2 level to a set-value (2.0 mg/L). A total of 120 clams with similar size were divided into six groups in separate tanks. Hypoxia challenge group (H1, H2 and H3) was kept at low oxygen level (2.0 mg/L), and normoxia group (N1, N2 and N3) was at normal oxygen level (7.0 mg/L). The temperature and salinity of seawater were at 15 °C and 32 ppt, respectively. After 24 h, the gills from three Manila clams in each hypoxia challenge groups and control groups were sampled, respectively, and frozen in liquid nitrogen and stored at –80 °C until RNA extraction. Meanwhile, the gills and digestive gland of Manila clams were sampled at 0, 2, 5, 8 d with three clams per tanks for each of the hypoxia challenge groups and control groups (three biological replicates with three technical replicates).
The total RNA was extracted from 30 mg gill tissue sample from each individual using RNAprep pure Tissue Kit (TianGene, Beijing, China), according to the manufacturer’s protocol. RNA purity was determined with a NanoPhotometer® spectrophotometer (Implen, CA, USA). The RNA concentration was measured with the Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed with the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
About 100 mg sample of gills and digestive glands were taken from each randomly chosen individual from each tank in the hypoxia challenge groups and control groups. Glass homogenizers were used to homogenize this tissue on ice for 2 min. The frozen samples were treated with 5 volumes of 0.05 M PBS buffer (pH 7.5). The homogenates were centrifuged for 8 min at 3500 × g, and the supernatant was collected and assayed for enzyme activities. All samples were tested using analysis kits manufactured by Nanjing Jiancheng Bioengineering Institute (Nanjing, China). The activities of LDH (EC 1.1.1.28), SDH (EC 1.3.5.1), and AKP (EC 3.1.3.1) were measured according to the methods [30–32].
Six cDNA libraries were constructed from R. philippinarum in the hypoxia challenged group (H1, H2 and H3) and control group (N1, N2 and N3). About 3 μg of RNA of each sample was used for libraries preparation and Illumina paired-end sequencing. The NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) was used for sequencing libraries according to the manufacturer’s instruction. The AMPure XP system (Beckman Coulter, Beverly, USA) was used for purification of the library and 150–200 bp cDNA fragments were select. The adaptor-ligated cDNA was incubated with 3 μl of USER Enzyme (NEB, USA) at 37°C for 15 min, and then at 95 °C for 5 min. The Phusion High-Fidelity DNA Polymerase was used for PCR with universal PCR primers and Index (X) Primer. The AMPure XP system was used for the purification of PCR products and the library quality was evaluated with the 2100 Agilent Bioanalyzer.
Trinity software (version 2.3.2) was used for de novo transcriptome assembly without reference genome. The sequence reads have been cleaned by removing adaptor sequences, ambiguous ‘N’ nucleotides (with the ratio of ‘N’ to be more than 10%), and low quality sequences (with quality score less than 5) [33]. For annotation, non-redundant sequences were subjected to public databases, including the non-redundant protein (Nr) and non-redundant nucleotide (Nt) databases at NCBI (http://www.ncbi.nlm.nih.gov/), Swiss-Prot (http://www.ebi.ac.uk/uniprot/), Gene Ontology (GO) (http://www.geneontology.org/), Clusters of Orthologous Groups (COG) (http://www.ncbi.nlm.nih.gov/COG/), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/). The sequences were matched with an E-value cut-off of 10–10 in Nr, Nt, and Swiss-Prot using BLASTX. The KEGG classifications and COG were conducted with BLASTX with an E-value cutoff of 10–5.
Differentially expressed genes (DEGs) were measured by counting tags from hypoxia/thermal treated samples against the control, which are normalized using the RNA Sequence Expected Maximization (RSEM) method [34]. Initially, reads from control sample were mapped to reference transcriptome and subjected to check the differential expression using trinity utility scripts (align_and_estimate_abundance.pland abundance_estimates_to_matrix.pl) as instructed (http://trinityrnaseq.github.io/). The edgeR program package was applied and adjusted the read counts for each sequenced library using before the differential gene expression analysis [35]. To elucidate the differences between the hypoxia exposure and normoxiatranscriptomes, DEGseq was used to screen the differentially expressed genes in three replicate H1, H2 and H3 as the hypoxia challenged group transcriptome library and three replicate N1, N2 and N3 as the control group transcriptome library [36]. P values <0.05 and |log2 (fold change)|>1 was set as the threshold for significantly DEGs.
To corroborate the RNA-seq results, we select both upregulated and downregulated hypoxia-responsive differentially expressed genes to perform the quantitative RT–PCR (RT–qPCR) analysis. The expressions of fourteen hypoxia responsive genes in gills of R. philippinarum were conducted at different time 0, 2, 5, and 8 d in hypoxia challenge group (DO: 2.0 mg/L) and control group (DO: 7.0 mg/L), respectively (three biological replicate for each gene and time point). The Primer 5.0 software was used to design primers (Premier Biosoft International, Inc.) (Additional file 1). RNA (about 500 ng) samples were measured and treated with RQ1 RNase-Free DNase (Promega) to remove genomic DNA. The β-actin was used as the internal control for the qPCR analysis. The cDNA was synthesized using a reverse transcriptase reagent kit (PrimeScript™ RT reagent Kit, Takara). The qPCR was performed with SYBR Premix Ex Taq II (Takara). The reactions were carried out in a total volume of 25 μl containing 2.5 μl of diluted cDNA, 2.5 μl of each primer, and 12.5 μl of SYBR Green PCR Master Mix, with the following cycling profile: 95 °C for 15 min for polymerase activation, followed by 40 cycles at 95 °C for 15s, at 55 °C for 30s, and at 70 °C for 30s. Each sample was processed in triplicate in the Roche LightCycler 480 Real-Time PCR System (Roche). All data were analyzed using the 2 –△△Ct method [37].
Significant differences among sample data were analyzed using SPSS 20.0 for Windows (IBM SPSS Inc., Chicago, IL, USA) Software package. A one-way analysis of variance (ANOVA) was used to compare the effect of hypoxia on gene expression and three enzymes activities. A one-way ANOVA was used to compare the effects of hypoxia on gene expression and three enzymes activities including AKP, SDH, and LDH activities, followed by Duncan’s multiple comparison tests. Differences were considered to be significant at P < 0.05.
[Additional files are not available with this version.]
DEGs: Differentially expressed genes; ROS: Reactive oxygen species; GO: Gene ontology; KEGG: Kyoto encyclopedia of eenes and eenomes pathways; KOG: Eukaryotic orthologous groups; NF-κB: Nuclear factor kappa B; SDH: succinate dehydrogenase; LDH: lactate dehydrogenase; AKP: alkaline phosphatase; qPCR: quantitative real-time PCR; PEX: Peroxisomal membrane protein; SCP: Sterol carrier protein; GPx: Glutathione peroxidase; IAP: Inhibitor of apoptosis protein; ATP: Adenosine–5-triphosphate; COG: Clusters of Orthologous Groups; ANOVA: one-way analysis of variance
This work was supported by the National Key R&D Program of China (2018YFD0901400), the Modern Agro-industry Technology Research System (CARS–49), the Scientific Research project of Liaoning Education Department (QL201703), and Dalian high level talent innovation support program (Dalian Youth Science and Technology Star Project Support Program) (2016RQ065).
We would like to give thanks to Mr. Changwei Lu and Mr. Chenlin Chai for their help to prepare the samples and animals.
The raw sequences for R. philippinarum have been deposited in the NCBI PRJNA478917 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA478917). The RNA-seq datasets are available in the NCBI Sequence Read Archive (SRA) with the accession numbers SRP151767, SRS3494254, SRX4332302, SRR7462276, SRS3494255, SRX4332301 and SRR7462277.
HTN conducted preparation of the transciptome data, bioinformatics analysis and manuscript preparation. HMW and KYJ were involved in one or more processes of sample collection, RNA extraction or experimental and data analysis. HTN and XWY conceived, designed and guided the research. All authors read and approved the final manuscript.
Not Applicable.
Not Applicable.
The authors declare that they have no competing interests.
1College of Fisheries and Life Science, Dalian Ocean University, 116023 Dalian, China. 2Engineering Research Center of Shellfish Culture and Breeding in Liaoning Province, Dalian Ocean University, 116023 Dalian, China.
Table 1
Summary statistics of Ruditapes philippinarum transcriptome assembly using Trinity software.
Unigene |
Min. Length (bp) |
Mean Length (bp) |
Max. Length (bp) |
Median Length (bp) |
N50 |
N90 |
Total Length (bp) |
Number |
201 |
935 |
46,346 |
587 |
1389 |
408 |
271,479,222 |
Table 2 Summary statistics of functional annotation of R. philippinarum transcriptome.
Category |
All sequences |
300-999bp |
>=1000bp |
Total number of unigene |
148593 |
64105 |
26610 |
Unigene matches against Nr and Swiss-Prot |
20955 |
5787 |
12873 |
Unigene matches against GO |
32370 |
11094 |
15442 |
Unigene matches against KEGG |
10434 |
2985 |
6253 |
Table 3 GO classification of the differentially expressed genes from R. philippinarum under hypoxia stress.
Term |
Comparison |
Type |
Gene |
p-Value |
immune system process (GO:0002376) |
H1 vs N1 |
BP |
71 |
0.000218 |
antioxidant activity (GO:0016209) |
H1 vs N1 |
MF |
20 |
0.000309 |
immune response (GO:0006955) |
H1 vs N1 |
BP |
60 |
0.000311 |
tumor necrosis factor receptor binding (GO:0005164) |
H1 vs N1 |
MF |
20 |
0.000474 |
peroxidase activity (GO:0004601) |
H1 vs N1 |
MF |
7 |
0.003901 |
regulation of apoptotic process (GO:0042981) |
H1 vs N1 |
BP |
58 |
0.004872 |
oxidoreductase activity (GO:0016491) |
H2 vs N2 |
MF |
278 |
3.24E-07 |
antioxidant activity (GO:0016209) |
H2 vs N2 |
MF |
23 |
2.05E-05 |
peroxiredoxin activity (GO:0051920) |
H2 vs N2 |
MF |
7 |
0.000251 |
immune system process (GO:0002376) |
H3 vs N3 |
MF |
61 |
0.000900 |
antioxidant activity (GO:0016209) |
H3 vs N3 |
MF |
19 |
0.000159 |
immune response (GO:0006955) |
H3 vs N3 |
BP |
51 |
0.001554 |
tumor necrosis factor receptor binding (GO:0005164) |
H3 vs N3 |
MF |
21 |
2.92E-05 |
cytokine receptor binding (GO:0005126) |
H3 vs N3 |
MF |
34 |
4.99E-05 |
defense response (GO:0006952) |
H vs N |
BP |
5 |
0.000556 |
response to stress (GO:0006950) |
H vs N |
BP |
7 |
0.007614 |
oxidoreductase activity (GO:0016491) |
H vs N |
MF |
8 |
0.005337 |
peroxiredoxin activity (GO:0051920) |
H vs N |
MF |
1 |
0.017815 |
Table 4 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapping for R. philippinarum.
KEGG Pathways |
Sub-Pathways |
Number of Unigenes |
Cellular Processes |
Cell growth and death |
966 |
|
Cell motility |
451 |
|
Cellular commiunity |
1340 |
|
Transport and catabolism |
1802 |
|
Membrane transport |
118 |
|
Signal transduction |
3410 |
|
Signaling molecules and interaction |
680 |
Genetic Information Processing |
Folding, sorting and degradation Replication and repair |
1470 425 |
|
Transcription |
609 |
|
Translation |
1346 |
Metabolism |
Amino acid metabolism |
839 |
|
Biosynthesis of other secondary metabolites |
39 |
|
Carbohydrate metabolism |
960 |
|
Energy metabolism |
324 |
|
Glycan biosynthesis and metabolism |
598 |
|
Lipid metabolism |
964 |
|
Metabolism of cofactors and vitamins |
401 |
|
Metabolism of other amino acids |
393 |
|
Metabolism of terpenoids and polyketides |
71 |
|
Nucleotide metabolism |
596 |
|
Xenobiotics biodegradation and metabolism |
202 |
Organismal Systems |
Circulatory system |
621 |
|
Development |
614 |
|
Digestive system |
1014 |
|
Endocrine system |
1815 |
|
Environmental adaptation |
366 |
|
Excretory system |
375 |
|
Immune system |
1296 |
|
Nervous system |
1097 |
|
Sensory system |
420 |