Integrated SRNA-Seq and RNA-Seq Analysis Reveals the Regulatory Roles of miRNAs in the Low-Temperature Responses of Canarium album

: Chinese olive ( Canarium album ), a characteristic fruit tree in tropical and subtropical areas, suffers greatly from low-temperature stress (LTS). The regulatory roles of microRNA (miRNA) in plant LTS responses have been conﬁrmed in many plant species but not in C. album . In this study, a cold-tolerant cultivar ‘Rui’an 3 (cid:48) (RA) and a susceptible cultivar ‘Qinglan 1’ (QL) treated at 25 ◦ C (control, CK) and − 3 ◦ C (cold temperature treatment, CT) were subjected to small RNA (sRNA) and transcriptome sequencing for the exploration of the cold responses of C. album . Comparative sRNA sequencing analysis identiﬁed much fewer LTS-responsive, differentially expressed miRNAs (DEMs) in RA (4 DEMs) than in QL (23 DEMs). Cal-miR482-22 was found to be speciﬁcally induced by LTS in RA. Cal-miR397-3 was upregulated, while cal-miR398_2-3 and cal-undef-190 were downregulated after LTS only in QL. However, when compared with QL, a higher basic expression of cal-miR397-3, and lower expression of cal-miR398_2-3 and cal-undef-190 were found in RA, suggesting that they may contribute to the cold tolerance of RA. Comparative transcriptome analysis showed that the number of LTS-responsive differentially expressed genes (DEGs) identiﬁed in QL was larger than that in RA, and some DEGs were also predicted as the target genes of the identiﬁed DEMs, forming multiple differentially expressed miRNA–target gene pairs, such as cal-miR397-3_ laccase 2, 4, 17 , cal-miR482-22_ suppressor of npr1-1 , etc. Quantitative real time PCR results showed that the expression changes of DEGs and DEMs in different samples were generally consistent with the sequencing results. Our study indicated that the basic expression levels of some miRNAs (especially the cal-miR397-3, cal-miR398_2-3, and cal-miR482-22), and their target genes contribute greatly to the cold-tolerance characteristics of C. album . Our study is helpful for understanding the roles of miRNAs in the cold resistance and responses of C. album .


Introduction
Low-temperature stress (LTS) is one of the most common abiotic stresses threatening plant growth and development. It often causes extensive damages to crops, particularly in tropical and subtropical areas in late winter and early spring. Generally, chilling and freezing injuries will be caused in plants by above 0 • C and below 0 • C LTS, respectively. The plant's response to LTS is a complex physiological process involving in series of cellular, metabolic, transcriptional, and post-transcriptional changes, and changes in cell structure, antioxidant process, membrane structure, enzymatic activity systems and gene expression, etc., were widely studied in many plant species [1]. 2 of 20 MicroRNA (miRNA) is a kind of endogenous non-coding single-stranded small RNA. Generally, miRNAs display regulatory functions at the post-transcriptional level as negative regulators of gene expression [2]. Accumulated evidences have shown that miRNAs are involved in LTS response of many plants. In Ipomoea batatas, miRNAs function in the LTS response by regulating the expression of coding genes, including APETALA2 (AP2), auxin response factor (ARF), CNR transcription factor (CNR), dicer-like (DCL), MYB transcription factor (MYB), SQUAMOSA promoter binding protein like (SPL), teosinte branched1/Cincinnata/Proliferating cell factor (TCP), etc., and by mediating genes involved in salicylic acid, abscisic acid signaling, and reactive oxygen species (ROS) response pathways [3,4]. The LTS responsive miRNAs of Zea mays mainly include miR160, miR319, miR395, miR396, miR408, miR528, and miR1432 [5]. Some miRNAs, such as miR156k, miR159a, miR167a, miR169a, and miR172a, have been proven to contribute to the LTSinduced dormancy of Paeonia suffruticosa [6]. A large number of LTS-esponsive miRNAs have also been identified in Solanum aculeatissimum, including miR168a, miR2652a, miR812v, miR4414a-5p, miR5813, miR167c-3p, miR9478-3p, miR4221, and miR8577 [7]. Moreover, miR159, miR164, miR168, miR172, miR393, miR397, miR529, and miR1029 were reported to function in the cold stress responses of Triticum aestivum [8], and miR159, miR166, miR472, and miR482 function in the enhancing cold tolerance of S. lycopersicum [9] by interacting with target genes. Therefore, miRNAs play important roles in the LTS responses of plants.
Chinese olive (Canarium album) is a typical tropical and subtropical fruit tree primarily distributed in southern China and Southeast Asia. Its fruit is rich in nutrients and has important edible and medicinal values with great cultivation potential [10][11][12]. LTS is a serious natural threat for the C. album industry, often leading to the reduction of fruit yield and even plant mortality. It has been generally believed that −3 • C is the lowest temperature that C. album plants could withstand. In our previous study, we obtained a cold-tolerant C. album cultivar, named 'Rui'an 3' (RA), which could survive well when the temperature was below −3 • C. Additionally, we investigated the transcriptomic changes of the C. album cultivar 'Fulan 1', a special fresh edible cultivar, in response to LTS [13]. However, up to now, there is no report about the role of miRNAs in C. album cold response. Hence, an in-depth study on the functional mechanism of miRNAs and their target genes in response to LTS would be helpful for the systematical clarification of the cold resistance mechanism of C. album. In our present study, the expression changes of miRNAs and mRNAs of C. album in a cold-tolerant cultivar RA and a cold-susceptible cultivar 'Qinglan 1' (QL) under LTS were compared. The results obtained in this study will be helpful for revealing the function of miRNAs in the cold resistance and responses of C. album.

