Transcriptomic and Proteomic Analysis of Shaan2A Cytoplasmic Male Sterility and Its Maintainer Line in Brassica napus

Cytoplasmic male sterility (CMS) lines are widely used for hybrid production in Brassica napus. The Shaan2A CMS system is one of the most important in China and has been used for decades; however, the male sterility mechanism underlying Shaan2A CMS remains unknown. Here, we performed transcriptomic and proteomic analysis, combined with additional morphological observation, in the Shaan2A CMS. Sporogenous cells, endothecium, middle layer, and tapetum could not be clearly distinguished in Shaan2A anthers. Furthermore, Shaan2A anther chloroplasts contained fewer starch grains than those in Shaan2B (a near-isogenic line of Shaan2A), and the lamella structure of chloroplasts in Shaan2A anther wall cells was obviously aberrant. Transcriptomic analysis revealed differentially expressed genes (DEGs) mainly related to carbon metabolism, lipid and flavonoid metabolism, and the mitochondrial electron transport/ATP synthesis pathway. Proteomic results showed that differentially expressed proteins were mainly associated with carbohydrate metabolism, energy metabolism, and genetic information processing pathways. Importantly, nine gene ontology categories associated with anther and pollen development were enriched among down-regulated DEGs at the young bud (YB) stage, including microsporogenesis, sporopollenin biosynthetic process, and tapetal layer development. Additionally, 464 down-regulated transcription factor (TF) genes were identified at the YB stage, including some related to early anther differentiation such as SPOROCYTELESS (SPL, also named NOZZLE, NZZ), DYSFUNCTIONAL TAPETUM 1 (DYT1), MYB80 (formerly named MYB103), and ABORTED MICROSPORES (AMS). These results suggested that the sterility gene in the Shaan2A mitochondrion might suppress expression of these TF genes in the nucleus, affecting early anther development. Finally, we constructed an interaction network of candidate proteins based on integrative analysis. The present study provides new insights into the molecular mechanism of Shaan2A CMS in B. napus.


