Sexual Dimorphism Floral MicroRNA Profiling and Target Gene Expression in Andromonoecious Poplar (Populus tomentosa)

Although the molecular basis of poplar sex-specific flower development remains largely unknown, increasing evidence indicates an essential role for microRNAs (miRNAs). The specific miRNA types and precise miRNA expression patterns in dioecious plant flower development remain unclear. Here, we used andromonoecious poplar, an exceptional model system, to eliminate the confounding effects of genetic background of dioecious plants. This system, combined with high-throughput sequencing and computational analysis, allowed us to characterize sex-specific miRNAomes from female and male flowers. Comparative miRNAome analysis combined with quantitative real-time PCR revealed the expression patterns of 27 miRNAs in poplar flower and showed that the targets of these miRNAs are involved in flower organogenesis, Ca2+ transport, phytohormone synthesis and metabolism, and DNA methylation. This paper describes a complex regulatory network consisting of these miRNAs expressed in sex-specific flower development in a dioecious plant. The conserved and novel miRNA locations were annotated in the Populus trichocarpa genome. Among these, miRNA Pto-F70 and 4 targets are located in the sex-determination regions of chromosome XIX. Furthermore, two novel miRNAs, Pto-F47 and Pto-F68, were shown for the first time to be regulatory factors in phytohormone interactions. To our knowledge, this report is the first systematic investigation of sex-specific flower-related miRNAs and their targets in poplar, and it deepens our understanding of the important regulatory functions of miRNAs in female and male flower development in this dioecious plant.


Introduction
The genus Populus encompasses approximately 100 species divided into 5 sections and many of these species are important components of terrestrial ecosystems or important cultivated trees for wood, pulp, paper, and possibly biofuels uses. As dioecious trees, the genus Populus provides opportunities to study perennial woody plant biology and plant sexual dimorphism [1]. Moreover, the high-quality, annotated genome sequence of P. trichocarpa provides a key genomics tool, enabling genome-wide regulatory analyses [1]. Sexual dimorphism is a widely studied phenomenon in dioecious plants and likely results from different modes of selection operating in males and females [2]. As dioecious trees, poplar floral development is sexual dimorphism. The poplar flower is composed of two whorls, including a reduced perianth cup surrounding either the stamens or carpels, which begin to develop in spring (early May in Shandong province) of the year before the spring in which they will open [3,4]. Until late May or early June, structural development of the male inflorescence is virtually indistinguishable from development of the female inflorescence. After that, the flower morphology, bract morphology, and cell division patterns differ markedly between the sexes [5]. Pioneering experiments found that a series of genes involved in flowering pathways are differently expressed between female and male flowers [6,7]. Although the molecular basis of poplar sex-specific flower development remains largely mysterious, increasing evidence suggests that miRNAs may play essential roles in this process.
In plants, there are two classes of endogenous small RNAs: small interfering RNAs (siRNAs) and microRNAs (miRNAs). Both classes, ranging from 20-24 nt, are short non-coding RNA molecules that negatively regulate gene expression at the transcriptional and/or post-transcriptional levels [8]. siRNAs are generated from double-stranded RNA. By contrast, miRNAs are transcribed from a long precursor molecule that folds back on itself to form a hairpin. Dicer-Like1 (DCL1) protein cleaves this precursor molecule, resulting in a miRNA:miRNA* complex, which separates into miRNA and miRNA* after transport to the cytoplasm [9]. In the cytoplasm, the miRNA is bound by ARGONAUTE proteins to form part of the RNA-induced silencing complex (RISC), which interacts with target mRNAs to cleave the RNA of target genes at the paired region [10]. Since the mature miRNA and its complementary target sequence have almost perfect complementarity, identifying a miRNA usually leads to the prediction and/or identification of its target. miRNAs have been shown to target genes during plant organ development, stress tolerance, phytohormone signaling, growth phase change, and disease resistance [11,12]. miRNAs have been identified as regulators of flower development and other pathways. For example, in Arabidopsis, the miR164 family targets a subset of NAC transcription factors including CUP-SHAPED COTYLEDON1 (CUC1) and CUC2, which mediate organ boundary formation. Among these members, miRNA164c affects petal number in early flowers [13]. miRNA172 regulates APETALA2 (AP2) and several of its homologs that share two tandem AP2 DNA-binding domains and play an important role in regulating flower development [14][15][16]. In maize, TASSELSEED 4 (TS4) and the AP2 homolog indeterminate spikelet 1 are confirmed miRNA172 targets. In ts4 mutants, female organs develop in the male inflorescence. Moreover, branching is increased in ts4 mutants indicating a link between sex determination and meristem fate [17]. In tree species, the unknown functions of the majority of genes are still the greatest gap. miRNA studies provided a best chance for directing to identify functions for genes [18]. Small RNA libraries prepared from poplar leaves and vegetative buds have been analyzed by high-throughput sequencing, finding 48 Populus-specific miRNA families [19]. However, sex-specific miRNA differences between female and male flower organs remain unclear, partly because of the effects of different genetic backgrounds of dioecious plants. For this work, we first used andromonoecious poplar, an exceptional model system that has both male and hermaphrodite flowers on a single plant, to eliminate confounding effects of genetic background of dioecious plants ( Figure S1). In this study, our aim is to screen sex-specific potential functional miRNAs in poplar floral tissue, investigate their defined expression patterns and demonstrate a possible functional network of sex-specific miRNAs in poplar flower organs.

