The experimental design of spawning referred to our previous study (Li A. et al., 2017). In short, we collected two congeneric oyster species, C. gigas and C. angulata, from their natural habitat in Qingdao (35°44′ N) and Xiamen (24°33′ N), respectively. According to our measurements, the mean air and sea surface temperatures at the southern site were 8.7°C and 6.6°C higher than those at the northern site over the past 2 years, respectively (Li A. et al., 2017). To alleviate environmental affect, we conducted a one-generation common garden experiment (Sanford and Kelly, 2011). For each species, we mixed eggs from 30 mature females and divided them into 30 beakers. Sperm from each of the 30 mature males was crossed separately with the eggs in each beaker. Zygotes were fertilized and larvae were cultured in indoor hatchery of Laizhou city (37°18′N, Shandong province, China). The fertilization was conducted in a 100 L bucket at 20 ∼ 22°C and salinity for 28 ∼ 30‰. The larvae were cultured in 20 m³ indoor pond at 22 ∼ 24°C and 28 ∼ 30‰ for salinity, feeding with 1.5 L 30,000 ∼ 40,000 cells/ml mixed algae [golden algae (Chrysophyta) and Platymonas subcordiformis] five times a day. One month later, 3-month intermediate rearing at the sea of Laizhou city were conducted, during which the cultured density was decreased gradually. Then, 4-month-old F₁ progenies were transported to the sea of Muping city (37°39′N, Shandong province, China) for grow-out culture. After 2 months of acclimatization, we sampled 15 adult progeny oysters, gill tissues, and other tissues of each species in situ for subsequent experiments. In addition, we collected C. gigas and C. angulata from Muping and brought them back to Qingdao to measure their oxygen consumption rates.

The gill tissues were used for RNA extraction. To increase the accuracy of the sequencing, we increased the sample number to eight individuals of each species. RNA from the 16 oyster tissues were extracted using TRIzol (Invitrogen, United Kingdom) for RNA sequencing by BMK (Biomarker Technologies, China). Each sample was used to generate an independent library with NEBNext Ultra^TM RNA Library Prep Kit for Illumina (NEB, United States). RNA was enriched using Oligo(dT) beads. Then, the enriched mRNA was divided into short fragments using a fragmentation buffer, and it was reverse transcribed into cDNA using random primers. Second-strand cDNA was synthesized using DNA polymerase I, RNase H, dNTP, and buffer solution. Afterward, the cDNA fragments were purified using AMPure XP beads (Beckman Coulter, Beverly, MA, United States), the ends were repaired, poly(A) tails were added, and the fragments were ligated with Illumina adapter sequences. The ligation products were size-selected using AMPure XP beads, amplified by PCR, and sequenced using an Illumina NovaSeq 6000. To obtain high-quality clean reads, we filtered the reads by removing the adapter sequences, reads with >10% unknown nucleotides (N), and those with >50% low-quality bases (Q-value ≤ 20).

The HISAT2 v2.0.4 package (Kim et al., 2015) was used to map reads to the oyster genome (GenBank accession No., GCA_011032805.1) with a parameter setting of “–dta -p 6 –max-intronlen 5000000.” Then, the SAM files were converted to BAM files using SAMtools v0.1.18 (Li et al., 2009). The StringTie v1.3.4d package (Pertea et al., 2015) was used to process the read alignments and annotations with a parameter setting of “–merge -F 0.1 -T 0.1.” StringTie was then used to estimate abundance and create a new transcript table for input to the R package “ballgown” (Pertea et al., 2016). Gene expression levels were quantified using “ballgown” based on the expected number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) method (Florea et al., 2013).

Differential expression analysis of C. gigas and C. angulata was performed using the “DESeq2” R package (version 1.6.3) (Love et al., 2014). The Benjamini-Hochberg method was used to control the false discovery rate (FDR) by adjusting the P-values. Genes with FDR ≤ 0.01 and a log₂ fold change (FC) of |log₂FC| > 1 were identified as differentially expressed genes (DEGs). Principal component analysis (PCA) and hierarchical cluster analysis (HCA) were performed using the “ropls” R package (version 3.3.2). The DEGs were enriched by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the clusterProfiler R package (version 3.10.1) (Yu et al., 2012). GO is an international standardized gene function classification system consisting of terms that use controlled vocabulary to provide a comprehensive representation of gene function. KEGG is a major public pathway-related database (Kanehisa et al., 2008) that can identify significantly enriched metabolic pathways in DEGs from the background of the whole genome¹. DEGs with a Q-value ≤ 0.05 were considered to be significantly enriched pathways.