Materials
The cold-tolerant cultivar RA and susceptible cultivar QL plants used in this study were provided by the Fuzhou Germplasm Repository of Chinese Olive (Fuzhou, CHN). Seedlings with a uniform height of about 20 cm were cultured in an incubator at 25 • C, 60-80% relative humidity and 12 h (2000 ± 200 lx) photoperiod for 1 week (CK). Then half of the seedlings were moved to −3 • C (CT) artificial climate boxes for 24 h. The QL plants treated at 25 • C and −3 • C were named as QLCK and QLCT, and the RA plants treated at 25 • C and −3 • C were recorded as RACK and RACT, respectively. After phenotype observation, the leaves of plants from the four groups, i.e., QLCK, QLCT, RACK and RACT, were collected, frozen with liquid nitrogen and stored in a −80 • C ultra-low temperature freezer for future use. For each group, three biological replicates were made.

RNA Extraction and Quality Determination
The total RNA was extracted from C. album leaves by E.Z.N.A.TM Plant RNA Kit (Omega bio-tek, Norcross, GA, USA). The concentration and purity of the total RNA were measured using the Thermo Scientific NanoDrop 2000 (Thermo Scientific, Waltham, MA, USA) and the integrity was detected using the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA).

Library Construction and Sequencing
NEB Next Multiplex Small RNA Library Prep Set for Illumina and NEB Next ® Ultra TM RNA Library Prep Kit (NEB, Ipswich, MA, USA) was used for the construction of libraries for small RNA sequencing and mRNA sequencing. The library quality was measured using an Agilent 2100 Bioanalyzer and Agilent High Sensitivity DNA Kit (Agilent Tehcnologies Inc., Santa Clara, CA, USA). Pico green was used to detect the total library concentration (Quantifluor-ST fluorometer, Promega, Madison, WI, USA) and qPCR was used to quantitatively detect the effective library concentration. The multi-sample DNA libraries were homogenized, mixed in equal volume, diluted and quantified, and sequenced on the Illumina sequencer.