Small RNA Sequence Analyses
Two small RNA (sRNA) libraries were constructed from female (F) and male (M) flowers of andromonoecious P. tomentosa. A total of 58,517,150 (F) and 47,918,351 (M) raw reads were generated from these libraries, respectively, using the Solexa sequencing technology (Table 1). After removing contaminant reads, clean reads were obtained and screened against rRNA, tRNA, snRNA, snoRNA, and mRNA in the Rfam 9.1 and NCBI GenBank databases, resulting in 50,782,479 (F) and 42,837,470 (M) reads remaining for further analyses. Conserved miRNAs were identified by alignment to data in miRBase 19.0 with 62 nt mismatch; a total of 4,789,550 (F) and 4,306,015 (M) reads that could form characteristic hairpin structures were identified. 4,459,932 (F) and 3,810,615 (M) unannotated unique sRNA sequences, which were further analyzed to predict novel miRNAs. The sequenced sRNAs were then mapped to the Populus genome, miRBase19.0, the NCBI GenBank database, and the Rfam database, and classified into six categories: rRNA, known miRNA (miRNAs in miRBase 19.0), exon, intron, repeat associated RNA, and unknown sRNA. More than 60% of the sequences were unknown sRNAs ( Figure 1A). The majority of total sRNA reads were either 21 or 24 nt for both libraries ( Figure 1B). In female flower libraries, the 21-nt sRNAs were the most abundant, making up 41.4% of total sequence reads. The 24-nt sRNAs were the second most abundant class at 17.3% of total sequence reads. In male flower libraries, the 21 and 24 nt sRNAs were not significantly different, making up 29.9% and 29.7%, respectively.

Conserved miRNAs in Andromonoecious P. tomentosa
To date, 46 conserved and non-conserved miRNA families have been discovered in the Populus genome [20]. In the two libraries sequenced in our study, 134 unique miRNA sequences were identified belonging to 38 conserved miRNA families ( Table 2). The expression levels of a few miRNA families, such as miR166 and miR472, were extraordinarily high in both libraries (Table 2). miR166 was the most abundant, with 2,157,665 (F) and 142, 916 (M) reads accounting for 62.5% and 47.1% of all conserved miRNA reads, respectively. Several miRNA families, such as miR156, miR159, miR169, miR319, miR396, and miR1447, had moderate expression levels ( Table 2). By contrast, some miRNA families showed very low levels of expression, with fewer than 200 reads. Different members in the same miRNA family displayed drastically different expression levels. For example, miR166 members varied in abundance from 695 to 2,157,665 reads. Also, 38 unique miRNA sequences belonging to 38 conserved miRNA families were produced from 134 loci (Table S1). Of these miRNA sequences, 11 were encoded by a single locus in the Populus genome, whereas the other 27 sequences had multiple loci. Most of these had 2-6 loci in the genome, and only a few had more than eight loci. For example, Pto-miR166a-l and Pto-miR169a-m had 12 and 13 loci, respectively. Thus, the size of miRNA families varies in andromonoecious P. tomentosa. Among these conserved miRNAs, miR6459 was only expressed in the female library. By contrast, two miRNAs, miRNA397 and miRNA6462, were detected as having male-specific expression ( Figure 2).

