Integrated application of transcriptomics and metabolomics provides insights into glycogen content regulation in the Pacific oyster Crassostrea gigas (2024)

Experiment materials

Half-sibling families (thirty male oysters individually mated to an egg mix from thirty female oysters) were constructed and cultured on a farm in Jiaonan, Qingdao, China. No specific permissions were required for the present study and all experiments were conducted with approval from the Experimental Animal Ethics Committee, Institute of Oceanology, Chinese Academy of Sciences, China. At the age of sixteen months, 185 individuals were selected randomly and used for metabolome and transcriptome analysis. The adductor muscle was dissected from the oyster, and the remaining tissues were ground in the presence of liquid nitrogen and mixed together to form a powdered tissue extract which was then stored at −80°C before analysis. This powdered tissue extract was then used for the measurement of glycogen content, as well as for transcriptome and metabolome analysis. Thirty oysters, out of the initial 185 individuals, consisting of fifteen individuals with the highest glycogen content, and fifteen individuals with the lowest glycogen content, were selected for transcriptome analysis. The top ten oysters with the highest glycogen content and the ten oysters with the lowest glycogen content, were used for metabolome analysis. Shell height, length, and width, total weight, and tissue weights were measured before dissection.

Glycogen content assay

The glycogen content of the entire tissues without the adductor muscle, was detected using a kit for detecting liver and muscle glycogen content (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). The procedure used was as follows: after tissues were ground into a powder in the presence of liquid nitrogen, a 0.50 μg sample was added to a tube containing alkaline liquor. The samples were incubated at 100°C for 20min in a water bath. After that, the hydrolysate was diluted 16-fold by the addition of distilled water. Subsequently, 2mL of color reagent was added to the diluted hydrolysate and the samples were incubated for 5min at a 100°C in a water bath. Finally, the OD value of each sample was measured at 620nm using a Multiscan Spectrum with a path length of 1cm. Set the blank and standard for each group (Additional file 1: Table S1). The glycogen content was calculated according to the following formula: where 1.11 is the conversion coefficient of glucose content to glycogen content in this method. 0.01 is the glycogen content in the standard group. The coefficients 20 and 10 respectively represent the dilution ratio of the sample before and during detection.

* The OD value of each tube was measured at 620nm wavelength.

RNA sequencing data analysis

The same tissues used for glycogen content detection were also used for RNA extraction (i.e. the total tissue lacking the adductor muscle). To increase the accuracy, we increased the sample number in the transcriptome analysis to fifteen. RNA from the thirty oyster tissues (i.e. the fifteen high glycogen and fifteen low glycogen oysters) were exacted using TRIzol. Each of the thirty samples was used to generate an independent library. The thirty RNA-seq paired-end libraries were generated using NEBNext® Ultra™ RNA LibraryPrep Kit for Illumina® (NEB, USA) following the instruction’s procedure and index codes were added to attribute sequences to each sample. Briefly, poly-dT oligo-attached magnetic beads were used to extract mRNA from total RNA. mRNA fragmentation was performed under elevated temperature with divalent cations in NEBNext First Strand Synthesis Reaction Buffer. First strand cDNA synthesis was carried out using hexamer primers and M-MuLV reverse transcriptase (RNase H-) while the second strand was synthesis by DNA polymerase I and RNase H subsequently.

The remaining overhangs were converted into blunt ends using a combination of exonuclease/polymerase. To prepare for hybridization, NEBnext adapters with a hairpin loop structure were ligated to DNA fragments with the adenylated 3′ ends. After size selection (150~200bp) of the ligation products, the ligated cDNA fragments that contained the adapter sequences were enhanced via PCR using PCR primers and Index (X) Primer. At the end, purification of PCR products were performed (AMPure XP system) and Agilent Bioanalyzer 2100 system were used for library quality assessment. Sequencing of libraries was performed using an Illumina Hiseq 2000 platform and 150bp paired-end raw reads were generated. After removing adapter sequences, poly-N sequences, and low quality reads, clean reads were obtained. Meanwhile, the Q20, Q30, and GC content of the clean data were calculated.

Read mapping and gene expression calculation

The HISAT2 v2.0.4 package [24] was used to map RNA-seq reads to the oyster genome (GenBank accession No., GCA_000297895.1) using default settings. The SAM files produced by HISAT2 were converted to BAM files and sorted using SAMtools v.1.3 [25]. The StringTie v1.3.3 package [26] was used to process the read alignments and the reference annotation with a parameter setting of “-eB”. Using this input, StringTie estimates abundance and creates a new transcript table for input to Ballgown. Gene expression levels were assessed using FPKM calculated using the R package Ballgown [27]. The pipeline for transcriptome analysis is shown in Additional file 2: Figure S1.

Analysis of differentially expressed genes

Calculation of FPKM of each gene was performed by Cufflinks based on the length of gene read counts mapped to genome. Differential expression analysis was performed between high and low glycogen content groups using the DEGseq R package (1.10.1) [28]. Benjamini and Hochberg’s approach were used to control the false discovery rate (FDR) by resulting P-values adjustment. Genes with an adjusted P-value <0.001 found by DEGSeq were selected as differentially expressed. Principal component analysis (PCA) was performed using the SIMCA14 software package (Umetrics, Umea, Sweden).