Small RNA Data Analysis
The quality information of raw data was calculated and then the raw data were filtered. Clean data were obtained by removing adapter and low-quality sequences. The number of clean reads with a sequence length of 18-36 nt were counted. The clean reads were combined with GenBank (https://www.ncbi.nlm.nih.gov/genbank/ (accessed on 18 October 2020)) and Rfam 14.2 (http://rfam.xfam.org/ (accessed on 18 October 2020)) databases to compare and annotate all small RNAs. The Arabidopsis database in miRBase 22.1 (http://www.miRbase.org/ (accessed on 19 October 2020)) was selected as the reference database, and the clean reads were aligned to the reference sequences. SRNAs with no match in the miRBase were further subjected to novel miRNA prediction using miREvo and miRdeep2. psRNAtarget (https://www.zhaolab.org/psRNATarget/ (accessed on 19 October 2020)) was used for target gene prediction with default parameters. For the identification of differentially expressed miRNAs (DEMs), |log 2 Fold change (FC) | > 1.0 and p < 0.05 were used as criteria.

Transcriptome Data Analysis
Trinity [14] was used for the transcriptome assembly. Gene function annotation was performed using NR, NT, Pfam, KOG, Swissprot, KEGG, and GO databases. DEGseq was used to analyze DEGs where p < 0.05 and |log 2 FC | > 1.0 were set as the threshold of differentially expressed genes (DEGs). Then, GO and KEGG analyses of DEGs was performed based on the GOseq R [15] package and KOBAS [16], respectively.

qRT-PCR Verification of DEMs and DEGs
β-Actin 7 and tubulin 5 were used as reference genes for the qRT-PCR verification of mRNAs and 18s rRNA was used as a reference gene for miRNAs. The TransScript ® Uni Allin-One First-Strand cDNA Synthesis SuperMix for qPCR (One-Step gDNA Removal) and TransScript ® miRNA First-Strand cDNA Synthesis SuperMix (Transgen Biotech, Beijing, China) were used for the synthesis of cDNA used for qRT-PCR verification of mRNAs and miRNAs, respectively. By using PerfectStart ® Green qPCR SuperMix (Transgen Biotech, Beijing, China), qRT-PCR reactions were performed on a Roche LightCycler480 fluorescence quantitative machine (Roche, Rotkreuz, Switzerland). The primer sequences are presented in Tables 1 and 2.

Phenotypic Changes of Cold-Tolerant and Susceptible C. album cultivars under LTS
After treatment at −3°C for 24 h, QL plants showed obvious damage symptoms, including leaf wilting, water loss, and browning of the leaf edges ( Figure 1A,B). The RA plants, however, did not show obvious cold stress-related symptoms ( Figure 1C,D).
The distribution of the first base of miRNAs with different lengths is presented in Figure 2A. The first base of miRNAs with lengths of 18 nt was mainly A and G, while that of miRNAs with lengths of 19 nt was mainly G and U. The first base of miRNAs with lengths ranging from 20 to 23 nt was mainly U. For the miRNAs with lengths from 24 to 26 nt, A was the major first base. And the first base of miRNAs with lengths of 28 nt was mainly C. The base distribution of miRNA sequences within 28 nt at positions 1-25 nt was
The distribution of the first base of miRNAs with different lengths is presented in Figure 2A. The first base of miRNAs with lengths of 18 nt was mainly A and G, while that of miRNAs with lengths of 19 nt was mainly G and U. The first base of miRNAs with lengths ranging from 20 to 23 nt was mainly U. For the miRNAs with lengths from 24 to 26 nt, A was the major first base. And the first base of miRNAs with lengths of 28 nt was mainly C. The base distribution of miRNA sequences within 28 nt at positions 1-25 nt was
Among these DEMs, the expression level of cal-miR164-1 increased after LTS in both QL and RA. The expression levels of some miR159, miR164, and miR396 family members were also upregulated. When compared with CK, DEMs belonging to the miR396 and miR397 families were both upregulated in QL and RA, while miR166 family members were downregulated. Notably, cal-miR397-3 was upregulated, while cal-miR398_2-3 and cal-undef-190 were downregulated in QL after LTS. The basic expression of cal-miR397-3 was higher and the expression of cal-miR398_2-3 and cal-undef-190 was lower in RA when compared with QL.

RNA-Seq and De Novo Assembly of Transcriptome
Comparative transcriptome analysis of QL and RA leaves before and after LTS treatment was performed. In total, an average of 47,306,978 raw reads were obtained from the 12 cDNA libraries. After removing the low-quality reads and reads with joints, 43,631,097 clean reads (accounting for about 92.22% of raw reads) were retained. The transcriptome data size of 12 samples ranged from 6.69 to 7.56 Gb. The Q20 and Q30 values for each library were higher than 98 and 94%, respectively, while the proportion of undefined bases was less than 8.0 × 10 −5 . The total length of the transcripts and unigenes was 337,927,555 and 87,678,110 bp, with an average length of 1380 and 998 bp, respectively. The N50 and N90 value of contigs was 2112 and 590 bp, respectively. The N50 and N90 values of unigenes were 1549 and 415 bp, respectively. In addition, the GC content of transcripts and unigenes was 39.40 and 38.50%, respectively.
The annotation proportions of all unigenes obtained in NR, GO, KEGG, Pfam, egg-NOG, and Swissprot databases were 53.43, 29.40, 21.95, 23.52, 51.19, and 37.54% respectively, and among all these annotated unigenes, 7945 (accounting for 9.04%) were annotated by all databases. Among these DEMs, the expression level of cal-miR164-1 increased after LTS in both QL and RA. The expression levels of some miR159, miR164, and miR396 family members were also upregulated. When compared with CK, DEMs belonging to the miR396 and miR397 families were both upregulated in QL and RA, while miR166 family members were downregulated. Notably, cal-miR397-3 was upregulated, while cal-miR398_2-3 and cal-undef-190 were downregulated in QL after LTS. The basic expression of cal-miR397-3 was higher and the expression of cal-miR398_2-3 and cal-undef-190 was lower in RA when compared with QL.

RNA-Seq and De Novo Assembly of Transcriptome
Comparative transcriptome analysis of QL and RA leaves before and after LTS treatment was performed. In total, an average of 47,306,978 raw reads were obtained from the 12 cDNA libraries. After removing the low-quality reads and reads with joints, 43,631,097 clean reads (accounting for about 92.22% of raw reads) were retained. The transcriptome data size of 12 samples ranged from 6.69 to 7.56 Gb. The Q20 and Q30 values for each library were higher than 98 and 94%, respectively, while the proportion of undefined bases was less than 8.0 × 10 −5 . The total length of the transcripts and unigenes was 337,927,555 and 87,678,110 bp, with an average length of 1380 and 998 bp, respectively. The N50 and N90 value of contigs was 2112 and 590 bp, respectively. The N50 and N90 values of unigenes were 1549 and 415 bp, respectively. In addition, the GC content of transcripts and unigenes was 39.40 and 38.50%, respectively.

DEGs Identification and Enrichment Analysis Results
By using the criteria of |log 2 FC | > 1.0 and p < 0.05, DEGs were identified from four comparisons, including QLCK vs. QLCT, QLCK vs. RACK, QLCT vs. RACT, and RACK vs. RACT. After LTS, more DEGs were identified in QL (11,026 in total with 5991 upregulated and 5035 downregulated) when compared to RA (8548 in total with 4238 upregulated and 4310 downregulated) ( Figure 4A, supplemental data Tables S3 and S4). Among these DEGs, 75 were identified as being common DEGs among all the comparisons ( Figure 4B). When compared with CK, the number of common DEGs both in QL and RA after LTS was 5180, and the number of specific DEGs in QL and RA were 5846 and 3368, respectively. These results indicated that the overall transcriptome change in RA after LTS was milder than that in QL.    KEGG enrichment analysis was performed for the DEGs identified in different comparisons. When compared with CK, 126 KEGG metabolic pathways were enriched after LTS in QL (supplemental data Table S5). Of them, 15 pathways, including phenylpropanoid biosynthesis, carotenoid biosynthesis, plant hormone signal transduction, starch and sucrose metabolism, sesquiterpenoid and triterpenoid biosynthesis, plant-pathogen interaction, plant circadian rhythm, plant MAPK signaling pathway, cutin, suberin, and wax biosynthesis, biotin metabolism, glyoxylate and dicarboxylate metabolism, anthocyanin biosynthesis, glycine, serine, and threonine metabolism, flavonoid biosynthesis, and isoflavonoid biosynthesis were extremely significantly enriched (p < 0.01) (Table 3). Additionally, seven pathways including galactose metabolism, glutathione metabolism, amino sugar and nucleotide sugar metabolism, glycerolipid metabolism, steroid biosynthesis, other types of O-glycan biosynthesis, and pentose and glucuronate interconversions were significantly enriched (p < 0.05) ( Table 3).   Table S6). Among these, 11 pathways, including plant hormone signal transduction, carotenoid biosynthesis, plant MAPK signaling pathway, plant circadian rhythm, starch and sucrose metabolism, glyoxylate and dicarboxylate metabolism, glycine, serine, and threonine metabolism, galactose metabolism, phenylpropanoid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, and carbon fixation in photosynthetic organisms, were extremely significantly enriched (p < 0.01) ( Table 4). Moreover, six pathways, including anthocyanin biosynthesis, tyrosine metabolism, ABC transporters, alanine, aspartate, and glutamate metabolism, steroid biosynthesis, and glycerolipid metabolism, were significantly enriched (p < 0.05) ( Table 4).  Phenylpropanoid biosynthesis, carotenoid biosynthesis, plant hormone signal transduction, starch and sucrose metabolism, sesquiterpenoid and triterpenoid biosynthesis, plant circadian rhythm, plant MAPK signal pathway, glyoxylate and dicarboxylate metabolism, anthocyanin biosynthesis, glycine, serine, and threonine metabolism, galactose metabolism, glyceride metabolism, steroid biosynthesis were significantly or extremely significantly enriched both in the comparisons of QLCK vs. QLCT and RACK vs. RACT, suggesting that these pathways play important roles in LTS responses in both the two C. album cultivars. Moreover, carbon fixation in photosynthetic organisms, tyrosine metabolism, ABC transporters, and alanine, aspartate, and glutamate metabolism pathways were found to be specifically enriched in the RACK vs. RACT comparison, which might contribute to the high cold tolerance of RA.
Additionally, the DEGs identified in the QLCK vs. RACK comparison and the QLCT vs. RACT comparison were found to be enriched in 98 and 99 KEGG pathways, respectively (supplemental data Tables S7 and S8), and five and four pathways were identified to be significantly enriched in the QLCK vs. RACK comparison and QLCT vs. RACT comparison, respectively (Tables 5 and 6) (p < 0.05). Notably, the phenylpropanoid biosynthesis pathway was significantly enriched in all comparisons (QLCK vs. QLCT, RACK vs. RACT, QLCK vs. RACK and QLCT vs. RACT), indicating that this pathway contributed greatly to the cold stress responses of C. album. In our previous study, plant hormone signal transduction was verified to be freezingstress-responsive in C. album cv. 'Fulan 1' [13]. In this study, plant hormone signal transduction was significantly enriched by DEGs identified from the QLCK vs. QLCT and RACK vs. RACT comparisons. In QL, AUX1, TIR1, and AUX/IAA genes, involved in auxin signaling pathways, were identified to be upregulated DEGs by LTS, while differentially expressed ARF, GH3, and SAUR genes contained both up-and downregulated members. However, in RA, GH3 was upregulated, AUX/IAA, ARF, and SAUR DEGs both had up and downregulated members, while AUX1 and TIR1 did not display differential expression change after LTS ( Figure 5A). In QL, the NPR1 and PR-1 genes involved in the salicylic acid signaling pathway were downregulated by LTS, but TGA genes had both up-and downregulated members. However, in RA, PR-1 was upregulated by LTS, and both upand downregulated NPR1 and TGA genes members were identified ( Figure 5B).

