Phenotypic Characteristics of C303A and 303B stamen
The wheat CMS line C303A was developed from stable sterile lines by consecutive backcrossing with 303B, as the donor parent, over 20 times in Yangling, China. To get the abortive morphological feature of C303A, we collected from C303A stamens in five developmental stages. Field comparative observation showed that except for floral organs, there was no significant diversity in growth and overall morphology between the CMS line C303A and its isomaintainer line (303B) (Fig. 1). At BNS, the anthers of 303B were quite full and yellowish, and which was bright yellow and the upper and lower ends were bifurcated, normal cracking accompanied by a large number of pollination phenomenon at the trinucleate stage. However, the anther color of C303A was white in the tetrad stage. Stamen malformations occurred in uninucleate stage and individual stamens with curved folds. Some stamens have become a combination of stamens and pistils at the binuclear stage. The stamens at the trinucleate stage have been completely transformed into pistils (Fig. 2). Moreover, to accurately observe the cytological structure of pistillody stamens, we used paraffin section. As shown in Figure 3, there were apparent differences between fertile stamen and pistillody stamen. Compared with the tetrad stage of 303B, the four anther locules in C303A sequentially shrunk and shriveled. Meanwhile, there was a marked degradation in one anther locule without microspores, another anther locules were empty and only with a very thin tapetum. However, a small number of microspores were detected in the other two anther locules. At early uninucleate stage, the degradation anther locules have been increased. Besides, the microspores shrank and condensed obviously in another anther locule. During the later uninucleate stage, the outlines of the tapetal cell and microspores were wholly invisible and some vascular bundles were found on degradation anther locule, thereby contributing to transport of nutrients to the locule for ovule formation. At the binucleate stage, the two degenerated chambers began to merge, and we speculate that this may be the beginning of the formation of ovules. Up to the trinucleate stage, the structure of the ovule in the pistillody stamen was determined in transverse and longitudinal sections. The above results indicate that the pistillody stamens contained ovule structures, instead of pollen grains and tapetum, thereby the stamens of C303A are unable to produce mature pollen grains, which could be detected by potassium iodide staining. So, we deduced that pistillody might occur at the binucleate stage.
Sequence Analysis Using RNA-Seq
To understand the basic molecular mechanism responsible for the pistillody at the transcriptional level, we conducted employing an Illumina HiSeq PE1500 sequencer for a transcriptome sequencing analysis of binucleate stage stamen of CMS line C303A and its maintainer line 303B. Stamens were collected three times, for a total of three biological replicates and sequencing reads were 150 bp in length. After filtering out the reads with>10% ambiguous nucleotide, adapter sequences and low-quality regions, 270,683,956 clean reads were produced, with 131,000,548 reads from the maintainer line, and 139,683,408 from the CMS line. Additionally, the GC content ranged from 55.80% to 59.63%，and the Q20 percentage surpassed 88.93%. Each sample’s clean reads were matched to Triticum aestivum reference sequence, where the alignment efficiency was between 60.46% and 68.09% (Table 1). The throughput and sequencing quality show that the RNA-Seq data we obtained were adequately high quality to ensure further analysis.
Identification of DEGs by RNA-Seq
In total, RNA-Seq detected 179,898 genes. For determining the significant difference in gene expression levels, we used the FDR < 0.05 and log 2 Fold Change (|log 2 FC|)>1) as the threshold. To compare the DGEs at the binucleate stage for C303A and 303B based on significant differences (Fig. 4). A total of 20,444 genes were discovered as being differentially expressed between C303A and 303B stamen. These DEGs included 10,283 up-regulated and 10,161 down-regulated genes in the C303A stamen compared with the corresponding expression levels in the 303B stamen.
Gene Ontology Enrichment Analyses of the DEGs
GO is a universal standardized gene functional classification system that gives a dynamically-updated controlled vocabulary and a rigidly defined concept to comprehensively describe properties of genes and their products in any organism . After enrichment analysis, the DEGs found in C303A have been annotated 36 functional groups (Fig. 5). Among the biological process functions, the central DGEs were associated with cellular processes, single-organism process, localization and metabolic processes function. With respect to the cellular component, the DGEs associated with cell, cell parts, membranes, and organelles. Besides, binding, catalytic activity and transporter activity are closely related to molecular function. There are twenty significantly enriched GO terms in the biological process functions which revealed by hypergeometric tests. Among twenty GO terms (Table S1), there are two terms the q-value equals zero: (GO: 0005975) carbohydrate metabolic process and (GO: 0044710) single-organism metabolic process. The DEGs in different functional categories may provide a valuable resource for the study of stamen development in C303A.
To identify biological pathways, the DEGs in the binucleate stage were mapped to 129 pathways in the KEGG database, where the top 36 pathways (Table S2) were deemed significant at a cut-off FDR corrected q-value < 0.05. And the main DEGs mainly comprised starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, glycolysis, phenylpropanoid biosynthesis, pyruvate metabolism, citrate cycle (TCA cycle), pentose and glucuronate interconversions and cyanoamino acid metabolism (Fig. 6). Most of the genes mapped in the nine significantly enriched pathways had down-regulated expression trends except for the circadian and phenylpropanoid biosynthesis. Indicating that starch and sucrose metabolism, glycolysis and TCA cycle could be essential for pollen development. Previous studies have proved that the TCA cycle has a pivotal role in plant male sterility  and the normal energy metabolism process of plants can satisfy their growth and development. So, we pay more attention to the significant genes related to these pathways, which may explain why no microspores in pistillody stamen.
Identification of MADS-box transcription factor involved in pistilldoy
In plant, the MADS-box gene family play essential role in the ABCDE model of flowering. According to this mole, 3 (TaAGL40, WM19B, and WM9B), 1 (WM27A), 2 (WM24A, TaAGL36) and 5 (TaAGL14, WM25, MAD-Box factor 2A, MAD-Box factor 2D, and MAD-Box factor 2B) wheat genes were identified to be class A, class D, SVP and class B, respectively (Fig. 7). Meanwhile, previous studies have identified several class B MADS-box genes in wheat, such as WPI2 and WPI1. These genes are mainly related to the transformation of wheat from vegetative growth to reproductive growth. These genes change their expression patterns in heterogeneous wheat and common wheat, and which contribute to the transformation of their stamens into pistil-like structures. The conserved domains and motif compositions of class B MAD-box genes were analyzed. As shown in Figure 8, a conserved K-box is present in the most class B MAD-box genes, except for MAD-box factor 2 genes. A distinct MADS superfamily was exclusively found in OsMADS32 and TaAGL14. From the analysis of motif compositions, there are 10 conservative motifs identified in class B MAD-Box, which are named motif1 to motif 10. Motif_1, as a conservative motif of MAD-Box, appears in all MAD-boxes of Class B. Motif_2 only appears in AP1 subfamily, motif_10 only appears in GGM13 subfamily and motif-10 only appears in MAD-box factor transcription factor 2. Based on the above results, we find that WM25 and OsMAD29, TaAGL14 and OsMAD32 have the same conserved domains and motif compositions, which they might have the same functionality. Meanwhile, the results of amino acid sequence alignment showed that the similarity was more than 92% (Fig. s1, Fig. s2).
A Possible molecular basis for the energy deficiency model in C303A
According to the analysis of the above results, we get three metabolic pathways involved in carbohydrate metabolism and energy metabolism (glycolysis, TCA cycle and pyruvate metabolism) to regulate pollen development. Combined with previous studies, we suggested a probable regulatory network leading to no microspores in C303A (Fig. 9). During anther and pollen development, the amount of pyruvate, the final product of glycolysis, might be decreased due to the down-regulated expression of the enzyme involved in glycolysis, which affects the TCA cycle to some extent. Also, the down-regulated expression of the enzyme genes associated with the TCA cycle contributed to a reduction in the number of coenzymes (FADH2 and NADH), whereby decreasing coenzymes entering the electron transport chain. In the electron transport chain, the down-regulation of antioxidant enzymes and key complexes represses electron transfer, and then electrons directly transfer to molecular oxygen, which produces excess ROS. Besides, upregulation of the activities of the antioxidant enzyme, which ruins the balance of the antioxidant defense system. Because ROS cannot be eliminated adequately immediately, which increases the consumption of ATP and the production of H2O2 and this might be finally causing no microspores in C303A.
Validation of key DEGs by Real-time Quantitative PCR (RT-qPCR)
To validate the sequencing data and possible pathways, the transcriptional expression of fifteen key DEGs involved in energy metabolism and pistillody were further analyzed by qRT-PCR (Fig. 10). Although quantitative difference existed in the two methods, the tendencies were same. The difference in the gene expression level was reasonable because of different working principles. Overall, these results demonstrated that our sequencing results were accurate and reliable, and further confirmed possibility on pathways related to stamen and pollen development.
Analysis of DEGs in energy metabolism pathway and determination ATP content and related enzyme activity
The DEGs were annotated and analyzed in KEGG pathway database to obtain the DEGs related to energy metabolism. Out of which 82 DEGs were annotated as NADH dehydrogenase, ATP synthase, Citrate synthase, Isocitrate dehydrogenase, 2-oxoglutarate dehydrogenase and Malate dehydrogenase in the energy synthesis pathway. Then heat map was used to analyze the expression level of DEGs in binuclear stage (Fig. 11). The results showed that the expression level of energy synthesis pathway genes was significantly down-regulated at the binuclear stage. To further confirm the accuracy of results described above, we measured the amount of ATP at the binucleate stage (Fig. 12). The energy metabolism-related enzymes genes were significantly lower in C303A than that of 303B presented a heat map in Figure 8. Because of the altered of energy metabolism-related enzymes genes in the C303A pistilldoy stamen, we considered that the ATP in the pistilldoy stamen was also lower when compared with 303B stamen. According to the assays, we found that the ATP content in C303A stamen was also significantly lower than that in 303B. So, from these results we hypothesized that these genes involved in energy metabolism are associated with C303A pistilldoy.
ROS assay and activities of antioxidant enzymes
Compared with the physiological research of fertile wheat, the previous research results showed that sterile wheat higher accumulation of H2O2 and MDA and had a higher O2− generation rate [16, 17]. Thus, we determined the O2− generation rate as well as the H2O2 and MDA contents during all of the anther developmental stages (Fig. 13). The rate of ROS production was significantly higher in C303A than that of its maintainer line during all of the anther developmental stages. Also, the O2−, H2O2 and MAD contents were continuously raised in C303A with peak values at the uninucleate stage, whereas the MDA contents peak in the trinucleate stage. Therefore, these results imply that excessive accumulation of ROS may lead to abnormal of tapetal cell in C303A stamens. Furthermore, we also measured the activities of antioxidant enzymes such as CAT, SOD and POD. The SOD and POD activity in C303A remained at a high level throughout the whole course of pollen development, while the CAT activity remained high in the pistillody stamen only from the tetrad stage to the later uninucleate stage. According to these outcomes, it is possible that upregulation of the activities of antioxidant enzyme also reflects the extreme accumulation of ROS in the C303A, which upsets the balance of the antioxidant defense system and finally could contribute to pistillody in C303A.