Homozygous broccoli seeds named TNK-002 (stored in Tianjin Kernel Vegetable Research Institute, Tianjin, China) were planted in soil (Klasmann Deilmann GmbH, Germany) under controlled conditions with a 16/8 h light/dark cycle at 25° and 22°C, respectively. The 10-old day seedlings were transplanted in a greenhouse under growth condition as mentioned above. At approximately 70 days, the curds formed, which are mostly composed of flower buds with different sizes. The flower buds were isolated from the curds, and their sizes were measured by vernier caliper. The pollen in different developmental stages such as the uninucleate microspores, binucleate and trinucleate pollen grains, from differently sized flower buds were stained with DAPI fluorescent dye, and detected under a fluorescence microscope (Olympus BX53, Japan). Using the information obtained about the relationship between bud size and pollen developmental stage, three flower bud size ranges (3.0, 4.5, and 6.0 mm in length) were selected that would reproducibly yield uninucleate, binucleate pollen grains and trinucleate pollen grains, respectively.

Three groups of different-sized flower buds with uninucleate microspores, binucleate pollen grains and trinucleate pollen grains were collected. At least 50 flower buds per group were used. Then, stamens were isolated from the flower buds and deposited in 10-mL centrifuge tubes that contained 5 mL modified B5 culture medium. The stamens were then homogenized for 10 min to release single microspores, binucleate pollen grains or trinucleate pollen grains. In this process, the stamens were gently ground to avoid releasing single somatic cells. Subsequently, the mixtures that contained uninucleate microspores, binucleate pollen grains or trinucleate pollen grains were filtered with a 50-μm filter membrane to remove large impurities, and then filtered with a 40-μm filter membrane. The filtrate was transferred to other 10-mL centrifuge tubes and centrifuged at 1,000 g/min for 10 min. The remaining solid materials were resuspended with 5 mL modified B5 culture medium and centrifuged at 1,000 g/min for 10 min. Finally, the rinsed solid materials were collected, immediately frozen in liquid nitrogen, and stored at −80°C. To assess the purity of the isolated samples, approximately 50 μL of resuspended cells were stained with 5 μL DAPI. The stained cells were observed by fluorescence microscope (Olympus BX53, Japan) under UV light. In each sample, cells in at least 10 independent visual fields under the microscope were counted and analyzed based on their cytological characteristics.

Total RNAs were isolated using TRIzol reagent (Invitrogen, USA). RNA samples with high purity (OD₂₆₀/₂₈₀ = 1.8–2.2) and high integrity (RNA integrity number, RIN > 8.0) were used to construct small RNA libraries by TruSeq Small RNA Library Prep Kit (Illumina, USA) according to the manufacturer's protocols. In brief, total RNAs were subjected to 15% denaturing polyacrylamide gel electrophoresis. Small RNA fractions of 18–30 nt were isolated from the gel and purified. The 5′ and 3′ adapters were ligated to the isolated small RNAs by T₄ RNA ligase (TaKaRa, Japan), and then converted to cDNAs which were further used as samples to conduct RT-PCR amplification. The adapter sequences and RT-PCR primers were showed in Table S1. The PCR products were purified and subjected to deep sequencing by HiSeq™ 2000 (Illumina, USA) [Beijing Genomics Institute (BGI), China]. Three repeats were conducted in the sequencing process. 49-nt raw reads were produced by the sequencing system.

