Genome-Wide Analysis of the miRNA–mRNAs Network Involved in Cold Tolerance in Populus simonii × P. nigra

Background: Cold tolerance is important for plants’ geographical distribution and survival in extreme seasonal variations of climate. However, Populus simonii × P. nigra shows wide adaptability and strong cold resistance. Transcriptional and post-transcriptional regulation of cold-responsive genes is crucial for cold tolerance in plants. To understand the roles of regulatory RNAs under cold induction in Populus simonii × P. nigra, we constructed cDNA and small RNA libraries from leaf buds treated or not with −4 °C for 8 h for analysis. Results: Through high-throughput sequencing and differential expression analysis, 61 miRNAs and 1229 DEGs were identified under cold induction condition in Populus simonii × P. nigra. The result showed that miR167a, miR1450, miR319a, miR395b, miR393a-5p, miR408-5p, and miR168a-5p were downregulated, whereas transcription level of miR172a increased under the cold treatment. Thirty-one phased-siRNA were also obtained (reads ≥ 4) and some of them proceeded from TAS3 loci. Analysis of the differentially expressed genes (DEGs) showed that transcription factor genes such as Cluster-15451.2 (putative MYB), Cluster-16493.29872 (putative bZIP), Cluster-16493.29175 (putative SBP), and Cluster-1378.1 (putative ARF) were differentially expressed in cold treated and untreated plantlets of Populus simonii × P. nigra. Integrated analysis of miRNAs and transcriptome showed miR319, miR159, miR167, miR395, miR390, and miR172 and their target genes, including MYB, SBP, bZIP, ARF, LHW, and ATL, were predicted to be involved in ARF pathway, SPL pathway, DnaJ related photosystem II, and LRR receptor kinase, and many of them are known to resist chilling injury. Particularly, a sophisticated regulatory model including miRNAs, phasiRNAs, and targets of them was set up. Conclusions: Integrated analysis of miRNAs and transcriptome uncovered the complicated regulation of the tolerance of cold in Populus simonii × P. nigra. MiRNAs, phasiRNAs, and gene-encoded transcription factors were characterized at a whole genome level and their expression patterns were proved to be complementary. This work lays a foundation for further research of the pathway of sRNAs and regulatory factors involved in cold tolerance.


Background
Plants have evolved sophisticated mechanisms to adapt various abiotic injuries such as heat, drought, salinity, and low temperature during their life cycle [1]. Cold injury, including chilling (<20 • C) and frost (<0 • C) is one of the key environmental factors in cold temperate zones that affect plant development and growth [2] and decide their geographical distribution. Low temperature can

Plant Materials and Treatment
Twelve overwintering branchlets of Populus simonii × P. nigra were collected from the campus of Northeast Forestry University, Harbin, China, (126 • 37 E, 45 • 42 N) and treated with hydroponic culture under a constant relative humidity of 50-60% and a 14 h light/10 h dark cycle at 24 • C ± 2 • C. Two days later, the branchlets were divided into two groups and subjected to cold (−4 • C) or not for 8 h with light (15 W/m 2 ). Then, three leaf buds were sampled from each branchlet of two groups (total 18 leaf buds from each cold treated and untreated samples) and immediately frozen in liquid nitrogen for RNA preparation. For real time PCR analysis, twelve samples were divided into two groups treated under the same conditions in triplicate and then the total RNA of the leaf buds (six buds for each replication of cold treated and untreated samples) was prepared for experiment 2.6. Total RNA extraction and RNA library construction for sequencing.

Total RNA Extraction and RNA Library Construction for Sequencing
Total RNA was isolated from the cold-treated and untreated P. simonii × P. nigra plantlets using Universal Plant Total RNA Extraction Kit (Bioteke, Beijing, China) following the manufacturer's protocol. Total RNA was quantified and assessed for quality using a NanoPhotometer (IMPLEN, Westlake Village, CA, USA) and Agilent 2100 (Agilent Technologies, Folsom, CA, USA). Cold-treated (CTD) and untreated (UD) sRNA libraries and cDNA libraries were constructed using Small RNA Sample Pre Kit (Illumina, San Diego, CA, USA) and NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations. Then, the sRNA libraries and cDNA libraries were sequenced on an Illumina Hiseq platform (Novogene, Beijing, China).

