Skip to main content

Identification of miRNAs and their target genes in genic male sterility lines in Brassica napus by small RNA sequencing

Abstract

Background

Brassica napus is the third leading source of edible oil in the world. Genic male sterility (GMS) lines provide crucial material for harnessing heterosis for rapeseed. GMS lines have been used successfully for rapeseed hybrid production in China. MicroRNAs (miRNAs) play crucial regulatory roles in various plant growth, development, and stress response processes. However, reports on miRNAs that regulate the pollen development of GMS lines in B. napus are few.

Results

In this study, 12 small RNA and transcriptome libraries were constructed and sequenced for the flower buds from the fertile and sterile lines of two recessive GMS (RGMS) lines, namely, “6251AB” and “6284AB”. At the same time, 12 small RNA and transcriptome libraries were also constructed and sequenced for the flower buds from the fertile and sterile lines of two dominant GMS (DGMS) lines, namely, “4001AB” and “4006AB”. Based on the results, 46 known miRNAs, 27 novel miRNAs on the other arm of known pre-miRNAs, and 44 new conserved miRNAs were identified. Thirty-five pairs of novel miRNA-3p/miRNA-5p were found. Among all the identified miRNAs, fifteen differentially expressed miRNAs with over 1.5-fold change between flower buds of sterile and fertile lines were identified, including six differentially expressed miRNAs between “4001A” and “4001B”, two differentially expressed miRNAs between “4006A” and “4006B”, four differentially expressed miRNAs between “6251A” and “6251B”, and ten differentially expressed miRNAs between “6284A” and “6284B”. The correlation analysis of small RNA and transcriptome sequencing was conducted. And 257 candidate target genes were predicted for the 15 differentially expressed miRNAs. The results of 5′ modified RACE indicated that BnaA09g48720D, BnaA09g11120D, and BnaCnng51960D were cleaved by bna-miR398a-3p, bna-miR158-3p and bna-miR159a, respectively. Among the differentially expressed miRNAs, miR159 was chosen to analyze its function. Overexpression of bna-miR159 in Arabidopsis resulted in decreased seed setting rate, and shortened siliques, illustrating that miR159 may regulate the fertility and silique development in rapeseed.

Conclusions

Our findings provide an overview of miRNAs that are potentially involved in GMS and pollen development. New information on miRNAs and their related target genes are provided to exploit the GMS mechanism and reveal the miRNA networks in B. napus.

Peer Review reports

Background

Male sterility is largely due to the abnormal development of either the sporophytic or gametophytic anther tissues. According to its origin or genetic analysis, male sterility can be divided into cytoplasmic male sterility (CMS) and genic male sterility (GMS) [1]. CMS is caused by the interaction between the mitochondrial genome and nuclear genes [2, 3]. GMS is derived from natural mutations in nuclear genes that control stamen development. GMS mutants always show advantages, such as complete and stable male sterility and no potential negative cytoplasmic effect, compared with most CMS mutants. GMS can be further divided into dominant GMS (DGMS) and recessive GMS (RGMS). In rapeseed, DGMS and RGMS are widely used for hybrid rapeseed production [4]. Multiple-allele DGMS was generally accepted and was proposed by Song et al. in 2005. This model presented multiple alleles in one locus inheritance to explain the fertility heredity of a newly reported DGMS line 609AB. In this model, Mf, Ms, and ms are three alleles at the same locus, with a relationship of Mf dominant over Ms and Ms over ms. The recessive allele is for normal fertility. The maintainers and restorers are easily screened. Thus, multiple-allele DGMS are widely used for hybrid rapeseed seed production through the construction of a three-line hybrid system [5, 6].

Additionally, the recessive GMS (RGMS) systems have another distinct advantage. In RGMS systems, most inbred lines can restore their fertility, so hybrids with strong heterosis are easily bred. Several RGMS lines have been successfully commercialized in China, including S45AB [7], 117AB [8], and 9012AB [9]. For S45AB and 117AB, approximately 50% of the fertile plants are required to be artificially removed before hybridizing the rest (50% sterile ones) with restorer lines in hybrid production, because no complete maintainer line is available [10, 11]. However, a three-line hybrid production system was developed for 9012AB [12], and this system has been well documented [10, 13]. The RGMS line 9012AB has been used successfully for rapeseed hybrid production in China. This male sterility was previously thought to be controlled by three independent genes (BnMs3, BnMs4, and BnRf). In 2012, Dong et al. demonstrated a major modification of the sterility inheritance model in 9012A. The modified inheritance model indicated that the male sterility was essentially controlled by two loci (BnMs3 and BnRf). The previously designated BnMs4 locus was just one allele of BnRf; it was then designated as BnRfa, which was designated in addition to BnRfb (the allele from 9012A) and BnRfc (the allele from temporary maintainer). The dominance relationship of the three alleles is in the following order: BnRfa > BnRfb > BnRfc. The BnRf allele-specific molecular markers were identified; these markers would simplify the breeding process involving this RGMS line [14].

Oilseed rape (Brassica napus, 2n = 38, AACC), which has low erucic acid and glucosinolate contents, is the third leading source of edible oil worldwide. In recent years, some conserved and novel miRNAs associated with silique length [15], thickness of pod canopy [16], cadmium stress [17,18,19], flower organ development [20], seed maturation [21, 22], cold stress [23], seed development, and oil synthesis [24] have been widely identified in rapeseed. However, little information is available about the DGMS and RGMS occurrence at the post-transcriptional level in rapeseed.