INTRODUCTION
Rapeseed (Brassica napus) is an important oil crop producing both edible oil and industrial materials such as lubricants and biodiesel. Heterosis is widely observed in B. napus, with excellent hybrids having yields over 30% higher than their parents (Radoev et al., 2008). Hybrid breeding largely relies on male-sterile lines, mainly involving chemically induced male sterility (CIMS), genic male sterility (GMS), and cytoplasmic male sterility (CMS) (Chen and Liu, 2014). In higher plants, CMS is a maternally inherited trait, which is largely resulting from rearrangements of mitochondrial DNA, which results in an inability to generate pollen or in abnormal pollen (Schnable and Wise, 1998;Horn et al., 2014). The first CMS line was found in onion (Jones and Clarke, 1943). To date, CMS has been observed in more than 150 plant species (Laser and Larsten, 1972;Schnable and Wise, 1998). Four major CMS systems have been used in rapeseed production: nap CMS (Shiga and Baba, 1971), pol CMS (Fu, 1981), Ogu CMS (Ogura, 1968), and Shaan2A CMS (Li, 1980). pol CMS and Shaan2A CMS are the most commonly used CMS systems in B. napus in terms of number of three-line hybrids and the area planted with these hybrids (Fu, 1995).
In recent years, transcriptomic and proteomic analyses have been used in CMS studies. The transcriptome of CMS onion revealed three nuclear-related genes, AGAMOUS (AG), ABORTED MICROSPORES (AMS), and SOMATIC EMBRYOGENESIS RECEPTOR-LIKE KINASE 1 (SERK1) (Yuan et al., 2018). Transcriptomic analysis in Gossypium hirsutum indicated decreased ability to eliminate reactive oxygen species (ROS) in the CMS line, with ROS released from mitochondria acting as signal molecules in the nucleus, triggering formation of abnormal tapetum (Suzuki et al., 2013;Yang et al., 2018). In Ogu CMS, key genes participating in the secretion and translocation of sporopollenin precursors were significantly down-regulated in the cabbage R2P2CMS line . Differentially expressed genes (DEGs) involved in protein synthesis and metabolic pathways have been identified between sterile and maintainer lines in both pol CMS and SaNa-1A CMS in B. napus Du et al., 2016).
In flowering plants, anthers are important organs that generate pollen grains for propagation. In Arabidopsis, anther development has been systematically divided into 14 stages (Sander et al., 1999). First, the stamen primordium arises from the floral apex, then the archesporial cells from the four corners of the L2 cell layer undergo a series of differentiation and division events to promote formation of the four microsporangia of the butterfly shaped anther (Feng and Dickinson, 2007). These processes are controlled by a number of genes that function downstream of floral identity genes, most of them unknown. In Arabidopsis, AG has dual roles in limiting stem cell proliferation and determining floral organ identities (Lohmann et al., 2001). SPL, encoding a MADs transcription factor (TF), is a direct downstream gene of AG and the candidate gene for specifying the identity, position, and number of archesporial cells in the anther L2 layer; microsporogenous cells are not generated in the spl mutant (Schiefthaler et al., 1999;Yang et al., 1999). Another gene, HYPONASTIC LEAVES1 (HYL1), modulates the stamen architecture of the four microsporangia by coordinating expression of REVOLUTA (REV) and FILAMENTOUS FLOWER (FIL) through miR165/166 in Arabidopsis (Lian et al., 2013). Interestingly, FIL directly interacts with SPL (Sieber, 2004). Recently, GROWTH-REGULATING FACTOR (GRF) and GRF-INTERACTING FACTOR (GIF) were identified as novel positive regulators in specifying archesporial cells in Arabidopsis, and grf1/2 and gif1/2 mutants show a similar phenotype to the spl mutant (Lee et al., 2018). BARELY ANY MERISTEM (BAM1 and BAM2) encode CLAVATA1-related and leucine-rich repeat (LRR) receptor-like kinases, and anthers in the double mutant bam1 bam2 exhibit abnormalities at a very early stage, lacking the tapetum, middle layers, and endothecium (Hord, 2006). EXCESS MICROSPOROCYTES 1 (EMS1, also known as EXTRA SPOROGENOUS CELLS, EXS) also encodes a putative LRR receptor kinase, regulating tapetum identity and number of reproductive cells in Arabidopsis (Canales et al., 2002;Zhao et al., 2002). Interestingly, SPL is expressed in most early anther development mutants, for example, bam1/2 double mutants and ems1 (Canales et al., 2002;Zhao et al., 2002;Hord, 2006), indicating that SPL might be the first "reproductive gene" to be activated in anther development.
Shaan2A CMS was identified by Professor Dianrong Li in 1976, and the representative hybrid rapeseed cultivar 'Qinyou 2' was generated using the Shaan2A CMS line and its restorer line KC01. This was the first three-line hybrid cultivar in China, and was widely cultivated from 1985 to 2008 (Li and Tian, 2015). Although the Shaan2A CMS system has been widely used in the production of hybrids for decades, its male sterility mechanism remains unclear. In the present study, we performed transcriptomic and proteomic analyses, combined with additional morphological observation, to reveal the mechanism of Shaan2A CMS. We aimed to identify differences between the sterile line Shaan2A and its maintainer line Shaan2B at the transcriptional and protein level, and elucidate the regulative and metabolic pathways involved in the male sterility. The results will provide new insights into the molecular mechanism of Shaan2A CMS in B. napus.

Plant Materials
Male sterile plants were found among hybrid offspring by Professor Dianrong Li in 1976. The male sterile lines were conserved by self-fertilization with limited pollen, hybridization of sister plants, and test crossing of different rapeseed cultivars. The male sterile Shaan2A line and its maintainer line Shaan2B were then rapidly bred by individual pollination of a single plant, biparental crossing, and mutual selection of parents according to the female parent's phenotype as well as sterility of the female and agronomic characteristics of the parents. In the present study, Shaan2A and Shaan2B were cultivated in the test field of Huazhong University of Science and Technology (Wuhan, Hubei, China) under the same conditions. Plant samples were quickly removed from buds of different length on ice and frozen in liquid nitrogen, then kept at −80 • C for extraction of total RNA and protein.

Observation of Paraffin Sections and Transmission Electron Microscopy (TEM)
Anthers of different length were vacuum-infiltrated and fixed with 50%FAA (v/v) and 2.5% (w/v) glutaraldehyde in 0.1 M phosphate buffer (pH 7.4). For paraffin section analysis, fixed materials were dehydrated through a graded series of ethanol (70,85,90,95, 2 × 100%) and then cleared twice with xylene for 2 h. The material was infiltrated with Paraffin wax and subsequently embedded in paraffin wax. Finally, sections of approximately 10 µm were obtained using a KD-1508A microtome (China) and observed under a LEICA DMLB microscope (Germany) with a CCD camera (DFC420 FX). Procedures for TEM observation followed a previous study (Yi et al., 2010) with minor modification (Ning et al., 2018).

