The novel CMS line (SaNa-1A) generated from progenies of somatic hybrids between B. napus and S. alba and the maintainer line (SaNa-1B) were used in this study. Plant materials were cultivated in the experimental fields of Jiangsu Institute of Agricultural Science in the Lixiahe District (Yangzhou, Jiangsu Province). Different sizes of sterile and fertile anthers were fixed for morphological analysis, aiming to acknowledge the abortion stage. In accordance with the abortion stage (tetrad stage), floral buds (0.6–1.3 mm in length) were collected and frozen at 80°C prior to semithin sections and RNA-seq.

Different sizes of sterile and fertile anthers were fixed in 2.5% glutaraldehyde. After 24 h, the fixative solution was replaced with a fresh solution. After 1 day, the anther was washed three times with phosphate buffer (PB; pH 7.2; 15 min for each wash) and postfixed with 1% osmium tetroxide for 4 h. The postfixed tissue was again washed three times with PB. The anther was dehydrated with a graded ethanol series (50, 70, 80, 90, 95, and 100%) for 15 min each, infiltrated with acetone, and embedded in 812 resin. Semithin sections (1 μm) were cut from the polymerized blocks and stained with 1% toluidine blue for 3 min for light microscope observations (Olympus CX51).

Total RNA from buds (0.6–1.3 mm in length, at abortion stage) of four sterile and fertile individual plants was extracted using RNAiso Plus (Takara, China) in accordance with the protocol of the manufacturer (Wang et al., 2014). The quality of total RNA was analyzed using the Agilent 2100 BioAnalyzer with threshold values of RIN ≥ 8 and 28S:18S RNA≥ 1.5:1. Then, RNA (0.1–4 μg) was precisely quantified using the QUBIT RNA Assay Kit. mRNA was purified with oligo(dT) and fragmented into 120–210 bp. cDNA was synthesized using the SuperScript II Kit (Invitrogen), following end-filling and adding A in the 3′-end. After ligation with the Illumina paired-end adapter, cDNAs were purified twice with AMPure XP Beads (Beckman) to eliminate redundant adapters and amplified with PCR. Finally, a gel purification procedure was conducted to select the fragments ranging from 300 to 350 bp to produce the paired-end library. Fragment size was controlled by the High-Sensitivity ChIP Kit, and the precise concentration of each library was tested with the KAPA qPCR Kit. A 10 μL library (2 nM) was fixed onto cBot and sequenced with Illumina HiSeq 2000.

According to Cox et al. (2010), clean reads were obtained by filtering the adapter sequences and low-quality sequences in the raw data using the FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit), which were then assembled into contigs by the Velvet and Oases software (Zerbino and Birney, 2008). Afterwards, the clean reads were mapped onto the unigenes using Bowtie (Langmead et al., 2009). Function of unigenes was annotated by BLASTX searches against the NR, Swiss-Prot, and COG databases (E < 1e−5). The GO annotations for them were determined using Blast2GO (Conesa et al., 2005), which were then submitted to WEGO for the classification graph (Ye et al., 2006; An et al., 2014). For GO enrichment analysis, a corrected p ≤ 0.05 was chosen as the threshold value. The GO term (p ≤ 0.05) is defined as significantly enriched terms of DEGs. Pathway annotations of differential gene expression were conducted to understand gene function through BLASTX of the KEGG database. The raw reads were deposited in the NCBI Short Read Archive with the accession number SRP075203. The assembled unigenes are shown in Additional File 1.

After filtering the low-quality data, all the tags were aligned to the unigenes. Similar to credibility interval approaches reported for the analysis of SAGE data (Vêncio et al., 2003), we employed IDEG6 (Romualdi et al., 2003) to identify mRNAs showing statistically significant differences based on their relative abundance (as reflected by the total count of individual sequence reads) between the two libraries. The general chi-square test was conducted, as it has been proven to be one of the most efficient tests. Finally, genes with p ≤ 0.01 and fold change ≥ 2 or ≤ 0.5 were marked significantly different between the two libraries.

The DGE results were verified by real-time qRT-PCR analysis, using the same RNA samples for library construction. Two developing stages of buds (~2 mm in length, at tetrad stage; 3–5 mm in length, at or after uninucleate stage) were chosen for expression validation. First-strand cDNA was synthesized using the Revert Aid First Strand cDNA Synthesis Kit (Thermo, USA) from 1200 ng total RNA. Gene-specific primers were designed for the selected unigene sequences (Additional File 2). Reactions were conducted with the SYBR PrimeScript™ RT-PCR Kit (TaKaRa, China) in the Bio-Rad CFX96 instrument. The PCR cycling was denatured using a program of 95°C for 10 s and 40 cycles of 95°C for 5 s and 55°C for 30 s. Three biological replicates for each sample and three technical replicates were conducted, and the relative expression level was calculated using the 2^−ΔΔCt method. Then, B. napus actin was used to normalize gene expression.

