The experimental animal protocols in the present study were reviewed and approved by the Animal Experimental Ethics Committee of Guangdong Ocean University, China (approval number: 0901-2019). Experiment procedures were performed in accordance with the Provisions and Regulations for the National Experimental Animal Management Regulations (China, July 2013) and the Experimental Animal Policies and Regulations of Guangdong Province (China, October 2010).

The experimental process (Figure 1A) included temporary fish culture, cold stress treatment, liver tissue collection under anesthesia, and total RNA extraction from samples and was performed following the methods described in our previous transcriptome study (Miao et al., 2021). Briefly, 150 healthy juveniles of the pearl gentian grouper hatched from Hainan (northern South China Sea, China) breeding parents were equally divided into three parallel experimental tanks. In the pre-experiment, we found that a considerable number of experimental fishes died at 11^°C. Accordingly, we set 12^°C as the lowest stress temperature point to obtain rich gene expression information and sufficient sample size under the extreme low-temperature stress. The water temperature in the three experimental tanks was gradually decreased from 25 to 20^°C, 15, and 12^°C at a constant rate of 1^°C/h and then kept at 12^°C for 6 h using a cooling-water machine. Three samples were randomly sampled from each tank at each of the five temperature points, for a total of 45 samples. These sampling groups were named CT (control group), LT20, LT15, LT12, and LT12-6. Considering that the liver is an important metabolic organ in the process of fish temperature adaptation (Hazel and Williams, 1990), it was sensitive to the transcriptomic changes of various pathways under temperature stress (Gracey et al., 2004; Sun et al., 2019). So, the liver of pearl gentian grouper was used for the targeted analysis sample in this study. After testing the total RNA quality of all 45 liver samples, 15 better quality RNA samples (taking one sample per parallel tank of all groups) were selected for subsequent small RNA library construction and sequencing.

Small RNA fragments were separately selected from each of the 15 total RNA samples using TruSeq Small RNA Library Preparation Kits (Illumina, San Diego, CA, United States) following the manufacturer’s instructions. Illumina adapters were directly added to the 5′ and 3′ ends of small RNA, and first strand cDNA was synthesized by reverse transcription and purified to obtain a cDNA library of about 150 bp sequences. After quality control using the Agilent 2100, small RNA libraries were sequenced on the Illumina HiSeq 2000 platform. Finally, 50 bp single-end reads were generated and used for subsequent analyses (Figure 1B).

Clean reads were obtained from raw data by a stringent filtering process, including the removal of low-quality reads (containing Ns), 5′ adapter contaminants, reads without 3′ adapter sequences, and reads containing poly-N sequences. Known miRNA, rRNA, tRNA, snRNA, snoRNA sequences were then identified and annotated by mapping small RNA clean reads to miRbase v22 (Kozomara et al., 2019) and Rfam v12.0 (Nawrocki et al., 2015) databases using bowtie2 v2.2.2 (Langmead and Salzberg, 2012). The programs mirEvo v1.3 (Wen et al., 2012) and miRdeep v2.0 (Friedländer et al., 2012) were integrated to predict novel miRNAs by identifying the Dicer cleavage site and secondary structure and calculating the minimum free energy of unannotated miRNA tags.

Expression levels of known and novel miRNAs in each sample were standardized by Transcripts Per Million (TPM). Pearson correlation coefficients (PCC) for expression levels in different samples were evaluated based on a standardized expression data of all miRNAs. Differential expression between groups was evaluated using the DESeq2 v3.18.1 (Love et al., 2014) R package. The DE miRNAs were determined using threshold values of Q < 0.05 and | log2 (Fold Change)| > 1. Subsequently, the DE miRNAs were visualized by volcano diagrams using the HiPlot.¹ Functional annotation and hyper-geometric distribution enrichment analyses of the target genes of DE miRNAs were performed using DIAMOND v0.9.24 (Buchfink et al., 2014) based on the Gene Ontology (GO) (Ashburner et al., 2006) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) databases. The GO terms and KEGG pathways were considered significantly enriched when FDR ≤ 0.05.

