Long Non-Coding RNA and Its Regulatory Network Response to Cold Stress in Eucalyptus urophylla S.T.Blake

: Long non-coding RNA (lncRNA) plays an important regulatory role in plant growth and development, but its systematic identiﬁcation and analysis in Eucalyptus has not yet been reported. Cold stress has a huge impact on the survival and yield of Eucalyptus seedlings, but the regulatory mechanism of lncRNA in Eucalyptus in response to cold stress is still unclear. In this study, the transcriptomes of young leaves of Eucalyptus urophylla S.T.Blake under low-temperature treatment and restoration were analyzed by RNA-seq. A total of 11,394 lncRNAs and 46,276 mRNAs were identiﬁed, of which 300 were differentially expressed lncRNAs (DE_lncRNAs) and 5606 were differentially expressed target genes of lncRNAs under cold stress, with the total number of target genes of DE_lncRNAs being 1681. A total of 677 differentially expressed transcription factors (TFs) were also identiﬁed, mainly including ERF , MYB and the NAC transcription factor family. Gene ontology (GO) analysis of the differentially expressed genes (DEGs) and target genes of DE_lncRNAs was mostly related to the response to cold stress and external stimuli. Furthermore, lncRNA–miRNA– mRNA regulatory networks were constructed, and 22 DE_lncRNAs were predicted to be targets or targeting mimics of 20 miRNAs. A qRT-PCR was used to verify the relative expression of genes in the regulatory EuGBF3 - EUC_00002677 - MSTRG.7690 network, and it matched the transcriptome data, indicating that it may play an important role in the response to cold stress in E. urophylla . This study provides a new insight into lncRNA and its regulatory network under abiotic stress, especially cold stress in E. urophylla .


Introduction
Eucalyptus is the general term for Eucalyptus L'Herit in the Myrtaceae family. It is usually a perennial evergreen tree, mainly distributed naturally in Australia and its northern islands [1]. Eucalyptus is the most important plantation tree species in the world, providing a large amount of wood for industrial production [2,3] and sufficient raw material for pulping and papermaking [4][5][6]. Eucalyptus oil contained in leaves also has a variety of biological activities and is an important and effective ingredient in pesticides and herbicides [7,8]. As the most significant source of forest products in the world, Eucalyptus has been widely planted in more than 200 countries or regions, with a cumulative planting area of more than 20 million hectares [9,10]. Since it mostly originates in areas with hot climates, Eucalyptus does not tolerate low temperatures well. Cold stress limits the distribution and productivity of Eucalyptus plantations. When the temperature drops to 8 • C or below, Eucalyptus globulus exhibits various symptoms of chilling injury due to its inability to adapt to cold stress [11]. Especially after afforestation in mountainous regions in spring, occasional low temperatures may damage the buds and young leaves of seedlings, ultimately causing afforestation to fail. Therefore, research to reveal the cytological and molecular mechanisms of Eucalyptus in response to cold stress can provide important theoretical guidance for its future resistance breeding by molecular design or gene editing.
Cold stress is a kind of abiotic stress which is usually not conducive to plant growth and development [12]. In recent years, the signal pathways in plants that respond to external stress have been revealed. Many studies have shown that plants can respond to a short period of low but non-freezing temperatures through a complex transcriptional regulatory network which increases the plant's cold tolerance. This process is called cold acclimation [13]. It changes the physiological state by altering the membrane lipid composition, photosynthetic pigments, affinity permeants (such as soluble sugars, betaine and proline) and antioxidant levels and by causing extensive rearrangement of gene expression [14][15][16][17][18][19]. The physiological and biochemical changes that occur during plant cold acclimation are mainly due to the fluctuations in cold-responsive gene expression [19]. For example, atypical R2R3 MYB transcription factors positively regulate cold hardiness and cold-responsive gene expression under cold stress and then increase the cold tolerance of apples through CBF-dependent and CBF-independent pathways [20]. After low-temperature treatment of two pepper varieties with different cold tolerances, ERF family genes had a specific regulatory effect on the cold response pathways through transcriptome sequencing [21]. miRNAs and lncRNAs are both involved in cold stress responses. For instance, overexpressing miR397a in Arabidopsis thaliana (L.) Heynh. increases the tolerance to cold stress [22], and a cold-responsive lncRNA named SVALKA plays a role in regulating CBF1 expression [23]. However, there is no detailed study on how lncRNA interacts with mRNA and miRNA to respond to cold stress.
In the past decade, many studies have shown that the genomes of plants are generally transcribed, while some transcripts do not encode proteins, producing a large number of long non-coding RNAs [24]. Long non-coding RNAs (lncRNAs) characterize a highly heterogeneous class of RNA molecules defined as transcripts, with over 200 nucleotides not translated into proteins [25]. lncRNAs contain natural anti-sense RNA (NATs) of proteincoding genes, small RNA precursors (microRNA, small nucleolar RNA/SnoRNA and other small regulatory RNA) and structural RNAs (tRNA and ribosomal RNA) [26]. They can originate from exonic, intronic, intragenic, intergenic, promoter regions, 3 -and 5 -UTR or enhancer sequences transcribed in either the sense or the antisense strand [27]. Based on genomic locations, lncRNAs can generally be divided into sense lncRNAs, antisense lncRNAs, bidirectional lncRNAs, intergenic lncRNAs and intronic lncRNAs [27]. Long non-coding RNAs act as significant managers in gene expression networks by governing structure and transcription in the nucleus and modulating mRNA steadiness, translation and post-translational modifications in the cytoplasm [28]. Because of the continuous developments in high-throughput sequencing technology, lncRNAs of various species have been identified, such as Arabidopsis [29], rice [30], maize [31], poplar [32] and Eucalyptus grandis W.Hill [33]. lncRNA in plants are involved in responses to salinity, osmotic stress, drought, heat stress and other stress. For example, lncRNAs are differentially and specifically expressed in roots and leaves in response to salinity in A. thaliana [34]. There are also studies on the response of lncRNA to low-temperature stress. For example, some cold-responsive lncRNAs were identified in Medicago truncatula Gaertn. seedlings, and potential regulatory networks between lncRNA MtCIR1 and MtCBFs were discovered [35]. However, no study has provided a comprehensive and detailed explanation of the role of lncRNAs in regulating the response to cold stress in plants.
Recently, many studies have investigated the physiological and gene expression changes in Eucalyptus in a low-temperature environment. Oberschelp et al. [36] conducted metabolomics and proteomics studies to compare cold acclimation and freezing tolerance in three Eucalyptus species and found that E. benthamii showed a higher cold tolerance than E. dunnii and E. grandis, possibly due to its higher accumulation of phenolics, anthocyanins and soluble sugars, as well as lower levels of photosynthetic pigments and related proteins. Gaete-Loyola et al. [37] used transcriptome analysis to reveal the molecular mechanism of cold acclimation and deacclimation in E. nitens, providing a comprehensive view of the total gene expression and revealing meaningful information about the dynamic and complex nature of gene expression.
Eucalyptus urophylla S.T.Blake is widely distributed in disjunct populations that occupy the islands of Indonesia and Timor [38]. Eucalyptus urophylla is fast-growing and easy to cross-breed. It is not only a vital cultivated tree species in China and other Asia-Pacific regions, but it is also an important parent of high-yielding hybrids such as Eucalyptus urophylla x E. grandis and Eucalyptus grandis x E. urophylla. This species grows well and is widely planted in China, especially southern China. It provides an indispensable source for the production of paper, pulp, cellulose products and lignocellulosic biofuels. The study revealed that the cold resistance mechanism of E. urophylla is important for the genetic improvement of Eucalyptus.
However, lncRNAs and mRNAs under cold stress in E. urophylla have not been identified and characterized, and the mechanisms by which lncRNA and mRNA participate in cold tolerance are still unclear. Therefore, this study subjects E. urophylla to low-temperature treatment to simulate damage of its young leaves by cold stress and uses next-generation DNA sequencing (NGS) to perform RNA-seq, exploring the response of lncRNA and mRNA to cold stress. Furthermore, an lncRNA-miRNA-mRNA regulatory network under cold stress in E. urophylla is constructed, and it provides a new perspective for Eucalyptus to better adapt to abiotic stress.