Twenty gill tissues (i.e., ten from C. gigas and ten from C. angulata) were used for metabolite extraction and detection. Metabolites were extracted according to the protocol described by BMK (Biomarker Technologies, China). Approximately (50 ± 1 mg) of each oyster tissue was extracted using 1 mL extraction liquid (V_Methanol:V_Acetonitrile:V_Water = 2:2:1) and 2 μL 2-Chloro-L-phenylalanine, which was used as an internal standard. The extract was vortexed for 30 s. Subsequently, the extract was homogenized using a ball mill for 10 min at 45 Hz and then ultrasonicated for 5 min at 4°C. The solution was incubated for 1 h at −20°C. Then, the mixture was centrifuged for 15 min at 12,000 rpm and 4°C. Approximately 500 μL of the supernatant was transferred into 1.5 mL microcentrifuge tubes. The extracts were dried completely in a vacuum concentrator without heating. Then, 160 μL of extraction liquid (V_Acetonitrile:V_Water = 1:1) was added to the dried metabolites, vortexed for 30 s, and ultrasonicated for 15 min at 4°C. The solution was centrifuged for 15 min at 12,000 rpm and 4°C. Approximately 120 μL of the supernatant was transferred into 2 mL sample bottles, and then 10 μL of the extract from each sample was mixed with a QC sample for detection. All samples were analyzed using a liquid chromatograph coupled with a quadrupole time-of-flight mass spectrometer (LC-QTOF-MS).

LC-QTOF-MS analysis was performed using an Acquity UPLC I-Class PLUS liquid chromatograph system (Waters, Milford, MA, United States) and a Waters UPLC Xevo G2-XS QTof mass spectrometer. Each sample was injected into an Acquity UPLC HSS T3 column (1.8 μm, 2.1 mm × 100 mm; Waters, Milford, MA, United States). The column was maintained at 50°C at a flow rate of 0.4 mL/min. The mobile phases consisted of (A) formic acid:water (1:999, v/v) and (B) formic acid:acetonitrile (1:999, v/v). The chromatographic separation was obtained using the following gradient: 0–0.25 min, 98% A and 2% B; 0.25–10 min, 98% A and 2% B; 10–13 min, 2% A and 98% B; 13–13.1 min, 2% A and 98% B; and 13.1–15 min, 98% A and 2% B. The analytical setup was operated in both positive and negative ion modes with the following parameters: capillary voltage of 2,000 V for positive ion mode and −1,500 V for negative ion mode; cone voltage, 30 V; source temperature, 150°C; desolvation temperature, 500°C; desolvation gas flow, 500 L/h; cone gas flow, 50 L/h.

Progenesis QI software (Waters, Milford, MA, United States), the METLIN database (Want, 2005) and the BMK Self-built Database (Biomarker Technologies, China) were used for raw peak exacting, data baseline filtering, calibration of the baseline, peak alignment, deconvolution analysis, peak identification, and integration of the peak area.

Principal component analysis and orthogonal projections to latent structures–discriminate analysis (OPLS-DA) were conducted using the “ropls” R package (version 3.3.2). PCA was used to determine the distribution of the original data. OPLS-DA was performed to obtain group separation and identify the most discriminant variables. A permutation test was used to estimate the robustness and predictive ability of the model. Then, OPLS-DA was used to identify the significantly different metabolites between C. gigas and C. angulata. The variable importance for projection (VIP) values along the predictive component was obtained to explain the results. First, VIP values exceeding 1.0 were selected as changed metabolites. Then, a Student’s t-test was used to evaluate the remaining variables. The metabolites with P-values > 0.05 were discarded. Subsequently, the FC value of each metabolite was calculated to further filter the variables. Taken together, metabolites with VIP values > 1, P < 0.05, and | log₂FC| > 1 were identified as significantly different metabolites (SDMs). The SDMs were enriched with KEGG using the “clusterProfiler” R package (version 3.10.1) (Yu et al., 2012).