Gene Ontology (GO) as an international standardized gene-function classification system consists of terms that provide a more global representation of gene functions using a controlled vocabulary. We selected 742 genes that had gene names annotated in the UniProt database from the 1888 differentially-expressed genes and transformed these gene names to protein names according to ID mapping on the UniProt website. GO enrichment was then analyzed using DAVID 6.8 Beta (david.ncifcrf.gov) and the results were classified with GOTERM_BP, GOTERM_CC, and GOTERM_MF. Significantly enriched metabolic and signal transduction pathways represented by DEGs were determined using pathway-enrichment analysis, compared with the whole genome background. KOBAS software were used to test the statistical enrichment of DEGs in KEGG pathways. A Q value ≤0.001 was considered as significantly enriched pathways among DEGs.

Metabolite extraction and detection

0.05g of each oyster tissue was extracted with 0.4mL methanol–chloroform (3/1, v/v) and 20μL of 2-Chloro-L-phenylalanine was added as an internal standard. The extract was hom*ogenized using a ball mill for 4min at 45Hz, followed by centrifugation at 13,000rpm for 15min at 4°C. Supernatant fluid (about 0.35ml) was transferred to a new 2ml glass vial. The extracts were dried under the vacuum concentrator without heating, and added 60μL of methoxyamine hydrochloride (20mg/mL in pyridine) to the dried metabolites, and then incubated for 30min at 80°C, after mixing and sealing. BSFTA reagent (80μL in 1% TMCS, v/v) was then added to the samples and they were incubated for 2h at 70°C. Add FAMEs (5μL), consisting of a standard mix of fatty acid methyl esters (C8–C16: 1mg/mL; and C18–C24: 0.5mg/mL in chloroform), to samples. Prior to GC–time-of-flight (TOF)–MS analysis, the samples were cooled to room temperature and fully mixed. GC-TOF/MS is a method that is capable of detecting low abundance metabolites, and has the ability to identify secondary metabolites as well as metabolic intermediates [29].

The detection of the metabolites were referred to methods of Hairui Wang et al. in 2017 [30]. In brief, GC–TOF–MS analysis was carried out utilizing an Agilent 7890 gas chromatograph system with the addition of a Pegasus 4D time-of-flight mass spectrometer. A DB-5MS capillary column were applied in this system. The column with 30m×250μm inner diameter, 0.25μm film thickness were coated with 5% diphenyl cross-linked with 95% dimethylpolysiloxane (J&W Scientific, Folsom, CA, USA). Splitless mode were used for injection with Helium as the carrier gas. The front inlet purge and gas flow rate were 3mL/min and 1mL/min respectively. The temperature program was 80°C for 1min, followed by raised to 290°C at a rate of 10°C per min and finally held at 290°C for 5min. Temperature for injection, ion source, and transfer line were 280, 220, and 295°C respectively, with the −70eV energy in electron impact mode. Mass spectrometry data were collected in full-scan mode with a range of 30–600m/z at a rate of 20 spectra per second after a solvent delay of 7min [30].

Chroma TOF 4.3X software from LECO Corporation and the LECO-Fiehn Rtx5 database were used for raw peak exacting, data baseline filtering and calibration of the baseline, peak alignment, deconvolution analysis, peak identification, and integration of the peak area. Peak identification were performed by retention time index method with the 5000 RI tolerance.

Multivariate statistical analysis of metabolites

PCA (principal component analysis) and OPLS-DA (orthogonal projections to latent structures–discriminate analysis) were performed by SIMCA14 software package (Umetrics, Umea, Sweden) with the results of three-dimensional data including the peak numbers (which is a metabolite marker), sample names, and normalized peak areas. PCA was applied to distribute the original data while OPLS-DA were performed for obtaining of effective group separation and better understanding of factors responsible for the classification. The parameters of classification were R2Y=0.961 and Q2Y=−0.113 which indicated adequate goodness-of-fit and predictive ability.

Identification of significantly different metabolites and pathways between high and low glycogen content oysters