In this study, to systematically explore the roles of miRNAs and their targets involved in GMS occurrence during pollen development in rapeseed, 12 small RNA and transcriptome libraries were constructed and sequenced for the flower buds from the fertile and sterile lines of two RGMS lines (“6251AB” and “6284AB”). Meanwhile, 12 small RNA and transcriptome libraries were also constructed and sequenced for the flower buds from the fertile and sterile lines of two DGMS lines (“4001AB” and “4006AB”). The aims of this study were to identify known and potential novel miRNAs from the 24 libraries and to analyze the expression profiles of the miRNAs and their targets in relation to DGMS and RGMS during rapeseed microspore development. The results would provide a foundation for evaluating the important regulatory roles of miRNAs in pollen formation and GMS occurrence in rapeseed and other crops.

Results

Analysis of small RNA library data sets and the small RNA profile

To identify miRNAs related to DGMS and RGMS during pollen development, the flower buds were collected from the sterile (6251A, 6284A) and fertile (6251B and 6284B) lines of the RGMS lines. Meanwhile, the flower buds were also respectively collected from the sterile (4001A, 4006A) and fertile (4001B and 4006B) lines of the DGMS lines. Three biological replicates were conducted for each of the eight kinds of samples. Thus, total of 24 sRNA libraries were constructed and deep-sequenced. The raw reads of the 24 sRNA libraries ranged from 20.58 to 42.68 million (Table 1). The raw reads of the 24 sRNA libraries were uploaded to SRA database of NCBI and 24 accession numbers were obtained, including SRX11350295, SRX11350296, SRX11350307, SRX11350312, SRX11350313, SRX11350315, and SRX11350316 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA743414?reviewer=t674c02cj415380e8oldre4s5a). After removing the low-quality reads and contaminated adapter sequences, the clean reads of the 24 sRNA libraries ranged from 19.77 to 41.61 million. The mapped reads were further annotated against the Pfam database and subsequently divided into rRNAs, tRNAs, snRNAs, snoRNAs, ta-siRNA, and others. The endogenous sRNAs were identified as known and novel miRNAs. The average sRNA lengths of the three biological replicates for each sample were calculated, which showed the length distribution patterns of the sRNAs being similar to one another. In general, the majority of the small RNAs ranged from 21 nt to 24 nt in size. The 24 nt small RNAs were the most dominant, followed by 21 nt small RNAs (Fig. 1).

Table 1 Overview of sRNA sequencing reads in Brassica napus
Fig. 1
figure 1

Length distribution of sRNAs in flower buds of A line and B line libraries of DGMS and RGMS in Brassica napus. A Length distribution of sRNAs in 4001A and 4001B. B Length distribution of sRNAs in 4006A and 4006B. C Length distribution of sRNAs in 6251A and 6251B. D Length distribution of sRNAs in 6284A and 6284B

Identification of known and novel miRNAs in B. napus

To identify known miRNAs in B. napus, all mapped small RNA sequences were compared with the known mature bna-miRNA sequences deposited in the miRBase database 22.1. Forty-six small RNAs that have the same sequences with the known bna-miRNAs in miRBase were identified. The numbers of reads of the 46 known miRNAs in 24 libraries were listed in Additional Table S1. Among the 46 known miRNAs, bna-miR159a, bna-miR166a, and bna-miR167c showed very high expression levels.

To predict novel miRNAs in B. napus, BLAST analysis was conducted for all the mapped small RNAs to the B. napus genome sequence in Brassica database and known plant miRNAs in miRBase. The small RNAs that exactly map to the genome sequence but not the known plant miRNAs were classified as candidate novel miRNAs. Five criteria described in the Materials and Methods were used to search for novel miRNAs. As a result, 35 pairs of novel miRNA-3p/miRNA-5p were identified. The mature sequences, reads numbers, positions in chromosomes, precursor sequences and minimum free energy were listed in Table 2. The length distribution of the novel miRNAs was between 18 nt to 26 nt. The length of the novel miRNA precursors ranged from 51 nt to 300 nt with an average length of 154 nt. The minimum free energy ranged from − 240.26 to − 9.9 kcal mol− 1 with an average of − 70.14 kcal mol− 1. The precursor sequences and secondary structures of the novel miRNA were shown in Additional file 1 (Table S2) and Additional file 2 (Fig. S1).

Table 2 Novel miRNAs identified in the fertile and sterile lines of Brassica napus by high-throughput sequencing

Novel miRNA on the other arm of known pre-miRNA

Through sRNA high-throughput sequencing, miRNA-3p and miRNA-5p were found to always be simultaneously present on the 5′ arm and 3′ arm of pre-miRNA secondary structures. To identify novel miRNAs on the other arm of known pre-miRNAs, all mapped small RNA sequences were compared to the known precursor sequences of bna-miRNAs in the miRBase database 22.1. Finally, 27 novel miRNAs on the other arm of known B. napus pre-miRNAs were identified. The miRNA sequences and the number of reads in 24 libraries were listed in Table 3.

Table 3 Identification of novel miRNAs on the other arm of known pre-miRNAs in Brassica napus

Identification of new conserved miRNA families and new miRNA members

To identify new conserved miRNAs in B. napus, all mapped small RNAs were mapped to known plant miRNAs in miRBase and B. napus genome sequences. If the small RNAs can match known plant miRNAs with no more than three mismatches and can exactly map to B. napus genome sequences, then these small RNAs were initially classified as candidate new conserved miRNAs. Five criteria described in the Materials and Methods were used to strictly screen the candidate conserved miRNAs. As a result, 44 miRNAs (22 pairs of miRNAs) belonging to 15 miRNA families were identified (Table 4). Among them, bna-miR159b was a new miRNA member of bna-miR159 family. The rest of the 36 miRNAs (14 pairs of miRNAs) have not been previously reported as bna-miRNAs in miRBase; they show high sequence similarity to some of the known plant miRNAs. The bna-miR158a.1 and bna-miR158a.2 were identified for bna-miR158a member. The two pairs of bna-miR158a shared the same mature sequences. Their precursor sequences were highly similar with each other, and these sequences were from different loci of the B. napus genome. These two pairs of miRNAs were called sub-members. This type of sub-member was also observed for bna-miR159b and bna-miR408a. Four sub-members (bna-miR159b.1, bna-miR159b.2, bna-miR159b.3, and bna-miR159b.4) were identified for bna-miR159b, and two sub-members (bna-miR408a.1 and bna-miR408a.2) were identified for bna-miR408a. This phenomenon suggests that some MIRNA genes might be produced through a replication event from one origin to another one, which results in more copies of the miRNA group. Two members were identified for bna-miR319 and bna-miR398 families. Except the above mentioned five miRNA families, the rest of 10 miRNA families had only one miRNA member (Table 4). The secondary structures of these new conserved miRNAs were shown in Additional file 2 (Fig. S1).

