Transposon-associated small RNAs involved in plant defense in poplar


 Background: Small RNAs (sRNAs) are hypothesized to contribute to plant defense responses by increasing the overall genetic diversity and regulating gene expression; however, their origins and functional importance in plant defense remain unclear. Here, we use Illumina sequencing to assess how sRNA populations vary in the Chinese white poplar (Populus tomentosa) during a rust fungus (Melampsora larici-populina) infection. We sampled sRNAs from the biotrophic growth phase (T02; 48 h post infection) and the urediniospore formation and release phase (T03; 168 h), two essential stages associated with plant colonization and biotrophic growth in rust fungi. Results: The proportion of siRNA clusters located in pseudogenes and transposons was significantly larger than would be expected by chance and infection-stage-specific differences in siRNAs primarily originated from those in the transposon regions. We also found that the abundance of clusters comprising 24-nt siRNAs located in the transposon and intergenic regions underwent more substantial changes as the infection progressed. A target analysis revealed that 95% of fungal genes were predicted to be targets of Populus sRNAs. Pathogen effector genes were targeted by more sRNAs identified during the biotrophic growth and urediniospores formation and release phases than in the control plants, suggesting a clear selection for sRNA-target interactions. Compared with the miRNAs conserved between different plant species, a significantly higher proportion of Populus-specific miRNAs appeared to target NB-LRR genes. Conclusions: This integrated study on the plant colonization and biotrophic growth in rust fungi profiles could provide evolutionary insights into the origin and potential roles of the sRNAs in plant defense, coevolution with pathogens, and functional innovation.


Background
Small RNAs (sRNAs) are ~21-24-nucleotide (nt) RNAs that can regulate gene expression, maintain genome integrity and chromatin structure, and influence plant development, such as seed germination [1][2][3][4]. sRNAs typically have an inhibitory effect on gene expression, and the corresponding regulatory mechanisms are therefore subsumed under the heading of RNA silencing [5].
Two major classes of endogenous sRNAs have been identified in plants, namely microRNAs (miRNAs) and small interfering RNAs (siRNAs) [5,6]. Plant miRNAs are transcribed from imperfect hairpin-shaped precursors and are typically 21 nt long, while siRNAs are 21-24 nt in length and are generated from double-stranded RNA duplexes [7]. miRNAs reduce the protein abundances of their targets through the direct cleavage or posttranscriptional regulation of target mRNAs to reduce their translation [8], which has important functions in a wide range of biological processes [5]. In Arabidopsis, over 100 miRNAs have been shown to be important for plant development [9] and abiotic stress tolerance [10].
siRNAs have been sequenced in many plant species, such as Arabidopsis thaliana, rice (Oryza sativa), and wheat (Triticum aestivum) [11][12][13][14], and have the ability to silence their targets by directly complementing their transcripts to prevent their translation. Studies in maize (Zea mays) have revealed that sRNA populations vary between parents and their progeny, contributing to the dramatic vigor of maize hybrids [15].
Innate immunity is an evolutionarily ancient mechanism that protects plants from a wide range of pathogens [16,17]. Many lines of evidence have confirmed that sRNAs contribute to plant defenses against pathogens [17][18][19]; for example, changes in the accumulation of specific miRNAs, such as miR393, miR160, and miR773, have been reported during the infection of Arabidopsis by the bacterial pathogen Pseudomonas syringae [20,21] [28,29]. It is almost certain that a perennial plant will encounter a pathogen before it can reproduce, and the long generation intervals of trees makes it difficult for them to match the evolutionary rates of microbes [30]. These analyses reveal the evolutionary importance and functional novelty of the sRNAs.