As described above, OPLS-DA was used to identify significantly different metabolites between the high- and low-glycogen content groups. To clarify this analysis, the VIP (variable importance for the projection) values along the predictive component were obtained. VIP values exceeding 1.0 were first selected as being changed metabolites. T test (Student’s t test) were used to the evaluation of remaining variables. Discard the variables with P>0.05 between the two comparison groups. The fold change (FC) value of each metabolite was calculated by comparing the mean value of the peak area obtained from the high-glycogen content group to that from the low-glycogen content group. The metabolites were blasted against the metabolite database of Caenorhabditis elegans, which is evolutionarily close to C. gigas. Search the online databases involving the Kyoto Encyclopedia of Genes and Genomes (KEGG), and NIST (http://www.nist.gov/index.html) to further identify and confirm differential metabolites. Each differential metabolite and related pathways were cross-listed in KEGG. Additionally, the top enriched pathways were selected and finally constructed in accordance with the potential functional analysis.

Integrative analysis of metabolome and transcriptome

Pearson correlation coefficients were calculated for metabolome and transcriptome data integration. MetScape 3 was utilized to visualize and clarfy the metabolomic and transcriptome data. This allow us to build networks of genes and metabolites. It also contribute to visualize metabolites difference and identify enriched pathways from transcriptome data.

Metabolome analysis of the oysters with high glycogen and low glycogen

SPSS version 16.0 was used to build a descriptive statistical analysis of the glycogen content, as well as to construct a correlation analysis between glycogen content and shell height, shell length, shell width, total weight, and tissue weight. The glycogen content of the experimental oysters ranged from 3.1 to 35.7mg/g, with a mean±std. value of 16.3±5.3. Glycogen content was normally distributed (P>0.05), with more than ten-fold degree of variation. (Additional file 3: Figure S2). A correlation analysis showed that oyster glycogen content had a positive correlation with both total weight and tissue weight (p<0.01) (Additional file 4: Table S2).

Representative chromatograms of the oyster tissue samples analyzed by GC-MS are shown (Additional file 5: Figure S3). To compare the metabolite composition of the two oyster groups having either high- or low-glycogen content, both hierarchical cluster analysis (HCA) and principal component analysis (PCA) models were tested using the data obtained from GC-TOF-MS analysis. The percentage of explained value in the metabolome analysis of PC1 and PC2 was 13.2% and 9.5% respectively. Based on the PCA (Fig. ​(Fig.1a)1a) and the heat map analysis (Fig. ​(Fig.1b),1b), there was an obvious separation between the high- and low-glycogen content oysters.

In total, seventy-five differentially-abundant metabolites were identified from 767 peaks in the chromatograms, including amino acids, fatty acids, and sugars (Additional file 6: Table S3). In the high-glycogen content group, fifty-six metabolites were significantly up-regulated and nineteen were significantly down-regulated compared with the low-glycogen content group. Specific metabolite information and the concentrations in the high- and low- glycogen groups are shown in Additional file 7: Figure S4. The metabolome view map revealed that thirty pathways (Additional file 8: Table S4 and Fig. ​Fig.1c)1c) were enriched based on the seventy-five differentially-abundant metabolites originally identified between the two groups. However, only four of these thirty pathways displayed significant enrichment (P<0.05) namely, “alanine, aspartate and glutamate metabolism”(map00250), “d-glutamine and d-glutamate metabolism”(map00471), “arginine and proline metabolism”(map00330), and “citrate cycle (TCA cycle)”(map00020) (Fig. ​(Fig.1d).1d). These enriched pathways were related to amino acid metabolism and glycogen metabolism. Interestingly, five fatty acids (palmitic acid, myristic acid, arachidonic acid, pelargonic acid, and pentadecanoic acid) and nine amino acids (alanine, aspartic acid, asparagine, glutamic acid, glutamine, histidine, methionine, isoleucine, and valine) were found to be abundant in high-glycogen content oysters (Fig. 2a, b). A correlation analysis conducted using SPSS version 16.0 showed that glycogen content had a positive correlation with the individual amino acid and free fatty acid levels (Fig. ​(Fig.2c).2c). Besides this, pyruvic acid, oxaloacetic acid, and alpha-ketoglutaric acid, which are involved in the TCA cycle, were also found to be abundant in high-glycogen content oysters (Additional file 7: Figure S4).

Transcriptome analysis of oysters with high-glycogen and low-glycogen contents

cDNAs prepared from the thirty oysters representing the two group were sequenced individually using an Illumina Hiseq2000 platform. An average of 11.9 and 12.5 million high-quality 150bp paired-end reads (107GB reads in total) were mapped to the C. gigas genome from the fifteen high- and low-glycogen content oysters, respectively (Additional file 9: Table S5). A total of 25,911 genes had reads that mapped to the C. gigas genome and the average map rate was 73.55%. The transcriptomes of the two groups of oysters were compared using a principle component analysis (PCA), which allowed an unbiased analysis in a format in which the groups could be visually and quantitatively compared. The percentage of explained value of PC1, PC2, and PC3 in the transcriptome was 10.1%, 9.7%, and 7.2% respectively. Figure ​Figure3a3a shows that the high- and low-glycogen content oysters could be perfectly separated.

A total of 1888 differentially-expressed genes (DEGs) were identified between the two glycogen content groups (geometric algorithm, P-value <0.001). The expression level of the genes mentioned are shown in Additional file 10: Table S6. Compared to the low-glycogen content group, a total of 981 genes were up-regulated and 907 genes were down-regulated in the high-glycogen content group.

The DEGs were annotated according to Gene Ontology (GO) and KEGG. The majority (732) of the DEGs mapped to at least one GO term, and 181 GO terms were significantly enriched. Among these, “biological process”, “response to stress”, “single-organism metabolic process” and “negative regulation of metabolic process” were the most highly represented (Fig. ​(Fig.3b3b).

We used the method described by Kobas [31] to find pathways in the KEGG database that were significantly altered in the transcriptome. We mapped all DEGs to KEGG pathways and twenty-nine significant specific pathways (P<0.05) were represented (Additional file 11: Table S7 and Fig. ​Fig.3c).3c). Among these twenty-nine pathways “Ribosome”, “Phagosome”, “Glycolysis/Gluconeogenesis” were the three most-represented pathways (Fig. 3c and d). Glycolysis and gluconeogenesis are both important pathways in glucose metabolism. The rate-limiting enzymes hexokinase (HK) and pyruvate kinase (PK) in glycolysis, and phosphoenolpyruvate carboxykinase (PEPCK) in gluconeogenesis, were highly expressed in high-glycogen oysters (Additional file 10: Table S6).

Key pathways related to oyster glycogen content difference analysis by metabolome and transcriptome

