Comprehensive Analysis of Differentially Expressed mRNA, Non-coding RNA, and Their Competitive Endogenous RNA Network of Pacific Oyster Crassostrea gigas With Different Glycogen Content Between Different Environments

Glycogen content is a quantitative trait, its phenotype differences are found between individual oysters due to genetic effects and environmental factors which were including food, water temperature, salinity, and so on. In this study, a full sibling family of Pacific oyster Crassostrea gigas showed different phenotypes with high and low glycogen content between South Huanghai Sea (Rizhao offshore area, RZ) and North Huanghai Sea (Kongtong Dao area, KTD), respectively. At the same time, the content of 11 glucogenic amino acids and 13 fatty acids were also significant differences between RZ and KTD. RNA-seq and small RNA-seq technologies were used for transcriptome sequencing and functional enrichment analysis of differentially expressed RNA were used by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. A total of 2,084 mRNAs, 1,080 long non-coding RNAs (lncRNAs), 34 circular RNAs (circRNAs), and 7 microRNAs (miRNAs) were differentially expressed. Based on these differentially expressed genes (DEGs), miRNA target interactions (lncRNA/circRNA–miRNA pairs and miRNA–mRNA pairs) were predicted using the miRanda software. The differentially expressed mRNAs in this network were mainly shown to be involved in calcium signaling pathway and insulin signaling pathway. These findings could help to speculate that environmental factors may be epigenetically regulated by non-coding RNA in C. gigas, thereby further affecting glycogen content.


INTRODUCTION
Pacific oyster (Crassostrea gigas) has the characteristics of strong adaptability to the external environment, fast growth, high yield, rich nutrition, and delicious flavor. Therefore, it is widely cultivated worldwide as a high-economic value aquatic product. At present, the improvement of oyster meat quality is one of the important research topics for its industrial development, because the quality of meat is often measured by people's basic sensory evaluation of taste. Glycogen content accounts for 20-40% of the dry weight of oysters, which directly affects the taste of oysters and is an important quality trait. Previous studies have found that on the one hand, the content of glycogen metabolism-related enzymes directly affects the changes in glycogen content (Hata et al., 1993), on the other hand, the metabolic level of amino acids and fatty acids also affects glycogen metabolism (Randle et al., 1963;Fromentin et al., 2011). Therefore, genetic and environmental factors (temperature, nutrition, salinity, etc.) may affect the glycogen content (Li et al., 2017b). At present, several key genes that have been cloned in the oyster glycogen metabolism pathway, such as: glycogen phosphorylase (Hata et al., 1993), phosphoglucomutase (Tanguy et al., 2006), glycogen synthase kinase (Zeng et al., 2013), glycogenin (Li et al., 2017a), protein phosphatase 1 regulatory subunit 3B , etc. Glycogen content fluctuates seasonally with gonad development, and glycogenrelated genes express differences in different seasons and tissues (Bacca et al., 2005). And the increase of free fatty acid content could lead to reduce the consumption of glycogen. Similarly, high amino acid content also might increase glycogen content, because some amino acids could guide glycogen synthesis (Li et al., 2017b).
With the development of sequencing technology and the continuous improvement of the genome project, the human genome sequencing revealed that there are only about 20,000 protein-coding genes, and the proportion of the total genome sequence is <2%, and the remaining about 98% belong to noncoding RNA (ncRNA) (Chan and Tay, 2018). At present, ncRNA is divided into two types: housekeeping ncRNA and regulatory ncRNA. The common regulatory ncRNAs include microRNA (miRNA), long non-coding RNA (lncRNA), and circular RNA (circRNA). Studies have found that regulatory ncRNA is involved in the larval metamorphosis, shell pigmentation, and immunerelated mechanisms (Yu et al., 2016;Feng et al., 2018). However, whether ncRNA mediates the genetic mechanism that causes changes in glycogen content due to environmental factors is still unknown.
Therefore, in order to increase our understanding of the regulation of glycogen content, we performed ncRNA transcriptome analysis. During the research process, we placed a full sibling family of C. gigas in two offshore areas and produced two different phenotypes: high glycogen content and low glycogen content. The gonads of individuals with high and low glycogen content were selected for analysis. These findings may help us to further understand the regulatory effects of ncRNA on glycogen content in different environments, and provide new information for future genetic breeding of C. gigas in aquaculture.

