Ocean Acidification Triggers Cell Signaling, Suppress Immune and Calcification in the Pacific Oyster Larvae

Elevated carbon dioxide levels in ocean waters, an anthropogenic stressor, can alter the chemical equilibrium of seawater through a process called ocean acidification (OA). The resultant reduction of pH can be detrimental during the early developmental stages of the commercially important edible Pacific oyster Crassostrea gigas; the ability of larvae to join a population is likely to be compromised by declining ocean pH. Given this threat, it is important to study the molecular mechanisms that these organisms use to overcome OA stress at the gene expression level. Here, we performed transcriptome profiling in oyster larvae following exposure to ambient (8.1) and reduced (7.4) pH during the pre-settlement growth period (i.e., 18 d post fertilization) using RNA-seq with Illumina sequencing technology. In total, 1,808 differentially expressed genes (DEGs) were identified, 1,410 of which were matched by BLAST against the Swiss-Prot database. Gene ontology classification showed that most of these DEGs were related to ribosomal, calcium ion binding, cell adhesion and apoptotic processes. Pathway enrichment analysis revealed that low pH (7.4) enhanced energy production and organelle biogenesis but prominently suppressed several immune response pathways. Moreover, activation of the MAPK signaling pathway was observed along with inhibition of the Wnt, VEGF, and ErbB pathways, highlighting the fact that the initiation of stress responses is given priority over larval development or shell growth when the larvae cope with low pH. In conclusion, our study demonstrated a unique gene expression profiling approach in studying oyster larval responses to OA, which not only provides comprehensive insights into the mechanisms underlying oyster tolerance to CO2-driven decreases in ocean pH but also supplies a valuable genomic resource for further studies in this species.


