Comparative transcriptome analysis revealed the conversions of stamens into pistil-like structures in Aegilops crassa cytoplasmic male sterile wheat (Triticum aestivum)

Background: Aegilops crassa cytoplasm is an essential material for investigating the cytoplasm of cytoplasmic male sterility (CMS). Moreover, the stamens of C303A exhibit a high degree of pistillody, turning almost white. However, the underlying molecular mechanism of C303A pistillody remains unclear. Therefore, to gain a better understanding of C303A, the phenotypic and cytological features of C303A were observed to identify the key stage that the homeotic transformation of stamens into pistil-like structures, transcriptome profiles were determined by Illumina RNA sequencing technology (RNA-Seq) of stamen. Results: Through morphological observation, for CMS wheat with Aegilops crassa cytoplasm (CMS-C) line C303A, the pistil of which developed normally, but stamens were ultimately aborted and the stamens released no pollen when mature. According to the results of the paraffin section, stamens began to transform into pistils or pistil-like structures at binucleate stage (BNS). Therefore, the stamens of line C303A and its maintainer 303B at BNS were collected for transcriptome sequencing. A total of 20,444 wheat genes were detected as being differentially expressed between C303A and 303B stamens, included 10,283 up-regulated and 10,161 down-regulated genes. Gene Ontology Enrichment Analyses showed that most differentially expressed genes (DEGs) distributed on the metabolic process, cell, cellular process, catalytic activity and cell part. From KEGG, we knew that DEGs were mainly enriched to energy metabolism. We also found several essential genes that may contribute to pistillody in C303A. Based on the above analysis, we believe that due to the confusion of energy metabolism and reactive oxygen metabolism, thereby inducing the pistillody and eventually lead to the abortion of C303A. Conclusion:This study unravels the complex transcriptome profiles in C303A stamen, highlighting the energy metabolism and class B MADS-box genes related to pistillody. This

work should lay the foundations of future studies in the mechanical response to wheat stamen and pollen development in CMS.

Background
As a staple food for 35% of the world's population [1], wheat was the second largest staple food crop after rice and cultivated around 220 million hectares worldwide [2]. China was the largest producer and consumer of wheat around the world, with a cultivation area of about 24 million ha, an average yield of 4762 kg ha − 1 .Many significant challenges are facing in China, such as the increasing population and the reducing arable land area. Therefore, improving grain yield was an inevitable requirement to ensure food security. Therefore, improving grain yield was an inevitable requirement to ensure food security [3]. Overall, increasing wheat yield was a long-term goal of wheat breeding. And, utilization heterosis was a best method to increase yield and meet global food safety needs, such as maize, rape, sunflower, rice and sorghum [4]. Moreover, male sterility plants provide crucial breeding tools to harness hybrid vigor, or heterosis, in hybrid crops and also provide valuable material to study stamen and pollen development and nuclearcytoplasmic interactions [5].
Aegilops crassa cytoplasm is the vital source of the cytoplasm of cytoplasmic male sterility (CMS), which had no harmful effect on the agronomic characters of common wheat. C303A that conferred by the cytoplasm from Ae. Crassa, was an outstanding wheat germplasm resources in cytoplasmic male sterility. But there was a fatal flaw on C303A, just like complex restoration of fertility, poor outcrossing amount of wheat and fewer restorer lines. Due to these reasons, C303A was hard to use its heterosis but had great research potential as an extensive germplasm resource. For example, the stamens of C303A exhibit a high degree of pistillody, turning almost white because of the nuclear-cytoplasmic interaction when compared with its maintainer line 303B. In novel studies discovered an alloplasmic line of Norin 26 (N26) with Ae. crassa cytoplasm, which showed male sterility under long-day conditions (>15 h light period) because of pistillody [6]. So, we might assume that some factors in the Ae. crassa cytoplasm provoked pistillody, probably mitochondrial gene (s) [7].However, studies on the pistillody identity in wheat are still superficial. C303A is a significant material for studying the pistillody of wheat.
Therefore, this study describes some characteristics of C303A, including its floret morphology, cytological mechanism, physiological index, and molecular mechanism. At present, the molecular mechanism of stamen development and nuclear-cytoplasmic interactions were clearly explained in Arabidopsis thaliana, rice, maize and other model plants, but it is rarely studied in non-model plants. Wheat is a heterologous hexaploid plant with a complex genetic background, and it is inefficient and difficult to study a great quantity of genes involved in the stamen development and nuclear-cytoplasmic interactions network by traditional molecular biology methods. Therefore, high-throughput transcriptome sequencing (RNA-Seq) provides comprehensive and rapid access to almost all transcript information in wheat stamen during binucleate stage. So, it is easy to systematically analyze the DEGs which related to pistillody. Our research carried out RNA-Seq on the stamen of C303A and its maintainer line 303B and presented herein focused on the discovery of metabolic processes and transcription factors contributing to pistillody in C303A. A comprehensive understanding of the genetic network changes and expression patterns of C303A pistillody will enable the future utilization of molecular mechanism in stamen development and nuclear-cytoplasmic interactions.