Identification of Novel miRNAs in Andromonoecious P. tomentosa
The characteristic hairpin structure of miRNA precursors was used to predict novel miRNAs. In total, 78 unique sRNA sequences were identified including 8 sRNAs with complementary miRNA*s (Table 3, 4 and 5). The majority of these novel miRNA candidates had lengths of 21 and 22 nt, and started with a 59 U ( Figure 3). The 78 sRNA sequences were transcribed from 85 loci; most were only produced from one locus, but Pto-F2, Pto-F13, Pto-F31 and Pto-F31 were each generated from 3 loci. The length of the predicted novel miRNA precursors ranged from 76 to 303 nt, with an average of 140 nt. The average minimum free energy (MFE) value was 255.58 kcal/mol, with a range of 221.2 to 2177.4 kcal/mol. The structures of the 78 novel miRNA precursors are shown in Figure S2. Most of them showed differential expression in the two libraries. For instance, the mature miRNA counts ranged from 63 to 8,480,860. Among those 78 miRNAs, only 6 sRNAs were found both in female and male libraries. 61 female-specific sRNAs and 11 male-specific sRNAs were detected, suggesting that more kinds of sRNA are involved in female floral development than in male development. To investigate whether these novel miRNA sequences were conserved across plant species, we used them as query sequences to perform BLASTN searches against the Arabidopsis thaliana, Brassica rapa, Citrus sinensis, Medicago truncatula, Oryza sativa, Ricinus communis, Prunus persica, Vitis vinifera and Zea mays genome databases in Phytozome v7.0. A large proportion of the novel sRNAs identified in andromonoecious P. tomentosa did not have perfect matches in different genome databases. Only 6 male-specific sRNA and 5 female-specific sRNAs were detected in other species (Table S2).
Novel miRNA and Targets were Located in Gender Determination Regions.
It was reported that Chromosome XIX is an incipient sex chromosome in Populus [21,22]. In the present study, 5 novel miRNAs located in chromosome XIX were identified, including Pto-F6, Pto-F70, Pto-F71, Pto-F72 and Pto-F73 (Table 6; Figure 4A). Pto-F6 was expressed both in female and male flower libraries, but expressed at higher levels in female flowers than in male ( Table 3). The Pto-70, Pto-71, Pto-72 and Pto-73 female-specific miRNAs were weakly expressed in female flowers (Table 5). A total of 33 targets located on chromosome XIX were identified for these novel miRNAs. Among these, Pto-F6 and Pto-73 were predicted to match 14 targets respectively that are located on chromosome XIX (Table 6). By contrast, no targets were identified for Pto-F72. Putative functions of these targets genes were annotated by BLAST searches against the TAIR, NCBI and JGI protein datasets. The results revealed that target genes are involved in disease resistance, response to auxin stimulus, decreased DNA methylation and membrane transport. Yin et al (2008) indicated that a 2.45 Mb segment in the peritelomeric region of chromosome XIX and is associated  with gender determination [21]. In our study, Pto-70 and 4 targets genes were found in sex-determination regions ( Figure 4B). Among the 4 target genes, 3 targets were predicted to match Pto-73 miRNAs and one was predicted to match Pto-F6 miRNA.
Target Prediction for Pto-miR156, Pto-miR159, Pto-miR172 and Pto-miR319 and Novel miRNAs in Andromonoecious P. tomentosa To understand the functions of sex-specific flower development related miRNAs, the first step is to predict and experimentally validate their targets. Consequently, we predicted 25 target genes for Pto-miR156, Pto-miR159, Pto-miR172 and Pto-miR319 that were different from the previously defined Populus miRNA targets (miRBase 19.0), and a total of 464 targets were identified for the 78 novel miRNAs (Table S3). The number of predicted targets varied from 1 to 15 per miRNA. The categorization of female and male-specific miRNA target genes according to biological process is shown in Figure 5. After putative functional annotation by BLAST searches against the TAIR, NCBI and JGI protein datasets, all targets of female and male-specific miRNAs that had a greater than twofold difference in expression between the sexes (Ratio §2or!0.5) were identified as candidate genes for gene functional enrichment analysis. All the target genes between the sexes were functionally characterized using enrichment analysis of gene sets, mostly as described by Gene Ontology terms ( Figure 5). GO analysis showed that response to hormone stimulus, negative regulation of gene expression, regulation of flower development, and organ formation were significantly enriched terms for biological processes.