Differentially expressed genes (FDR ≤ 0.01 and | log₂FC| > 1) and SDMs (VIP > 1, P < 0.05, and | log₂FC| > 1) were used for the integrative analysis of C. gigas and C. angulata. Pearson correlation coefficients were used to integrate metabolome and transcriptome data. A heatmap was used to show the connections between genes and metabolites.

Growth parameters was measured before we sampled in situ. The shell length, width, and height of each sample were measured using a digital caliper (0.02 mm accuracy). The wet weight of each sample was measured using an electronic balance (0.01 g accuracy).

The glycogen content in the gill tissues of three oysters was determined using a kit for detecting liver and muscle glycogen content (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Soxhlet’s method was used to measure the crude fat content of all tissues excluding the gills from five oysters, and the fatty acid content and composition of these samples were measured using gas chromatography. The protocol was as follows: C19:0 fatty acid methyl ester (Sigma, United States) was added as an internal standard to the sample, and 0.01% butylhydroxytoluene (BHT) methanol solution was added as an antioxidant. After extracting the total fat using dichloromethane:methanol, the sample was dried using high-purity nitrogen. Then, 1 mL 0.5 M KOH methanol solution was added and the mixture was saponified in a water bath at 80°C for 2 h under the protection of nitrogen. After cooling, 1 mL of 14% BF3 methanol solution was added to the sample, and a methyl esterification reaction was carried out in a water bath at 80°C for 1 h. The fatty acid methyl esters were extracted using n-hexane. The sample volume was adjusted to 0.5 mL, and then the fatty acid composition and content was analyzed using an Agilent 7890A gas chromatograph (Agilent Technologies, United States). The chromatographic conditions were as follows: capillary column: DB-FFAP (30 m × 0.32 mm × 0.25 μm); inlet temperature: 220°C; detector temperature: 280°C; column temperature: Program heating 150°C (1 min), 3°C/min, 220°C (33 min).

Ten 6-month-old individuals for each species were used for the normalized oxygen consumption rate measurements. To avoid the impact of periphyton, all of twenty oysters were cleaned and reared temporarily in a 500 L tank for 2 weeks. Averagely, 4 g/m³ commercial spirulina powder was added once per day as a food source (Zhao et al., 2021). An air pump (Super Silent Power Air Pump, China, V-10) was used to ensure a sufficient and stable oxygen supply.

The oysters were placed in sealed bottles (0.35 L) filled with air-saturated seawater. The temperature was regulated using a water bath at 20°C. The seawater was refreshed every 2 h and mixed using a magnetic stirring bar placed beneath the chamber to ensure adequate and uniform distribution of oxygen. The dissolved oxygen levels were measured using a ten-channel optical fiber oxygen meter (OXY-10 ST Prototype, PreSens Precision Sensing GmbH, Germany). The dissolved oxygen content was measured every 3 s using fiber optic probes. Due to the initial violent fluctuation of the dissolved oxygen and unstable water temperature, the data points obtained in the first 30 min were ignored. The subsequent 1,200 data points were used to obtain the fitted line and the slope coefficient K. The normalized oxygen consumption rate of the oysters was calculated as (K× V)/W, where V was the volume of the sealed bottle (L), K was the slope of the fitted line, and W was the dry weight (g) of the oysters (Zhao et al., 2021).

Gills from 15 C. gigas and 15 C. angulata were sampled in situ, frozen in liquid nitrogen, and stored at −80°C for subsequent RNA extraction. RNA extraction was conducted using the Trelief^TM RNAprep FastPure Tissue and Cell Kit (Tsingke Biotechnology, China). cDNA was synthesized using 1 μg of mixed RNA from three oysters, and qRT-PCR was performed using three technical replicates for each cDNA and five biological replicates for each group. The primers that were used are shown in Supplementary Table 1. The relative expression of antioxidant genes [glutathione S-transferase (GST), catalase (CAT), superoxide dismutase (SOD), glutathione (GSH), and xanthine dehydrogenase (XOD)] were determined using the 2^–ΔΔCT method (Livak and Schmittgen, 2001). A comparison between the two species has been described in detail in a previous study (Li et al., 2019). The total antioxidant capacity in the gill tissues of five oysters was measured using a total antioxidant capacity assay kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China).