Results
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 [1415].
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)  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 H 2 O 2 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 metabolismrelated 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. and MDA contents during all of the anther developmental stages (Fig. 13). Meanwhile, in order to show the content of O 2− or H 2 O 2 directly, we stained the the normal anther and pistillody anther by NBT and DAB (Fig. S3, S4). 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 O 2− , H 2 O 2 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.

Discussion
Wheat is an allohexaploid, which contained three genomes (A, B, and D), and its genome approximately 17G. The formation of common wheat involves three primitive ancestral species and two natural hybrids. Triticum urartu, Aegilops speltoides, Aegilops tauschii were the progenitor species genome of the wheat. Compared with rice, corn and other crops, wheat molecular basic research is still weak, and the genetic background is relatively poor. In recent years, RNA-Seq, which based on Illumina sequencing platform, has shown itself to be a reliable tool with an excellently diverse range of applications, from detailed researches of biological processes at the cell type-specific level, to give penetrations into fundamental questions in plant biology on an evolutionary time scale.
In this study, we performed comprehensive RNA-seq analyses of the stamen from C303A and 303B at the binucleate stage. Furthermore, compared with 303B, we detected many DEGs that involved the process of carbohydrate and energy metabolism of stamen in C303A and numerous of these DGEs decreased at the binucleate stage. Overall, we could suggest that the decrease or non-expression of these DEGs due to pistillody. And these DEGS may result in insufficient cellular energy supply and abnormal starch synthesis, which also disturbs the balance of material and energy metabolism in stamens and eventually leads to no pollen grain. showed that they were functionally interdependent and lacked a single gene. In our research found that AGPase was decreased at the binucleate stage in C303A. Moreover, exopolygalacturonase and Alpha-glucan phosphorylase were involved in starch synthesis and accumulation, which activated the corresponding metabolic pathways to improve the accumulation and synthesis of starch. Thus, these results showed that starch and sucrose metabolism in C303A anther was probably weakened that in 303B anther, which was closely associated with no pollen in C303A.

DEGs involved in Energy Metabolism
The normal energy metabolism process of plants can satisfy their growth and development. Energy metabolism mainly includes these major processes: such as photosynthesis, glycolysis, oxidative phosphorylation, and the TCA cycle. Energy metabolism process requires a lot of related genes or proteins to be adjusted together. A large number of investigations have revealed that the abnormal performance of these proteins or genes in stamen will inevitably interfere with the energy supply of pollen development-for example, using the fluorescein-luciferase way to determine the content ATP of the anther, which found that the content of ATP in AS was much lower than that in normal anther from the early stage of pollen grains to maturation phase [2324]. Similarly, Xia and Liu [2425] used the same method to determine the ATP content of stamen at different developmental stages using the corn CMS line and its maintainer system. The results showed that the content of ATP of the sterile line was significantly lower than that of the maintainer in the anther, and it was considered that a large amount of energy was consumed during the pollen abortion process because the formation of microspores requires a large amount of energy and nutrients. In addition, the absence of microspores in C303A may be due to the down-regulation of energy metabolism in C303A during pistillody stamens formation. Moreover, Hexokinase-5, 6-phosphofructokinase 5 and pyruvate kinase played a vital role in the glycolysis pathway, which was significantly down-regulated in C303A. These genes may have an impact on the development of stamen. The down-regulation of these significant genes directly led to the number of respiratory substrates decreased and can also interfere with the electron transmit chain (ETC) on mitochondria. Also, many DEGs that related to ETC and TCA decreased in this stage and our examinations of the ATP contents also uphold this standpoint ( Figure Fig.  12). So, we could suggest that the energy of C303A pistillody stamens was lower than that of 303B stamens, which can not meet the basic energy requirement of microspores in C303A. As a result, the pistillody stamens of C303A cannot form microspores.

DEGs Related to Pistillody
When the stamen of flower showed the homeotic transformation of stamens into pistil-like structures, the phenomenon is called pistillody, which can cause stable and complete male sterility, and existed in some species, such as in Arabidopsis [2526] and Antirrhinum [1314,2627]. In present study, pistillody advents were occurred in CMS wheat line C303A with Aegilops crassa cytoplasm, but the pistillody phenomenon of C303A was not clear, so we used RNA-seq to elucidate the molecular mechanism of pistillody induction. Therefore, we selected some reported genes that affected flower development. wheat. In the present study, we find a paralogous gene (TaAGL14) may have an analogous role to OsMADS32. In the research of cfo1-1, the strong alleles of OsMADS32, discovered that the number of pistil was increased and OsMADS32 was required for pistil development and floral meristem determinacy. Furthermore, in C303A gene TaAGL14 was upregulated.
Thus, class-B MAD-box genes may contribute to pistil development.
developmental processes, such as ovule development, vegetative growth, flower morphogenesis and fruit formation [31, 32]. The MADS-box transcription factor, WM27A, the expression pattern was compatible with its function as the class D gene and regulated ovule identity specification on the basis of the "ABCDE" model in flower development [33].
Previous studies found that WM27A was deeply expressed in pistils and caryopses, with weak expression in stamens and extremely high expression during late spike development. Moreover, the phylogenetic tree of the Nucleotide sequences of a MADS-box transcription factor TaMS-MADSbox (GenBank accession number: 36925702) which linked to fertility conversion in male sterile wheat lines [34], along with MADS-box transcription factor from our sequencing data, showed that WM27A and TaMS-MADSbox in the same branch. In our study, WM27A expression was unregulated in C303A. Thereby, it is obvious to find that this result is very similar to the results of previous studies. EPIDERMAL PATTERNING FACTOR-LIKE (EPFL) family functions in various plant growth and development aspects, such as the guidance of inflorescence architecture and pedicel length. And in wheat TaEPFL1 secreted peptide gene is essential for stamen development [35]. Moreover, TaEPFL1 gene is expressed at an abnormally high level in pistillody stamens compared with that in pistils and stamens. In our study, we found that genes were coding the EPFL protein in C303A. On the whole, these results could speculate that these genes brought about pistillody. However, the mechanisms by which these genes control wheat stamen and pistil development remains to be widely explored.

