Natural populations of adult pacific oysters were acclimated in natural seawater with salinity 24 psu and a temperature of 24°C for 2 days prior to experimental spawning. D-shaped larvae collected after strip spawning of 8 male and 23 female adult oysters, C. gigas were used for low pH exposure experiments (Breese and Malouf, 1975). Larvae of the C. gigas cultured from the early D-shaped developmental stage (i.e., 12 h post-fertilization) to the pediveliger stage (i.e., 15–17 day post-fertilization) under OA condition (Ko et al., 2014). For detailed information on the culture conditions, please refer to the “Materials and Methods” section of Ko et al. (2014) and Dineshram et al. (2016). The experiment was done in three replicates, using 40-liter tanks for the larval culture. Each larval culture tanks were mixed with air to obtain the ambient CO₂ levels for the control (pH 8.1) and CO₂ pre-mixed air for the low pH (7.4), which was within environmentally relevant range predicted for the Pacific oysters in the fluctuating coastal marine systems and reflects the futuristic atmospheric CO₂ emission scenarios (Zeebe, 2008), for the year 2300 as predicted by the Intergovernmental Panel on Climate Change IPCC (Caldeira and Wickett, 2005). The low pH 7.4 was chosen to reflect the futuristic CO₂ emissions as per the IPCC projected scenarios for the study region with respect to decreased pH to represent 2300 years considering far future acidification (Zhai et al., 2013). The CO₂ concentration in the experimental tanks was regulated through flowmeters (Cole-Parmer Inc., United States) and pH of the experimental system was constantly monitored thrice a day using a pH meter (Orion Star™, Thermo Electron Co., United States), respectively.

Regular laboratory oyster larval rearing procedures was followed throughout the experiment that includes filtered seawater (0.22μm), ambient temperature: 22–24^°C maintained using a continuous flow through system water bath. The carbonate chemistry of seawater used in the experiment was determined by total alkalinity (TA) measurements using Alkalinity Titrator (AC-A2, Apollo SciTech’s Inc., United States), calibrated against certified seawater reference material (Batch 103, A.G. Dickson, Scripps Institution of Oceanography). Using the measured temperature, salinity, pH and TA, other parameters of carbonate chemistry were calculated using co2sys.xls spreadsheet (Pelletier et al., 2007) with the equilibrium constants K₁, K₂ by Mehrbach (Mehrbach et al., 1973) and KSO₄ constants by Dickson (Dickson and Millero, 1987). The carbonate chemistry conditions for the pH levels were described in Table 1. Samples from pediveliger larvae with a size of ∼280 μm were collected after 18 days and stored at –80°C prior to quantitative transcriptome analysis. This study aimed to provide a set of differential uni-genes for the oyster pediveliger larvae exposed to seawater at elevated CO₂ levels, low pH 7.4 using next generation high-throughput transcriptome sequencing.

Total cellular RNA was extracted from larval samples (each group in triplicates) exposed to pH 8.1 and 7.4 using RNeasy mini kit (Qiagen). Briefly, the experiment started with four biologically independent replicates (replicated culture tanks), however, we have used only three of these replicates for transcriptomic analysis. This pooling of samples was done with purpose, i.e., to reduce biological variability between replicate samples obtained from multiple parents and to deduct the treatment effect without baseline variation (Dineshram et al., 2015, 2016; Griffith et al., 2019). RNA concentration and integrity were measured by NanoDrop (Thermo Fisher Scientific) and Agilent 2100 BioAnalyzer (Agilent Technologies), respectively. After treatment with DNase I, mRNA isolated with Oligotex mRNA Mini Kit (Qiagen) was quantified using NanoDrop 2000 spectrophotometer (Thermo Scientific). The mRNA pooled from triplicate samples in equal quantities was ready for cDNA synthesis.