All statistical analyses were performed with GraphPad Prism version 8.0.2 for Windows. Statistical analyses were conducted after confirming normality of the distributions using the Shapiro–Wilk test and homogeneity of variance using Bartlett’s test. The shell height, shell length, and shell width of C. gigas and C. angulata were compared using the Mann–Whitney U test. The comparison of wet weight, relative gene expression, fatty acid, crude fat, and glycogen content, oxygen consumption rate, and total antioxidant capacity between the two oyster congeners was carried out using an unpaired t-test. Significant differences between groups were marked with “^∗” for P < 0.05, “^∗∗” for P < 0.01, and “^∗∗∗” for P < 0.001.

Based on the sequencing results of the cDNA from the 16 oysters representing the two species, an average of 41.3 and 44.7 million high-quality 150 bp paired-end reads (107.83 Gb reads in total) were mapped to the genomes of C. gigas (Supplementary Table 2). Quality analysis showed that C. gigas and C. angulata had an average guanine-cytosine content of 42.36 and 42.36% of clean reads, a Q30 of 93.52 and 93.43%, and a map rate of 72.63 and 69.63%, respectively (Supplementary Table 2). The data showed a high sequencing quality, indicating that the subsequent transcriptome analysis results were reliable. According to the HCA of gene expression, outlier samples of C. gigas-1 and C. gigas-7 were removed to obtain accurate results for differential genes (Supplementary Figure 1). The transcriptomes of the two oyster species were compared using PCA, which allowed the species to be visually and quantitatively compared in an unbiased format. The percentages of explained values of PC1 and PC2 in the transcriptome were 24.14 and 9.53%, respectively. Figure 1A shows that C. gigas and C. angulata were perfectly separated.

A total of 2,006 DEGs were identified between the two oyster species (FDR ≤ 0.01 and | log₂FC| > 1). The expression levels of these genes are shown in Supplementary Table 3. Compared with C. angulata, a total of 935 genes were upregulated and 1,071 genes were downregulated in C. gigas.

The DEGs were annotated using GO and KEGG. The majority (1,231) of the DEGs were mapped to at least one GO term, and 52 GO terms were significantly enriched. Among these, “cellular process,” “metabolic process,” “membrane,” “membrane part,” “binding,” and “catalytic activity,” were the most strongly represented (Figure 1B).

Pathway enrichment analysis identified significantly enriched metabolic pathways in the DEGs and showed significantly altered pathways in the transcriptome. We mapped all DEGs to KEGG pathways, and 17 significant specific pathways (P < 0.05) were found (Supplementary Table 5 and Figure 1C). Among these, “peroxisome,” “pyrimidine metabolism,” “glutathione metabolism,” “glycine, serine, and threonine metabolism,” “phenylalanine metabolism,” and “alanine, aspartate, and glutamate metabolism” were the most represented pathways. Enrichment of the antioxidant and amino acid metabolism pathways showed that the two oyster species may exhibit differences in stress resistance and protein synthesis.

A total of 3,149 metabolites (2,464 in positive and 685 in negative ionization modes) were identified using the LC-QTOF-MS platform. Gross changes in metabolite composition were determined using PCA and HCA. The percentages of explained values in the metabolome analysis of PC1 and PC2 were 62.52 and 11.3% in positive ionization mode and 47.51 and 18.06% in negative ionization mode, respectively. Based on the PCA (Supplementary Figures 2A,B) and the heatmap analysis (Supplementary Figures 2C,D), there was a slight separation between C. gigas and C. angulata. To maximize the discrimination between the two species, we used OPLS-DA to elucidate the different metabolic patterns. The parameters of the models were R2X = 0.492 and 0.637, R2Y = 0.965 and 0.961, and Q2 = 0.784 and 0.86 for positive and negative ionization modes, respectively (Figures 2A,B). All of them displayed an adequate goodness-of-fit and reliable prediction. In addition, pQ2 = 0.005 and pR2Y = 0.005 for the two modes, which demonstrated the robustness of the models, low overfitting, and low reliability risks (Figures 2C,D). All samples in the score plots were within the 95% Hotelling’s T-squared ellipse, and effective discrimination and separation were found between C. gigas and C. angulata. The OPLS-DA model was determined to be useful for subsequent analyses.