Conclusion
For CMS-C line C303A, the pistil of which developed normally, but stamens were ultimately aborted and the stamens released no pollen when mature. According to the results of the paraffin section, stamens began to transform into pistils or pistil-like structures at binucleate stage (BNS). And we employed RNA-Seq results showed that mainly via effects on carbohydrate and energy metabolism, and downregulating the expression of enzyme genes that involved in carbohydrate and energy metabolism, which were closely related to pollen development. Thus, we suggested a possible regulatory network and did deem that the downregulation of key genes in this regulatory network, which is the leading cause of Because C303A had no microspore, we determined the developmental stage of C303A by observing the characteristics of 303B. And previous studies have shown that the spike development stages of 303B and C303A are the same [36,37]. Therefore, the average of the spike length and the external morphological characteristics of each period of 303B were determined. In the tetrad stage, the spike length of the 303B was 3.7 cm and the spike apical was in the middle of the second top leaves and the third top leaves. At the early uninucleate stage, the spike length of the 303B was 10.5 cm and the spike apical was in the middle of the second top leaves and the flag leaf. At the later uninucleate stage, the spike length of the 303B was 12.9 cm, and the awn was just wholly extracted.
At the binucleate stage, the spike length of the 303B was 11.4 cm, and wheat earing was about 2/3. At the trinucleate stage, the spike length of the 303B was14.5 cm and all spike have been pulled out, which also found some of glume opening and anther exposed. Total RNA extraction, cDNA library preparation and

Illumina sequencing
Total RNA was extracted into C303A and 303B stamen during BNS stage using Trizol reagent (Takara Biotechnology, Dalian, China). The integrity of total RNA was identified by 1% agarose gel electrophoresis, the sample OD260 and OD230 were detected by Nanodrop spectrophotometer, then the purity of RNA was analyzed, and the sample concentration, RIN value and 28S/18S value were determined by Agitent2100 bioanalyzer Agilent Technologies, Santa Clara, CA, USA).
All samples were given to Sagene Biotech Co. Ltd (Guangzhou, China) for library and sequencing, and the eukaryotic mRNA was enriched with Oligo (dT) beads. And then the fragmentation buffer was used to interrupt the mRNA into short fragments to mRNA. cDNA is synthesized by using mRNA as a template by using random primers. DNA polymerase I, RNase H, dNTP and buffer were added to synthesize Second-strand Then those were purified with QiaQuick PCR extraction kit (GeneStar Co. Ltd, Beijing, China). The purified cDNA was needed to end repaired, poly (A) added. Finally, PCR amplification was performed and the PCR product was purified with AMPure XP beads Beckman Coulter, Brea, CA, USA to obtain a final library. After the library is qualified, different libraries are pooled according to the requirements of the effective concentration and the target data volume, and then the Illumina HiSeq sequence is performed.

Raw data filtering and transcripts splicing
The process of raw data filtering was as follows: (1) removing reads with adapter; (2) eliminating reads containing more than 10% of unknown nucleotides (N); (3) excluding low-quality reads. In order to ensure the quality of information analysis, raw reads were filtered to obtain clean reads, and subsequent analysis is based on clean reads.

Physiological indexes
The content of ATP was determined by spectrophotometric method according to the protocol of ATP kit (Comin Biotechnology Co., Ltd., Suzhou, China). Three ground C303A anther samples (0.1 g each) were placed in a centrifuge tube. Added 1 mL lysate and stirred well. The mixture was well homogenized for 10 min, centrifuged at 8000 r/min for 10 min at 4 °C, and the supernatant was taken to another centrifuge tube. Added 500 μL of chloroform and mixed well by shaking. Centrifuge at 10000 r/min for 3 min at 4 °C, take the supernatant, and place on ice for the ATP testing. And then though spectrophotometric method determined the content of ATP. In the same way, we measured the ATP content in the anther of 303B.