The STEM (Ernst and Bar-Joseph, 2006) program was used to characterize the expression profiles of miRNA that were consistently expressed with gradient temperature changes in all groups. First, the miRNA expression data standardized by the TPM algorithm were used as the input and then converted into preprocessed data by the log2 algorithm. Second, the maximum unit of change between time nodes and the minimum Pearson correlation coefficient were set to 2 and 0.7, respectively, for the construction of expression profiles. Finally, the miRNA expression profiles with statistical significance (P < 0.05) were screened for visualization and further analysis. The miRNA profiles that were significantly correlated with the temperature gradient were used for subsequent GO and KEGG enrichment analyses conducted using DIAMOND.

All full-length transcripts and DEGs were derived from our previous transcriptome data obtained by Illumina and PacBio sequencing. We applied MiRanda v3.3 (Riffo-Campos et al., 2016) and TargetScan v7.2 (Agarwal et al., 2015) to predict potential target genes of DE miRNAs based on DEG sequences, and overlapping target genes were selected for further analysis. The high-quality full-length transcript sequences were mapped to Animal TFDB v3.0 (Hu et al., 2019) to identify potential TFs using DIAMOND. A differential expression analysis of TFs among groups was performed using the DESeq2 v3.18.1 R package, with P < 0.05 and | log2 (Fold Change)| > 1 as thresholds to determine the DE TFs.

Pairwise relationships between DE miRNAs, DE TFs, and DEGs were evaluated using PCC > 0.7 and P < 0.05 as thresholds to construct a potential miRNA-TF-mRNA regulatory network. A network diagram was drawn using Cytoscape v3.7.2 (Shannon et al., 2003). Nodes with different shapes in the network represent different molecule types (i.e., miRNA, TF, and gene), different colors represent up- and down-regulation, and arcs and arrows indicate regulatory relationships. Further, by integrating pathways involved in each branch of the miRNA-TF-mRNA regulatory network, the cold-induced reaction process was revealed.

To verify the accuracy of sequencing data, five DE miRNAs and six target genes were randomly selected for qRT-PCR to detect expression trends in all groups. Specific primers were designed based on the miRNA and target gene sequences using Primer 5 (Table 1). Total RNA (1 μg) from each sample was used to synthesize cDNA using the TUREscript 1st Stand cDNA Synthesis Kit (AidLab, Beijing, China). Since 18S rRNA was verified by qPCR as stably expressed in all groups, it was used as an internal reference gene for miRNAs and target genes. The qRT-PCR system included 20 μL of reaction mixture consisting of 7.8 μL of RNase-free water, 10 μL of 2 × SYBR^® Green Supermix reagent (SYBR, Dalian, Liaoning, China), 1 μL of cDNA, and 0.6 μL of each primer. Amplification was performed using an Analytik Jena QTower 2.2 (Analytik Jena AG, Jena, Germany) under the following conditions: preheating at 95^°C for 5 min, followed by 40 cycles of 95^°C for 15 s and 60^°C for 30 s. All samples were evaluated in three repeated reactions and a no template control was included for each verified miRNA and target gene. Then, an amplification curve was generated to evaluate the range of cycle threshold (CT) values, and the melting curve was used to determine the specificity of primers and primer dimers. Lastly, the relative expression levels of miRNAs and target genes were calculated using the comparative CT method (2^−ΔΔCT).

A total of 4.82 Gbp of clean small RNA reads (Q30 > 97.24%) were obtained after filtering 11.28 Gbp of raw reads. As illustrated in Table 2, the GC content ranged from 49.16 to 49.80%, and the error rate per sample was < 0.01%. Among the 15 samples, the base numbers of clean reads and unique reads ranged from 284,153,887 to 284,153,887 and from 133,080 to 218,551, respectively. Sequence lengths were mainly 18–35 nt. The most frequent lengths were 22 nt (43.97%), 21 nt (21.48%), and 23 nt (15.12%).

The abundances of known miRNAs, novel miRNAs, rRNAs, tRNAs, snRNAs, snoRNAs, and other sequences in each sample are listed in Table 3. Among small RNA reads, the most frequent type was known miRNAs, followed by rRNAs, snoRNAs, snRNAs, and tRNAs.