Plant Materials
Eucalyptus urophylla plants preserved in the Guangxi Dongmen Forest Farm were used. In vitro plantlets produced by tissue culture were transplanted to pots and cultivated for about 2 months in a greenhouse under natural light conditions (1250 µmol m −2 s −1 of photosynthetically active radiation) at 25 ± 1 • C and 50 ± 1% relative humidity, with a 16/8 h light/dark cycle. For cold treatment, the seedlings were transferred to a chamber at 4 • C for different time periods: LT6, LT12, LT24 and LT48, referring to low-temperature treatment for 6, 12, 24 and 48 h, respectively, and RE24, referring to returning to a normal temperature of 25 • C for 24 h after 48 h of low-temperature treatment. Some seedlings were moved out of the chamber and used as the non-treatment control group (CK). After each treatment, young leaves were collected, immediately frozen in liquid nitrogen and stored at −80 • C until use. Three biological replicates were set for each treatment and sampling.

Measurement of Fv/Fm and Relative Electrolyte Leakage
After 20 min of dark adaptation, the first through third fully expanded young leaves were used to measure the maximal quantum efficiency of photosystem II (PSII) (Fv/Fm) [39]. These values were automatically monitored and calculated using a chlorophyll fluorometer (MultispeQ, East Lansing, MI, USA). Subsequently, 0.1 g of a young leaf cut into small pieces of 0.5 cm 2 was used to measure the relative electrolyte leakage (REL) using a conductivity meter (Leici-DDS-307A, Shanghai, China). After 6 h of shock at 160 rpm, leaf samples packed in the centrifuge tube containing 10 mL of deionized water were used to determine the initial conductivity (REL1). Then, these leaf samples were boiled for about 20 min, and the final conductivity (REL2) was measured after cooling to a normal temperature. The REL value was calculated as (REL1/REL2) × 100% [40].

RNA Extraction and Sequencing
The total RNA was extracted and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), following the manufacturer's instructions. The total RNA quantity and purity were analyzed using Bioanalyzer 2100 and the RNA 1000 Nano LabChip Kit (Agilent, CA, USA) with an RIN >7.0. About 5 µg of the total RNA was used to deplete the ribosomal RNA according to the manual of the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, USA). The enriched circRNAs were broken into small pieces with divalent cations in a high-temperature environment. After removing ribosomal RNAs, the remaining RNAs were reverse-transcribed to create cDNA, which was used to synthesize U-labeled, secondstranded DNAs with E. coli DNA polymerase I, RNase H and dUTP. An A-base was then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the Atailed fragmented DNA. Single-or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme treatment of the U-labeled, second-stranded DNAs, the ligated products were amplified by PCR under the following conditions: initial denaturation at 95 • C for 3 min, 8 cycles of denaturation at 98 • C for 15 s, annealing at 60 • C for 15 s, extension at 72 • C for 30 s and final extension at 72 • C for 5 min. The average insert size for the final cDNA library was 300 bp (±50 bp). Finally, paired-end sequencing was performed on an Illumina Novaseq™ 6000, following the vendor's recommended protocol.

Transcript Assembly
Cutadapt [41] was used to remove the reads that contained adapter contamination, low-quality bases and undetermined bases. The clean reads obtained were aligned to the Eucalyptus genome using Hisat2 [42]. Transcriptome assembly and expression quantification were performed using StringTie [43].