MetScape 3 was utilized to visualize and clarfy the metabolomic and transcriptome data [32]. This allow us to build networks of genes and metabolites. It also contribute to visualize metabolites difference and identify enriched pathways from transcriptome data. An overview of the differentially expressed genes and metabolites found between high glycogen oysters and low glycogen oysters is shown in Fig. ​Fig.4.4. The pathways enriched between the high-glycogen content group and the low-glycogen content group are shown in Additional file 12: Figure S5. Specific pathways, metabolites, and gene information are shown in Additional file 13: Table S8. A total of twenty-seven signaling pathways played important roles between the two groups, including thirty-nine key metabolites and seventeen genes known to be involved in “Glycolysis and Gluconeogenesis”, “the TCA cycle”, “fatty de novo acid biosynthesis and beta-oxidation”, and “amino acid metabolism”. For instance, in the TCA cycle, abundant levels of oxaloacetic acid, alpha-ketoglutaric acid, and pyruvic acid (Fig. ​(Fig.1c1c and Fig. ​Fig.4)4) were detected in high-glycogen oysters together with high levels of malate dehydrogenase (MDH) and pyruvate carboxylase (PYC) (Additional file 10: Table S6).

In the above analysis, we found that three differentially-expressed genes and metabolites were directly related and were a cross validation (Fig. 5a and b). Oleoyl-ACP hydrolase (OLAH) and palmitoyl-protein thioesterase (PPT2) which are both hydrolases that participate in fatty acid biosynthesis are highly expressed in high-glycogen content oysters (Additional file 10: Table S6). Consistent with this, their products palmitic acid (hexadecanoic acid) and myristic acid (tetradecanoic acid) were also abundant in high-glycogen content oysters. PPT is a thioesterase that normally functions to remove long-chain fatty acids like palmitate from modified cysteine residues in proteins [33]. PPT releases hexadecanoic acid according to the R01274 reaction. Oleoyl-ACP hydrolase (OLAH) hydrolyzes acyl-ACP into the free fatty acid and ACP, thereby releasing the fatty acid from the fatty acid synthase complex [34]. OLAH releases hexadecanoic acid and tetradecanoic acid according to reactions R01706 and RE1576. The expression of these DEGs and the metabolites are shown in Fig. ​Fig.5c5c.

Notably, three pathways that share a common enzyme, referred to as glutathione peroxidase 5 (GPX5), were identified based on an integrated analysis of both metabolome and transcriptome data (Additional file 14: Figure S6). These are the urea cycle, along with the metabolism of arginine, proline, glutamate, aspartate and asparagine, arachidonic acid metabolism, and linoleate metabolism. GPX5, which was highly expressed in high-glycogen content oysters, belongs to the glutathione peroxidase family and has the ability to reduce hydroperoxides [35].

The glycogen content is close related to the fatty acid degradation capacity

There is experimental evidence to suggest that a high FFA concentration leads to the accumulation of glycogen through the suppression of glycolysis and TCA cycle [16, 40] but the molecular mechanism is not understood. Based on the data presented here, we hypothesize that the high glycogen content in oysters may be related to a decrease in expression of CPT2 gene and a decrease in the conversion of free fatty acids to acyl-CoAs. In our study, the “de novo fatty acid biosynthesis” pathway was highly expressed in high-glycogen content oysters along with abundant levels of myristic acid and palmitic acid. There was also high expression of OLAH and PPT, which participate in the release of fatty acids from the fatty acid synthase complex [33, 34]. In contrast to this, in low glycogen-content oysters, we showed that the rate-limiting enzyme of fatty acid degradation, namely carnitine O-palmitoyltransferase 2 (CPT2), was highly expressed and there was low expression of fatty acids like palmitic acid and myristic acid (Fig. ​(Fig.1c).1c). CPT I and CPT II (Carnitine palmitoyltransferases I and II), together with the acyl-carnitine translocator, mediate the transfer of acyl-groups into mitochondria [41], and are the rate-limiting enzymes in fatty acid degradation and mitochondrial β-oxidation [42]. It has been reported in rats that a CPT inhibitor could inhibit FFA oxidation and thus affect glycogen levels [43]. This study is consistent with our result that FFA concentration was related to glycogen accumulation in the Pacific oyster. High levels of CPT expression are indicative of robust fatty acid degradation that would result in a decrease in the levels of free fatty acids along with an accumulation of acyl-CoA in low-glycogen oysters. This accumulation of long chain fatty acyl-CoA may inhibit glycogen synthase activity suppressing glycogen synthesis in low-glycogen oysters [44, 45]. In short, a high level of CPT2 expression increases flux through the fatty acid degradation pathway, causes a lowering of the free fatty acid levels and a reciprocal increase in long chain fatty acyl-CoA levels that then suppress glycogen accumulation.

Amino acid content and conversion efficiency to glucose play crucial role in glycogen content

