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.).

MicroRNAs (miRNAs) are a type of endogenous non-coding small RNAs that play critical roles in plant growth and developmental processes. Cytoplasmic male sterility (CMS) is typically a maternally inherited trait and widely used in plant heterosis utilization. However, the miRNA-mediated regulatory network of CMS occurrence during anther development remains largely unknown in radish. In this study, a comparative small RNAome sequencing was conducted in floral buds of CMS line 'WA' and its maintainer line 'WB' by high-throughput sequencing. A total of 162 known miRNAs belonging to 25 conserved and 24 non-conserved miRNA families were isolated and 27 potential novel miRNA families were identified for the first time in floral buds of radish. Of these miRNAs, 28 known and 14 potential novel miRNAs were differentially expressed during anther development. Several target genes for CMS occurrence-related miRNAs encode important transcription factors and functional proteins, which might be involved in multiple biological processes including auxin signaling pathways, signal transduction, miRNA target silencing, floral organ development, and organellar gene expression. Moreover, the expression patterns of several CMS occurrence-related miRNAs and their targets during three stages of anther development were validated by qRT-PCR. In addition, a potential miRNA-mediated regulatory network of CMS occurrence during anther development was firstly proposed in radish. These findings could contribute new insights into complex miRNA-mediated genetic regulatory network of CMS occurrence and advance our understanding of the roles of miRNAs during CMS occurrence and microspore formation in radish and other crops.

Introduction

MicroRNAs (miRNAs) are a type of endogenous non-coding small RNAs of ~21–24 nucleotides that are known to be important negative regulators of gene expression at transcriptional and post-transcriptional level by mediating mRNA degradation or translational repression (Voinnet, 2009). In plants, primary miRNAs (pri-miRNAs) are transcribed from nuclear-encoded MIR genes by RNA polymerase II and cleaved by Dicer-like1 (DCL1) assisted by the dsRNA binding protein HYL1 to generate miRNA:miRNA* duplexes called pre-miRNAs (Jones-Rhoades et al., 2006; Kurihara et al., 2006; Ruiz-Ferrer and Voinnet, 2009). The duplexes are then methylated by HEN1 and one of the strands combines with the argonaute protein1 (AGO1) to form the RNA-induced silencing complex (RISC), which regulates gene expression through mRNA degradation with nearly perfect complementarity or translational repression with partial complementarity (Yu et al., 2005; Jones-Rhoades et al., 2006; Bodersen et al., 2008).

Cytoplasmic male sterility (CMS) is a maternally inherited trait in plant, which is unable to produce functional pollen, and is a widely observed phenomenon in nearly 200 species (Brown et al., 2003; Hu et al., 2012). CMS lines have been widely used for the production of F1 hybrid seeds and utilization of heterosis in many crops, such as cotton, maize, sorghum, wheat, rice, beet, and rapeseed (Schnable and Wise, 1998; Bentolila et al., 2002; Kubo et al., 2011). In addition to its crucial breeding tools, CMS lines also provide important materials for studying anther and pollen development, and cytoplasmic-nuclear interactions (Chen and Liu, 2014). CMS is usually due to the effect of sterilizing factors found in the mitochondrial genome (Touzet and Meyer, 2014). In most cases, CMS can be restored by nuclear-encoded fertility restorer (Rf) gene(s), which relies on Rf suppressing cytoplasmic dysfunction caused by mitochondrial genes (Eckardt, 2006). High-throughput sequencing now is widely used and has been proven to be an excellent application for the identification of plant miRNAs. As a class of negative regulators, miRNAs have also been identified and characterized during anther development in several plant species, including Arabidopsis (Chambers and Shuai, 2009), Oryza sativa (Wei et al., 2011; Yan et al., 2015), Gossypium hirsutum (Wei et al., 2013), Brassica juncea (Yang et al., 2013), and B. rapa (Jiang et al., 2014). In G. hirsutum, 16 conserved miRNA families were identified during anther development between the Genetic male sterility (GMS) mutant and its wild type. In O. sativa, Wei et al. (2011) identified 292 known miRNAs and 75 novel miRNAs from sporophytic tissues and pollen at three developmental stages. Additionally, many CMS occurrence-associated miRNAs have also been identified in some vegetable crops. In B. rapa, 54 conserved and eight novel miRNA families involved in pollen development were identified (Jiang et al., 2014). In B. juncea, 197 known and 93 new candidate miRNAs during pollen development between CMS line and its maintainer line were also identified (Yang et al., 2013). Although a large number of miRNAs during anther development have been isolated and identified in many crop species, the miRNA-mediated regulatory network of CMS occurrence during anther development remain to be clarified in root vegetable crops.

Radish (Raphanus sativus L. 2n = 2x = 18) is an important annual or biennial root vegetable crop of Brassicaceae family. In recent years, some conserved and novel miRNAs associated with taproot thickening, embryogenesis, flowering-time, and heavy metal stresses had been widely identified in radish (Xu et al., 2013a; Zhai et al., 2014; Nie et al., 2015; Wang et al., 2015; Yu et al., 2015). However, there is little information about the CMS occurrence at the post-transcriptional level in radish. To systematically explore the roles of miRNAs and their targets involved in CMS occurrence during anther development in radish, two small RNA libraries from ‘WA’ (male sterile line), and ‘WB’ (maintainer, fertile line) floral buds of radish were constructed. The aims of this study were to identify known and potential novel miRNAs from the two libraries and investigate the dynamic expression patterns of the CMS occurrence-related miRNAs and their targets during anther development in radish plant. Furthermore, the miRNA-mediated regulatory network of CMS occurrence during anther development was constructed in radish. These results would lay a valuable foundation for elucidating the regulatory roles of CMS occurrence-related miRNAs in radish and facilitate further dissection of the molecular mechanisms underlying microspore formation and CMS occurrence in other crops.