Proteomic Analysis
Different-sized anthers were collected from flower buds and mixed together. Preparation of total protein, two-dimensional electrophoresis (2-DE), analysis of differentially expressed proteins (DEPs), and protein identification by MALDI-TOF-MS-MS were conducted as described previously (Gan et al., 2013;Ning et al., 2018).

Transcriptomic Analysis
B. napus anthers from buds <0.5, 0.5-1, 1-1.5, and 1.5-2.0 mm are at the pollen mother cell stage, meiosis stage, tetrad stage, and stage of microspore release from the tetrad, respectively (Zhou et al., 2012). Young buds (YB) shorter than 1 mm represented the stage before meiosis, and small anthers (SA) in 1-2 mm buds represented the tetrad stage to the microspore release stage. Total RNA extraction, assembly, and mapping of clean reads, and annotation of the transcriptome were conducted according to a previous report (Ning et al., 2018). For Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis 1 , KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways.

qRT-PCR
Total RNA was reverse-transcribed using ReverTra Ace R qPCR-RT Master Mix with gDNA Remover (TOYOBO) according to the protocol of manufacturer, with actin as an internal reference gene . qRT-PCR experiments were conducted and calculations performed as described previously (Miao et al., 2016). All primer sequences are given in Supplementary Table S1.

Interaction Analysis
Interaction analysis of DEGs and DEPs was based on the STRING database 2 , including known and predicted proteinprotein interactions. Blastx (v2.2.28) was used to align target gene sequences to reference protein sequences, and a network was constructed according to known interactions in Arabidopsis thaliana.

Morphological and Cytological Comparison Between Shaan2A and Shaan2B
The petals of Shaan2B flowers were significantly wider than those of Shaan2A flowers, and the anthers and filaments in Shaan2B were obviously longer than those in Shaan2A (Figures 1A,B). The anthers in Shaan2A appeared to stop developing at an early stage and hardly split, releasing pollen eventually, while the pistil in Shaan2B was normal. Cytological differences between the anthers of Shaan2A and Shaan2B are shown in Figure 2. The anthers in 1 mm buds of Shaan2B were developed normally, and it was easy to distinguish the sporogenous cells and the middle layer, endothecium, and tapetum ( Figure 2E). However, most of the anther locules in 1 mm buds of Shaan2A were differentiated abnormally, and the anther wall and sporogenous cells could not be clearly distinguished (Figure 2A). In 2-4 mm buds, the anthers of Shaan2B were well developed, while the anthers in Shaan2A were obviously aberrant (Figure 2).
We conducted TEM observation to further understand ultrastructural differences in the anther wall between Shaan 2A and Shaan2B. As shown in Figure 3A, epidermis cells in 1 mm buds of Shaan2A had few starch grains in chloroplasts and a large vacuole. Except for the large vacuole, Shaan2B epidermis cells had starch grains in most of the chloroplasts ( Figure 3D). Tapetum cells in 1 mm buds of Shaan2A had a few mitochondria (2.8 ± 0.98) (Figure 3G), while those in Shaan2B were enriched with mitochondria (19 ± 1.55) ( Figure 3J). The lamella structure of chloroplasts in both epidermis and tapetum cells were obviously abnormal in 2 mm buds of Shaan2A (Figures 3B,H). However, lamella structure in Shaan2B was normal ( Figure 3E), and the structure of tapetum cells in Shaan2B ( Figure 3K) were distant not alike with those in Shaan2A ( Figure 3H). Moreover, epidermis and tapetum cells in 3 mm buds of Shaan2A had the same structure ( Figures 3C,I), while those of Shaan2B exhibited different differentiation (Figures 3F,L). The different cell layers in anther walls during early anther differentiation are roughly shown in Supplementary Figure S1. These phenomena indicated that the differentiation of anther wall cells in Shaan2A was disordered.

