Regulation of the Cell Cycle, Apoptosis, and Proline Accumulation Plays an Important Role in the Stress Response of the Eastern Oyster Crassostrea Virginica

Background Understanding how organisms respond and adapt to environmental changes is central to evolutionary biology. As a sessile organism that has adapted to life in estuaries and intertidal zones, the eastern oyster Crassostrea virginica can tolerate wide fluctuations in temperature and salinity and survive for weeks out of water. To understand the molecular mechanisms underlying the remarkable stress tolerance of the eastern oyster, we studied the transcriptomic changes induced by exposure to air and cold stress. Eastern oysters were maintained for 7 days under four conditions, namely, in seawater (normal) at 22°C, in air at 22°C, in seawater at 5°C and in air at 5°C, and then sampled for RNA sequencing. Results Transcriptomic analysis revealed that many genes involved in cell cycle progression and DNA replication were downregulated in oysters exposed to air and cold, which indicates that stress inhibits cell division. Exposure to air at 22°C induced a concerted inhibition of apoptosis through the upregulation of expanded inhibitors of apoptosis and the downregulation of caspases. Interactions between TNF and NF-κB signalling implied a reduction in the inflammatory response and immune functions. Key genes for proline production, fatty acid synthesis and chromosomal proteins were upregulated during exposure to low temperatures, which suggested that proline accumulation, energy conservation, and epigenetic modification of chromosomes are important for coping with cold stress. The upregulation of melatonin, FMRFamide, and neural acetylcholine receptors indicate the significance of the neurohormonal regulation of homeostasis. Conclusion These results show that air exposure and cold stress alter the expression of key genes for cell division, apoptosis, proline accumulation, fatty acid metabolism, neurohormonal signalling, and epigenetic modifications, suggesting regulation of these processes plays an important role in the stress response of the eastern oyster and possibly other marine molluscs. This study provides new insights into molecular mechanisms of stress response that are essential for understanding the adaptive potential of marine organisms under climate change.