Expression Profiling of miRNAs and Targets in Andromonoecious P. tomentosa
Andromonoecious poplar, an exceptional model system, eliminates confounding effects of genetic background of dioecious plants and thus is a powerful system in which to identify accurate sex-specific factors related to flower development. Therefore, the comparison of the small RNA profiles of female and male flowers of andromonoecious poplar was informative for characterization of functional miRNAs involved in sex-specific flower development. We sequenced 134 conserved miRNA sequences and 78 predicted novel miRNA sequences. 70.1% (94/134) of conserved miRNA sequences were up-regulated in female flower, and only 40 conserved miRNAs were up-regulated in male flower. A total of 104 conserved and 2 novel miRNAs were verified as differentially expressed by greater than two-fold between the two libraries (ratio .2 or ,0.5 and p,0.01), including 94 conserved and 2 novel miRNAs that were up-regulated in female flower and 37 conserved and 4 novel miRNAs that were up-regulated in male flower ( Table 2 and 3). Moreover, 3 conserved miRNAs (miRNA 397, miRNA 6459 and miRNA6462) and 72 novel miRNAs, including 11 male-specific and 61 female-specific, were expressed in a sex-specific manner ( Table 4 and 5). To examine the detailed expression pattern of these potential miRNAs in andromonoecious poplar flowers, the 27 differentially expressed miRNAs, including four conserved and 23 novel miRNAs, were verified by real-time quantitative PCR ( Figure 6A). Overall, comparisons of these data correspond to the deep sequencing results, indicating that the miRNA expression levels determined by high-throughput sequencing are reliable.
Base on the results of gene enrichment analysis, 31 putative target genes were selected as candidate genes for gene expression analysis (Table 7 and Figure 6B). The expression of miRNA target genes was detected using real-time quantitative PCR. For most target genes, expression was negatively correlated with the levels of a given miRNA, which is in accordance with the gene silencing  function of miRNAs. It should be noted that expression of CONSTANS-LIKE 6 (COL6), the target of PtoF25, was not correlated with the abundance of PtoF25 ( Figure 6B), indicating that COL6 expression might be regulated by other mechanisms.

A Possible Functional Network of Sex-specific Flower Development Related miRNAs in Andromonoecious P. tomentosa
The idea that miRNAs could play important roles in flower development has been suggested in previous studies. However, no systematic investigation has been conducted into how flower development may be regulated via miRNA-mediated processes. In the present study, a possible functional network, which connects flower development related miRNAs with their targets were first reported (Table 7 and Figure 7). Base on the gene functional annotation, all target genes were divided into four subgroups as shown in Figure 7.
Three miRNAs are involved in regulation of the second subgroup, genes that relate to Ca 2+ transport. Two female flower specific miRNAs, Pto-F28 and Pto-F45, were found to be negative factors for the expression of CALMODULIN-BINDING FAMILY PROTEIN and CALCIUM-BINDING EF HAND FAMILY PRO-TEIN. By contrast, CALCIUM-DEPENDENT PROTEIN KINASE FAMILY PROTEIN (CDPK) gene expression was repressed by Pto-F16, a male-specific miRNA.
In the third subgroup, there were 10 miRNAs involved in hormone related gene regulation. All of these miRNAs were female flower specific, except Pto-F6, which had expression significantly higher in female flower than in male. Pto-F68 has three target genes, including AUXIN EFFLUX CARRIER COM-PONENT 1, AUXIN-INDUCED PROTEIN (IAA4) and GIBBER-   In the last subgroup, genes that relate to DNA methylation were regulated by 3 miRNAs. MAINTENANCE OF METHYLATION1 (MET1) was negatively regulated by miRNA167, which had expression significant higher in male flower than in female. By contrast, two members of DECREASED DNA METHYLATION 1 (DDM1) were negatively regulated by Pto-F6, which had higher expression in female flower and by female flower specific Pto-F56, respectively.

