A Characterization of the RNA Modification Response to Starvation under Low Temperatures in Large Yellow Croaker (Larimichthys crocea)

Abstract

Emerging evidence shows that N6-methyladenosine (m6A) is a post-transcriptional RNA modification that plays a vital role in regulation of gene expression, fundamental biological processes, and physiological functions. To explore the effect of starvation on m6A methylation modification in the liver of Larimichthys crocea (L. crocea) under low temperatures, the livers of L. crocea from cold and cold + fasting groups were subjected to MeRIP-seq and RNA-seq using the NovaSeq 6000 platform. Compared to the cryogenic group, the expression of RNA methyltransferases mettl3 and mettl14 was upregulated, whereas that of demethylase fto and alkbh5 was downregulated in the starved cryogenic group. A Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis showed that the differentially m6A-modified genes were mainly enriched in steroid biosynthesis, DNA replication, ribosome biogenesis in eukaryotes, PPAR, ECM-receptor interaction, lysine degradation, phosphatidylinositol, and the MAPK signaling pathway, suggesting that L. crocea responds to starvation under low-temperature stress through m6A methylation modification-mediated cell growth, proliferation, innate immunity, and the maintenance of lipid homeostasis. This study advances understanding of the physiological response mechanism exerted by m6A methylation modification in starved L. crocea at low temperatures.

1. Introduction

Water temperature affects almost all the life activities of aquatic organisms, including growth, development, reproduction, metabolism, behavior, and geographical distribution [1]. The ability to cope with cold stress is critical for the survival of fish. Extremely low temperatures in winter have led to the massive mortality of cultured fish, seriously affecting the economic returns of farmers [2]. Studies have shown that Nile tilapia (Oreochromis niloticus) and Chanos chanos (Chanos chanos) are highly sensitive to cold stress and that low temperatures in winter frequently lead to massive mortality [3,4]. The ability of fish to tolerate low temperatures is related to their genetics and developmental stage, and juvenile fish are more sensitive to low-temperature stress than adult fish [5]. Therefore, exploring strategies and related mechanisms to better understanding of the cold tolerance of fish is necessary.

Large yellow croaker (Larimichthys crocea [L. crocea]) is a significant marine economic fish endemic to China; in marine fish farming, aquaculture production has ranked first for several years [6]. L. crocea is a warm-temperate fish, and its optimal growth temperature range is approximately 18–25 °C. It is less tolerant to low temperatures, will not feed below 8 °C, and cannot survive for a long time when it is lower than 5 °C [7]. Suitable starvation can increase the survival rate of zebrafish (Danio rerio) under cold stress and improve cold-stress resistance [8]. In our previous study, we reported that the cold + fasting group (CF) of L. crocea had significantly higher survival than the cold group (CC) [9]. Previous studies have shown that starvation stress helps L. crocea adapt to low temperature but is likely to affect their immune system [10]. However, the molecular mechanism by which starvation improves cold tolerance in L. crocea remains unclear.

RNA modification is critical for development and regeneration and has recently attracted attention [11]. With the advancement of RNA detection technology, there have been advances in the understanding of m6A methylation modification. Methyltransferases and demethylases including mettl3, mettl14, fto, and alkbh5 dynamically and reversibly regulate m6A methylation modification. mettl3 and mettl14 are methyltransferases, and fto and alkbh5 are two demethylases [12]. These RNA modification regulators play a significant role in several cellular processes, including mRNA splicing, translation, and degradation [13], and are involved in the regulation of various biological processes, including organismal development, cellular differentiation, and stress response [11,14,15]. Studies have shown that the expression of m6A methylation modification-related enzymes changes under heat stress in sheep, demonstrating that the m6A methylation modification of RNA is involved in heat-stress regulation [16]. Furthermore, it has been shown that m6A methylation modification can regulate the UV-induced DNA damage response in humans (Homo sapiens) [17]. Starvation stress can cause changes in epigenetic modifications that regulate differences in organismal phenotypes in zebrafish [18]. Studies related to m6A RNA methylation in starved L. crocea under low temperatures have not been reported. Therefore, this study aimed to investigate the expression pattern and biological role of m6A modification in starved and low-temperature exposed L. crocea using MeRIP-seq and RNA-seq techniques to understand the physiological response mechanism of L. crocea. The results of this study will provide a theoretical basis for enhancing the response mechanism of L. crocea to starvation stress at low temperatures.