In total, 1,312 (1,094 in positive and 218 in negative ionization modes) SDMs were identified between the two oyster species (P-value ≤ 0.05, VIP > 1, and | log₂FC| > 1) (Supplementary Table 6). In the C. gigas group, 157 and 29 SDMs were significantly upregulated in the positive and negative ionization modes, respectively, and 937 and 189 SDMs were downregulated compared with the C. angulata group (Supplementary Table 6). A total of 189 pathways (135 in positive and 54 in negative ionization modes) were identified when the SDMs between C. gigas and C. angulata were introduced into KEGG (Supplementary Table 7). Based on the P-values and metabolite ratio, we identified that the most significantly enriched pathways were “carbon metabolism,” “purine metabolism,” “bile secretion,” “arachidonic acid metabolism,” “alpha-linolenic acid metabolism,” and “glutathione metabolism” (Figures 3A,B).

Pearson correlation coefficients were used to display the correlation between the transcriptome and the metabolomic data. An overview of the DEGs and metabolites found in C. gigas and C. angulata is shown in Figure 4. The pathways that were enriched between the two groups are shown in Supplementary Figures 3A,B.

We found that the glycolysis pathway was downregulated in C. gigas (Figure 4). Hexokinase (HK) and glucose-6-phosphate isomerase (GPI) were significantly downregulated in C. gigas (Supplementary Table 4). Additionally, acetyl-CoA, an important oxidation product of pyruvic acid that enters the tricarboxylic acid (TCA) cycle (the final product in the glycolysis pathway), was also downregulated in C. gigas (Supplementary Table 6). HK phosphorylates glucose to produce glucose-6-phosphate, which is the first essential step in most glucose metabolism pathways. This phosphorylation event directly associates extramitochondrial glycolysis with intramolecular oxidative phosphorylation. GPI interconverts glucose-6-phosphate and fructose-6-phosphate, which is the second step in glycolysis.

We also found that the gene expression of enoyl-CoA hydratase and 3-hydroxyacyl CoA dehydrogenase (EHHADH) and butanoyl-CoA in the fatty acid degradation pathway were lower in C. gigas than in C. angulata (Figure 4 and Supplementary Tables 4, 6). In addition to the fatty acid degradation pathway, the fatty acid synthesis pathway was upregulated along with high expression of fatty acid synthase (FASN) and high contents of behenic acid (C22:0) and pentadecanoic acid (C15:0) (Figures 4, 5C and Supplementary Table 4). Moreover, stearoyl-CoA desaturase (SCD) involved in fatty acid desaturation pathway was also higher in C. gigas (Figure 4 and Supplementary Table 4).

The metabolome data showed that malate and succinate, which are intermediate products of the TCA cycle, were significantly lower in C. gigas (Figure 4 and Supplementary Table 6). And there was also a lower content of flavin mononucleotide (FMN) which is an important component of electron transport chain in C. gigas (Supplementary Table 6).

According to the transcriptome and metabolome data, we found that many amino acid metabolic pathways were significantly enriched, including glycine, serine, and threonine metabolism; alanine, aspartate, and glutamate metabolism; tyrosine metabolism; cysteine and methionine metabolism; phenylalanine metabolism; and arginine and proline metabolism (Supplementary Figure 3). In the glycine, serine, and threonine metabolism pathways, we observed higher gene expression of betaine-homocysteine S-methyltransferase (BHMT), sarcosine oxidase (PIPOX), choline dehydrogenase (CHDH), and higher content of phosphatidylserine in C. angulata relative to those in C. gigas (Figure 4 and Supplementary Table 4). Interestingly, the expressions of gamma-glutamyl transpeptidase and aminopeptidase N, which catalyze the conversion of glutathione to glycine, were higher in C. angulata (Figure 4 and Supplementary Table 4). In the alanine, aspartate, and glutamate metabolism pathways, we found that aspartate aminotransferase (ASPB) and d-aspartate oxidase (DDO) which mediate aspartic acid into the TCA cycle to supply ATP were upregulated in C. angulata (Figure 4 and Supplementary Table 4). The same result was also identified for cysteine and methionine metabolism; there was higher expression of ASPB and lower expression of cystathionine gamma-lyase in C. angulata, which also mediates the conversion of cysteine and pyruvate (Figure 4 and Supplementary Table 4).

