Comprehensive analysis of Ogura cytoplasmic male sterility-related genes in turnip (Brassica rapa ssp. rapifera) using RNA sequencing analysis and bioinformatics

Ogura-type cytoplasmic male sterility (Ogura-CMS) has been widely used in the hybrid breeding industry for cruciferous vegetables. Turnip (Brassica rapa ssp. rapifera) is one of the most important local cruciferous vegetables in China, cultivated for its fleshy root as a flat disc. Here, morphological characteristics of an Ogura-CMS line ‘BY10-2A’ and its maintainer fertile (MF) line ‘BY10-2B’ of turnip were investigated. Ogura-CMS turnip showed a reduction in the size of the fleshy root, and had distinct defects in microspore development and tapetum degeneration during the transition from microspore mother cells to tetrads. Defective microspore production and premature tapetum degeneration during microgametogenesis resulted in short filaments and withered white anthers, leading to complete male sterility of the Ogura-CMS line. Additionally, the mechanism regulating Ogura-CMS in turnip was investigated using inflorescence transcriptome analyses of the Ogura-CMS and MF lines. The de novo assembly resulted in a total of 84,132 unigenes. Among them, 5,117 differentially expressed genes (DEGs) were identified, including 1,339 up- and 3,778 down-regulated genes in the Ogura-CMS line compared to the MF line. A number of functionally known members involved in anther development and microspore formation were addressed in our DEG pool, particularly genes regulating tapetum programmed cell death (PCD), and associated with pollen wall formation. Additionally, 185 novel genes were proposed to function in male organ development based on GO analyses, of which 26 DEGs were genotype-specifically expressed. Our research provides a comprehensive foundation for understanding anther development and the CMS mechanism in turnip.


Introduction
As an important and valuable resource, male-sterile varieties are extensively exploited in crop hybrid breeding. Cytoplasmic male-sterility (CMS) is a category of male-sterility resulted from PLOS ONE | https://doi.org/10.1371/journal.pone.0218029 June 14, 2019 1 / 29 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 a genomic conflict between the mitochondrial and nuclear genomes, and has been extensively utilized [1]. Various types of CMS have been developed and adopted in plant breeding [2]. It has been proposed that normal microsporogenesis needs appropriate timing of tapetum degeneration and specific gene expression [3]. In CMS system, this elaborate process is complex because of the mitochondrial retrograde signaling pathway and the interaction of nuclear and organelle genomes [4][5][6][7]. Considerable variations in morphological phenotype of anther development, particularly of microspore and tapetum behaviors, arise with different nuclear backgrounds and/or cytoplasmic genotype [8]. Mostly, for example, CMS causes the premature degradation of the tapetal cells [9]. However, CMS sorghum and CMS-T maize show a persistent tapetum which likely inhibits nutrient delivery, resulting in the failure of microspore development [8,10].
Ogura-type CMS was first discovered in Japanese radish (Raphanus sativus) and is now widely applied in the breeding of Brassicaceae crops, such as Brassica napus, B. juncea, and B. oleracea, providing a classic model to probe the role of nuclear-cytoplasmic genome interactions [11]. Although interactions of specific mitochondrial gene orf138 with different nuclear backgrounds have been reported to be responsible for Ogura-CMS, diverse floral behaviors attributed to the same cytoplasm in different species have not been fully investigated [11]. For example, anther morphology in Ogura-CMS B. napus is normal, whereas pollen development is impaired and sensitive to temperature [12][13][14][15][16]. However, Ogura-CMS Chinese cabbage (B. rapa ssp. pekinensis) shows reduced plant height and delayed flowering, has shorter filaments, and produces few and infertile pollen grains in indehiscent anthers [17]. All these differences suggest the presence of various regulatory mechanisms and/or multiple regulatory pathways in Brassica spp.
Recently, numerous candidate genes involved in CMS have been discovered in different species, such as onion (Allium cepa), cabbage (B. olerace var. capitata), rice (Oryza sativa), and pepper (Capsicum annuum) [7,9,18,19]. In addition, participation of miRNAs and non-coding RNAs is becoming increasingly evident in retrograde regulation of CMS [20][21][22][23]. Despite previous extensive work, no specific retrograde pathway has been reported for CMS to date and the regulatory mechanism of CMS is still largely unknown.
Exploring the molecular mechanisms underlying CMS is of great importance for improving seed yield in many crop species, especially in crucifers. As a Brassica root crop, turnip (Brassica rapa ssp. rapifera) has been important for human consumption for thousands of years [24]. In this study, the morphological characteristics of an Ogura-CMS line 'BY10-2A' and its maintainer fertile (MF) line 'BY10-2B' of turnip were investigated, and a detailed RNA sequencing (RNA-Seq) analysis for inflorescences in turnip was conducted. These data provide a comprehensive view on the dynamic gene expression networks and their potential roles in controlling anther development. Using pairwise comparisons, we identified 5,117 DEGs, which might respond to the mutation of the mitochondrial ORF138 locus. Among them, 185 novel genes were proposed to function in male organ development based on GO analyses. These findings provide a comprehensive insight into the regulatory networks responsible for Ogura-CMS tapetum abnormality and pollen abortion in turnip, and demonstrate that cytoplasmic retrograde regulation is probably a principal molecular mechanism for CMS in turnip.