Identification and Differential Expression Analysis of miRNAs
The raw sequences were firstly processed to filter out the 5 adapter (5 -GTTCAGAGTTCT ACAGTCCGACGATC-3 ) and 3 adapter (5 -AGATCGGAAGAGCACACGTCT-3 ) sequences, N% > 10%, low-quality as well as polyA/T/G/C sequences, and to obtain the final clean reads. Then, the small RNA sequences with 18-30 nt in length were subjected to Populous tomentosa mRNAs, Rfam (http://rfam.xfam.org/), and NCBI GenBank databases (http://ftp.ncbi.nlm.nih.gov) to identify mRNA, rRNA, tRNA, snRNA, snoRNA, and repeat sequences. Finally, the remaining unique sequences were analyzed by BLAST against miRBase (Release 21, http://www.mirbase.org/) and sequences with identical or related (1-3 bases mismatch) sequences from mature miRNAs were identified as known miRNAs [24]. Potentially novel miRNAs were identified by folding the flanking genome sequence of unique small RNAs using MIREAP https://sourceforge.net/projects/mireap/. The RNA secondary structure of the predicted miRNAs was checked using Mfold http://www.bioinfo.rpi. edu/applications/mfold/cgi-bin/rna-form1.cgi [25]. The potential targets of conserved and novel miRNAs were predicted using the psRobot program [26] http://omicslab.genetics.ac.cn/psRobot/ with default parameters. The miRNA expression levels between the two libraries (CTD and UD) were analyzed by the ratio of normalized read count in CTD and normalized read count in UD. The read counts from each library were normalized to transcripts per million reads (TPM) as follows: Normalized expression = (actual miRNA count/total count of clean reads) × 1,000,000. Then DEGseq was used to differential expression analysis and the differential expression miRNAs were screened with q value (q value < 0.01) and foldchange (|log 2 (foldchange)| > 1). Hierarchical clustering of differentially expressed miRNAs was performed using the Cluster (v3.0)/Treeview software (v3.0) (http://bonsai.hgc.jp/~{}mdehoon/software/cluster/software.htm). The TPM data adjusted by log 10 (TPM + 1) was used as the value of hierarchical clustering.

Analysis of Differential Expression Genes Based on Transcriptome Sequencing
Raw data (raw reads) of fastq format were firstly processed to remove adaptor sequence, low-quality reads (more than 20% nucleotides with quality value ≤10), and shorter reads (<50 bp) and to obtain clean reads. After that, Q20, Q30, GC-content, and sequence duplication level of the clean data were calculated. Then contigs were assembled and obtained by using Trinity [27] with min_kmer_cov set to two by default. Next, the reads were mapped back to contigs to obtain longer sequences (transcripts), and the main transcript was defined as a unigene. After that, ORFs were predicted for all unigenes using Getorf (http://emboss.sourceforge.net/apps/cvs/emboss/apps/getorf.html) and Bowtie program was used to map the reads from two libraries to unigenes [28]. Differential expression analysis of CTD and UD data was performed using the DEGSeq R package (1.12.0) [29] based on the negative binomial distribution [30]. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate [31]. Genes with q value < 0.005 & |log 2 (foldchange)| > 1 found by DEGSeq were assigned as differentially expressed.

GO and KEGG Pathway Enrichment Analysis
Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R packages based Wallenius noncentral hyper-geometric distribution [32], which can adjust for gene length bias in DEGs. KEGG [33] was utilized for understanding high-level functions and utilities of the biological system (http://www.genome.jp/kegg/). We used KOBAS [34] software to test the statistical enrichment of differential expression genes in KEGG pathways.

Quantitative Real-Time PCR Analysis
Total RNA was obtained from cold-treated and untreated P. simonii × P. nigra cuttages. For cDNA synthesis, 1 µg of total RNA of each sample was reverse transcribed with stem-loop primers (miRNA expression analysis) or oligo dT primers (target gene expression analysis) using the PrimeScript Reagent Kit with gDNA Eraser (Takara, Dalian, China). Several cold-responsive miRNAs and genes were validated and quantified using qRT-PCR. Primers were designed according to the miRNA sequences (Table S1) and assembly data using the online Primer BLAST program (https://www.ncbi.nlm.nih.gov/ tools/primer-blast) ( Table S2). The qRT-PCR was performed using the LightCycler 480 SYBR GREEN I Master (Roche Applied Systems, Mannheim, Germany) and ROCHE LightCycler 480 real-time system (Roche Applied Systems). Each reaction with a 20 µL volume was carried out containing 10 µL of PCR Master Mix, 1 µL of the first-strand cDNA, and 0.5 µM of each primer. The PCR amplification program consisted of 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s and 60 • C for 1 min. A thermal denaturing cycle of 95 • C for 1 min and 65 • C for 30 s was applied to determine the dissociation curves for verifying the specificity of the PCR amplifications. Three biological replicates and three technical replicates were performed for each of the analyzed genes. The actin gene was used as an internal reference for normalization. Relative transcript levels of each gene were calculated with the comparative cycle threshold (2-dCt) method [35].

Identification of miRNAs in Populus simonii × P. nigra
Two RNA-seq libraries from cold-treated (CTD) and untreated (UD) plantlets were constructed for high-throughput sequencing. After discarding the low-quality sequences and adapter sequences, 6,356,896 (CTD) and 6,125,123 (UD) clean reads were obtained. Of 18-30 nt sRNAs, 24 nt miRNAs were the most abundant in CTD and UD samples ( Figure 1) which is consistent with the distribution patterns of sRNAs in Brassica napus and tea plant [15,36]. Also, 10,089 (CTD) and 45,241 (UD) small RNA reads were found to be similar to known miRNAs after comparing with mature miRNA sequence or miRNA precursor in miRBase 21.0. The rest of the sequences were rRNA, snRNA, snoRNA, novel miRNA, or ta-siRNA (Table 1). Among the classified sRNAs, the proportion of sRNAs reduced in CTD compared to sRNAs in UD sample implying that miRNAs might be involved in the regulation of cold tolerance. Except known miRNAs, the function of unknown cold responsive sRNAs remained to be identified. A total of 87 (UD) and 63 (CTD) conserved miRNAs belonging to 34 miRNA families were identified in the two libraries (Table S3). Among these miRNA families, the miR166 family is the most abundant, followed by miR159 and miR472 family in both libraries. Other miRNA families such as miR482, miR319, and miR398 were also detected at high copy number. According to the predicting hairpin structures, 16 novel miRNAs candidates with hairpin sequence were identified (Table S4). Additionally, we identified the candidate phased-siRNA in Populus simonii × P. nigra from CTD and UD sRNA libraries. In total, 31 phased-siRNA were obtained (reads ≥4) and mainly derived from four predicted phased loci (Table S5). The Cluster-12280.0 loci overlapped by siRNAs is homologous to TAS4 which is predicted to be triggered by miR828, and Cluster-16493.19914 loci is similar to TAS3 which is predicted to be bound by miR390 and initiated phasiRNA ( Figure 2). Among the candidate ta-siRNA, most were identified from UD sRNA libraries.

Differential Analysis of Candidate miRNAs Involved In Cold Tolerance in Populus simonii × P. nigra
To accurately measure the abundance of miRNA, the reads of each miRNA was normalized to generate RPM (reads per million) and analyzed by DEGseq [29]. The threshold of log 2 Ratio ≥1 was applied to identify the differential expression miRNAs. In total, 61 miRNAs including 57 known miRNAs and four novel miRNAs were significantly up-or downregulated under cold-treated or untreated conditions ( Figure 3). Among the sRNAs detected in both CTD and UD libraries, 11 miRNAs, including miR164a, miR168a-5p, miR160e-5p, and miR482d-3p, were upregulated with the treatment of cold, while 50 miRNAs such as miR159a, miR1447, miR167e, miR168a-3p, miR1450, and miR162a were downregulated after cold treatment. Also, the greatest changes in abundance were miR171e and miR482d-3p, which showed more than four-fold differences in expression with cold treatment compared to without cold treatment. In addition, 25 miRNAs, including miR399f, miR6427-3p, miR160b-3p, miR482b-3p, and miR7817a, were only detected in the UD library, but certain miRNAs such as miR171c, miR172h-5p, miR390d-3p, miR393a-3p, miR482b-5b, and miR6421 were transcribed at low levels (Table S6).

Identification of Candidate Cold-Responsive Genes by Transcriptome Analysis
The cold-responsive differential expression genes (DEGs) were also analyzed and identified through sequencing CTD and UD plantlet cDNA libraries. After removing low-quality reads and adapter sequences, 54,357,974 and 48,912,024 clean reads were obtained from the UD and CTD raw data separately. The Trinity program [27] was used to assemble transcripts and among 178,416 transcripts, 92,755 transcripts were annotated to be encoding genes ( Table 2). Among the databases of NR, NT, PFAM, GO, and KOG, 78,569 (84.7%) genes could be annotated in NT and 63,389 (68.34%) in NR, and only about 50% (45,817) genes could be annotated by GO ( Figure S1). In addition, most of genes had significant matches with Populus trichocarpa (60.4%) and Populus euphratica (25.4%) ( Figure S2). The FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced) was calculated and only those genes that had FPKM values > 0.3 were kept for analysis. After being normalized by TMM, a total of 1229 DEGs were detected by DEGseq analysis (q value < 0.005 & |log 2 (fold change)| > 1) including 667 upregulated and 562 downregulated genes under cold-induced conditions ( Figure 4 and Table S7). Differential expression analysis showed that during the cold treatment, many significantly differentially expressed genes are involved in the physiological and biochemical process of cold response. Further analysis by GO showed that among the 45 subsets of GO categories, the functional categories of metabolic process, carbohydrate metabolic process, oxidation-reduction process and oxidoreductase activity, metal ion binding, and cation binding are predominant ( Figure S3). The cluster of GO: 0009409 represents response to cold and had four genes  (Table S9).

Expression Analysis of miRNAs by qRT-PCR
To validate the obtained differentially expressed miRNAs in response to cold, eight miRNAs were selected for further analysis by stem-loop qRT-PCR. Our result showed that miR167a, miR1450, miR319a, miR395b, miR393a-5p, miR408-5p, and miR168a-5p were downregulated after cold treatment ( Figure 5), whereas transcription level of miR172a increased after cold treatment. The expression of most detected miRNAs was consistent with those by high throughput sequencing, except for miR172a and miR168a-5p. However, the transcription level of miR168a-5p was higher in both CTD and UD plantlets than that of other detected miRNAs, which was consistent with high throughput detection result. We also detected the expression of miR395b, miR393a-5p, miR408-5p, and miR168a-5p in CTD and UD roots and the results indicated that the expression of miR395b, miR393a-5p, and miR168a-5p decreased while the transcription level of miR408-5p was high in roots of cold-treated plants ( Figure 6).

qRT-PCR Validation of the Transcription Level of Target Genes
Through high throughput sequencing, many candidate genes have been screened from cold-treated and untreated plantlets of Populus simonii × P. nigra. To validate the results obtained by high-throughput sequencing, the transcription level of nine genes including Cluster-15451. (annotated squamosa promoter-binding-like protein 1, SBP1) and Cluster-16493.31604 decreased after cold treatment. The most abundant transcription level of these detected genes was Cluster-16493.31604 which was annotated coatomer subunit alpha-1-like isoform X1 (COPA1) and the least expression was Cluster-1378.1 which was annotated auxin response factor 2 (ARF2). In addition, transcription factor GAMYB was predicted to have the binding sites of ptc-miR319i, ptc-miR319a, ptc-miR319e, novel_63, ptc-miR159a, and ptc-miR159d. Also, the RFK1 (with binding site of ptc-miR164a), APS1 (with binding site of ptc-miR395b), SBP1 (with binding site of ptc-miR167a), bZIP44 (with binding site of ptc-miR172), LHW (with binding site of ptc-miR395b), ATL38 (with binding site of ptc-miR393a-5p), COPA1 (with binding site of ptc-miR482a.1, ptc-miR168a-5p, novel_48), ARF2 (with binding site of ptc-miR319a) were all predicted to be regulated by isolated miRNAs. Moreover, the expression patterns of ptc-miR395b with its candidate target LHW, and ptc-miR393a-5p with its candidate target ATL4, are opposite in CTD and UD plantlets (Figures 5-7).

Discussion
In plants, cold injury severely affects geographic distribution. Recent reports have suggested that miRNAs play important roles in various abiotic response during development of plants through negatively controlling the expression of their target gene [2,37,38]. Although many genes have been identified in Populus under cold treatment [2,13], few studies have reported on Populus simonii × P. nigra, which is widely distributed in Northeast China. Further study is needed to elucidate the functions of miRNAs at a genome-wide level in response to cold tolerance. There is always a cold spell at the end of April and early May in which the temperature can drop below zero degree with snow in Heilongjiang area of China. Through cold treatment, we found plants treated with −4 • C for 12 h show obvious damage by freezing and can recover normal growth, but plants treated with −8 • C for 12 h cannot grow normally and plants treated with −4 • C for 4 h show light damage ( Figure S4). Then, we used high throughput sequencing to analyze the transcription level of miRNAs and transcriptome in Populus simonii × P. nigra with or without treatment with −4 • C for 8 h. The results showed that the expression of most conserved miRNAs decreased and the transcription level of their candidate targets significantly increased in response to cold treatment, suggesting that cold induction lead to increased cold-responsive gene expression through decreasing the level of miRNAs. Thus, genome-wide analysis of miRNAs and targets response to cold can provide crucial information for our understanding of the molecular regulation of cold tolerance in Populus simonii × P. nigra.

Cold-Responsive miRNAs in Populus simonii × P. nigra
In our study, the transcription levels of miRNAs were analyzed in Populus simonii × P. nigra with and without cold treatment and the results showed that the expression of about 70% of miRNAs was affected by cold treatment. Among them, about 80% of miRNAs were downregulated, indicating that decreased expression of miRNAs might play more crucial roles than increased expressions of miRNAs in Populus simonii × P. nigra. In Populus tomentosa, downregulated miRNAs were more common than upregulated miRNAs in response to cold treatment [2]. Previous research indicated that psu-miR475 in Populous suaveolens [39], miR160, miR168, miR390, and miR396 in Populus trichocarpa [40], and pto-miR171, pto-miR319, and pto-miR395 in Populus tomentosa [2] respond to cold. Similarly, in the present study, these miRNAs showed obvious changes under cold treatment. Additionally, conserved miRNAs, such as ptc-miR1450, ptc-miR159a, ptc-miR164a, ptc-miR167a, ptc-miR172a, and candidate novel ptc-novel_48, ptc-novel_50, ptc-novel_59, and ptc-novel_76, were differentially expressed under cold treatment. These obtained miRNAs might downregulate the expression of their target genes which encode regulatory and functional proteins involved in cold tolerance.
Phased, secondary small interfering RNA (phasiRNA) are derived from phasiRNA-producing loci which were originally discovered in Arabidopsis and were also called trans-acting small interfering RNA (ta-siRNA) [41]. Four known TAS genes (TAS1 and TAS2 are targeted by miR173, and TAS3, TAS4 are targeted by miR390 and miR828 separately) are conserved in plants and their functions relate to each target of ta-siRNA and the expression pattern of targets [41,42]. Our results showed the siRNA Nt_403935_×34, with a similar sequence to 5' D8 (+) originated from Cluster-16493.19914, was predicted to be ptc-TAS3 loci and was bound by ptc-miR390. The targets of siRNA (Nt_403935_×34) were predicted to be auxin response factors and their transcription level increased in CTD plantlets. Also, miR390-TAS3 tasiRNA-ARF2/3/4 has been reported to be integrated with auxin signaling to regulate lateral root growth of Arabidopsis [43]. This implied that miR390-TAS3 tasiRNA-ARF2/3/4 might also be involved in the process of cold resistance of Populus simonii × P. nigra. TAS4 loci, where the siRNAs were triggered by miR828 and targeted a set of MYB transcription factors [44] was also be predicted, but ptc-miR828 was not detected in either CTD or UD plantlets of Populus simonii × P. nigra. We postulated that the reason for undetected ptc-miR828 might be related to the low transcription level of ptc-miR828 verified by stem-loop RT-PCR according to the sequence of pto-miR828a in Populus tomentosa.
The predicted targets of differentially expressed miRNAs identified in our research have also been reported to be involved in cold tolerance. Our results showed the target gene of cold-responsive ptc-miR164a was LRR (leucine-rich repeat) receptor kinase, which is reported to be cold-inducible and has Ser/Thr protein kinase activity [45]. The conserved miR160 with potential target of ARF (auxin response factor) also showed a cold-resistance response in wheat [20]. The cold responsive miR167 has been identified to have LTR (low temperature response) element in the upstream regions [46] and candidate low temperature-induced protein was also predicted in our research. Other targets of miR167, such as ARF and SPL12 (Squamosa promoter-binding-like protein 12) were also predicted and reported to response to cold treatment [17] or thermo tolerance [47]. DnaJ, the candidate target of miR162a was also reported to contribute to maintenance of photosystem II under chilling induction [48]. In addition, Myb and TCP, candidate targets of ptc-miR319, have also been reported to be subjected to cold resistance in sugarcane [49]. Our results showed that cold induced the changes in transcription level of ptc-miRNAs and therefore led to the expression of target genes involved in cold resistance. Certainly, many of the cold responsive miRNAs have also been reported to be involved in plant growth and development, such as miR159 in growth and programmed cell death [50], miR172 in floret development [51], and miR319 in leaf development [52]. These results suggest that miRNAs have complex and widely regulatory roles in various plant stress and growth responses.

Differentially Expressed Transcripts under Cold Treatment
The transcriptome analysis showed many candidate genes differentially expressed in CTD and UD plantlets of Populus simonii × P. nigra which might be cold responsive genes. The functions of orthologs for these candidate genes, especially those encoded transcription factors involved in cold tolerance, have been clarified in other plants. DREB, the orthologs of which have been identified from Populus simonii × P. nigra in our research, is well known to be involved in abiotic stress. For example, AaDREB1, isolated from cold tolerant plant Adonis amurensis, showed enhanced tolerance to low temperature in transgeneic Arabidopsis and rice [53]. DREB/CBF (C-repeat binding factor) can recognize and bind to the CRT/DRE (C-repeat/Dehydration Responsive Element) motif to regulate cold-responsive (COR) genes [54]. The orthologs of leucine zipper protein (Cluster-16493.40514, Cluster-16493.27643) were also identified in Populus simonii × P. nigra and leucine zipper protein was reported to interact with CBF1B and involved in low temperature response in hot pepper [55]. MYB-like transcription factor (Cluster-16493.30912, Cluster-16493.18975) belongs to the R2R3-MYB family and recently MdMYB23 has been reported to be involved in cold tolerance in apple [56]. In rice, WRKY71 had positive function in cold tolerance by regulating downstream target genes [57] and Cluster-16493.45530 and Cluster-16493.18576 were also annotated to be WRKY family genes. MfERF1, one of ethylene-responsive transcription factor (ERF) was isolated from Medicago falcata and conferred cold tolerance [58] and the orthologs of ethylene-responsive transcription factor (Cluster-16493.20988, Cluster-16493.19581) were also isolated in our results. Moreover, the MADS genes [59], bHLH [60], bZIP [61], were also reported to be involved in cold resistance.

miRNAs May Be Involved in Cold Tolerance by Negatively Regulating Target mRNAs
The downregulated miRNAs in response to cold treatment in Populus simonii × P. nigra would positively regulate the expression of their targets for cold resistance. Integrated analysis of miRNAs and their targets expression helps to reveal the regulatory pathways of functional miRNA-mRNA modules ( Figure 8). Based on data analysis, the putative cold responsive regulation model showed sophisticated up-and downregulated network. Firstly, the cold signal was perceived and transmitted to transcription factor genes such as MYB, SBP, bZIP, ARF, LHW, ATL, MIR genes such as miR319, miR159, miR167, miR172, miR395, miR393, miR390, novel_63, and ncRNA genes which can be targeted by miRNA for generating siRNAs. The expression of putative MIR showed downregulation, whereas the expression of candidate transcription factor genes were upregulated. Secondly, the negative regulation between miRNAs and their candidate targets, such as miR390 and noncoding RNA TAS3, miR319, and MYB led to the expression of cold responsive transcription factors and accumulation of cold responsive metabolic pathway genes. Thirdly, the metabolic genes related to the LRR receptor kinase pathway, ARF pathway, SPL pathway, and DnaJ-mediated photosystem II are taken as candidate genes involved in cold tolerance of Populus simonii × P. nigra. Although cold responsive miRNAs and candidate genes were isolated and characterized in Populus simonii × P. nigra and several other species, such as wheat [20], sugarcane [62], tea plant [15], soybean [63], and tomato plant [64], the mechanism of molecular regulation at the transcriptional and post-transcriptional level involved in the network of cold tolerance is still largely unknown. Further work is needed to elucidate the function of these important genes in cold resistance-related regulatory pathways.

Conclusions
Cold tolerance miRNAs and candidate target genes were identified through integrated sRNA and transcriptome analysis in cold treatment of Populus simonii × P. nigra. MIR genes such as miR319, miR159, miR167, miR172, miR395, miR393, miR390, and novel_63 and transcriptional factors including MYB, SBP, bZIP, ARF, LHW, and ATL showed differential expression and they might be the main contributors related to LRR receptor kinase, the ARF pathway, the SPL pathway, and DnaJ-related photosystem II involved in the cold tolerance of Populus simonii × P. nigra. These results not only increase our knowledge of sRNAs involved in the post-transcriptional regulation of cold tolerance, but also provide candidate genes for future functional analysis of the cold tolerance-related signaling pathways in Populus simonii × P. nigra.