Ethics Statement
In this study, all experiments were performed according to local and central government regulations. No endangered species were used in the study. No special permission was required to collect oysters or to perform the described experiments.

Experimental Animals and Sample Collection
In this study, we used a near-infrared (NIR) model  which constructed in our laboratory to determine the phenotypic value of the glycogen content of C. gigas. The families with high breeding value were selected through EBV calculation, and comprehensively consider the genetic relationship between families (the genetic relationship is less than 0.2) which could prevent the occurrence of inbreeding decline. During the oyster breeding season, the phenotypic value of glycogen content of the parents of each candidate family was determined by a NIR spectrometer (MicroNIR 1700 JDSU, United States) on site, and the male and female individuals with high glycogen content between each family were selected. Through mating between families, male and female parents from different families were fertilized to construct a full-sib family. Later, a full sibling family of C. gigas was farmed in the South Huanghai Sea (Rizhao offshore area, RZ) and North Huanghai Sea (Kongtong Dao area, KTD). In March 2018, we collected nine two-year-old oysters from RZ and KTD areas, and measured the glycogen content. In the RZ group and KTD group, three samples with similar growth characteristics were selected for follow-up experiments. We calculated the average glycogen content of RZ and KTD oysters, and each group selected three oyster individuals whose glycogen content was close to the average glycogen content of each group, and collected the gonadal tissues of these samples for RNA-seq and small RNA-seq, and determine its fatty acid content and amino acid composition. All samples were quick-frozen in liquid nitrogen and stored at −80 • C. by gas chromatography (GC). Simply, 2 ml methanol solution of hydrochloric acid, 1 mL n-hexanewere added into 100 mg powder sample and then hydrolyzed at 80 • C for 2 h. After cooling down, 3 mL 6% K 2 CO 3 was added to the samples. After filtered using a 0.22−µm membrane, the samples were analyzed on a high performance gas chromatography analyzer (GC-2010, Shimadzu, Japan), each injection volume was 1 µL. The types and relative contents of fatty acids were calculated according to the standard peak time screening and the corrected area normalization method.

Total RNA Isolation and Illumina Sequencing
RNA-seq and small RNA-seq were completed by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). The total RNA from gonadal tissue was isolated using TRIzol reagent (Invitrogen, CA, United States) according to the manufacturer's instructions. The RNA degradation and contamination was monitored on 1% agarose gels, purity was checked using the NanoPhotometer Spectrophotometer R (IMPLEN, CA, United States), concentration was measured using Qubit R RNA Assay Kit in Qubit R 2.0 Flurometer (Life Technologies, CA, United States), integrity was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States).
According to the manufacturer's instructions, a total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations. Firstly, ribosomal RNA was removed by Epicentre Ribo-zero TM rRNA Removal Kit (Epicentre, United States), and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, sequencing libraries were generated using the rRNA depleted RNA by NEBNext R Ultra TM Directional RNA Library Prep Kit for Illumina R (NEB, United States) following manufacturer's recommendations, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster GenerationSystem using TruSeq PE Cluster Kit v3-cBot-HS (Illumia). After cluster generation, the libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated.
A total amount of 3 µg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext R Multiplex Small RNA Library Prep Set for Illumina R (NEB, United States), library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. The clustering using TruSeq SR Cluster Kit v3-cBot-HS (Illumia) After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500/2000 platform and 50 bp single-end reads were generated.