Comparative Transcriptomic Analysis Between Shaan2A and Shaan2B
We used comparative transcriptomics to identify genes associated with male sterility in Shaan2A CMS. RNA-seq of three biological replicates at the YB (<1 mm buds, before meiosis) and SA (1-2 mm buds, representing tetrad to microspore release) stages generated 637,533,762 raw reads. The raw reads were deposited in the SRA (NCBI Sequence Read Archive) database under accession number PRJNA502996. These raw data were filtered to obtain clean reads, of which over 90% perfectly matched to the reference genome (Chalhoub et al., 2014). More detailed information for RNA-seq data is shown in Table 1. We performed hierarchical cluster analysis of all DEGs to examine global gene expression patterns in different samples ( Figure 4A). Shaan2A_YB and Shaan2B_YB showed similar global expression patterns. When compared with Shaan2B, a total of 3,509 and 9,636 up-regulated DEGs and 5,634 and 12,291 down-regulated DEGs (P < 0.05) were identified at the YB and SA stages in Shaan2A, respectively ( Figure 4B). The number of DEGs shared at the two stages was 4,328, only 1/5 of the total DEGs at the SA stage ( Figure 4C). The number of DEGs, and the densities of single-nucleotide polymorphisms (SNPs) and insertions and deletions (indels) in different chromosomal regions at the YB and SA stages, respectively, are shown in Figure 4D. High-density areas of both SNPs and indels were focused on chromosomes A01, A02, A03, and A05, and DEGs were more concentrated in genome A than genome C.
KEGG analysis (P < 0.05) of DEGs between Shaan2A and Shaan2B at the two stages revealed that pathways associated with protein processing in endoplasmic reticulum, carbon metabolism, and fatty acid degradation were enriched at the YB stage (Supplementary Table S2), as some genes in these pathways, such as HEAT SHOCK PROTEIN (Hsp70), SHEPHERD (SHD), and LONG CHAIN ACYL-COA SYNTHETASE 9 (LACS9) were significantly differentially expressed. In Arabidopsis, Hsp70 is a molecular chaperone that plays numerous important roles in protein folding (Wisén and Gestwicki, 2008). At the SA stage, genes associated with starch and sucrose metabolism, porphyrin and chlorophyll metabolism, and carbon metabolism pathways were annotated (Supplementary Table S2). Among these genes, FERROCHELATASE (FC) and CHLOROPHYLLASE (CLH) were identified. FC is the last enzyme for heme biosynthesis in the tetrapyrrole pathway and is essential for stress tolerance in Arabidopsis (Chow et al., 1998;Fan et al., 2018).
To obtain an overview of the pathways in which the DEGs participated, we further analyzed the DEGs (| log 2 ratio| > 1) at the YB and SA stages using MapMan software. The metabolism pathway visualization is shown in Figure 5, and detailed information is listed in Supplementary Table S3. Although the DEGs were dispersed in different primary and secondary metabolism pathways, most involved in lipid and flavonoid metabolism and mitochondrial electron transport/ATP synthesis pathways were down-regulated at  both the YB ( Figure 5A) and SA (Figure 5B) stages; for example, BnaA02g21480D, which encodes a member of the SMO1 family of sterol 4-alpha-methyl oxidases in the lipid metabolism pathway, and BnaA02g33800D, encoding flavonol synthase 3 (FLS3) in the flavonoid metabolism pathway. Interestingly, a large number of DEGs were up-regulated in light reactions, photorespiration, and the Calvin cycle ( Figures 5A,B). These included BnaA03g22080D, which encodes a PsbP domain-OEC23 like protein localized in thylakoids (peripheral-lumenal side) in the light reaction pathway, BnaA01g00170D, which encodes a protein with mitochondrial serine hydroxymethyl transferase activity and functions in the photorespiration pathway, and BnaC04g05700D, which encodes a Rubisco activase and functions in the Calvin cycle.  Furthermore, 11 DEGs related to anther development were selected for qRT-PCR validation (Supplementary Figure S2), for example, BnaA01g09760D (AG), BnaA10g23720D (EMS1), and BnaA01g16350 (SPL). Most of the expression patterns of these DEGs at the two stages were basically consistent with the RNA-seq data, indicating that the RNA-seq data were reliable.

Comparative Proteomic Analysis Between Shaan2A and Shaan2B
To detect DEPs between Shaan2A and Shaan2B, we conducted 2-DE with three biological triplicates using total protein from anthers of Shaan2A and Shaan2B. Approximately 1,000 protein spots were detected (Figures 6A,B). DEPs exceeding the least significant difference (p < 0.05) and showing more than twofold change in abundance were chosen for MALDITOF-MS-MS analysis.
To further facilitate the biological interpretation of the identified DEPs, we conducted hierarchical clustering analysis of the DEPs showing quantitative changes in expression. We identified four clusters of DEPs ( Figure 7A) with dramatic changes in expression patterns between fertile and sterile anthers. Cluster A contained two proteins, one belonging to the "carbohydrate metabolism" group and the other unknown (Figure 7B), and both were dramatically down-regulated in Shaan2A compared with Shaan2B. Cluster B, the largest cluster, included 51 proteins down-regulated in Shaan2A, with the "genetic information processing" and "carbohydrate metabolism" categories dominating ( Figure 7B); for example, acetyl-CoA C-acetyltransferase protein (spot 428; AT5G48230) and adenosine kinase 2 (spot 7412; AT5G03300). Cluster C proteins showed increased expression in Shaan2A compared with Shaan2B, and included 12 proteins with "energy metabolism" as the dominant classification ( Figure 7B); for instance, ATP synthase beta subunit (spot 6523; ATCG00480) and ribulose-1,5bisphosphate carboxylase/oxygenase large subunit (spot 2502; ATCG00490). Cluster D contained only one protein with high expression in Shaan2A and low expression in Shaan2B ( Figure 7B). This was a V-type proton-ATPase protein related to "energy metabolism."