In our study, we observed an increase in the abundance of glucogenic amino acids, including alanine, aspartate, glutamic acid, as well as others, in high-glycogen oysters. There was also an increase in the level of pyruvic acid, as well as increased expression of PEPCK which catalyzes the first committed step in gluconeogenesis [46] in these same oysters. In rats, Nomura et al. found that Co²⁺-induced PEPCK suppression causes a reduction in glycogen content [47]. It has also been reported that in Litopenaeus vannamei, the gluconeogenic capacity of the hepatopancreas can be modulated by PEPCK activity which affects the glycogen content in shrimp [48]. It has also been reported that in humans, the glycogen content in liver was dependent on the activity of PEPCK [49]. These studies further validate our data whereby high levels of PEPCK expression correlate with high glycogen levels in the Pacific oyster. Taken together, we speculate that this increase in abundance of glucogenic amino acids (i.e. gluconeogenic substrates) coupled with an increase in gluconeogenic capacity leads to increased de novo production of glucose and ultimately to increased glycogen synthesis. In this model, the glucogenic amino acids are first transformed into the corresponding ketonic acid via a deamination reaction. Through the TCA cycle, these ketonic acids are transformed into pyruvic acid, which is then converted to oxaloacetate [50] and then to PEP by the action of PEPCK [46], leading to the generation of glucose. The high levels of PEPCK expression and the high pyruvic acid content imply a high level of gluconeogenesis in the high-glycogen content oysters. Together with the high levels of glucogenic amino acids, this will result in a high level of de novo glucose production leading to increased glycogen storage. A recent study that examined the effect of high protein feeding in rats also supported our hypothesis. Wistar rats fed a high protein diet also increased their glycogen content because a significant portion of the dietary amino acids was channeled to glycogen synthesis [17].

Oysters with higher glycogen content exhibited increased energy metabolism

Based on an analysis of key metabolic pathway using the seventy-five differentially-abundant metabolites and the 1888 differentially-expressed genes, we observed that both the glycolytic pathway and the TCA cycle were enhanced in oysters with high-glycogen content, which implied increased energy metabolism occurs in the high-glycogen oysters (Fig. ​(Fig.4).4). HK and PK, which are both rate limiting enzymes in glycolysis and catalyze irreversible reactions, were found to be highly expressed in high-glycogen content oysters. Up-regulation of the TCA cycle was inferred by the increased abundance of oxaloacetic acid, alpha-ketoglutaric acid, and pyruvic acid, as well as increases in the expression of MDH and PYC. MDH and PYC play crucial roles in the citrate cycle wherein MDH catalyzes the interconversion of oxaloacetate and malate [51], while PYC catalyzes the formation of oxaloacetate from pyruvate with concomitant ATP cleavage [52]. Since glycogen is a glucose polysaccharide that serves as a form of energy storage [3], these data imply that there is an concurrent upregulation of the mobilization of glycogen through glycogen phosphorylase, with the resultant liberated glucose then completely metabolized to CO₂ and H₂O through glycolysis and oxidative metabolism (TCA cycle), along with the release of the energy.

The high level of expression of the HK and PK genes involved in glycolysis and the MDH and PYC genes that participate in the TCA cycle in high-glycogen oysters reveals that these oysters have a higher energy metabolism compared to low-glycogen content oysters and implies that they have a greater resistance to stress.

Higher glycogen content oyster manifests as a higher anti-adversity ability

Interestingly, we also showed that the antioxidant enzyme GPX5 was highly expressed in the high-glycogen content oysters. The enzyme GPX5 was found to be in common with three pathways that were identified based on an integrated analysis of both the metabolome and the transcriptome. Glutathione peroxidase (GPX) is an antioxidant that has the ability to reduce hydroperoxides to alcohols and hydrogen peroxide to water in a glutathione-dependent reaction [35]. The high expression of GPX5 implies there is an increased antioxidant capacity, which is expected to contribute to an ability of the high-glycogen oysters to cope with stress. Glycogen may therefore play an important role in antioxidant defense and stress resistance. It has been reported that glycogen levels are massively reduced under hypo-osmotic stress [53], and that when the glycogen is exhausted, the osmotic imbalance caused by salinity stress is more severe [54]. It has also been reported that glycogen could be a potential biomarker for Cu-induced toxic effects in the Pacific oyster [55]. Although it has been widely recognized that glycogen participates in stress resistance in the oyster, the molecular mechanism behind this has not been identified. Making use of integrated analysis of transcriptome and metabolome, we suggest that GPX5 may be an important enzyme in the mechanism by which glycogen contributes to stress-resistance in oysters.

The authors thank all members of the laboratory for valuable discussions.

Ethics approval and consent to participate

The Pacific oysters, C. gigas, used in this study were marine-cultured animal and were collected from a local farm in Qingdao, China. No specific permissions are required for the collection of oysters. All the experiments were conducted with approval from Experimental Animal Ethics Committee, Institute of Oceanology, Chinese Academy of Sciences, China and did not involve any endangered or protected species.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

The online version of this article (doi:10.1186/s12864-017-4069-8) contains supplementary material, which is available to authorized users.

1. Samain JF. Review and perspectives of physiological mechanisms underlying genetically-based resistance of the Pacific oyster Crassostrea gigas to summer mortality. Aquat Living Resour. 2011;24(3):227–236. doi:10.1051/alr/2011144. [CrossRef] [Google Scholar]

2. Maguire GB, Gardner NC, Nell JA, Kent GN, Kent AS. Studies on triploid oysters in Australia .2. Growth, condition index gonad area, and glycogen content of triploid and diploid Pacific oysters, Crassostrea gigas, from oyster leases in Tasmania, Australia. Aquaculture. 1995;137(1–4):357. doi:10.1016/0044-8486(96)83594-0. [CrossRef] [Google Scholar]

3. Roach PJ, Depaoli-Roach AA, Hurley TD, Tagliabracci VS. Glycogen and its metabolism: some new developments and old themes. Biochem J. 2012;441:763–787. doi:10.1042/BJ20111416. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