Table 4 Identification of new conserved miRNA families in Brassica napus

Expression profiling of differentially expressed miRNAs in sterile and fertile lines

The normalized expression levels of miRNAs were used for identifying differentially expressed miRNAs between the sterile line and the corresponding fertile line, such as “4001A” and “4001B”, “4006A” and “4006B”, “6251A” and “6251B”, and “6284A” and “6284B”. The known, identified conserved and novel miRNAs were followed to differential expression analysis criteria (qvalue < 0.01 and |log2 (fold change)| > 0.73). As a result, 6, 2, 4, and 10 differentially expressed miRNAs were obtained between the flower buds of “4001A” and “4001B”, “4006A” and “4006B”, “6251A” and “6251B”, and “6284A” and “6284B”, respectively (Fig. 2, Table 5). To further explore the miRNAs involved in the two DGMS lines, a Venn diagram analysis was conducted. The results indicated that two differentially expressed miRNAs (bna-novel_34-5p and bna-novel_31-5p) were shared between the DGMS lines “4001AB” and “4006AB” (Fig. 3A). In addition, three differentially expressed miRNAs (bna-novel_34-5p, bna-miR408a-5p, and bna-miR398a-3p) were shared between the RGMS lines “6251AB” and “6284AB” (Fig. 3B). The novel miRNA “bna-novel_34-5p” was the only miRNA that was simultaneously differentially expressed in the DGMS and RGMS lines.

Fig. 2
figure 2

Relative expression analysis of miRNAs in the flower buds of the A line and B line of RGMS and DGMS in Brassica napus by high-throughput sequencing. Six miRNAs (bule columns) showed significant expression differences between ‘4001A’ and ‘4001B’ flower buds. Two miRNAs (orange columns) showed significant expression differences between ‘4006A’ and ‘4006B’ flower buds. Four miRNAs (red columns) showed significant expression differences between ‘6251A’ and ‘6251B’ flower buds. Ten miRNAs (gray columns) showed significant expression differences between ‘6284A’ and ‘6284B’ flower buds. The Y-axis represents the differentially expressed miRNAs. The X-axis represents the log2(FoldChange). Left indicates the miRNAs enriched in B line flower buds. Right indicates the miRNAs enriched in A line flower buds

Table 5 Differentially expressed miRNAs in “4001AB”, “4006AB”, “6251AB” and “6284AB” libraries and their candidate targets by sRNA sequencing and transcriptome analysis in Brassica napus
Fig. 3
figure 3

Venn diagram showing overlaps of significantly differential expressed miRNAs between different genic male sterility (GMS) lines of Brassica napus. A Two differentially expressed miRNAs were shared between 4001AB and 4006AB lines. B Three differentially expressed miRNAs were shared between 6251AB and 6284AB lines

qRT-PCR was conducted to verify the expression profiles of these differentially expressed miRNAs in deep sequencing. Five differentially expressed miRNAs were chosen for qRT-PCR analysis. The results were consistent with those of deep sequencing. In qRT-PCR, miR158 was up-regulated in “4001B” flower buds compared with that in “4001A”. Novel_34 was greatly up-regulated in “4006B” flower buds (730-fold) compared with that in “4006A”. MiR159 and miR827 were both up-regulated in “6284A” flower buds compared with that in “6284B”. miR398 was up-regulated in “6284B” flower buds compared with that in “6284A” (Fig. 4).

Fig. 4
figure 4

The qRT-PCR analysis of differentially expressed miRNAs between the flower buds of A lines and B lines. The flower buds used for qRT-PCR analysis were collected from corresponding A lines or B lines

Target prediction and identification of differentially expressed miRNAs in sterile and fertile lines

