337S phenotype of different growth conditions from different sowing dates
The SL and NN plants of 337S were planting at different plots of the same experimental field (Additional file 2: Figure S1a and Additional file 3: Figure S2a). Significant difference among the 337S, HZ09 and Huamai 2566 plants at heading and flowering in SL early-sowing time was observed (Additional file 2: Figure S1b-d). 337S plants flowered about twenty days later than HZ09 and Huamai 2566 plants. When the 337S plants were at early booting stage, HZ09 and Huamai 2566 plants were at heading stage in SL condition (Additional file 2: Figure S1e). However, there was no significant difference on booting stage, heading date and flowering period among the 337S, HZ09 and Huamai 2566 plants in NN sowing-time (Additional file 3: Figure S2b-e).
Our previous study indicated that degeneration of microspore mother cells in meiosis I was the main reason resulting in 337S line sterile at SL condition [27]. The growth of anthers prior to the microspore stage was divided into 3 phases: early anther development stage (stage 1, from spikes differentiation to formation of sporogenous cell, spike length < 3.5 cm), microspore mother cell stage (stage 2, 3.5 cm ≤ spike length < 4.5 cm) and meiosis stage (stage 3, 4.5 cm ≤ spike length ≤ 5.5 cm) according to the length of the wheat spikes planted in the field. Morphologically, the vegetative development of spikes appears no obvious difference at different stages between SL and NN conditions (Fig. 1).
Compared with NN plants of 337S, spikelets of SL plants were abnormal dehiscence (Fig. 2a-d). The glumes and lemmas of SL spikelets had a bigger angle than that of NN plants at full-bloom stage (Fig. 2e, f). The anthers of SL were shorter than those of NN plants, and SL anthers showed a slight yellow discolouration at mononuclear pollen stage and late binuclear pollen stage, while the NN anthers remained green (Fig. 2g, h). However, there was no significant difference of the pistil filaments between NN and SL plants (Fig. 2h). Unlike the normally grain filling on NN spikes, no visibly filled grain on SL plants at the filling stage was observed (Fig. 2i). To determine whether the SL plants could produce abnormal anthers leading to male sterility as reported previously [27], spike lengths about 3.9–4 cm, 4.5–4.6 cm and 5–5.1 cm corresponding to microspore mother cell stage (NN1/SL1) and meiosis stage (NN2/SL2, NN3/SL3) were collected from NN and SL plants for performing semi-thin sections, respectively (Fig. 1). Anthers transverse cross-section observation indicated that the anthers of SL1 plants developed normally until the microspore mother cell (MMC) stage (Fig. 2j, k). Histological results showed that the obvious aberration of SL2 anther occurred during meiosis, around the prophase stage, in which meiotic cells degenerated and collapsed rapidly, whereas the meiotic cells of NN2 anthers remained normally development (Fig. 2l, m). At the stage when the meiotic cells in NN3 anther produced dyads, the SL3 plant appeared similar phenotype to itself at prophase stage without any intact meiotic cells in anther chambers (Fig. 2n, o). These observations indicated that the first detectable sign of male sterility occurred at the meiotic prophase (MP) stage with degenerated meiocyte cell in SL2.
Overview of the small RNA sequencing
In order to reveal the regulatory mechanism of SL abortion and explore new regulators for the anther development in wheat, small RNA and degradome sequencing was applied to identify critical genes leading to male sterility for 337S line. A total of 60,119,469 raw reads were obtained from four libraries; Of which 14,947,750 were generated from NN1 library, 14,290,542 from NN2, 14,998,252 from SL1, and 15,882,925 from SL2 library. After filter out reads of N% > 10%, reads with polyA/T/G/C, low-quality and adapter sequences, 14,597,462 (97.66%), 13,981,982 (97.84%), 14,635,455 (97.58%) and 15,584,410 (98.12%) clean reads were obtained from NN1, NN2, SL1 and SL2 libraries, respectively (Additional file 4: Table S2). The lengths of the majority of the sRNAs were 21~ 25 nt, in which the 24 nt sequences showed the highest enrichment among them, accounting for about half of the total reads, followed by 21 nt. A similar size distribution was observed for the length of 22 nt, 23 nt and 25 nt sRNAs in each library. Moreover, the frequency of 21 nt sRNAs in SL plants was consistent with the scale in NN plants at the same stage, such as 15.01 and 14.8% for NN1 and SL1, 12.04 and 12.55% for NN2 and SL2, respectively. However, the libraries of SL had higher abundance (49.6 and 51.02%) in 24 nt sequences than that of NN plants (44.38 and 46.06%) at both stages (Fig. 3a). sRNAs were annotated and grouped into several classes by aligning with known non-coding RNAs in the Rfam and NCBI database. Overall, the percentage of sRNA reads that matched for a specific group was similar among the four libraries. Nevertheless, interesting data was found from the unannotated reads between NN and SL plants. SL1 samples showed lower accumulation of unique unannotated reads than NN1 samples (54 and 58% of mapped sRNA for SL1 and NN1, respectively), likewise between SL2 and NN2 samples (Fig. 3c-f).
Differentially expressed miRNAs between the SL and NN plants
The sRNAs with classic miRNA secondary structure for a length of at least 18 nt were aligned to all known plant miRNA sequences in miRBase21.0, 94 annotated known miRNAs were identified from the NN and SL libraries (Additional file 5: Table S3); Of which, 83 were from NN1 library, 84 from NN2, 89 from SL1, and 87 from SL2 library. Analysis of nucleotide bias at each position of miRNAs showed that the first nucleotide tended to be uracil (U) for each sample (Additional file 6: Figure S3). As shown in Fig. 3b, 62 known miRNAs were assigned to 35 known families, in which three fifths of the families contained a single member. Ten miRNA families were represented by only two members. Both MiR9657 and MiR9666 family contained 3 members. The two large numbers of these miRNA families were seven and eight, represented by MiR1120 and MiR1122 family, respectively (Fig. 3b). Prediction of novel miRNAs from the 4 sRNA libraries discovered the precursors for eight of the putative novel miRNAs, which were designated as Novel_01-Novel_08. All novel miRNAs were detected in the each of the four libraries. The precursor length of the novel miRNAs ranged from 57 to 292 nt, and were processed into 20 nt to 24 nt miRNAs (Additional file 7: Table S4). The expression level of each miRNA was normalized by TPM. More than 76% of TPMs for each of the four libraries were more than 60 (Additional file 8: Table S5).
Comparative analysis of 102 unique miRNAs (94 known and 8 novel miRNAs) showed that 88 and 89 miRNAs overlapped between MMC stage and MP stage, respectively (Fig. 4a, b). Compared to NN samples, 9 and 6 specific miRNAs were detected in the SL samples at stages of MMC and MP, respectively (Fig. 4a, b). Two microRNAs, tae-miR399 and tae-miR9778, were specifically expressed in SL1 (Fig. 4c). One microRNA, tae-miR9674a-5p was only identified in NN1. However, these three sample-specific miRNAs were expressed at low levels (Additional file 9: Table S6). In addition, 83 miRNAs were common in all NN and SL samples (Fig. 4c). The heatmap showed 53 miRNAs differentially expressed among the four samples, including 52.1% (49 of 94) conserved miRNAs and 50% (4 of 8) novel miRNAs (Fig. 4d). 27 miRNAs were significantly differentially expressed in the sample of SL1 in comparison with the NN1 (control) at MMC stage, and 25 miRNAs were significantly differentially expressed between the SL2 and the NN2 (control) at MP stage (Additional file 10: Figure S4, Additional file 11: Table S7 and Additional file 12: Table S8). There were 10 common miRNAs across these two comparisons. Thus, total 42 significantly differentially expressed miRNAs were found in both MMC and MP stages together (SL1 vs NN1 and SL2 vs NN2) (Additional file 10: Figure S4), 60 miRNAs showed no differential expression at both MMC and MP stages (Additional file 13: Table S9–1). The 42 differentially expressed miRNAs were divided into 5 categories based on expression pattern. Only tae-miR1122a was upregulated in the both two comparisons (Additional file 13: Table S9–2), while tae-miR1122c-3p and tae-miR5200 were downregulated in these two comparisons (Additional file 13: Table S9–3). Seven miRNAs showed opposite expression patterns between the two comparisons (Additional file 13: Table S9–4). 17 miRNAs exhibited significantly differential expression at MMC stage (Additional file 13: Table S9–5), 15 miRNAs showed significantly differentially at MP stage (Additional file 13: Table S9–6).
qRT-PCR was performed to validate the sequencing data. Based on the data analysis for each category above, the miRNAs which had been confirmed to be conservative in reproductive development of plant and two novel miRNAs were selected for qRT-PCR analysis. The results showed that the expression patterns of qRT-PCR for these differentially expressed miRNAs were consistent with the expression trend of RNA-Seq data, indicating high reliability of the sequencing (Fig. 5 and Additional file 14: Figure S5). According to previous studies, miR156, miR159, miR160, miR164, miR167, miR319, miR396 and miR5200 were mainly involved in floral development [14, 28]. No obvious expression changes of tae-miR159a, tae-miR160, tae-miR164, tae-miR167a and tae-miR167c were found between SL1 and NN1, and between SL2 and NN2 samples. tae-miR156 and tae-miR319 did not show two-fold change as RNA-seq data did. Sequencing data indicated the expression levels of tae-miR2275-3p, miR396-5p and tae-miR5200 were significant difference between NN1 and SL1 plants, and downregulated in SL1 plants. tae-miR2275-3p and tae-miR5200 were also extremely suppressed in SL2 compared to NN2 plants (Fig. 5).
Identification of miRNA targets through degradome sequencing
From degradome sequencing, more than 1*10 million raw reads from each library were obtained. After removing the reads < 15 nt and 3′ adaptor, 9,829,031 and 9,421,684 transcript reads from SL library and NN library were mapped to wheat genome, respectively (Additional file 15: Table S10). The predicted targets were classified into five categories (0, 1, 2, 3 and 4). The targets with category 0 were evaluated as the most significant.
According to the expression characteristics of the wheat transcriptome data, t-plots were built for the high-efficiency analysis of the identified miRNA targets. As a result, a total of 643 transcripts targeted by 53 known miRNAs were obtained from transcriptome analysis. The tae-miR1127b-3p targets the highest number of 124 annotated and unknown transcripts. Nevertheless, no candidate target was validated for the novel miRNAs identified here. Further, total of 257 transcripts targeted by 20 miRNAs that were significant differentially expressed at MMC or MP stages were screened out from degradome data (Additional file 16: Table S11). Then, 12 miRNAs with TPMs > 25 in at least one sample in Table S9–5 and Table S9–6 together with tae-miR1122a in Table S9–2, tae-miR1122c-3p in Table S9–3 and tae-miR1122b-3p in Table S9–4 (marked by *) were validated by qRT-PCR. qRT-PCR showed that tae-miR1122a, tae-miR1122b-3p, tae-miR1122c-3p, tae-miR1124, tae-miR1127a, tae-miR9655-3p, tae-miR9668-5p, tae-miR9675-3p and tae-miR9679–5p were two-fold differentially expressed between SL1 and NN1, or between SL2 and NN2 (Fig. 6). However, the targets of tae-miR1124, tae-miR9655-3p and tae-miR9668-5p were unknown (Additional file 17: Table S12). Therefore, the 6 candidate miRNAs tae-miR1122a, tae-miR1122b-3p, tae-miR1122c-3p, tae-miR1127a, tae-miR9675-3p, tae-miR9679–5p together with tae-miR2275-3p were selected to further dissect their functions in wheat reproductive development. In addition, no target was identified for both tae-miR396-5p and tae-miR5200 from degradome data.
Our degradome data showed a new and remarkable counterpart gene of mammalian SWI/SNF chromatin-remodeling complex, which may be associated with meiotic DSB in 337S. SMARCA3L3 was the only target belonging to the cleavage sites classified for category = 0 (Fig. 7a) from all of the 14 genes targeted by tae-miR1127a (Additional file 16: Table S11). Noteworthy, the expression levels of SMARCA3L3 was negatively correlated with that of tae-miR1127a, which was significantly downregulated in SL1 compared with NN1 at MMC stage (Fig. 7c). The expression of XPB2 [29], a DNA repair helicase targeted by tae-miR1122c-3p, was highly upregulated in SL samples at both MMC and MP stages (Fig. 7c), implying that XPB2 as a DNA damage detection recognizer may be induced to repair DNA in sterile plants of 337S. EME1, another target gene of tae-miR1122c-3p, is required for generating meiotic crossovers by resolving the double Holliday junction in fission yeast [30], its expression pattern was also changed in SL1 (Additional file 18: Figure S6). Moreover, previous studies showed that the CAF1 encoding a conserved subunit of the CCR4-NOT complex was vital to meiotic progression [10]. Our degradome data revealed eight genes targeted by tae-miR2275-3p, three are CAF1 homologs in the category 0 (Fig. 7b). Interestingly, the CAF1 was dramatically suppressed in SL compared with NN at both MMC and MP stages (Fig. 7c). In addition, cyclin A3;2 (CYCA3;2) [31], the cell cycle regulator, targeted by tae-miR1122b-3p, showed significantly differential expression at MP stage (Fig. 7c). The expression levels of DEMETER (DME) with DNA glycosylase activity for removing 5mC [32] targeted by tae-miR9679–5p were increased in SL plants (Fig. 7c). Methylenetetrahydrofolate reductase (MTHFR) that is the methyl donor for numerous cellular reactions [33] targeted by tae-miR1122a, was also upregulated in SL compared to NN plants. Furthermore, 60S ribosomal protein L23 (RPL23) as the unique target gene of tae-miR9675-3p [34] and Defective in Anther Dehiscence1 (DAD1) [35], one of the targets of tae-miR1127b-3p, showed no obvious differential expression changes in SL plants (Additional file 19: Figure S7). The conserved miRNA-target interaction previously proved to be essential for floral development in plant, like tae-miR156-SPL17, tae-miR159a-GAMYB, tae-miR160-ARF18, tae-miR164-CUC2 and tae-miR167a-ARF12, showed very little changes in expression between SL and NN plants at MMC and MP stages (Additional file 19: Figure S7).
To further determine which miRNA-target interaction of tae-miR2275-CAF1, tae-miR1127a-SMARCA3L3, tae-miR1122c-XPB2, tae-miR1122c-EME1, tae-miR1122b-CYCA3;2, tae-miR9679-DME or tae-miR1122a-MTHFR is the vital and dominant determiner for regulating male sterility in 337S. Their expression patterns were analyzed in the stages prior to MMC. Spike lengths ranging from 2 to 2.1 cm and 3–3.1 cm corresponding to the early anther developmental stages were collected from normal day-length/normal temperature (NN-A/NN-B) and short day-length/low temperature (SL-A/SL-B) environment conditions (Fig. 1). Only one pair tae-miR1127a-SMARCA3L3 showed significant difference in both microRNAs and targets between SL-B and NN-B at early anther developmental stage. In addition, only the expression levels of SMARCA3L3 and CAF1 were negatively correlated with expression levels of tae-miR1127a and tae-miR2275-3p between SL-B and NN-B plants, respectively (Fig. 8). Therefore, the interactions of tae-miR1127a-SMARCA3L3 and tae-miR2275-CAF1 might be involved in regulating the male reproductive development in the 337S.