cDNA libraries were constructed following TruSeq™ RNA sample preparation guide (Illumina). Briefly, the mRNA was fragmented to sizes of 120–200 bp. Then, first and second strand were synthesized using SuperScript II Reverse Transcriptase (Invitrogen) and Second Strand Master Mix, respectively. The double strand (ds) cDNA were isolated using Ampure XP beads. The adapter containing the A-Tailing oligonucleotide was ligated to dsDNA, and then amplified the library using 12 cycles of PCR. After purification, libraries were quantified by Qubit^® 2.0 Fluorometer and validated by Agilent 2100 bioanalyzer to confirm the insert size and calculate the mole concentration. Clusters were generated by cBot with the library dilution of 10 pM, and finally were sequenced on the Illumina Genome Analyzer (Hiseq 2500).

Quality of raw data generated from Illumina sequencing was checked by FastQC.¹ The raw reads were preprocessed using the FASTX² to remove the adapter sequences, low-quality reads, and reads less than 20 bp in length. Then, the clean sequence reads obtained from preprocessing were aligned to the C. gigas genome³ using Tophat 2.0 (Trapnell et al., 2009), allowing maximum 2 mismatch per read with genome sequence. Final results were stored in BAM format files.

The mapped reads were counted and then normalized with transcript length using cufflink (Trapnell et al., 2010), resulting in FPKM (Fragments Per Kilo bases per Million reads) for each transcript from two libraries. Thus, FPKM value facilitated to calculate the expression abundance of each gene that was on behalf of one transcript within two libraries. The differentially expressed genes (DEGs) were identified using DEGseq with each pair-wise comparison (Wang et al., 2010). Furthermore, the DEGs were filtered using the threshold of false discovery rate (FDR ≤ 0.1) and the absolute value of log₂ (FC) (≥ 1).

The DEGs were homology-searched using BLASTX against the Swiss-Prot database with an E-value of 1e^–5. Gene names were assigned to each sequence based on the best BLAST hit. To annotate sequences, the Gene Ontology (GO) terms were assigned to query sequences by BLAST2GO (Götz et al., 2008), producing group of genes classified in the three ontology categories, cellular component, biological process and molecular function. Furthermore, pathways were analyzed using the online KEGG Automatic Annotation server (Kyoto encyclopedia of genes and genomes⁴) with bi-directional best-hit setting. Finally, functional networks of enriched KEGG pathways were visualized using ClueGO.⁵

Expression of five up-regulated genes such as ADP-ribosylation factor 1, Hemicentin-1, Iron (II)-dependent oxidoreductase EgtB, NADH dehydrogenase and Glutathione S-transferase omega-1 and three down-regulated genes such as Zinc finger protein, FMRF- amide neuropeptides and Calcium/Calmodulin-dependent protein kinase type I revealed as differentially expressed due to low pH exposure by RNA sequencing was determined by RT-qPCR (Table 2). PCR amplifications was performed in triplicate biological samples on a Roche 454 sequencer (United States) in a final volume of 20 μl containing 1μl of template RNA mixed with 10 μl of 2X Master Mix (Roche) and 0.8 μl of each forward and reverse primer in each RT- PCR reaction for cDNA synthesis. Amplifications was carried out for 40 cycles on a 7500 Real Time PCR system (Applied Biosystems) under conditions of 30 s at 95°C for 1 cycle, 5 s at 95°C, 30 s at 60°C for 40 cycles. Thermal denaturing step to generate a dissociation curve to verify amplification efficiency and primer dimer was performed by increasing temperature from 65–99°C with a 5 s hold for every one degree temperature increment. RT-qPCR products were also further checked on 1.5% agarose gels to check the specificity of amplicons. The mRNA expression of each gene was quantified and calculated using their normalized expression level relative to the housekeeping gene by the 2^–ΔΔCt method (Livak and Schmittgen, 2001).