KEGG Enrichment Analysis of DEM Target Genes
Target gene prediction analysis using psRNAtarget identified 2743 target genes of 26 DEMs, including 48 and 4 target genes with expectation values less than 2.0 and equal to 0, respectively. Pathway enrichment analysis of these predicted DEM target genes showed that, in QL, the target genes of LTS-responsive DEMs were enriched in 26 pathways. The top ten enriched pathways included pyrimidine metabolism, amino sugar and nucleotide sugar metabolism, biotin metabolism, cutin, suberin, and wax biosynthesis, carbon fixation in photosynthetic organisms, pyruvate metabolism, carotenoid biosynthesis, plant-pathogen interaction, phenylalanine metabolism, and ubiquinone and other terpenoid-quinone biosynthesis (Table 7). Additionally, peroxisome, phenylpropane biosyn-thesis, RNA degradation, and others related to plant stress resistance pathways were also enriched. However, only pyrimidine metabolism was significantly enriched (p < 0.05).

KEGG Enrichment Analysis of DEM Target Genes
Target gene prediction analysis using psRNAtarget identified 2743 target genes of 26 DEMs, including 48 and 4 target genes with expectation values less than 2.0 and equal to 0, respectively. Pathway enrichment analysis of these predicted DEM target genes showed that, in QL, the target genes of LTS-responsive DEMs were enriched in 26 pathways. The top ten enriched pathways included pyrimidine metabolism, amino sugar and nucleotide sugar metabolism, biotin metabolism, cutin, suberin, and wax biosynthesis, carbon fixation in photosynthetic organisms, pyruvate metabolism, carotenoid biosynthesis, plantpathogen interaction, phenylalanine metabolism, and ubiquinone and other terpenoidquinone biosynthesis (Table 7). Additionally, peroxisome, phenylpropane biosynthesis, RNA degradation, and others related to plant stress resistance pathways were also enriched. However, only pyrimidine metabolism was significantly enriched (p < 0.05).
When compared with CK, the predicted DEM target genes in RA were enriched in only five pathways (Table 8), including propanoate metabolism, glyoxylate and dicarboxylate metabolism, pyruvate metabolism, amino sugar and nucleotide sugar metabolism, and glycolysis/gluconeogenesis. Among them, propanoate metabolism was significantly enriched (p < 0.05). Moreover, the target genes of LTS-responsive DEMs involved in propanoate metabolism, glyoxylate and dicarboxylate metabolism, and glycolysis/gluconeogenesis were identified to be enriched in both QL and RA.  When compared with CK, the predicted DEM target genes in RA were enriched in only five pathways (Table 8), including propanoate metabolism, glyoxylate and dicarboxylate metabolism, pyruvate metabolism, amino sugar and nucleotide sugar metabolism, and glycolysis/gluconeogenesis. Among them, propanoate metabolism was significantly enriched (p < 0.05). Moreover, the target genes of LTS-responsive DEMs involved in propanoate metabolism, glyoxylate and dicarboxylate metabolism, and glycolysis/gluconeogenesis were identified to be enriched in both QL and RA.