Materials and methods

Plant materials

The radish cytoplasmic male sterile line ‘WA’ and its maintainer line ‘WB’ were used as materials in this study. The ‘WB’ was advanced inbred line through multiple self-pollination for more than 10 generations, while CMS line ‘WA’ was developed through continuously backcrossing with ‘WB’ for more than 10 generations. ‘WA’ had completely aborted anthers without pollen, whereas ‘WB’ had normal anthers with fertile pollen (Figure S1). The materials were planted under normal conditions at Jiangpu Breeding Station of Nanjing Agricultural University, China. According to the cytological characterization of the developmental stages identified with paraffin section technique, the longitudinal length of floral buds reaching 1–1.5, 2–2.5, and 4–5 mm corresponds to the stage of meiosis, tetrad, and early microscope, respectively (Figure S1), which was in highly accordance with results of previous studies (Sun et al., 2012, 2013). Floral buds at three stages were independently collected from the two lines with three biological replicates. Each sample was collected from three randomly selected individual plants and immediately frozen in liquid nitrogen and stored at −80°C for further use.

High-throughput sequencing of small RNAs

Total RNAs were extracted from three stages of floral buds of ‘WA’ and ‘WB’ using Trizol® Reagent (Invitrogen, USA) according to the manufacturer's protocols, respectively. RNAs from the three different stages were equally pooled and used for two small RNA libraries (WA and WB) construction according to previously described procedures (Hafner et al., 2008; Xu et al., 2013b). In brief, small RNA fractions of 18–30 nt were separated and purified from total RNA by 15% denaturing polyacrylamide gel electrophoresis. Then the isolated sRNAs were ligated to 5′ and 3′ adaptors and reverse transcribed to cDNA through SuperScript II Reverse Transcriptase (Invitrogen) and amplified by PCR. Finally, sRNA libraries were sequenced by the Solexa sequencer (Illumina) HiSeqTM 2500.