Sequencing Results
In RNA-seq, raw data were obtained by removing reads containing adapter, reads on containing ploy-N and low quality reads. At the same time, Q20, Q30, and GC contents of the clean data were calculated. Raw sequencing reads were submitted to Sequence Read Archive in NCBI; the SRA accession numbers were SRR15293269, SRR15293270, SRR15293271, SRR15293272, SRR15293273, and SRR15293274 1 .
The reference genome and gene model annotation files are downloaded from the genome website. 2 Index of the reference genome was built using bowtie2 v2.2.8 and paired-end clean reads were aligned to the reference genome using HISAT2 (Langmead et al., 2009) v2.0.4. HISAT2 was run with "-rnastrandness RF, " other parameters were setas default. The mapped reads of each sample were assembled by StringTie (v1.3.1) (Pertea et al., 2016) in a reference-based approach.
In small RNA-seq, raw data were obtained by removing reads containing ploy-N, with 5 adapter contaminants, without 3 adapter or the insert tag, containing ploy A or T or G or C and low quality reads. At the same time, Q20, Q30, and GC contents of the clean data were calculated. Raw sequencing reads were submitted to Sequence Read Archive in NCBI; the SRA accession numbers were SRR15293267 and SRR15293268 (see text footnote 1).
The reference genome and gene model annotation files are downloaded from the genome website (see text footnote 2). The small RNA tags were mapped to reference sequence by Langmead et al. (2009) without mismatch to analyze their expression and distribution on the reference.

Identification of Long Non-coding RNA and Circular RNA
The spliced transcripts were screened for lncRNA according to the following steps: first, selected transcripts with exon number ≥2 and length >200 bp; second, Cuffcompare software was used to screen annotation lncRNA (Trapnell et al., 2010); third, selected the transcripts with FPKM ≥0.5; finally, CPC (Coding Potential Calculator) and PfamScan were used to screen for coding potential (Kong et al., 2007;Mistry et al., 2007). CircRNA was screened and identified in the constructed lncRNA library. Because of the high false positives in circRNA identification, circRNA was screened using find_circ and CIRI2 (Gao et al., 2018), the results were merged and the intersection of the two was taken.

Identification of MicroRNA
The mapped small RNA tags were compared with the Danio rerio sequence in miRBase 20.0 to find known miRNAs. The modified software mirdeep2 (Friedlander et al., 2012) and srna-tools-cli 3 were used to obtain potential miRNAs and draw secondary structures. To remove tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, small RNA tags were mapped to RepeatMasker, Rfam database. The software miREvo (Wen et al., 2012) and mirdeep2 were integrated to predict novel miRNA. To make every unique small RNA mapped to only one annotation, we follow the following priority rule: known

Target Gene Prediction
The lncRNA cis-acting target genes were screened for protein coding genes within 100 kb upstream and downstream; the lncRNA trans-acting target genes were screened for the Pearson correlation coefficient between lncRNA and mRNA than 0.95. Used miRanda (Enright et al., 2003) to predict the target gene of miRNA. The prediction was based on the miRNA binding site in the 3 UTR region of the gene. If the gene was uncertain or there was no 3 UTR region, then the miRNA binding site was determined according to the gene's stop codon within 1,000 bp downstream point.

Quantification of Gene Expression Level
Cuffdiff (v2.1.1) was used to calculate FPKMs of both lncRNAs and coding genes in each sample (Trapnell et al., 2010). Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group. MiRNA expression levels were estimated by TPM (transcript per million) through the following criteria (Zhou et al., 2010).

Differentially Expressed RNA Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis
During the analysis, Ballgown was used for differential expression of lncRNA and mRNA analysis (Pertea et al., 2016), DESeq R package (1.10.1) was used for differential expression circRNA analysis, and DESeq R package (1.8.3) was used for differential expression miRNA analysis. The p < 0.05 and | log 2 FoldChange| > 1 as the threshold for differential screening to identify differential genes.
The GOseq R software package (Young et al., 2010) and KOBAS software (Mao et al., 2005) were used to perform Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on differentially expressed genes (DEGs) or ncRNA target genes, respectively.

Construction of the Long Non-coding RNA/Circular RNA-MicroRNA-mRNA Network
Used miRanda to predict the interaction between DEMis and DEG, DELs and DECs (Friedlander et al., 2012) and Cytoscape V3.2 software to visualize it, and choosed to express the trend as "up-down-up" or "down-up-down" network for further research (Saito et al., 2012).

Quantitative Real-Time PCR Validation
To verify the results of RNA-seq and small RNA-seq, quantitative real-time PCR (qRT-PCR) was performed. The sample used for qRT-PCR was the same as the sequencing sample. The verification primers involved were designed using Primer Premier 5 software (Table 1), and were specifically compared by NCBI Primer Blast, synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). First, total RNA was used as the initial  (Renault et al., 2011). Finally, the Light Cycler 480 real-time PCR instrument (Roche Diagnostics, Burgess Hill, United Kingdom) was used for amplification, and 2 − ct method was used to calculate relative gene expression levels (Ballester et al., 2013).