Expression Analysis of DEMs and Their Corresponding Target Genes
The differentially expressed target genes of DEMs included scarecrow-like protein (SCL), laccase (LAC), auxin response factor (ARF), U-box domain-containing protein (U-box protein), transketolase (TKT), heavy metal-associated isoprenylated plant protein (HIPP), resistance to uncinula necator protein (RUN), 4-coumarate-CoA ligase (4CL), suppressor of npr1-1, constitutive 1 (SNC1), E3 ubiquitin-protein ligase, and other genes with unknown functions (Table 9). According to the FPKM values, the expression of several miRNAs and their target genes varied greatly in different samples. These included cal-miR171_2-1 and SCL6, cal-miR171_2-6 and SCL6, cal-miR394-1 and U-box protein, cal-miR396-17 and TKT, cal-miR397-3 and LAC2, LAC4, LAC17, cal-miR482-22 and SNC1, cal-undef-161 and 4CL, and cal-undef-8 and E3 ubiquitin protein ligase. Noteworthily, after LTS, the expression of cal-miR482-22 (log 2 FC = 1.668) in RA was significantly negatively correlated with its target gene SNC1 (log 2 FC = −3.987). To reveal putative interactions between miRNAs and mRNA regulating LTS responses in C. album, the Pearson correlation coefficient (PCC) between DEMs and target genes was calculated based on the FPKM values and the condition of miRNAs or mRNAs where PCC ≥ 0.8 was selected to construct the regulatory network diagram ( Figure 6). The DEMs and the differentially expressed target genes after LTS were verified by qRT-PCR (Figure 7). The expression patterns of miRNAs and their corresponding target genes in different samples were almost all consistent with the sequencing results, suggesting that these miRNAs jointly regulate the LTS response in C. album via target gene interaction. Additionally, the expression of cal-miR488-22 and its predicted SNC1 was confirmed to be negatively correlated in RA after LTS by using qRT-PCR. Table 9. Differentially expressed target genes of DEMs after low temperature stress treatments. *-a significant difference of gene or miRNA in these samples. Inf-the FPKM value was zero in CK.   Table 9. Differentially expressed target genes of DEMs after low temperature stress treatments. *a significant difference of gene or miRNA in these samples. Inf-the FPKM value was zero in CK.