2. Materials and Methods

2.1. Ethics Statement

All procedures for treating L. crocea were approved according to the relevant guidelines and regulations of the Committee for Laboratory Animal Research at the East China Sea Fisheries Research Institute of the Chinese Academy of Fisheries Sciences.

2.2. Animals and Sample Collection

The experimental juvenile L. crocea (mean body weight, 21.38 ± 2.46 g) were obtained from Fujian Fuding Research Center, East China Sea Fisheries Research Institute, Chinese Academy of Fisheries Sciences. The fish were temporarily reared in six tempered glass buckets with a 600 L volume for 2 weeks (30 fish/bucket). During the temporary rearing period, commercial feed was fed twice daily (08:00 and 16:00), and the residual feed was removed 1 h after feeding. The following were the mean water quality parameters: water temperature, 20 °C ± 0.42 °C; salinity, 24.78 ± 0.37; pH, 7.29 ± 0.34; dissolved oxygen, 7.31 ± 0.42 mg/L; and total ammonia nitrogen, 0.14 ± 0.02–0.19 ± 0.03 mg/L. At the end of the temporary rearing, the water temperature in the six buckets was reduced from 20 °C to 8 °C at a rate of 1.5 °C/day. To stop feeding, three buckets were randomly taken among the six buckets. At this time, the experimental fish were divided into a cold (CC) group and a cold + fasting (CF) group, with three replicates for each treatment group. Following 30-day cold exposure, the liver of one fish from each group was used for MeRIP-seq and RNA-seq.

2.3. qRT–PCR Analysis

methyltransferase-like 3 (mettl3), methyltransferase-like 14 (mettl14), fat mass and obesity-associated protein (fto), and alk b homolog 5 (alkbh5), which are genes related to m6A methylation modification regulation, were selected for real-time PCR. Total liver RNA was extracted using a TRIzol kit (Invitrogen, Carlsbad, CA, USA), and high-quality RNA was synthesized using a first-strand cDNA synthesis kit (Fermentas, Carlsbad, CA, USA). Real-time PCR primers were designed using Primer Premier 6.0 (Premier Biosoft International, Carlsbad, CA, USA) on the basis of the gene sequences in the L. crocea genome database (

Table 1. Primers used for real-time PCR analysis.

Primer Name	Forward Primer	Reverse Primer
β-actin	TCCTCGGTATGGAATCTTGC	GAGTATTTACGCTCAGGTGGG
mettl3	GCAGAGCAAGAAGGTGAGCAGAG	TCCCGTGGTCGCAGAACTCC
mettl14	TGCCCAACTCCCACCTATCCG	CCTCTGTCTCCGCCTCTTCCTG
fto	CCTCATGCTGCCACTGGAATCTG	TGAACTGAAGCGTGCGATGTCTC
alkbh5	CGGTGTTTGTCCTGCCTGTGAG	ACGCCGCTGCTTGATGTCTTG

). A 2 × QuantiFast SYBRGreen Master Mix kit (Roche, Leverkusen, Germany) was used to perform qRT–PCR on a fluorescence quantitative PCR instrument (Applied Biosystems Prism 7500 Sequence Detection System, Foster City, CA, USA). Three biological and three technical replicates were performed for each gene in each individual using β-actin as the internal reference gene. The 2^−ΔΔCt method was used to statistically analyze the expression levels of the target genes. One-way analysis of variance (ANOVA) (p < 0.05) was performed using Statistical Package for the Social Sciences (SPSS version 25, IBM, Armonk, NY, USA).

2.4. Total RNA Extraction, cDNA Library Construction, and Sequencing