Statistical Analyses
The Student's t-test analysis in SPSS statistical software v.22 (SPSS, Inc., Chicago, IL, United States) was used to analyze the data of glycogen content, fatty acid content and amino acid content. Quantitative data were expressed as means ± standard deviation (SD), and comparisons between different groups were completed using t-tests. The significance level for the statistical analyses was p < 0.05.

Phenotypic Statistics
The average glycogen content of the selected samples (fresh samples) in the RZ group and the KTD group were 7.173 ± 0.165 and 2.973 ± 0.055 mg/g, respectively. The differences in glycogen content between the two comparison groups were significantly different (p < 0.001). Therefore, in the follow-up analysis, we divided the two groups of oysters into high glycogen content of RZ (RZGH) group and low glycogen content of KTD (KTDGL) group. As shown in Table 2, a total of 16 amino acid compositions were measured. 11 glucogenic amino acids were identified that were more abundant in KTDGL oysters (p < 0.05), such as aspartic acid, threonine, serine, glycine, alanine, valine, methionine, isoleucine, tyrosine, phenylalanine, and arginine. As shown in Table 3, there are differences in the content of 13 fatty acids in the RZGH and KTDGL groups (p < 0.05), and the content of 9 fatty acids in the RZGH group was significantly higher than that in the KTDGL group.

Functional Analysis of Differentially Expressed mRNAs, Long Non-coding RNAs, and Circular RNAs Gene Ontology Enrichment Analysis on Differentially Expressed Genes
It mainly includes three categories: biological process (BP) category, cellular component (CC) category, and molecular Bold fonts were differentially expressed of glycogenic amino acids. function (MF) category. GO enrichment analysis found that 120 GO terms were significantly enriched (p < 0.05). Among the three categories, microtubule-based movement, proteinaceous extracellular matrix and intramolecular oxidoreductase activity, and interconverting aldoses and ketoses were the most significant ones, which were selected separately the top 10 terms of significance are displayed (Figure 2A). KEGG enrichment analysis was performed on these DEGs to further determine the metabolic process pathway. The results showed that DEGs were significantly enriched in 13 pathways, such as calcium signaling pathway, ECM-receptor interaction, and insulin signaling pathway, among which 29 genes were enriched into the calcium signaling pathway, and 17 genes were enriched into the insulin signaling pathway (Figure 2B). Including protein phosphatase 1 regulatory subunit 3B (LOC105340237, PPP1R3B), calmodulin (LOC105327344, CAM), protein tyrosine phosphatase (LOC105329458, PTP), insulin-like peptide receptor (LOC105348544, INSR), acetyl coenzyme A carboxylase (LOC105343342, ACCA), inositol Triphosphate 3-kinase B (LOC105336889, IP3KB), and so on.

The Cis and Trans Role of Long Non-coding RNAs in Target Genes
To investigate the possible functions of lncRNAs, the cis and trans roles of DELs in target genes were identified, and the filtered mRNAs were analyzed in combination with DEGs. A total of 471 lncRNA-mRNA interactions were found to have cis roles and 1,327 lncRNA-mRNA interactions were found to have trans roles. The target genes were analyzed by GO and KEGG pathway enrichment. In the results of GO enrichment analysis, it was found that fucose metabolic process, mRNA cap binding complex, and L-fucose isomerase activity were the most significant among the three categories, respectively. The top 10 terms of significance were selected for display ( Figure 3A). KEGG path analysis of lncRNA target genes to explore its function. Similar to the KEGG pathway of DEGs, it was found that 22 genes were enriched in the calcium signaling pathway, and 14 genes were enriched in the insulin signaling pathway, including PPP1R3B, CAM, PTP, INSR, DDO, etc. (Figure 3B).