Correlation Analysis Between DEGs and DEPs
B. napus is an allotetraploid plant, and its genome contains many repeats and homologous sequences (Chalhoub et al., 2014). The 66 DEPs identified in the proteomic analysis corresponded to 344 genes in B. napus. We conducted correlation analysis between DEGs and DEPs at both the YB and SA stage. DEPs and DEGs were divided into four groups according to their expression in Shaan2A compared to Shaan2B (Supplementary Table S4). In Group I, 43 and 115 DEGs showed the same trend in expression as their corresponding DEPs at the YB and SA stage, respectively. For example, actin 3 (AT3G53750) and chalcone and stilbene synthase family protein (AT4G00040), which are associated with protein processing in the endoplasmic reticulum, and SAM. In Group II, 21 and 24 DEGs showed the opposite trend in expression from their corresponding DEPs at the YB and SA stage, respectively, such as MAT4 (AT3G17390) and GLN1.3 (AT3G17820), associated with amino acid and carbohydrate metabolism, respectively. In Group III, the corresponding genes of 280 and 205 DEPs at the YB and SA stage, respectively, were not differentially expressed. In Group IV, the genes corresponding to 33 DEPs at the YB stage and 10 DEPs at the SA stage did not show differential expression in the transcriptomic analysis. These results showed that both positive and negative correlations exist between the mRNA and protein expression profiles at different stages, indicating that there might be a complex post-transcriptional regulatory network in Shaan2A male sterility.

Gene Ontology (GO) Analysis of DEGs and DEPs
We further analyzed the up-and down-regulated DEGs (| log 2 ratio| > 1) at the YB and SA stages and total DEPs between Shaan2A and Shaan2B using Cytoscape BiNGO and visualized the functional enrichment with Cytoscape EnrichmentMap. As shown in Figure 8, and DEPs generated 102,59,300,117, and 85 nodes, respectively. These nodes were classified into different categories. Interestingly, the common GO terms among these stages were "response to stimulus" and "metabolic process, " which might largely be due to the toxic action of the sterility gene in mitochondria disrupting the balance of metabolic processes in Shaan2A. Very importantly, nine categories associated with anther and pollen development were enriched among down-regulated DEGs at the YB stage: "androecium development, " "anther development, " "microsporogenesis, " "pollen development, " "pollen exine formation, " "pollen wall assembly, " "sporopollenin biosynthetic process, " "stamen development, " and "tapetal layer development." This indicates that anther and pollen development were disordered at the YB stage in Shaan2A, consistent with the cytological observations mentioned above. Expression levels of the genes in these GO categories at the YB stage are given in Supplementary Table S5. For example, three copies of CYP703A2 were remarkably down-regulated (BnaA10g00600D, log 2 ratio = −13.53; BnaC05g00670D, log 2 FIGURE 8 | Biological process analysis of DEGs and DEPs. GO modules enriched by up-and down-regulated DEGs and DEPs visualized using Cytoscape EnrichmentMap. Yellow and red circles show different and common biological processes between up-and down-regulated DEGs, respectively. The related GO categories are circled with a black dashed ring. ratio = −13.30; BnaCnng08630D, log 2 ratio = −10.29); this gene is expressed in developing anthers, and knockout lines display a partial male sterile phenotype in Arabidopsis (Morant et al., 2007). Two copies of LESS ADHESIVE POLLEN 6 (LAP6) were also significantly down-regulated (BnaA04g21720D, log 2 ratio = −9.50; BnaC04g45570D, log 2 ratio = −12.63); this gene is transiently and specifically expressed in tapetal cells during pollen development in Arabidopsis, and the mutants show defective pollen exine (Kim et al., 2011).