INTRODUCTION
The mechanism through which organisms respond and adapt to environmental shifts is a fundamental question of evolutionary biology. As climate change causes significant shifts in the oceanic environment, understanding its effects on marine organisms and ecosystems is critical. The eastern oyster (Crassostrea virginica) is a marine bivalve widely distributed along the Atlantic coast of America, ranging from the Gulf of St. Lawrence in Canada to the Gulf of Mexico and extending to West Indies, and an important fishery and aquaculture species in the US and Canada. Due to its adaptation to estuarine and intertidal environments across a wide latitudinal range, the eastern oyster is remarkably resilient and can tolerate wide fluctuations in salinity and temperature as well as long periods of exposure to air. In the intertidal zone, the eastern oyster may face icy conditions in the winter and 45°C in the summer (Galtsoff, 1964;Shumway, 1996). This species can survive well for at least a month during exposure to air at 4°C (Kawabe et al., 2010). This remarkable resilience makes the eastern oyster an interesting model for studying stress adaptation.
Among all environmental challenges, prolonged exposure to air may be very stressful for oysters and can cause heavy mortalities during summer months (Malek, 2010;Clements et al., 2018). As observed in field experiments, increases in intertidal air exposure increase disease prevalence and mortality and reduce growth (Malek, 2010). Genomic analyses have provided some insights into the molecular mechanisms of the stress response and adaptation of oysters. Transcriptomic analyses of the eastern oyster have identified genes and pathways related to the immune responses of the species and its responses to salinity and crude oil (Eierman and Hare, 2014;McDowell et al., 2014;Zhang et al., 2014a;Loṕez-Landavery et al., 2019;Modak and Gomez-Chiarri, 2020). Transcriptomic studies of Pacific oysters Crassostrea gigas subjected to thermal, salinity, air exposure and heavy metal stresses have shown that air exposure induces transcriptional changes in the largest number of genes (Zhang et al., 2012). These changes include strong upregulation of the expression of heat shock proteins (HSPs) and inhibitor of apoptosis (IAPs), which suggests that these proteins play key roles in stress response and adaptation (Zhang et al., 2012;Guo et al., 2015). A proteomic study of Pacific oysters under desiccation stress revealed increases in ubiquinone synthesis, reductions in reactive oxygen species (ROS), and downregulation of the expression of proteins involved in metabolism, ion transport, DNA replication and protein synthesis (Zhang et al., 2014b). In the eastern oyster, air exposure induces the upregulation and alternative splicing of alternative oxidases that participate in ROS removal (Liu and Guo, 2017).
Few studies have examined the relationship between air exposure and the innate immune system. Allen and Burnett (2008) found that hypoxia and high temperatures reduce the ability of haemocytes to kill bacteria in C. gigas under air exposure. Air exposure may result in reduced oxygen levels in oysters or hypoxia. Alvarez and Friedl (1992) reported that hypoxia and anoxia have no effects on phagocytosis and particle clearance in the eastern oyster, whereas another study (Boyd and Burnett, 1999) in the same species showed that the respiratory burst function (ROS production) of haemocytes is reduced under hypoxic conditions. Environmental stressors exert immunosuppressive effects in Pinctada imbricata (Kuchel et al., 2010). The eastern oyster populations are seriously impacted by two lethal diseases, MSX caused by Haplosporidium nelsoni and Dermo caused by Perkinsus marinus (Guo and Ford, 2016). The effects of air exposure on the immune function of the eastern oyster are unclear.
Compared with air exposure, cold exposure is a relatively mild stressor for oysters. The LT50 (median lethal time to reach 50% mortality) values for oysters exposed to air at temperatures of 4°C, 15°C and 20°C are 47.8, 15.9 and 12.2 days, respectively (Kawabe et al., 2010). The regulation of energy storage and consumption is crucial for the survival of animals during long periods of cold conditions (Storey and Storey, 1990). During cold periods, bivalves can enter a state of metabolic rate depression (MRD) (Saucedo et al., 2004;Xiao et al., 2014). An analysis of the physiological states of oysters at low temperatures has suggested that the metabolic pathways in oysters change after one week of exposure (Kawabe et al., 2010), although most of the regulatory networks at the protein and transcriptomic levels are largely unknown.
The availability of the C. virginica genome and highthroughput next-generation sequencing technology has allowed accurate measurement of transcriptional levels of all genes, which provided a valuable tool for studying molecular mechanisms of stress response in oysters (Zhang et al., 2012;Goḿez-Chiarria et al., 2015;Guo et al., 2015;Guo and Ford, 2016). In this study, we used the reference genome (GCA_002022765.4 C_virginica-3.0) and transcriptomic sequencing to characterize transcriptional changes in eastern oysters subjected to air exposure and cold stress and thus elucidate the molecular mechanisms or regulatory pathways underlying stress response in this species.

Air Exposure and Cold Stress
Hatchery-produced eastern oysters of the same age (18 months) and similar shell lengths (5-8 cm) were used in this study. The oysters were obtained from the Haskin Shellfish Research Laboratory, Rutgers University, Port Norris, New Jersey, USA. The oysters were maintained in ambient seawater at 9°C~12°C and then acclimated to 22°C for 3 days for this experiment in April. To elicit a strong response to chronic air exposure and cold stress, oysters were exposed to air and cold (5°C) for 7 days. The low temperature of 5°C was selected because it is the lowest possible temperature above freezing that can be stably maintained in the lab. The duration of 7 days was selected based on a previous study (Zhang et al., 2012), where oysters exposed to 5°C for 7 days and to air for 7-9 days exhibited the strongest responses. Although 7-day exposure is rare in the field (only experienced by oysters at extremely high tidal positions), this extreme length helps to identify key genes responding to air exposure stress. Also, eastern oysters are often stored in cold rooms without water for 1-2 weeks before being consumed live. Taking all these factors into account, we subjected oysters to 7 days of air exposure. In winter, the coldest month in New Jersey is January when the average temperature ranged from -6.7°C to 3.3°C. The oysters may be covered by ice for several days at a time in the field. The low temperate tested in this study is quite common for eastern oysters in NJ and the northeastern US. The average ambient temperature in April was 10°C. The control oysters were maintained in seawater at 22°C with heaters, an optimal temperature that is routinely used for oyster conditioning in our lab. Accordingly, the oysters were divided into 4 groups (7 per group): one group was maintained in seawater at 22°C with aeration and served as the control group (hereafter the 22SW group), one group was exposed to air in an incubator with temperature set at 22°C (the 22AE group), one group was maintained in seawater at 5°C (the 5SW group), and one group was exposed to air at 5°C (the 5AE group). The cold treatment was conducted in a cold storage room with both water and air temperature set at 5°C. The temperature decreased from 22°C to 5°C over approximately one day. The oysters exposed to air were covered with a wet towel to keep them moist, and the towel was rewetted daily. The oysters kept in seawater were not fed to ensure consistency with those exposed to air. The salinity of the seawater was 32 parts per thousand.

