The MicroRNAs-Transcription Factors-mRNA Regulatory Network Plays an Important Role in Resistance to Cold Stress in the Pearl Gentian Grouper

Pearl gentian grouper (Epinephelus fuscoguttatus ♀ × E. lanceolatus ♂) is a hybrid fish with high commercial value. It is widely cultured on the Asian coast; however, it is not cold-tolerant. Although we have previously characterized the liver transcriptomic responses of this grouper to cold stress, the roles of miRNAs and transcription factors (TFs) in cold resistance and the underlying regulatory mechanisms are still unclear. In this study, we integrated miRNA and mRNA sequencing data for pearl gentian grouper under cold stress and constructed a miRNA-TF-mRNA regulatory network. Furthermore, we screened seven key miRNAs (i.e., gmo-miR-221-5p, ssa-miR-7132b-5p, ola-let-7c, ssa-miR-25-3-5p, ccr-miR-489, gmo-miR-10545-5p, ccr-miR-122) that regulated target genes (including TF ACSS2, TF PPARD, TF PPP4CB; CYP2J2, EHHADH, RXRs, NR1D2, PPP1CC-A, PPP2R1A, FOXK2, etc.). These miRNAs participated in several important pathways and biological processes by the direct or indirect regulation of target genes, such as antioxidation and membrane fluidity, glucose and lipid metabolism, circadian rhythm, DNA repair, and apoptosis. The key cold-related miRNAs, TFs, and genes and their potential regulatory relationships identified in this study provide a deeper understanding of the complex molecular basis of the response to low-temperature environments in the grouper. In particular, our results provide the first identification for the role of NR1D2 gene in the cold tolerance of fish via the regulation of circadian rhythm. Furthermore, the key miRNAs and genes provide a basis for the molecular breeding of new cold-tolerant varieties of the pearl gentian grouper.


