High-Throughput Sequencing of Small RNAs from Pollen and Silk and Characterization of miRNAs as Candidate Factors Involved in Pollen-Silk Interactions in Maize

In angiosperms, successful pollen-pistil interactions are the prerequisite and guarantee of subsequent fertilization and seed production. Recent profile analyses have helped elucidate molecular mechanisms underlying these processes at both transcriptomic and proteomic levels, but the involvement of miRNAs in pollen-pistil interactions is still speculative. In this study, we sequenced four small RNA libraries derived from mature pollen, in vitro germinated pollen, mature silks, and pollinated silks of maize (Zea mays L.). We identified 161 known miRNAs belonging to 27 families and 82 novel miRNAs. Of these, 40 conserved and 16 novel miRNAs showed different expression levels between mature and germinated pollen, and 30 conserved and eight novel miRNAs were differentially expressed between mature and pollinated silks. As candidates for factors associated with pollen-silk (pistil) interactions, expression patterns of the two sets of differentially expressed miRNAs were confirmed by stem-loop real-time RT-PCR. Transcript levels of 22 predicted target genes were also validated using real-time RT-PCR; most of these exhibited expression patterns contrasting with those of their corresponding miRNAs. In addition, GO analysis of target genes of differentially expressed miRNAs revealed that functional categories related to auxin signal transduction and gene expression regulation were overrepresented. These results suggest that miRNA-mediated auxin signal transduction and transcriptional regulation have roles in pollen-silk interactions. The results of our study provide novel information for understanding miRNA regulatory roles in pollen-pistil interactions.