Function Analysis of Differentially Expressed Circular RNAs
The functional enrichment of circRNA host genes was explored, 33 GO terms were significantly enriched (p < 0.05). KEGG pathway analysis was performed on the host genes of DECs to explore its function. The two pathways were significantly enriched (p < 0.05), other glycan degradation and propanoate metabolism. And it was found that novelcirc-0004211 were transcribed from acetyl-CoA synthase 2 (LOC105319828, ACAS2).

Construction of a Potential Long
Non-coding RNA/Circular RNA-MicroRNA-mRNA Regulatory Network A total of seven DEMis (four up-and three down regulated) were screened out with small RNA-seq. By analyzing the predicted target genes and DEGs, a total of 1,077 target genes were found. Analysis found that three down-regulated DEMis, novel-miR-34, novel-miR-78, and novel-miR-13 target a total of 283 up-regulated DELs; four up-regulated DEMis, novel-miR-154, novel-miR-44, novel-miR-225, and novel-miR-263 target a total   In addition, two up-regulated DEMis were found, novel-miR-154 targeted novel-circ-0001866; novel-miR-263 targeted novel-circ-0004857 and novel-circ-0003587. One down-regulated DEMi, novel-miR-78 targets novel-circ-0003765 and novel-circ-0000916. In the interacting network, the other novel-circ-0000916 and INSR gene were all negatively regulated by novel-miR-78. Some circRNA-miRNA-mRNA network interactions were shown in Figure 4B.
Gene Ontology and KEGG pathway analyses of the functions of DEGs were included in the network. A total of 106 GO terms were significantly enriched (p < 0.05), and the top 10 terms of significance in the three categories were selected to display, as shown in Figure 5A. The BP mainly involves the response of microtubule-based motion, fucose metabolism, etc. The MFs mainly involved in deoxyribonuclease activity, mercury ion transmembrane transport protein activity, and intramolecular involvement of aldose, etc. A total of 14 pathways were significantly enriched (p < 0.05), of which 25 genes were enriched into calcium signaling pathway, and 15 genes were enriched into insulin signaling pathway ( Figure 5B).

Quantitative Real-Time PCR Validation
The expression of twelve different genes was selected to confirm the utility of RNA-seq in qPCR analysis. The results showed that these transcripts were differentially expressed in oysters with different glycogen contents and were generally consistent with RNA sequencing data (Figure 6).