The flower organ in rapeseed comprises the sepal, petal, stamen, and pistil. The male reproductive organ in normally developed rapeseed flower has four stamens longer than the other two stamens. Each stamen comprises the filament and anther, which are located at the top of each flower (Figure 1B). However, sterile flowers were visually smaller than fertile flowers, and wrinkles were clearly observed on sterile petals. During the developmental process, the anthers and filaments in sterile flowers were shorter than that in fertile flowers from the tetrad stage. The sterile anthers produced only a few or no pollen, but the pistil was normal (Figure 1A).

Anther development in the fertile and sterile lines were cytologically observed in different developmental stages to determine the abortion stage of the sterile line. At the PMC stage, primary parietal cells were differentiated into the anther epidermis, endothecium, middle wall layers, and tapetum after periclinal division and surrounded the microspore mother cells, which were in the middle of each pollensac (Figure 2B1). At the tetrad stage, PMC cells were separated and formed tetrads after meiosis and were surrounded with callose. The middle wall layers disappeared gradually, and the tapetum cells were abundant with cytoplasm (Figure 2B2). At the uninucleate stage, the microspores were mostly occupied with central vacuole, with the nucleus squeezed to the edge, and tapetum cells started to degrade (Figure 2B3). Afterward, the microspores were differentiated into pollen grains, which were released after the breakdown of pollen sac wall (Figures 2B4,B5). Compared with the maintainer line, tapetum cells in the sterile line were observed with vacuolation during the microsporocyte and tetrad stages (Figures 2A1,A2), although they could form pollensac. In the uninucleate stage, tapetum cells were dissolved and moved to the center of the pollensac, and the cytoplasm was condensed and located at the center of microspores released from the tetrad (Figure 2A3). After the late uninucleate stage, microspores failed to separate and no mature pollen was formed, which were degraded together with tapetum cells (Figures 2A4,A5). The main reason for sterility was the obstacle in the degradation of tapetum cells, which hindered the necessary nutrition that microspores needed. We observed that the abortion of this sterile line occurred during the tetrad stage.

After removing low-quality tags in the raw data, 36,539,702 and 43,903,006 clean tags were obtained in the sterile and fertile lines, respectively (Additional File 3). A total of 195,568 unigenes (1326.8 bp on average) were obtained from clean data (average length = 96 bp) by de novo assembly with the Velvet and Oases software (Additional File 4). Considering the N₅₀ value and average length, we select k-mer = 51 for de novo assembly. Meanwhile, 184,146 unigenes (94.2%) were equal to or greater than 200 bp and 109,168 unigenes (55.8%) were equal to or greater than 1000 bp (Table 1).

For functional annotation, the 195,568 unigenes were subjected to BLASTX searches against the sequences in the Swiss-Prot, non-redundant protein sequences in the National Center for Biotechnology Information (NCBI; NR), and Clusters of Orthologous Groups (COG) databases (E ≤ 1e−5). A total of 122,380 (62.6%), 186,245 (95.2%), and 70,315 (36.0%) unigenes were annotated in the Swiss-Prot, NR, and COG databases, respectively. Meanwhile, 68,429 unigenes (35.0%) were well annotated in all of the databases and 186,752 unigenes (95.5%) were identified with annotation at least in one databases (Figure 3). In the COG database, these unigenes were classified into 24 functional groups, including general function prediction only (23.61% of unigenes), posttranslational modification, protein turnover, chaperones (9.4% of unigenes), and nuclear structure (0.04% of unigenes) (Figure 4). All the unigenes were displayed by searching GO database and classified into three hierarchies, cellular location, molecular function, and biological process. Of the 45 GO groups, unigenes were mainly classified in cell and cell part (33.90% of unigenes), binding (48.37% of unigenes), and metabolic process (33.99% of unigenes). Only a few unigenes were classified into virion, locomotion, and viral reproduction (Figure 5).

The expression abundance of tag-mapped genes was analyzed by counting the number of transcripts per million (TPM). We noted that unigenes in the fertile and sterile buds showed high expression levels in several basic bioprocesses. For example, unigene304 (TPM = 663.76 in Fer and TPM = 1517.84 in Ste) encodes ribulose bisphosphate carboxylase that participates in energy production and conversion; unigene24 (TPM = 610.47 in Fer and TPM = 496.96 in Ste) encodes ubiquitin, which functions in cellular processes and signaling; and unigene78636, which is involved in RNA transport, showed high expression levels in fertile and sterile buds. With the threshold p ≤ 0.01 and fold change ≥2 or ≤0.5, 7811 unigenes (3.92%) were identified with expression differences between sterile and fertile lines, of which 1736 unigenes were upregulated and 6.075 unigenes were downregulated in the sterile line compared with the fertile line (Table 2 and Additional File 5). These DEGs with at least two-fold differences in the two libraries are shown in Figure 6. The red and green dots represented transcripts higher or lower in abundance for more than two-fold in the sterile library, respectively. The blue dots represented transcripts that differed less than two-fold between the two libraries, which were arbitrarily designated as “no difference in expression.”