INTRODUCTION
As aquatic ectothermic animals, the survival and growth of fish are closely related to the water temperature (Yoneda and Wright, 2005). When the water temperature reaches or drops below the lower tolerance limit, a cascade of complex physiological, biochemical, and behavior responses occur in fishes (Donaldson et al., 2008). Many physiological, biochemical, and transcriptome studies have showed that in the process of continuous water temperature decline, the metabolic rate of fishes decreases gradually, accompanied by decreases in oxygen uptake and ammonia excretion (Li et al., 2017); however, the blood glucose content could increase via increased glucose and fat metabolism to supplement this energy loss (Cheng et al., 2017). When low-temperature stress exceeds the tolerance range in fish, hydrogen peroxide and free radicals produced by oxidative phosphorylation cannot be eliminated in a timely manner, which could cause irreversible damage to cells and lead to increased DNA repair and even apoptosis (Gabbianelli et al., 2003). These cold-induced responses involve a large number of genes and biological pathways. The underlying causes of transcriptomic changes are closely related to the regulatory roles of both microRNAs (miRNAs) and transcription factors (TFs) (Bartel, 2009). Only a few studies have identified the functional roles of miRNAs and TFs in the cold stress response in fish (Yang et al., 2011;Nie et al., 2019). Moreover, the main pathways involved in the response to cold stress differ substantially among studies.
MiRNAs are a class of small non-coding RNAs (18-25 nt) with pivotal roles in the post-transcriptional regulation of genes (Bushati and Cohen, 2007). The wide application of miRNA microarray and next-generation sequencing technologies has led to miRNA expression profiling studies of many species, especially mammals and plants (Lv et al., 2010). Despite several studies focused on the identification and functional roles of miRNAs in fish development, osmotic pressure regulation, and immune responses (Yan et al., 2012), relatively little is known about the functions of miRNAs in low-temperature stress in fish. Hung et al. (2016) confirmed that dre-miR-29b of zebrafish larvae can accurately regulate the core clock gene PER2 (Period Circadian Regulator 2) under cold shock stress, although weak correlations between cold-affected miRNAs and mRNAs were reported in an earlier study of cold-acclimated zebrafish (Yang et al., 2011). Subsequently, miRNAs associated with enhanced cold tolerance in the turbot Scophthalmus maximus have been shown to be related to circadian pathways and various signaling pathways (Nie et al., 2019). During cold stress in the common carp Cyprinus carpio, miRNAs regulate genes actively involved in energy metabolism and cell membrane components (Sun et al., 2020). Moreover, clear interactions among miRNAs and signaling pathways under cold stress have been revealed in the farmed tilapia Oreochromis niloticus (Qiang et al., 2018a). Collectively, numerous genes are directly and precisely regulated by miRNAs in several important biological processes involved in energy metabolism, antioxidation, the biological clock, and other processes, and miRNAs play key roles in genomic plasticity against low-temperature stress in fishes.
MiRNAs indirectly regulate a number of target genes via the regulation of TFs (Zhou et al., 2008). In particular, TFs can regulate the expression of target genes at a certain time point by combining with a specific sequence upstream of the 5 end of the gene (Lovering et al., 2021). Numerous studies have demonstrated the biological functions of TFs in the specific physiological responses of a wide variety of animals and plants (Liu et al., 2019). TFs also play essential roles in the regulation of fish embryonic development, cell differentiation, and hypoxia stress (Joseph et al., 2017). However, little is known about the functional roles of fish TFs and their effects under lowtemperature conditions. Nie et al. (2019) reported that TFs of the turbot S. maximus are mainly involved the cell membrane, signal transduction, and circadian rhythm pathways via regulatory genes under cold stress. TFs and miRNAs can accurately regulate gene expression pre-and post-transcription, respectively. Thus, the construction of a miRNA-TF-mRNA regulatory network will clarify the complex regulatory relationships involved in the response to low-temperature stress in fish.
Pearl gentian grouper (Epinephelus fuscoguttatus ♀ × E. lanceolatus ♂) is an F1 hybrid fish characterized by fast growth and strong disease resistance (Fan et al., 2014). It has become a popular grouper variety with high commercial value in the aquaculture industry along the Asian coast (Shapawi et al., 2014), particularly in southern China (Lin et al., 2017). However, this hybrid grouper is not cold-tolerant. Within the water temperature range of 20-13 • C, the feeding and activity of the pearl gentian grouper decrease significantly as the temperature decreases, and its lower critical temperature is 11 • (Liang et al., 2013). This obviously limits the breeding scope and environment of this grouper, which is not conducive to the expansion of its aquaculture industry. For this reason, we previously analyzed changes in the liver transcriptome of pearl gentian grouper under cold stress by integrating PacBio SMRT-Seq and Illumina RNA-Seq technologies and successfully screened various cold-related genes and biological pathways; our results indicated that energy homeostasis may play a crucial role in physiological responses to cold stress in the species (Miao et al., 2021). However, the regulatory relationships and dynamics of these genes are still unclear, and the construction of miRNA-TF-mRNA network is necessary to clarify regulatory mechanisms at the pre-and post-transcriptional levels underlying the transient response to low-temperature stress in this grouper and to narrow the screening of cold-related candidate genes. A clear understanding of these molecular mechanisms will contribute to the breeding of cold-tolerant fish varieties by the editing or knockout of cold-sensitive genes .
In this study, we investigated the complex regulatory relationships among miRNAs, TFs, and mRNAs under cold stress in the pearl gentian grouper and constructed a miRNA-TF-mRNA regulatory network. Based on high-quality sequencing data, we first identified miRNAs and screened differentially expressed (DE) miRNAs in cold-stressed fish, followed by analyses of the molecular functions and biological pathways associate with their target genes. Second, we used the STEM (Short Time-series Expression Miner) trend analysis to obtain expression profiles for miRNAs that were consistently expressed with decreasing temperatures. Finally, we integrated all DE miRNAs, DE TFs, and differentially expressed genes (DEGs) to construct a miRNA-TF-mRNA regulatory network and clarify important pathways and biological reactions. Our findings provide new insights into the molecular regulatory mechanism underlying adaptation to low-temperature environments in this grouper and are expected to facilitate the molecular breeding of cold-tolerant varieties.