DISCUSSION
Oyster glycogen content is a quantitative trait, which varies greatly among individuals (Bacca et al., 2005;Aagesen and Haese, 2014). It is a complex and important trait that is affected by genetic and environmental factors, including seasons, nutrients, salinity, and temperature (Bacca et al., 2005;Varis et al., 2016). Liu's (2019) previous research showed that March is the proliferation period of oyster gonads, and the glycogen content of oysters is higher than the mature period of oyster gonads. Therefore, in this study, oysters of a full sibling family in two different sea areas RZ and KTD were farmed, to eliminate any potential genetic influence. In March 2018, we collected samples separately and found that the oyster glycogen in the RZ area was significantly higher than that of oysters in the KTD area. Since carbohydrates such as glycogen, amino acids, and fatty acids are the three major nutrients in organisms, their metabolic pathways are interrelated. An increase in free fatty acid levels may lead to an increase in glycogen content by reducing glycogen consumption (Randle et al., 1963). Similarly, a high-protein diet can increase glycogen content, because a large part of dietary amino acids are directed to glycogen synthesis (Fromentin et al., 2011). Therefore, we also measured the amino acid composition and fatty acid composition of the oysters in the RZGH and KTDGL groups, and the results found that there were significant differences in the amino acid composition and fatty acid composition of the oyster gonads between the two groups.
In this study, in order to get insight into the molecular regulation mechanism of glycogen content changes in oysters under different environments, we compared the differences in the expression profiles of mRNA, miRNA, lncRNA, and circRNA in the oyster gonads in the RZGH and KTDGL groups of a full sibling family. To determine the key factors involved in the regulation of glycogen content. Functional analysis of 2,084 differentially expressed genes showed that the insulin-like family genes were significantly enriched. Previous studies had shown that insulin could promote glycogen synthesis in mammals . Lymnaea stagnalis was the first shellfish to discover insulin-like family genes, and for invertebrates, the insulinlike family was widespread and highly conserved (Li et al., 1992). It was mainly involved in growth, development, and metabolism (Brogiolo et al., 2001). However, whether insulin-like in C. gigas was related to glycogen content had not been FIGURE 6 | RNA sequence data was verified using qPCR. (A) RNA-seq and q-PCR data for DEGs and DELs expression levels. (B) RNA-seq and q-PCR data of DECs and DEMis expression levels. Using eEF-1 and U6 as reference genes, the data represent the mean ± SD of three biological replicates, and each measurement was repeated three times. The fold change from the qPCR was calculated using the 2 − ct method.
reported. INSR is an insulin-like receptor protein tyrosine kinase substrate (Payankaulam et al., 2019), PTP dephosphorylates the tyrosine residues on the insulin receptor or its substrate, which prevents the insulin receptor from binding to insulin and negatively regulates insulin signal transduction (Haeusler et al., 2018). During the research, it was found that PTP was lowly expressed in RZGH oyster gonads, and INSR was highly expressed in RZGH oyster gonads. Therefore, it is speculated that the low expression of PTP in RZGH oyster gonads reduces the negative regulation of insulin-like signal transduction and further promotes glycogen synthesis. There are several genes related to glycogen synthesis in the insulin signaling pathway. For example, PPP1R3B was highly expressed in RZGH oyster gonad, it could cause glycogen synthase dephosphorylation, promote glycogen synthase activity, and further promote glycogen synthesis . LncRNA performs its function mainly by regulating mRNA, and analysis of the potential function of lncRNA is mainly predicted by analyzing their protein-coding target genes (Zhang et al., 2016). In the RZGH vs KTDGL group, it was found that lncRNA cis targets 471 mRNAs and trans targets 1,327 mRNAs. Functional enrichment analysis found that it was similar to differentially expressed mRNA and was related to insulin-like synthesis. For example, some of differentially expressed lncRNAs, XLOC-019254, XLOC-003125, XLOC-097516, and XLOC-020486 regulate CAM, PTP, INSR, IP3KB, etc., respectively. Therefore, lncRNA may play a certain regulatory role in insulin-like metabolism.
In eukaryotes, miRNAs have been reported to affect a variety of BPs, including glycogen metabolism (Dou et al., 2015), and their role in regulation is very important. With the development of high-throughput sequencing technology, studies on the miRNA transcriptome profiles of various mollusks have been reported (Xu et al., 2014;Zhou et al., 2014), and a new theory competitive endogenous RNA (ceRNA) has been proposed. However, in oysters, there is no report about miRNAs targeting glycogen metabolism in different environments. In this study, functional analysis of the seven differentially expressed miRNA target genes revealed that the differentially expressed miRNA might also be related to the insulin signaling pathway. The analysis found that in the network interaction analysis, INSR genes and lncRNAs XLOC-116536, XLOC-104062, etc., were all negatively regulated by novel-miR-13, novel-miR-34, and novelcirc-0000916 and INSR gene were all negatively regulated by novel-miR-78. These reports and our research results indicate that the ceRNA network formed by the interaction of lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA related to the insulin-like signaling pathway may play a role in the regulation of glycogen content in different environments, but the specific function needs further study.

CONCLUSION
In summary, we compared the expression characteristics of mRNA and ncRNA in the same high-glycogen content sibling family produced two different phenotypes of C. gigas with high and low glycogen content. The results show that there are a large number of ncRNAs in gonadal tissues, and bioinformatics analysis shows that ncRNAs are involved in insulin-like metabolism. These include novel-miR-78, novel-miR-154, novel-miR-225, and novel-circ-0000916 and XLOC-116536, XLOC-105344, and XLOC-104062, etc. These findings may help us to further understand the regulatory effects of ncRNA on glycogen content in different environments, and provide new information for selective breeding of C. gigas in future.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih. gov/Traces/study/?acc=SRP330386&o=acc_s%3Aa.

AUTHOR CONTRIBUTIONS
WW and JY contributed to the conception and design of the study. WW, XW, ZL, QL, and BL completed the material preparation, and data collection and analysis. GS, TX, XX, and YF supervised the data collection and analysis. WW and XW wrote the first draft. All authors read and approved the final manuscript.