Overview of the Deep Sequencing Datasets
Although an increasing number of studies about miRNAs and their expression patterns from different tissues have been reported, there are still few studies on sex-specific miRNAs between female and male flowers of poplar. In this study, andromonoecious poplar, an exceptional model system, was used to eliminate confounding effects of genetic background of dioecious plants. We generated sRNA libraries from female and male flower of andromonoecious poplar and obtained over 40 million reads of 18-30 nt per library, substantially increasing the available data on poplar flower sRNAs [23,24]. Mapping analysis of unique sRNAs revealed that over 60% of the identified sRNAs were previously unknown, suggesting that even more sRNAs remain to be identified. The size distribution of small RNAs indicated that the numbers of 21-nt sRNAs were significantly higher than 24-nt sRNAs in female flower, suggesting that female flower development may require more sRNAs in this size class. By contrast, the numbers of 21-nt sRNAs and 24-nt sRNAs were not significantly different in male flower libraries. The 24-nt sRNAs mainly comprised siRNAs associated with repeats and transposons [25,26]. 24-nt sRNA numbers in male flower were significantly higher than in female, suggesting that a number of siRNAs related to repeats and transposons are involved in male flower development.

Conserved miRNAs Expressed in Flower
We detected dramatically different expression levels among different members of the same miRNA family, suggesting that expression of conserved miRNAs varies in andromonoecious P. tomentosa. These conserved miRNAs have been identified in multiple flowering pathways. For example, three miRNA families (miR172, miR159/miR319 and miR156) are involved in flowering-time regulation [27]. In this study, miRNA159 and miRNA319 expression in female flower was 38.8-fold and 51.42fold higher, respectively, than in male flower. miRNA172 was weakly expressed in female and male flower, with no significant difference. By contrast, miRNA156 was expressed in male flower at higher levels than in female flower. This result suggested that different miRNAs are involved in female and male flowering time regulation. miRNA159 and miRNA319 also have been reported to act with miRNA164 and miRNA167 in specifying particular cell types during the later stages of flower development [28]. The present study revealed that miRNA164 was not expressed in female or male flower, indicating miRNA164 might show stagespecific expression in flower development. miRNA167 was upregulated dramatically in male flower, indicating that it may play an important role in specifying male flower cell types. miRNA156, miRNA157 and miRNA172 may be components of a regulatory pathway mediating the transition between the vegetative and reproductive phases in plants [28]. miRNA166, targeting APETALA2 and type III homeodomain-leucine zipper (HD-Zip) family genes, was dramatically up-regulated in female flower, suggesting AP2 and HD-ZiP genes might be negatively regulated in female flower [29]. miR167 and miR160 are thought to control transcription in response to the phytohormone auxin by targeting mRNAs coding for ARF DNA-binding proteins [30,31]. These miRNAs were detected in this study as having opposite expression patterns between female and male flower, suggesting that two miRNAs involved in regulation of auxin responses in female flower are different compare with in male flower.

Novel miRNAs Displayed Different Expression Patterns between Female and Male Flower
The 78 novel miRNAs identified here spanned a size range of 20-22 nt, and the majority started with a 59 U, consistent with results in Japanese apricot [32]. The majority of novel miRNAs showed relatively low expression levels, which was consistent with previous studies predicting that novel miRNAs are often expressed at lower levels than conserved miRNAs [33]. The miRNA*s were generally less abundant than the corresponding mature miRNAs. It has been reported that miRNA*s have important functions through regulation of the miRNA precursors themselves when miRNA*s were more abundant than their corresponding miRNAs [34]. miRNA*s identified were used as extremely strict criteria to identify novel miRNAs [35]. In our study, only 8 unique miRNAs were detected with complementary miRNA*s. Most novel miRNAs do not have identified miRNAs*, likely because of the depth of sequencing, which results in a low probability of identifying low-frequency miRNA*s. The numbers of femalespecific novel miRNAs was 6-fold higher than male-specific novel miRNAs, suggesting that abundant novel miRNAs are involved in female flower development. A large proportion of the novel sRNAs identified in andromonoecious P. tomentosa did not have perfect matches in different genome databases, providing more evidence to support our conclusion that these novel miRNAs might be poplar-specific and were first reported in present study. Among all of the novel miRNAs, 11 miRNAs were found with perfect matches in other species, suggesting that some of these newly identified miRNAs may be broadly conserved within angiosperms.