Raw reads of small RNAs were obtained via Solexa sequencing. Low-quality reads and contaminants in the raw reads, for example, reads that were less than 18 nt in length, adaptor-only and polyA reads, were discarded. The adaptor sequence of each read was trimmed. Filtered high-quality reads that were 18–30 nt in length were designated as clean reads and further assembled into unique reads. The clean reads and unique reads were used in the following analysis. The length distributions of clean reads in uninucleate microspores, binucleate and trinucleate pollen grains were analyzed. Pairwise comparisons of common/specific clean and unique reads were also conducted. Then, all unique reads were aligned with structural non-coding RNAs (rRNA, tRNA, snRNA, snoRNA, etc.) deposited in the GenBank and Rfam databases. Reads that matched with the sequences in these two databases were excluded in the further analysis. Subsequently, the unique reads were subjected to Blastn analysis against plant miRNAs, particularly Arabidopsis and other Brassicaceae plant miRNAs deposited in miRBase 18.0 (http://www.mirbase.org/) to identify the conserved miRNAs. The unique reads detected in each small RNA library were annotated in accordance with the criteria: ribosomal RNA > known miRNA > repeat > extron > intron.

The unannotated reads and reads from the introns were aligned to broccoli EST data (Accession number: PRJNA361430) to identify potentially novel miRNAs. Sequences >100 bp surrounding the matched region were extracted, and utilized to predict pre-miRNA candidates using the Mireap program (https://sourceforge.net/projects/mireap/). Parameters were as follows: Minimal miRNA sequence length (18 nt); Maximal miRNA sequence length (25 nt); Minimal miRNA reference sequence length (20 nt); Maximal miRNA reference sequence length (23 nt); Maximal copy number of miRNAs on reference (20 nt); Maximal free energy allowed for a miRNA precursor (−18 kcal/mol); Maximal space between miRNA and miRNA^* (300 nt); Minimal base pairs of miRNA and miRNA^* (16 nt); Maximal bulge of miRNA and miRNA^* (4 nt); Maximal asymmetry of miRNA/miRNA^* duplex (4 nt); Flank sequence length of miRNA precursor (20 nt). The predicted pre-miRNA sequences were further assessed by M-fold (Zuker, 2003), and only structures with the lowest free energies were selected. In addition, to separate novel miRNA candidates from possible siRNAs, small RNA read distribution was analyzed using Blastall and Omega 2.0 softwares. Small RNAs with wide distribution on the precursor sequences and having reads that almost equally map to both plus and minus strands were excluded.

In small RNA deep sequencing, the count of clean reads originating from each miRNA represents the expression abundance or level of the corresponding miRNA. At least 16-nt overlap was required to confirm a read that generate from a certain miRNA. To explore the expression patterns of miRNAs in three different developmental phases of broccoli pollen, the frequency of each miRNA was normalized to the same order of magnitude according to the formula: Normalized expression = actual miRNA count/total count of clean reads × 1,000,000. The normalized read count in miRNAs from which no reads were detected, was set at 0.01. The differential expression level of each miRNA in two arbitrary developmental phases was evaluated by the fold change of the normalized expression level. The log base 2-fold change was calculated as follows: the log base 2-fold_change = log₂ (A/ B) (A and B represent the normalized expression level of miRNA in any two developmental phases). Then, statistical analysis was performed according to Poisson distribution. The P-value was calculated and further corrected by Bonferroni Correction. For the identification of significantly expressed miRNAs, the criteria was used as if (log base 2-fold change ≥1 or ≤–1) and (corrected P < 0.05). A hierarchical cluster analysis of all known miRNAs was done by the package “gplots” of the R project according to the log base 2-fold_change of the normalized expression levels in two arbitrary developmental phases (http://www.r-project.org/).

Conserved and novel miRNAs were aligned with the broccoli ESTs to predict the potential targets of these miRNAs. If one EST sequence satisfied the parameters that were suggested by Allen et al. (2005) and Schwab et al. (2005), the EST was considered as a putative target of a specific miRNA. In addition, the Gene Ontology (GO) analysis of putative targets was conducted. Target gene candidates from broccoli uninucleate microspores, binucleate pollen grains and trinucleate pollen grains were mapped to GO terms in the database (http://www.geneontology.org/), respectively. Then, the gene numbers per term were calculated. The hypergeometric test was conducted to identify the GO terms in target gene candidates that were significantly enriched compared with the reference gene background (broccoli ESTs). The calculating formula is:

In the formula, N is the number of all genes with GO annotation; n is the number of target gene candidates in N; M is the number of all genes that are annotated to a certain GO term; m is the number of target gene candidates in M. Then, the Bonferroni Correction is used for the P-value to obtain a corrected P-value. GO terms with corrected P < 0.05 are defined as significantly enriched in target gene candidates. To further describe the enrichment level of each GO term, the formula [Enrichment = (n/N)/(m/M)] is used.

Total RNA was isolated from microspores, binucleate and trinucleate pollen grains by TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. RNAs free of contaminated genomic DNAs were reverse-transcribed using miRNA-specific stem-loop primers (Table S1). The subsequent PCR amplification and re-sequencing analysis of miRNAs were conducted as described by Geng et al. (2014).

The expression levels of miRNAs and miRNA targets were further validated by real-time quantitative RT-PCR (qRT-PCR). Specific primers were designed based on the sequences of miRNAs and corresponding targets (Table S1), respectively. Faststart Universal SYBR Green Master (Roche, Germany) was used in all qRT-PCR experiments. Small nuclear RNA U6 was used as an internal reference in miRNA expression level analysis, while the actin gene of broccoli was used as an internal control in target expression level analysis. The relative expression levels of miRNAs and targets were calculated by the comparative 2^−ΔΔCT method according to the manufacturer's recommendations. To ensure the reliability of quantitative analysis, three batches of independently isolated RNAs from the three different developmental phases of broccoli pollen were used, and three technological replicates were performed.

To validate miRNA targets, a modified 5′ RLM-RACE (RNA ligase-mediated rapid amplification of cDNA ends) was conducted using a 5′-Full RACE Kit (TaKaRa, Japan) in accordance with the manufacturer's instructions with some modifications. Total RNA that was isolated from uninucleate microspores, binucleate and trinucleate pollen grains were ligated to the 5′ RNA adaptors by T₄ RNA ligase (TaKaRa, Japan). A reverse transcription reaction was performed with 9-nt random primers. Subsequently, a nested PCR reaction was performed. In the first PCR reaction, the reverse transcription product was amplified using the 5′ RNA adaptor outer primers and gene-specific outer primers (Table S1). In the second PCR reaction, the first PCR products (1 μL) were amplified as templates with an inter primer combination (5′ RNA adaptor inter primers and gene-specific inter primers) (Table S1). After amplification, 5′ RACE products were gel-purified and cloned using the T-A cloning method by Mighty TA-cloning Kit (TaKaRa, Japan). A minimum of 10 independent clones from each PCR reaction were randomly selected and sequenced.

Microspores, binucleate pollen grains and trinucleate pollen grains, were isolated from broccoli flower buds (Figure 1). Statistical analysis indicated that the purity of the isolated uninucleate microspores and unmatured pollen grains in binucleate and trinucleate phases exceeded 95% (Table S2). This result indicated that the quality of the isolated samples was sufficient to conduct the following analysis.

Three small RNA libraries from the corresponding isolated samples were constructed and sequenced by Solexa technology (Accession number: PRJNA361414). In total, 21,659,809, 20,439,333, and 21,429,804 raw reads generated from the RNA of uninucleate microspores, binucleate pollen grains and trinucleate pollen grains, respectively. After removing low-quality reads, contaminants and adaptors, 21,429,050 clean reads representing 7,160,154 unique sequences in uninucleate microspores, 20,173,912 clean reads representing 6,299,777 unique sequences in binucleate pollen grains and 21,209,081 clean reads representing 6,952,258 unique sequences in trinucleate pollen grains were obtained (Table S3). Pairwise comparisons of the small RNAs from the three small RNA libraries were conducted. Consisted with reported investigations, the vast majority of unique reads are non-overlapping. In uninucleate microspores and binucleate pollen grains, only 14.79% unique reads were shared in the two developmental phases. In uninucleate microspores and trinucleate pollen grains, 43.28 and 41.59% unique reads were specific in uninucleate and trinucleate phase, respectively. Similarly, the majority of the unique reads, which were detected in binucleate and trinucleate pollen grains, were specific to the binucleate (39.95%) pollen grains or trinucleate (45.59%) pollen grains. Only a small proportion of these reads (14.46%) was shared in these two developmental phases (Figure S1). All unique reads were then further annotated. Only a few reads were annotated as ribosomal RNAs, tRNAs, snRNAs, snoRNAs, miRNAs, and other previously known RNAs. More than 97% of the reads were unannotated. The length distributions of all of the clean reads were also explored. The majority of the reads were 18–24 nt in length. The highest abundance was found for reads that were 24 nt in length followed by 21 nt, 22 and 23 nt in length (Figure S2).

In total, 333 known miRNAs were detected in the three early stages of broccoli pollen development. These miRNAs originated from 235 different miRNA families. The majority of these conserved miRNAs (325/333) were simultaneously detected in three samples. However, the expression abundance of each miRNA varied (Tables S4, S5). The members in the miRNA families, such as miR167, miR166, miR156/157, miR397, miR165, miR158, miR168, miR2199, miR2911, miR5782, and miR6300 families, had a high expression abundance (clean read count of miRNA > 10,000) in all three developmental phases. MiRNAs, such as miR159, miR162, miR172, miR390, miR2916, and miR391, showed a moderate expression abundance (clean read count of miRNA > 1,000), whereas most miRNAs, especially, miR169i, miR169j, miR4221, miR399d, miR5026, miR5642a, miR5642b, miR775, and miR828 displayed a low expression abundance (clean read count of miRNA < 1,000) (Figure S3, Table S4, S5). Further analysis indicated that only one gene was detected in most miRNA families. At least two genes were detected in the other 25 miRNA families. The miR156/157, miR166, miR167, miR169, miR172, miR395, and miR399 families had more than four members, although high or moderate expression abundance of several these miRNA families was not detected (Figure S4).

A total of 55 novel miRNA candidates were observed. Among these novel miRNAs, 31, 28, and 26 predicted novel miRNAs were detected in the uninucleate microspores, binucleate pollen grains and trinucleate pollen grains, respectively. The lengths of these newly predicted miRNAs ranged from 20 to 23 nt, among which 21-nt miRNAs were dominant. The negative minimum folding free energy (MFE) of these miRNA precursors ranged from −20 kcal/mol to −85.3 kcal/mol. The average free energies were approximately −40.37 kcal/mol (Table 1). Each novel miRNA precursor can form regular stem-loop structures. The star sequences (miRNA^*) of 28 novel miRNA candidates were observed, which further confirmed the existence of these miRNAs. Further analysis indicated that the expression patterns of most of the predicted novel miRNA varied. Forty-one of the 55 predicted novel miRNAs were specifically expressed in the uninucleate microspores, binucleate pollen grains or trinucleate pollen grains. Among them, 15 predicted novel miRNAs were only detected in the uninucleate microspores, 13 predicted novel miRNAs were specifically expressed in the binucleate pollen grains. Similarly, 13 predicted novel miRNAs were only detected in the trinucleate pollen grains (Table S6). Four predicted novel miRNAs (bol-miR01, bol-miR03, bol-miR24, and bol-miR25) were detected in the uninucleate microspores and binucleate pollen grains but not in the trinucleate pollen grains. Bol-miR27, bol-miR30, bol-miR33, and bol-miR38 were only detected in the binucleate and trinucleate pollen grains. However, no predicted novel miRNAs were specifically expressed only in uninucleate microspores and trinucleate pollen grains. Nevertheless, six predicted novel miRNAs (bol-miR02, bol-miR04, bol-miR7, bol-miR9, bol-miR18, and bol-miR22) were simultaneously detected in the three developmental phases (Table S6). In total, 49 of the 55 predicted novel miRNAs displayed specific expression patterns in the three different developmental phases of broccoli pollen.

miRNAs with similar expression patterns may have similar functions. Thus, the expression patterns of 333 conserved known miRNAs were further analyzed. According to the hierarchical clustering data, these 333 miRNAs could be roughly classified into 11 groups (Figure 2 and Table S7). miRNAs from the same family, such as the members of miR156/157, miR160, miR165, miR167, miR319, miR3955, and miR858, always belonged to the same group, implying that they exhibited similar expression patterns in the development of broccoli pollen.

Based on the expression levels and patterns of each miRNA, sixty of the 333 known conserved miRNAs were confirmed to exhibit significantly differential expression levels in the three different developmental phases of broccoli pollen (Figure 3 and Table S8). Pairwise comparison analysis indicated that 19 known miRNAs in the uninucleate microspores and binucleate pollen grains showed significantly differential expression levels. Among these known miRNAs, 7 miRNAs displayed higher expression levels, whereas the remaining other 12 known miRNAs displayed lower expression levels in the uninucleate microspores than those in the binucleate pollen grains. Compared with those in the trinucleate pollen grains, 35 known miRNAs with significantly differential expression levels were detected in the uninucleate microspores. The expression levels of miRNAs in the binucleate and trinucleate pollen grains were also compared. Twenty-six known miRNAs with significantly differential expression levels were detected in the two phases. Of these miRNAs, 12 had up-regulated expression in the binucleate phase and 14 showed down-regulated expression (Figure 3, Tables S8, S9). Further analysis indicated that the expression levels of miR156h, miR391, miR3954, miR472, miR5665, and miR5716 in the uninucleate microspores significantly differed from those in the binucleate and trinucleate pollen grains. MiR156 h and miR5716 showed considerably low expression levels in the uninucleate microspores, while four other miRNAs showed particularly high expression levels. In binucleate pollen grains, miR5767 showed particularly low expression levels, whereas miR172, miR827a, and miR862b showed specifically high expression levels. In the trinucleate pollen grains, the expression levels of 11 known miRNAs (miR164c, miR169b, miR169c, miR172c, miR172d, miR858a, miR164a, miR6034, miR6108f, miR858b, and miR3434), were significantly different from those in uninucleate microspores and binucleate pollen grains. MiR169b, miR169c, miR858a, miR6034, and miR858b showed considerably high expression levels in the trinucleate pollen grains, while the six other known miRNAs demonstrated evidently low expression levels (Figure 4, Table S8).

The expression trends of miRNAs with significantly differential expression levels were also analyzed. Nine of the 60 known miRNAs increased their expression levels with the development of pollen, while the expression levels of 22 miRNAs decreased. Approximately half of these miRNAs showed “V” or reversed “V” expression patterns. In addition, 22 of the 60 differentially expressed miRNAs belonged to seven miRNA families (miR164, miR165, miR169, miR172, miR2111, miR159/319, and miR858). Similarly, miRNAs from the same miRNA family exhibited similar expression trends (Figures 2, 4 and Table S7).

In 55 predicted novel miRNAs, the expression profiles of predicted novel miRNAs with read counts >20 were further observed. The results indicated that, except for seven novel miRNAs only detected in a specific developmental phase, four predicted novel miRNAs showed significantly differential expression levels in the three samples (Figure S5). Bol-miR01 showed similar expression levels in the uninucleate microspores and binucleate pollen grains but was not detected in the trinucleate pollen grains. Bol-miR03 was specifically and highly expressed in the binucleate pollen grains but was lowly expressed in both uninucleate microspores and trinucleate pollen grains. Bol-miR33 and bol-miR38 were undetected in the uninucleate microspores, while both were highly expressed in the binucleate and trinucleate pollen grains. The other four novel miRNAs all detected in three developmental phases did not exhibit significantly differential expression levels (Figure S5).

A total of 552 target sites were identified for 127 conserved miRNAs. More than 10 putative genes were targeted by miR156, miR157, miR414, miR4993, miR5021, miR5137, miR5139, miR5215, miR5227, miR5631, miR5782, miR6443, and miR854, whereas only a few targets were detected for most miRNAs (Table S10). Homologous analysis and functional annotation of the putative targets were performed. The annotated targets were mainly associated with the regulation of plant growth and development (Table S10). All of the targets were subjected to GO analysis. The GO terms of the conserved miRNA targets from the uninucleate microspores and the two following developmental products of microspores were significantly enriched in 23 biological processes, 11 cellular components and 15 molecular functions. In biological processes, 17 GO terms, including GO:0003156 (regulation of organ formation) and GO:2000027 (regulation of organ morphogenesis), were enriched for the targets predicted in all three developmental phases. GO:0009943 (adaxial/abaxial axis specification), GO:0009955 (adaxial/abaxial pattern specification) and GO:0050793 (regulation of developmental process) were significantly enriched for the targets predicted in the uninucleate microspores and binucleate pollen grains. Three other GO terms were only significantly enriched for the targets predicted only in the trinucleate pollen grains. In cellular components, nine of the 11 GO terms were significantly enriched for the targets predicted in the three developmental phases. In molecular function, eight of the 15 GO terms were especially enriched for the targets predicted in the trinucleate pollen grains. Other GO terms from this section were specifically enriched for the targets predicted in the uninucleate microspores and binucleate pollen grains (Figure 5 and Table S11).

A total of 69 putative targets were predicted to be cleaved by 40 predicted novel miRNAs (Table S12). The annotated targets also underwent GO analysis. No significantly enriched GO terms were detected for the targets of novel miRNAs found in the uninucleate microspores. The GO terms associated with floral organ development, pollen development, postembryonic organ development and response to nitrate only significantly enriched for the targets predicted in the binucleate pollen grains. For the predicted miRNA targets in the trinucleate pollen grains, 10 significantly enriched GO terms mainly involved in biological processes were detected. The hormone binding GO term (GO:0042562) exhibited significant enrichment for the targets predicted in both the binucleate and trinucleate pollen grains (Figure S6).

Twelve conserved miRNAs and 25 of the 55 predicted novel miRNAs were randomly selected for stem-loop RT-PCR and re-sequencing validation. Most of these detected miRNAs could be amplified except for bol-miR10 and bol-miR21. Nevertheless, the precursor sequences of the two novel miRNAs were successfully amplified. Sequencing analysis confirmed that the mature miRNA and precursor sequences of bol-miR10 and bol-miR21 were consistent with those predicted by the bioinformatics method (Table 1). The expression patterns of six conserved miRNAs (miR165a, miR169b, miR159a, miR319c, miR391, and miR858a) and eight novel miRNAs (bol-miR04, bol-miR06, bol-miR17, bol-miR21, bol-miR23, bol-miR35, bol-miR47, and bol-miR51) were further explored by real-time quantitative stem-loop RT-PCR. The expression patterns of six conserved miRNAs and six of the eight novel miRNAs that were detected by qRT-PCR were highly consistent with the data produced by high-throughput sequencing (the coefficient of determination R² > 0.86) (Figure 6 and Figure S7). In addition, the expression patterns of the 12 putative targets for bol-miR04, bol-miR17, bol-miR23, bol-miR34, and bol-miR51 were observed (Table S12). The expression trends of at least one putative target for each novel miRNA were confirmed to be negatively correlated with the expression levels of the corresponding miRNA (Figure S8).

The modified 5′ RLM-RACE method was conducted to observe the cleavage sites of 12 putative targets of several novel miRNAs. Consistent with the expression trend analysis, the targets, whose expression trend negatively correlated with their miRNAs, were successfully cleaved by the corresponding miRNAs mainly in the complementary region of the miRNA and mRNA sequences. The cleavage sites were predominant in the 10^th and 11^th nucleotides from the 5′ end of miRNAs (Figure 7).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.