Two libraries were sequenced with Illumina (GEO accession numbers SRR960964 and SRR967068) and characteristics of transcriptome were summarized below. There were 16,602,870 and 14,795,142 raw reads generated from the control and one low pH library, respectively. After low quality reads removal and adaptor trimming, a total of 15,081,854 effective reads were produced from the control group, while 13,455,384 effective reads were obtained from the treated group, whose effective reads ratio were 90.8 and 90.9%, respectively. These effective reads were then mapped into C. gigas genome (Zhang et al., 2012), revealing 11,515,356 reads in control group and 9,746,368 reads in low pH group were able to be mapped. The mapping ratios for the control and low-pH libraries mapping ratios were 83.4 and 81.1%, respectively. Finally, these mapped reads were collected and used for quantification and annotation analyses.

To investigate the mechanisms underlying oyster larval tolerance to low pH, transcriptomic expression profiling was performed by comparing the control (ambient pH) library and low pH libraries using the FPKM method. Scatter plot showed a total of 1,808 DEGs obtained from the larval oysters exposed to low pH treatment, including 1,176 up-regulated and 632 down-regulated genes, respectively (Figure 1). Among these genes, 1410 DEGs were significantly matched against Swiss-Prot database by BLASTX, accounting for 78% of all DEGs (Supplementary Table 1). After functional annotation of these DEGs, all assigned GO terms were classified into three GO categories: biological process (1,138, 26.4%), cellular component (1,380, 32.1%) and molecular function (1,787, 41.5%) (Figure 2 and Supplementary Table 2). Within molecular function, structural constituents of ribosome (GO: 0003735) and calcium ion binding (GO: 0005509) were highly represented (p < 0.001), highlighting the effect of protein translation and shell formation at low pH. Regarding the cellular component, the most represented (p < 0.001) GO categories were nucleus (GO: 0005634) and extracellular region (GO: 0005576). As in the case of biological process, translation (GO: 0006412) was the most significantly represented (p < 0.001) followed by proteolysis (GO: 0006508), cell adhesion (GO: 0007155) and apoptotic process (GO: 0006915). Almost two-thirds of the differentially regulated genes identified in the current study were up-regulated. The proportion of down-regulated genes was lower, indicating that pH stress does not suppress the activities of major cellular pathways but rather predominantly induces cellular homeostatic mechanisms that up-regulate energy metabolism, cell signaling and ribosomal biogenesis. This gene regulation pattern confirms that low pH alters post-transcriptional regulatory mechanisms irrespective of gene expression.

To explore the possible biological pathways associated with low pH DEGs, this study mapped these DEGs into the KEGG database. Among them, 315 DEGs received KEGG annotations and were assigned to 27 pathways (Figure 3 and Supplementary Table 3). Most of the DEGs were found to be enriched in a few pathways, such as Ribosome (ko03010, 42 DEGs), MAPK signaling pathway (ko04010, 37 DEGs), Viral carcinogenesis (ko05203, 35 DEGs), Peroxisome (ko04146, 34 DEGs) and Oxidative phosphorylation (ko00190, 29 DEGs) with high significance (p < 0.001). Based on the intersection of biological functions and pathways identified in this transcriptomic study, ClueGO Cytoscape pathway analysis was performed to provide integrated and comprehensive insights (Bindea et al., 2009, 2013; Figure 4).