Novel miRNAs Located in Sex Determination Regions
Although studies on determination of gender indicate that dioecy in Populus is genetically controlled, the precise genderdetermining systems remain unclear. Based on genetic mapping results, a gender determination locus mapped to different positions on chromosome XIX, depending upon the Populus species [36][37][38]. Klevebring et al. (2009) found that the proposed sexdetermining peritelomeric region of chromosome XIX showed a distinctive pattern of sRNA occurrence that differed significantly from the rest of the genome [39]. In our study, 5 novel miRNAs and 33 targets genes were located on chromosome XIX, indicating that these miRNAs might be involved in sex determination of poplar. The annotated results revealed that over one-third of targets were annotated as disease resistance genes. It has been reported that a NBS-LRR gene supercluster was located on Populus chromosome XIX. These results are consistent with the hypothesis that resistance to a floral pathogen and regulation of gender determination coevolved in a region of the genome that experiences reduced recombination [22]. Interestingly, Pto-F6 was predicted to match a target that encodes DNA (cytosine-5)methyltransferase implying that DNA methylation might act in flower development. To better understand the functions of identified miRNAs, we sought to identify potential targets in andromonoecious P. tomentosa using computational analyses. Many targets of conserved miRNAs have been predicted and their functions annotated [40,41]. Integrating target genes and miRNA expression data identified sub-networks for further investigation. COL1 encodes a zinc finger protein belonging to the first groupings of the flowering-time gene CONSTANS (CO) [42]. Over-expression of COL1 can shorten the period of two distinct circadian rhythms by affecting a light input pathway [43]. In our study, COL1 was targeted by female flower specific Pto-F66, which was expressed in female flowers suggesting that the novel miRNA Pto-F66 likely plays a negative role in the light input pathway in female flower. By contrast, COL6 is part of the third groupings of CO, and is a target of Pto-F19 and Pto-F25. However, these two miRNAs were not correlated with COL6 transcript abundance, indicating that COL6 expression might be regulated by other mechanisms. The target of male specific miRNA Pto-F11, EPR1 encodes a nuclear-localized MYB protein showing rapid induction in response to red light irradiation. EPR1 was significantly down-regulated in male flower, suggesting that this primary phytochrome-responsive factor is likely repressed by the male flower specific novel miRNA Pto-F11 [44]. Another gene involved in the photoperiodic pathway, PFT1 transcript abundances were significantly decreased in male flower. PFT1 is a target of PtoF14 and functions downstream of phyB to regulate the expression of FLOWERING LOCUS T (FT), providing evidence that Pto-F14 might be involved in regulation of flowering time by affecting the light quality pathway [45].
The AP2 gene family contains 144 members encoding at least one AP2 DNA-binding domain and the varied biological functions of AP2 family members range from development to stress and defense responses [46][47][48]. Conserved miRNA172 has been identified as an important factor regulating AP2 transcription factors [17]. In our study, miRNA172 was equally expressed in female and male flower, but two novel miRNAs, Pto-F36 and Pto- F51, were detected with female-specific expression. AP2 transcript abundance decreased in female flower, providing evidence that these two miRNAs act in regulation of AP2. UFO encodes an Fbox protein required for the determination of floral-organ and floral-meristem identity [49]. PM36 is involved in reproductive system development. These genes were targeted by the femalespecific novel miRNA Pto-F54 and UFO and PM36 expression were down-regulated in female flower, indicating the broad biological functions of Pto-F54. Calcium is a ubiquitous second messenger in plant signal transduction cascades [50]. Three novel miRNAs were identified as related to Ca 2+ transport, including one male flower specific and two female flower specific miRNAs, suggesting that Pto-F16, Pto-28 and Pto-45 are likely play different roles in regulating calcium transport in female and male flower development.
Hormones play a central role in the coordination of internal developmental processes with external environmental signals [51]. A series of genes related to phytohormone synthesis and metabolism are differentially expressed between female and male during flower development (unpublished data). In our study, 13 genes related to five important phytohormones were found to be targeted by 12 novel miRNAs, suggesting that these miRNAs act broadly in regulation of phytohormone synthesis and metabolism. Previous work reported that auxin is necessary for GA-mediated control of root growth, and that attenuation of auxin transport or signaling delays the GA-induced disappearance of RGA (repressor of gal-3) [52]. Among these miRNAs, Pto-F68 targets three genes involved in auxin transport and gibberellin synthesis suggesting that Pto-F68 may regulate the interaction of auxin and gibberellin. Also, a series of physiological studies indicate that ethylene and auxin interact in a number of different ways depending on the cell type, developmental stage, and environmental conditions [51]. EREBF4 and SAUR2 expression were downregulated in female flower and targeted by Pto-F47 miRNA. This result may provide a clue that deepens our understanding of the important regulatory functions of miRNA Pto-F47 in phytohormone interactions. Poplar male flowers abscise earlier than female flowers, suggesting the ABA accumulation and regulation of related genes might be different between female and male flowers. ABA3 encodes an enzyme that catalyzes the generation of the sulfurylated form of MoCo, a cofactor required by aldehyde oxidase, which functions in the last step of ABA biosynthesis in plants [53]. ABI3 encodes a DNA binding factor targeted to promoters responsive to ABA and auxin [54]. These two genes related to ABA biosynthesis and responses are targeted by novel miRNAs Pto-F55 and Pto-F65, which have female flower specific expression, suggesting these miRNAs have negative roles in flower abscission through regulation of ABA metabolism.
MET1, encoding DNA (cytosine-5) -methyltransferase, which maintains CG methylation, was down-regulated in male flower [55]. DDM1, encoding a chromatin remodeling complex with ATPase activity that can remodel nucleosomes, was up-regulated in male flower [56]. MET1 and DDM1 play opposite roles in the regulation of DNA methylation levels; therefore, the different expression of these genes in female and male flowers might alter DNA methylation levels, which are higher in female flower than male (unpublished). Our study indicated that miRNA167, Pto-F6 and Pto-F56 might play a positive regulatory role in maintaining DNA methylation levels of female flower through regulation of MET 1 and DDM1.