Moreover, we observed that many antioxidant and resistance-related pathways were significantly enriched. In the glutathione metabolism pathway, C. gigas had a higher content of cysteine, greater expressions of glutathione peroxidase and glutathione S-transferase, and a lower content of oxidized glutathione (Figure 4 and Supplementary Tables 4, 6). In the arginine and proline metabolism pathways, we observed that C. gigas had a higher content of gamma aminobutyric acid (GABA) and a higher level of amidase and spermine oxidase (Figure 4 and Supplementary Tables 4, 6), which catalyze arginine and spermine to GABA (Wang et al., 2003). In the tyrosine metabolism pathway, increased levels of aromatic-L-amino-acid/L-tryptophan decarboxylase, catechol O-methyltransferase, and phenylalanine-4-hydroxylase in C. gigas implied that catecholamines such as dopamine and epinephrine were upregulated (Figure 4 and Supplementary Table 4). Moreover, the integration analysis showed that the linoleic acid metabolism pathways were significantly enriched. In the linoleic acid metabolism pathway, we found that linoleic acid and linolenic acid were higher in C. gigas (Figure 4 and Supplementary Table 6).

In morphological data, our results showed that the mean shell height, shell width, and wet weight of C. angulata were significantly higher than those of C. gigas (P < 0.01; Figure 5A and Supplementary Table 8).

In nutritional data, C. gigas had higher glycogen and crude fat contents (Figure 5B and Supplementary Table 9). Additionally, we found that C14:1 (tetradecanoic acid), C15:0 (pentadecanoic acid), C18:1 (oleic acid), C18:2ω6 (linoleic acid), C18:3ω3 (linolenic acid), C20:1ω9 (eicosenoic acid), C20:3ω6 (eicosatrienoic acid), and C22:0 (behenic acid) were significantly higher in C. gigas (Figure 5C). And C20:3ω3 (icosatrienoic acid) was found to be higher in C. angulata (Figure 5C).

In physiological data, the normalized oxygen consumption rate in C. angulata under normal conditions was also higher than that in C. gigas (Figure 6A). Additionally, total antioxidant capacity (Figure 6B) and five anti-oxidation genes (glutathione S-transferase, catalase, superoxide dismutase, glutathione, and xanthine dehydrogenase) expression (1.97–11.44 fold change; P < 0.05; Figure 6C) were higher in C. gigas.

On the basis of the pathway enrichment analysis of DEGs and SDMs, we found that the glycolysis pathway and fatty acids degradation pathway were downregulated in C. gigas. And fatty acids synthesis pathway was upregulated comparatively.

In the glycolysis pathway, HK and GPI were lower in C. gigas. It has been reported that HK expression can be affected by environmental stressors, such as hypoxia (Hall et al., 2009). The glycogen content assays further validated that C. gigas inhibited the glycolysis pathway to increase glycogen content. As for lipids, lower EHHADH expression and butanoyl-CoA content suggested that C. gigas suppressed the fatty acid degradation pathway, which was also confirmed by the crude fat assay. Additionally, higher expression of fatty acid synthase and an increase in the contents of behenic acid and pentadecanoic acid suggested that the fatty acid synthesis pathway was upregulated. Fatty acid synthase is a multifunctional protein, and its main function is to catalyze the synthesis of palmitate from acetyl-CoA and malonyl-CoA into long-chain saturated fatty acids in the presence of NADPH. These results suggest that C. gigas allocates more energy to store glycogen and lipids than does C. angulata.