Plant material and fungal treatments
The Chinese white poplar (Populus tomentosa) was planted in pots under natural light conditions (12 h of 1,250 μmol m -2 s -1 photosynthetically active radiation) at 25±1°C (day and night) and 50±1% relative humidity (day and night) in an air-conditioned glasshouse using soilless culture technology. The infection experiments were performed using M. larici-populina strain G5. Six inoculation spots per poplar leaf were each inoculated with 5 μl of 10 6 condia ml -1 , with three biological replicates per treatment. The inoculated leaves were incubated in an artificial climate incubator (LT36VL; Percival Scientific, Inc; Perry, U.S.A.) under 95% relative humidity and 25±1°C, and were harvested at 48 (T02) and 168 (T03) hours post inoculation (hpi).
In the control group (T01), the leaves were sprayed with sterile tap water and harvested at 48 hpi. The samples from all leaves exposed to the same treatment were pooled and treated as one biological repeat, and three independent experimental repeats were performed for each treatment (control, 48 hpi, 168 hpi). All samples were immediately frozen in liquid nitrogen and stored at -80°C for RNA extraction.

miRNA library construction
Total RNA was extracted with Trizol (Thermo Fisher Scientific) according to the manufacturer's protocol [35]. The total RNAs generated by Dicer processing are efficiently targeted by the included modified adapters. RNA-seq libraries and sRNA libraries were prepared using the TruSeq RNA Library Preparation Kit and TruSeq Small RNA Library Preparation Kit (Illumina), respectively [36]. sRNAs with adapters on both ends were used as templates to create cDNA constructs using reverse transcription PCR. After being purified and quantified using a Qubit dsDNA HS (Qubit 2.0 Fluorometer) and Agilent 2100, the PCR products were used for cluster generation and sequencing on an Illumina HiSeq 2500, according to the cBot and Hiseq 2500 user guides, respectively. sRNA sequencing was performed with two biological replicates for each treatment and transcriptome sequencing was performed with three biological replicates for each treatment.

mRNA sequencing, alignment, and normalization
Quality control of raw reads was carried out using FastQC to remove reads containing >50% low-quality nucleotides (sQ <= 5) and >10% ambiguous nucleotides.
The qualified reads that passed the filters above were then aligned using TopHat2 [37] with the parameters --read-mismatches 2 -p -G to generate read alignments for each sample. Up to two mismatches were permitted in each read alignment. The transcript abundance was calculated, and differential transcript expression was computed using CuffDiff with the parameters -p -b [38]. Run bias detection and correction algorithms to improve the accuracy of transcript abundance calculations.

Effector prediction
Effector prediction was performed using a machine learning program, EffectorP, with default parameters [39]. In this analysis, all of the protein sequences in the M.