Strikingly, all 23 DEGs assigned to the Oxidative phosphorylation (ko00190) showed to be increased after low pH treatment (two fold change) in genes NADH dehydrogenase, isoforms of ATP synthase, cytochrome c oxidase and V-type ATPase involved in electron transport chain. Oxidative phosphorylation is a metabolic pathway to produce ATP and energy release by the oxidation of nutrients. Thus, oyster larvae likely increase their energy production by up-regulating electron transport chain genes in response to low pH exposure. These results were consistent with previous studies in which adult oysters exposed to elevated temperature exhibited increased mitochondrial respiration (Chamberlin, 2004; Ivanina et al., 2012). In contrast to previous studies that showed metabolic depression in response to low-pH stress in adult blue mussel (Hüning et al., 2013), sea urchins (Padilla-Gamiño et al., 2013) and corals (Todgham and Hofmann, 2009; Kaniewska et al., 2012, 2015; Moya et al., 2012; Rocker et al., 2015), the current study found that genes corresponding to the basal energy metabolic pathways were up-regulated by two fold. The metabolic pathways involved in energy production were significantly down-regulated in PDI stage of pacific oyster larvae (Liu et al., 2020). Also studies on adult pearl oyster Pinctada fucata and razor clams reported metabolic depression along with up-regulation of V-type ATPase and energy producing genes in response to ocean acidification (Li et al., 2016a; Peng et al., 2017) showing that the metabolic depression can also happen without the suppression of energy producing genes. However, other studies in mussels (Thomsen and Melzner, 2010; Fernández-Reiriz et al., 2012), sea-urchins (Stumpp et al., 2011), in arctic pteropods (Comeau et al., 2010) and Corals (Vidal-Dupiol et al., 2013; Davies et al., 2016; González-Pech et al., 2017; Rivest et al., 2018) showed either up-regulation or no change in the expression of energy metabolism genes in response to moderate changes in the ocean acidification. The expression pattern observed in the current study suggests that aerobic pathways are activated to generate more energy to compensate for low pH-induced energy depletion.

Immune responses are energy-intensive processes that require large amounts of metabolic energy which involves physiological trade-offs in energy balance to meet energy demands upon activation of the immune response (French et al., 2009). Several pathways involved in immune response were identified by the pathway enrichment analysis, including the Chemokine signaling pathway (ko04062), Natural Killer cell mediated cytotoxicity (ko04650) and B cell receptor signaling pathway (ko04662). Interestingly, all three immune pathways were dramatically inhibited in the low pH group, comprising of 90, 83.3, and 85.7% of the down-regulated DEGs, respectively. Immune related genes in the hemocytes were significantly down-regulated in adult Crassostrea gigas when exposed to ocean acidification (Wang et al., 2016) but it was found to be increased in other oyster species (Goncalves et al., 2016, 2017). Most of the genes corresponding to the immune pathways encoded kinases responsible for immune defense in many bivalves. Kinases are usually activated by stressors such as UV, pollutants, heat and osmotic stress in most bivalves (Canesi et al., 2006). The roles of kinases in defense mechanisms have been illustrated in mussel hemocytes; mussels exposed to a bacterial challenge showed induction of kinase-mediated signal transduction pathways (Canesi et al., 2006). The down-regulation of genes encoding kinases such as mitogen activated protein kinase 3 (–2.8-fold), cAMP-dependent protein kinase catalytic subunit (–3-fold) and glycogen synthase kinase-3 beta (–2.6-fold), in addition to the up-regulation of phosphatidylinositol 3-kinase regulatory subunit alpha (+3.1-fold), clearly suggests that seawater acidification could inhibit immune responses in pediveliger larvae through impaired signal transduction. Notably, down-regulation of fucolectin-1, an immune protein that binds to bacterial surfaces and contributes to defense, was also observed, consistent with the effect of OA on the innate immune system (Honda et al., 2000). Stimulation of TLR pathway and immune factors is known to occur in Crassostrea gigas when exposed to low pH (Cao et al., 2018). In humans, B-cell lymphoma/leukemia 11A (Bcl11a) functions as a key transcriptional factor that controls the normal immune cell development, and in the present study its homolog in pediveliger larvae was also suppressed upon exposure to low pH seawater conditions. Therefore, OA may impair immune cell development and immune signaling in oysters. However, this speculation requires further verification (Liu et al., 2003).