Transcription Factor Analysis
Transcription factors play key roles in multiple biological processes. We identified 349 up-regulated and 464 downregulated TFs in Shaan2A compared to Shaan2B at the YB stage, and 867 up-regulated and 596 down-regulated TFs at the SA stage (P < 0.05, Supplementary Table S6). Most genes belonging to "bHLH, " "C2H2, " "MADS, " and "MYB" TF families were downregulated in Shaan2A at the YB stage. Interestingly, 110 TFs were only expressed in Shaan2B at this stage. These included two copies of DYT1 (a putative bHLH TF), predicted to play a key role in a putative regulatory network model for Arabidopsis anther development, and three copies of MYB80, the downstream gene of DYT1 (Feng et al., 2012;Cui et al., 2016). Further, some TFs reported to be key regulators of anther differentiation were down-regulated at the YB stage.

Genes Related to Early Anther Differentiation
During early anther differentiation, SPL, which is activated by AG, together with its down-stream genes, specifies archesporial cells in Arabidopsis before meiosis (Feng and Dickinson, 2007). The expression levels of 23 of these genes at the YB stage are shown in Supplementary Table S7, with 12 of them significantly differentially expressed between Shaan2A and Shaan2B. Among these, three genes were slightly up-regulated in Shaan2A: ER, GRF, and ERL2. Interestingly, the eight down-regulated genes showed large log 2 ratios for expression level; for example, two copies of SPL (BnaA01g16350D, log 2 ratio = −2.72; BnaC01g19500D, log 2 ratio = −2.66), ROXY2 (BnaA02g01990D, log 2 ratio = −4.72), four copies of AMS (BnaA03g39400D, log 2 ratio = −13.00; BnaA07g03340D, log 2 ratio = −9.00; BnaC03g46740D, log 2 ratio = −10.33; BnaC07g05950D, log 2 ratio = −11.92), and two copies of MS1 were only expressed in Shaan2B. There were 11 genes with no significant difference in expression between Shaan2A and Shaan2B or no expression in either the sterile or maintainer lines, for instance, AG, SPL8, TPD1, and SERK1/2, indicating that they might not participate in Shaan2A CMS.

Interaction Analysis for Candidate Genes
Transcriptomic and proteomic analyses suggested that CMS in Shaan2A is controlled by a complex regulatory network. To further elucidate this mechanism, we investigated known and predicted interactions among candidate proteins corresponding to our DEGs and DEPs; for example, TFs associated with anther differentiation and genes belonging to the anther and pollen development GO categories mentioned above. Protein-protein interactions based on Arabidopsis orthologs were identified using STRING 10.0 and visualized using Cytoscape 3.6.1, and an interaction network associated with anther development was constructed. As shown in Figure 9, this network contained 236 nodes. Interestingly, TFs down-regulated at the YB stage showed interactions with many other proteins. For instance, SPL interacted with 26 targets including MYB65, TDF1, ROXY2, AMS, and MS1, DYT1 interacted with 19 targets including ABCG26, LAP3, LAP5, AT1G33430, and CYP704B1, and MYB65 interacted with seven targets including SPL, AMS, MS1, and DYT1. Therefore, these TFs and their target proteins might form a complex network for regulating anther differentiation and pollen development in Shaan2A CMS.

DISCUSSION
In recent years, CMS has been widely studied by transcriptomic and proteomic analysis in different species, including onion (Yuan et al., 2018), cotton (Suzuki et al., 2013;Yang et al., 2018), maize (Li et al., 2017), and rapeseed Du et al., 2016). These studies revealed that cotton CMS-D8 lines have reduced ability to eliminate ROS , secretion and translocation of sporopollenin precursors is significantly down-regulated in cabbage R2P2CMS , and significant DEGs in rapeseed SaNa-1A CMS and pol CMS mostly participate in protein synthesis and metabolic pathways Du et al., 2016). In the present study, DEGs were enriched in "protein processing in endoplasmic reticulum, " "carbon metabolism, " and "fatty acid degradation" categories at the YB stage, and "starch and sucrose metabolism, " "porphyrin and chlorophyll metabolism, " and "carbon metabolism" categories at the SA stage. Moreover, DEPs indicated that functional processes of "genetic information processing, " "carbohydrate metabolism, " "energy metabolism, " and "amino acid metabolism" were disordered in Shaan2A.