sRNAome analysis
The sequences generated from the leaves exposed to the three infection treatments were used to detect the transcript abundance of mature sRNAs. Figure S1 shows the procedures used for sRNA library construction. Over 16 million raw reads were produced per library, of which 33.5-44.0% aligned perfectly (no more than one mismatch) to version 3.0 of the P. trichocarpa genome. These perfectly aligned sequences were annotated by BLAST-searching them against the GenBank and Rfam databases (version 13; http://rfam.xfam.org/), allowing one mismatch. The tRNAs, rRNAs, snRNAs, snoRNAs, and scRNAs were removed from the sequencing reads.
The remaining unannotated sRNAs were searched against the known miRNAs from miRBase version 22.0 (http://www.mirbase.org/), allowing a maximum of two mismatches. To predict the novel miRNAs, miRDeep2 [40] was used.

Detection of sRNA clusters
To identify the siRNA clusters, the strategy used by Johnson et al. (2009) [41] was applied, with some modifications. The 21-24-nt sRNAs from all the samples that did not match the miRNAs, ribosomal rRNA, or tRNAs and that could map to the P. trichocarpa genome were processed together. siRNAs within 100 bp of each other were merged into blocks referred to as siRNA clusters. The coordinates of the clusters were defined by the first and last siRNAs of the overlapping sequences. siRNA clusters with a transcription abundance of at least 5 tpm were included in the downstream analysis. The clusters were labeled according to the most abundant length of siRNA they contained (21, 22, or 24 nt), and the type of genomic region in which they were located.

Genomic locations of miRNAs in the P. trichocarpa genome
The siRNA clusters were screened for their localization within the transposons, introns, pseudogenes, intergenic regions, 5′ UTR, 3′ UTR, and upstream and downstream (2 kb) of coding genes of the P. trichocarpa genome (Phytozome version 12), with an overlapping rate of above 80%. To assess the overrepresentation of siRNA clusters in each genomic element, deviations from expected proportions were assessed using randomization tests. The genomic locations containing the same length of siRNA clusters were randomly sampled 1,000 times, after which the overlap of the siRNA clusters with these known genomic elements was assessed. A one-tailed P-value was estimated based on the Z score of the difference between the actual count minus the mean count of 1,000 surrogates relative to the standard deviation from the 1,000 surrogates. The P-value was the probability to the right of the observed count, calculated using a one-tailed z test.

TE annotation
The TE annotations used in this study were obtained from the outputs of the These RM outputs were filtered to remove non-TE elements, such as satellites, simple repeats, low complexity sequences, and rRNA.

Pseudogene annotation
The intergenic sequences of the P. trichocarpa genome were used to identify the putative Pseudogenes. The overall pipeline used for this identification was generally based on the PlantPseudo workflow [43], and consisted of four major steps: (1) identify the masked intergenic regions with sequence similarity to known proteins using BLAST; (2) eliminate redundant and overlapping BLAST hits in places where a given chromosomal segment has multiple hits; (3) link homologous segments into contigs; and (4) realign sequences using tfasty to identify features that disrupt contiguous protein sequences.

Poplar sRNAs expressed during a rust fungus infection
To study the expression profile of the Chinese white poplar sRNAs during the biotrophic growth (T02) and urediniospore formation and release (T03) stages of infection by Melampsora larici-populina [44], sRNA libraries were prepared and subjected to high-throughput sequencing. Specifically, we prepared a library of RNAs approximately 18-30 nt in length from the three samples and sequenced them using an Illumina Hiseq 2500, generating 81.3 million clean reads in total (Table S1).
After processing the sRNA data and combining the sequences across the libraries, we identified a set of distinct sRNAs for the two infection stages and the mock-inoculated control plants. The resulting data were evaluated to ensure they were consistent with the expected size range of functional siRNAs and miRNAs, with 56.6-65.4% of the generated sRNAs having an expected sequence length of 21-24 nt ( Figure S2). In total, 33.4-44.0% of the total reads could be aligned perfectly (with no more than one mismatch) to the Populus trichocarpa genome (version 3.0)[30].

Pseudogenes and transposons act as catalysts for the formation of siRNA clusters
We next studied the locations of the siRNAs to elucidate how they might have evolved. The sRNA reads matching miRNA, tRNA, or ribosomal DNA sequences were removed from the datasets, and only those that mapped to the genome and were of the characteristic siRNA length (21-24 nt) were retained. The vast majority of siRNAs expressed during the rust fungus infection were found to be unique to one stage, while few siRNAs were common to all three samples ( Figure 1A). Thus, infections could therefore create stage-specific siRNAs, producing a more complex and targeted siRNA population.
We used an approach similar to that reported by Barber et al. (2012) [15] to identify the types of genetic features associated with the stage-specific differences in siRNAs. We grouped overlapping 21-24-nt siRNAs within a ≤ 100-bp window on the P. trichocarpa genome into clusters based on their location within the genome.
Totally, we identified 301,155 siRNA clusters. The relative abundances of the siRNAs clusters were calculated (tpm: transcripts per million mapped reads). We then classified and annotated the siRNA clusters according to their common siRNA length (21, 22, or 24 nt) and their location in the genome. The 22-nt clusters accounted for less than 10% of all the clusters, and the 24-nt clusters accounted the largest proportion (> 50%). The siRNA clusters with an abundance of at least 5 tpm in at least one stage were retained for the downstream analyses (Data S1). To investigate the clusters exhibiting large stage-specific differences, we ordered the clusters by their absolute value of fold change values (T02/T01 or T03/T01).The degree of difference in the 24-nt siRNA cluster abundances between stages increased for those located in the transposon and intergenic regions ( Figure 1C and 1D), with the strongest trend observed in the transposon regions.
To further investigate the global differences in siRNA abundance in the retrotransposon families between the different infection stages, we mapped the 21-nt, 22-nt, and 24-nt siRNAs that perfectly matched the Populus genomes onto the retrotransposons. Figure 2A shows the retrotransposon families that showed a siRNA abundance of at least 100 tpm in at least one of the stages. The production of 21-nt siRNA lengths and the abundance of these siRNAs located in the retrotransposon families are more similar across the three stages than those of the 22-nt and 24-nt siRNAs.
The stage-specific differences in the abundance of siRNAs located in the transposons were examined using a X 2 test for each siRNA length. The differences in the siRNA abundance for most transposon families between the stages (T02 vs T01; T03 vs T01) was contingent on the 24-nt length sequences (P < 0.05; Chi-square test; Data S3; Figure 2A). The stage-specific differences in siRNA abundance for these retrotransposon families thus primarily resulted from the 24-nt siRNAs.

Populus siRNAs mediate plant-pathogen interactions, as revealed using a shotgun strategy
To explore the potential role of the Populus siRNAs in plant-pathogen interactions, we profiled the siRNA sequences that mapped to the P. trichocarpa genome which could regulate the genes of the fungal pathogen Melampsora larici-populina. In the biotrophic growth (T02) and urediniospore formation and release (T03) stages, more than 95.0% of the fungal genes (T02: 15,954; T03: 15,472; total fungal genes: 16,372) were predicted to be targeted by siRNA sequences (Data S4). Nearly one-third of the expression of Melampsora larici-populina genes detected in our transcriptome were altered at least by twofold between sample pairs (T02/T01; T03/T01; T02/T01) (Data S5). We hypothesize that some siRNAs could suppress pathogen pathogenicity by targeting pathogen virulence genes using a shotgun strategy. Importantly, if the siRNA-pathogene targets not in a random match manner, we would expect to see the targets of sRNAs may bias some pathogen genes. As expected, compared with the control sample (T01), pathogen effectors were targeted by more siRNAs at the biotrophic growth and urediniospore formation and release phases (P < 0.001; permutation test; Figure 2B), suggesting there is a clear selection for the siRNA-target interactions.
The Populus-specific miRNAs are more involved in the regulation of the

disease-resistance (DR) genes
We next analyzed the presence of miRNAs in the sequencing data, as miRNAs have been shown to play critical roles in plant defense systems [45,46]. We annotated the sRNAome by aligning the sRNA reads to Rfam (version 13). After removing the scRNAs, tRNAs, rRNAs, snRNAs, and snoRNAs, we identified 423 sRNAs precursors, which grouped into 211 miRNA families containing an average of 2.0 genes per family for each of the three stages. We found that 67.7% of the miRNAs could be annotated in miRBase (version 22); however, most of the putative novel miRNAs did not remain in the final dataset. These excluded miRNA sequences may not represent bona fide miRNAs or may be expressed at levels too low to reliably show that they have the characteristics of miRNAs (at least one read matching their star sequence; Data S6). The expression values of each miRNA gene were also calculated and normalized (Data S7).
Using the miRNAs identified in this study and the miRNAs of the closely related species annotated in miRBase (version 22), we first sought to identify the Populus-specific miRNAs. Of the 211 miRNA families, 181 were Populus-specific, and 70 were supported by the findings of a previous study ( Figure 3A) [47]. We then predicted the potential miRNA targets using psRNATarget [42]. A total of 65,777 predicted miRNA-target gene pair matches with an expectation score of 5 were identified (Data S8). An examination of the expression of the miRNA/target genes revealed 1,940 pairs with negatively correlated expression patterns (P < 0.05; Pearson correlation; Data S9). Using the publicly available degradome library [48], we identified 553 miRNA/target pairs representing 165 miRNAs (Data S10), of which several conserved miRNA/target pairs (miR156-SPL, miR172-AP2, miR397-LAC) were observed in our dataset [49][50][51]. When the predicted DR gene targets of miRNAs in both the conserved and Populus-specific datasets were excluded, a total of 161 DR targets remained for the conserved miRNAs and 711 DR targets for Populus-specific miRNAs. More than 66% (120 of 181) of the Populus-specific miRNAs target DR genes, a significantly larger proportion than that of the conserved miRNAs (23.3%, P = 0.009; Fisher's exact test; Table 2, Data S11; Figure 3B).
To gain a better understanding of the functional roles of the Populus-specific miRNA families, we performed a pfam domain analysis of their predicted targets. The top 20 domains of the Populus-specific miRNA targets were highly correlated with the domains of the plant DR genes ( Figure S1); for example, NB-ARC and LRR domains, which are known to dissociate upon effector-dependent activation [52], were overrepresented in the samples.

Populus-specific miR6579 and miR6590 target numerous NB-LRR genes
Plant NB-LRRs (encoded by NB-LRR genes) are commonly targeted by miRNAs, typically generating phasiRNAs (phased small-interfering RNAs), which can reduce the levels of these targets in cis or in trans [48]. In the two stages of the rust fungus infection, a conserved miRNA superfamily composed of miR482 targeted the NB-LRRs at the encoded and conserved P-loop motifs. This family was found to be highly conserved among divergent plant species, including rice, Arabidopsis, grape (Vitis vinifera), and soybean (Glycine max) [23, 53,54]. A second conserved family, miR2111, was identified in the eudicotyledons (Figure 3C), indicating the existence of this miRNA family prior to the divergence of these species from a common ancestor around 103 Mya [55]. The highly conserved 20-nt miR482s were predicted to target 221 NB-LRRs. Two related Populus-specific families, miR6579 and miR6590, which both contain 21-nt miRNAs, were predicted to target 132 and 49 NB-LRRs, respectively ( Figure 3D-F). These two Populus-specific miRNA families had the highest transcript abundance during the mock-inoculated control plants, decreasing to a low level during the biotrophic growth (T02) and urediniospore formation and release (T03) stages (Data S7). By contrast, the highly conserved 22-nt miR2111s were predicted to target only nine NB-LRRs, three of which were also predicted targets of miR482 while two were predicted to be targeted by miR6579 ( Figure 3D-F).
These results indicate that the NB-LRRs in Populus are primarily targeted by conserved miR482s and Populus-specific miR6579s.
A careful examination of their precursor sequences revealed that miR6579 and miR6590 were intron-derived and pseudogene-derived miRNAs, respectively ( Figure   3E). The parent of the pseudogene Chr03|4847214-4847421 was a NB-LRR gene (Potri.001G426600), suggesting that pseudogenes could act as catalysts for the formation of miRNAs. In poplar, the 21-nt miR6590 family showed characteristics canonical of the typical miRNA triggers associated with ARGONAUTE1 (AGO1) for the initiation of the PHAS loci, such as "U" in the 5′ position, "A" in position 10, and "C" in the 3′ position, indicating that this miRNA may trigger the biogenesis of secondary siRNAs [56]. It seems that species-specific miRNAs could become fixed in the regulatory network once a number of targets have evolved for which the presence of a miRNA-binding site would be beneficial.

Transposons and posttranscriptional regulation
Despite previously being referred to as junk DNA, transposons are now known to be essential elements of most eukaryotic genomes, making important contributions to their structure, diversity, capacity, and adaptation [57]. The widely distributed transposable components are the most rapidly evolving part of the genome because they are more susceptible to sequence change than coding sequences [58]. Here, while exploring the distribution of poplar siRNAs in different parts of genome regions, including the 5′ UTRs, 3′ UTRs, introns, upstream (2 kb) and downstream (2 kb) of the coding genes, transposons, pseudogenes, and intergenic regions, we determined that transposons contributed a large proportion of the siRNAs (48%). In particular, the siRNA clusters were enriched in genomic locations comprising high-copy retrotransposon families. We identified significant variation in the 24-nt siRNAs derived from specific retrotransposon families between the different infection stages.
Based on their length and sequence similarity to the retrotransposon families from which they are derived, these siRNAs may act posttranscriptionally to mediate repressive histone modification, thus contributing to the reinforcement of the silenced state of TE chromatin [59,60].
Our findings further illustrate that it is difficult to assess the connection between siRNAs and transposons in the process of plant defense. In general, the mobilization of transposons can be mutagenic [61], so host genomes have evolved elaborate mechanisms to suppress their activities. Most of the 24-nt siRNAs are likely to be involved in the transcriptional regulation of transposons through RNA-directed DNA methylation [62], but many TE-associated siRNAs were also found to influence non-TE transcripts exerting certain regulatory roles [63]; for example, 24-nt siRNAs were found to target pathogen genes [61]. Transposons are rapidly activated by biotic and abiotic stresses however [64][65][66], potentially complicating the regulatory loop between the siRNAs and TEs.
siRNAs are an ideal tool for use in plant defense because of their target specificity and highly efficient responses [67]. Ohtsu et al. (2007) [68] proposed that the retrotransposon-associated siRNAs may target genes with homologous sequences.
Also, direct evidence in Arabidopsis suggests that 21-nt siRNAs derived from an Athila retrotransposon could posttranscriptionally regulate the UBP1b (Upstream binding protein 1b) gene to mediate the stress response [69]. In rice, transposon-derived 24-nt siRNAs could target the genes involved in hormone homeostasis, suggesting that transposon-associated regulatory modules could play a variety of roles in plant development and defense [70]. It is therefore possible that the 21-24-nt siRNAs may have another purpose besides the maintenance of genome stability, and may enhance plant defense by targeting the pathogen genes. Considering the high genetic diversity and mobilization of transposons, the regulatory system mediated by the TE-derived siRNAs could play a major role in the performance of plant defense.

Cross-kingdom RNAi and sRNAs mediate plant-pathogen interactions
An endless arms race drives host-pathogen coevolution [71]. In plants, pathogens can invoke multiple layers of immune responses, and many defense genes were found to play key roles in this process [72,73]. Most potential pathogens are blocked by nonhost resistance and pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) [67]; however, many pathogens deliver a variety of effectors into plant cells to suppress host immunity. To counteract this attack, plants have evolved resistance proteins and immune mechanisms to sense these effectors and trigger a robust resistance response [74]. Plant sRNAs are essential regulatory noncoding RNAs that can posttranscriptionally silence their target genes [75,76]. Recent progress in profiling sRNAs, especially advances in next-generation sequencing technology, has revealed their extensive and complicated involvement in interactions between plants and viruses, bacteria, and fungi [75].
Fungal pathogens deliver sRNAs into plant cells to induce cross-kingdom RNAi of the plant immunity genes [77]. A similar strategy was also observed in animal-parasite systems [78]. Furthermore, sRNAs generated from the host plant cells can also be transferred into pathogen cells [71], as fungal pathogens are known to be capable of taking up external sRNAs and long dsRNAs [71,79]. Pathogen genomes encode variable effector repertoires for suppressing host immunity [80]; specifically, these effectors are often encoded in genomic regions of high plasticity [81], such as TEs [82]. Given that effectors are often clustered in these regions, TEs must function as an evolutionary force shaping these virulence factors by contributing to their diversification [82]. Sequence diversification is particularly important for pathogen effectors, which are essential components facilitating their coevolution with, and outmaneuvering of, the host.
Considering the high diversification of virulence factors generated by pathogens, it is reasonable to assume that transposon-derived sRNAs are also the best sources of new regulators for plant defense resistance, facilitating cross-kingdom RNAi. TEs, the most abundant component of most eukaryote genomes, are hotspots for generating sRNAs, which is particularly true of the retrotransposons [83]. Here, we found that pathogen effectors were targeted by more poplar siRNAs in the biotrophic growth and urediniospore formation and release stages of infection by a rust fungus, in comparison with control. It is therefore reasonable to conclude that the activation of the siRNA repertoire does not occur randomly, and that their targets are instead involved in more central roles in the plant-pathogen interactions. Supporting this hypothesis, a previous study reported that the transfer of host sRNAs could silence virulence-related genes and suppress fungal pathogenicity in Arabidopsis [84].

Disease resistance associated with Populus-specific miRNAs
The plant NB-LRR genes mediate effector-triggered immunity by acting as key receptors during the innate immunity response against a wide range of pests and diseases [85]. NB-LRR genes are generally grouped into two subclasses: the Toll/Interleukin-1 receptor-like group (TIR-NB-LRRs) and a coiled-coil domain-containing group (CC-NB-LRRs) [86]. Both classes can be triggered by miRNAs to generate phasiRNAs, which can reduce the levels of the transcripts of their targets in cis and in trans [87]. In contrast to low-copy genes, many NB-LRR genes have undergone dramatic duplications and losses, domain architecture variations, the partitioning of subfamilies, and copy number variation among species [85,88,89]. Thus, the NB-LRR genes are highly variable, lineage-specific, and associated with the plant immune response, providing material to allow rapid adaptive evolution.
The target sites of conserved miRNAs are often located within the highly conserved domains of the target genes [90]. Some of these conserved domains, such as the P-loop motifs, are targeted by several miRNA families [56]. Unlike conserved miRNAs, newly evolved miRNAs tend not to target these conserved functional domains [90], and may instead target mRNAs simply by chance [91]. The recent salicoid whole genome duplication (WGD) event, which was shared among the Salicaceae, resulting in overrepresentation of genes associated with special functional categories [92]. Indeed, the target genes of the newly emerged miRNAs in Populus were found to have various functions, including a large number of DR NB-LRR genes, ABC transporters, and a few transcription factors, consistent with the retained salicoid duplicates observed in this genus [30]. Our analysis demonstrated that more Populus-specific miRNAs than conserved miRNAs target the NB-LRR genes, although the species-specific miRNAs have fewer targets than the conserved miRNAs.
Considering that this de novo diversity may be associated with plant defense, these new miRNAs have a greater potential to target newly evolved plant DR genes, further contributing to the phenotypic innovation of the host. Once the newly emerged miRNAs become fixed in the regulatory modules, they may gradually evolve to target more genes linked to their specific function [48,91]. Our study therefore provides insights into the origination of siRNAs and the regulation of siRNAs in the defense process.

Conclusions
In this study, we sequenced the sRNAome of Chinese white poplar during the two essential stages associated with plant colonization and biotrophic growth in rust fungi (Melampsora larici-populina). We identified a set of distinct sRNAs for the two infection stages and the mock-inoculated control plants. Compared with the control plants, more poplar siRNAs in the biotrophic growth and urediniospore formation and release stages targeted pathogen effector genes. Furthermore, there are more Populus-specific miRNAs than conserved miRNAs targeting the NB-LRR genes, although the species-specific miRNAs have fewer targets than the conserved miRNAs.
The results reveal the evolutionary importance and functional novelty of the sRNAs, would provide valuable information on comprehensive understanding of the origin, evolution, and plant defense effects of the poplar sRNAs.