Interestingly, our data showed that C. gigas had less fatty acid saturation than did C. angulata. It has been widely recognized that membrane lipid unsaturation is closely related to membrane fluidity, which plays an essential role in the cold resistance of organisms (Nishida and Murata, 1996; Guderley, 2004). Higher expression of SCD in C. gigas further illustrated that the unsaturated fatty acid synthesis pathway was upregulated. When on the membrane of the endoplasmic reticulum, SCD is the rate-limiting enzyme that catalyzes the formation of monounsaturated fatty acids from saturated fatty acids. This process is especially important for palmitoyl-CoA (C16:0) and stearoyl-CoA (C18:0), which form palmitoleate (C16:1) and oleic acid (C18:1), respectively, as they are important targets for fatty acid metabolism (Stukey et al., 1990; Wang et al., 2005). Palmitoleate and oleate are the main components of monounsaturated fatty acids in membrane phospholipids, triglycerides, cholesterol esters, and ester waxes, which are closely related to membrane fluidity (Nakamura and Nara, 2004). It has been reported that SCD-1 was significantly induced by low temperature stress in orange-spotted grouper (Epinephelus coioides) (Sun et al., 2019) and large yellow croaker (Larimichthys crocea) (Xu et al., 2015). This is consistent with the relatively low temperatures of the natural habitats in Northern China. Therefore, we speculate that C. gigas upregulated fatty acid unsaturation to maintain membrane fluidity and adapt to the relatively cold environment in Northern China.

Moreover, lower content of malate and succinate reflect that the TCA cycle is downregulated in C. gigas (Figure 6). The TCA cycle is the final oxidative pathway in aerobic metabolism, which links carbohydrates, fats, and amino acids (Akram, 2014). Oxidation of acetyl-CoA by the TCA cycle accounts for 2/3 of total oxygen consumption and ATP production. In addition, low content of flavin mononucleotide (FMN) in C. gigas indicates downregulated oxidative phosphorylation (Supplementary Table 8). FMN is the active form of the flavoprotein prosthetic group in complex I of the mitochondrial electron transport chain and it participates in the electron transfer process during biological oxidation (Liu et al., 2002). The suppression of the TCA cycle and oxidative phosphorylation indicates that C. gigas reduces the ATP supply via aerobic metabolism, which is consistent with previous finding that C. angulata had higher ATP content than C. gigas (Li A. et al., 2021). And the reduction of ATP supply may influence the costs of reproduction, growth, activity, and maintenance (Sokolova et al., 2012).

There were many significantly enriched amino acid metabolism pathways in the integrated analysis. Glycine is an important structural amino acid in the glycine, serine, and threonine metabolism pathways, and its synthesis pathway was upregulated in C. angulata. It has been reported that glycine plays an important role in the metabolism and nutrition of organisms and 80% of all glycine is used for protein synthesis, which means that it has a crucial function in body growth (Razak et al., 2017). Based on a growth improvement experiment, it has also been reported that the glycine, serine, and threonine metabolism pathways are enriched in Pacific white shrimp (Litopenaeus vannamei) (Yu et al., 2021). In C. angulata, the genes that mediate the alanine, aspartate, and glutamate metabolism pathways and the cysteine and methionine metabolism pathways (such as ASPB and DDO) into the TCA cycle to supply ATP were upregulated. Furthermore, we found that C. angulata had higher contents of arginine and histidine than did C. gigas (Supplementary Table 8). The oxygen consumption rate is a physiological and metabolic indicator frequently used to reflect aerobic respiration and energy metabolism levels (Prins et al., 1997), and as a basic physiological activity of energy metabolism, it indicates the physiological state of organisms and the influence of environmental factors (Lannig et al., 2006; Sokolova et al., 2012; Juárez et al., 2018). Therefore, we suggest that the higher normalized oxygen consumption rate in C. angulata than in C. gigas reflects its stronger aerobic metabolism capacity and higher allocation of energy to growth. We infer that C. angulata upregulates amino acid synthesis (such as glycine) and amino acids oxidation (such as aspartate and cysteine) to supply ATP for other life-history components.