A plant small RNA target analysis server (psRNATarget)-based analysis was performed to predict miRNA target genes with default parameters and a maximum expectation value of 3.5 (https://www.zhaolab.org/psRNATarget/). A total of 560 transcripts were predicted to be targets of the 15 miRNAs (Additional Table S3). In addition, transcriptome sequencing was conducted using the same samples as sRNA sequencing (unpublished data). The differentially expressed and up-regulated mRNAs were predicted as the candidate targets for the differentially expressed and down-regulated miRNAs. At the same time, the differentially expressed and down-regulated mRNAs were predicted as the candidate targets for the differentially expressed and up-regulated miRNAs. As shown in Table 5. Thirty-eight candidate target genes were predicted for the six differentially expressed miRNAs between “4001A” and “4001B”. Eleven candidate genes were predicted for the two differentially expressed miRNAs between “4006A” and “4006B”. Twenty-seven candidate genes were predicted for the four differentially expressed miRNAs between “6251A” and “6251B”. One hundred and eighty-one candidate genes were predicted for the ten differentially expressed miRNAs between “6284A” and “6284B”.

To further demonstrate the potential target genes, 5′ modified RACE was performed using mixed samples from flower buds of the fertile lines (“6284B” and “4001B”). Three target genes were validated using 5′ modified RACE (Fig. 5). Bn.A09.CSD1 (BnaA09g48720D) was cleaved by bna-miR398a-3p. Bn.A09.PPR (BnaA09g11120D) was cleaved by bna-miR158-3p. Bn.Cnn.MYB (BnaCnng51960D) was cleaved by bna-miR159a.

Fig. 5
figure 5

Identification of miRNA target genes in B. napus by 5′ modified RACE. The arrows show the cleavage sites of target mRNA. Watson-Crick pairing and mismatches are indicated by vertical dashes and asterisks

Overexpression of bna-miR159a affected seeds and siliques development in Arabidopsis

Among all the differentially expressed miRNAs in the two DGMS and RGMS lines, bna-miR159a had the highest expression level. To reveal miR159 potential function, two constructs containing pre-miR159a-C6 and pre-miR159a-A7 were transformed to Arabidopsis, and corresponding MIR159OE-1 and MIR159OE-2 transgenic plants were obtained. Five and four lines were obtained for MIR159OE-1 and MIR159OE-2 in T1, respectively. In wild-type plants of Arabidopsis, the transcript of mature miR159a and its targets (AtMYB33 and AtMYB65) were detected in root, stem, rosette leaf, stem leaf, flower, and silique through qRT-PCR. The expression level of mature miR159a was the highest in silique (5875-fold), relatively lower in stem, stem leaf, and flower compared with that in root. The expression levels of AtMYB33 and AtMYB65 were very low and almost undetectable in silique, whereas they were relatively high in stem, stem leaf, and flower compared with that in the root (Fig. 6). In T2 transgenic plants, the transcripts of miR159a and its targets were detected in stem leaf from line 1 of MIR159OE-1 and mixed stem leaf from line 3 and line 4 of MIR159OE-2. The results indicated mature miR159a was overexpressed in MIR159OE-1 (4.04-fold) and MIR159OE-2 (13.6-fold) compared with that in WT. Meanwhile, the transcripts of AtMYB33 and AtMYB65 were suppressed in MIR159OE-1 and MIR159OE-2, especially in MIR159OE-2, compared with that in WT (Fig. 7).

Fig. 6
figure 6

The relative expression levels of miR159 and its targets in different tissues from wild type plants of Arabidopsis. R, root; St, stem; Ro-L, rosette leaf; St-L, stem leaf; F, flower; Si, silique

Fig. 7
figure 7

The relative expression levels of miR159 and its targets in stem leaf of MIR159OE transgenic lines of Arabidopsis. In the T2 transgenic plants, stem leaf from line 1 of MIR159OE-1 and mixed stem leaf from line 3 and line 4 of MIR159OE-2 were used for qRT-PCR analysis

The morphological characters of MIR159OE-1, MIR159OE-2, and WT were observed along with their development processes, especially in the flowering and fruiting periods. No significant difference was observed between transgenic and WT plants in the vegetative growth phase. However, during the reproductive growth period, in the T1 and T2 transgenic plants of MIR159OE-1 and MIR159OE-2, the seed setting rate decreased, and siliques became shorter compared with that in the WT (Fig. 8). The length of siliques from WT, MIR159OE-1, and MIR159OE-2 transgenic plants were measured. In T1, the silique length of WT was approximately 13.4 mm, while in the MIR159OE-1 transgenic plants, the silique lengths of line 1, line 4, and line 5 were 6.7 mm, 4.6 mm and 6.4 mm, respectively. And the silique lengths of line 1, line 3, and line 4 of MIR159OE-2 were 5.4 mm, 5.6 mm, and 6.4 mm, respectively (Fig. 8E). In T2, the silique length of WT was approximately 10.6 mm, while in the MIR159OE-1 plants, the silique length of line 1 was 4.3 mm. The silique lengths of line 3 and line 4 of MIR159OE-2 were 4.4 mm and 4.0 mm, respectively (Fig. 8F). These results indicated that overexpression of MIR159 resulted in significantly shorter siliques and reduced seed setting rate.

Fig. 8
figure 8

The phenotypes of MIR159OE transgenic lines. The T2 transgenic plants (A) and siliques (B-D) of MIR159OE-1 and MIR159OE-2 were observed. The length of siliques was measured in T1 (E) and T2 (F) transgenic plants of MIR159OE-1 and MIR159OE-2. L1, L3, L4, L5 indicate line 1, line 3, line 4, and line 5. More than 30 siliques for each plant were measured. Data represent means of three replicates ± SD. *, P < 0.05, **, P < 0.01. Student’s t-test

Discussion

MiRNAs, as the key post-transcriptional regulators, participate in various biological processes in plant. Recently, an increasing number of studies showed that plant miRNAs were also involved in pollen and fertility development. In 2014, Jiang et al. identified 54 new conserved miRNAs and 25 pairs of novel miRNA/miRNA* in a GMS system of B. campestris ssp. chinensis. Eighteen differentially expressed miRNAs with over two fold change between flower buds of male sterile and fertile lines were identified; they might be involved in the pollen development process [25]. In 2017, Ma et al. verified that the overexpression of MIR158 caused pollen abortion and reduced pollen vitality, which were caused by the degradation of pollen content from the binuclear microspore stage [26]. Dong et al. identified 85 known miRNAs and 120 novel miRNAs, which were expressed during rapeseed anther development in a novel recessive GMS system “CN12AB.” Moreover, 19 and 18 known miRNAs were found to be differentially expressed in 0.5–1.0 mm buds and in 2.5–3.0 mm buds between CN12A and CN12B, respectively. Among these, 14 miRNAs were up-regulated, and 23 miRNAs were down-regulated expressed in CN12A compared with those in CN12B [27]. In this study, to identify miRNAs and their targets involved in pollen development and GMS occurrence in rapeseed, 24 small RNA libraries and transcriptome libraries were constructed and sequenced for the flower buds from the fertile and sterile lines of two RGMS lines (“6251AB” and “6284AB”) and two DGMS lines (“4001AB” and “4006AB”). Based on the sequencing results, fifteen differentially expressed miRNAs with over 1.5-fold change between flower buds of male sterile and fertile lines were identified, including six differentially expressed miRNAs between “4001A” and “4001B”, two differentially expressed miRNAs between “4006A” and “4006B”, four differentially expressed miRNAs between “6251A” and “6251B”, and ten differentially expressed miRNAs between “6284A” and “6284B”. Among them, bna-novel_34-5p was common and differentially expressed between the fertile and sterile lines of “4001AB”, “4006AB”, “6251AB”, and “6284AB”. The results of previous studies and this study verified that miRNAs may play important regulatory roles in rapeseed pollen development and GMS occurrence.

MiR159 is conserved in many plants and is involved in multiple growth and development processes of plants. Allen et al. obtained a mir159ab double mutant in Arabidopsis, which showed pleiotropic morphological defects, including altered growth habit, curled leaves, small siliques, and small seeds [28]. Millar and Gubler verified that miR159 regulated anther development by regulating the expression of its targets, such as MYB33 and MYB65 [29]. Overexpression of miR159 caused the down-regulation of MYB103 transcripts and earlier degeneration of the tapetum and aberrant pollen formation during anther development [30]. In radish, differential expression level of miR159 during anther development was observed among male sterile and maintainer lines. Increased levels of miR159 transcripts decreased the expression of MYB101, thereby inhibiting tapetum development and exine formation [31]. Chen et al. identified 17 differentially expressed miRNAs between long and short siliques of rapeseed, including miR159. Correlation analysis of miR159 and its targets suggested that miR159 repressed cell proliferation to control silique length [15]. Hu et al. found that the overexpression of Bra-MIR159a caused pollen abortion and abnormal pollen germination [32]. In this study, the differentially expressed miR159 in “6284A” and “6284B” was chosen to analyze its function during anther development. The bna-miR159 was overexpressed in Arabidopsis and resulted in decreased seed setting rate, and shortened siliques, illustrating that miR159 may regulate the fertility and silique development of rapeseed. The results of previous reports and the results of the present study verified that miR159 and its target genes might be involved in the regulatory network of pollen development and male sterility and this module is conserved in plants.

Conclusion

A large number of miRNAs were identified during pollen development in the two DGMS and two RGMS lines by deep sequencing. These identified miRNAs included 27 novel miRNAs on the other arm of known pre-miRNAs, 44 new conserved miRNAs, and 35 pairs of novel miRNA-3p/miRNA-5p. Among all the identified miRNAs, 15 differentially expressed miRNAs with over 1.5-fold change between flower buds of male sterile and fertile lines were identified, including six differentially expressed miRNAs between “4001A” and “4001B”, two differentially expressed miRNAs between “4006A” and “4006B”, four differentially expressed miRNAs between “6251A” and “6251B”, and ten differentially expressed miRNAs between “6284A” and “6284B”. The qRT-PCR results of 5 differentially expressed miRNAs (miR158, novel_34, miR159, miR827, and miR398) were consistent with deep sequencing results. The association analysis of small RNA and transcriptome sequencing was conducted, and the analysis results indicated that 257 genes were predicted to be the candidate targets of 15 differentially expressed miRNAs. The results of 5′ modified RACE verified that Bn.A09.CSD1 (BnaA09g48720D), Bn.A09.PPR (BnaA09g11120D), and Bn.Cnn.MYB (BnaCnng51960D) were cleaved by bna-miR398a-3p, bna-miR158-3p, and bna-miR159a. Additionally, overexpression of bna-miR159 in Arabidopsis resulted in decreased seed setting rate, and shortened siliques, illustrating that miR159 may regulate rapeseed fertility and silique development. All the results in our study would provide valuable clues for exploring miRNA-mediated regulatory networks in fertility development of GMS lines in B. napus.

Methods

Plant materials

“6251AB” and “6284AB” are two recessive genic male sterile (RGMS) lines of B. napus. The “6251A” and “6284A” are male sterile lines and the “6251B” and “6284B” are the fertile lines. Their male sterility is controlled by two loci (BnMs3 and BnRf) [14]. The original source of the two RGMS lines was “9012A”, which was identified by Chen et al. [9] and Sun et al. [33]. “4001AB” and “4006AB” are two dominant genic male sterile (DGMS) lines of B. napus. The “4001A” and “4006A” are male sterile lines and the “4001B” and “4006B” are the fertile lines. Their male sterility is controlled by three alleles (Mf, Ms, and ms) at the same locus [5]. The original source of the two DGMS lines was “Yi3A”, which was identified by Li et al. [34]. The above RGMS lines and DGMS lines used in this study have been grown for several generations. The progenies of the two RGMS lines are both segregated into sterile and fertile types during reproduction at a ratio of 1:1. While the two DGMS lines are segregated into sterile and fertile types at a ratio of 3:1. These plant materials were planted in the experimental farm of Zhuanghang comprehensive experimental station of Shanghai Academy of Agricultural Sciences. During flowering stage, mixed flower buds were respectively harvested from more than ten plants of the eight lines. The eight lines were “6251A”, “6251B”, “6284A”, “6284B”, “4001A”, “4001B”, “4006A”, and “4006B”. The eight kinds of samples were quickly frozen in liquid nitrogen and stored at − 80 °C. Three independent biological replicates were collected for each kind of sample.

Small RNA library construction and sequencing

Total RNA for each kind of sample was extracted in three biological replicates using Trizol reagent (Invitrogen, USA). RNA samples with an OD260/OD280 ratio of 2.0 and a total content of more than 2 μg were qualified for small RNA library construction. Then, 24 sequencing libraries (three biological replicates respectively for “6251A”, “6251B”, “6284A”, “6284B”, “4001A”, “4001B”, “4006A”, and “4006B”) were constructed using TruSeqSmall RNA Sample Preparation Kit (Illumina, USA) and then sequenced using Illumina Hiseq 2500/Miseq at Beijing Novogene Bioinformatics Technology Co. Ltd.

Data analysis

Clean reads were obtained by removing low-quality reads, N-containing fragments, and adapters. Then, the length of 18 to 30 nt clean reads were mapped to the B. napus genome sequence (http://brassicadb.agridata.cn/brad/) using Bowtie2 with no mismatches allowed, more details were described in Niu et al. [35]. Unmapped sequences were removed.

The mapped small RNAs reads were aligned to known miRNAs in miRBase22.1. Modified software mirdeep2 [36] were used to predict the potential miRNAs and secondary structures. The ncRNAs includes rRNAs, tRNAs, snRNAs, snoRNAs, and small genome repeat sequences were removed. The rest of sRNA sequences were aligned to B. napus NAT-siRNAs in PlantNATsDB to remove NAT-siRNAs. Then miREvo [37] and mirdeep2 [36] were used to predict novel miRNAs in B. napus. To reveal the differentially expressed miRNAs, the miRNAs expression was analyzed using the DESeq2 [38]. The online software of Venny 2.1.0 was used to draw Venn diagrams.

Identification of conserved and novel miRNAs

The candidate sRNAs were mapped to all the known plant miRNA sequences from the miRBase database (http://www.mirbase.org/). The matched sRNAs with no more than three mismatches were considered as candidate conserved miRNAs, while the unmatched sRNAs were considered as candidate novel miRNAs. In addition, Mfold software was used to predict the secondary structures of pre-miRNAs with the flanking sequences of the candidate small RNAs in the genome [39]. Five criteria must be met for identifying conserved and novel miRNAs [40, 41]. The sequences and lengths, read counts, and positions in chromosome were further analyzed for conserved and novel miRNAs.

qRT-PCR

Total RNA was treated with DNase I (Takara, Japan) to remove residual genomic DNA. For the qRT-PCR analysis of AtMYB33 and MYB65, AtACTIN2 gene was used as internal control. PrimeScript™ II reverse transcriptase (Takara) and oligo (dT) primers were used for first-strand cDNA synthesis. For the qRT-PCR analysis of mature miRNAs, U6 was used as internal control gene, and the first cDNA was synthetized using a miRNA First Strand cDNA Synthesis kit (Stem-loop Method) (Sangon Biotech). The qRT-PCR reactions were performed in a MyiQ2 qRT-PCR detection system (Bio-Rad, www.bio-rad.com/) using iQ SYBR Green supermix (Bio-Rad) [42]. Each experiment was conducted in three biological replicates, and the same sample was performed in three technical replicates. Relative expression levels of miRNAs and their target genes were quantified by using the 2-ΔΔCt method [43]. All the primers used for qRT-PCR analysis are listed in Additional file 1: Table S4.

5′ modified RACE analysis

A mixture of flower buds from the fertile lines (“6284B” and “4001B”) was used for total RNA isolation. The 5′ modified RACE was performed using a FirstChoice™ RLM-RACE Kit (Invitrogen, USA). Total RNA was directly ligated to the 5’RACE oligo. The first-strand cDNA synthesis and the two rounds of PCR reactions were conducted following the manufacturer’s instructions. The PCR products containing the target gene bands were ligated into pGEM-T Easy Vector (Promega, USA) for sequencing [25]. The primers are listed in Additional file 1: Table S4.

Vector construction and plant transformation

Two precursor sequences of bna-miR159a, which were located in A7 and C6 chromosomes of B. napus were designated as pre-miR159a-C6 and pre-miR159a-A7, respectively (Additional file 1: Table S5). The 411 nt and 470 nt genomic fragments containing pre-miR159a-C6 and pre-miR159a-A7 were amplified from B. napus using gene specific primers with endonuclease cleavage sites Sma I and Sal I. Then the fragments were cloned into pCAMBIA1301 binary vector with CaMV 35S promoter. The two vectors were designated as p35S::MIR159a-C6 and p35S::MIR159a-A7, which were introduced into Agrobacterium tumefaciens strains GV3101 and further transformed into Arabidopsis by floral dip. The inflorescences were dipped in the Agrobacterium solution containing sucrose and Silwet-77 for 2 min. The infected plants were cultured for 48 h in the dark environment and then transferred to greenhouse [44]. The seeds were harvested and screened by germination on MS medium containing 25 mg/L hygromycin. The T1 and T2 hygromycin-resistant seedlings were transplanted and grown in greenhouse. Seeds from each transgenic plant were harvested separately. The corresponding Arabidopsis transgenic plants were designated as MIR159OE-1 and MIR159OE-2. The transgenic and wild-type plants were cultivated in the greenhouse under the same environment. The phenotypes of T1 and T2 transgenic plants were observed and recorded. Their silique length data were collected from more than three plants for each transgenic line. More than 30 siliques were measured for each plant.

Availability of data and materials

The raw reads of the 24 sRNA libraries were uploaded to SRA database of NCBI and 24 accession numbers were obtained, including SRX11350295, SRX11350296, SRX11350307, SRX11350312, SRX11350313, SRX11350315, and SRX11350316 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA743414?reviewer=t674c02cj415380e8oldre4s5a). The raw transcriptome data will also be further utilized to excavate differentially expressed mRNAs in the fertile and sterile lines of RGMS and DGMS line.

Abbreviations

GMS:

Genic male sterility

RGMS:

Recessive genic male sterility

DGMS:

Dominant genic male sterility

CMS:

Cytoplasmic male sterility

References

  1. Kaul MLH. Male sterility in higher plants. Berlin: Springer; 1988.

    Google Scholar 

  2. Hanson MR, Bentolila S. Interactions of mitochondrial and nuclear genes that affect male gametophyte development. Plant Cell. 2004;16:154–69.

    Google Scholar 

  3. Chase CD. Cytoplasmic male sterility: a window to the world of plant mitochondrial-nuclear interactions. Trends Genet. 2007;23:81–90.

    CAS  PubMed  Google Scholar 

  4. Fu TD. Breeding and utilization of rapeseed hybrid. 2nd ed. Wuhan: Hubei Science and Technology Press; 2000.

    Google Scholar 

  5. Song LQ, Fu TD, Tu JX, Ma CZ, Yang GS. Molecular validation of multiple allele inheritance for dominant genic male sterility gene in Brassica napus L. Theor Appl Genet. 2006;113:55–62.

    CAS  PubMed  Google Scholar 

  6. Liu J, Lu W, He QB. Genetic analysis and molecular mapping of gene associated with dominant genic male sterility in rapeseed (Brassica napus L.). Gene Genomics. 2008;30(6):523–32.

    CAS  Google Scholar 

  7. Pan T, Zeng FY, Wu SH, Zhao Y. A study on breeding and application GMS line of low eruci acid in rapeseed (B. napus). Oil Crop China. 1988;3:5–8.

    Google Scholar 

  8. Hou GZ, Wang H, Zhang RM. Genetic study on genic male sterility (GMS) material no.117A in Brassica napus L. (in Chinese). Oil Crop China. 1990;2:7–10.

    Google Scholar 

  9. Chen FX, Hu BC, Li QS. Discovery and study of genic male sterility (GMS) material 9012A in Brassica napus L. (in Chinese). Acta Agric Univ Pekinensis. 1993;19(Suppl):57–61.

    CAS  Google Scholar 

  10. Ke LP, Sun YQ, Hong DF, Liu PW, Yang GS. Identification of AFLP markers linked to one recessive genic male sterility gene in oilseed rape, Brassica napus. Plant Breed. 2005;124:367–70.

    CAS  Google Scholar 

  11. Yi B, Chen Y, Lei S, Tu J, Fu T. Fine mapping of the recessive genic male—sterile gene (Bnms1) in Brassica napus L. Theor Appl Genet. 2006;113:643–50.

    CAS  PubMed  Google Scholar 

  12. Chen FX, Hu BC, Li C, Li QS, Chen WS, Zhang ML. Genetic studies on GMS in Brassica napus L.: I. inheritance of recessive GMS line 9012A (in Chinese). Acta Agron Sin. 1998;24:431–8.

    Google Scholar 

  13. Xiao L, Yi B, Chen YF, Huang Z, Chen W, Ma CZ, et al. Molecular markers linked to Bn;rf: a recessive epistatic inhibitor gene of recessive genic male sterility in Brassica napus L. Euphytica. 2008;164:377–84.

    CAS  Google Scholar 

  14. Dong FM, Hong DF, Xie YZ, Wen YP, Dong L, Liu PW, et al. Molecular validation of a multiple-allele recessive genic male sterility locus (BnRf) in Brassica napus L. Mol Breed. 2012;30:1193–205.

    CAS  Google Scholar 

  15. Chen L, Chen L, Zhang X, Liu T, Niu S, Wen J, et al. Identification of miRNAs that regulate silique development in Brassica napus. Plant Sci. 2018;269:106–17.

    CAS  PubMed  Google Scholar 

  16. Chen ZY, Huo Q, Yang H, Jian HJ, Qu CM, Lu K, et al. Joint RNA-Seq and miRNA profiling analyses to reveal molecular mechanisms in regulating thickness of pod canopy in Brassica napus. Genes (Basel). 2019;10(8):591–607.

    CAS  Google Scholar 

  17. Zhou ZS, Song JB, Yang ZM. Genome-wide identification of Brassica napus microRNAs and their targets in response to cadmium. J Exp Bot. 2012;63(12):4597–613.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Jian HJ, Yang B, Zhang AX, Ma JQ, Ding YR, Chen ZY, et al. Genome-wide identification of microRNAs in response to cadmium stress in oilseed rape (Brassica napus L.) using high-throughput sequencing. Int J Mol Sci. 2018;19(5):1431–47.

    PubMed Central  Google Scholar 

  19. Fu Y, Mason AS, Zhang YF, Lin BG, Xiao ML, Fu DH, et al. MicroRNA-mRNA expression profiles and their potential role in cadmium stress response in Brassica napus. BMC Plant Biol. 2019;19(1):570–89.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Wang TY, Ping XK, Cao YR, Jian HJ, Gao YM, Wang J, et al. Genome-wide exploration and characterization of miR172/euAP2 genes in Brassica napus L. for likely role in flower organ development. BMC Plant Biol. 2019;19(1):336–50.

    PubMed  PubMed Central  Google Scholar 

  21. Körbes AP, Machado RD, Guzman F, Almerao MP, de Oliveira LFV, Loss-Morais G, et al. Identifying conserved and novel microRNAs in developing seeds of Brassica napus using deep sequencing. PLoS One. 2012;7(11):e50663. https://doi.org/10.1371/journal.pone.0050663.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Huang D, Koh C, Feurtado JA, Tsang EW, Cutler AJ. MicroRNAs and their putative targets in Brassica napus seed maturation. BMC Genomics. 2013;14:140.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Megha S, Basu U, Joshi RK, Kav NNV. Physiological studies and genome-wide microRNA profiling of cold-stressed Brassica napus. Plant Physiol Biochem. 2018;132:1–17.

    CAS  PubMed  Google Scholar 

  24. Wei W, Li G, Jiang X, Wang Y, Ma Z, Niu Z, et al. Small RNA and degradome profiling involved in seed development and oil synthesis of Brassica napus. PLoS One. 2018;13(10):e0204998. https://doi.org/10.1371/journal.pone.0204998.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Jiang JX, Jiang JJ, Yang YF, Cao JS. Identification of microRNAs potentially involved in male sterility of Brassica campestris ssp. chinensis using microRNAs array and quantitative RT-PCR assays. Cell Mol Biol Lett. 2013;18(3):416–32.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Ma ZM, Jiang JX, Hu ZW, Lv TQ, Yang Y, Jiang JJ, et al. Over-expression of miR158 causes pollen abortion in Brassica campestris ssp.chinensis. Plant Mol Biol. 2017;93:313–26.

    CAS  PubMed  Google Scholar 

  27. Dong Y, Wang Y, Jin FW, Xing LJ, Fang Y, Zhang ZY, et al. Differentially expressed miRNAs in anthers may contribute to the fertility of a novel Brassica napus genic male sterile line CN12A. J Integr Agric. 2020;19(7):1731–42.

    CAS  Google Scholar 

  28. Allen RS, Li J, Stahle MI, Dubroué A, Gubler F, Millar AA. Genetic analysis reveals functional redundancy and the major target genes of the Arabidopsis miR159 family. Proc Natl Acad Sci U S A. 2007;104(41):16371–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Millar AA, Gubler F. The Arabidopsis GAMYB-like genes, MYB33 and MYB65, are microRNA-regulated genes that redundantly facilitate anther development. Plant Cell. 2005;17:705–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Zhang ZB, Zhu J, Gao JF, Wang C, Li H, Li H, et al. Transcription factor AtMYB103 is required for anther development by regulating tapetum development, callose dissolution and exine formation in Arabidopsis. Plant J. 2007;52:528–38.

    CAS  PubMed  Google Scholar 

  31. Zhang W, Xie Y, Xu L, Wang Y, Zhu XW, Wang RH, et al. Identification of microRNAs and their target genes explores miRNA-mediated regulatory network of cytoplasmic male sterility occurrence during anther development in radish (Raphanus sativus L.). Front Plant Sci. 2016;7:1054–69.

    PubMed  PubMed Central  Google Scholar 

  32. Hu Z, Shen X, Xiang X, Cao J. Evolution of MIR159/319 genes in Brassica campestris and their function in pollen development. Plant Mol Biol. 2019;101(6):537–50.

    CAS  PubMed  Google Scholar 

  33. Sun CC, Zhao H, Wang WR, Li YL, Qian XF, Fang GH. Inheritance and utilization of a recessive genic male sterile two-type line 20118AB in Brassica napus L. (in Chinese). Acta Agric Shanghai. 2004;20(1):30–2.

    Google Scholar 

  34. Li SL, Qian YX, Wu ZH. Inheritance of genic male sterility in Brassica napus and its application to commercial production (in Chinese). Acta Agric Shanghai. 1985;1(2):1–12.

    Google Scholar 

  35. Niu CD, Li HY, Jiang LJ, Yan MJ, Li CY, Geng DL, et al. Genome-wide identification of drought-responsive microRNAs in two sets of Malus from interspecific hybrid progenies. Hortic Res. 2019;6(1):75–92.

    PubMed  PubMed Central  Google Scholar 

  36. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52.

    PubMed  Google Scholar 

  37. Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012;13:140–9.

  38. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550–700.

    PubMed  PubMed Central  Google Scholar 

  39. Zuker M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003;31:3406–34.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, Bowman JL, et al. Criteria for annotation of plant MicroRNAs. Plant Cell. 2008;20:3186–90.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. Jiang JX, Lv ML, Liang Y, Ma ZM, Cao JS. Identification of novel and conserved miRNAs involved in pollen development in Brassica campestris ssp. chinensis by high-throughput sequencing and degradome analysis. BMC Genomics. 2014;15:146–58.

    PubMed  PubMed Central  Google Scholar 

  42. Jiang JX, Bai JJ, Li SX, Li XR, Yang LY, He YK. HTT2 promotes plant thermotolerance in Brassica rapa. BMC Plant Biol. 2018;18:127–35.

    PubMed  PubMed Central  Google Scholar 

  43. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25(4):402–8.

    CAS  PubMed  Google Scholar 

  44. Clough SJ, Bent AF. Floral dip: a simplified method for Agrobacterium -mediated transformation of Arabidopsis thaliana. Plant J. 1998;16:735–43.

    CAS  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by grants from the Natural Science Foundation of China (Grant No. 31901502), and Shanghai Agriculture Applied Technology Development Program, China (Grant No.Z20160105). The former fund supported the design of the study and collection, including deep sequencing and data analysis. The second fund supported the writing and publishing of this manuscript.

Author information

Authors and Affiliations

Authors

Contributions

J-XJ and L-YY conceived the project and research plan. J-XJ finished the manuscript. J-XJ, P-FX and Y-JL analyzed the sequencing data, performed the experiments and drafted the manuscript. X-RZ, J-YZ, and J-FZ planted and observed the DGMS plants. Y-LL, M-YJ, and W-RW planted and observed the RGMS plants. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Liyong Yang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Identification of known miRNAs. Table S2. Precursors of 35 novel miRNAs in Brassica napus. Table S3. Target genes were predicted by psRNATarget. Table S4. The primers were used for qRT-PCR and 5, modified RACE analysis. Table S5. The sequences used for constructing MIR159 over-expressed vectors.

Additional file 2: Figure S1.

The secondary structures of new conserved miRNAs and novel miRNAs identified in Brassica napus.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jiang, J., Xu, P., Li, Y. et al. Identification of miRNAs and their target genes in genic male sterility lines in Brassica napus by small RNA sequencing. BMC Plant Biol 21, 520 (2021). https://doi.org/10.1186/s12870-021-03306-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-021-03306-w

Keywords