In total, 34 and 18 DEGs involved in formation of Peroxisome (ko04146) and Lysosome (ko04142) formation, respectively, were mapped and shown to be activated, and 32 and 17 of these DEGs were up-regulated at the mRNA level. Additionally, 42 DEGs were assigned to the ribosome pathway (ko03010), and all of the associated DEGs were up-regulated after exposure to low pH. Most of the genes involved in the peroxisomal pathway were oxidases involved in the production of reactive oxygen species (ROS) in cells (Cancio and Cajaraville, 1997). Similarly induction of oxidative stress and subsequent increase in ROS and apoptosis has been observed in the adult Crassostrea gigas (Wang et al., 2016; Cao et al., 2018) while there was no significant increase ROS production in adult Crassostrea gigas (Wang et al., 2020), Saccostrea glomerata (Goncalves et al., 2016, 2017) under long term exposure to low pH. The up-regulation of these genes in low-pH seawater suggests that oxidative stress is induced in mollusks, especially in oysters, by intracellular acidosis and inefficient electron transport (Tomanek et al., 2011). Increased oxidative damage caused by low pH also induced the expression of ubiquitin-related genes in this study, including Ubiquitin-like modifier-activating enzyme 1, Ubiquitin-conjugating enzyme E2, E3 ubiquitin-protein ligase SHPRH, UBR4, HRD1, and de-ubiquitinase BRCC36. On contrary, NADH dehydrogenase ubiquinone and oxidoreductase activity were down-regulated in Crassostrea gigas larvae prior to shell formation (Liu et al., 2020). Considering the crucial role of ubiquitin system in survival against many environmental stressors, including oxidation, activation of ubiquitin system may contribute to host protection from OA-derived stress (Hanna et al., 2007). Activation of ubiquitination protein in Crassostrea hongkongensis (Lim et al., 2021) and up-regulation of oxidoreductase activity was observed in Crassostrea virgininca (Downey-Wall et al., 2020). Lysosomal genes such as Cathepsin L1, F and B, also were found up-regulated (≥2fold) in this study. Studies using different forms of stressors including temperature (Zhang et al., 2006), pollutants (Sheehan and Power, 1999) and salinity (Stickle et al., 1985) have observed lysosomal membrane destabilization leading to the release of lysosomal enzymes into the cytoplasm (Moore et al., 1982). In this study, the oyster larvae displayed lysosomal destabilization in response to low pH stress. Increased activity of lysosomal enzymes plays a vital role in the mobilization of energy reserves to ensure normal metabolic function during periods of stress in bivalves (Tremblay et al., 1998). However, excessive lysosomal activity may cause drastic changes in the cellular activities leading to apoptosis (Zhao et al., 2003; Sokolova, 2009). Similarly, the highly active, autolytic, arylsulfatases, were found to be up-regulated by fourfold in this study. Spinal cord injury in mice causes local tissue acidification and stimulates increased arylsulfatase B activity to improve locomotor function (Yoo et al., 2013). Therefore, it can be assumed that increased catabolism triggered by lysosomal enzymes may help to maintain the cellular metabolism in low pH-stressed oyster larvae. The increased expression of the 40S and 60S ribosomal subunits revealed in the study also appears to suggest an increase in energy requirements during low pH stress.

In total, 79 DEGs were assigned into four main signaling pathways, including the MAPK signaling pathway (ko04010), the Wnt signaling pathway (ko04310), the vascular endothelial growth factor (VEGF) signaling pathway (ko04370) and the ErbB signaling pathway (ko04012). Among these pathways, the MAPK signaling pathway was significantly activated with 67.6% of total mapped DEGs expressed after low pH treatment. Several MAPK genes are highly expressed in the oyster larvae with a strong coverage in the transcriptome dataset.

The activation of MAPK genes may disrupt many signal transduction pathways in cells including immune response, cell death and other kinase-mediated signaling cascades. Induction of MAPK gene pathways has proven to be a useful biomarker of the responses of marine invertebrates to a variety of environmental stressors (Canesi et al., 2006; Hardege et al., 2011), and MAPK pathways could thus provide information on the responses of early stage oyster larvae to low-pH exposure. However, the Wnt, VEGF and ErbB signaling pathways appeared to be suppressed, with 86.7, 83.3, and 83.3% of the constituent genes down-regulated.