Plant Materials
Four-hundred and sixty unrelated individuals, representing almost the whole geographic distribution of P. tomentosa, were collected in the national nursery of Guan Xian County (Shandong Province, China). After observation of flower morphology, three 29-year-old unrelated andromonoecious clones ('2-14', '3605', and '5103') were found.
All three individuals were used in our study. Sampling that female and male flowers were dissected from the andromonoecious flowers was performed at the last phase of flower development, before pollination. All materials were collected from poplar individuals and immediately frozen in liquid nitrogen, and stored at 280uC for isolation of RNA. To eliminate genetic backgroud, the two bulks were made by mixing equal amounts of RNA extracted from three andromonoecious poplar female and male flower, respectively.
This study was carried out in strict accordance with the recommendations in the Guide for Observational and field studies (http://www.plosone.org/static/publication). All necessary permits were obtained for the described field studies. The sampling of all individuals of P. tomentosa was approved by Youhui Zhang, director of National Garden of P. tomentosa in Guan Xian County, Shandong Province.
Total RNA Isolation, sRNA Library Construction, and Solexa Sequencing Total RNA was isolated from female and male floral tissue by a modified CTAB method [57] with isopropanol instead of lithium chloride for RNA precipitation. Extracted RNA was used for sRNA library construction. Small RNAs were sequenced using an Illumina HiSeq 2000 at the Shanghai Bio Institute.
Sequence data from this article have been deposited in the GenBank Data Library under the accession Nos: KC477266-KC477295.

Bioinformatic Analyses of Sequencing Data and Target Prediction of miRNAs
Clean reads were screened from raw reads by removing contaminating reads, including sequences with adapters, without the insert tag, with poly-A tails, shorter than 18 nt, or longer than 30 nt. Clean reads were then aligned with the Populus genome (http://www.phytozome.net/poplar) using SOAP [58], with no mismatches allowed. Small RNA tags matching exons and introns of mRNAs were excluded from further analysis. Tags matching non-coding RNAs, including rRNAs, scRNA, tRNAs, snRNAs, and snoRNAs, deposited at the Rfam (http://www.sanger.ac.uk/ Software/Rfam) and NCBI GenBank databases (http://www. ncbi.nlm.nih.gov/blast/Blast.cgi) [59] were also excluded from further analysis. The remaining unannotated sRNAs were searched against miRBase 19.0 with a maximum of two mismatches allowed [20] to identify conserved miRNAs in P. beijingensis and then the resulting sequences were screened for the presence of the characteristic hairpin structures using the program RNAfold (http://www.tbi.univie.ac.at/*ivo/RNA/ViennaRNA-1.8.1.tar.gz) [60]. We used the prediction software Mireap (https://sourceforge.net/projects/mireap/) to predict novel miR-NAs by exploring the secondary structures, DL1 cleavage sites, and minimum free energies of the unannotated sRNA tags that could be mapped to the Populus genome. Basic criteria [35] were used for selecting the potential novel miRNAs. Target genes were predicted in accordance with previous research [61,62] but with higher stringency; only three mismatches were allowed in the miRNA/target duplex. Target searches were performed using Populus transcripts deposited in the Phytozome v7.0 (http://www. phytozome. net/) and Pfam databases (http://pfam.sanger.ac.uk/), and the results were used to annotate the functions of potential targets. We compared miRNA expression between the two libraries (control and treatment) to determine which miRNAs were differentially expressed. We normalized the expression levels of miRNAs within each library to get the expression of transcripts per million.

59-RACE
RNA Ligase-Mediated 59-RACE (RLM-RACE) was performed with the First Choice RLM-RACE Kit (Ambion) according to the manufacturer's instructions, with slight modifications. Briefly, 10 mg of total RNA was directly ligated to the 59 adaptor followed by reverse transcription with an oligo (dT) primer. PCR was performed with 59 adaptor primers and 39 gene-specific primers (Table S4) using cDNA as the template. The RACE products were gel-purified, cloned, and sequenced.

qRT-PCR
Small RNAs (,200 nt) were isolated using the mirVana miRNA Isolation Kit (Ambion, USA) following the manufacturer's instructions. Quantitative real-time PCR (qRT-PCR) was carried out as previously described [63] in an ABI PRISM 7500 Fast Realtime PCR System (Ambion, USA) using the SYBR Premix Ex TaqTM Kit (TaKaRa, Japan). The reactions were carried out in a 20 ml volume containing 2 ml of diluted cDNA, 200 nM of each primer, and 19 PCR Master Mix with the following conditions: 95uC for 30 s, and 45 cycles of 95uC for 5 s, 58uC for 15 s, and 72uC for 20 s. Then, a thermal denaturing cycle of 95uC for 15 s and 60uC for 1 min was applied to determine the dissociation curves, which were used to verify the specificity of PCR amplifications. All reactions were run in triplicate for each sample. Twelve miRNAs, including six conserved and six novel miRNAs were validated, and 5.8S rRNA was selected as a reference gene for normalization [41] (Table S5 and Table S6).

Network Analysis
The details of sex-specific flower development related miRNAs functional network methods were done according to previously described [64]. Comparative miRNAome analysis combined with realtime quantitative PCR and RACE-PCR revealed 23 miRNAs expressed in andromonoecious poplar flower and each of them has functional targets. Negative correlation between the miRNA and target expression profiles was validated using real-time quantitative PCR. These results indicated that miRNA-mediated regulation does occur in flower development and that these miRNAs with dynamic changes ultimately down-regulate their targets. All miRNA and functional targets were connected for network analysis.