Plant materials and growth conditions
Previously, the Ogura-CMS line 'BY10-2A' of B. rapa ssp. rapifera was developed by inter-specific hybridization between B. rapa ssp. chinensis as the Ogura-CMS cytoplasm donor and fertile B. rapa L. ssp. rapifera, followed by 10 recurrent generations of back-crossing. The Ogura-

RNA sample collection and total RNA isolation
After flowering, all floral buds of an inflorescence from the Ogura-CMS and MF lines of B. rapa ssp. rapifera were collected. In each case, samples were harvested and pooled from ten individual plants with transcriptome profiles representing 'f' difference, then immediately frozen in liquid nitrogen and stored at -70˚C until RNA isolation. For biological repetitions, RNA was extracted from three samples using the EASYspin Plant RNA kit (Aidlab Biotechnologies Corporation, Bejing, China). RNA quality and quantity were characterized on a 1% agarose gel, and determined with a NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and a Qubit RNA Assay Kit in Qubit2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Illumina sequencing and de novo transcriptome assembly
A total of 3 μg RNA per sample was used for library preparation using a NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA). The library preparations were then sequenced on an Illumina Hiseq 2000 platform by the Biomarker Biotechnology Corporation (Beijing, China). Clean reads were filtered from the raw reads by removing reads containing adapter, poly-N and low-quality reads. All the downstream analyses were based on clean data with quality determined by Q20, Q30, GC-content and sequence duplication levels. Then de novo transcriptome assembly was accomplished using the Trinity platform [25] with min-kmer-cov set to 2 by default and all other parameters set to default. The RNA-Seq data were uploaded to the Sequence Read Archive of the National Center for Biotechnology Information (NCBI) (DOI: https://www.ncbi.nlm.nih.gov/sra/PRJNA505114; accession number: PRJNA505114).

Annotation of differentially expressed genes
Differentially expressed genes (DEGs) were screened out using the DESeq (2010) R Package. Genes with a Benjamini-Hochberg false discovery rate (FDR) < 0.05 and a log 2 fold change (FC) � 1 or � -1 in each pairwise comparison were assigned as differential expressed. DEG sequences were blast in the NCBI non-redundant protein (Nr) database, Swiss-Prot database and orthologous groups of genes (eggNOG) database, and also aligned to the Clusters of Orthologous Group (

Real-time reverse transcription polymerase chain reaction (RT-PCR) validation
Sixteen DEGs were randomly selected for validation using real-time RT-PCR. Residual RNA samples for DEG analysis were transcribed into cDNA with a HiScrip II 1 st Strand cDNA Synthesis Kit (Vazyme Biotech, Nanjing, China) as the template. ACT7 was used as the normalization control [29]. Real-time RT-PCR was performed on an Applied Biosystems 7500 Real-time PCR System (ThermoFisher, MA, USA) and the relative expression levels were analyzed. Three technical repeats were performed. All primers used were listed in S1 Table.

The Ogura-CMS turnip displays complete male sterility
Despite 10 generations of back-crossings, Ogura-CMS plants showed a reduction in fleshy root size compared with its MF line (Fig 1A-1C). Furthermore, during the reproductive growth phase, the Ogura-CMS line was distinguishable from the MF line by its short filaments and withered white anthers (Fig 1D and 1E). The floral buds from the Ogura-CMS line showed the same developmental pattern as those of its MF line (Fig 2A and 2C). However, compared with the yellow and plump anthers of the MF line (Fig 2B), no pollen was observed in the mature Ogura-CMS anthers (Fig 2D). Other than pollen absence and stamen abnormality, floral organs presented normal morphologies (S1 Fig), including female gametophytes as pollinating the MF pollen grains on the Ogura-CMS stigma led to normal silique growth and full seed set (S2 Fig).

Aberrant anther development occurs during transition from microspore mother cells to tetrads
To determine the precise stage at which Ogura-CMS anther shrinkage begins, semi-thin sections of anthers from various developmental stages were prepared and further analyzed by microscopy. No differences in microspore development were observed inside the anther locules during the process of meiotic division up to the tetrad stage (Fig 3A, 3B, 3F and 3G). Each Ogura-CMS tetrad ( Fig 3G) contained four microspores, similar to those of the MF line ( Fig  3B), indicating that meiosis is normal in the Ogura-CMS line. Both Ogura-CMS and MF tetrads were surrounded by a tapetum, a middle layer, an endothecium, and an epidermis from the inside out at the tetrad stage (Fig 3B and 3G). However, it was noticeable that the Ogura-CMS tapetum (Fig 3G) swelled at the center of the locule, and the cytoplasm was distinguishably clear from that of the MF line ( Fig 3B). After the tetrad stage, the MF middle layer degenerated, then disappeared by the uninucleate stage, and the MF anthers released microspores that developed into mature pollen grains (Fig 3C-3E). However, the collapse of the Ogura-CMS microspores started at the uninucleate microspore stage and was accompanied by extensive degeneration (Fig 3H). Moreover, the middle layer persisted at this point, together with significantly enlarged tapetum, crushing the free microspores to the locule center. At later stages of development, the collapse of the microspore was even more remarkable due to the degeneration of its entire contents ( Fig 3I). Eventually, defective microspore development and early clearing of tapetal cytoplasm led to shrunken anthers with collapsed locules and remnants of pollen grains adhered to the inner face of the epidermis (Fig 3J).   Abnormal Ogura-CMS microspore development was further confirmed by TEM. Ogura-CMS microspores underwent similar development to those of MF from the microspore mother cell stage to the tetrad stage ( Fig 4A, 4B, 4F and 4G) and, at the uninucleate stage, were distinguishable from MF microspores (Fig 4C and 4H). At the uninucleate stage, the MF microspores had almost complete basic intine and exine structure ( Fig 4C). Exine was comprised of inner nexine and outer sexine. Sexine further possessed a three-dimensional structure composed of baculae and a roof-like tectum, whereas the bilayer nexine, consisting of nexine I and nexine II, was laid down on the intine layer ( Fig 4K). Subsequently, bicellular microspores were generated concurrently with the size increase of exine, and the mature exine structure was visually completed at this point (Fig 4D and 4L). Finally, the mature pollen grain was complete with tryphine ( Fig 4E and 4M). However, after dissolution of the callose wall, free Ogura-CMS microspores were deformed with a thin and incomplete exine layer (Fig 4H). The nexine I was the last layer overlying the microspore plasma membrane but not the intine layer ( Fig 4N). By the bicellular stage, with an empty body, pollen grains were prepared with incomplete-developed layers of exine and tryphine (Fig 4O), which developed into remnants at the dehiscent stage ( Fig 4J).
Coordinated with microspore development, visible changes occurred to the surrounding walls in the anther locules and were also observed when using TEM. No defects on Ogura-CMS surrounding walls were detected at the microspore mother cell stage ( Fig 5F) compared with those of the MF line ( Fig 5A). Anther primordia which were enclosed in an epidermis, differentiated inwardly into the endothecium, middle layer, tapetum, and microspore mother cells. At the tetrad stage, MF tapetal cells became mature and vacuolated, with a heterogeneous density in the cytoplasm ( Fig 5B). However, the anther pattern was significantly altered in Ogura-CMS tapetum at this stage. Ogura-CMS tapetal cells enlarged and swelled to expand to the center of the locules, with larger vacuoles and a clearing cytoplasm (Fig 5G). When the callose wall totally dissolved and free individual microspores released into the anther locules, a large number of elaioplasts emerged in the MF tapetum, and remnants of the middle layer were absent in the MF anther ( Fig 5C), but still clearly present in the Ogura-CMS anther ( Fig  5H). In addition, Ogura-CMS tapetal cells were full of shedding materials such as tapetosomes and significantly swelled, crushing the free microspores to the locule center ( Fig 5H). This premature degradation of the tapetum was obvious, which fulfilled the programmed cell death ahead of schedule. All that remained of the Ogura-CMS tapetum was an almost empty shell at the bicellular stage ( Fig 5I), but integral tapetal cells with a large amount of elaioplasts were still observed in the MF anther ( Fig 5D). Moreover, the endothecium also appeared abnormal in the Ogura-CMS anther, devoid of any content in the cytoplasm at the bicellular stage ( Fig  5I), and had disappeared completely at the mature pollen stage (Fig 5J).
Overall, the Ogura-CMS anthers showed two distinct defects that occurred during the transition from microspore mother cells to tetrads: the failure of microspore development and the swollen tapetum layer. Defective microspore production and premature tapetum degeneration during microgametogenesis led to complete male sterility of the Ogura-CMS line.

RNA-Seq analysis on inflorescences of the Ogura-CMS line and its maintainer fertile line in turnip
To explore the molecular basis for the morphological differences in anther and microspore development described above, RNA-Seq analyses were conducted to generate transcriptome profiles of the whole inflorescences from the Ogura-CMS and MF lines. RNA-Seq analysis was performed with three biological replicates for each. After removing low-quality reads, an average of 24.2 × 10 6 clean reads per library were generated (S2 Table). The de novo assembly resulted in a total of 84,132 unigenes (S3 Fig and S3 Table), which were proposed to be expressed during floral bud development in turnip. To verify the quality of the RNA-Seq data, real-time RT-PCR analysis was conducted on 16 randomly selected genes (Fig 6). The strong correlation between the RNA-Seq and real-time RT-PCR results indicated high reliability of our transcriptomic profiling data.
To determine whether the genes regulating pollen and tapetum development are defective or expressed abnormally, differential expression of genes during thereproductive development of the MF and CMS lines was analyzed. The reads were mapped to unigenes and quantified to show gene expression abundance by Fragments Per Kilobase of transcript per Million mapped reads (FPKM) [30]. The FPKM expression data were tested by correlation analysis to evaluate sampling between biological replicates, and all correlation coefficients were � 0.82. Using the DESeq (2010) R Package with a FDR < 0.05 and a log 2 FC � 1 or � -1, pairwise comparisons of Ogura-CMS versus MF showed that 5,117 genes were significantly differentially expressed, of which 1,339 genes were significantly up-regulated and 3,778 genes significantly down-regulated in the Ogura-CMS line relative to the MF line (Fig 7). Representative genes for the upand down-regulated DEGs are listed in Tables 1 and 2, respectively, according to their

Functional annotation of differentially expressed genes during reproductive development of the Ogura-CMS line and its maintainer fertile line of turnip
To perform annotation analysis of the DEGs, eight public databases including COG, GO, KEGG, Swiss-prot database, KOG, Pfam, eggNOG, and Nr were searched. In total, 4,864 DEGs were found and annotated in detail in at least one of these databases (S4 Table). To forecast functional classifications of annotated DEGs, the GO and KEGG pathway analyses were performed to provide a clue. We utilized all DEGs for GO analysis and found that 76% of Transcriptome of Ogura-CMS turnip DEGs (3,889 out of 5,117) have at least one GO term assigned and were categorized into 46 functional groups ( Fig 8A). Among these, the top three dominant categories were involved in cellular, metabolic, and single-organism processes ( Fig 8A). The KEGG pathway analysis manifested that 50 pathways were significantly enriched, particularly metabolic pathways, plantpathogen interaction, and plant hormone signal transduction ( Fig 8B). The DEGs were found to be mostly enriched in ether lipid metabolism ( Fig 8C). Furthermore, many genes involved in fatty acid metabolism which is essential for the assembly of exine and tryphine [3,31] were dysregulated ( Fig 8D).
In all annotated DEGs, 1,289 genes were significantly up-regulated and 3,575 genes were significantly down-regulated in the Ogura-CMS relative to MF inflorescences (S4 Table). Among these, 610 DEGs were specifically expressed in the MF inflorescences and only 31 DEGs in the Ogura-CMS inflorescences (S5 and S6 Tables), implying that there is a considerable scope for further research to discover novel CMS-associated genes in turnip. In addition, several genes that were classified as function unknown but specifically expressed in the MF inflorescences, such as TRINITY_DN22922_c0_g2, TRINITY_DN37525_c0_g1, and TRINI-TY_DN23102_c0_g1, could be good candidates for CMS-related genes.

Genes related to anther development and microspore formation
Based on genetic and transcriptomic studies, it has been long assumed that Arabidopsis pollen development involves precise spatial and temporal cooperation of the tapetum and the gametophyte itself, and relies on the functions of numerous genes and their dynamic regulatory network [32]. We compared the expressive alteration of homologs of Arabidopsis genes previously reported to be associated with anther and pollen development, to unravel whether Transcriptome of Ogura-CMS turnip cytoplasmic retro-regulated counterparts of those genes from the nucleus. In addition, some functionally known genes involved in this unique process in species other than Arabidopsis were also compared. As the innermost layer surrounding the sporogenous cells in the anther, the tapetum provides not only energy, but also nutrients, metabolites, and sporopollenin precursors for microspore development [33]. Coordinated with the defective tapetum, homologs of some extensively demonstrated genes and enzymes associated with tapetum development in Arabidopsis (Table 3) and other species (Table 4) exhibited altered expression in Ogura-CMS. For example, AMS is a basic helix-loop-helix (bHLH) transcription factor, one of the master regulators for tapetum and microspore development in Arabidopsis [3]. Expression of a counterpart of AMS (TRINITY_DN27860_c1_g1) was down-regulated in Ogura-CMS. Some turnip homologs of AMS-dependent genes including QRT3 (TRINITY_DN22468_c2_g3), CYP98A8 (TRINITY_DN26854_c1_g2), CHS (TRINITY_DN27063_c0_g1), EXL6 (TRINI-TY_DN24807_c1_g1, and TRINITY_DN26730_c0_g3), and PAB5 (TRINITY_DN22996_ c2_g9, TRINITY_DN22996_c2_g5, and TRINITY_DN22996_c2_g10), showed reduced expression, but some genes did not, such as CYP703A2, CYP704B1, KCS7, LAP5 and LTP12 (Table 3). It was noteworthy that counterparts of two AMS-dependent genes, TKPR1 (TRINI-TY_DN23571_c0_g2, TRINITY_DN25196_c0_g1, and TRINITY_DN23571_c0_g1) and A6 (TRINITY_DN25127_c0_g3), suggested to be directly regulated by AMS in Arabidopsis, had increased expression in Ogura-CMS (Table 3). In addition, AMS regulatory pathway oriented analyses showed that turnip homologs of ATA20 (TRINITY_DN25943_c1_g3 and TRINI-TY_DN25943_c1_g1), bHLH89 (TRINITY_DN24940_c2_g5), and bHLH91 (TRINITY_ Table 2. Functional categories of representative genes significantly down-regulated in inflorescences of the Ogura-CMS line relative to its maintainer fertile (MF) line.

Gene ID log 2 fold change (Ogura-CMS line/MF line) Description
Transcription factors Transcriptome of Ogura-CMS turnip DN22984_c1_g2 and TRINITY_DN22984_c1_g1) were down-regulated (Table 3), but DYT1, TDF1, MS188, and MS1 were not in the pool of DEGs. Based on our cytological results in pollen wall development of turnip, a dramatically altered expression of numerous counterparts of known genes associated with pollen wall formation in CMS florets (Tables 3 and 4) also seemed to be responsible for obvious defects in male gametophyte development. Accordingly, we suggest a somewhat similar regulatory network underlying microspore development in Arabidopsis and turnip. Transcriptome of Ogura-CMS turnip Apart from the functionally known genes presented in Tables 3 and 4, 185 other novel DEGs were addressed to participate in male organ development in our GO analyses (S7 Table). Among these DEGs, 26 were specifically expressed in the MF inflorescences and none in the Ogura-CMS inflorescences (Table 5). Of the down-regulated genes, two were predicted to encode uncharacterized protein and four were functionally unknown, implying that they could be good candidates for further research to discover novel anther and microspore development-associated genes in turnip. Overall, these results indicated that mutation of cytoplasmic ORF138 retro-regulated genes from the nucleus, and interactions between them led to male sterility in Ogura-CMS turnip.

Morphological characteristics of Ogura-CMS turnip
Ogura-CMS originated from Japanese radish and has been extensively applied in hybrid breeding of crop species [1]. The morphological changes occurring during anther Table 4. List of known anther and microspore development-related genes in turnip and species other than Arabidopsis. Here, we emphasized the value of detecting variation in apparently uniform turnip. In Ogura-CMS flowers of turnip, the length of the stamens was greatly reduced, with shortened filaments and thin withered white anthers (Figs 1 and 2). It has been proposed that filament elongation concerns the polar auxin transport from pollen grains/anther to filaments, and hormone-biosynthetic enzymes are also essential for stamen development, including filament elongation [79][80][81]. When RNA-Seq results were classified by function, we assumed that altered expression of the auxin-, Transcriptome of Ogura-CMS turnip jasmonate-, and gibberellin-related genes (Fig 8 and S4 Table) may be the reason for short Ogura-CMS filaments. Further investigation showed that microspores exhibited aberrant development in Ogura-CMS anthers and this was initially observed at the uninucleate stage, leading to remnants of the aborted microspores with incomplete-developed layers of exine and tryphine (Fig 4). Defects in Ogura-CMS pollen wall formation were consistent with the observation that apart from the abortion of microspore development, the tapetum in Ogura-CMS anthers was swollen and abnormally vacuolated, initially at the tetrad stage (Figs 3 and 5). Tapetal organelles, such as elaioplasts and tapetosomes, were abnormally formed. In Brassicaceae plants, such as Arabidopsis, elaioplasts and tapetosomes are the carriers in the traffic of lipid molecules from the tapetum to the microspore [82]. Proper lipid accumulation in these organelles is critical for normal tapetum development and pollen wall formation, as disruptions of lipid metabolismrelated genes affect tapetum morphology and pollen wall patterning [31]. These results showed that defective microspore production and early tapetum degeneration during microgametogenesis led to complete male sterility of the Ogura-CMS line.

Anther development/Tapetum development and degeneration
It was also noteworthy that in sterile anthers, the middle layer persisted for a long time ( Fig  5). The question is still poorly documented for that if other non-tapetal anther layers, such as the middle layer, endothecium, and epidermis, are also essential for microspore development. Previously, a delay in PCD of the middle layer was reported to induce male sterility in kiwifruit (Actinidia deliciosa) [83]. The absence of well-defined genes in which mutation leads to delayed degeneration of the middle layer has enriched our knowledge on microspore development. For instance, Oryza sative tdr mutant shows a delay in middle layer degeneration as well as tapetum PCD [73]. In our RNA-seq data, the turnip homologue (TRINITY_DN27860_ c1_g1) of this bHLH transcription factor-encoding gene showed a significant decrease in transcript level in the Ogura-CMS line relative to the MF line (Table 4). Mutation of a leucine-rich repeat receptor protein kinase EXS/EMS1 in Arabidopsis causes abnormal persistence of the middle layer [84,85]. It is not clear what genes regulate the middle layer degeneration. Although non-tapetal anther layers argue against an essential role in pollen development, the molecular basis and effort of non-tapetal anther layers in microspore development still need further exploration in the future.
Morphology of all other floral organs was normal, suggesting that genes involved in floral organ identity were normal. Apart from affecting stamen development, several negative effects other than heterosis are resulted from the introduction of Ogura-CMS from radish, via B. napus into B. rapa [17]. It has been reported that the CMS mitochondrial signaling can also cause aberrations in other aspects including vegetative growth [2]. Though 10 generations of back-crossing, protoplast fusion eliminated the negative effects on crown diameter and plant height of Ogura-CMS plants, the fleshy root still had reduced size (Fig 1 and S2 Fig).

Cytoplasmic retro-regulates nuclear gene expression
Based on findings on Arabidopsis microsporogenesis [86], the early stage of this process in the Ogura-CMS line should be normal (Fig 3). Instead, genes associated with tapetal development and meiotic tapetal function are defective in the Ogura-CMS line. Moreover, the expression of genes associated with microspore development could be altered. It has been reported in CMS rice, carrots, wheat, and B. napus that mitochondrial mutations affect nuclear gene expression, through a communication pathway from mitochondria to the nucleus called the mitochondrial retrograde signaling pathway [4][5][6][7]. In addition, variation in anther behavior among Brassica species suggests the presence of different regulatory mechanisms and/or multiple regulatory pathways [17]. Therefore, we explored the differential gene expression patterns of transcripts between the Ogura-CMS and its MF inflorescences in turnip. RNA-Seq identified prodigious DEGs (S4 Table), supporting an extensive involvement of mRNAs in anther development. The complex transcript network was simplified by modularization. A lower number of dysfunctional DEGs in the Ogura-CMS inflorescences may also reflect the underdevelopment of stamens. Overall, these results indicate that the regulatory network of anther development is altered in Ogura-CMS turnip, the nuclear genes involved in this complex process are retro-regulated by mitochondrial orf138, and mitochondrial-to-nucleus signaling in Ogura-CMS plants perhaps causes aberrant male sporophyte and gametophyte development.

Genes involved in tapetum development and microspore formation show dysregulated expression
Pollen development involves delicate and coordinated spatio-temporal regulation of numerous biosynthetic genes by specific transcriptional regulators, and requires coordinated activity of different cell types and tissues of both gametophytic and sporophytic origin [3,32,87,88]. It is estimated that about 14,000 genes give a positive expression signal during microspore development, of which only a few have been shown to be pollen-specific with a defined requirement at different stages of pollen development [32]. In the CMS system, it may be more complicated because of the presence of the mitochondrial retrograde signaling pathway and the interaction between the nucleus and mitochondria [89]. We performed a transcriptomic comparison between the expression of some well-known anther and tapetum development-associated genes and their homologs in turnip. The results showed that cytoplasmic retro-regulated counterparts of some genes from the nucleus (Tables 3 and 4), such as GPAT6, which is essential for tapetum development, when disrupted causes defective exine deposition and partial pollen degradation [39]; RBOHE, encoding a NADPH oxidase and required for tapetal PCD, if interfered results in aborted male gametophytes [34]; and BAM1, a Leu-rich repeat receptor-like kinase encoding gene, is an important regulator for anther development [36].
During early tapetum development in Arabidopsis, EXS/EMS1 is a trigger to the signaling pathway for tapetal fate determination [84,85]. Here, the counterparts of EXS/EMS1 were not distinguishable in turnip (S3 and S4 Tables), indicating that retrograde signaling did not affect gene expression at the early stage of tapetum development and function, consistent with normal early tapetum behavior we observed (Fig 5). Later, core genetic networks involving several transcription factors regulate tapetum and microspore development, including a bHLH transcription factor Arabidopsis AMS/rice TDR, one of the master regulators [3]. AMS triggers at least 23 pollen wall development-related genes, such as A6, TKPR1 [33]. Mutation in AMS causes tapetal hypertrophy extending into the locule and subsequently leads to abnormality of microspore development [33,35]. The turnip ortholog of AMS was found to be down-regulated in the Ogura-CMS inflorescence (Table 3). Some AMS-dependent genes were also dysregulated in the Ogura-CMS inflorescence, but others, such as CYP703A2, CYP703B1, WBC27, TEK, and LAP5, were not. In addition, a transcriptional regulatory pathway involving AMS with DYT1-TDF1-AMS-MS188-MS1 has been proposed for tapetal and microspore development in Arabidopsis [3]. DYT1 is a bHLH transcriptional factor and TDF1 is a putative R2R3 MYB transcriptional factor, both of which are responsible for tapetum development and microspore formation [90,91]. However, the turnip orthologs of DYT1 and TDF1 were absent in our RNA-Seq data. Instead, the counterparts of the other two Arabidopsis bHLH transcriptional factors, bHLH89 and bHLH91 [37], and another R2R3 MYB transcription factor, MYB32 [48], had decreased expression in the Ogura-CMS inflorescence. It has been suggested that BrbHLH89 might replace DYT1 function and transcription factor genes BrTDF1, BrAMS, BrMS188, and BrMS1 in Chinese cabbage, however, it exhibited delayed or extended expression in the Ogura-CMS line [17]. Therefore, we suggest that the function of AMS may be retained but the regulatory AMS pathway may have some alteration between turnip and Arabidopsis, or other Brassica spp.
Recently, based on research on the inflorescence transcriptome of Ogura-CMS cabbage, 22 DEGs were assigned to pollen development [18]. Among these DEGs, some functionally known key members involving pollen development, such as BoLAP5, BoLAP6, BoACOS5, BoMS2, BoCYP703A2, and BoCYP704B1, were found to be down-regulated in Ogura-CMS anthers compared to the fertile line. In addition, some new members, such as BobHLH1, BoMYB1, and BoMF2, were also identified [92][93][94]. However, some of these well-defined genes were not part of the pool of our DEGs. These differences in gene expression may be the cause of different anther behaviors that Ogura-CMS cabbage has an abnormally fat tapetum with a delayed tapetum PCD, compared to our Ogura-CMS turnip which exhibited premature tapetum degeneration during microgametogenesis. These results provide new clues for a better understanding of different Ogura-CMS Brassica spp regulatory mechanisms.
Apart from those well-defined genes, 185 other novel DEGs were predicted to be involved in male organ development based on GO analyses, of which 26 DEGs were genotype-specific (S7 Table and Table 5), implying that they are good candidates for further research on novel anther and microspore development-associated genes in turnip. It is interesting that some members were previously reported to function in other developmental processes, not just anther development. For example, putative indole-3-acetic acid-amido synthetase GH3.9, a member of the GH3 family of auxin-responsive genes, was previously reported to contribute to primary root length [95]; GAMETE EXPRESSED 2 encodes a trans-membrance domain containing protein which is targeted by DUO1 in Arabidopsis [96]. These genes may be endowed with new putative roles in male organ development.

Genes related to pollen wall development
The indispensable biological function of pollen is somehow reflected by its unique surrounding wall. Pollen wall synthesis is regulated by both the sporophytic tapetum and the microspore [97]. As the innermost layer of the pollen wall, intine not only plays an important role in aperture configuration and pollen stability, but also functions in pollen tube growth into the stigma [77,87,98,99]. Intine comprises of pectin, cellulose, hemicellulose, hydrolytic enzymes and hydrophobic proteins, and is controlled by the gametophyte [97,100]. Several genes and proteins have been identified in intine formation [66,74,75,97,[101][102][103][104][105][106][107]. Here, we found that the expression of homologs of RGP1, RGP2, and PLL9 in the Ogura-CMS inflorescences of turnip was reduced (Tables 3 and 4), disruptions of which contribute to the premature degradation of developing microspores [66,75]. Arabinogalactan proteins (AGPs) are extensively glycosylated hydroxyproline-rich glycoproteins and are also proved to be indispensable for intine formation [78,108]. Currently, two classical AGP-encoding genes AGP6 and AGP11 and, their B. campestris ortholog BcMF18 and BcMF8, and a fasciclin-like AGP (FLA) gene FLA3, have been isolated and identified [77,78,109,110]. The orthologs of all these AGPencoding genes were found in turnip and showed down-regulated expression in the Ogura-CMS inflorescences (Tables 3 and 4). In addition, counterparts of another two FLA genes, FLA5 (TRINITY_DN22476_c2_g15 and TRINITY_DN22476_c2_g16) and FLA14 (TRINITY_ DN23986_c3_g2, TRINITY_DN23986_c3_g5, and TRINITY_DN23986_c3_g6), showed down-regulation in the Ogura-CMS inflorescences (S4 Table). Similarly, in Chinese cabbage, BrFLA14 was highly expressed in fertile buds and could be associated with pollen development [72], implying that this FLA gene could be important for intine formation in turnip.
In addition to the extensively reviewed genes and enzymes associated with pollen intine formation summarized above, the expression of some genes acting at different steps of exine formation was also dysregulated (Table 3). Exine formation requires proper pollen wall pattern determination and sporopollenin deposition. Pollen wall pattern determination is dependent on the formation and dissolution of the callose wall, which provides a structural basis for the primexine deposition [3]. Previous attempts to explore the roles of callose wall-related genes have indicated that overexpression of β-1,3-glucanase in the tapetum results in defective exine pattern formation [111]. Brassica napus and Arabidopsis A6 and its rice ortholog Osg1 encode β-1,3-glucanase that are essential for timely callose degradation [70,71]. As for the abnormality of tapetum, their counterpart in turnip (TRINITY_DN25127_c0_g3) was found to be up-regulated in the Ogura-CMS inflorescences (Table 4). In addition, the counterparts of CalS5 (TRI-NITY_DN23912_c2_g2 and TRINITY_DN24310_c1_g1), in which mutation causes abnormal exine wall formation [49], had a decreased expression in the Ogura-CMS inflorescences ( Table 3). Although the appearance of the Ogura-CMS callose wall at the tetrad stage was normal, the dysregulated expression of homologous genes of A6/Osg1 and CalS5 was suggested to be positively correlated with efficient microspore development. Sporopollenin deposition is dependent on the synthesis, secretion and translocation of sporopollenin precursors, which is closely associated with fatty acid metabolism [3,31]. In our RNA-Seq data, many genes involved in fatty acid metabolism pathway were dysregulated in the Ogura-CMS inflorescences (Fig 8D), indicating that they might also be responsible for exine defects. Notably, a number of functionally known key members of the complex biochemical biosynthesis of sporopollenin precursors in the tapetum were not in the pool of DEGs. ACOS5, a medium-and long-chain fatty acyl CoA-synthetase [112], are condensed with fatty acyl ACP reductase MS2 [113], cytochrome P450 family members CYP703As and CYP704Bs [114][115][116], flavonoid synthases LAP5 (PKSA) and LAP6 (PKSB) [117], together with tetraketide α-pyrone reductases TKPR1 and TKPR2 [57], yield precursors of sporopollenin building units. Disruption of these genes drastically hampers biosynthesis of sporopollenin precursors in the tapetum. Among all these genes, only the counterparts of TKPR1 were found up-regulated in our RNA-Seq results. We suggest that the reduced expression of ATP-binding cassette subfamily G (ABCG) transporters contributes to the transport of sporopollenin precursors across anther tissues [67]. Here, we only found ABCG1 had decreased expressed counterparts in the Ogura-CMS inflorescences (Table 3), which functions redundantly with ABCG16 in exine layer formation. These results suggest that the synthesis and secretion of sporopollenin precursors may be little affected by mitochondrial retrograde signaling. In addition, RNA-Seq analyses showed that a large number of lipid transfer proteins (LTPs)-encoding genes were dysregulated (S4 Table). Presently, no experimental evidence was obtained for the mechanism by which LTPs sporopollenin precursors are translocated, but LTPs are likely candidates for their delivery from the tapetum to the microspores [100]. Overall, these results indicated that during the process of sporopollenin deposition in Ogura-CMS turnip, the focal point of the effect of mitochondrial orf138 may be sporopollenin precursor translocation.
Genes related to pollen coat formation, which is also closely associated with fatty acid metabolism, are also indispensable in Ougra-CMS turnip. FLP1/CER3, is a member of the ECERIFERUM family, which appears to participate in the synthesis of wax, components of tryphine and sporopollenin of exine in Arabidopsis [61]. In the flp-1 mutant, excessive tryphine fills the cavity of the exine leading to a smooth surface that is sensitive to acetolysis. We found the homolog of FLP1/CER3 (TRINITY_DN27258_c1_g1) was down-regulated (Table 3). In addition, the expression of the counterparts of ABCG9 (TRINITY_DN25053_c4_g1), ABCG31 (TRINITY_DN27312_c0_g2 and TRINITY_DN27312_c0_g1) and GRP17 (TRINI-TY_DN22659_c0_g1, TRINITY_DN23861_c0_g2, and TRINITY_DN22389_c1_g4), which contribute to the accumulation of pollen coat [62,63], was also decreased ( Table 3). The altered expression of genes associated with the pollen coat may be related to male sterility in Ogura-CMS turnip.

Conclusions
In this study, detailed morphological characteristics of an Ogura-CMS line and its MF line of turnip were described. We employed high-throughput sequencing approaches to identify candidate nuclear genes that responded to mitochondrial retrograde signaling, and focused on the anther and microspore development-associated genes. Comparison of the gene expression between CMS and MF lines supported the notion that CMS is attributed to the mitochondrial retrograde signaling pathway and the interaction of nuclear and organelle genomes. Furthermore, our results suggest that different regulatory mechanisms and/or multiple regulatory pathways may exist in Brassica spp.