Identification of lncRNA
Assembled transcripts with overlapping of known mRNAs or that were shorter than 200 bp were discarded. CPC and CNCI were used to predict the coding potential of the remaining transcripts with default parameters. The transcripts with a CPC score ≥−1 and a CNCI score ≥0 were considered lncRNAs.

Different Expression Analysis of mRNAs and lncRNAs
Edger R was used to identify differentially expressed mRNA or lncRNA genes. Genes with |log 2 (fold change)| ≥ 1 and p < 0.05 were considered differentially expressed genes.

Target Gene Prediction of lncRNAs with Protein-Coding Genes
In this study, coding genes in 100,000 upstream and downstream components were selected by a Python script.

Functional Analysis of DEGs and DE_lncRNAs:
Functional annotation of differentially expressed lncRNA loci was based on their target genes. The target genes and differentially expressed genes were annotated based on the GO database (http://www.geneontology.org; accessed on 12 July 2020).

Identification of Transcription Factors
The Plant Transcription Factor Database was used to identify and classify TFs from differentially expressed genes.

Prediction of lncRNAs as miRNA Target Mimics
Eucalyptus grandis miRNA sequences were downloaded from a previous study [44]. The target mimics were predicted by submitting the lncRNAs and miRNAs to psRNATarget (http://plantgrn.noble.org/psRNATarget; accessed on 30 July 2020). Less than three mismatches and G/U pairs were allowed within the lncRNA and miRNA pairing regions [45]. Cytoscape was used to establish the lncRNA-miRNA-mRNA interaction network.

Quantitative Real-Time PCR
The total RNA was isolated from the young leaves of E. urophylla using an RNAprep Pure Plant Plus Kit (Tiangen, China, cat DP441). About 1.5 µg of RNA was reversetranscribed into first-stand cDNA using a cDNA synthesis kit (Tiangen, China, cat KR106). qPCR was subsequently performed using a TransStart ® Tip Green qPCR SuperMix (Trans-Gen Biotech, Beijing, China) with EuGADPH as a reference gene. qPCR on the Applied Biosystems 7500 FAST real-time PCR system according to the manufacturer's instructions. The primers used for qRT-PCR analysis are listed in Table S1.

Statistical Analysis
The Statistical Product and Service Solutions (SPSS) 20 program was used to conduct all statistical analyses of data, and data were presented as the mean ± SD (standard deviation). The significant differences between means (p < 0.05) were examined by a t-test and analysis of variance (ANOVA). The Waller-Ducan test was used for post hoc comparison after ANOVA.

The Effect of Cold Stress on the Phenotype of Eucalyptus urophylla
To analyze the effect of cold stress on E. urophylla, the leaf morphology, relative electrolyte leakage (REL) and maximum photochemical efficiency of PSII (Fv/Fm) under different durations of low-temperature treatment (4 • C) were carefully observed. Compared with the control group (CK), after 6, 12, 24 and 48 h of low-temperature treatment (LT6, LT12, LT24 and LT48, respectively), the young leaves had a different extent of wilting ( Figure 1a). In the first 6 h of low-temperature treatment (LT6), the morphological changes in the young leaves were not obvious. As the duration of low-temperature treatment increased from 12 h to 48 h (LT12, LT24 and LT48), the young leaves showed different levels of curling and wilting. After 48 h of treatment, the plants were placed in a normaltemperature environment for 24 h of recovery (RE24), but the leaves did not return to normal. The REL value is a key indicator of the severity of cell penetration. Compared with low temperature treatment for 6 h, the REL value when treated for 12 h did not change significantly ( Figure 1c). After that, they gradually increased from 12.67% to 57.79%, especially in the 48 h and later treatment groups (LT48 and RE24), and they were significantly higher than in other treatments. Chlorophyll fluorescence is a key parameter for evaluating photosynthesis. The Fv/Fm value is used to represent the maximum photochemical efficiency of PSII, which is often used to determine the response degree of plants to various stress treatments. Compared with the control group (CK), Fv/Fm decreased from 0.68 to 0.56 in the initial 6 h of low-temperature treatment (LT6) ( Figure 1b) and then rose slightly after adaptation at 12 h of treatment (LT12). With the increase in the cold stress duration, Fv/Fm further gradually decreased to 0.53 at 48 h of treatment (LT48). After that, with 24 h of recovery (RE24), Fv/Fm directly decreased to 0.09, which was significantly lower than in other treatments. The phenotype of E. urophylla underwent major changes under different durations of low-temperature treatment.

Genome-Wide Identification and Characterization of lncRNAs and mRNAs of E. urophylla under Cold Stress
Young leaves of E. urophylla from six groups were collected, including one control group (CK), four low-temperature treatment groups chilled at 4 • C for 6, 12, 24 and 48 h (LT6, LT12, LT24 and LT48, respectively) and one recovery group (recovered to a normal-temperature environment for 24 h after 48 h of a low-temperature treatment (RE24)). RNA sequencing from 18 cDNA libraries with three biological replicates from these six groups was conducted to systematically identify the lncRNA response to cold stress. Low-quality reads were abandoned from the RNA-seq data, and adapter sequences from residual reads were trimmed to obtain clean reads (Table 1)

Genome-Wide Identification and Characterization of lncRNAs and mRNAs of E. urophylla under Cold Stress
Young leaves of E. urophylla from six groups were collected, including one control group (CK), four low-temperature treatment groups chilled at 4 °C for 6, 12, 24 and 48 h (LT6, LT12, LT24 and LT48, respectively) and one recovery group (recovered to a normaltemperature environment for 24 h after 48 h of a low-temperature treatment (RE24)). RNA sequencing from 18 cDNA libraries with three biological replicates from these six groups Values shown in (b,c) are means ± SD with three biological replicates. The letters indicate significant differences (one way ANOVA followed by post hoc Waller-Duncan test, p < 0.05).
A series of screening processes was followed to identify the lncRNAs from CK, LT6-LT48 and RE24 groups ( Figure 2). Numerous candidate lncRNAs were predicted using the Coding Potential Calculator (CPC) and Non-Coding Index (CNCI) databases, including 433 intronic lncRNAs, 188 bidirectional lncRNAs, 1429 sense lncRNAs, 7407 intergenic lncRNAs and 1937 antisense lncRNAs (Figure 3a). The mRNAs and lncRNAs involved in cold stress were characterized according to five aspects: the exon number, ORF length, isoform number per gene, transcript length and differential expression level. Most of the mRNAs and lncRNAs encompassed one or two exons, accounting for 34% and 91%, respectively. The maximum number of exons of lncRNAs was five, while the number of exons in mRNAs was in a wider distribution range ( Figure 3b). The whole transcript length and ORF length of the lncRNAs were apparently much shorter than those of the mRNAs. The number of lncRNAs with ORF lengths less than 100 nt was 2081, occupying more than 83%, while the number of mRNAs with ORF lengths less than 100 nt was 2902. The ORF lengths of the mRNAs were primarily distributed between 100 and 500 nt. The longest ORF length of the mRNAs was more than 1200 nt, while the longest ORF length of the lncRNAs was shorter than 600 nt (Figure 3c). In terms of the transcript length, most of the mRNA transcripts were longer than 1000 nt, accounting for about 69%, while lncRNA transcripts not exceeding 300 nt accounted for a maximum ratio of 46% (Figure 3e). The expression levels of the mRNAs and lncRNAs were also compared using the fragments per kilobase of exon per million fragments mapped (FPKM). The lncRNAs and mRNAs were expressed at similar levels in all six groups (Figure 3f). However, the overall expression levels of the mRNAs were generally higher than those of lncRNAs. Around 8% of the lncRNAs were alternatively spliced, while about 23% of the mRNAs were spliced (Figure 3d).       Three lncRNAs and three mRNAs from the treatment and recovery groups were randomly selected and analyzed by qRT-PCR in order to explore their expression patterns. The results showed that the expression patterns of lncRNAs and mRNAs involved in cold stress were relatively consistent and comparable (Figure 4), suggesting that their expression based on RNA-seq data was reliable.

Identification of DEGs and DE_lncRNAs
To understand gene expression under cold stress, the transcript levels of the lncRNAs and mRNAs in terms of fragments per kilobase of exon per million fragments mapped (FPKM) based on RNA-Seq data were measured ( Figure S1, Tables S2 and S3). The conditions to select differentially expressed mRNAs (DEGs) and differentially expressed lncRNAs (DE_lncRNAs) were set as |log2 (fold change)| ≥ 1 and p < 0.05. The heat maps of the mRNAs and lncRNAs demonstrated gene expression after normalizing the data of all six groups. The results showed significant differences in gene expression in E. urophylla under cold stress, especially comparing the CK group with the LT6, LT12, LT24, LT48 and RE24 groups, indicating a rapid response to low-temperature stimulus of gene regulation. The number of DEGs (DE_lncRNAs) of five comparison groups (LT6 vs. CK, LT12 vs. CK, LT24 vs. CK, LT48 vs. CK and RE24 vs. CK) were 1893 (57), 2108 (85), 3076 (159), 4707 (232) and 3123 (45), respectively. The expression patterns of the mRNAs and lncRNAs in each treatment group were generally similar, except in the LT6 and LT12 groups (Figure 5c,d). Compared with the control group, with the increase in cold stress time, the number of DEGs and DE_lncRNAs gradually increased. After 24 h of recovery, the number of DEGs and DE_lncRNAs decreased to a certain extent. Pairwise comparisons were performed among the treatment groups and the recovery group, and a large number of DEGs and DE_lncRNAs were obtained. Compared with the LT48 group, the RE24 group had 4359 DEGs and 74 DE_lncRNAs, suggesting the recovery process had a vital influence on the biological process in E. urophylla after low-temperature treatment. Detailed analysis of differentially expressed genes showed that the number downregulated DEGs was more than the upregulated DEGs with the increase in treatment time (Figure 5a), while the upregulated DE_lncRNA were mostly greater than the downregulated ones (Figure 5b). These results show that cold stress caused different changes in the mRNAs and lncRNAs in E. urophylla. Figure 4. Confirmation of the expression patterns of three lncRNAs and three mRNAs of Eucalyptus urophylla under cold stress using qRT-PCR. The first layer is three randomly selected mRNAs, and the second layer is three lncRNAs. The values shown are the means ± SD of three replicates. EuGADPH was used as the reference gene.

Identification of DEGs and DE_lncRNAs
To understand gene expression under cold stress, the transcript levels of the lncRNAs and mRNAs in terms of fragments per kilobase of exon per million fragments mapped (FPKM) based on RNA-Seq data were measured ( Figure S1, Tables S2 and S3). The conditions to select differentially expressed mRNAs (DEGs) and differentially expressed lncRNAs (DE_lncRNAs) were set as |log 2 (fold change)| ≥ 1 and p < 0.05. The heat maps of the mRNAs and lncRNAs demonstrated gene expression after normalizing the data of all six groups. The results showed significant differences in gene expression in E. urophylla under cold stress, especially comparing the CK group with the LT6, LT12, LT24, LT48 and RE24 groups, indicating a rapid response to low-temperature stimulus of gene regulation. Compared with the control group, with the increase in cold stress time, the number of DEGs and DE_lncRNAs gradually increased. After 24 h of recovery, the number of DEGs and DE_lncRNAs decreased to a certain extent. Pairwise comparisons were performed among the treatment groups and the recovery group, and a large number of DEGs and DE_lncRNAs were obtained. Compared with the LT48 group, the RE24 group had 4359 DEGs and 74 DE_lncRNAs, suggesting the recovery process had a vital influence on the biological process in E. urophylla after low-temperature treatment. Detailed analysis of differentially expressed genes showed that the number downregulated DEGs was more than the upregulated DEGs with the increase in treatment time (Figure 5a), while the upregulated DE_lncRNA were mostly greater than the downregulated ones (Figure 5b). These results show that cold stress caused different changes in the mRNAs and lncRNAs in E. urophylla.

Prediction of Target Genes of lncRNAs
Since lncRNA plays an important role in regulating gene function, identification and analysis of their target genes could be helpful to explore the function of lncRNAs. In this study, the number of target genes was 19,418, of which 5606 were differentially expressed target genes, while 300 were DE_lncRNAs and 579 were differentially expressed target genes of DE_lncRNA (Table 2).

Identification of Transcription Factors in E. urophylla
A total of 677 differentially expressed TFs under cold stress were identified in each group (Table S4). These TFs belonged to 48 families, of which the top 12 TF families were ERF, MYB, NAC, bHLH, C2H2, WRKY, GRAS, HD-ZIP, G2-like, bZIP, Dof and MYB_related, consisting of 77, 66, 64, 53, 42, 41, 26, 23, 22, 21, 21 and 21 genes, respectively. Most of these gene families are related to the response to cold stress [46] (Figure 6). Among them, the TF family with the maximum differentially expressed genes was ERF, another subfamily of the AP2/ERF TF family that plays a vital role in the biological regulation of the response to stress [47,48]. Among these differentially expressed TFs, 14 were the target genes of DE_lncRNAs, belonging to the ERF, MYB_related, NAC, C3H, CO_like and WRKY TF families (Table S4). Of these 14, 4 differentially expressed TFs were the target genes of MSTRG.6543, namely Eucgr.B02398, Eucgr.B02399, Eucgr.B02396 and Eucgr.B02393, which all belonged to the C2H2 TF family. target genes, while 300 were DE_lncRNAs and 579 were differentially expressed target genes of DE_lncRNA (Table 2). Table 2. Statistical data of annotation lncRNAs and target genes.

Identification of Transcription Factors in E. urophylla
A total of 677 differentially expressed TFs under cold stress were identified in each group (Table S4) (Figure 6). Among them, the TF family with the maximum differentially expressed genes was ERF, another subfamily of the AP2/ERF TF family that plays a vital role in the biological regulation of the response to stress [47,48]. Among these differentially expressed TFs, 14 were the target genes of DE_lncRNAs, belonging to the ERF, MYB_related, NAC, C3H, CO_like and WRKY TF families (Table S4)

Functional Analysis of DEGs and Target Genes of DE_lncRNAs
Gene ontology (GO) was used to evaluate the functions of DEGs and was divided into three parts: biological process, cell component and molecular function. To further analyze the effect of cold stress on the gene expression of E. urophylla, GO enrichment analysis of the DEGs and target genes of DE_lncRNAs was performed (Figure 7). The DEGs were mainly enriched in four GO terms in the cell component process under lowtemperature stress in E. urophylla: anchored component of membrane, kinesin complex, nucleosome and cell surface. In general, the DEGs were mainly enriched in parts of the

Functional Analysis of DEGs and Target Genes of DE_lncRNAs
Gene ontology (GO) was used to evaluate the functions of DEGs and was divided into three parts: biological process, cell component and molecular function. To further analyze the effect of cold stress on the gene expression of E. urophylla, GO enrichment analysis of the DEGs and target genes of DE_lncRNAs was performed (Figure 7). The DEGs were mainly enriched in four GO terms in the cell component process under low-temperature stress in E. urophylla: anchored component of membrane, kinesin complex, nucleosome and cell surface. In general, the DEGs were mainly enriched in parts of the biological processes and molecular functions. Most of the enrichment pathways of the biological processes were related to the cold stress response, such as the response to cold being enriched within 134 genes. The pathway with most DEGs enriched was the response to salt stress, which was enriched within 164 genes. The results also show that the pathways including the response to absicic acid, response tosalicylic acid, enthylene-activated signaling pathway, response to ethylene and response to jasmonic acid were enriched within 160, 100, 99, 62 and 61 genes, respectively. In the GO terms of molecular functions, the most enriched pathways were reansferase activity, transferring glycosyl groups, protein kinase activity,  The GO enrichment of the target genes of DE_lncRNAs was similar to that of DEGs, with only one pathway, the cortical actin cytoskeleton, enriched in the GO terms of cellular components, which were also mainly enriched in biological processes and molecular functions. Most of the biological process enrichment pathways were related to cold stress and external stimuli, including the response to cold, cellular response to extracellular stimulus, cellular response to external stimulus, response to extracellular stimulus, cellular response to abiotic stimulus, cellular response to environmental stimulus and cellular response to cold. The cold-related and drought-and water-deficiency pathways were also enriched, such as the response to water deprivation, response to water and cellular response to water deprivation. The target genes of DE_lncRNAs were also enriched in the area of molecular functions, and the number of pathways enriched within more than 20 genes was 6. The most enriched pathway was protein homodimerization activity, containing 68 genes, while the second-and third-most enriched pathways were transition metal ion binding and lipase activity, including 41 and 38 genes, respectively.

DEG, DE_lncRNA, Differentially Expressed TF and miRNA Interaction Network Construction
A large number of lncRNAs that could be used as targets or target mimics of miRNA were recognized. Among these 300 cold stress responsive DE_lncRNAs, 22 were identified as targets or target mimics of 20 miRNAs, belonging to the miR166, miR167 and miR482 families containing 27 targeting relationships (Table S5). In the prediction of the miRNA targeting relationship, the results also show the interaction between the target genes of the DE_lncRNAs and miRNAs. A total of 11 target genes of the DE_lncRNAs were identified as target mimics of 14 miRNAs.
To analyze the interaction between lncRNA and miRNA, transcriptional regulatory networks among the DEGs, DE_lncRNAs, differentially expressed TFs (DE_TFs) and miRNAs under low-temperature stress were constructed in E. urophylla (Figure 8). One of the regulatory networks showed that EUC_00002677 from the miR167d family had four target genes ( Figure 8I), one of which was a DEG named Eucgr.J01413, while the others were all TF, including EuGNC (Eucgr.C03895), EuFRSII (Eucgr.K03522) and EuGBF3 (Eucgr.E000078), which belong to the GATA, FAR1 and bZIP TF families, respectively. EUC_00002677 can also form a targeting relationship with two lncRNAs, MSTRG.48389.1 and MSTRG.7690.1, which respond to cold stress. Another regulatory network showed that ( Figure 8VI) miRNA EUC_00008890, belonging to the miR167d family, can target the Dof zinc finger gene EuDAG2 (Eucgr.B0245) and lncRNA MSTRG.46425, which respond to cold stress. The third regulatory network showed that MSTRG.26750 has two target DEGs, and it forms a targeting relationship with two miRNAs, EUC_00005453 (miR482e) and EUC_00006406 (miR447B-3P). Furthermore, EUC_00005453 (miR482e) could also bind to four target mRNAs, all of which were TFs, including three members from the NAC family and one from the SRS family ( Figure 8II).
In addition, qRT-PCR was used to determine the relative gene expression in the mRNA-miRNA-lncRNA regulatory network of E. urophylla under low-temperature treatment to verify their interaction. EUC_00002677 from the miR167d family could target an lncRNA named MSTRG7690 and a transcription factor named EuGBF3 (Eucgr.E00078) from the bZIP family. The relative expression of the above three genes was measured to explore the relevance of their changing trends. The expression of miRNA EUC_00002677 was downregulated under low-temperature treatment, which promoted the upregulation of lncRNAs MSTRG.7690 and EuGBF3. The change was particularly obvious after 6 h of low-temperature treatment (LT6). The upregulated expression of MSTRG.7690 and EuGBF3 was nearly 9 and 1.5 times that of the control group, respectively, which indicates that lncRNAs may be more sensitive to cold stress. The expression of lncRNAs MSTRG.7690 and EuGBF3 was downregulated after 12 h of treatment (LT12), but overall, it was higher in the LT6, LT12, LT24 and LT48 groups than in the control group, while the expression of miRNA EUC_00002677 was lower than in the control group. In addition, at 24 h of recovery after the low-temperature treatment, EUC_00002677 was upregulated, and its expression was approximately twice that of the control group, while the expression of MSTRG.7690 and EuGBF3 was downregulated and was lower than that of the control group. In general, under cold stress, the expression of EUC_00002677 in E. urophylla was downregulated, which promoted the upregulation of MSTRG.7690 and EuGBF3, showing a negative correlation, which was in line with the ceRNA theory.

Discussion
Many studies have shown that lncRNAs can interact with DNA, RNA or TFs to alter the expression of target genes in order to regulate different biological processes [49]. With the continuous development of high-throughput sequencing technology, lncRNAs of many species of plants have been identified, such as Arabidopsis thaliana [29], rice [30], corn [31] and poplar [32]. Although some lncRNAs have been reported in Eucalyptus [33], the lncRNA response to cold stress and its mRNA-miRNA-lncRNA regulatory network have not been identified in E. urophylla. In this study, E. urophylla seedlings were subjected to lowtemperature treatments (4 • C) for different durations. Through transcriptome sequencing, 11,394 lncRNAs and 46,276 mRNAs were identified, of which 300 were DE_lncRNAs responding to low temperatures, and 5606 were garget genes of lncRNAs. The DE_lncRNAs and DEGs were mainly enriched in the abiotic stress response, hormone signal response and carbohydrate biosynthetic processes, among others. In the newly constructed mRNA-miRNA-lncRNA regulatory network, 22 DE_lncRNAs were identified as miRNA targets. The miRNAs in E. urophylla negatively regulate the expression of lncRNA and mRNA to re-spond to a low-temperature environment. In general, the regulatory network in E. urophylla under low-temperature stress was explored, which could provide a molecular basis for the research on resistance to cold stress and improve its adaptation to an external environment.
Plant responses to cold stress involve multifaceted processes, such as changes in physiology, metabolism and macromolecular processes [50]. Cold stress can cause obvious damage to plants, including wilting, drying of leaf edges, discoloration, incomplete ripening, accelerated aging and even death [51]. The phenotype changes significantly after cold treatment. Young leaves are more sensitive to cold stress than mature leaves [52]. Phenotypic changes in young leaves began to appear in the LT12 group, and the leaves began to wither in the LT24 group, indicating that low temperatures damage the cell membrane, especially in young leaves, and further cause water shortage [53]. In addition, the REL value increased after the cold treatment, which confirmed that low temperatures destroy the structures of the cell membranes and further cause a large amount of extravasation of soluble substances. The Fv/Fm value significantly reduced after cold treatment, indicating that low temperature can destroy photosynthetic components, thereby affecting the growth and development of plants. Both REL and Fv/Fm showed obvious changes after 6 h of cold treatment and a certain degree of recovery after 12 h of cold treatment, which may be related to the balance of oxidative stress and the antioxidant system in the plant [53,54]. Plants can induce reactive oxygen species (ROS) in a low-temperature environment [55][56][57][58][59][60]. To cope with the massive production of ROS, plants can run an efficient purification system, including increasing the activity of non-enzymatic antioxidants and antioxidant enzymes [61,62]. When the temperature was restored to room temperature for 24 h after low-temperature treatment, neither the leaf state nor the REL or Fv/Fm recovered, which indicated that low temperatures cause irreversible damage to plant leaves, especially young leaves.
To understand the effect of cold stress on the regulation of gene expression in E. urophylla, GO enrichment analysis of DEGs and the target genes of DE_lncRNAs was performed. In the cell components, the GO term with the maximum DEG enrichment was "anchored component of membrane", suggesting that the main effect of cold stress on the cells in E. urophylla is a change in the composition of the cell membrane, especially the saturation of constituent lipids [63], eventually causing young leaves to lose water and wither. The response of DE_lncRNAs to cold stress at the cellular level is mainly through the influence of the cortical actin cytoskeleton. In the biological process, in addition to the "response to cold" term, the DEGs and the target genes of DE_lncRNAs were mainly enriched in response to salt, drought, osmotic stress and plant hormones, indicating crosstalk in the response to different abiotic stresses in E. urophylla under cold stress. This is because when plants are under abiotic stresses, such as cold stress, signal transduction, such as the calcium ion signal transduction pathway, could play a corresponding role to create various responses to external stimuli [50,64]. Plant hormones also play an important role in the response to stimuli and regulate signal transduction pathways [65], such as jasmonic acid [66], ethylene [67] and abscisic acid [68]. Some DEGs were also enriched in the flavonoid biosynthetic process, because cold stress can regulate flavonoid synthesis pathways through transcriptome modification, including anthocyanin synthesis and metabolism pathways [69]. In addition, some DE_lncRNAs were enriched in the carbohydrate synthesis pathway. Notably, DE_lncRNAs related to the carbohydrate biosynthetic process under cold stress are likely beneficial to homeostasis between cellular membranes and ROS [70,71], which could partly explain the increase in Fv/Fm and the decrease in REL in the LT12 group compared with the LT6 group. Regarding molecular function, most of the DEGs were enriched in pathways related to enzyme activity, such as transferase activity, and glycosyl group transfer, which may be because E. urophylla accumulates carbohydrates through the activation of specific enzymes under cold stress [72,73], and these carbohydrates further act as osmotic regulators and signaling molecules [74]. The DE_lncRNAs were mostly enriched in protein homodimerization activity, which shows that cold stress had a significant impact on protein structure formation. These analyses indicate that large amounts of DEGs and DE_lncRNAs could respond to cold stress to regulate biological activities, which ultimately led to changes in the plant phenotype and physiological indicators.
TFs are proteins that affect many biological processes, including growth, development, cell division and responses to environment stimuli such as cold stress [75]. Transcriptome data showed that 677 transcription factors were differentially expressed under cold stress in E. urophylla. Of these, 14 TFs were identified as the target genes of DE_lncRNAs. Further investigation of the DE_TFs showed that 77 DEGs belonged to the ERF transcription factor family, which was the most enriched TF family, followed by the MYB and NAC TF families. APETALA2/ETHYLENE RESPONSIVE FACTOR (AP2/ERFs) are key regulators of several abiotic stresses [76,77], and CBF/dehydration-responsive element-binding factor (DREB1) from the AP2/ERF family (CBF1/DREB1B, CBF2/DREB1C and CBF3/DREB1A) in Arabidopsis plays a prominent role in cold acclimation [78,79]. Analyzing changes in the expression of EuCBF1/EuCBF2/EuCBF3 after low-temperature treatment showed that, compared with the CK group, the expression of EuCBF1-3 in E. urophylla after lowtemperature treatment significantly increased in the treatment groups, while its expression in the recovery group reduced compared with the treatment groups (Table S6). This indicates that the CBF transcription factor in E. urophylla could respond to cold stress. Low temperatures induce its expression, while normal conditions inhibit its expression.
In addition to the ERF family of TFs, MYB transcription factors also responded significantly to cold stress. We found 66 DEGs belonging to the MYB transcription factor family. MYB proteins are one of the largest TF families and are considered unique to plants. They are classified into four major subfamilies based on their structural features: 1R-MYB proteins, R2R3-MYB proteins, 3R-MYB proteins and 4R-MYB proteins [80]. They can promote or reduce the cold tolerance of plants by regulating the expression of downstream genes [81,82]. The NAC transcription factor family also showed significant differences in expression under cold stress in E. urophylla. All proteins in this family have a conserved Nterminal DNA-binding domain called the NAC structure domain [83]. Although most of the NAC genes are involved in regulating the salt and drought stress tolerance in plants, some are also involved in regulating temperature tolerance [46]. For example, overexpression of GmNAC20 and GmNAC11 in transgenic Arabidopsis enhances its salt and cold tolerance, respectively, due to the interactive effects of these genes with the DREB/CBF pathway and other stress-related genes [83]. Transcription factors themselves can be regulated at the transcription level by upstream components and also be subjected to various tiers of modifications at the post-transcription level, such as ubiquitination and sumoylation, thus forming a complex regulatory network to modulate the expression of stress-responsive genes [84,85]. Therefore, analyzing the roles of these TFs in various signal transduction pathways could help us better understand the mechanism of response to cold stress in E. urophylla.
Studies have shown that miRNA recognizes target mRNA through base pairing and participates in the degradation or inhibition of the silent complex that regulates the translation of target mRNA [86]. The competitive endogenous RNA (ceRNA) hypothesis proposed by Salmena et al. [87] is considered a new regulatory mechanism between non-coding RNA and mRNA. Based on the previous research on genome-wide miRNA identification, we analyzed the targeting relationship among miRNAs, DE_lncRNAs and mRNAs under cold stress to construct an mRNA-miRNA-lncRNA regulatory network in E. urophylla. In addition to the DE_lncRNAs identified in this study, miRNAs in E. urophylla also responded to cold stress. In the regulatory network constructed, miRNAs and lncRNAs constitute a targeting relationship, and 22 DE_lncRNAs are predicted to be targets or targeting mimics of 20 miRNAs. In this study, miRNA also showed differential expression under cold stress, which affected the expression of mRNA and lncRNA, indicating that cold stress can change the expression of miRNA and further regulate mRNA and lncRNA to respond to cold temperatures in E. urophylla (Figure 9). Cold stress can cause significant changes in the expression of these miRNAs. For example, miR167 in tomatoes [88] and miR193, miR169 and miR168 in Arabidopsis thaliana [89][90][91] are upregulated under low-temperature treatment, and the expression of miR166 in cabbage [92], miR482 in Hevea brasiliensis and miR858 in Citrullus lanatus change under cold stress [93,94]. miRNAs do not directly play a role in plant growth and development or in the response to abiotic stress but do participate in a plant's response to abiotic stress by regulating key components of a complex regulatory network. Both lncRNAs and mRNAs can act as target genes of miRNAs. miR167 is a highly conserved family acting in a wide range of biological processes in plants [95]. ARF6 and ARF8 are the targets of miR167 that promote adventitious rooting. Cold stress induces the significant differential expression of miR167 in many plant species [96]. In this study, we found new regulatory networks that involved miR167 in E. urophylla under cold stress. For example, miRNAs EUC_00002677 and EUC_00008890, both belonging to the miR167 family, responded to cold stress by interacting with the target genes in E. urophylla. EUC_00002677 targeted novel lncRNAs MSTRG.7690 and EuGBF3, and EUC_00008890 targeted novel lncRNAs MSTRG.46425 and EuDAG2. EuGBF3 and lncRNA MSTRG.7690 were both involved in the response to cold stress regulated by miRNA EUC_00002677, inducing the gene expression difference in E. urophylla ( Figure 9). DOF AFFECTING GERMINATION 2 (DAG2), which encodes Dof zinc finger transcription factors, is a positive regulator of light-mediated seed germination and is repressed by DAG1 [97]. This interaction network suggested that EuDAG2 and novel lncRNA MSTRG.46425 were both responsive to a low-temperature environment, which may be due to the regulation of miRNA EUC_00008890. In general, the mRNA-miRNA-lncRNA regulatory network we constructed revealed a new molecular mechanism of the response to cold stress in E. urophylla and may play a significant role in understanding the genetic improvement of cold tolerance in this tree species.

Conclusions
This study identified the lncRNAs' response to cold stress in E. urophylla and found that lncRNAs exhibit different functions. We also constructed mRNA-miRNA-lncRNA regulatory networks in E. urophylla under cold stress, which may provide new insights into enhancing the cold tolerance in this tree species.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1. Figure S1: Number of DE_lncRNAs and DEGs in the control group, four cold treatment groups and one recovery group; Table S1: List of primers for genes used for qPCR; Table S2: List of DEGs in Eucalyptus urophylla; Table S3: List of DE_lncRNAs and target genes in Eucalyptus urophylla; Table S4: Annotation of DE_TFs; Table S5: List of putative targets and target mimics of DE_lncRNAs for miRNA; Table S6: EuCBF gene expression under cold stress. EuGADPH was used as the reference gene. The values shown are the means ± SD of three replicates. Asterisks indicate statistically significant differences compared with the control group by a t-test (* p < 0.05; ** p < 0.01; *** p < 0.001). miRNAs do not directly play a role in plant growth and development or in the response to abiotic stress but do participate in a plant's response to abiotic stress by regulating key components of a complex regulatory network. Both lncRNAs and mRNAs can act as target genes of miRNAs. miR167 is a highly conserved family acting in a wide range of biological processes in plants [95]. ARF6 and ARF8 are the targets of miR167 that promote adventitious rooting. Cold stress induces the significant differential expression of miR167 in many plant species [96]. In this study, we found new regulatory networks that involved miR167 in E. urophylla under cold stress. For example, miRNAs EUC_00002677 and EUC_00008890, both belonging to the miR167 family, responded to cold stress by interacting with the target genes in E. urophylla. EUC_00002677 targeted novel lncRNAs MSTRG.7690 and EuGBF3, and EUC_00008890 targeted novel lncRNAs MSTRG.46425 and EuDAG2. EuGBF3 and lncRNA MSTRG.7690 were both involved in the response to cold stress regulated by miRNA EUC_00002677, inducing the gene expression difference in E. urophylla ( Figure 9). DOF AFFECTING GERMINATION 2 (DAG2), which encodes Dof zinc finger transcription factors, is a positive regulator of light-mediated seed germination and is repressed by DAG1 [97]. This interaction network suggested that EuDAG2 and novel lncRNA MSTRG.46425 were both responsive to a low-temperature environment, which may be due to the regulation of miRNA EUC_00008890. In general, the mRNA-miRNA-lncRNA regulatory network we constructed revealed a new molecular mechanism of the response to cold stress in E. urophylla and may play a significant role in understanding the genetic improvement of cold tolerance in this tree species.

Conclusions
This study identified the lncRNAs' response to cold stress in E. urophylla and found that lncRNAs exhibit different functions. We also constructed mRNA-miRNA-lncRNA regulatory networks in E. urophylla under cold stress, which may provide new insights into enhancing the cold tolerance in this tree species.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12070836/s1. Figure S1: Number of DE_lncRNAs and DEGs in the control group, four cold treatment groups and one recovery group; Table S1: List of primers for genes used for qPCR; Author Contributions: Conceptualization, H.C., J.L. and J.Y.; methodology, H.C., J.L. and X.K.; formal analysis, H.C., J.L. and J.Y.; investigation, J.L., B.Q., Y.Z. and Z.L.; writing-original draft preparation, H.C. and J.L.; writing-review and editing, J.Y. and X.K.; funding acquisition, J.Y. and X.K. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural Science Foundation of China (31901337).
Data Availability Statement: The data sets generated or analyzed during this study are included in this published article and its Supplementary Materials.

Acknowledgments:
The authors are grateful to Jianzhong Wang from Guangxi Dongmen Forest Farm for collecting materials and additional help.

Conflicts of Interest:
The authors declare no conflict of interest.