INTRODUCTION
Carbon dioxide levels in the atmosphere and oceans have been increasing at unprecedented rates, since the beginning of the industrial revolution (Falkowski et al., 2000). One-third of atmospheric anthropogenic CO 2 emission dissolves in the ocean, raising acidity of the oceans that in turn leads to changes in the carbonate chemistry, the process called as ocean acidification (OA) (Raven et al., 2005). Addition of excess amount of CO 2 in the ocean leads to overall shift in the seawater acid-base chemistry tending toward more acidic condition, reducing pH and saturation state of carbonate minerals (Hurd et al., 2019;Downey-Wall et al., 2020). Subsequent increase in carbonic acid reduces the availability of carbonate ions in oceans, thus affecting marine organisms that depend on ocean carbonate levels to form their calcium carbonate shells Guinotte and Fabry, 2008;Cooley et al., 2009;Thomsen et al., 2018;Hurd et al., 2019). Ocean acidification studies revealed both positive and negative effects on marine organisms (Kroeker et al., 2010(Kroeker et al., , 2011. Though some marine organisms are resilient and can acclimate to changes (Kroeker et al., 2010;Munday et al., 2011;McCulloch et al., 2012), most shell-forming invertebrates, such as sea urchins, corals and mollusks , find it difficult to calcify in the presence of elevated CO 2 levels and reduced pH in seawater (Gazeau et al., 2007(Gazeau et al., , 2010(Gazeau et al., , 2013Kurihara et al., 2007;Miller et al., 2009;Talmage and Gobler, 2009;Todgham and Hofmann, 2009;Sheppard Brennand et al., 2010;Parker et al., 2012;Doney et al., 2020). Ocean acidification also impacts homeostasis and behavioral response in marine fishes (Esbaugh et al., 2012;Jutfelt et al., 2013;Hamilton et al., 2014) and leads to metabolic suppression to conserve energy for survival and immune responses in many marine organisms (Bibby et al., 2007(Bibby et al., , 2008Lannig et al., 2010;Hernroth et al., 2011).
A thorough understanding of the mechanisms that control the physiological processes is needed to uncover the sensitivity of the organism exposed to OA . In this context, studies investigating the molecular mechanisms of responses of marine organisms to ocean acidification assume importance, in evaluating whether they can withstand the deleterious effects associated with it. Although currently available transcriptomic analyses and genomic resources are useful in determining mRNA plasticity in response to OA, the availability of these techniques is still limited to only a few species (Todgham and Hofmann, 2009;Todgham and Stillman, 2013). For example, gene expression studies of early developmental stages in coral identified a number of genes affected by OA (Moya et al., 2012;Rivest et al., 2018). Studies on sea urchins also revealed the overall molecular responses of sea urchin larvae exposed to CO 2 (Todgham and Hofmann, 2009;Devens et al., 2020). OA is known to interfere with several key physiological process in marine organisms, transcriptomics has proven to be effective method to elucidate the molecular level response and when coupled with phenotypic data confers a better understanding of whole organismal response (Davies et al., 2016;Melzner et al., 2020;Strader et al., 2020).
To explore the adaptation and tolerance mechanisms of oyster larvae to low pH, this study performed transcriptome sequencing of early stage (pediveliger) oysters exposed to seawater at a low pH 7.4 (projected to reflect the atmospheric CO 2 emission scenario in 2300, by the Intergovernmental Panel on Climate Change IPCC) (Caldeira and Wickett, 2005), using RNA-seq with Illumina sequencing, to identify underlying adaptation and tolerance mechanisms of oyster larvae coping with environmental pH variations and to provide an integrated view of the cellular pathways affected; to emphasize the physiological responses of ocean acidification.

Larval Culture and CO 2 Perturbation Set Up
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 postfertilization) to the pediveliger stage (i.e., 15-17 day postfertilization) 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 2 levels for the control (pH 8.1) and CO 2 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 2 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 2 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 2 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 TM , 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 1 , K 2 by Mehrbach (Mehrbach et al., 1973) and KSO 4 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 2 levels, low pH 7.4 using next generation high-throughput transcriptome sequencing.

RNA Extraction
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(Dineshram et al., , 2016Griffith 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 Library Construction and RNA-Seq cDNA libraries were constructed following TruSeq TM 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 R 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 Control, Pre-processing, and Genome Mapping
Quality of raw data generated from Illumina sequencing was checked by FastQC. 1 The raw reads were preprocessed using the FASTX 2 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 3 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.

Differential Gene Expression Analysis
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 2 (FC) (≥ 1).

Differentially Expressed Gene Identification, Classification, and Pathway Analysis
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 4 ) with bi-directional best-hit setting. Finally, functional networks of enriched KEGG pathways were visualized using ClueGO. 5

Validation of Transcriptomic Profiles by Real Time-Quantitative Polymerase Chain Reaction
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 Measurements of temperature (Temp), salinity, pH and total alkalinity (TA) were used to derive the other calculated carbonate system components using CO2SYS.xls spreadsheet. Data are presented as the mean ± SD of four measurements. 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).

General Characteristics of Transcriptome Assembly
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.

Differential Expression Genes and Gene Annotation
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)  followed by proteolysis (GO: 0006508), cell adhesion (GO: 0007155) and apoptotic process (GO: 0006915). Almost twothirds of the differentially regulated genes identified in the current study were up-regulated. The proportion of downregulated 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.

Biological Pathways Responsive to Ocean Acidification Acclimation
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 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(Bindea et al., , 2013; Figure 4).

Genes for Energy Production and Cellular Homeostasis
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., 2012Kaniewska et al., , 2015Moya 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 . 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 ( 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 Responsive Genes
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(Goncalves et al., , 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.1fold), 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).

Oxidative Stress and Ribosome Biogenesis Genes
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(Goncalves et al., , 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 FIGURE 4 | ClueGo pathway analysis of enriched differentially expressed genes. Functional group networks of enriched KEGG pathways induced in oyster larvae to low pH are represented using a cerebral layout. GO terms are shown by nodes and are linked based on their kappa score (> 0.3). The node size represents the term enrichment significance, and the labels shown are the most significant leading group terms. Terms with up-regulated and down-regulated genes are shown in red and green, respectively. The color gradient represents the extent of regulation. Terms and pathways with no significance are shown in gray.
Frontiers in Marine Science | www.frontiersin.org 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 . 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 upregulated (≥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.

Signal Transduction and Calcification Regulatory Genes
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 downregulated.
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  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  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 .
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 downregulation 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).

Validation of RNA-Seq by RT-qPCR
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).

CONCLUSION
This study described a unique transcriptome profiling analysis of pediveliger larvae exposed to low pH and also provided comprehensive insights into the mechanism underlying their tolerance to extremely decreased pH. We identified a number of OA-responsive genes that were enriched in various pathways involved in controlling basic physiological processes. In particular, our results revealed that oyster larvae may adopt an energy "trade-off " strategy through the suppression of energydemanding physiological processes, such as immune responses, calcification or shell growth, to overcome the stress induced by OA. In addition, our transcriptome data analyses could also provide a substantial and valuable genomic resource for further studies of this ecologically and economically important aquatic species.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: SRR960964 and SRR967068).

AUTHOR CONTRIBUTIONS
RD and YZ designed the research, analyzed the data, and drafted the final version of the manuscript. RD performed the larval culture experiments with SX, GK, and VT. RD and JL performed transcriptome analysis with the help of YZ and ZY.
RD and KS performed the literature review and discussion of the results. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the National Science Foundation of China (No. 32073002). The larval rearing and ocean acidification facility of this project was funded by a three GRF grants from the HKSAR-RGC (Grant Nos. 17303517, 17304619, and 17120120).

ACKNOWLEDGMENTS
We thank Mr. Tuo Yao for helping in the animal collection from Tsingdao.