Introduction
Appropriate pollen-pistil interactions are essential to subsequent successful reproduction and seed formation in angiosperms [1]. At the beginning of the reproductive process, mature pollen grains are released from the dehisced anthers and land on the receptive stigmatic surface. After pollen capture, compatible pollen grains hydrate and germinate to extrude pollen tubes. Elongated by tip growth, pollen tubes penetrate the cell layers of the stigma, navigate within the style tissues, and eventually reach the ovule where fertilization occurs [2]. Incompatible pollen grains, on the other hand, are arrested at certain points in their journey toward the ovule [3]. Pollen-pistil interactions are a series of intensively regulated processes that ensure delivery of appropriate male gametes to female ones [4]. During these processes, pollen controls germination and tube growth, while the pistil provides the foundation and directional guidance for pollen tube navigation [5]. Because pollen-pistil interactions are essential to seed and fruit formation, revealing that their mechanisms is of great importance, both for understanding the completion of the plant life cycle and for enhancing agricultural production. Recent transcriptomic and proteomic studies have extended our knowledge of pollen/pistil gene and protein expression and revealed genes possibly involved in pollen-pistil interactions [2,[6][7][8]. The post-transcriptional regulation underlying these interactions remains largely unknown, however.
MicroRNAs (miRNAs), which are negative regulators of gene expression, are involved in many biological processes in eukaryotic species [9]. Most miRNA genes are transcribed by RNA polymerase II into primary miRNAs (pri-miRNAs), which are converted in turn into single-strand stem-loop precursors (pre-miRNAs) and then miRNA:miRNA* duplexes by proteins such as Dicer like 1, HYPONASTIC LEAVES1, and SERRATE in plants [10]. miRNAs are then incorporated into RNA-induced silencing complexes (RISCs) and cleave their target mRNAs with the assistance of ARGONAUTE proteins, while miRNA*s are usually degraded [11]. To elucidate the post-transcriptional regulation underlying various biological processes, miRNA expression profiles during these processes have recently been analyzed in different plant species [12][13][14][15]. No studies identifying miRNAs with potential roles in pollen-pistil interactions are available, however.
Maize (Zea mays L.) is a model species for studying pollen-pistil interactions, and is one of the most important crops in the world [16]. In this study using maize as the material, we sequenced four small RNA libraries obtained from mature pollen (MP), in vitro germinated pollen (GP), mature silks (pistils; MS), and pollinated silks (PS). Two sets of differentially expressed miRNAs were identified and considered as candidate factors related to pollen-silk interactions. The first set included miRNAs differentially expressed during pollen germination and tube growth, and the second set included miRNAs responsive to pollination in silks. GO analysis of their target genes revealed the enrichment of terms related to auxin signaling and transcriptional regulation, suggesting the potential regulation of these processes by miRNAs during pollen-pistil interactions.

Results and Discussion
High-throughput sequencing of small RNAs from maize pollen and silk tissues To profile small RNAs expressed in maize pollen and silk, and to identify miRNAs potentially involved in regulation of pollen-silk interactions, we sequenced four small RNA libraries from mature pollen (MP), in vitro germinated pollen (GP), mature silks (MS), and pollinated silks (PS) using Solexa technology. We obtained 6,433,541, 11,133,002, 9,173,416, and 8,327,428 raw reads from MP, GP, MS, and PS, respectively. After removing adaptors, contaminants, and low-quality sequences, 3,638,941, 10,575,042, 5,775,303, and 5,205,129 high-quality clean reads 18-30-nt long were generated from MP, GP, MS and PS, respectively. Length distribution based on total number of clean reads from the four sequenced libraries is illustrated in Figure 1. Small RNAs with 24 and 22 nt were the two major classes, similar to previous findings in many tissues of maize and rice [13,14,17].
The clean reads were mapped to the maize genome (B73 RefGen_2v, released 5a.59) using SOAP, revealing 3,107,037 (85.38%), 9,527,191 (Table 1). Based on a Blastn search against Rfam base, sequences corresponding to rRNAs, tRNAs, snoRNAs, and other small RNAs were identified (Table 1). A large number of  sequences (22,969,291 representing 9,816,635 unique sequences) could not be annotated. The well-represented unannotated small RNAs provided an opportunity to identify novel miRNAs.

Identification of conserved miRNAs
To identify known maize miRNAs, candidate sequences were Blastn-searched against known mature miRNAs and their precursors in miRBase (version 18.0) to recover perfectly matching sequences from maize. We identified 161 conserved miRNAs belonging to 27 miRNA families from the four libraries: 97 in MP, 105 in GP, 143 in MS, and 147 in PS (Table S1). Total read numbers of conserved miRNAs in the four libraries were 94,487 (MP), 168,406 (GP), 643,011 (MS), and 580,305 (PS) ( Table S1).
We further analyzed conserved miRNAs with relatively high (.1,000) read numbers. As shown in Figure 2, zma-miR156a/b/ c/d/e/f/g/h/i/l and miR168a/b were highly abundant in all four libraries. miR156 is involved in diverse processes in various model species, including anther development, juvenile-to-adult phase transition, gibberellin signal transduction, and flower and leaf development [18][19][20][21]. In maize, however, the accumulation of zma-miR156 in pollen and silk samples suggests its possible role in reproductive processes. miR168 is a known feedback regulator of the miRNA pathway, where it targets ARGONAUTE1 (AGO1), which is important for miRNA function [22,23]. The modulation of AGO1 mediated by miR168 is essential to proper development in Arabidopsis [22]. We suggest that miR168-controlled feedback regulation is also required for plant reproduction. Members of several other miRNA families, including zma-miR164, zma-miR167, zma-miR169, zma-miR171, and zma-miR827, were also highly expressed in MS and PS; they were barely detected in MP and GP, however, indicating their specific roles in maize female reproductive tissues.
To further understand the potential biological roles of identified conserved miRNAs, we predicted their target genes using a previously described approach [24]. Based on this method, 144 target genes were identified (Table S3).

Prediction of novel miRNAs
In plants, miRNAs are excised from the stem of their singlestranded precursors' stem-loop structures [25]. Based on this fundamental defining miRNA characteristic, we predicted novel miRNAs using the following criteria. First, flanking sequences of the unannotated small RNAs were analyzed with the Einverted program of the Emboss software package, and secondary structures were predicted using RNAfold [26][27][28]. Second, MirCheck analysis was applied to filter out sequences mapping to the loop region. Sequences passing MirCheck were further inspected manually [27]. Third, we removed putative miRNAs with minimal folding free energy indexes # 0.85, based on previous studies [29]. Finally, to minimize noise, only small RNAs with more than five reads in at least one library were selected. A total of 82 novel miRNA sequences classifiable into 61 families were identified (Table S2). The star sequences of some miRNAs were also identified. The lengths of these newly identified miRNAs ranged from 20 to 24 nt. Their negative folding free energies varied from 219.9 to 2171.2 kcal mol 21 , with an average of 265 kcal mol 21 , similar to values uncovered in other recent studies on maize [17,30]. Unlike conserved miRNAs, novel miRNAs were nearly equally distributed across the four libraries, with 51 in MP, 58 in GP, 69 in MS, and 52 in PS. Target genes of these novel miRNAs were predicted using the same approach employed for the conserved miRNAs; 44 possible targets were identified (Table S4).
For some novel miRNAs, both sequences from 3p and 5p arms were identified (Table S2). Most read counts were similar between the two arms, indicating accumulation of both miRNAs and miRNA*s. Similar observations have been made in many other tissues of maize [14,17]. It is currently believed that some miRNA*s may be able to accumulate in large quantities and exert the same functions as miRNAs-regulating their target genes-in both plants and animals [31][32][33][34]. Because such miRNA*s have been referred to as de facto miRNAs [14], we considered all the newly identified 3p-and 5p-sequences in our dataset to be novel miRNAs.

Identification of miRNAs differentially expressed between mature and germinated pollen
Identification and characterization of miRNAs differentially expressed between MP and GP may facilitate our understanding of post-transcriptional regulation during pollen germination and tube growth. For this purpose, we compared normalized expression levels of miRNAs from MP and GP, selecting miRNAs with fold-changes greater than 1.5 (log 2 ratio) and P-values less than 0.01 (Chi 2 test). As shown in Table 2, 40 conserved and 16 novel miRNAs were found to be differentially expressed between MP and GP.
Five conserved and 12 novel miRNAs, including zma-miR159h/i, zma-miR390a/b and novel-a-5p, were up-regulated in GP compared with MP (Table 2). Previous studies have shown that spatial and temporal regulation of miR159a/b/c by APC/C   is required for mitotic progression in male gametophytes in Arabidopsis [35]. On the other hand, a mir159ab double mutant did not exhibit defects in pollen germination [36]. These results suggest that miR159a/b/c participates in pollen development, but not in pollen germination and tube growth. In our study, different members of the miRNA159 family displayed diverse expression patterns. miR159a/b/c expression levels varied moderately between MP and GP, whereas miR159h/i exhibited significant changes, indicating the probable involvement of miR159h/i in pollen germination and tube growth. These results suggest that individual members of the same miRNA family may execute their respective functions through different expression patterns. Another miRNA up-regulated in GP, miR390a/b, has been shown to respond to the concentration of auxin, a stimulator of pollen tube growth, and to mediate its signal transduction [37,38]. In addition, novel-a-5p was predicted to target SNARE-associated Golgi proteins involved in vesicle trafficking during pollen tube elongation [39,40]. It is thus possible that zma-miR390a/b and novel-a-5p participate in regulation of maize pollen tube elongation. Thirty-five conserved and four novel miRNAs were downregulated in GP, of which 11 conserved and three novel miRNA sequences were detected with read counts greater than 10 in at least one library (Table 2; Tables S1 and S2). Several of these down-regulated members may be involved in pollen tube growth. For instance, miR167 controls the expression of ARF6 and ARF8, which are essential for pollen tube elongation [41]. Overexpression of MIR167b in Arabidopsis leads to defective pollen germination [42]. We suggest that the observed down-regulation of miR167 in GP ensures normal ARF6 and ARF8 expression and promotes pollen germination and tube growth. Another miRNA down-regulated in GP, miR166, influences establishment of adaxial/abaxial leaf polarity in maize, and regulates shoot apical meristem and floral development in Arabidopsis [43][44][45]. Finally, miR396 mediates leaf development in Arabidopsis by controlling GROWTH-REGULATING FACTOR [46]. In maize, however, the obvious down-regulation of miR166f/g/h and miR396a/b in GP compared with MP implies that their targets are probably involved in pollen tube growth.

Characterization of miRNAs differentially expressed preand post-pollination in silks
Although pollinated silks were harvested after shaking off unadhered pollen grains, PS was still a mixture of pollinated silks, germinated pollen, and a few pollen grains. In spite of their low percentages in PS, residual pollen grains and tubes could interfere with identification of miRNAs differentially expressed between MS and PS before and after pollination. More specifically, miRNAs accumulating at high levels in MP and/or GP might be considered to be up-regulated in PS compared with MS even though their expression levels were unchanged between mature and pollinated silks. In addition, miRNAs varying moderately in silks pre-and post-pollination might be identified as significantly down-regulated in PS because of their lower levels in PS compared with those in MS. To avoid these errors, we used the following criteria to select miRNAs differentially expressed between MS and PS. First, miRNAs displaying expression level increases more than 1.5-fold (log 2 ratio) higher in PS than in MS, and with more than 10-fold higher read counts in PS or MS than in GP or MP were considered to be up-regulated in PS. Second, miRNAs with more than 2-fold (log 2 ratio) decreased expression in PS compared with MS were considered to be down-regulated in PS. On this basis, 30 conserved and eight novel miRNAs were identified as differentially expressed in silks pre-and post-pollination ( Table 3).
Previous studies have found that pollen rehydration processes on dry stigmatc surfaces may cause water deficit stress and alter expression of genes responsive to drought stress [7,59]. Moreover, pollen tube reception and fungal invasion have been shown to share conserved components, indicating similar molecular mechanisms in these two processes [61]. Consistent with these observations, all conserved miRNAs identified as differentially expressed pre-vs. post-pollination in this study were involved in drought and/or fungal invasion responses. We thus assume that these miRNAs play a role in pollen-silk interactions; however, dissecting the underlying molecular mechanisms requires further investigation.

Analysis of miRNAs potentially involved in pollen-silk interactions
miRNAs associated with pollen-silk interactions could be derived from either pollen or silk. In this study, we identified miRNAs differentially expressed during in vitro pollen germination and tube growth, as well as from silks pre-and post-pollination. Most current information about pollen germination and tube growth has been obtained from in vitro studies [62]. Many common cellular and molecular events, as well as genes regulating pollen germination and tube growth, are shared between in vivo and in vitro processes [8,62]. We therefore assumed that most of the differentially expressed miRNAs identified from in vitro germinated pollen also play a role in in vivo processes. It was thus likely that the two sets of differentially expressed miRNAs, including 72 conserved and 24 novel ones, contained miRNAs functioning in pollen-silk interactions (Tables 2 and 3). These differentially expressed miRNAs were subsequently considered as potential factors involved in pollen-silk interactions.
We used stem-loop real-time RT-PCR to confirm expression patterns of the two sets of differentially expressed miRNAs uncovered via Solexa sequencing. According to our RT-PCR analysis, expression patterns of these differentially expressed miRNAs were all similar to those detected by sequencing (Figures 3, 4A and B). Expression levels of a few memberszma-miR156k, zma-miR166m, zma-miR167a/b/c/d, novel-23-5p, and novel-27-3p-were slightly different ( Figure 3). This difference may be due to variations in MP and GP sampling times, or differences in sensitivity and specificity of the two technologies.
In general, results obtained by stem-loop real-time RT-PCR were in accordance with the sequencing data. In addition, we further validated the expression level of zma-miR171b/f using northern blot analysis ( Figure 4C). The result was consistent with those of stem-loop real-time RT-PCR and sequencing data. We also analyzed expression of 22 predicted target genes of seven selected miRNAs through real-time RT-PCR. As shown in Figure 5, expressions of most of these putative target genes were negatively correlated with expression of their corresponding miRNAs, indicating that these genes might be repressed by miRNAs through transcriptional regulation. GRMZM2G081406, GRMZM2G153233, and GRMZM2G475882, however, were positively correlated with their miRNAs ( Figure 5). It is possible that these three genes were regulated through translational repression by miRNAs.
To investigate their functions and study their potential roles in pollen-silk interactions, we performed Gene Ontology (GO) analysis on predicted target genes of the differentially expressed miRNAs using singular enrichment analysis (SEA) [63]. Six GO terms were found to be enriched (Table 4). Within the GO biological process category, the terms ''response to hormone stimulus'', ''regulation of gene expression'', and ''macromolecule biosynthetic process'' were overrepresented (Table 4).
Notably, all seven genes enriched in the ''response to hormone stimulus'' category-GRMZM2G005284, GRMZM2G078274, GRMZM2G081406, GRMZM2G153233, GRMZM2G159399, GRMZM2G390641, and GRMZM2G475882-were auxin response factors (Table S3), indicating that this hormone might be involved in pollen-silk interactions. The functions of miR167regulated ARF6 and ARF8 in pollen germination and tube growth have been discussed above. Real-time RT-PCR analysis revealed that the expression level of GRMZM2G078274, which encodes ARF6, was negatively correlated with that of zma-miR167a/b/c/ d ( Figure 5; Table S3). In addition, zma-miR160a/b/c/d/e/g and zma-miR393a/b/c were observed to be down-regulated in silks after pollination ( Figure 4B; Table 3). Recent studies in model species have provided evidence that miR160 targets ARF10, ARF16, and ARF17, and that miR393b targets TIR1/AFB2, the gene encoding an auxin receptor [64,65]. Our real-time RT-PCR results revealed that expression levels of GRMZM2G005284, GRMZM2G159399, GRMZM2G390641, and GRMZM5G808366, which encode ARF10, ARF16, ARF16, and ARF17, respectively, were elevated when zma-miR160a/b/ c/d/e/g was down-regulated ( Figure 5; Table S3). Additionally, three targets (GRMZM2G135978, GRMZM2G137451, and GRMZM5G848945) encoding auxin receptors exhibited opposite expression patterns to zma-miR393 ( Figure 5; Table S3). Observed down-regulation of zma-miR160 and zma-miR393 therefore implies increased accumulation of ARF10, ARF16, ARF17, and TIR1/AFB2 during pollen-silk interactions, indicating activation of auxin signal transduction. miR393 is suggested to function as an integrator of complex environmental cues through regulation of auxin receptors [66]. We propose that during pollination processes, zma-miR393 responds to signals released from pollen grains and tubes by activating auxin signal transduction, which enables silk cells to achieve suitable states for acceptance and guidance of elongating pollen tubes. Reduced fertility observed in the mutant miR393b-1 supports the notion that zma-miR393 plays an important role in pollen-silk interactions [66].
Moreover, the terms ''regulation of gene expression'', ''DNA binding'', and ''nucleus'' were overrepresented in GO biological processes, molecular functions, and cellular components categories, respectively (Table 4). Taken together, these results indicate that miRNAs may mediate transcriptional regulation through their target genes during pollen-silk interactions.
In conclusion, we obtained 22,065,744 genome-matched small RNA sequences representing 9,859,314 unique sequences from maize pollen and silk samples. Among these sequences, 161 were identified as known miRNAs belonging to 27 families, and 82 were predicted as novel miRNAs. We characterized identified miRNAs differentially expressed between MP and GP, and between MS and PS. Expression patterns of these differentially expressed  miRNAs, candidates for involvement in pollen-silk interactions, were confirmed by stem-loop real-time PCR. Regarding miRNAs differentially expressed between MP and GP, several members showed potential involvement in pollen germination or tube growth. In silks, miRNAs responding to pollination have been found in previous studies to be responsive to stresses, especially drought and fungal invasion. GO analysis of putative target genes of the two sets of differentially expressed miRNAs indicated that auxin signal transduction and transcriptional regulation controlled by miRNAs may be involved in the modulation of pollen-silk interactions.

Plant materials
Maize (Zea mays L.) inbred line B73 was used in this study. Plants were grown in the greenhouse at 22uC and 60% relative humidity under a 16h/8h light/dark cycle. Mature pollen (MP) was collected by shaking tassels between 9:00-10:00 a.m. To ensure viability of collected pollen, we tested their in vitro germination percentages. Only pollen grains with germination percentages greater than 90% were used. Germinated pollen (GP) was obtained after culturing MP for 30 min in liquid medium [67], with germination percentages greater than 90% obtained (data not shown). Silks extending approximately 5 cm beyond the bracts were harvested as mature silks (MS). To obtain pollinated silks (PS), mature silks were pollinated by an excess of active pollen. After 30 min, most pollen grains had germinated and their tubes had penetrated into silk tissues (data not shown). The pollinated silks were then harvested after shaking off unadhered pollen grains.

Small RNA library construction and Solexa sequencing
Library construction and sequencing were performed as described previously [68]. Briefly, total RNA was isolated using Trizol reagent (Invitrogen, Carlsbad, CA, USA). Small RNAs ranging from 18-30 nt were gel-purified and ligated to proprietary adaptors at 59 and 39 termini. RT-PCR was performed for 18 cycles, and the products were purified and subjected to Solexa sequencing. Sequencing data is available in the GEO database under accession number GSE44787 (http://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc = GSE44787).

Analysis of sequencing data
Low quality reads were produced by an Illumina 1G Genome Analyzer (Illumina, San Diego, CA, USA). Clean reads were obtained after removing adaptor/acceptor sequences, contaminants formed by adaptor-adaptor ligation, and low-quality tags. The clean reads were then mapped to the B73 genome (http:// www.maizesequence.org, release 5a.59) using SOAP, and sequences derived from rRNAs, tRNAs, and snoRNAs annotated in the National Center for Biotechnology Information (NCBI) and Rfam database were eliminated [69]. The remaining small RNA sequences were considered as candidates for identification of conserved and novel miRNAs.

Identification of conserved and novel miRNAs
To identify conserved miRNAs, all candidate sequences were Blastn-searched against miRBase 18.0. Sequences perfectly matched to mature and precursor sequences in the database were identified as conserved miRNAs.
For novel miRNA identification, we used the Einverted program in the Emboss software package to located inverted repeats (stem-loop structures) [26] with the following parameters: threshold = 30, match score = 3, mismatch score = 3, gap penalty = 6, and maximum repeat length = 240 [27]. Each inverted repeat was extended 10 nt on either side, and its secondary structure was predicted by RNAfold [28]. Unique reads in inverted repeats were evaluated by MirCheck using default parameters [27]. miRNA sequences passing MirCheck were then inspected manually. Of the remaining sequences, only those with more than five read counts in at least one library and minimal folding free energy indexes (MFEI) higher than 0.85 were selected as novel miRNAs (MFEI = [(MFE/RNA sequence length) 6 100]/(G+C)%) [29].

Stem-loop real-time PCR
Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and treated by RNase-free DNase I treatment (Promega, Madison, WI, USA) to eliminate genomic DNA contamination. Then, 2 mg total RNA was used for first strand cDNA synthesis with a One Step PrimeScript miRNA cDNA synthesis kit (Takara, Dalian, China). Real-time RT-PCR was carried out with 0.5 mL cDNA on CFX96 real-time system (Bio-Rad, Hercules, CA, USA) using a SYBR Premix Ex Taq kit (Takara, Dalian, China). Protocol of the stem-loop real-time PCR was as follows: initial denaturation at 95uC for 10 s, followed by 40 cycles at 95uC for 10 s, 60uC for 30 s. At the end of the PCR reaction, a melting curve was obtained by holding at 95uC for 5 s, cooling to 60uC for 5 s, and then heating slowly at 0.5uC/s until 95uC. Only the results of miRNA amplification that produced one-peak melting curve at the correct annealing temperature were used for next analysis. Simultaneously, the threshold was automatically set and the threshold cycle (Ct value) was recorded. Expression levels of small nuclear RNA U6 were used as a reference for normalization [70]. The DDCt method was used to determine relative expression levels [71]. Four replicates were performed for each experiment. Primers used in real-time PCR analysis are listed in Table S5.

Northern blot analysis
Northern blot hybridization was performed using the MiRNA Northern Blot Assay Kit (Signosis, Sunnyvale, CA, USA) according to manufacturer's protocol. For each sample, 20 mg total RNA was mixed with RNA loading buffer and load onto one well of 15% polyacrylamide gel. The gel was then run at 60 V until bromophenol blue reaches approximately 3 cm away from the bottom of the gel. RNA was transferred to membrane at 60 V at for l hour using a BioRad Trans-Blot Cell (Bio-Rad, Hercules, CA, USA). Biotin labeled miRNA probe (Signosis, Sunnyvale, CA, USA) was used for hybridization. U6 bands were shown as a loading control.

Target gene prediction
Potential target genes of miRNAs were predicted using Patscan as previously described [24] with default parameters. Targets were selected according to the following criteria: 1) less than three mismatches between miRNA and target; 2) no insertions and deletions; and 3) no mismatches in positions 10 and 11 of mature miRNAs.