Shaan2A Is Defective in Early Anther Differentiation
Tapetum cells in a sterile onion CMS line begin to degrade from the tetrad stage and microspores start to abort, while degradation in its maintainer line is not observed until the microspore stage (Yuan et al., 2018). Fertile and sterile cotton anthers show differences at the meiotic stage, but no obvious variations before the stage when the microsporocytes and four anther wall cell layers are formed . In the C-type CMS of maize, no significant differences in microspores are observed between fertile and sterile plants until the mononuclear stage (Li et al., 2017). Although the pollen sac in the SaNa-1A CMS line forms normally, tapetum cells develop vacuolation at the microsporocyte and tetrad stages compared with the maintainer line (Du et al., 2016). In this study, the anther epidermis, endothecium, middle wall layers, tapetum, and sporogenous cells could not be clearly distinguished in Shaan2A by paraffin section observation (Figure 2). TEM analysis also revealed fewer mitochondria in the anther wall cells of Shaan2A than of Shaan2B, and the lamella structure of chloroplasts in Shaan2A anther wall cells was obviously abnormal (Figure 3). In contrast to Shaan2B, there was no difference between epidermis and tapetum cells in Shaan2A at the late stage of anther differentiation (Figure 3). Hence, we deduce that defective differentiation of anther wall cells and sporogenous cells in Shaan2A results in male sterility.