Discussion
Low-temperature stress is an important environmental factor influencing plant production. The plant responses to LTS are very complex and diverse [17,18]. In our present study, we compared the expression changes of sRNAomes and transcriptomes, screened the regulatory miRNAs, and obtained several potential miRNA-mRNA pairs in cold-tolerant and -susceptible C. album cultivars under cold stress. The preliminary analysis of their interactions and potential roles in cold resistance response will provide new insights for revealing the underlying response mechanisms of C. album to LTS.

Some miRNAs Contribute to the Cold-Tolerance of C. album
MiRNAs are an important class of non-coding RNA that are commonly involved in plant responses to abiotic stress via interacting with their corresponding target mRNAs or other non-coding RNA molecules (such as long non-coding RNA) [19,20]. As a characteristic crop of tropical and subtropical areas, C. album is sensitive to low temperature. In this research, we identified 547 known miRNAs and 237 novel miRNAs from C. album leaves. Among these miRNAs, 21 known and five novel miRNAs showed differential expression in response to LTS. These differentially expressed known miRNAs belong to 11 families, most of these miRNA families, including miR159 [8], miR164 [8], miR166 [21], miR171 [22], miR319 [21][22][23], miR394 [24], miR396 [21,25,26], miR397 [8], miR398 [21], and miR482 [9], which were all reported to play roles in regulating cold or freezing stress responses in plants such as Arabidopsis thaliana, Populus simonii×, Medicago sativa, Z. mays, and so on. Consistently, some miR159, miR164, and miR396 family members were found to be upregulated after LTS in both QL and RA. Previous studies have shown that miR159 acts as a small RNA-regulating calcium-dependent protein kinase (CDPK) and is up-regulated during cold stress in S. lycopersicum and S. habrochaites [27]. MiR164 targeting NAC transcription factors and phytosulfokines play important functions in signal transduction and oxidative stress responses, and are upregulated by cold stress in T. aestivum [8]. The miR396 that targets growth-regulating factors (GRFs) is upregulated in cold stress responses in Hylocereus polyrhizus [25].
The overexpression of miR397a may modulate the lignification of plant cell walls to improve the endurance of low-temperature stresses in plants [28], whereas miR397b is slightly upregulated by cold stress treatments [29]. In response to cold stress, the expression level of miR398 decreased in T. aestivum [30]. In our study, we found that cal-miR397-3 was upregulated, whereas cal-miR398_2-3 was downregulated in QL, but both of them were not significant changed in RA after LTS. Notably, the basic expression of cal-miR397-3 was higher and the expression of cal-miR398_2-3 was lower in RA, suggesting that their basic expression levels might related with the cold tolerance of C. album cultivars. Besides, miR482 has also been reported to be closely related to plant stress resistance [9]. In our present study, cal-miR482 was identified to be specifically upregulated only in RA after LTS, suggesting that this miRNA might contribute to cold-tolerance in C. album. These results indicated that the basic expression levels of some miRNAs, including cal-miR397-3, cal-miR398_2-3, and cal-miR482-22, contribute greatly to the cold-tolerance characteristics of C. album. Additionally, we also identified several novel miRNAs that were differentially expressed in different samples, including cal-undef- 8, 17, 161, 190, and 204. Further researches are needed to identify the roles of these miRNAs in the low temperature responses of C. album.