Validation of predicted target genes by real-time RT-PCR
Total RNA from MP, GP, MS and PS was isolated as described above, and 2 mg total RNA was reverse-transcribed using oligo (dT) primer and M-MLV reverse transcriptase (Promega, Madison, WI, USA) according to manufacturer's protocol. The real-time PCR was performed in triplicate on a CFX96 real-time system (Bio-Rad, Hercules, CA, USA) using SYBR Green Real-time PCR Master Mix (Toyobo, Osaka, Japan). Briefly, 1 mL cDNA was amplified with the following protocol: 94uC for 5 min, followed by 41 cycles of 94uC for 15 s, 65uC for 30 s. Finally, the melting curve was adjusted as 94uC for 5 s and 60uC for 1 min, and then rising to 94uC by increment of 0.5uC/s. 18S rRNA was used for each sample as an internal control. Relative expression levels of the genes were calculated using the DDCt method [71]. All primer sequences were listed in Table S5.

Identification of differentially expressed miRNAs
Transcripts per million reads (TPM) were used to represent normalized miRNA expression levels. TPM of miRNAs in the four libraries were calculated as: actual miRNA count/total count of clean reads 6 1,000,000. For the identification of miRNAs differentially expressed between mature and germinated pollen, we used fold-change = log 2 6 (GP/MP). Only miRNAs with foldchanges greater than 1.5 and P-values less than 0.01 were selected (Chi 2 test) [72]. For the identification of miRNAs differentially expressed before and after pollination in silks, we used fold-change = log 2 6(PS/MS). miRNAs with expression level increases 1.5-fold higher in PS, and read counts more than 10-fold greater in PS or MS than in GP and MP were considered to be up-regulated in PS. miRNAs with expression level decreases 2-fold greater in PS compared with MS were considered to be down-regulated in PS.

Gene Ontology (GO) analysis
We performed GO analysis on target genes of miRNAs involved in pollen-silk interactions using Singular enrichment analysis (SEA) [63] (http://bioinfo.cau.edu.cn/agriGO/analysis.php). Overrepresented gene terms in biological processes, molecular functions, and cellular components categories were selected using the thresholds of P-value ,0.001 and false discovery rate (FDR) ,0.05.

Author Contributions
Conceived and designed the experiments: XSZ. Performed the experiments: XML XYZ. Analyzed the data: YLS XML. Wrote the paper: YLS XSZ.