Total RNA Extraction
After 7 days of treatment, 6 oysters were sampled from each group. Gill tissue from each oyster was flash-frozen in liquid nitrogen and used for transcriptome sequencing. Gills were used for this study because gills are the organ needed for respiration and are important for maintaining homeostasis. Total RNA was extracted from the gill tissue using TRIzol reagent (Invitrogen, USA) following the manufacturer's protocol. The quantity and quality of the extracted total RNA were assessed with a NanoDrop 2000 (Thermo Scientific, USA) and agarose gel electrophoresis. For each group, RNA was extracted from six oysters, and equal amounts of RNA from each oyster were combined into one sample. A total of 3 mg of RNA per group was used for library preparation.

Library Preparation and Bulk RNA Sequencing
Sequencing libraries were generated using an NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) following the manufacturer's recommendations, and unique index codes were utilised to distinguish the samples. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was performed using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H -). Second-strand cDNA synthesis was performed using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends by exonuclease/ polymerase. After adenylation of the 3' ends, NEBNext adaptors with hairpin loop structures were ligated to the fragments to prepare them for hybridization. To preferentially select cDNA fragments with a length of 150~200 bp, the library fragments were purified with an AMPure XP System (Beckman Coulter, USA). The size-selected, adaptor-ligated cDNA was then incubated with 3 ml of USER Enzyme at 37°C for 15 min and then at 95°C for 5 min before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and the Index (X) Primer. The PCR products were purified (AMPure XP System), and the library quality was assessed with an Agilent Bioanalyzer 2100 system. Clustering of the indexcoded samples was performed using a cBot Cluster Generation System with a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to manufacturer's instructions. After cluster generation, the libraries were sequenced on an Illumina HiSeq 2000 platform to generate 100-bp paired-end reads.