Shaan2A Anther Cells Have Damaged Mitochondria
Cytoplasmic male sterility is a maternally inherited trait resulting from the interaction of a nuclear fertility restoring gene (Rf ) and a mitochondrial CMS gene (Horn et al., 2014;Yamagishi and Bhat, 2014). A few CMS-associated genes have been reported, for example, orf224 in pol CMS (L'Homme and Brown, 1993;Menassa et al., 1999) and orf138 in Ogu CMS (Duroc et al., 2005). However, the sterility gene in the Shaan2A mitochondrial genome has not been identified. The mitochondrion is a crucial organelle for metabolic pathways such as respiratory electron transfer, ATP synthesis, and the TCA cycle (Reichert and Neupert, 2004;Logan, 2006). As well as the number of mitochondria in Shaan2A anthers being less than that in Shaan2B, most DEGs associated with the "mitochondrial electron transport/ATP synthesis pathway" and "TCA cycle" were obviously down-regulated at the YB and SA stages in Shaan2A anthers (Figure 5), including some important genes such as Bnac04g22440D (CITRATE SYNTHASE 5, log 2 ratio = −7.69 at YB stage), BnaC03g08060D (NAD(P)H DEHYDROGENASE B4, log 2 ratio = −5.55 at YB stage), and BnaA05g17260D (AOX1a, log 2 ratio = −2.34 at YB stage). Additionally, 15.2% DEPs were grouped into the "energy metabolism" functional category. It is interesting that AT4G11150 (a V-type proton-ATPase) FIGURE 9 | Interaction analysis of candidate proteins. Red nodes represent TFs related to early anther differentiation. Rose red nodes represent proteins in the nine GO categories of anther and pollen development. Light blue nodes represent proteins in the significantly enriched pathways identified by KEGG or MapMan analysis. and AT1G76030 (a nucleotide-binding subunit of vacuolar ATPase) were both up-regulated in Shaan2A anthers compared to Shaan2B anthers; if defective mitochondria in Shaan2A do not produce enough energy to meet the needs for rapid anther growth, these ATPases might be up-regulated to rescue the energy metabolism in the sterile anthers. As the sterility gene in Shaan2A is unknown, the mechanism by which it functions and induces a series of response in mitochondria will need further study.

Shaan2A Anther Cells Have Injured Chloroplasts
The abundance of starch in chloroplasts is necessary for microspore development and is a crucial feature of fertile pollen (Du et al., 2016). Most studies of male sterility in plants have identified "carbon metabolism" as disordered (Cheng et al., 2013;Li et al., , 2016Du et al., 2016;Liu et al., 2018;Yuan et al., 2018). Interestingly, the lamella structure of chloroplasts in Shaan2A anther wall cells was obviously abnormal. Most DEGs involved in "carbon metabolism, " "light reactions, " "photorespiration, " and "Calvin cycle" at the YB and SA stages (Figure 5 and Supplementary Table S3) and in "starch and sucrose metabolism" and "porphyrin and chlorophyll metabolism" at the SA stage were up-regulated (Supplementary  Table S2). However, proteomic analysis revealed that the 12 DEPs grouped into "carbohydrate metabolism" were mostly downregulated, except for ATCG00490, (ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit, encoded by a chloroplast gene) and AT3G62030 (a nuclear-encoded chloroplast stromal cyclophilin CYP20-3, also known as ROC4). Despite many carbon metabolism related pathways being identified by upregulated genes, chloroplasts of Shaan2A contained fewer starch grains than those in Shaan2B (Figure 3). Therefore, we speculate that high expression of many genes related to carbon metabolism is triggered by the injured chloroplast structure in Shaan2A. However, this up-regulation cannot balance the carbon metabolism in cells, which might be affected by the sterility gene. The direct or indirect mechanism by which the sterility gene impacts chloroplasts requires further study.

Genes Involved in Early Anther
Differentiation Are Down-Regulated at the YB Stage SPL functions downstream of AG, specifying the identity, position, and number of archesporial cells in the anther L2 layer (Schiefthaler et al., 1999;Yang et al., 1999). We identified two copies of SPL that were down-regulated in Shaan2A compared with Shaan2B, while three copies of AG were not significantly expressed (Supplementary Table S7). We therefore presume that SPL is the first gene repressed by the sterility gene in the mitochondrion. Furthermore, SPL was predicted to interact with 26 targets (Figure 9). These include TFs such as MYB65, ROXY2, AMS, MS1, and DYT1, which were also obviously downregulated in Shaan2A (Supplementary Table S7). In Arabidopsis, MYB33 and MYB65 share high sequence similarity. The myb33 myb65 double mutant has a deficiency in anther development, with hypertrophied tapetum at the pollen mother cell stage and aborted microspores before meiosis; the single mutant shows no phenotype, implying that MYB65 and MYB33 are functionally redundant (Millar, 2005). In B. napus, we only identified putative orthologs of AtMYB65; however, three copies of BnMYB65 were significantly down-regulated. ROXY1 and ROXY2 together control anther development in Arabidopsis, and single mutants produce normal anthers; however, double mutants are sterile owing to defects in early lobe differentiation (Xing and Zachgo, 2008). In the present study, two copies of ROXY2 showed significant down-regulation in Shaan2A compared with Shaan2B (Supplementary Table S7). In addition, four copies of AMS, two copies of MS1, and three copies of DYT1 exhibited large log 2 ratios of down-regulation in Shaan2A (Supplementary Table S7); all of these are located downstream of SPL and essential for anther and pollen development in Arabidopsis (Wilson et al., 2001;Feng and Dickinson, 2007;Feng et al., 2012;Ferguson et al., 2017).
SPL was also predicted to interact with ABCG26, MEE48, LAP5, LAP6, ACOS5, EMS1, and other proteins (Figure 9). In Arabidopsis, ABCG26 is an ATP binding cassette (ABC) transporter localized to the plasma membrane and endoplasmic reticulum that mediates trafficking of hydroxycinnamoyl spermidines and polyketides needed for sporopollenin formation (Quilichini et al., 2010(Quilichini et al., , 2014. Maternal Effect Embryo Arrest 48 (MEE48) is related to morphogenesis and development (Chakraborty et al., 2015). LAP5 and LAP6 catalyze condensation of malonyl-CoA units with CoA ester starter molecules to generate various natural products, and are required for sporopollenin biosynthesis and pollen development (Kim et al., 2011). The acyl-CoA synthetase 5 (ACOS5) gene encodes a fatty acyl-CoA synthetase that plays a key role in sporopollenin biosynthesis and exine formation (de Azevedo Souza et al., 2009). EMS1 encodes a putative LRR receptor protein kinase (LRR-RPK) controlling reproductive cell fates in anthers (Canales et al., 2002;Zhao et al., 2002). Interestingly, in B. napus, six copies of MEE48, seven copies of LAP5, two copies of LAP6, and two copies of ACOS5 showed high log 2 ratios of down-regulation in Shaan2A or were only expressed in Shaan2B, and two copies of EMS1 also displayed lower expression in Shaan2A than in Shaan2B (Supplementary Tables S5, S7). These genes were also predicted to interact with many other genes (Figure 9). Taken together, we propose that after the male sterility signal from the mitochondrion is transferred to the nucleus, it initiates suppression of SPL expression resulting in down-regulation of its downstream genes, which form a huge and complicated regulatory network for anther differentiation in Shaan2A. However, the mechanism underlying this network requires further research.

DATA AVAILABILITY
Transcriptomic raw reads were submitted to the SRA (Sequence Read Archive of NCBI) database with accession number PRJNA502996.

AUTHOR CONTRIBUTIONS
ML and HW conceived and designed the experiments. LN and ZL performed the experiments and analyzed the data. DL provided the Shaan2A and Shaan2B materials. YL, WZ, HC, and LM helped to analyze the data. LN and ML wrote the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We thank the Research Core Facilities for Life Science for qRT-PCR experiments.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00252/ full#supplementary-material Result from RNA-seq. Numbers are log 2 (X)-normalized ratio values. Red represents higher gene expression levels. Green corresponds to lower gene expression levels. Blocks without a numerical value indicate that gene expression was not significantly detected by RNA-seq.