The Low Temperature Responses of C. ablum Involved Multiple Metabolic Pathways
A large number of DEGs involved in LTS responses were identified to be enriched in phenylpropanoid biosynthesis, carotenoid biosynthesis, plant hormone signal transduction, starch and sucrose metabolism, sesquiterpenoid and triterpenoid biosynthesis, and plant circadian rhythm in both QL and RA, which were also found to be enriched in cold-treated C. album cv. 'Fulan 1' [13]. It was thus suggested that the genes involved in these pathways may be important for LTS responses in C. album. However, sesquiterpenoid and triterpenoid biosynthesis, plant MAPK signal pathway, and glyoxylate and dicarboxylate metabolism were specifically enriched in our study. Sesquiterpenoid and triterpenoid biosynthesis could reduce the lipid oxidation under stress conditions in S. lycopersicum [31,32]. The plant MAPK signal pathway has been reported to be closely related to low-temperature stress responses in plants [33]. Brassica rapa adapted to LTS by changing protein function in the glyoxylate and dicarboxylate metabolism [34]. The significant enrichment of these pathways indicated that they were all necessary for the LTS response of C. album. Carbon fixation in photosynthetic organisms, tyrosine metabolism, ABC transporters, alanine, and aspartate and glutamate metabolism have been confirmed to participate in the LTS responses of several plant species, such as Lilium davidii [35], Oryza sativa [36], Panax ginseng [37], and Lolium perenne [38]. Noteworthily, the DEGs involved in these pathways were also found to be enriched in the cold-tolerant cultivar RA. Additionally, when compared with CK, the number of DEGs in RA was significantly lower than in QL after LTS. The enrichment of these pathways in cold-treated RA and much milder transcriptome changes of RA after LTS support well the high cold-tolerance of this cultivar.
The accumulated evidence has proved that phytohormones play important roles in plant's low-temperature responses. In our previous study, the auxin signal pathway in C. album cv. 'Fulan 1' significantly changed under −3 • C stress but not under chilling stress [13]. In this study, after LTS, the auxin signaling pathway was significantly induced, many involved genes exhibited differential expression, including AUX1, TIR1, AUX/IAA, ARF, GH3, SAUR, etc. Besides, using sRNA-seq, a series of auxin signaling or inductionrelated DEMs, including miR159 [39], miR164 [40], miR398 [40], and miR482 [41], were also identified. These results suggest that the auxin signaling might play important role in the response to the low-temperature stress of C. album. Salicylic acid promotes cell heat production, inhibits ROS accumulation, and alleviates cell dehydration after LTS [42]; the significant enrichment of salicylic acid signal transduction after LTS in C. album in our research suggests that this might be close related with its high cold tolerance.