In total, 668 known miRNAs were identified by mapping small RNA tags to the miRbase database. These miRNAs belonged to various miRNA families, such as mir0124, mir0190, mir0205, mir0218, and mir0455. In an analysis of the first base of mature miRNAs, adenine (A) or uracil (U) was the primary first bases. In addition, 48 novel miRNAs were predicted, of which 12 novel miRNAs were expressed in all samples, including novel_2, novel_7, novel_10, novel_33, novel_34, novel_38, novel_44, novel_47, novel_63, novel_69, novel_77, and novel_85. The secondary structures of these novel miRNAs were drawn to show the detailed stem-loop structure and mature miRNA sequences (indicated by red lines in Figure 2).

Expression data for all 716 miRNAs were standardized by the TPM algorithm and were used for analyses of correlations and differential expression among groups. PCC for pairwise comparisons between samples were in the range of 0.944–0.995 (Figure 3).

We detected 27 (5 up- and 22 down-regulated), 23 (5 up- and 18 down-regulated), 15 (1 up- and 14 down-regulated), and 38 (13 up- and 25 down-regulated) DE miRNAs in the CT vs. LT20, CT vs. LT15, CT vs. LT12, and CT vs. LT12-6 comparisons, respectively (Figure 4 and Supplementary Tables 1–4). Notably, 4 DE miRNAs were shared among the four comparisons. Moreover, 4 (1 up- and 3 down-regulated), 11 (4 up- and 7 down-regulated), 28 (12 up- and 16 down-regulated), 5 (2 up- and 3 down-regulated), 24 (12 up- and 12 down-regulated), and 17 (6 up- and 11 down-regulated) DE miRNAs were observed in the LT20 vs. LT15, LT20 vs. LT12, LT20 vs. LT12-6, LT15 vs. LT12, LT15 vs. LT12-6, and LT12 vs. LT12-6 comparisons, respectively. Obviously, the most DE miRNAs (38) were screened in the comparison between the CT group and the LT12-6 group. In addition, 89 unique miRNAs were screened in the pairwise comparisons among all groups, and 72 unique miRNAs were detected in the pairwise comparisons between the control group and the four low-temperature groups.

GO and KEGG pathway enrichment analyses of the DE miRNAs in the pairwise comparisons among groups were performed to further evaluate the functions of target genes regulated by miRNAs under cold stress. The main GO enrichment terms in all pairwise comparisons were generally similar (Figures 5A,C), including metabolic process, single-organism process, and cellular process in the biological process category, cell, cell part, and membrane in the cellular component category, and binding, catalytic activity, and transporter activity in the molecular function category. Several common pathways, such as vitamin digestion and absorption, biotin metabolism, AMPK signaling pathway, fat digestion and absorption, and phagosome pathways were within the top 20 significantly enriched pathways in the KEGG enrichment analysis (Figures 5B,D).