Different genes cooperate to achieve their biological functions. Pathway-based analysis helps in further understanding the biological functions of DEGs. The upregulated and downregulated unigenes were submitted to the Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatics Annotation Server and classified using the single-directional best hit method (Kanehisa et al., 2008). Most of the upregulated unigenes functioned as transcription factors (TFs). By contrast, a large number of the downregulated unigenes functioned as lipid biosynthesis proteins, chaperones and folding catalysts, protein kinases, prenyltransferases, peptidases, ubiquitin systems, translation factors, and TFs (Table 3). GO annotation of DEGs revealed that, compared with the fertile line, the upregulated unigenes in the sterile line mainly functioned in the binding, cellular process, cell, and cell part. In addition, downregulated unigenes mainly functioned in the catalytic activity, binding, and metabolic process (Figure 7).

GO enrichment analysis of functional significance was conducted to reveal significantly enriched GO terms of DEGs. The GO term with p ≤ 0.05 was defined as significantly enriched GO term of DEGs. For molecular function, significantly enriched GO terms of DEGs included catalytic activity, ATP binding, oxidoreductase activity, transferase activity (transferring phosphorus-containing groups), and zinc ion binding (Additional File 6). For enriched biological processes, 10 GO terms of DEGs were available, which participated in the oxidation–reduction process, pollen development, transmembrane transport, carbohydrate metabolic process, cellular metabolic process, lipid metabolic process, fatty acid biosynthetic process, response to oxidative stress, cell wall modification, and cellular lipid metabolic process (Additional File 7).

The development of pollen is a fundamental and complex process in flowering plants. Pollen is essential for propagation and evolution. As a model plant, Arabidopsis thaliana has been well investigated for this bioprocess. Therefore, all of the unigenes identified in this study were annotated to the TAIR database. Unigenes annotated to 15 genes regulating pollen development (from AG to MS2) were analysized (Higginson et al., 2003; Koizuka et al., 2003; Yang et al., 2003, 2005). One gene encoding tapetum determinant 1 (BnTPD1, unigene40544) was significantly upregulated in the sterile line than in the fertile line (log2ratio = 2.128). We determined that 13 unigenes were downregulated in the sterile line, including BnAG (unigene560, log2ratio = −5.517), BnNZZ/SPL (unigene15428, log2ratio = −1.951), BnEMS1 (unigene4387, log2ratio = −2.569), BnSERK1 (unigene11458, log2ratio = −2.206), BnDYT1 (unigene31429, log2ratio = −3.135), BnMYB99 (unigene2438, log2ratio = −7.081), and BnMYB26 (unigene14326, log2ratio = −2.147). These genes mainly participated in the formation and differentiation of tapetum. Moreover, BnTDF1 (unigene5833, log2ratio = −5.56), BnAMS (unigene472, log2ratio = −6.74), BnMYB103 (unigene25542, log2ratio = −6.639), BnMS1 (unigene17214, log2ratio = −5.98), BnMS2 (unigene1036, log2ratio = −8.279), and BnA6 (unigene90, log2ratio = 8.251) relate to the functional regulation of tapetum at the later developmental stage, such as degradation of callose and formation of pollen wall. We also identified unigenes with no significant expressional difference, and several other genes were unidentified in the sterile and fertile lines (Additional File 8).

Accumulation of starch is necessary for the development of microspore. In other words, abundance of starch is an important characteristic of fertile pollen. Genes participated in carbon metabolism were differentially expressed between fertile and sterile lines. Altered expression of numerous genes involved in pentose and glucuronate interconversions, starch and sucrose metabolism, and amino sugar and nucleotide sugar metabolism were observed (Additional File 5). With regard to pentose and glucuronate interconversions, 13 genes were significantly downregulated, including 6 genes encoding aldehyde dehydrogenase (unigene242, unigene1535, unigene2791, unigene5839, unigene7055, and unigene7131), 4 genes encoding pectinesterase family protein (unigene7776, unigene8408, unigene17104, and unigene25928), and 3 genes encoding pectate lyase family protein (unigene13148, unigene15151, and unigene31020). Moreover, 9 genes related to starch and sucrose metabolism were significantly downregulated, including 5 genes encoding pfkB-type carbohydrate kinase family protein (unigene981, unigene2366, unigene4077, unigene16521, and unigene20485), 1 gene encoding beta-glucosidase (unigene4620), 2 genes encoding glycosyl hydrolase (unigene13516 and unigene17104), and 1 gene encoding cell wall invertase (unigene13728). Of the genes mentioned above, unigene8408, unigene17104 and unigene13148 were not identified in sterile line. Two genes participating in amino sugar and nucleotide sugar metabolism were downregulated, encoding mannose-6-phosphate isomerase (unigene25396, log2ratio = −2.859) and 6-phosphofructokinase (unigene28654, log2ratio = −2.907), respectively.