The Interactions of Some MiRNAs and Their Target mRNAs Might Contributed Greatly to the Cold Resistance Response of C. album
In our study, we found that the predicted target genes of some DEMs were simultaneously identified as DEGs, indicating that these miRNAs and their corresponding target genes jointly regulate the response process through interaction. The differentially expressed target genes of DEMs (i.e., SCL, LAC, ARF, U-box protein, SNC, TKT, HIPP, 4CL, E3 ubiquitinprotein ligase, etc.) were involved in plant stress-resistance, auxin, and phenylpropane and flavonoid metabolism. Previous studies have shown that the knockout of miR171i or overexpression of its target gene SCL26.1 can improve the abiotic stress resistance of Malus domestica by regulating the expression of antioxidant enzyme genes [43]. Laccase regulates cell wall synthesis [44]; miR397-3 and its several predicted target LAC genes were differentially expressed in RA after LTS in this research, indicating that these genes may be regulating cell wall formation in RA to response to low temperatures via thickening or improving the compression. ARF members respond to the induction of LTS in T. aestivum [45] and Sorghum bicolor [46]. Sixteen MtPUB members of the U-box protein family of M. sativa were reported to be LTS responsive [47]. TKT is an important enzyme in the pentose phosphate pathway that was also found to be significantly related to response to abiotic stress; the LTS tolerance of transketolase antisense transgenic cucumber plants are decreased [48]. HIPP is a key protein for the safe transport of metal ions in cells and the expression of O. sativa HIPP41 increased greatly under cold stress [49]. 4CL is a key enzyme in phenylpropane metabolism, and research on Ocimum basilicum, Amygdalus persica, and B. rapa shows that it significantly affects plants' response to low temperature [50][51][52]. Specifically, SNC1 was very important for plants' adaptation to anomalous temperature stress. Previous research has shown that the activities of superoxide dismutase and other enzymes of the snc1 Arabidopsis mutant increase to clear the ROS and help the plant to adapt the oxidative stress environment. In our present study, a DEG encoding SNC1 was predicted to be one of the targets of low-temperature-inducible cal-miR482-22 in C. album. It was thus predicted that the downregulation of SNC1 gene might be helpful in improving the antioxidant reductase activity to eliminate the accumulation of ROS caused by low temperature in RA and improve its cold resistance [53]. Our results provide new insights into the underlying mechanisms of cold-resistance in C. album.

Conclusions
In conclusion, the basic expression levels of some miRNAs (especially the high expression of cal-miR397-3, the low expression of cal-miR398_2-3, and the specifically-induced expression of cal-miR482-22) in cold-tolerant cultivar of C. album, as well as their differentially expressed target genes under low temperature stress, contribute greatly to its cold-tolerance characteristics.