Ethics Statement
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).

Experimental Process and Sample Collection
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.
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).

Identification and Differential Expression Analysis of MicroRNAs
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. 1 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.

Short Time-Series Expression Miner Trend Analysis and Enrichment Analysis
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.

Target Gene Prediction and Transcription Factors Identification
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.

Regulatory Network of MicroRNAs, Transcription Factors, and Target Genes
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 downregulation, 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.

qRT-PCR Validation of Key Differentially Expressed MicroRNAs and Target Genes
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 R 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 ).

Illumina Sequence Data for Small RNAs
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%).

Classification of Small RNAs and Identification of MicroRNAs
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).
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, singleorganism 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).

Expression Profiles of MicroRNAs With Temperatures
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  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 downregulated 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)].

illustrated in
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 lowtemperature 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).

Target Differentially Expressed Genes Prediction and Transcription Factor Identification
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  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.

Regulatory Network Analysis of MicroRNAs, Transcription Factors, and Target Genes
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).

Validation of the Differentially Expressed MicroRNAs and Target Genes
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).

DISCUSSION
Based on an integrated analysis of small RNA sequencing data and previous RNA-Seq data, we obtained 89 DE miRNAs, 368 DE TFs, and 3271 DEGs involved in the response to cold stress in the pearl gentian grouper. A miRNA-TF-mRNA network was constructed based on 23 DE miRNAs, 24 DE TFs, and 91 DEGs (Figure 9). By further exploring the functions and pathways related to these molecules, we finally confirmed that 28 key molecules, including 7 DE miRNAs, 3 DE TFs, and 18 DEGs, were involved in several important pathways and biological reactions associated with cold stress in pearl gentian grouper by both direct and indirect regulatory relationships (Figure 11). The main functions identified in this analysis could be summarized into four general categories (Figure 11): antioxidation and membrane fluidity, glucose and lipid metabolisms, circadian rhythm, and DNA repair and apoptosis.

Antioxidation and Membrane Fluidity to Resist Damage
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 downregulation 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 downregulation 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.

Glucose Metabolism and Lipid Metabolism
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 miRNA-TF-mRNA Regulatory Network 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 FIGURE 11 | Several important pathways and biological reactions associated with the miRNA-TF-mRNA regulatory network involved in the response to cold stress in pearl gentian grouper. 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.

Circadian Rhythm Prevents Cold Shock
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.

DNA Repair and Apoptosis
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 upregulation 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.

CONCLUSION
In this study, DE miRNAs, mRNAs, and TFs in pearl gentian grouper were screened, and a candidate miRNA-TF-mRNA regulatory network related to cold stress was constructed. The network involved two regulatory mechanisms, i.e., four miRNAs directly regulated eight target genes and three miRNAs indirectly regulated ten target genes via the regulation of three TFs. These miRNAs, TFs, and mRNAs played essential roles in several important pathways and biological reactions in pearl gentian grouper under cold stress, such as antioxidant and membrane fluidity, glucose and lipid metabolism, circadian rhythm regulation, DNA repair, and autophagy. In short, the two regulatory mechanisms revealed in pearl gentian grouper protected against damage caused by cold stress and promoted the maintenance of normal physiological activity. Of note, the potential function of NR1D2 gene in cold tolerance in fish via the regulation of the circadian rhythm was identified for the first time in this study. The complex regulatory relationships revealed in the present study not only provide a deeper understanding of the molecular mechanism underlying the adaptation of pearl gentian grouper to low-temperature environments but also provide clear directions for cultivating cold-tolerant varieties in the future.
We detected a number of TFs with important roles as bridges in our regulatory network. However, our prediction analyses revealed only a few TFs involved in miRNA regulation for further pathways analyses. Since the miRNA-mRNA regulatory relationship was the focus of our work, further studies of the detailed regulatory roles of additional TFs are needed, including ATAC-Seq studies of TFs in pearl gentian grouper and analyses of regulatory effects on target genes.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: PRJNA779537).

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Experimental Ethics Committee of Guangdong Ocean University.