Moreover, we found that 40S ribosomal protein S19 (RPS19) and 60S acidic ribosomal protein P1 gene expression was higher in C. angulata (Supplementary Table 4). These proteins participate in the ribosome pathway. Ribosomes, which are vital organelles in all living organisms, are involved in the translation process. The abundance of ribosomes depends on the level of ribosomal protein-encoding mRNAs (RP-mRNAs) and the subsequent translation process (Jesús et al., 2015; Dermit et al., 2019). Therefore, higher expression of RP-mRNAs suggests that more ribosomes are present and that protein synthesis is upregulated. There is experimental evidence to suggest that chitinase negatively regulates biomineralization (Wang et al., 2017; Zhou et al., 2020). Our data showed a significantly lower level of chitinase (Supplementary Table 4) in C. angulata, indicating an upregulated biomineralization process.

In addition, we found some growth performance-related differential metabolites, including uridine diphosphate-glucose in the pyrimidine metabolism pathway, creatinine in the arginine and proline metabolism pathway, and taurocholate and taurocyamine in the taurine and hypotaurine metabolism pathway (Supplementary Table 8). In a study on the effects of temperature on the growth of European seabass (Dicentrarchus labrax), Zhang Z. et al. (2021) found that uridine diphosphate-glucose and taurocholate are potential markers of growth (Zhang Z. et al., 2021). It has been reported that creatinine levels are positively correlated with muscle quality and bone mineral density (Schiffrin et al., 2007; Yamada et al., 2017).

Our data showed that C. angulata has better growth performance than C. gigas, demonstrated by significantly greater shell height, shell weight, and wet weight. Based on our results and those of previous studies, we propose that C. angulata upregulates the amount of available energy by increasing the aerobic metabolism of carbohydrates and lipids and then distributing more energy to synthesize amino acids and proteins.

Glutathione is essential for antioxidant defense, energy metabolism, and regulation of cellular events (Wu et al., 2004; Jones, 2008; Forman et al., 2009). It has been reported that cysteine is the amino acid that is limiting for glutathione synthesis (Lyons et al., 2000; Lu, 2009), and its antioxidant function is related to its sulfhydryl group (Valko et al., 2016). Therefore, a higher cysteine content and higher expression of glutathione peroxidase and glutathione-S-transferase in C. gigas indicated a stronger antioxidant ability. Furthermore, the lower oxidized glutathione content further validated our hypothesis. It has been reported that the glutathione oxidation-reduction product is a signal molecule that participates in the defense responses of organisms (Jones, 2006; Locato et al., 2017). Importantly, antioxidant gene expression and the total antioxidant capacity data were consistent with our results.

Gamma aminobutyric acid is a four-carbon amino acid that is an important inhibitory neurotransmitter in animals (Kinnersley and Turano, 2000; Bouche et al., 2003). Previous studies have mostly focused on the functions of GABA in the growth, development, and resistance of plants (Guo et al., 2021; Wang P. K. et al., 2021). However, it has also been reported that GABA content is higher in more hypoxia-tolerant species such as C. gigas (Haider et al., 2020) and it is helpful for mitigating heat stress in the green-lipped mussel Perna canaliculus Gmelin (Dunphy et al., 2015) and soft corals (Farag et al., 2018). In the tyrosine metabolism pathway of C. gigas, the catecholamine synthesis pathway is upregulated. Under environmental stress, catecholamines and adrenaline can be induced in mammals (Gatzke-Kopp, 2011; Vindas et al., 2017) and Drosophila melanogaster (Adamo, 2017). Linolenic acid (C18:3n-3) is an eighteen-carbon polyunsaturated fatty acid that is a precursor of synthetic eicosanoids, such as prostaglandins and leukotrienes, and it can activate immune function and improve the ability of antioxidants (Geay et al., 2015; Chen et al., 2016).

Previous research has found that C. gigas has more glycogen and lipids than C. angulata. Under environmental stress, glycogen and lipids can be used to quickly supply energy to maintain the basic functions of organisms (Rajamanickam and Muthuswamy, 2008; Ferreira et al., 2011; Wu et al., 2013). Therefore, based on the increased levels of antioxidant-related genes and metabolites, as well as the content of resistant metabolites, we speculate that C. gigas allocates more energy to basal maintenance of some cellular processes related to defense than does C. angulata. Considering the common garden experiment was conducted in the region of the native inhabit of C. gigas (Northern of China), our result was consistent with previous study that C. gigas has higher adaptive capacity than C. angulata in its native inhabit (Li A. et al., 2021).