Using STEM trend analysis, 20 expression profiles were obtained for all miRNAs in the five groups, and expression changes in profile 0, profile 2, and profile 18 were statistically significant. As illustrated in Figure 6A, profile 0 was characterized by the overall down-regulation of miRNAs and profile 18 involved an overall up-regulation of miRNAs. The expression trend in profile 2 was unclear. In particular, 79 consistently down-regulated miRNAs were assigned to profile 0, most of which showed slight expression changes (Figure 6B). In profile 0, the normalized range of fold change values [log2(V(i)/V(0)] for most miRNAs with down-regulated expression was 0 to –4 [| Fold Change| ∈ (1, 16)], and a few miRNAs had values from –4 to –8 [| Fold Change| ∈ (16, 256)]. In contrast, 32 miRNAs were consistently up-regulated in profile 18 (Figure 6C), most of which were up-regulated in the LT20 group. Expression trends in profile 18 varied. Some miRNAs maintained their original expression levels after an initial down- or up-regulation, while others were continuously down- or up-regulated at all temperatures. The normalized range of values for up-regulated expression for the 32 miRNAs was 0–8 [| Fold Change| ∈(1, 256)].

Further, the GO and KEGG enrichment analyses of miRNA target genes assigned to profiles 0 and 18 were performed to reveal the key biological processes and pathways related to low-temperature responses. As demonstrated in Figure 7, the main GO terms associated with profiles 0 and 18 were similar, while the KEGG enrichment pathways involving miRNA target genes varied between the two profiles. In the GO enrichment analysis (Figures 7A,C), the main terms for the two profiles included metabolic process, cellular process, and single-organism process in the biological process category, cell, cell part, and membrane in the cellular component category, and binding, catalytic activity, and molecular function regulator in the molecular function category. Enriched KEGG pathways involving genes in profile 0 mainly included the AMPK signaling pathway, amino sugar and nucleotide sugar metabolism, phagosome, regulation of actin cytoskeleton, and arginine biosynthesis (Figure 7B). Profile 18 was mainly enriched in RNA transport, protein processing in endoplasmic reticulum, terpenoid backbone biosynthesis, platelet activation, and PI3K-Akt signaling pathways (Figure 7D).

Predicting target DEGs of all DE miRNAs provides a basis for annotation and enrichment analyses as well as the development of regulatory networks. Based on our 3,271 previously identified DEGs, 89 DE miRNA target genes were predicted in this study, and 9,637 miRNA–mRNA interacting pairs were obtained.

In the present study, 917 TFs were predicted against the Animal TFDB, and these could be assigned to 46 TF families. The top three TF families with the largest number of TFs were TF_bZIP (131 TFs), zf-C2H2 (124 TFs), and bHLH (103 TFs). In total, 63 (40 up- and 23 down-regulated), 138 (84 up- and 54 down-regulated), 148 (110 up- and 38 down-regulated), and 266 (143 up- and 123 down-regulated) TFs were differentially expressed in CT vs. LT20, CT vs. LT15, CT vs. LT12, and CT vs. LT12-6, respectively (Figure 8 and Supplementary Tables 5–8). Moreover, 25 (15 up- and 10 down-regulated), 91 (70 up- and 21 down-regulated), 186 (111 up- and 75 down-regulated), 56 (35 up- and 21 down-regulated), 129 (72 up- and 57 down-regulated), and 90 (35 up- and 55 down-regulated) DE TFs were identified in LT20 vs. LT15, LT20 vs. LT12, LT20 vs. LT12-6, LT15 vs. LT12, LT15 vs. LT12-6, and LT12 vs. LT12-6, respectively. The most DE TFs (270) were screened between the CT group and the LT12-6 group. Further analyses revealed that 368 unique DE TFs were obtained in pairwise comparisons among all groups, most of which (328) were obtained in pairwise comparisons between only the control group and the four low-temperature groups.

The factors related to low-temperature stress were screened by further exploring the interactions and regulatory relationships among DE miRNAs, DE TFs, and DEGs. As a result, 161 regulatory interactions involving the above-described 9637 miRNA-mRNA pairs were selected to construct a potential miRNA-TF-mRNA regulatory network (Figure 9), using thresholds of PCC > 0.7 and P < 0.05. Although there were fewer up-regulated miRNAs than down-regulated miRNAs, the up-regulated miRNAs regulated more target genes and TFs (Figure 9).

The direct regulation of target genes by miRNAs and indirect regulation of genes via TFs were both observed in the miRNA-TF-mRNA regulatory network (Figure 9). For example, gmo-miR-221-5p directly regulated CYP2J2 (Cytochrome P450 2J2), CYP4B1 (Cytochrome P450 4B1), TM7SF3 (Transmembrane 7 Superfamily Member 3), and GYS2 (Glycogen Synthase 2), ssa-miR-25-3-5p regulated TSKU (Tsukushi, Small Leucine Rich Proteoglycan) and NR1D2 (Nuclear Receptor Subfamily 1 Group D Member 2), and ccr-miR-489 and ccr-miR-122 regulated PPP1CC-A (Protein Phosphatase 1 Catalytic Subunit Gamma) and FOXK2 (Forkhead Box K2), respectively. The miRNA ssa-miR-7132b-5p indirectly regulated EHHADH (Enoyl-CoA Hydratase and 3-Hydroxyacyl CoA Dehydrogenase), ACAT2 (Acetyl-CoA Acetyltransferase 2), HADHA (Hydroxyacyl-CoA Dehydrogenase Trifunctional Multienzyme Complex Subunit Alpha), PCK1, PCK2 (Phosphoenolpyruvate Carboxykinase 1, 2) via the TF ACSS2 (Acyl-CoA Synthetase Short Chain Family Member 2), ola-let-7c indirectly regulated RXRAA, RXRBA, RXRGB (Retinoid X Receptor Alpha, Beta, Gamma) by the regulation of the TF PPARD (Peroxisome Proliferator Activated Receptor Delta), and gmo-miR-10545-5p had an indirect effect on PPP2R1A (Protein Phosphatase 2 Regulatory Subunit 1) and PTPA (Protein Phosphatase 2 Phosphatase Activator) by regulating the TF PPP4CB (Protein Phosphatase 4 Catalytic Subunit). In brief, these regulatory mechanisms and relationships may contribute to the response to cold stress in the pearl gentian grouper.

The relative expression levels of five miRNAs and six target genes in all groups were validated by qRT-PCR. As demonstrated in Figure 10A, all five miRNAs (i.e., ola-let-7c, ccr-miR-122, drc-let-7c-5p, gmo-miR-221-5p, and dre-miR-221-5p) were up-regulated. Among the six target genes, NPR2, FOXK2, NIPSNAP1, and MAF1 were down-regulated and PPP1CC-A and C7 were up-regulated. These results for miRNAs and target genes were consistent with those obtained by quantitative RNA-Seq (Figure 10) and regulatory network construction (Figure 9).

The oxidation of long-chain fatty acids is closely related to the stability of cell membrane components, and long-chain unsaturated fatty acids can maintain the fluidity of membrane fat and normal physiological functions (Jobin et al., 2013). Low temperatures could induce the abnormal oxidation of lipids and proteins, thus destroying the structure and function of the cell membrane (Malek et al., 2004). Hence, the antioxidant system plays a crucial role in the early stage of cold stress in animals.

In this study, we observed that the up-regulation of gmo-miR-221-5p inhibits the expression of CYP2J2 and CYP4B1; these two genes were involved in linoleic acid metabolism and biological oxidation, respectively. As cytochrome P450 monooxygenases, CYP2J2 and CYP4B1 are involved in the catabolism of polyunsaturated fatty acids (PUFAs) in various tissues (Wu et al., 1996), and the oxidation reaction involving CYP2J2 could reduce the fluidity of the cell membrane (Cheng et al., 2017). Accordingly, the down-regulation of CYP2J2 and CYP4B1 reduced the catabolism and oxidation of unsaturated fatty acids, which was conducive to the maintenance of the fluidity and normal function of cell membrane in the pearl gentian grouper under cold stress. On the other hand, the down-regulation of ssa-miR-7132b-5p weakened the inhibitory effect on the TF ACSS2, which reduced the expression of EHHADH. This gene was mainly involved in the fatty acid β-oxidation (unsaturated) pathway of pearl gentian grouper. As a bifunctional peroxisomal enzyme, EHHADH contains enoyl-CoA hydratase and 3-hydroxyacylcoa dehydrogenase, catalyzing the second and third steps of fatty acid β-oxidation of ultra-long chain fatty acid (Ranea-Robles et al., 2021). Therefore, the down-regulation of EHHADH reduced the oxidation of long-chain unsaturated fatty acids. This is beneficial for the resistance of oxidative damage in the cell membrane under cold stress and for the maintenance of the stability of unsaturated fatty acids. Taken together, these miRNAs and genes may reduce the oxidation of unsaturated fatty acids in the cell membrane.

In addition, we found that the down-regulation of ssa-miR-7132b-5p indirectly up-regulated the expression of ACAT2 and HADHA by regulating the TF ACSS2. The ACAT2 gene encodes cytosolic acetoacetyl-CoA thiolase, which participates in cholesterol biosynthesis (Kursula et al., 2005). Cholesterol is indispensable in animal cell membranes; it can effectively regulate the permeability, molecular order, and elasticity of lipid membranes (Xu et al., 2018). Thus, the up-regulation of the ACAT2 detected in our study could help to maintain the steady state of cell membrane cholesterol and ensure the basic fluidity and biological function of the cell membrane under cold stress. The HADHA gene acylates monomer dicarboxylic acid into cardiolipin, which is a major membrane phospholipid in eukaryotic mitochondria and plays a key role in cell apoptosis and ATP production (Taylor et al., 2012). This suggests that the up-regulation of HADHA was conducive to maintaining mitochondrial membrane function in pearl gentian grouper, thereby ensuring the progression of the electron transport chain in mitochondrial intima and the generation of ATP to provide energy during cold stress. We further discovered that the down-regulation of ssa-miR-25-3-5p weakened the inhibitory effect on TSKU. This gene was mainly involved in the ectoderm differentiation pathway of pearl gentian grouper. TSKU can stabilize the cholesterol content by reducing the circulation of high-density lipoprotein cholesterol, the outflow capacity of cholesterol, and the conversion of cholesterol into bile acids in the liver (Niimori et al., 2014). Therefore, we infer that the up-regulation of TSKU would maintain the stability of the cholesterol content in the cell membrane of pearl gentian grouper. This ensures the fluidity of the cell membrane under cold stress conditions. In brief, these miRNAs and genes may contribute to the fluidity of the cell membrane by maintaining essential lipid components, such as cholesterol and cardiolipin.

Collectively, we conclude that the three miRNAs gmo-miR-221-5p, ssa-miR-25-3-5, and ssa-mir-7132b-5p played pivotal roles in the processes of antioxidation and the maintenance of cell membrane fluidity under cold stress in pearl gentian grouper.

Energy is a necessary and important material source during adaptation to environmental temperature changes. Under continuous low-temperature conditions, a temperature decrease will directly cause the rapid loss of energy and further affect the normal biological and behavioral processes of animals (Sun et al., 2020). Carbohydrate and lipid metabolism are the main metabolic pathways involved in resistance to cold stress in fish (Wen et al., 2018).

In our regulatory network, the up-regulation of gmo-miR-221-5p inhibited the expression of both TM7SF3 and GYS2. These two genes were involved in the positive regulation of insulin secretion and glycogen biosynthesis, respectively. TM7SF3 promotes insulin secretion, thereby contributing to the inhibition of cytokine-induced pancreatic β cell death (Beck et al., 2011). Insulin reduces the blood glucose content to normal levels by reducing the production of glucose in the liver and promoting the storage of glycogen in both the liver and muscle and triglyceride fuel in adipose tissue (Najjar, 2003). By inhibiting TM7SF3, the gmo-miR-221-5p may weaken the hypoglycemic effect of insulin in response to cold stress, thereby stabilizing the blood glucose level in pearl gentian grouper. Glycogen synthase encoded by the GYS2 gene catalyzes glycogen synthesis. This synthase promotes the transfer of glucose molecules from UDP-glucose to the terminal branch of glycogen molecules (Orho et al., 1998). Consequently, the inhibition of GYS2 expression by gmo-miR-221-5p may prevent glycogen synthesis to maintain the blood glucose level and provide necessary energy for pearl gentian grouper under cold stress via glucose metabolism. We further discovered that the down-regulation of ssa-miR-7132b-5p indirectly up-regulated both PCK1 and PCK2 by regulating the TF ACSS2. The PCK1 and PCK2 genes participate in glycolysis and gluconeogenesis and in glucose/energy metabolism, respectively. At low glucose levels, PCK1 catalyzes the conversion of oxaloacetate to phosphoenolpyruvate (PEP) and generates glucose from other precursors derived from the lactic acid and citric acid cycle (Latorre-Muro et al., 2018). PCK2 catalyzes the compensatory transformation from phosphoenolpyruvate to oxaloacetate at high glucose levels (Latorre-Muro et al., 2018). Thus, the indirect up-regulation of both PCK1 and PCK2 could jointly regulate the stability of blood glucose levels, preventing energy loss caused by cold stress in pearl gentian grouper. In summary, the above miRNAs and genes may maintain the stability of blood glucose by participating in glucose metabolism under low-temperature conditions, e.g., by reducing insulin and glycogen synthesis.

We also found that the up-regulation of ola-let-7c inhibited the expression of the TF PPARD, thus decreasing the expression of RXRAA, RXRGB, and RXRBA genes, involved in fatty acid, triacylglycerol, and ketone-body metabolism pathway. The protein encoded by TF PPARD is considered an integratory of transcription inhibition and nuclear receptor signaling, which may inhibit the transcription activities of PPARA/PPARG (Peroxisome Proliferator-Activated Receptors Alpha and Gamma) induced by RXR genes (Van Der Veen et al., 2005). PPARA/PPARG are the key regulators of lipid homeostasis (Kliewer et al., 1997) and have been implicated in the regulation of lipid metabolism and adipocyte differentiation (Gorla-Bajszczak et al., 1999). Hence, the regulatory axis of ola-let-7c, TF PPARD, and three RXR genes promoted the maintenance of PPARA/PPARG activity in pearl gentian grouper, thereby maintaining lipid homeostasis in the liver and adipocyte differentiation in low-temperature environments. Furthermore, our results showed that the down-regulation of ssa-miR-25-3-5p weakened the inhibition of NR1D2 gene, involved in the nuclear receptor transcription pathway. NR1D2 encodes a member of the nuclear hormone receptor family and participates in a lipid metabolic pathway in a heme-dependent manner as a transcriptional inhibitor (Yu et al., 2018). NR1D2 regulates lipid and energy homeostasis in the skeletal muscle by inhibiting genes associated with lipid metabolism and regulates lipid metabolism in the liver by inhibiting APOC3 (Apolipoprotein C3) (Ramakrishnan et al., 2005; Bugge et al., 2012). Therefore, the up-regulation of NR1D2 detected in pearl gentian grouper was beneficial for the regulation of lipid metabolic homeostasis in the muscle and liver and for the maintenance of basic energy to resist cold stress. In short, the miRNAs and genes described are involved in lipid metabolism-related pathways and therefore may provide energy stored in lipids for necessary biological processes.

To sum up, the four miRNAs gmo-miR-221-5p, ssa-miR-7132b-5p, ola-let-7c, and ssa-miR-25-3-5p contributed to the regulation of blood glucose levels and energy in pearl gentian grouper under cold stress by directly or indirectly regulating the expression of their corresponding target genes and TFs.

In addition to its role in lipid metabolism, NR1D2 gene also regulates the circadian rhythm of pearl gentian grouper in a heme-dependent manner. The NR1D2 could directly inhibit the expression of the core clock components ARNTL (Aryl Hydrocarbon Receptor Nuclear Translocator Like) and CRY1 (Cryptochrome Circadian Regulator 1) and form a key negative regulatory branch of the circadian rhythm (Crumbley and Burris, 2011). NR1D2 plays a pivotal role in the regulation of the circadian rhythm in mammals (Laing et al., 2015) and has been a focus of human cancer research (Tong et al., 2020). In studies of fish, NR1D2 has only been shown to play a functional role in the circadian rhythm of muscle cells, which is beneficial for resistance to cell damage (Amaral and Johnston, 2012; Wu et al., 2016). However, this gene has not been reported in studies of fish responses to cold temperatures. In this study, we found that the down-regulation of ssa-miR-25-3-5p reduced the inhibition of NR1D2 in the pearl gentian grouper. Accordingly, we concluded that the up-regulation of NR1D2 may negatively regulate the core clock and its key elements in the circadian rhythm, thereby improving resistance to cold stress damage as well as cold tolerance and preventing cold shock in extremely cold environments. To our knowledge, our results provide the first identification that NR1D2 gene improves the cold tolerance of fish via the regulation of the circadian rhythm.

Moreover, we found that the down-regulation of ccr-miR-489 reduced the inhibition of PPP1CC-A gene. PPP1CC-A cooperates with CSNK1D (Casein Kinase 1 Delta) and CSNK1E (Casein Kinase 1 Epsilon) to jointly regulate the activity and stability of circadian clock elements and to determine the length of the circadian cycle (Song et al., 2015). PPP1CC-A also regulates the normal operation of the circadian clock by regulating PER1 (Period Circadian Regulator 1) and PER2 phosphorylation. Hung et al. (2016) confirmed that the core clock gene PER2 could enhance the cold tolerance of zebrafish larvae via its regulatory association with dre-miR-29b, although this gene was not identified in our study. Therefore, we inferred that PPP1CC-A was involved in the regulation of core clock genes in the circadian rhythm of pearl gentian grouper.

The miRNAs and genes described above may prevent pearl gentian grouper from entering cold shock in the later stage of lower-temperature stress by regulating the circadian rhythm, thereby improving the ability to cope with cold environments.

Cold stress causes a free radical imbalance, and free radical accumulation leads to DNA damage and cell function damage (Qiang et al., 2018b). When low-temperature stress exceeds the tolerance range of cells, cell signaling disruptions, extensive DNA damage, and apoptosis can occur (Elmore, 2007), affecting normal physiological activities and survival.

In this study, the down-regulation of gmo-miR-10545-5p reduced the inhibitory effect on the expression of TF PPP4CB and increased the expression of PPP2R1A gene. As an important protein phosphatase, the TF PPP4CB is involved in DNA repair, DNA damage checkpoint signaling, and apoptosis (Cohen et al., 2005). In particular, it catalyzes the dephosphorylation of RPA2 (Replication Protein A2) to increase expression and mediates the effective recruitment of RAD51 (RAD51 Recombinase) to chromatin, which is an important step in DNA repair (Lee et al., 2010). PPP2R1A encodes the scaffold A subunit of PPP2CA (protein phosphatase 2 catalytic subunit α, also known as PP2A) (Yang et al., 2013). As a substrate of RIDD (IRE1α-dependent decay), PPP2R1A participates in the RIDD catalytic reaction under DNA damage, which could promote DNA repair (Xu et al., 2018). Consequently, the up-regulation of PPP2R1A could provide sufficient substrate for RIDD, conducive to the repair of DNA damage. These results suggest that TF PPP4CB and PPP2R1A gene may protect against DNA instability and damage caused by cold stress in pearl gentian grouper.

Additionally, we found that the down-regulation of gmo-miR-10545-5p indirectly inhibited the expression of PTPA gene by regulating the TF PPP4CB. PTPA is mainly involved in cyclins and cell cycle regulation pathways in pearl gentian grouper. PTPA with peptidyl prolyl isomerase activity can stimulate weak phosphoserine phosphatase activity of the serine/threonine phosphatase PP2A (Chao et al., 2006). Luo et al. (2014) reported that the down-regulation of PTPA reduced the mitochondrial membrane potential, induces Box translocation into mitochondria, and releases Cyt C, which would induce apoptosis. Accordingly, the down-regulation of PTPA detected here may promote apoptosis via mitochondrial pathways, which is useful in maintaining tissue stability in pearl gentian grouper under cold stress. In our study, the up-regulation of ccr-miR-122 was observed to inhibit the expression of FOXK2 gene, involved in apoptosis and autophagy pathway. As a transcription inhibitor of autophagy in muscle cells and fibroblasts, FOXK2 inhibits the induction of autophagy via the TF FOXO3 (Bowman et al., 2014). This suggests that the down-regulation of FOXK2 could reduce the inhibitory effect on the TF FOXO3, which may enhance autophagy (i.e., the clearance of apoptotic cells and harmful substances) of cold stress. Taken together, the above miRNAs and genes may contribute to the stress response and reduce damage in pearl gentian grouper under extreme cold stress by regulating apoptosis and autophagy, thereby maintaining the basic functions of tissues and organs.

In conclusion, miRNA gmo-miR-10545-5p and ccr-miR-122 may play important roles in DNA repair, apoptosis, and autophagy by regulating their corresponding TF and target genes, thereby improving the survival of pearl gentian grouper in the later stage of low-temperature stress.