Total RNA was extracted and assayed for quality using TRIzol reagent (Invitrogen) from liver samples of one fish from each of the CC and CF groups. Experiments were performed using 32 μg total RNA, and the fragmented mRNA was enriched with magnetic beads bearing anti-m6A, an antibody that captures mRNA fragments containing m6A methylation. Additionally, fragmented RNA was diluted 1:9 as the input control. Using the input data from MeRIP-seq, transcriptome analysis was performed. High-throughput sequencing was performed using the NovaSeq 6000 (Illumina, San Diego, CA, USA) platform.

2.5. MeRIP-seq Data Analysis

Splice removal and quality control were performed on raw data. FastQC [19] was used to analyze the quality of sequencing data to obtain information on the sequencing quality distribution, base content distribution, and proportion of duplicated sequencing fragments. HISAT2 software [20] was used to compare the filtered clean reads with the reference genomes to obtain unique-mapped reads for the next step of analysis (GCF_000972845.2). The enrichment of clean reads on the genome was determined using macs software [21], which is the classic ChIP-seq analysis software that can be used for Peak Calling for MeRIP-seq after some post-processing. Peak annotation analysis was performed after obtaining the peaks; subsequently, the peaks were analyzed for their location on the gene structure and overall distributional features. Motif analysis of the peaks was performed using DREME software [22] (https://meme-suite.org/meme/doc/dreme.html, accessed on 20 May 2023). The circlize package in the R software [23] (http://circos.ca/, accessed on 21 May 2023) was used for mapping the distribution of differential peaks across the genome. Moreover, the selected genes were visualized using integrated genomics viewer (IGV 2.17.0) software [24] (UC San Diego and the Broad Institute of MIT and Harvard, CA, USA) (http://software.broadinstitute.org/software/igv/, accessed on 22 May 2023).

Gene Ontology (GO) classification was applied to describe the function of genes, categorizing them into the following three main groups: biological process (BP), molecular function (MF), and cellular component (CC). The m6A differentially modified genes obtained from the analysis were based on the GO annotations from the database from the three levels of BP, MF, and CC, respectively. To calculate the significance level (p-value) of each GO, Fisher’s test was used and significant GO was associated with a p-value of <0.05. The top 20 terms were selected to draw the histogram [25]. The differential modified m6A genes were annotated on the basis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [26]; the significance level (p-value) of the pathway was calculated using Fisher’s test, and the significant pathway term for m6A gene enrichment was screened at a p-value of <0.05. To draw the enrichment circle diagram, the top 10 upregulated and downregulated pathways and their associated gene names were selected.

2.6. RNA-seq Data Analysis

Gene expression values were normalized by RNA sequencing using Reads Per Kilo Base Million Reads. The DEseq algorithm [27] (http://bioconductor.org/packages/release/bioc/html/DESeq.html, accessed on 6 June 2023) was used for the differential screening of fragments with log₂fold-Change >1 or <−1.

2.7. MeRIP-seq and RNA-seq Association Analysis

Plots were made using the annotation results of all peaks used for differential computation by MeRIP-seq and the differential analysis results of all genes by RNA-seq. log₂FC was used as a variable in RNA-seq for the horizontal coordinate, and log₁₀FC columns were used as a variable in MeRIP-seq for the vertical coordinate. The screening criteria considered only log₂FC and log₁₀FC values and did not consider the p-value from the two sequences. Where the criterion for RNA-seq was |log₂FC| > 1, the screening criterion for MeRIP-seq was |log₁₀FC| > 1. Four-quadrant plots show mRNA with differential m6A modification levels in the liver of L. crocea and their mRNA expression distribution. The following are the four types: the upper-left quadrant represents hypermethylated-down (hypermethylated-RNA level downregulation), the upper-right quadrant represents hypermethylated-up (hypermethylated-RNA level upregulation), the lower-left quadrant represents hypomethylated-down (demethylated-RNA level downregulation), and the lower-right quadrant represents hypomethylated-up (demethylated-RNA level upregulation). The number of transcripts in each quadrant was recorded.

2.8. Statistical Analysis

SPSS 25 software and Excel 2019 were used for statistical analysis, and one-way ANOVA was performed for comparison between groups. The results were expressed as means ± standard deviations, with p < 0.05 indicating a statistically significant difference.

3. Results

3.1. Starvation at Low Temperatures Enhanced the Expression of Genes Related to the Modification of m6A RNA Methylation

We first measured the expression of genes related to m6A modification by fluorescence quantitative PCR to explore whether m6A epigenetic modification is involved in the process of starvation mediating the enhancement of cold tolerance in L. crocea. The results (Figure 1, Table S1) showed that the expression of RNA methyltransferases mettl3 and mettl14 was upregulated in the CF group compared with that in the CC group (p < 0.05), whereas the expression of demethylases fto and alkbh5 was downregulated (p < 0.05). These results suggest that m6A RNA methylation modifications can be altered by starvation at low temperatures. Therefore, the MeRIP-seq technique was used for obtaining a specific modification landscape.

3.2. Statistics and Quality Assessment of MeRIP-seq

We determined the m6A data and the corresponding input data for two sets of liver samples using the MeRIP-seq technology. The results (

Table 2. Statistics for sequencing data.

Sample	Raw Reads	Q30	Clean Reads	Clean Ratio	Mapped Reads
CF.IP	42,980,372	90.97%	42,939,500	99.90%	27,123,510
CC.IP	47,704,922	90.96%	47,656,038	99.90%	32,767,664
CF.Input	77,079,832	93.42%	75,403,704	97.83%	57,529,263
CC.Input	90,898,260	93.74%	89,057,952	97.98%	65,467,101

Clean ratio refers to the ratio of clean reads to raw reads. Mapped reads refers to the number of reads matched to the genome.

) showed that for the immunoprecipitation library, approximately 4.30–4.77 × 10⁷ raw reads were obtained; for the input library, approximately 7.71–9.09 × 10⁷ raw reads were obtained. The Q30 values were all above 90.96%, indicating that the sequencing results are reliable and can be used for subsequent analysis. After filtering the raw reads, the clean reads obtained were more than 4.30 × 10⁷, which were all above 97.83% of the raw reads. The valid data were compared with the L. crocea genome sequence; it was noted that the number of mapped reads was between approximately 2.71 and 6.55 × 10⁷.

3.3. The Identification of the m6A Peak and Motif Analysis

As shown in Figure 2A,B (Table S2), 8419 and 9100 m6A peaks were screened in the CF and CC groups, corresponding to 4765 and 5068 genes, respectively, and 6441 (4016 genes) m6A peaks were the same in the two groups, accounting for more than 70% of all transcripts. The distribution of the m6A peak in different gene functional elements was analyzed, and the results are shown in Figure 2C,D. The m6A peaks in the CF and CC groups were mainly distributed in the protein-encoding (CDS) region, the start codon, and the stop codon region, with less distribution in the 5′-UTR and 3′-UTR regions. As can be observed in Figure 2E, the peak density in the CF and CC groups is significantly higher in the CDS region, near the start codon, and near the stop codon than in the 5′UTR and 3′UTR regions, verifying that the m6A peak is mainly distributed in these regions in both groups. An analysis of the motif features of m6A modifications in the CF and CC groups (Figure 2F,G) showed that the most significant m6A motif sequence in both groups comprised GGAGG.

3.4. Difference in m6A Peak Analysis

Overall, 919 transcripts with differentially methylated sites (peaks) were observed among the CF and CC groups, corresponding to 765 genes (Figure 3A, Table S3), of which 442 and 477 peaks were significantly upregulated and downregulated (Fold Change > 10 and p < 0.05), respectively, corresponding to 358 and 407 genes (p < 0.05). Differential peaks were distributed on every chromosome, particularly on chromosomes 3, 10, 12, and 13 (Figure 3B).

3.5. Differential m6A Methylation Gene GO Analysis

We performed functional enrichment analysis on all differential m6A modification genes from the BP, CC, and MF GO ontology functional groups to reveal further m6A modification epitranscriptomic information. The GO enrichment results (Table S4) showed that 358 significantly upregulated m6A modifier genes were significantly enriched in 300 GO entries (p < 0.05), including 128 BP, 97 CC, and 75 MF entries. Significantly upregulated m6A modifier genes were significantly enriched for DNA replication, rRNA metabolic process, and ncRNA processing in the BP category; the anchored component of the plasma membrane, the intrinsic component of the plasma membrane, and P-body in the CC category; and catalytic activity, acting on a nucleic acid, helicase activity, and ATP-dependent activity in the MF category (Figure 4A). A total of 407 significantly downregulated m6A-modified genes (Table S5) were significantly enriched in 248 GO entries (p < 0.05), including 130 BP entries, 74 CC entries, and 44 MF entries. Significantly downregulated m6A modifier genes were mainly enriched in the BP category, including cell adhesion, basement membrane organization, and blood coagulation regulation; extracellular space, T-tubules, and a Rad17 RFC-like complex in the CC category; and glycosaminoglycan binding, metalloendopeptidase activity, and inositol 1,4,5 trisphosphate binding in the MF category (Figure 4B).

3.6. An Analysis of the KEGG Functional Enrichment Pathways

Enrichment analysis was performed for all differential modified genes; the pathway information, corresponding to gene enrichment and ranked in the top 10 for transcripts with significantly upregulated or downregulated m6A modifications, (p < 0.05) was identified. The results (Table S6) showed that significantly upregulated m6A modification genes were mainly enriched in signaling pathways, such as steroid biosynthesis (ebp, krit1, hsd17b7, nsdhl, msmo1, and dhcr7), ribosome biogenesis in eukaryotes (xpo1a, utp4, tbl3, utp15, nol6, POP1, and nop58), DNA replication (pold2, pold1, dna2, rnaseh2c, and mcm2), and the peroxisome proliferator-activated receptor (PPAR) signaling pathway (acsl3b, SCD, acsl4a, plin2, SCD1, acsbg2, and pck1) (Figure 5A). Significantly downregulated m6A modification genes (Table S7) were mainly enriched in signaling pathways, including the extracellular matrix (ECM)–receptor interaction (itga4, thbs1b, RELN, fras1, and col1a1b), lysine degradation (kmt5ab, ehmt1a, mecom, and COLGALT1), the phosphatidylinositol signaling system (itpkcb, dgkaa, dgkzb, itpr2, and ip6k2a), and the MAPK signaling pathway (rasgrf2a, angpt1, MAPK8IP1, CHUK, mapta, mecom, tgfb3, cacna1g, and mapkapk2b) (Figure 5B).

3.7. An Integrated Analysis of MeRIP and RNA Sequencing

We performed a correlation analysis of MeRIP-seq and RNA-seq data (Table S8) to reveal the potential regulatory role of m6A modifications upon gene expression during starvation and low-temperature stress. As shown in Figure 6A (Table S9), RNA-seq data were filtered by |log₂FC| > 1, and MeRIP-seq is filtered by |log₁₀FC| > 1; 1527 eligible differentially methylated peaks (DMGs) were noted, of which 651 and 876 differentially methylated peaks showed upregulated and downregulated methylation modifications, respectively. Of the 651 methylation modification upregulated differential methylation peaks, 494 and 157 genes were synchronously upregulated (hyper-up) and reverse downregulated (hyper-down), respectively. Of the 876 methylation modification downregulated differential methylation peaks, 622 and 254 genes were synchronously downregulated (hypo-down) and reverse upregulated (hypo-up), respectively. These results mainly indicated that the level of m6A methylation modification was positively correlated with gene expression in starved L. crocea under low-temperature stress. As shown by the KEGG enrichment results of MeRIP-seq, we inferred that pck1 plays a role in the insulin, PPAR, and FoxO signaling pathways, whereas tnfsf10 plays a role in the FoxO signaling pathway. Therefore, we used an Integrative Genomics Viewer (IGV) to examine the m6A peak profiles of the pck1 and tnfsf10 genes. The results showed that, compared with the low-temperature group, the m6A level of the pck1 gene in the starvation and low-temperature group increased, whereas the m6A level of the tnfsf10 gene decreased (Figure 6B,C). These results indicate that RNA methylation is reversible and plays a dynamic role in the response of starved L. crocea under low temperatures.

4. Discussion

m6A is one of the most common epigenetic modifications in eukaryotic RNA and plays a fundamental regulatory role for gene expression [28]. We performed MeRIP-seq on samples from the liver tissues of starved L. crocea under low-temperature stress: 8419 and 9100 m6A peaks were screened in the CF and CC groups, corresponding to 4765 and 5068 genes, respectively. An analysis of the distribution of m6A peaks in different gene functional elements revealed that m6A peaks were mainly distributed in the protein-encoding (CDS) region, start codon, and stop codon region in the CF and CC groups, implying that RNA methylation modifications in these regions play a significant role in gene expression regulation in starved L. crocea under low-temperature stress. The m6A modifications were particularly enriched and highly conserved near the termination codon of CDS and the first one-fourth of the 3′UTR region in humans and mice [29]. The results of this study are consistent with those of these studies. m6A modification has been shown to delay the extension of mRNA translation, prevent mRNA degradation, and regulate translation with the help of the mRNA secondary structure in the CDS region, in addition to regulating the expression of genes in the 5′-UTR and 3′-UTR regions [30]. The presence of more liver methylation sites in the CDS region under starvation and low-temperature stress in L. crocea is presumed to be because the methylation of the CDS region can increase the translation opportunity, which is beneficial for starved L. crocea to cope with the environmental stress caused by low temperatures. A motif analysis showed that the most significant m6A motif in both groups was GGAGG, which was similar to the RGACH motif (R = G/A; H = A/C/U) observed in zebrafish (D. rerio) [31], suggesting that the m6A methylation sequence motif has certain similarities and conservation among vertebrates.

Compared with the CC group, 765 m6A differentially methylated genes were observed in the CF group, which contained 358 and 407 m6A differentially modified genes with upregulated and downregulated expression, respectively (p < 0.05). During the process of temperature decrease, metabolic rate, the secretion of digestive enzymes, and the binding capacity of iron ions are downregulated to different degrees, which may indirectly inhibit gene expression [32]. Hence, we suggest that when starved L. crocea was faced with low-temperature stress, the expression of some genes would be upregulated or downregulated to adapt to the low-temperature condition, and the fact that there were more downregulated than upregulated genes may be because low temperatures inhibit some physiological functions of L. crocea, which in turn leads to the downregulation of the expression of some genes. A GO functional enrichment analysis of m6A differentially modified genes showed that significantly upregulated m6A-modified genes were mainly enriched in DNA replication, rRNA metabolic processes, and ncRNA processing, which was consistent with the extensive and critical role of m6A methylation in the mRNA life cycle [33]. Therefore, we hypothesize that m6A methylation modification plays a significant role in regulating the post-transcriptional levels of genes involved in key metabolic and biological processes in response to low-temperature stress in starved L. crocea. Significantly downregulated m6A-modified genes were mainly enriched in cell adhesion, glycosaminoglycan binding, and metalloendopeptidase activity. One of the most essential biological features of cells and a key step in cell survival or death, membrane homeostasis maintenance, and signaling/metabolic pathway induction is cell adhesion [34]. An important first step in bacterial pathogenesis is epithelial cell adhesion involving glycosaminoglycans; the involvement of glycosaminoglycans can facilitate host–pathogen interactions during infection [35]. Therefore, we hypothesize that the m6A-modified cell adhesion pathway plays a vital role in starved L. crocea under low-temperature stress.

A KEGG pathway enrichment analysis revealed that significantly upregulated m6A modified genes were mainly enriched in steroid biosynthesis, DNA replication, ribosome biogenesis in eukaryotes, and the PPAR signaling pathway. Steroid is a general term for a large group of cyclopentanoperhydro-phenanthrene fully hydrogenated phenanthrene derivatives; steroids are essential for fish to maintain growth and survival under harsh environmental conditions. Different environmental and physiological conditions can affect the steroid biosynthetic capacity of cultured fishes [36]. The most significant signaling pathway for KEGG enrichment under oxidative damage conditions in grass carp (Ctenopharyngodon idella) and yellow catfish (Pelteobagrus fulvidraco) is the steroid biosynthesis signaling pathway [37,38]. Therefore, we hypothesized that starved L. crocea would respond to low-temperature stress through the m6A-modified steroid biosynthesis pathway. Ribosome biogenesis in eukaryotes is directly related to gene expression and plays a significant role in biological processes, including cell growth and proliferation [39]. The PPAR signaling pathway is a class of pathways related to lipid metabolism; PPARs are steroid hormone receptors that belong to the steroid/thyroid/retinoic acid receptor superfamily of transducing proteins, which are lipid-activated transcription factors that play an important pivotal role in maintaining lipid homeostasis [40]. Therefore, we hypothesize that in starved L. crocea under low-temperature stress, the m6A-modified PPAR signaling pathway plays a significant role in maintenance of lipid homeostasis. A significant downregulation of m6A modifier genes was mainly enriched in the ECM–receptor interaction, lysine degradation, phosphatidylinositol signaling pathway, and MAPK signaling pathway. The ECM–receptor interaction is an essential process for maintaining cell survival and performing its functions; the ECM can affect cell proliferation, migration, and differentiation by regulating the receptor signaling pathways. Studies have shown that ECM–receptor interactions play a significant role in the gut regeneration process of sea cucumber (Apostichopus japonicus), which accelerates cell migration movements [41]. Therefore, we suggest that some genes in the m6A-modified ECM–receptor interaction signaling pathway in starved L. crocea may be involved in the regeneration process of body-wall tissues under low-temperature stress. Lysine is an important metabolite in organisms; lysine residues produced by lysine degradation affect the signaling pathway of innate immunity, thereby inhibiting and negatively regulating it [42]. Therefore, the m6A-modified lysine degradation signaling pathway in starved L. crocea may affect the innate immunity of the organism under low-temperature stress. In the phosphatidylinositol signaling pathway, two second messengers are generated, which can open calcium channels when combined with ligands on the endoplasmic reticulum, thereby resulting in an increase in the intracellular metal ion concentration [43]. We hypothesize that the starved L. crocea maintains the intracellular metal ion concentration through the m6A-modified phosphatidylinositol signaling pathway, thereby leading to normal cellular excitation in response to cold stress. The MAPK signaling pathway is one of the highly conserved signaling pathways in eukaryotic cells, which transduces extracellular signals to the cytoplasm and nucleus, thereby controlling physiological processes in the organism. In vertebrates, the MAPK signaling pathway plays a key role in stress resistance, reproduction, cell development, differentiation, and inflammation [44]. This pathway is activated in orange-spotted grouper (Epinephelus coioides), yellow drum (Nibea albiflora), and olive flounder (Paralichthys olivaceus), among other species, in response to low-temperature stress [45,46,47]. Therefore, we hypothesize that the m6A-modified MAPK signaling pathway plays a significant role in starved L. crocea under low-temperature stress.

We further investigated the relationship between m6A differentially methylated peaks and the expression levels of differentially expressed genes to assess the role of modified m6A methylation in starved L. crocea under low-temperature stress. We observed that the gene expression levels of most of the peaks with altered methylation modifications were positively correlated with the presence of their m6A peaks. The expression levels of most differentially expressed genes in A. japonicus are positively correlated with m6A methylation modifications during disease development [28]. Moreover, in human genes, the gene expression levels of transcription start site fragments are positively correlated with m6A modifications of methylation [48]. However, some studies have shown that the expression levels of differentially expressed genes in porcine liver and chicken ovaries are negatively correlated with m6A methylation modifications [49,50]. On the basis of these studies, we hypothesize that the regulatory mechanisms of m6A methylation on target genes vary among species.

5. Conclusions

Our results demonstrate that m6A modification is a significant RNA epigenetic modification involved in normal physiology in the liver of starved L. crocea, which was observed to respond to low-temperature stress through m6A-modified steroid biosynthesis, DNA replication, ribosome biogenesis in eukaryotes, PPAR, the ECM–receptor interaction, lysine degradation, phosphatidylinositol, and MAPK signaling pathways. Gene expression levels were positively correlated with the modification of m6A methylation in the liver of L. crocea starved at low temperatures. Our results provide a basis for understanding RNA modification in L. crocea and extend our knowledge of the function of RNA modification.