The clean reads were obtained after removing low quality reads, reads with 5′ primer contaminants or poly-A tails, trimming reads smaller than 18 nt or longer than 30 nt. The remaining unique sequences were then mapped to the radish reference sequences including genomic survey sequences (GSS), expressed sequence tag (EST) sequences and the radish mRNA transcriptome sequences (accession number: SRX1671013) using the SOAP2 program (Li et al., 2009; Xu et al., 2013a). Only perfect matched sequences with no more than two mismatches were retained for proceeding analysis. After using BLAST in GenBank (http://www.ncbi.nlm.nih.gov/genbank/) and Rfam 12.0 (http://rfam.xfam.org/) database, the clean reads compared with the non-coding RNAs (rRNAs, tRNAs, snRNAs, and snoRNA) were removed for further analysis. The remaining matched reads were aligned with known miRNAs in miRBase 21 (http://www.mirbase.org/index.shtml) for radish known miRNAs identification. Then, the unannotated reads were used to predict potential novel miRNAs using Mireap software (https://sourceforge.net/projects/mireap/) according to the previous criteria (Meyers et al., 2008). The stem-loop structure of miRNA precursors were folded by Mfold (http://unafold.rna.albany.edu/?q=mfold/RNA-Folding-Form) (Zuker, 2003).

Differential expression analysis of miRNAs between CMS line and its maintainer line

The frequency of miRNAs from two libraries was normalized to one million by total number of miRNAs per sample (Gao et al., 2012). If the normalized read of a given miRNA is zero, the expression value was set to 0.01 for further use. The differential expression of miRNAs between the two libraries was calculated as: Fold-change = log2 (WA/WB). The P-value was calculated following the previously reported methods (Li et al., 2009; Zhai et al., 2014). The miRNAs with P ≤ 0.05 and fold-change ≥ 1 or ≤ −1 were considered as up- or down-regulated miRNAs between the two libraries during anther development, respectively.

Prediction and annotation of potential targets for CMS occurrence-related miRNAs

The potential target genes of the identified miRNAs were predicted by the plant small RNA target analysis server (psRNATarget; http://plantgrn.noble.org/psRNATarget/) (Dai and Zhao, 2011). The criteria used for target prediction in plants were performed following previous methods (Allen et al., 2005). To understand the biological functions of the targets, gene ontology (GO) analysis were performed by Blast2GO program on the basis of the BLAST searching against the available Nr database in NCBI. In addition, KEGG Orthology Based Annotation System (KOBAS2.0; http://kobas.cbi.pku.edu.cn/home.do) was applied to predict the biological functions of target genes (Xie et al., 2011). Based on the differentially expressed miRNAs and their corresponding targets, the miRNA-targets regulatory network was constructed using Cytoscape_v3.2.1 software (Smoot et al., 2011).

qRT–PCR validation

Quantitative reverse transcription-PCR (qRT–PCR) was employed to evaluate the validity of small RNA sequencing and also to analyze the expression patterns of miRNAs and their targets during different stages. miRNAs and total RNAs were extracted from samples and reverse-transcribed to cDNA using the One Step Primer Script® miRNA cDNA Synthesis Kit (Takara Bio Inc., Dalian, China) and SuperScript® III Reverse Transcriptase (Invitrogen, USA) following the manufacturer's instructions, respectively. All reactions were performed on a BioRad iQ5 sequence detection system (BIO-RAD) and carried out in a total volume of 20 μl including 0.2 μM primer pairs, 2 μl diluted cDNA, and 10 μl 2 × SYBR Green PCR Master Mix (TaKaRa). The PCR amplification reaction was performed following the previous reports (Zhai et al., 2014). The 5.8S ribosomal RNA (rRNA) was used as the reference gene for normalization. All reactions were done in triplicate, the 2−ΔΔ method was used to calculate the relative expression data (Livak and Schmittgen, 2001). The statistical analysis was performed using SPSS 20 software (SPSS Inc., USA) with Duncan's multiple range test at the 5% level of significance. The primers for qRT–PCR were showed in Table S1.

Results

Overview analysis of sequences from small RNA libraries

To identify known and potential novel miRNAs involved in anther development and CMS occurrence, we constructed two small RNA libraries from the floral buds of ‘WA’ and ‘WB’ line. A total of 43,068,458 raw reads were obtained from the two sRNA libraries. After filtering low quality reads, adapter contaminants, and reads smaller than 18 nucleotides, we obtained 20,287,225 (representing 5,528,061 unique sequences), and 21,989,236 (representing 5,682,107 unique sequences) clean reads from WA and WB library, respectively (Table S2). Of these reads, 13.84% were WA library-specific with 42.68% unique sRNAs, 13.88% were WB library-specific with 44.24% unique sRNAs, and 72.28% were present in both with 13.08% unique sRNAs (Table S3).

By comparing with the NCBI GenBank and Rfam databases, these clean reads that matched non-coding sRNAs including rRNAs, snoRNAs, snRNAs, and tRNAs were eliminated. After that, 27,092 (WA) and 27,764 (WB) unique sequences were acquired by querying the unique reads against miRBase 21 (Table 1). The remaining 5,388,388 (WA) and 5,511,728 (WB) unannotated unique reads were used for identification of potential novel miRNAs (Table 1). The length distribution of sRNA reads ranged from 18 to 30 nt in both libraries (Figure 1), and the most abundant sequences in the two libraries ranged from 20 to 24 nt, which is the representative size range of products cleaved by DCLs (Henderson et al., 2006). The most abundant sRNAs in WA and WB library was 21 and 24 nt long, which accounted for 28.97 and 31.47%, respectively.

Table 1

Distribution of small RNAs among different categories in radish.

Category	WA		WB
	Unique sRNAs	Total sRNAs	Unique sRNAs	Total sRNAs
Total	5528061 (100%)	20287225 (100%)	5682107 (100%)	21989236 (100%)
miRNA	27092 (0.49%)	1453994 (7.17%)	27764 (0.49%)	1521955 (6.92%)
rRNA	92089 (1.67%)	2685339 (13.24%)	117885 (2.07%)	3446274 (15.67%)
snRNA	6099 (0.11%)	22072 (0.11%)	6602 (0.12%)	27969 (0.13%)
snoRNA	3383 (0.06%)	10011 (0.05%)	4120 (0.07%)	11837 (0.05%)
tRNA	11010 (0.20%)	1157890 (5.71%)	14008 (0.25%)	514526 (2.34%)
unannotated	5388388 (97.47%)	14957919 (73.73%)	5511728 (97%)	16466675 (74.89%)

Figure 1

Length distribution of small RNAs in WA and WB library.

Identification of known miRNAs during anther development

To identify known miRNAs from the two libraries, the unique sRNA reads were aligned to known miRNA precursors, and mature miRNA sequences in miRBase 21, allowing a maximum of two mismatches. A total of 124 unique reads belonging to 25 conserved miRNA families were identified in the two libraries (Table 2). The distribution of conserved miRNA family members was analyzed (Figure S2). A large part of conserved miRNA families had members of more than three, and miR165/166 family possessed the largest member of 17, followed by miR156/157, and miR169 with 14 and 11 members, respectively. However, some conserved miRNA families including miR158, miR161, miR391, miR395, miR397 miR398, and miR403 had only one or two members. In addition, 38 unique reads belonging to 24 non-conserved miRNA families were also discovered in these two libraries, which contained fewer members as compared with conserved miRNAs (Figure S2).

Table 2

Known miRNA families and their expression abundance in WA and WB library.

Family	Number of members	miRNA reads		Total reads	Ratio (WA/WB)
		WA	WB
CONSERVED miRNA
miR156/157	14	410237	304695	714932	1.35
miR158	2	15565	7641	23206	2.04
miR159	4	1726	1157	2883	1.49
miR160	6	13463	18893	32356	0.71
miR161	1	152	0	152	–
miR162	3	1211	1379	2590	0.88
miR164	5	9831	9100	18931	1.08
miR165/166	17	321161	405255	726416	0.79
miR167	8	318811	276963	595774	1.15
miR168	3	137677	135624	273301	1.02
miR169	11	36871	37179	74050	0.99
miR171	6	3241	4824	8065	0.67
miR172	7	9866	13246	23112	0.74
miR319	4	508	976	1484	0.52
miR390	5	7594	10560	18154	0.72
miR391	2	3769	3177	6946	1.19
miR393	3	1789	1643	3432	1.09
miR394	3	358	348	706	1.03
miR395	2	45236	44	45280	1028.09
miR396	6	3751	2853	6604	1.31
miR397	1	13	60	73	0.22
miR398	2	223	285	508	0.78
miR399	3	108	80	188	1.35
miR403	2	1455	1306	2761	1.11
miR408	4	3431	7780	11211	0.44
NON-CONSERVED miRNA
miR400	2	43	33	76	1.30
miR447	1	189	245	434	0.77
miR482	2	2073	35	2108	59.23
miR529	3	68	55	123	1.24
miR535	2	2229	83	2312	26.86
miR824	2	651	681	1332	0.96
miR825	1	354	479	833	0.74
miR827	3	1932	2558	4490	0.76
miR828	1	5	71	76	0.07
miR829	1	36	98	134	0.37
miR831	1	87	58	145	1.50
miR845	3	5094	5262	10356	0.97
miR854	1	640	468	1108	1.37
miR858	1	0	26	26	0.00
miR859	1	671	0	671	–
miR860	1	248	215	463	1.15
miR948	1	0	336	336	0.00
miR1878	1	0	1584	1584	0.00
miR1885	2	1961	4902	6863	0.40
miR2111	4	91	206	297	0.44
miR2118	1	42998	60406	103404	0.71
miR2199	1	42085	44471	86556	0.95
miR3444	1	0	452	452	0.00
miR5654	1	1298	1705	3003	0.76

The number of miRNA reads differed greatly in the two libraries (Figure S3). For instance, miR156/157 presented the highest expression abundance with 410,237 in WA library, while miRNA165/166 displayed the highest expression of 405,255 copies in WB library. Several miRNA families such as miR167, miR168, miR2118, and miR2199 also displayed extraordinarily high abundance in both libraries, while some other miRNA families (miR400, miR828, miR829, miR831, and miR858) were expressed with relatively low levels of expression with no more than 100 reads in WA and WB library. In addition, the expression levels of different members of the same miRNA family varied drastically (Table S4).

Identification of potential novel miRNAs in floral buds

A total of 30 precursor sequences and 27 novel miRNA families were identified in the two libraries (Table S5). The secondary structures of these predicted novel miRNA precursors were displayed in Figure S4. In addition to secondary structure prediction, identification of complementary sequences of the mature miRNAs is another way to provide forceful evidences for these predicted novel miRNAs (Meyers et al., 2008). Out of these potential novel miRNAs, only seven miRNAs with mature and complementary miRNA*s were detected as the novel miRNA candidates (Table 3). In the present study, the length of these mature miRNAs ranged from 20 to 23 nt, with a distribution peak at 21 nt (60.0%). Furthermore, the length of these potential novel miRNA precursors ranged from 72 to 255 nt with the average length of 142.6 nt. The minimum free energy (MFE) value ranged from −97.83 to −22.9 kcal/mol with an average value of −48.32 kcal/mol. In addition, nine potential novel miRNAs were expressed in both libraries, while a total of 14 and 8 potential novel miRNAs were WA library-specific, and WB library-specific, respectively (Table S5). Most of these potential novel miRNAs had relatively low expression levels when compared with known miRNAs, and the expression levels of miRNA* sequences were obviously less than those of their corresponding mature miRNAs, which was consistent with the viewpoint that miRNA* strands degraded quickly during the biogenesis of mature miRNAs (Rajagopalan et al., 2006).

Table 3

Novel miRNAs with their complementary miRNA.

miRNA name	Reads		Mature sequence (5′–3′)	Arm	Size	LP (nt)	MFE (kcal/mol)	miRNA location
	WA	WB
rsa-miRn3	443	0	TATTCCGACGACAATTCCGACG	5′	22	100	–56.31	CL2831.Contig4,Contig6
rsa-miRn3*	8	0	TCGGAATTCCGTCGGAATATA	3′	21	100	–56.31	CL2831.Contig4,Contig6
rsa-miRn4	519	637	AATGTATGTAGTCCAATCTAT	5′	21	117	–66	CL2870.Contig2
rsa-miRn4*	18	6	ACATTGGACTACATATATTAC	3′	21	117	–66	CL2870.Contig2
rsa-miRn5	1141	1452	GCTTCCATATCTAGCAGTAGG	5′	21	184	–75.8	CL2916.Contig2
rsa-miRn5*	6	12	TACCGATAGATGTGGAAGCGT	3′	21	184	–75.8	CL2916.Contig2
rsa-miRn7	21141	24827	TTTGCGTGAGTATGTGGATGT	5′	21	119	–49	CL4600.Contig2
rsa-miRn7*	32	42	ATCCACATACTCACGAAAATC	3′	21	119	–49	CL4600.Contig2
rsa-miRn9a	771	1062	CGTTCAGTTCTCCTTTTGCTTC	5′	22	106	–47.24	Rsa#S43006900
rsa-miRn9a*	7	28	AGCAAACGAGAATTGAACGGA	3′	21	106	–47.24	Rsa#S43006900
rsa-miRn19	0	186	GAACGATATAAAAGATCATGGA	5′	22	105	–30.2	CL6156.Contig1, Contig2
rsa-miRn19*	0	44	TATGGCCTTTATATCGTATTCG	3′	22	105	–30.2	CL6156.Contig1, Contig2
rsa-miRn24	0	25	GGTGCAGTTCGGGACTGATTG	5′	21	110	–48.8	FD955742
rsa-miRn24*	0	10	ATTGGCTCCCGCCTTGCATCAA	3′	22	110	–48.8	FD955742

LP (nt), The length of precursor; MFE (kcal/mol), Minimum free energy.

Identification of CMS occurrence-related miRNAs during anther development in radish

To identify miRNAs involved in CMS occurrence during anther development in radish, the differential expression of miRNAs in WA, and WB library was analyzed. Based on these rigorous set of criteria above, a total of 28 known and 14 potential novel miRNAs were differentially expressed during anther development (Figure 2, Table S6). Among them, 17 miRNAs including 11 known and 6 novel ones were up-regulated, and 25 miRNAs including 17 known and 8 novel ones were down-regulated. Of these, 15 miRNAs were differentially expressed at a ratio greater than 10-fold, including 13 known, and two novel miRNAs. Especially, two miRNAs, miR395x (17.77-fold) and rsa-miRn3 (11.09-fold) were the most significantly up-regulated known and novel miRNA, respectively (Figure 2). In addition, many of these CMS occurrence-related miRNAs including miR169m, miR171b-3p, miR396b, miR482c-5p, miR1878-3p, and miR3444a-5p were confined to be expressed only in the WA library, whereas miR171a-3p, miR396a, miR482a-5p, and miR859 were only detected in the WB library. The findings suggested that these miRNAs may play critical roles during anther development in radish.

Figure 2

Comparative relative expression of differentially expressed known (A) and potential novel (B) miRNAs between WA and WB library from radish floral buds.

Target prediction of CMS occurrence-related miRNAs in radish

Target prediction is a prerequisite to understand the biological functions of miRNAs during anther development. In this study, a total of 489 target transcripts were predicted for all the identified miRNAs in radish (Tables S7, S8). To further understand the biological functions of miRNAs, the annotation of these target transcripts were classified into three GO ontologies using the Blast2GO program (http://www.blast2go.com), including 21 biological processes, 12 cellular components, and 10 molecular functions (Figure 3). The main terms in biological processes were “cellular process” (GO: 0009987), “metabolic process” (GO: 0008152), “single-organism process” (GO: 0044699), and “biological regulation” (GO: 0065007). In regard to cellular components, “cell” (GO: 0005623), “cell part” (GO: 0044464), and “organelle” (GO: 0043226) were the three most abundant terms. In addition, “binding” (GO: 0005488) and “catalytic activity” (GO: 0003824) were the most abundant subcategories in the molecular functions.

Figure 3

Gene ontology classification of the predicted targets for differentially expressed miRNAs.

To understand the biological functions of the isolated miRNAs in radish, the miRNA-cleaved mRNAs during anther development were identified. In this study, 489 potential target sequences for 53 known, 16 potential novel and 84 unclassified non-conserved miRNAs from the transcripts of WA and WB library were further annotated by BLAST search against Arabidopsis sequences using KOBAS 2.0 program (Tables S7, S8). Among these predicted targets, a large proportion of them are known transcription factor families such as auxin response factors (ARFs), basic-leucine zippers (bZIPs), myb domain proteins (MYBs), and squamosa promoter-binding proteins (SPLs), which could play essential roles in anther development and CMS occurrence of radish. Moreover, several target genes encoding functional proteins play roles in a broad range of biological processes including agamous-like MADS-box protein 16 (AGL16), argonaute (AGO), F-box protein (F-box), NAC domain containing protein 96 (NAC096), pentatricopeptide repeat-containing protein (PPR), and protein TRANSPORT INHIBITOR RESPONSE 1 (TIR1) (Tables 4, S7). To gain further insight into the correlations between miRNAs and their targets, the miRNA-targets regulatory network was constructed (Figure S5, Table S9). Among them, 26 miRNAs including 19 known and 7 potential novel ones, and 87 unique targets formed a total of 93 miRNA–targets pairs with negatively correlated expression during anther development. In general, these results suggested that the differentially expressed miRNAs may play fundamental regulatory roles in diverse aspects of biological processes during anther development of radish.

Table 4

Identified candidate targets for some known and potential novel miRNAs during anther development.

miRNA	Target sequence	Target gene	Target gene annotation
miR156a	Rsa#S41982434	CYP705A15	Cytochrome P450, family 705, subfamily A, polypeptide 15
	Rsa#S43017568	SPL3	Squamosa promoter-binding-like protein 3
	CL289.Contig1	SPL5	Squamosa promoter-binding-like protein 5
	CL2234.Contig1	SPL13	Squamosa promoter-binding-like protein 13
	Unigene3780	OTP82	Chloroplast RNA editing factor
miR158b-3p	Rsa#S42049270	PPR	Pentatricopeptide repeat-containing protein
miR159a	Rsa#S42037487	MYB101	Myb domain protein 101
	Rsa#S41979156	SPL	Putative transcription factor SPL
	CL8717.Contig1	SPL	Putative transcription factor SPL
miR160a	Rsa#S42581764	ARF16	Auxin response factor 16
	Unigene466		Hydroxymethylglutaryl-CoA lyase
miR161	CL1282.Contig1	MSL10	Mechanosensitive channel of small conductance-like 10
	Unigene28541		Transcription initiation factor TFIIE alpha subunit
miR169b	CL2169.Contig1		26S proteasome non-ATPase regulatory subunit 14
miR169m	CL8331.Contig1		Sulfite exporter TauE/SafE family protein
miR171a-3p	FD953436		Peroxisomal nicotinamide adenine dinucleotide carrier
	CL271.Contig2	LTP4	Non-specific lipid-transfer protein 4
miR393a	FD955493	TIR1	Protein TRANSPORT INHIBITOR RESPONSE 1
	Unigene359	EMB2726	Elongation factor Ts family protein
miR395a	Unigene14836		Putative F-box/kelch-repeat protein
miR396a	CL879.Contig1		Transducin/WD-40 repeat-containing protein
	CL6202.Contig1	SIP2;1	Putative aquaporin SIP2-1
miR396b	Unigene20881	AGO5	Argonaute 5
	Unigene22800	PPR	Pentatricopeptide repeat-containing protein
miR396b-3p	Rsa#S41989522		Transducin/WD-40 repeat-containing protein
miR397a	CL379.Contig2		D-glycerate 3-kinase
	Unigene14031		Syntaxin/t-SNARE family protein
miR403	Rsa#S41987411	AGO2	Argonaute 2
	CL3585.Contig3	AGO2	Argonaute 2
miR482c-5p	CL561.Contig2		Serine/threonine protein kinase
	CL561.Contig4		Protein kinase family protein
miR854	Rsa#S42041817		Carboxylate clamp-tetratricopeptide repeat protein HOP2
	Rsa#S41978503	F-box	F-box protein
	CL5880.Contig1	bZIP	bZIP transcription factor
miR1878-3p	Rsa#S42043459		Putative metal tolerance protein C3
miR1885b	CL9579.Contig1	CSLE1	Cellulose synthase-like protein E1
	Unigene1615		Probable 26S proteasome non-ATPase regulatory subunit 3b
miR2111a-3p	CL444.Contig1	AHL19	AT-hook motif nuclear-localized protein
miR2111b-3p	CL4779.Contig1	GLN1.3	Glutamine synthetase cytosolic isozyme 1-3
miR2199	Rsa#S42571099	AGL16	Agamous-like MADS-box protein AGL16
miR3444a-5p	Rsa#S42004586	PRX Q	Peroxiredoxin Q
	Rsa#S42563276	AGP16	Arabinogalactan protein 16
	Rsa#S43011644	FLA3	Fasciclin-like arabinogalactan protein 3
	CL2205.Contig1	F-box	F-box protein
rsa-miRn2	Rsa#S42571626		Putative pectate lyase 18
rsa-miRn10	CL6997.Contig1	GR-RBP2	Glycine-rich RNA-binding protein 2
	CL7005.Contig1	ATR1	NADPH–cytochrome P450 reductase 1
rsa-miRn13	FD957134	HB20	Homeobox-leucine zipper protein ATHB-20
rsa-miRn15	Unigene22510	NAC096	NAC domain containing protein 96
rsa-miRn16	Unigene25057		Bromo-adjacent homology (BAH) domain-containing protein
rsa-miRn20	CL9688.Contig1	PKP-BETA1	Plastidial pyruvate kinase 2
	Unigene18839		Anticodon-binding domain-containing protein

qRT–PCR validation of miRNAs and their targets during anther development

To verify the quality of small RNA sequencing and analyze the expression patterns of CMS occurrence-related miRNAs in radish, a total of 15 miRNAs were randomly selected for qRT-PCR analysis. It was shown that the expression patterns of these miRNAs from qRT-PCR displayed a similar tendency with those from small RNA sequencing (Figure 4). To further study the dynamic expression patterns of CMS occurrence-related miRNAs and their corresponding targets during anther development, a total of 12 predicted target genes, SPL3 (Rsa#S43017568 targeted by miR156a), PPR (Rsa#S42049270 targeted by miR158b-3p), ARF16 (Rsa#S42581764 targeted by miR160a), HRE1 (Rsa#S43010415 targeted by miR164b-3P), TIR1 (FD955493 targeted by miR393a), AGO5 (Unigene20881 targeted by miR396b), Transducin/WD-40 (Rsa#S41989522 targeted by miR396b-3p), F-box (CL2205.Contig1 targeted by miR3444a-5p), HB20 (Rsa#S43028702 targeted by rsa-miRn13), NAC096 (Unigene22510 targeted by rsa-miRn15), RDM4 (CL8993.Contig1 targeted by rsa-miRn17), and UBQ1 (Rsa#S42012413 targeted by rsa-miRn27), were examined by qRT-PCR at three different stages of meiosis, tetrad, and early microspore. As shown in Figure 5, miR158b-3p, miR160a, miR164b-3p, and miR396b-3p were up-regulated and the expression levels maximized at meiosis stage, and then decreased at tetrad and early microspore stage. In addition, miR156a, miR393a, and miR3444a-5p were down-regulated at meiosis stage, and the expression levels then peaked at tetrad stage, but rapidly decreased at early microspore stage. miR396b showed an up-regulated expression pattern and peaked at tetrad stage, and then slightly decreased at early microspore stage. For the novel miRNAs, the expression levels of rsa-miRn13 and rsa-miRn27 were up-regulated at meiosis and tetrad stage, but dramatically decreased to the minimum at early microspore stage. Moreover, rsa-miRn15 was down-regulated at meiosis and tetrad stage, but rapidly increased to the maximum at early microspore stage. Transcripts of rsa-miRn17 reached its maximum at meiosis stage, but sharply declined at tetrad and early microspore stage. Furthermore, some negative correlations could be found between the expression levels of miRNAs and their corresponding target genes during various anther development stages, suggesting that miRNA-mediated mRNA silencing may be involved in CMS occurrence during anther development in ‘WA’ and ‘WB’ line (Figure 5).

Figure 4

Comparison of relative expression levels of miRNAs between qRT-PCR and small RNA sequencing in radish. Data are means ± SD from triplicate assays.

Figure 5

qRT-PCR validation of differentially expressed miRNAs and its corresponding target genes in floral buds of ‘WA’ (in green) and ‘WB’ (in blue). (S1) meiosis stage, (S2) tetrad stage, (S3) early microspore stage. Data are means ± SD from triplicate assays.

Discussion

High-throughput sequencing technology helps identify a large number of miRNAs and targets associated with CMS occurrence during anther development in several plant species (Wei et al., 2011, 2013; Fang et al., 2014; Yan et al., 2015), and provide an effective way to evaluate the expression profiles of miRNAs and targets in different tissues at different developmental stages. The production of functional pollen grain is a prerequisite for the propagation in flowering plants, and the tapetum cell plays a critical role in microspore and pollen formation (Goetz et al., 2001). Unlike the radish CMS line ‘WA’ having no pollen in aborted anthers, its maintainer line ‘WB’ has normal anthers with fertile pollen (Figure S1). Cytological studies show that there is no visible difference between these two lines during the meiosis and tetrad stage (Figure S1). Thereafter, as compared with ‘WB’, the expanded, and vacuolated tapetum cells of ‘WA’ resulted in microspore degeneration and finally aborted anther with no pollen grains (Figure S1). However, few studies on the relationships between miRNAs and CMS occurrence during anther development in radish were conducted. The lack of CMS occurrence-related genes seriously hampered our understanding of molecular mechanism in CMS occurrence, which became an obstacle to utilize the heterosis of radish. To uncover the miRNA-mediated regulatory network of CMS occurrence during anther development, a comparative small RNAome sequencing was conducted in ‘WA’ and ‘WB’ line in this study. To our current knowledge, this study is the first investigation on identification and characterization of miRNAs, and their targets during anther development in radish.

With the application of high-throughput sequencing technology, it has provided an efficient tool to identify a quite comprehensive set of miRNAs at different stages and to reveal the miRNA-mediated regulatory network of CMS occurrence during anther development in plant. In this study, the length distribution of sRNAs suggested that the 24 nt sRNAs were the most abundant, followed by 21 nt sRNAs, which has been reported in Arabidopsis (Voinnet, 2009), Prunus mume (Gao et al., 2012), O. sativa (Ma et al., 2013), and Medicago truncatula (Eyles et al., 2013). The whole frequent percent of 21 and 24 nt small RNAs (28.33 and 30.07%, respectively) in radish was strikingly different from that of B. juncea, which 21 nt RNAs had high abundance (> 95%), and 24 nt RNAs possessed low frequency (1.1%) (Yang et al., 2013). Interestingly, the same tendency also existed when compared with B. rapa in which 24 nt sRNAs were the most dominant, followed by 21, 22, and 23 nt small RNAs (Jiang et al., 2014), it could be speculated that the genetic relationship between radish and B. rapa is closer than that between radish and B. juncea in the process of evolution.

Identification a set of miRNAs is a crucial step to promote our understanding of miRNA-mediated regulatory network of anther development and CMS occurrence. Recently, numerous studies have presented that the majority of known miRNAs in plantae are evolutionarily conserved (Chen et al., 2012; Barvkar et al., 2013). The diversity of known miRNA families in radish might be decided by the abundance and number of members. In the present study, a large number of conserved miRNAs expressed relatively higher levels compared with non-conserved ones, which was in agreement with previous researches in other species (Gao et al., 2012; Wang F. D. et al., 2012; Wang Z. J. et al., 2012; Fang et al., 2014). In addition, several studies have reported a number of known and potential novel miRNAs involved in anther development and CMS occurrence in B. juncea (Yang et al., 2013), B. rapa (Jiang et al., 2014), Citrus reticulata (Fang et al., 2014), G. hirsutum (Wei et al., 2013), and O. sativa (Yan et al., 2015), which greatly enhanced our knowledge of the regulatory roles of miRNAs in CMS occurrence. In this study, 28 known miRNAs were differentially expressed and the majority of these miRNAs were down-regulated during anther development. The differential expression patterns of rsa-miR160a and rsa-miR169b were consistent with those observed in O. sativa (Yan et al., 2015). Moreover, the expression pattern of rsa-miR396b and rsa-miR171a-3p was similar to that identified in G. hirsutum and B. rapa, respectively (Wei et al., 2013; Jiang et al., 2014). Interestingly, the targets of the miR160 contain three critical regulators, ARF10, ARF16, and ARF17, which are important in mediating gene expression response to the plant hormone auxin and regulating floral organ formation (Mallory et al., 2005; Wang et al., 2005; Chapman and Estelle, 2009; Liu et al., 2010). The expression level of rsa-miR160a was down-regulated in ‘WA’ and validated by qRT-PCR (Figures 5, S5). Thus, it could be speculated that the decreased abundance of rsa-miR160a may partially increase the expression of ARF16, finally resulting in abnormal pollen development in the sterile line ‘WA’.

According to the negative correlation between differentially expressed miRNAs and their corresponding targets (Figure S5), a hypothetical schematic model of miRNA-mediated regulatory network of CMS occurrence during anther development in radish was put forward (Figure 6). As shown in the regulatory network, targets of these differentially expressed miRNAs containing important transcription factors (TFs) and functional proteins are involved in many biological processes, including auxin signaling pathways, signal transduction, miRNA target silencing, floral organ development, and organellar gene expression. For instance, SBP-box genes targeted by miR156, a group of TFs with significant regulatory functions controlling the transition from the vegetative phase to the floral phase in Arabidopsis, O. sativa, and Zea mays (Chuck et al., 2007; Gandikota et al., 2007; Jiao et al., 2010). It was reported that three genes, LEAFY, FRUITFULL, and APETALA1, are directly activated by SPL3 to regulate the timing of flower formation (Yamaguchi et al., 2009). Additionally, multiple SPL genes can lead to fully fertile flowers and regulate cell division and differentiation in Arabidopsis (Xing et al., 2010). In the present study, up-regulation of the rsa-miR156a decreased the expression of SPL3 in ‘WA’ compared to ‘WB’ (Figure 6, Table S7), leading to disordered floral organ development, cell division, and differentiation in radish. MiR159 is required for normal anther development, in which it regulates the expression of genes encoding MYB TFs (Achard et al., 2004; Tsuji et al., 2006). MYB TFs are involved in the control of plant development, determination of cell fate and identity, primary, and secondary metabolism (Stracke et al., 2007; Gonzalez et al., 2008; Kang et al., 2009). AtMYB103, specifically expressed in tapetums and middle layers of anthers, is important for pollen development, especially the pollen exine formation (Zhang et al., 2007; Chen et al., 2014). Down-regulation of the AtMYB103 resulted in earlier degeneration of tapetum and pollen grains aberration during anther development in A. thaliana (Zhang et al., 2007). In rice, anther and pollen defect in floral organ development are also found in the loss-of-function mutations of MYB (Kaneko et al., 2004). In the present study, the rsa-miR159a was found to be up-regulated in ‘WA’ line compared to ‘WB’ line (Figures 6, S5), indicating that the increased abundance of rsa-miR159a partially decreased the expression of MYB101, hampering normal tapetum development, callose dissolution, and exine formation in radish anthers. Moreover, AGL16, belonging to MADS-box transcription factors, was identified to be targeted by rsa-miR2199. The MADS-box TFs are essential regulators of the development of the floral meristems and floral organs in plants (Chen et al., 2014). These evidences indicated that rsa-miR2199 might be an essential component of gene regulatory network that involved in radish CMS occurrence (Figure 6).

Figure 6

The hypothetical schematic model of miRNA-mediated regulatory network of CMS occurrence during anther development in radish. The up- and down-regulated miRNAs are in red and green, respectively. Agamous-like MADS-box protein AGL16 (AGL16), argonaute 2 (AGO2), argonaute 5 (AGO5), auxin response factor 16 (ARF16), basic leucine-zipper (bZIP), F-box protein (F-box), myb domain protein 101 (MYB101), NAC domain containing protein 96 (NAC096), pentatricopeptide repeat-containing protein (PPR), protein TRANSPORT INHIBITOR RESPONSE 1 (TIR1), squamosa promoter-binding like protein 3 (SPL3).

Apart from key TFs, a variety of genes which encode important functional proteins, such as PPR proteins, F-box proteins, AGO proteins, and protein TRANSPORT INHIBITOR RESPONSE 1 (TIR1), were also considered to play important roles in CMS occurrence during anther development. PPR protein genes were identified as targets of miR158 (Lurin et al., 2004; Sunkar and Zhu, 2004). Previous studies indicated that PPR proteins are mostly located in the mitochondria and chloroplast and play crucial roles in pollen development, specific RNA sequence binding, post-transcriptional splicing and mRNA stability regulating (Okuda et al., 2006; Wang et al., 2006; Saze and Kakutani, 2007; Fujii and Small, 2011). In addition, some PPR proteins have also been identified as fertility-restoring genes (Rf) for CMS occurrence (Desloire et al., 2003; Wang et al., 2008; Yasumoto et al., 2009). In this study, rsa-miR158b-3p targeting the gene encoding PPR protein was up-regulated and the PPR gene was suppressed in ‘WA’ line compared with ‘WB’ line, and it could be suggested that the regular expression of CMS-associated mitochondrial genes and suppression of PPR gene result in sterility in radish ‘WA’ line (Figures 5, S5). Moreover, F-box proteins are involved in the regulation of various developmental processes in plants, including floral meristem, floral organ identity determination, and photomorphogenesis (Jain et al., 2007). The expression of rsa-miR3444a-5p was down-regulated at meiosis stage, and then peaked at tetrad stage, but rapidly decreased at early microspore stage, and a negative correlation was found between the expression levels of rsa-miR3444a-5p and its target gene which encoding F-box protein at three different stages according to the qRT-PCR analysis (Figure 5). In addition, F-box gene targeted by osa-miR528 was found to be involved in the regulation of the abortion process in male sterile line of rice. Moreover, the other 23 genes including APG2, AGP16, FIO1, FLA3, FLA5, NAC083, NSP5, TRP1, and VIP1 were also the targets of rsa-miR3444a-5p, indicating that the miRNA has multiple effects on the targets (Figure S5). All of these genes targeted by rsa-miR3444a-5p might function together to regulate the CMS occurrence during anther development in radish. Additionally, AGO proteins were reported to be involved in diverse biological processes including hormone response, developmental regulation, and stress adaptation (Yang et al., 2013). Up-regulation of TIR1 enhances auxin sensitivity, and causes altered leave phenotype and delayed flowering (Chen et al., 2011). In this study, AGO2 and AGO5 was targeted by miR403 and miR396b, respectively, and TIR1 was targeted by miR393a, indicating miR403, miR396b, and miR393a might modulate the hormone response to play roles in the microspore development and CMS occurrence.

In summary, CMS occurrence-associated miRNAs and their targets between the male sterile line ‘WA’ and its maintainer line ‘WB’ were firstly identified and characterized in radish. These results provide a valuable foundation for unraveling the complex miRNA-mediated regulatory network of CMS occurrence and facilitate further dissection of roles of miRNAs during CMS occurrence and microspore formation in radish and other crops.

Author contributions

WZ, YX, and LL designed the research. WZ and YX conducted experiments. LX and YW participated in the design of the study and performed the statistical analysis. WZ and YX analyzed data and wrote the manuscript. XZ, RW, YZ, and EM helped with the revision of manuscript. All authors read and approved the manuscript.

Conflict of interest statement

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.