4. Berri C, Debut M, Sante-Lhoutellier V, Arnould C, Boutten B, Sellier N, Baeza E, Jehl N, Jego Y, Duclos MJ, et al. Variations in chicken breast meat quality: implications of struggle and muscle glycogen content at death. Brit Poultry Sci. 2005;46(5):572–579. doi:10.1080/00071660500303099. [PubMed] [CrossRef] [Google Scholar]

5. Chaosap C, Parr T, Wiseman J. Effect of compensatory growth on forms of glycogen, postmortem proteolysis, and meat quality in pigs. J Anim Sci. 2011;89(7):2231–2242. doi:10.2527/jas.2010-2953. [PubMed] [CrossRef] [Google Scholar]

6. Gierus M, Rocha JB. Forage substitution in a grain-based diet affects pH and glycogen content of semimembranosus and semitendinosus rabbit muscles. J Anim Sci. 1997;75(11):2920–2923. doi:10.2527/1997.75112920x. [PubMed] [CrossRef] [Google Scholar]

7. Lekatz LA, Ward MA, Borowicz PP, Taylor JB, Redmer DA, Grazul-Bilska AT, Reynolds LP, Caton JS, Vonnahme KA. Cotyledonary responses to maternal selenium and dietary restriction may influence alterations in fetal weight and fetal liver glycogen in sheep. Anim Reprod Sci. 2010;117(3–4):216–225. doi:10.1016/j.anireprosci.2009.05.009. [PubMed] [CrossRef] [Google Scholar]

8. Hata K, Hata M, Matsuda K. Purification and properties of glycogen-Phosphorylase from the adductor muscle of the oyster, Crassostrea gigas. Comp Biochem Phys B. 1993;105(3–4):481–486. doi:10.1016/0305-0491(93)90077-I. [CrossRef] [Google Scholar]

9. Bacca H, Huvet A, Fabioux C, Daniel JY, Delaporte M, Pouvreau S, Van Wormhoudt A, Moal J. Molecular cloning and seasonal expression of oyster glycogen phosphorylase and glycogen synthase genes. Comp Biochem Phys B. 2005;140(4):635–646. doi:10.1016/j.cbpc.2005.01.005. [PubMed] [CrossRef] [Google Scholar]

10. Zeng Z, Ni J, Ke C. Expression of glycogen synthase (GYS) and glycogen synthase kinase 3beta (GSK3beta) of the Fujian oyster, Crassostrea angulata, in relation to glycogen content in gonad development. Comp Biochem Phys B. 2013;166(3–4):203–214. doi:10.1016/j.cbpb.2013.09.003. [PubMed] [CrossRef] [Google Scholar]

11. She Z, Li L, Qi H, Song K, Que H, Zhang G. Candidate gene polymorphisms and their association with glycogen content in the Pacific oyster Crassostrea Gigas. PLoS One. 2015;10(5):e0124401. doi:10.1371/journal.pone.0124401. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

12. Liu S, Li Q, Yu H, Kong L. Relationship between single nucleotide polymorphism of glycogen synthase gene of Pacific oyster Crassostrea Gigas and its glycogen content. J Ocean Univ China. 2017;16(1):168–174. doi:10.1007/s11802-017-3141-y. [CrossRef] [Google Scholar]

13. Zhong XX, Li Q, Guo X, Yu H, Kong LF. QTL mapping for glycogen content and shell pigmentation in the Pacific oyster Crassostrea Gigas using microsatellites and SNPs. Aquacult Int. 2014;22(6):1877–1889. doi:10.1007/s10499-014-9789-z. [CrossRef] [Google Scholar]

14. Ma J, Yang J, Zhou L, Ren J, Liu X, Zhang H, Yang B, Zhang Z, Ma H, Xie X, et al. A splice mutation in the PHKG1 gene causes high glycogen content and low meat quality in pig skeletal muscle. PLoS Genet. 2014;10(10):e1004710. doi:10.1371/journal.pgen.1004710. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

15. Duan YY, Ma JW, Yuan F, Huang LB, Yang KX, Xie JP, Wu GZ, Huang LS. Genome-wide identification of quantitative trait loci for pork temperature, pH decline, and glycolytic potential in a large-scale white Duroc x Chinese Erhualian resource population. J Anim Sci. 2009;87(1):9–16. doi:10.2527/jas.2008-1128. [PubMed] [CrossRef] [Google Scholar]

16. Randle PJ, Garland PB, Hales CN, Newsholme EA. The glucose fatty-acid cycle. Its role in insulin sensitivity and the metabolic disturbances of diabetes mellitus. Lancet. 1963;1(7285):785–789. doi:10.1016/S0140-6736(63)91500-9. [PubMed] [CrossRef] [Google Scholar]

17. Fromentin C, Azzout-Marniche D, Tomé D, Even P, Luengo-Guyonnot C, Piedcoq J, Fromentin G, Gaudichon C. The postprandial use of dietary amino acids as an energy substrate is delayed after the deamination process in rats adapted for 2 weeks to a high protein diet. Amino Acids. 2011;40(5):1461–1472. doi:10.1007/s00726-010-0756-3. [PubMed] [CrossRef] [Google Scholar]

18. Rondon R, Akcha F, Alonso P, Menard D, Rouxel J, Montagnani C, Mitta G, Cosseau C, Grunau C. Transcriptional changes in Crassostrea gigas oyster spat following a parental exposure to the herbicide diuron. Aquat Toxicol. 2016;175:47–55. doi:10.1016/j.aquatox.2016.03.007. [PubMed] [CrossRef] [Google Scholar]