The mitochondria is an important site for numerous metabolic pathways, including the TCA cycle, respiratory electron transfer, and ATP synthesis (Barrientos, 2003; Reichert and Neupert, 2004; Logan, 2006). The recombination and rearrangement of mitochondrial genomes could cause dysfunction and lead to CMS. Flower phenotypic variation or defection in pollen formation are presumed as secondary effects of mitochondrial mutation, and the primary defect may be a reduction in the efficiency of respiration or the impairment of other mitochondrial functions. We determined that six unigenes participating in the TCA cycle were significantly downregulated, with expressions increased by 3.43–53.47-folds in Fer, such as unigene12291 encoding a subunit of aconitate hydratase, unigene5423 encoding a subunit of citrate synthase, unigene232 encoding phosphoenolpyruvate carboxylase 2, unigene2462 encoding peroxisomal nad-malate dehydrogenase 1, unigene579 encoding 2-oxoglutarate dehydrogenase E1 component, and unigene40435 encoding the subunit B2 of ATP citrate lyase (Additional File 5). Citrate synthase, as a key enzyme of the TCA cycle, functions in the formation of citrate using acetyl-CoA and oxaloacetic acid as substrate. We also identified 11 DEGs related to the electron transport chain and oxidative phosphorylation, which were downregulated in the sterile line compared with the fertile line. For instance, four genes (unigene1243, unigene5691, unigene46531, and unigene31126) encoding subunits of ATPase, three genes (unigene7780, unigene56318, and unigene27486) encoding subunits of pyrophosphorylase, unigene26033 encoding subunits of cytochrome oxidase, unigene15953 and unigene952 encoding NADH dehydrogenase, and unigene1560 encoding subunit of ubiquinol oxidase were included (Additional File 9).

In the present study, oxidoreductase activity is one of the most enriched GO terms. A total of 145 DEGs were involved in oxidoreductase activity, including 36 members acting on the aldehyde or oxo group of donors, with NAD or NADP as acceptor. For instance, fatty acyl-CoA reductase (alcohol-forming)/oxidoreductase (unigene1036 and unigene7520), also called male sterility 2 (MS2), acting on the CH–CH group of donors, with NAD or NADP as acceptor, was more abundant in the fertile line. Other DEGs were 3-chloroallyl aldehyde dehydrogenase/aldehyde dehydrogenase (unigene1535, unigene5839, and unigene7055), aldehyde oxidase (unigene3230, log2ratio = −2.686), NAD or NADH binding/catalytic/glyceraldehyde-3-phosphate dehydrogenase (unigene1766, log2ratio = −1.898), cytosolic factor family protein (unigene12736 and unigene27821), and several unnamed oxidoreductase members (Additional File 5).

Proteins encoding the PPR motif are predicted as site-specific, RNA-binding adaptor proteins that mediate the interactions between RNA substrates and relative enzymes (Lurin et al., 2004; Stern et al., 2004; Shikanai, 2006). PPR proteins play important role in plant mitochondrial biogenesis (Wen et al., 2003). Most of the cloned restorer genes encode mitochondria-targeted PPR proteins (Brown et al., 2003; Desloire et al., 2003; Kazama and Toriyama, 2003; Koizuka et al., 2003; Akagi et al., 2004; Wang et al., 2006). In this study, 12 genes encoding PPR proteins (unigene1911, unigene16513, unigene16881, unigene17672, unigene20560, unigene21040, unigene23771, unigene27924, unigene28403, unigene40662, unigene42447, and unigene44799) were downregulated and one gene (unigene57904) was upregulated (Additional File 5). These PPR proteins were candidates for analyzing the link between nuclear and mitochondrial gene expression.

We randomly selected 16 genes for qRT-PCR assays to confirm the reliability of the Illumina/Solexa sequencing technology. The corresponding primers are listed in Additional File 2. The results showed that the expression of 14 genes were consistent between qRT-PCR and RNA-seq analyses (Figure 8). In qRT-PCR, unigene39076 and unigene36457 showed no difference in the fertile and sterile buds, whereas RNA-Seq analysis indicated a significant difference. The inconsistency between qRT-PCR and RNA-seq analyses of certain genes can be likely attributable to the fact that RNA-seq was more sensitive in the detection of low-abundance transcripts and small expressional changes than qRT-PCR.

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.