Genes related to calcium transport encoding cAMP-dependent protein kinase catalytic subunit, calcium/calmodulin-dependent protein kinase type II subunit delta and calcineurin subunit B type 1 and Wnt genes were found to be down-regulated during acidification stress in the present study. Wnt genes encode highly conserved signaling molecules that play an important role in embryogenesis, cell differentiation and osteogenesis in both vertebrates and invertebrates (Angers and Moon, 2009; De Robertis, 2010; Holstein, 2012; Karako-Lampert et al., 2014). Decreased calcification and retarded growth might result from the differential expression pattern of Wnt genes under ocean acidification as evidenced from the studies in coral holobiont exposed to low-pH conditions (Karako-Lampert et al., 2014). Recent studies reported hypermethylation of Wnt receptor catabolic process, calcium ion binding, calcium dependent phosphorylation in pediveliger larvae of Crassostrea hongkongensis when exposed to ocean acidification (Lim et al., 2021). Pacific oyster larvae at 14–16 h post fertilization, showed a delay in shell formation with an increase in the expression of tyrosinase gene, carbonic anhydrase, calcium and chitin binding proteins when exposed to pH of 7.5 and 7.6 (De Wit et al., 2018). The expression of genes related to key process in shell formation including calcium transportation, bicarbonate transport and organic matrix where significantly inhibited at the stage of initial shell formation (PDI) of pacific oyster larvae, when exposed to low pH 7.4 (Liu et al., 2020) disrupting the mobilization of calcium ions. Supplementary Table 4 provides complete literature of genes and pathways regulated in response to ocean acidification in marine bivalves across life stages. The expression levels of genes related to chitin synthesis was significantly decreased while the expression of genes related to chitin degradation was found to be increased in D shaped larvae of Crassostrea gigas at pH 7.4, negatively effecting chitin biosynthesis and hence biomineralization (Zhang et al., 2019). In adult pacific oyster, a series of calcium binding proteins and calcium signaling pathways were up-regulated in response to OA indicating regulation of calcium homeostasis under OA (Wang et al., 2020), also in larval sea-urchins (Evans et al., 2013) on the contrary, calcium dependent signaling pathway were suppressed in pearl oyster Pinctada fucata under low pH (Li et al., 2016b) but whereas in case of Crassostrea virginica, no significant expression of these genes was observed upon exposure to ocean acidification (Downey-Wall et al., 2020). The stimulation of calcium related genes as reported in some might be helping the larvae to survive and grow under low pH but in a slower rate (De Wit et al., 2018). However, the initial shell formation requires endogenous energy from metabolic process, the metabolic suppression can lead to alteration in the energy allocation to shell formation delaying or inhibiting the initial shell synthesis (Liu et al., 2020).

Similarly, VEGF and ErbB pathway genes involved in invertebrate cell proliferation and angiogenesis were also down-regulated, suggesting that decreased calcification affects building and repair of skeletal tissues (Green et al., 2013). Most molluscan species have evolved to cope with severe environmental conditions by suppressing their metabolic rate. The degree of metabolic suppression is controlled by the down-regulation of phosphorylation inside the cell (Hochachka et al., 1993). Among the ErbB genes, the negative regulation of the serine/threonine-protein phosphatase 2B catalytic subunit alpha isoform in the low-pH environment may affect the relative activities of phosphatases initiating metabolic suppression (Hand and Hardewig, 1996).

Eight differentially expressed genes obtained by RNA sequencing in response to OA exposure were selected to confirm the accuracy of gene expression by RT-qPCR approach. Further, the fold changes for the validated genes were shown to have a slight variation in their relative gene expression levels when compared between RNA-Seq and RT-qPCR methods (Figure 5).