Mapping of Reads Against the C. Virginica Genome
Raw sequence data in fastq format were first processed through inhouse Perl scripts. Clean reads were obtained by removing reads containing adapters or poly-N sequences (for which the proportion of N bases was > 10%) and reads containing more than 50% of low quality (Q-value ≤ 20) bases. The C. virginica reference genome and gene model annotation files were downloaded from NCBI (https:// www.ncbi.nlm.nih.gov/genome/398). An index of the reference genome was built, and clean paired-end reads were mapped to the reference genome using HISAT2.2.4 (Kim et al., 2015) with the "-rna-strandness RF" setting and default settings for other parameters.

Gene Expression Analysis
The mapped reads of each sample were assembled with StringTie v1.3.1 (Pertea et al., 2015;Pertea et al., 2016) using a reference-based approach. For each transcription region, fragments per kilobase of transcript per million mapped reads (FPKM) values were calculated to quantify transcript abundance of each gene using StringTie software. Differential expression between two samples was analysed with edgeR (Robinson et al., 2010;McCarthy et al., 2012). P-values were obtained using the exactTest function with normalized counts and assuming a negative binomial distribution. Correction for multiple testing was performed by applying the Benjamini-Hochberg method to the p-values to control the false discovery rate (FDR) (Benjamini and Hochberg, 1995). To reduce false positives due to the lack of replication, a fold-change threshold was also applied. Genes with a false discovery rate (FDRs) < 0.001 and absolute fold change ≥ 2 were considered DEGs. A hierarchical cluster analysis of the top 1000 DEGs was performed to assess transcriptional variations among oysters exposed to different stresses using in-house R scripts.

Screening of Stress Response Genes
Gene Ontology (GO) enrichment analysis of the DEGs was performed with the GOseq R package (Young et al., 2010). We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology-Based Annotation System (KOBAS) (Mao et al., 2005) to test the significance of the enrichment of DEGs in KEGG (Kanehisa et al., 2008) pathways. KOBAS uses the whole genome as the default background distribution. For each pathway that occurs in a geneset, it counts the number of genes in the geneset that belong to the pathway and compares it against the whole-genome background. The calculating formula of P-value is: where N is the number of all genes that with KEGG annotation, n is the number of DEGs in N, M is the number of all genes annotated to specific pathways, and m is number of DEGs in M. The calculated p-values were corrected for FDR (Benjamini and Hochberg, 1995) with FDR ≤ 0.05 as the threshold. Pathways meeting this condition were defined as significantly enriched pathways in DEGs.

qPCR Validation of the RNA-seq Results
Before combining RNA from 6 oysters for library construction, 2 µg of total RNA from each oyster was removed and separately reverse-transcribed into cDNA for qPCR validation. Complementary DNA (cDNA) was prepared using a FastQuant RT kit (with gDNase) (Tiangen, Beijing, China).
We selected eight genes for qPCR validation rather arbitrarily to cover genes with different expression levels. Most of the eight genes were significantly expressed in at least one treatment group, but not all of them were DEGs, such as EC-SOD (Extracellular superoxide dismutaseare). They were chosen to show the correlation between the read counts obtained by RNAseq and the relative expression levels obtained by qPCR. Some genes were abundantly expressed (KLF5) and some were not (EEF2). At the same time, they may or may not be key genes or in important pathways analyzed. Primers for the eight genes were designed with Primer-BLAST online (http://www.ncbi.nlm.nih. gov/tools/primer-blast/) (Supplementary Table 2). The expression of the selected genes and the housekeeping gene bactin (Ivanina et al., 2010) was quantified by qPCR using a 7500 Fast Real-time PCR System (Applied Biosystems, USA). The cycle threshold (CT) values of each gene were normalized to those of b-actin, and the relative amount of mRNA specific to each of the selected genes was calculated using the 2− DDCt method (Livak and Schmittgen, 2001). Validation analyses were performed by comparing the read count of each gene obtained from the RNA-seq data (adjusted by edgeR software using one scaling normalized factor) with the mean 2− DDCt value from the qPCR results.

Transcriptome Mapping and General Patterns of Gene Expression
A total of 247,921,390 clean paired-end reads were generated from the four groups of oysters, and these numbers ranged from 55.2 to 67.9 million per group (Supplementary Table 1 Table 2), a highly significant correlation (p < 0.0001, r = 0.8553) was found between the read counts obtained by RNA-seq and the relative expression levels obtained by qPCR ( Supplementary Figure 1), which indicated that the RNA-seq data are reliable due to the pooling of six individuals and a high read coverage. Further, we applied more stringent selection criteria (FDR < 0.001 plus foldchange ≥ 2) to DEG identification to compensate for the lack of sequencing replication, and identified affected pathways based on concerted regulation of multiple protagonist and antagonist genes (see below) so that the conclusions are not affected by random false positives. Nevertheless, it should be noted that the inclusion of independent samples or sample pools would increase the statistical power and detect more DEGs with small fold changes that could not be identified in this study. Additionally, findings of this study should be viewed with caution due to the lack of treatment replicates. Finally, the response to stress is a temporal process, and our observations are limited to the response at day 7, which may differ from the responses at other time points.
Differential expression analysis revealed large numbers of DEGs among the four groups. Compared with the control condition (22SW), exposure to air at 22°C resulted in the largest number of DEGs (2488), which suggested that air at 22°C is the harshest of the tested conditions for the eastern oyster ( Figure 1A). Exposure to air at 5°C produced 1969 DEGs, and exposure to seawater at 5°C produced 1802 DEGs, which indicated that these conditions are relatively less stressful. The two conditions at 5°C shared 987 DEGs ( Figure 1B). Under all three stress conditions, a higher number of genes were downregulated than upregulated. These results indicate that all three stress conditions depressed the transcriptional activities of the eastern oyster and that air exposure at 22°C was the most stressful. This finding is consistent with the results obtained for the Pacific oyster, in which air exposure also induces the highest number (4420) of DEGs (Zhang et al., 2012).
The three stress conditions differed among the types of DEGs induced. Although the three conditions shared 636 DEGs, they also had large numbers of unique DEGs: 1296 in 22AE, 713 in 5AE and 528 in 5SW ( Figure 1B). Among the four groups, the expression profiles obtained for 5AE and 5SW were most similar but differed significantly from that obtained for 22AE, as indicated by the clustering heatmap ( Figure 1C). As low temperature tends to reduce metabolic activity, oysters might be less susceptible to air exposure at low temperature than at normal or high temperatures. All DEGs revealed by pairwise comparisons are shown in Supplementary Table 3.
Enrichment analysis revealed that the following KEGG pathways were overrepresented in the DEGs between the 22AE and 22SW control groups: DNA replication, cell cycle control, retinol metabolism, apoptosis, nuclear factor kappa B (NF-kB) signalling, lipid metabolism, peroxisome and tumour necrosis factor (TNF) signalling (P-value < 0.05, Table 1). The main processes affected in the 5SW and 5AE groups were lipid metabolism and DNA replication. The only GO term significantly enriched in DEGs between 22AE and 5AE was DNA metabolism (0006259), and these DEGs included genes involved in DNA replication and cell cycle control. Some of the significantly enriched pathways and processes are highlighted below. For brevity, the gene IDs from the reference genome were shortened by omitting the string "OC1111", which is shared by all the genes (i.e., L00001 for LOC111100001).

Cell Cycle Control Under Stress
One of the most significantly affected processes in the eastern oyster under stress is cell division, as evidenced by the identification of large number of DEGs related to DNA replication and cell cycle control under all three stressors at day 7, particularly in the 22AE group (Table 1). Cell cycle progression, a process that is conserved in metazoans, is regulated by cyclin-dependent kinases (CDKs) and their specific cyclin partners, positive and negative regulators, and downstream targets (Morgan, 1997). In this study, cyclin A (L22860, L22939), cyclin D (L03611), cyclin E (L32839, L99092), and CDK2 (L21435), which are expressed at the G1 phase and determine whether cells proceed through the G1 phase into the S phase, were all downregulated under air normalized (read count+1) values from high (red, the maximum value is 8.14) to low (green, the minimum value is -6.82). Table 4 and  Table 4). All six genes under cold stress showed reduced expression, although the reduction in the expression of some of the genes was not significant. Once replication complexes are built and activated, a series of replicating forks form along the chromosome, generating bubbles that eventually fuse together to complete DNA replication (Alberts et al., 2010). Genes associated with the replication protein A complex, such as replication protein A (L20888), DNA polymerase (L12959, L30346, L37485, L30057, L26184, L28010), DNA primase (L26192, L21347) and DNA clamp proteins (RFC and PCNA, L08767, L31355, L32381, L34900, L22897, L35799, L20673), were downregulated (Supplementary Table 4). Genes encoding cyclins and CDKs that promote the progression of cells from the G2 phase to the M phase, such as G2/mitotic-specific cyclin-A (L22860, L22939) and M-phase inducer phosphatase (CDC25A, L32091, L32548, L35070), which function as dosage-dependent inducers of mitotic progression (Nilsson and Hoffmann, 2000), were also downregulated ( Figure 2).

exposure and cold stress (Supplementary
Other genes, such as CDK7 (L25401), involved in transcription initiation and DNA repair were also downregulated. These results strongly suggest that exposure of the eastern oyster to air and cold stress suppresses the transcription of key genes needed for cell cycle progression or cell division, which is crucial for growth and immune functions. The discovery of the downregulation of key genes needed for cell division provides molecular support for the well-known negative effects of stress on growth and the immune response (Malek, 2010;Guo et al., 2015;Guo and Ford, 2016;Clements et al., 2018;La Peyre et al., 2018). Cell cycle arrest may be induced by DNA damage, and cells attempt to repair this damage during cell cycle arrest (Pucci et al., 2000). In this study, genes that are induced by DNA damage, such as the growth arrest and DNA damage-inducible protein GADD45 (L18431, L18444, L18682, L21109), which is induced by p53 and arrests the cell cycle at the G2/M checkpoint (Papa et al., 2004), were upregulated. High-mobility group proteins  (L09677, L09703, L12806, L17359, L25762, L26226), which respond to DNA damage and oxidative stress in yeast (Reeves and Adair, 2005;Lange et al., 2008), were also upregulated. In addition, sestrin1 (L03315), which is induced by the p53 tumour suppressor protein and has the ability to arrest cell growth and decrease ROS production, was also upregulated (Sun et al., 2014). It is widely accepted that p53 mediates cell cycle arrest, apoptotic cell death, and/or cellular senescence . Although the transcription level of p53 did not change significantly under stress in this study, it is possible that p53 activity, a key event in the cell cycle arrest process, is controlled by a diverse array of posttranslational modifications (Kruse and Gu, 2009). Additionally, because this study only provides a snapshot in time, p53 and its regulators may be differentially expressed at other time points. It should be noted that in addition to DNA damage repair, stress and nutritional limitations may also lead to cell cycle arrest and growth arrest (Yanagida et al., 2011;Rossi and Antonangeli, 2015). Although the underlying mechanisms are not completely defined, cell cycle arrest is an active response to stress that enables organisms to survive under challenging environmental conditions. The results of this study demonstrate that the exposure of eastern oyster to air and cold stress suppresses the transcription of genes promoting cell cycle progression, which may lead to cell cycle arrest and in turn affect growth, immune responses and other biological processes. Growth retardation inferred from cell cycle arrest is consistent with the results from field experiments conducted by Malek (2010), who found that the growth of eastern oysters decreased with increases in the duration of intertidal air exposure, and La Peyre et al. (2018), who found that a lower growth rate among oysters that experienced daily air exposure. Air exposure in the field may also reduce the time dedicated to feeding and increase the use of glycogen storage, which contribute to reduced growth.

Evasion of Cell Death during Exposure to Air at 22°C
Exposure to air at 22°C triggers the inhibition of apoptosis and related regulatory pathways. Inhibitors of apoptosis (IAPs) play important roles in promoting cell survival. IAP family members inhibit programmed cell death by binding to and/or directly inhibiting caspases (LaCasse et al., 1998). The C. gigas genome exhibits an expanded IAP family (48 genes), members of which are highly expressed after air exposure and pathogen infection, which indicates that IAPs play crucial roles in anti-apoptotic processes in oysters (Zhang et al., 2012;Green et al., 2015). We identified 58 IAPs in the eastern oyster genome; 8 and 4 of these IAPs were highly and significantly upregulated and slightly downregulated, respectively, in the 22AE group, resulting in an overall upregulation of IAPs (Supplementary Table 4). Furthermore, genes encoding caspases, CASP7 (L3317, L2489, L24891) and CASP10 (L04823, L33185, L33189), were slightly downregulated in the 22AE group. These results suggest the existence of a concerted effort to inhibit apoptosis in oysters under stress (Figure 3 and Supplementary Table 4).

Xenobiotic and Lipid Metabolism During Exposure to Air at 22°C
The key genes and pathways involved in xenobiotic biotransformation and lipid metabolism were affected during exposure to air at 22°C ( Table 1). The xenobiotic biotransformation system helps oysters cope with environmental contaminants (Zhang et al., 2016). Genes that are known to be important for xenobiotic biotransformation, such as cytochrome P450s (CYP450s, 11 genes), myeloperoxidase (MPO, L099927, L01491, L01639, L01744), GSTs (L02975, L04070), and ATPbinding cassette subfamily D member (L30560) (Supplementary Table 4), were signifi cantly downregulated under 22AE conditions. Genes related to fatty acid synthesis, particularly the key gene fatty acid synthase (FAS, L18924, L21707), were highly upregulated under 22AE conditions. Moreover, genes playing key roles in fatty acid oxidation in peroxisomes, such as peroxisomal acyl-coenzyme A oxidase 1 (ACOX1, L13840, L15085, L32545, MSTRG.2967), which is the first enzyme in the fatty acid beta-oxidation pathway that produces hydrogen peroxide (Inestrosa et al., 1979), were downregulated. The 3ketoacyl-CoA thiolase (ACAA1, L15744, L15745), which catalyses the final step in the b-oxidation cycle, was also downregulated. The upregulation and downregulation of these key genes were also observed in the groups exposed to cold stress, albeit at a small magnitude, which suggested similar regulation of fatty acid metabolism under cold stress (Supplementary Table 4). Maintaining the energy balance is critical for stress tolerance (Sokolova et al., 2012), and our results indicate the involvement of lipid metabolism in coping with air exposure and cold stress. Additionally, organisms in cold environments often have more unsaturated fatty acids for membrane fluidity (Dey et al., 1993).
The downregulation of xenobiotic biotransformation and peroxisome-related genes may decrease fatty acid oxidation (Mayer et al., 1995), whereas the upregulation of fatty acid synthase increases the de novo synthesis of fatty acids; altogether, these effects alter the fatty acid composition. These results all indicate that fatty acid metabolism plays a critical role in mediating the responses of oysters to air exposure.

Proline Accumulation and Histone Upregulation Under Cold Stress
The proline accumulation pathway was activated in oysters under cold stress (Figure 4). The synthesis of proline is mediated by two pathways: the L-glutamic acid pathway and the ornithine pathway (Verbruggen and Hermans, 2008). Delta-1-pyrroline-5-carboxylate synthase (P5CS, L36221), a key enzyme responsible for the synthesis of L-glutamate 5semialdehyde from L-glutamate in steps 1 and 2 of the subpathway, was significantly upregulated (Figure 4 and Supplementary Table 4). Arginase (ARG, L30979), which transforms arginine to ornithine, and ornithine aminotransferase (OAT, L36290), which catalyses the transformation of ornithine into L-glutamate 5-semialdehyde, were both upregulated by cold stress. The enzyme pyrroline-5carboxylate reductase (PYCR2, L28747, PYCR3, L28753), which catalyses the last step in proline biosynthesis, was also upregulated. Together, these results suggest that the proline accumulation pathway is activated or upregulated under cold stress. The accumulation of proline has been observed in oysters under hyperosmotic stress (Meng et al., 2013) and is a key mechanism for cold tolerance in plants (Szabados and Savoure, 2010). The activation of the proline accumulation pathway observed in this study implies that a large proline pool is important for the survival of oysters under cold stress conditions. Nine genes for histone H4 and four genes for histone H2B were highly upregulated under cold stress, particularly under 5SW conditions. The upregulation of some but not all histone subunits was also observed in the Pacific oyster under stress (Zhang et al., 2012). Canonical histones are synthesized during the S phase for the packaging of newly synthesized DNA. Because genes that promote cell division were downregulated under stress, the upregulation of some histone genes is intriguing and begs other explanations. Moving, ejecting or restructuring nucleosomes is one of the strategies for maintaining chromatin stability (Cedar and Bergman, 2009;Clapier and Cairns, 2009). The high-mobility groups (HMGs) of proteins (L09677, L09703, L12806, L17359, L25762, L26226), another group of major chromatin architectural components of eukaryotes, were also highly upregulated under cold stress (Supplementary Table 4).
Although preliminary and to be confirmed through further studies, these results indicate that genes encoding chromosomal proteins are upregulated under cold stress, probably for repairing or stabilizing DNA under stress or as epigenetic modifications that regulate transcription. In plants, stress induces dramatic epigenetic changes, such as DNA methylation, histone modifications and changes in the composition of different histone variants, which regulate gene expression and enable stress responses and memory (Chinnusamy and Zhu, 2009). We hypothesize that the upregulation of the expression of certain chromosomal structural protein genes observed in this study is involved in the epigenetic modifications of chromatin that are needed for regulating the transcriptional response to stress and facilitating epigenetic memory in oysters.

Upregulation of Hormone and Neurotransmitter Receptors
Several genes encoding hormone and neurotransmitter receptors were upregulated in oysters under air exposure and cold stress conditions. The melatonin receptor (MTNR, L03709, L04987, L03710) and the FMRFamide receptor (FMRFAR, L22789, L22537, L21307) were upregulated under all three conditions, particularly under 22AE conditions (Supplementary Table 4). FMRFamide and FMRFamide-related neuropeptides are abundant in invertebrates and have several important functions. FMRFamide-related peptides (FaRPs) in gastropods are involved in movement, feeding, defecation, and reproduction (Loṕez-Vera et al., 2008). Melatonin is involved in circadian and seasonal timing (Reppert et al., 1994;Dubocovich et al., 2003) and was highly upregulated under all stresses in this study. In addition, genes encoding atrial natriuretic peptide receptor A and neuronal acetylcholine receptor were also upregulated under 5AE conditions. The up regulation of these hormone and neurotransmitter receptors may be critical for activating physiological responses or inducing physiological state transitions as in a circadian system shift (Ottaviani and Franceschi, 1996;Fuentes-Pardo et al., 2013;Lenz et al., 2015). In C. gigas, many GPCR genes, including those encoding prostaglandin E2 receptors, FMRFamide receptors, cholecystokinin receptors, melatonin receptors, prolactin-releasing peptide receptors, adenosine receptors and dopamine D2-like receptors, are upregulated by OsHV-1 infection . It has been shown that neuronal acetylcholine receptors are greatly expanded in bivalve FIGURE 4 | Effects of cold stress on the proline accumulation pathway in eastern oysters. The genes in red exhibit upregulated expression. The gene names are the following: P5CS, delta-1-pyrroline-5-carboxylate synthase; PYCR, pyrroline-5-carboxylate reductase; ARG, arginase; OAT, ornithine aminotransferase. molluscs and are important for adaptation to stationary life under variable environments (Jiao et al., 2019).

CONCLUSIONS
The remarkable resilience of the eastern oyster is not well understood at the molecular level. The results obtained in this study clearly demonstrate that exposure to air and cold stress induce significant changes in the transcription of large numbers of genes that are likely crucial for the response and tolerance of the eastern oyster to stress. The observed suppression of genes promoting cell cycle progression suggests that cell division may be arrested in oysters under stress. The upregulation of antiapoptosis genes and the downregulation of pro-apoptosis genes and regulatory pathways indicate a concerted effort to evade cell death. The inhibition of key signalling pathways that activate the immune response suggests that immune defence is impaired by stress. Receptors for hormones and neurotransmitters were found to be upregulated in response to stress. The upregulation of genes for proline synthesis and chromosome structural proteins indicates that proline accumulation and epigenetic modifications of chromosomes are important for coping with cold stress. Although it is well known that stress depresses growth and immune responses, which often lead to disease outbreaks and mass summer mortalities in molluscs (Guo and Ford, 2016), the underlying molecular mechanisms are not well understood. By revealing the downregulation of genes involved in cell cycle progression and the inflammatory response under stress, this study provides transcriptomic evidence for the negative effects of stress on growth and immune functions. This study reveals several previously undescribed responses to stress in oysters and other molluscs. The genes and pathways implicated in this study provide useful resources for further studies on the stress response and adaptation in molluscs and other marine invertebrates. Knowing how stress affects growth and immune response at molecular levels is essential for understanding stress impact and assessing the adaptive potential of marine organisms under climate change. Oysters are economically important and support major aquaculture industries worldwide. The sustainable development of the oyster aquaculture industry depends on genetic improvement of cultured stocks (Guo, 2021). The candidate genes and their polymorphisms can be used for marker-assisted selection for the breeding of more resilient oysters.

AUTHORS CONTRIBUTIONS
XG and CL conceived and designed the study. CL performed the experiments and analysed the data. CL, XG and HW interpreted the data and wrote the manuscript. All authors have read and approved the manuscript. All authors contributed to the article and approved the submitted version.