19. Lim HJ, Kim BM, Hwang IJ, Lee JS, Choi IY, Kim YJ, Rhee JS. Thermal stress induces a distinct transcriptome profile in the Pacific oyster Crassostrea Gigas. Comp Biochem Physiol Part D Genomics Proteomics. 2016;19:62–70. doi:10.1016/j.cbd.2016.06.006. [PubMed] [CrossRef] [Google Scholar]

20. Zhang L, Li L, Zhu Y, Zhang G, Guo X. Transcriptome analysis reveals a rich gene set related to innate immunity in the eastern oyster (Crassostrea Virginica) Mar Biotechnol (NY) 2014;16(1):17–33. doi:10.1007/s10126-013-9526-z. [PubMed] [CrossRef] [Google Scholar]

21. Zhao S, Zhao J, Bu D, Sun P, Wang J, Dong Z. Metabolomics analysis reveals large effect of roughage types on rumen microbial metabolic profile in dairy cows. Lett Appl Microbiol. 2014;59(1):79–85. doi:10.1111/lam.12247. [PubMed] [CrossRef] [Google Scholar]

22. Gracey AY, Connor K. Transcriptional and metabolomic characterization of spontaneous metabolic cycles in Mytilus Californianus under subtidal conditions. Mar Genomics. 2016;30:35–41. doi:10.1016/j.margen.2016.07.004. [PubMed] [CrossRef] [Google Scholar]

23. Xu G, Du F, Li Y, Nie Z, Xu P. Integrated application of transcriptomics and metabolomics yields insights into population-asynchronous ovary development in Coilia Nasus. Sci Rep. 2016;6:31835. doi:10.1038/srep31835. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

24. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–360. doi:10.1038/nmeth.3317. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

25. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Proc GPD The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–2079. doi:10.1093/bioinformatics/btp352. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

26. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–295. doi:10.1038/nbt.3122. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

27. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–1667. doi:10.1038/nprot.2016.095. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

28. Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–138. doi:10.1093/bioinformatics/btp612. [PubMed] [CrossRef] [Google Scholar]

29. Scano P, Murgia A, Pirisi FM, Caboni P. A gas chromatography-mass spectrometry-based metabolomic approach for the characterization of goat milk compared with cow milk. J Dairy Sci. 2014;97(10):6057–6066. doi:10.3168/jds.2014-8247. [PubMed] [CrossRef] [Google Scholar]

30. Wang HR, Zhong HJ, Hou R, Ayala J, Liu GM, Yuan SB, Yan Z, Zhang WP, Liu YL, Cai KL, et al. A diet diverse in bamboo parts is important for Giant panda (Ailuropoda Melanoleuca) metabolism and health. Sci Rep-Uk. 2017;7(1):3377. doi:10.1038/s41598-017-03216-8. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

31. Wu J, Mao X, Cai T, Luo J, Wei L. KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 2006;34(Web Server issue):W720–W724. doi:10.1093/nar/gkl167. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

32. Karnovsky A, Weymouth T, Hull T, Tarcea VG, Scardoni G, Laudanna C, Sartor MA, Stringer KA, Jagadish HV, Burant C, et al. Metscape 2 bioinformatics tool for the analysis and visualization of metabolomics and gene expression data. Bioinformatics. 2012;28(3):373–380. doi:10.1093/bioinformatics/btr661. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

33. Segal-Salto M, Sapir T, Reiner O. Reversible Cysteine Acylation regulates the activity of human Palmitoyl-protein Thioesterase 1 (PPT1) PLoS One. 2016;11(1):e0146466. doi:10.1371/journal.pone.0146466. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

34. Imai H, Nishida I, Murata N. Acyl-(acyl-carrier protein) hydrolase from squash cotyledons specific to long-chain fatty acids: purification and characterization. Plant Mol Biol. 1992;20(2):199–206. doi:10.1007/BF00014488. [PubMed] [CrossRef] [Google Scholar]

35. Vieira RC, Silva CMS, de Araujo MB, Garcia A, Voltarelli VA, dos Reis AD, Voltarelli FA. Aerobic swimming training increases the Actitivy of antioxidant enzymes and the glycogen content in the skeletal muscle of rats. Rev Bras Med Esporte. 2013;19(3):204–208. doi:10.1590/S1517-86922013000300012. [CrossRef] [Google Scholar]

36. Aagesen AM, Hase CC. Seasonal effects of heat shock on bacterial populations, including artificial Vibrio Parahaemolyticus exposure, in the Pacific oyster, Crassostrea Gigas. Food Microbiol. 2014;38:93–103. doi:10.1016/j.fm.2013.08.008. [PubMed] [CrossRef] [Google Scholar]

37. Varis J, Haverinen J, Vornanen M. Lowering temperature is the trigger for glycogen build-up and winter fasting in Crucian carp (Carassius Carassius) Zool Sci. 2016;33(1):83–91. doi:10.2108/zs150072. [PubMed] [CrossRef] [Google Scholar]

38. LØFstedt MB. The effect of food quality on glycogen content, the fatty acid profile and winter mortality in cultivated oyster spat (Ostrea Edulis) Aquac Nutr. 2010;16(6):625–636. doi:10.1111/j.1365-2095.2009.00700.x. [CrossRef] [Google Scholar]

39. Wang Y-R, Li E-C, Chen L-Q, Wang X-D, Zhang F-Y, Gao L-J, Long L-N. Effect of acute salinity stress on soluble protein, Hemocyanin, Haemolymph glucose and Hepatopancreas glycogen of Eriocheir Sinensis. Acta Hydrobiologica Sinica. 2013;36(6):1056–1062. doi:10.3724/SP.J.1035.2012.01056. [CrossRef] [Google Scholar]

40. Park JY, Kim CH, Hong SK, Suh KI, Lee KU. Effects of FFA on insulin-stimulated glucose fluxes and muscle glycogen synthase activity in rats. Am J Phys. 1998;275(2 Pt 1):E338–E344. [PubMed] [Google Scholar]

41. McGarry JD, Brown NF. The mitochondrial carnitine palmitoyltransferase system. From concept to molecular analysis. Eur J Biochem. 1997;244(1):1–14. doi:10.1111/j.1432-1033.1997.00001.x. [PubMed] [CrossRef] [Google Scholar]

42. Kerner J, Hoppel C. Fatty acid import into mitochondria. Bba-Mol Cell Biol L. 2000;1486(1):1–17. doi:10.1016/S1388-1981(00)00044-5. [PubMed] [CrossRef] [Google Scholar]

43. Li J, Stillman JS, Clore JN, Blackard WG. Skeletal muscle lipids and glycogen mask substrate competition (Randle cycle) Metabolism. 1993;42(4):451–456. doi:10.1016/0026-0495(93)90102-T. [PubMed] [CrossRef] [Google Scholar]

44. Argyraki M, Wright PD, Venables CW, Proud G, Taylor R. In vitro study of human skeletal muscle strips: effect of nonesterified fatty acid supply on glucose storage. Metabolism. 1989;38(12):1183–1187. doi:10.1016/0026-0495(89)90157-1. [PubMed] [CrossRef] [Google Scholar]

45. Witit*uwannakul D, Kim KH. Mechanism of palmityl coenzyme a inhibition of liver glycogen synthase. J Biol Chem. 1977;252(21):7812–7817. [PubMed] [Google Scholar]

46. Rognstad R. Rate-limiting steps in metabolic pathways. J Biol Chem. 1979;254(6):1875–1878. [PubMed] [Google Scholar]

47. Nomura Y, Okamoto S, Sakamoto M, Feng Z, Nakamura T. Effect of cobalt on the liver glycogen content in the streptozotocin-induced diabetic rats. Mol Cell Biochem. 2005;277(1–2):127–130. doi:10.1007/s11010-005-5777-y. [PubMed] [CrossRef] [Google Scholar]

48. Rosas C, Cuzon G, Gaxiola G, Le Priol Y, Pascual C, Rossignyol J, Contreras F, Sanchez A, Van Wormhoudt A. Metabolism and growth of juveniles of Litopenaeus Vannamei: effect of salinity and dietary carbohydrate levels. J Exp Mar Bio Ecol. 2001;259(1):1–22. doi:10.1016/S0022-0981(01)00222-2. [PubMed] [CrossRef] [Google Scholar]

49. Sanchez-Medina F, Sanchez-Urrutia L, Medina JM, Mayor F. Effect of muscular exercise and glycogen depletion on rat liver and kidney phosphoenolpyruvate carboxykinase. FEBS Lett. 1971;19(2):128–130. doi:10.1016/0014-5793(71)80495-7. [PubMed] [CrossRef] [Google Scholar]

50. Mendez-Lucas A, Duarte JA, Sunny NE, Satapati S, He T, Fu X, Bermudez J, Burgess SC, Perales JC. PEPCK-M expression in mouse liver potentiates, not replaces, PEPCK-C mediated gluconeogenesis. J Hepatol. 2013;59(1):105–113. doi:10.1016/j.jhep.2013.02.020. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

51. Goward CR, Nicholls DJ. Malate dehydrogenase: a model for structure, evolution, and catalysis. Protein Sci. 1994;3(10):1883–1888. doi:10.1002/pro.5560031027. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

52. Attwood PV. The structure and the mechanism of action of pyruvate carboxylase. Int J Biochem Cell Biol. 1995;27(3):231–249. doi:10.1016/1357-2725(94)00087-R. [PubMed] [CrossRef] [Google Scholar]

53. Wu H, Zhang X, Wang Q, Li L, Ji C, Liu X, Zhao J, Yin X. A metabolomic investigation on arsenic-induced toxicological effects in the clam Ruditapes philippinarum under different salinities. Ecotoxicol Environ Saf. 2013;90:1–6. doi:10.1016/j.ecoenv.2012.02.022. [PubMed] [CrossRef] [Google Scholar]

54. Coughlan BM, Moroney GA, van Pelt FNAM, O’Brien NM, Davenport J, O’Halloran J. The effects of salinity on the manila clam (Ruditapes philippinarum) using the neutral red retention assay with adapted physiological saline solutions. Mar Pollut Bull. 2009;58(11):1680–1684. doi:10.1016/j.marpolbul.2009.06.020. [PubMed] [CrossRef] [Google Scholar]

55. Cao C, Wang WX. Copper-induced metabolic variation of oysters overwhelmed by salinity effects. Chemosphere. 2017;174:331–341. doi:10.1016/j.chemosphere.2017.01.150. [PubMed] [CrossRef] [Google Scholar]