Comparative Metabolomic Analysis of Multi-Ovary Wheat under Heterogeneous Cytoplasm Suppression

: The multi-ovary trait of wheat inbred line DUOII is controlled by a dominant gene whose expression can be suppressed by the heterogeneous cytoplasm of TeZhiI (TZI), another inbred line with the nucleus of common wheat and the cytoplasm of Aegilops . DUOII ( ♀ ) × TZI ( ♂ ) shows multi-ovary trait, while TZI ( ♀ ) × DUOII ( ♂ ) shows mono-ovary. To elucidate the molecular mechanism regulating heterogeneous cytoplasmic suppression of the multi-ovary gene, we performed an untargeted metabolomic analysis of 2–6 mm young spikes of reciprocal crosses between DUOII and TZI at the critical stage of additional pistil primordium development. We identiﬁed 198 annotated differentially expressed metabolites and analyzed them according to their biological functions. The results showed that these metabolites had obvious functional pathways mainly implicated in amino acid, carbohydrate, nicotinate and nicotinamide, and purine metabolism and isoquinoline alkaloid biosynthesis. We also found that shikimate, phosphoglycolic acid, nicotinamide, guanine, and xanthine might play essential roles in cytoplasmic suppression of multi-ovary trait. Chloroplast metabolism was also implicated in the nuclear-cytoplasmic effect of the multi-ovary gene. The ﬁndings provide solid theoretical and empirical foundations for future studies elucidating the mechanisms controlling heterogeneous cytoplasmic suppression of the nuclear multi-ovary gene in


Introduction
Wheat (Triticum aestivum L.) is grown worldwide and is vital to global food security. It comprises~26% of all cereal production and provides~20% of the protein, 18% of the calories, and 3% of the fat consumed by humans [1]. Wheat production is expected to increase in response to growing food demand associated with expanding population growth [2]. Given that the cultivated land is limited, research continues to focus on maximizing wheat yield per unit area. Generally, each wheat floret forms only one grain. However, multi-ovary wheat is a new variety that has 1-3 ovaries and three stamens and can set 1-3 grains per floret [1,3]. This trait is genetically stable and can produce a relatively greater number of grains per spike. For example, the grain number of the main stem spike of DUOII was 130.4 [1]. Therefore, multi-ovary wheat could potentially increase wheat yield and facilitate the study of floral development mechanisms.
Since the wheat multi-ovary trait was first reported in 1983 [3], studies have mainly focused on floral organ development [4], biochemical basis of seed germination [5], discovery of molecular markers [6], and multi-ovary development mechanisms [7]. Previously, we found that the multi-ovary trait of DUOII was controlled by a dominant gene whose expression could be suppressed by the heterogeneous cytoplasm of TeZhiI (TZI), an alloplasmic line with the nucleus of common wheat and the cytoplasm of Aegilops [8]. The genetic basis of the multi-ovary trait of DUOII has been clarified. However, little is known about the mechanism by which TZI suppresses the expression of the multi-ovary gene in common wheat, although we have studied this process through genomic DNA methylation [9], transcriptomic [10] and proteomic [11] analyses.
Metabolomics, one of the most recent of the "-omics" technologies, is used to characterize the final products of the genetic, gene expression, and protein activities in a given moment and condition [12]. Compared to genomics, transcriptomics, and proteomics, the metabolomics studies are easier to relate to phenotypes. This is because the functions of genomic region and expressed sequences and proteins can undergo epigenetic regulation, post-translational modifications, and activity impairments depending on substrate and cofactor availability [13,14]. Metabolomics support high throughput and rapid detection of changes in a wide range of metabolites and provide in-depth analyses of the total metabolome of biological processes in plants [15][16][17]. In wheat, metabolomics has been successfully applied for seed germination [18], tiller [19], anther [20], pistil [21], and grain development [22], grain filling [23], and responses to various biotic and abiotic stresses [24][25][26][27][28].
In the present study, we conducted an untargeted metabolomic analysis to compare the alterations in the metabolite profiles of F 1 s derived from the reciprocal crosses between DUOII and TZI, and a considerable number of metabolites were identified to show differences in their relative abundances in these two crosses. The metabolomic results provided new insight into the underlying molecular mechanisms that control the heterogeneous cytoplasmic suppression of the multi-ovary gene in wheat. The findings will provide theoretical foundation for understanding the development of multi-ovary in wheat, and help to identify key metabolites or genes involved in this process. Further, these key metabolites or genes can be developed to biomarkers of multi-ovary trait, which can be implemented in future breeding programs focusing on the development of high yield wheat cultivars using multi-ovary wheat.

Plant Materials
The two inbred wheat lines DUOII (a common multi-ovary line) and TZI (an alloplasmic line with the nucleus of the "Chris" wheat variety and an Aegilops cytoplasm) were used here. Seeds of both lines were maintained at the Key Laboratory of Crop Heterosis of Shaanxi Province, Shaanxi, China. Before the onset of this experiment, the lines had been selfed for more than twenty years and there was no segregation. Hence, their genetic backgrounds were pure. In October 2016, DUOII and TZI were sown in the experimental field of Northwest A & F University, Yangling, Shaanxi, China (34 • 91 N, 106 • 86 E). In May 2017, reciprocal crosses between DUOII and TZI were performed and the F 1 seeds were sown in October 2017. In March 2018, 2-6-mm young spikes were hand-dissected from~100 plants in each F 1 population (DUOII (♀) × TZI (♂) and TZI (♀) × DUOII (♂)), immediately frozen in liquid nitrogen, and stored at −80 • C until metabolite extraction.

Morphological Analysis
The multi-ovary trait was examined during the flowering season and maturation stage, and the numbers of ovaries and grains in the florets of each spike were investigated. Spikes and grains were photographed with a Nikon D600 digital camera (Nikon, Tokyo, Japan). Pistils were photographed with a Nikon E995 digital camera (Nikon, Tokyo, Japan) mounted on a Motic K400 dissecting microscope (Preiser Scientific, Louisville, KY, USA).

Young Spike Metabolite Extraction
Young spikes (50 mg) were placed in a 2 mL microcentrifuge tube (Eppendorf, Hamburg, Germany) containing 1 mL extraction solvent (methanol:acetonitrile:water (v/v/v) = 2:2:1) and 2 µg/mL 2-chloro-L-phenylalanine (internal standard, Hengbai Biotech Co., Ltd., Shanghai, China). The mixture was vortexed for 30 s, homogenized in a ball mill for 4 min at 45 Hz, ultrasound-treated for 5 min (incubated in ice water). The homogenization was repeated thrice and the samples were then incubated at 20 • C for 1 h and centrifuged at 12,000 rpm for 15 min at 4 • C. The resulting supernatant (500 µL) was transferred to a fresh microcentrifuge tube and dried in a vacuum concentrator without heating. The dried pellet was redissolved in 100 µL extraction solvent (acetonitrile:water (v/v) = 1:1), vortexed for 30 s, sonicated for 10 min (ice water bath), and then centrifuged at 12,000 rpm for 15 min at 4 • C. The supernatant (60 µL) was carefully transferred to a fresh 2 mL LC/MS glass vial for untargeted metabolomic analysis. Ten-microliter aliquots were taken per sample. The aliquots were pooled and the combination was used as a quality control (QC) sample. Six technical replicates were performed in the metabolomic analysis.

Metabolomic Analysis
Reconstituted samples were analyzed in an ultra-high performance liquid chromatography coupled to quadrupole time of fight mass spectrometry with an electrospray ionization source (UHPLC-ESI/QTOF-MS) system using an untargeted metabolomics method. A 1290 UHPLC system (Agilent Technologies, Santa Clara, CA, USA) and a Triple TOF 6600 mass spectrometer (AB Sciex, Framingham, MA, USA) were used.
The UHLPC-ESI/QTOF-MS untargeted metabolomics analysis was conducted as described by [29], with certain modifications. Briefly, UHPLC chromatographic separation was performed with a binary mobile phase consisting of 25 mM NH 4 Ac and 25 mM NH 4 OH in water (pH = 9.75) (A) and acetonitrile (B) and an ACQUITY UPLC BEH Amide column (100 mm × 2.0 mm; 1.7 µm; Waters Corp., Milford, MA, USA). The linear gradient was as follows: 0 min, 95% B; 0.5 min, 95% B; 7 min, 65% B; 8 min, 40% B; 9 min, 40% B; 9.1 min, 95% B; and 12 min, 95% B. The flow rate of the mobile phase was 0.5 mL/min. All samples were analyzed in positive and negative mode and the injection volume was 1.5 µL. A triple time-of-flight mass spectrometer was used because it could acquire MS/MS spectra on an information-dependent basis (IDA) during LC/MS. In this mode, the acquisition software Analyst TF v. 1.7 (AB Sciex, Framingham, MA, USA) continuously evaluated the full-scan survey MS data as they were collected and triggered the acquisition of MS/MS spectra depending on the preselected criteria. In each cycle, 12 precursor ions with intensities > 100 were selected for fragmentation at collision energy (CE) of 30 eV (15 MS/MS events with product ion accumulation time of 50 msec each). The ESI source conditions were as follows: ion source gas 1, 60 psi; ion source gas 2, 60 psi; curtain gas, 35 psi; source temperature, 650 • C; ion spray voltage floating (ISVF), 5000 V or −4000 V in positive or negative mode, respectively.

Data Preprocessing and Annotation
Raw data files in wiff format were converted to mzXML files with ProteoWizard [30] and processed by R package XCMS (version 3.2) [31]. The preprocessing results generated a data matrix comprising retention time, mass-to-charge ratio (m/z), and peak intensity. The "CAMERA" package in R [32] was used for peak annotation after XCMS data processing. The in-house MS2 database and KEGG [33], PubChem [34], ChEBI [35], and HMDB [36] were used for metabolite identification. Functional annotations of the identified metabolites were extracted from HMDB, PubChem, ChEBI, and KEGG.

Statistical Analysis
Total ion chromatograms (TICs), Pearson's correlations, and principal component analysis (PCA) for different samples were used to assess the metabolomics data. Supervised orthogonal projections to latent structures-discriminate analysis (OPLS-DA) was applied to clarify the global metabolic changes between TZI × DUOII and DUOII × TZI according to the method of [37]. The evaluation parameters for the OPLS-DA model were R 2 X, R 2 Y, and Q 2 Y. The model was considered reliable when R 2 X, R 2 Y, and Q 2 Y neared 1, adequate when Q 2 Y > 0.5, and excellent when Q 2 Y > 0.9. The robustness and predictive ability of our model were evaluated by 7-fold cross validation. The corresponding variable importance in the projection (VIP) were assigned to the metabolites and expressed their contribution to OPLS-DA model construction. Only metabolites with VIP > 1 and Student's t-test p < 0.05, were considered as differentially expressed metabolites (DEMs). Commercial databases including KEGG and MetaboAnalyst [38] were utilized to search for the metabolic pathways of DEMs.

Morphological Characterization of F 1 Plants
As a rule, there is only one ovary per wheat floret. However, DUOII is a multi-ovary wheat variety with 1-3 pistils and three stamens per floret. It can set 1-3 grains per floret. The reciprocal crosses (F 1 s) between multi-ovary wheat DUOII and alloplasmic wheat TZI had similar spike phenotypes. Nevertheless, they differed in terms of their multi-ovary trait ( Figure 1). All 50 DUOII (♀) × TZI (♂) plants expressed the multi-ovary trait, whereas all 50 TZI (♀) × DUOII (♂) plants expressed the mono-ovary trait. Further, there were no other significant differences between the crosses. Thus, the heterogeneous cytoplasm only influenced the expression of multi-ovary-related traits but affected no other morphological traits.

Statistical Analysis
Total ion chromatograms (TICs), Pearson's correlations, and principal component analysis (PCA) for different samples were used to assess the metabolomics data. Supervised orthogonal projections to latent structures-discriminate analysis (OPLS-DA) was applied to clarify the global metabolic changes between TZI × DUOII and DUOII × TZI according to the method of [37]. The evaluation parameters for the OPLS-DA model were R 2 X, R 2 Y, and Q 2 Y. The model was considered reliable when R 2 X, R 2 Y, and Q 2 Y neared 1, adequate when Q 2 Y > 0.5, and excellent when Q 2 Y > 0.9. The robustness and predictive ability of our model were evaluated by 7-fold cross validation. The corresponding variable importance in the projection (VIP) were assigned to the metabolites and expressed their contribution to OPLS-DA model construction. Only metabolites with VIP > 1 and Student's t-test p < 0.05, were considered as differentially expressed metabolites (DEMs). Commercial databases including KEGG and MetaboAnalyst [38] were utilized to search for the metabolic pathways of DEMs.

Morphological Characterization of F1 Plants
As a rule, there is only one ovary per wheat floret. However, DUOII is a multi-ovary wheat variety with 1-3 pistils and three stamens per floret. It can set 1-3 grains per floret. The reciprocal crosses (F1s) between multi-ovary wheat DUOII and alloplasmic wheat TZI had similar spike phenotypes. Nevertheless, they differed in terms of their multi-ovary trait ( Figure 1). All 50 DUOII (♀ ) × TZI (♂ ) plants expressed the multi-ovary trait, whereas all 50 TZI (♀ ) × DUOII (♂ ) plants expressed the mono-ovary trait. Further, there were no other significant differences between the crosses. Thus, the heterogeneous cytoplasm only influenced the expression of multi-ovary-related traits but affected no other morphological traits.

Assessment of Metabolomic Data
Metabolic profiling of 2-6 mm young spikes from TZI × DUOII and DUOII × TZI plants was performed using an untargeted UHPLC-ESI/QTOF-MS analysis. After relative standard data pre-processing, we identified 3325 and 2783 metabolites in the positive and negative ion mode, respectively. For positive or negative ion mode, TICs of QC samples showed no significant differences, and the retention times and peak intensities were reproducible and stable (Supplementary Figure S1). Hence, the used UHPLC-ESI/QTOF-MS system was reliable and effective. Additionally, Pearson's correlation among the three QC samples were higher than 0.992 (Figure 2a,c). Further, principal component analysis was performed to detect metabolite distributions for the metabolomic data. The QC samples

Assessment of Metabolomic Data
Metabolic profiling of 2-6 mm young spikes from TZI × DUOII and DUOII × TZI plants was performed using an untargeted UHPLC-ESI/QTOF-MS analysis. After relative standard data pre-processing, we identified 3325 and 2783 metabolites in the positive and negative ion mode, respectively. For positive or negative ion mode, TICs of QC samples showed no significant differences, and the retention times and peak intensities were reproducible and stable (Supplementary Figure S1). Hence, the used UHPLC-ESI/QTOF-MS system was reliable and effective. Additionally, Pearson's correlation among the three QC samples were higher than 0.992 (Figure 2a,c). Further, principal component analysis was performed to detect metabolite distributions for the metabolomic data. The QC samples were tightly clustered together, and the samples were located in Hotelling's T-squared ellipses (95% confidence interval) (Figure 2b,d). Therefore, the metabolomic data had a sufficiently high quality to support subsequent analyses.
were tightly clustered together, and the samples were located in Hotelling's T-squared ellipses (95% confidence interval) (Figure 2b,d). Therefore, the metabolomic data had a sufficiently high quality to support subsequent analyses.

Identification of DEMs
To identify the metabolites describing changes between TZI × DUOII and DUOII × TZI plants, the supervised OPLS-DA were performed. The OPLS-DA score scatter plots indicated that the modeling highlighted a clear separation between TZI × DUOII and DUOII × TZI in positive or negative ion mode (Figure 3a,c). The OPLS-DA model fitting parameters were effective; the goodness-of-fit R 2 Y > 0.9 and the goodness-of-prediction Q 2 Y > 0.5. In positive ion mode, R 2 Y = 0.962 and Q 2 Y = 0.589 while in negative ion mode, R 2 Y = 0.969 and Q 2 Y = 0.634. Moreover, the permutation test confirmed the validation of OPLS-DA models (Figure 3b,d).

Identification of DEMs
To identify the metabolites describing changes between TZI × DUOII and DUOII × TZI plants, the supervised OPLS-DA were performed. The OPLS-DA score scatter plots indicated that the modeling highlighted a clear separation between TZI × DUOII and DUOII × TZI in positive or negative ion mode (Figure 3a,c). The OPLS-DA model fitting parameters were effective; the goodness-of-fit R 2 Y > 0.9 and the goodness-of-prediction Q 2 Y > 0.5. In positive ion mode, R 2 Y = 0.962 and Q 2 Y = 0.589 while in negative ion mode, R 2 Y = 0.969 and Q 2 Y = 0.634. Moreover, the permutation test confirmed the validation of OPLS-DA models (Figure 3b,d). were tightly clustered together, and the samples were located in Hotelling's T-squared ellipses (95% confidence interval) (Figure 2b,d). Therefore, the metabolomic data had a sufficiently high quality to support subsequent analyses.

Metabolic Pathway Analysis of DEMs
To elucidate the metabolic pathways that are involved in the suppression of multiovary gene in TZI × DUOII, a detailed metabolic pathway analysis was performed by using Oryza sativa japonica as the pathway libraries to associate the biological functions of identified known DEMs to different pathways. Sixty-six of the 198 known DEMs were placed in 50 KEGG pathways introduced to MetaboAnalyst (Supplementary Table S4).
We identified 11 significant pathways by combining the p and pathway impact values ( Figure 5). They comprised five upregulated and 30 downregulated DEMs. These pathways were divided into three subsets by KEGG classification (Table 1). Subset Ⅰ included arginine biosynthesis (five DEMs), glutathione metabolism (six DEMs), phenylalanine, tyrosine, and tryptophan biosynthesis (five DEMs), alanine, aspartate, and glutamate metabolism (four DEMs), tryptophan metabolism (four DEMs), and arginine and proline metabolism (four DEMs). Hence, all these six pathways were related to amino acid metabolism and involved 21 different DEMs. Of these, 18 were downregulated in TZI × DUOII. Therefore, amino acid metabolism might have been weakened by the heterogeneous cytoplasm. Subset Ⅱ consisted of butanoate metabolism (four DEMs) and glyoxylate and dicarboxylate metabolism (five DEMs). These pathways were classified as carbohydrate metabolism and all DEMs were downregulated. Thus, these metabolic processes were downregulated in TZI × DUOII. Subset III involved isoquinoline alkaloid biosynthesis (two DEMs), nicotinate and nicotinamide metabolism (three DEMs), and purine metabolism (nine DEMs). They were separately classified into biosynthesis of other secondary metabolism, metabolism of cofactors and vitamins, and nucleotide metabolism, respectively. Both two DEMs in isoquinoline alkaloid biosynthesis were upregulated while all three DEMs in metabolism of nicotinate and nicotinamide were downregulated. In

Metabolic Pathway Analysis of DEMs
To elucidate the metabolic pathways that are involved in the suppression of multiovary gene in TZI × DUOII, a detailed metabolic pathway analysis was performed by using Oryza sativa japonica as the pathway libraries to associate the biological functions of identified known DEMs to different pathways. Sixty-six of the 198 known DEMs were placed in 50 KEGG pathways introduced to MetaboAnalyst (Supplementary Table S4).
We identified 11 significant pathways by combining the p and pathway impact values ( Figure 5). They comprised five upregulated and 30 downregulated DEMs. These pathways were divided into three subsets by KEGG classification (Table 1). Subset I included arginine biosynthesis (five DEMs), glutathione metabolism (six DEMs), phenylalanine, tyrosine, and tryptophan biosynthesis (five DEMs), alanine, aspartate, and glutamate metabolism (four DEMs), tryptophan metabolism (four DEMs), and arginine and proline metabolism (four DEMs). Hence, all these six pathways were related to amino acid metabolism and involved 21 different DEMs. Of these, 18 were downregulated in TZI × DUOII. Therefore, amino acid metabolism might have been weakened by the heterogeneous cytoplasm. Subset II consisted of butanoate metabolism (four DEMs) and glyoxylate and dicarboxylate metabolism (five DEMs). These pathways were classified as carbohydrate metabolism and all DEMs were downregulated. Thus, these metabolic processes were downregulated in TZI × DUOII. Subset III involved isoquinoline alkaloid biosynthesis (two DEMs), nicotinate and nicotinamide metabolism (three DEMs), and purine metabolism (nine DEMs). They were separately classified into biosynthesis of other secondary metabolism, metabolism of cofactors and vitamins, and nucleotide metabolism, respectively. Both two DEMs in isoquinoline alkaloid biosynthesis were upregulated while all three DEMs in metabolism of nicotinate and nicotinamide were downregulated. In contrast, there were eight downregulated and one upregulated DEMs in purine metabolism. Therefore, all three pathways were altered in TZI × DUOII. Moreover, there were relatively clear relationships among the 11 pathways (Supplementary Figure S3). They formed a pathway network and might play essential roles in heterogeneous cytoplasmic suppression of the multi-ovary trait.
contrast, there were eight downregulated and one upregulated DEMs in purine metabolism. Therefore, all three pathways were altered in TZI × DUOII. Moreover, there were relatively clear relationships among the 11 pathways (Supplementary Figure S3). They formed a pathway network and might play essential roles in heterogeneous cytoplasmic suppression of the multi-ovary trait. Figure 5. KEGG pathway analysis of DEMs. The circles represent different KEGG pathway. The Impact is the value calculated from pathway topology analysis by MetaboAnalyst, and the value is bigger, the pathway is more important. The size of circle represents the number of DEMs in this pathway, and the color represents the p value, which is the significance level in enrichment analysis statistics. Details of DEMs in each pathway are listed in supplementary Table S4.  The Impact is the value calculated from pathway topology analysis by MetaboAnalyst, and the value is bigger, the pathway is more important. The size of circle represents the number of DEMs in this pathway, and the color represents the p value, which is the significance level in enrichment analysis statistics. Details of DEMs in each pathway are listed in Supplementary Table S4.

Identification of Important DEMs
Stricter criteria were used with VIP > 1.5, p < 0.01, and Fold Change >1.50 or <0.66 to identify important DEMs within the 35 DEMs involved in the 11 significant metabolic pathways. As a result, we identified shikimate, phosphoglycolic acid, nicotinamide, guanine, and xanthine as the most important DEMs (Table 2). They were all downregulated, enriched in various metabolic pathways, and might play critical roles in suppressing multi-ovary in TZI × DUOII.

Discussion
In general, each flower has a constant ovary number. Under specific conditions, however, the number of ovaries may vary. Increases in the number of ovaries have been reported for Arabidopsis thaliana [39], rice [40], and wheat [7]. Multi-ovary wheat has 1-3 pistils and three stamens per floret and can set more than one grain per floret. These traits increase the grain number per spike in multi-ovary wheat. Hence, multi-ovary wheat is useful for studying the mechanisms underlying multi-ovary trait and floral development in wheat and improving potential yield and propagation coefficients in wheat breeding.
DUOII is a multi-ovary wheat variety. Its multi-ovary trait is controlled by a dominant gene that is suppressed by the heterogeneous cytoplasm of TZI [8]. This study provided a detailed comparable metabolites analysis of 2-6 mm young spikes obtained from TZI × DUOII and DUOII × TZI at the critical stage of additional pistil development. We identified 198 annotated DEMs and found that 82.83% of them were downregulated. Hence, suppression of the additional pistil was associated with a relatively weak metabolite network. We also identified 11 significant metabolic pathways and five important DEMs possibly playing vital roles in heterogeneous cytoplasm suppression of multi-ovary trait in TZI × DUOII.
Normal plant development depends on the complex and precise interaction of various metabolic processes, in which proteins are the primary executors of life functions and reflect complex molecular and physiological processes [41,42]. Amino acids are bioactive macromolecules that form the basic units of proteins. Previous studies showed that proteins and amino acids play crucial roles in plant development [20,43]. Flowers utilize amino acids as building blocks for enzyme and structural protein biosynthesis and precursors of N-containing secondary metabolites and signaling molecules [44]. The concentrations of several amino acid metabolites vary significantly during anther development in maize [45] and wheat [20]. Cysteine and methionine metabolism play pivotal roles in the development of wheat pistillody, another phenomenon of ovary number increasing [21]. Here, 49 DEMs were classified into amino acids and their derivatives. This category was the largest and accounted for 24.7% of all annotated DEMs. The DEMs of free amino acids included glycine, L-arginine, L-glutamate, L-tyrosine, and L-tryptophan. In addition, six of the 11 significant metabolic pathways were related to amino acid metabolism, including arginine biosynthesis, arginine and proline metabolism, alanine, aspartate, and glutamate metabolism, glutathione metabolism, phenylalanine, tyrosine, and tryptophan biosynthesis, and tryptophan metabolism. Therefore, the levels of these amino acids and the activity of their corresponding metabolic pathways were altered and they might have contributed to the suppression of additional pistils in TZI × DUOII wheat. This conclusion aligns with our previous study [11].
Carbohydrate molecules are critical for energy storage and organ formation in plants [20,46]. The changes in carbohydrates that occur during flower-inducing conditions influence flower bud differentiation and flower development [47][48][49]. Under carbohydratelimited conditions, Arabidopsis inhibited flower primordium initiation and aborted flowers and very young siliques [50]. Carbohydrate content and metabolism are associated with anther and pistil development [51,52]. Pectin is essential for pistil development [53,54]. At low pectin levels, tissue morphogenesis is abnormal in rice pistils [55]. In the present study, there were 11 downregulated and four upregulated DEMs in the carbohydrate classifica-tion, and pectin was downregulated in the TZI × DUOII plants. Among the significant KEGG pathways of DEMs, two terms of butanoate metabolism and glyoxylate and dicarboxylate metabolism were downregulated in TZI × DUOII. All these processes are related to carbohydrate metabolism according to KEGG pathway hierarchy. Thus, carbohydrate content and metabolism were disrupted in the TZI × DUOII plants. For this reason, these hybrids might have had insufficient energy and metabolites to support the development of additional pistil of multi-ovary trait.
Along with the eight aforementioned KEGG pathways, the three other significant metabolic pathways identified in this study were isoquinoline alkaloid biosynthesis, nicotinate and nicotinamide metabolism, and purine metabolism. They were separately classified into biosynthesis of other secondary metabolites, metabolism of cofactors and vitamins, and nucleotide metabolism, respectively. These three pathways participate in various biological processes in plants including flower development [56][57][58][59]. Five major isoquinoline alkaloids were detected in the carpeloid but not the normal stamens of Papaver somniferum L. [60,61]. Nicotinate and nicotinamide metabolism contribute to early ear inflorescence development in maize [57]. Purine metabolism is associated with anther and pollen development in rice [59] and wheat [58]. Consequently, changes in these pathways and their associated metabolites might suppress additional pistil formation in TZI × DUOII plants.
Of the 35 DEMs involved in these 11 significant KEGG pathways, five were detected by using more stringent criteria. Shikimate, phosphoglycolic acid, nicotinamide, guanine, and xanthine were all downregulated in TZI × DUOII. In plants, shikimate is a precursor of the aromatic amino acids phenylalanine, tyrosine, and tryptophan, indole and its derivatives, vitamins K 1 and B 9 , several alkaloids, tannins, and lignin. The shikimate pathway for aromatic amino acid biosynthesis is localized to the chloroplasts [62][63][64]. Shikimate pathway downregulation might result in insufficient precursors for the biosynthesis of some substances and alter normal chloroplast metabolism in TZI × DUOII. Phosphoglycolic acid (2-phosphoglycolate) is synthesized by the action of oxygenase on ribulose 1,5-bisphosphate in the chloroplast [65] and inhibits at least two key enzymes in chloroplast carbon metabolism including triosephosphate isomerase (Calvin cycle) and phosphofructokinase (glycolysis) [66,67]. Therefore, 2-phosphoglycolate downregulation could affect certain metabolic processes in TZI × DUOII chloroplasts.
Nicotinamide, the amide of nicotinic acid, is an active component of the coenzymes NAD(P) and NADP(H), participates in numerous enzymatic redox reactions, and is required by all organisms for energy metabolism [68]. Nicotinamide and nicotinic acid help protect plants against various stressors [69] and are involved in pyridine nucleotide biosynthesis [70]. They also exert epigenetic effects by decreasing DNA methylation [71]. Our previous study showed that DNA methylation is associated with multi-ovary trait suppression [9]. Hence, nicotinamide downregulation may affect biosynthesis and energy metabolism, increase DNA methylation, and indirectly suppress multi-ovary gene expression.
Guanine and xanthine are in the purine metabolic pathway. Guanine is a nucleobase in purine nucleotides which are building blocks of nucleic acids and certain cofactors [72]. Purine nucleotides are precursors of purine alkaloids and cytokinins in plants [73]. Plants may use purine nucleotides as nutrient sources in biodegradation reactions [74][75][76]. In purine catabolism, xanthine is the first common intermediate and free purine nucleobases converge on it. Subsequent degradation generates glyoxylate, carbon dioxide, and ammonia through a series of oxidation and hydrolysis reactions [72,77]. In this study, guanine and xanthine were downregulated. Thus, purine nucleotide biosynthesis and degradation were altered in TZI × DUOII and could have influenced the expression of multi-ovary trait.

Conclusions
Overall, the present research showed that the suppression of multi-ovary trait was associated with DEMs during the critical development stage of additional pistil primordia (2-6 mm young spikes). Thirty-five DEMs were significantly enriched in 11 KEGG path-Agronomy 2021, 11, 658 10 of 13 ways and five key DEMs were identified. We demonstrated that chloroplast metabolism plays a vital role in cytoplasmic suppression of multi-ovary trait. This discovery corroborated our previous transcriptomic and proteomic analyses [10,11]. In both studies, we found that chloroplast metabolism was altered in TZI × DUOII, which might transmit signal to nucleus and trigger some changes, and these changes then influenced the development of additional pistil primordia and resulted in the suppression of multi-ovary trait. The results of this study provide solid theoretical and empirical foundations for the elucidation of the mechanism controlling heterogeneous cytoplasmic suppression of the nuclear multi-ovary gene in wheat.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11040658/s1, Figure S1: Total ion chromatograms (TICs) of quality control (QC) samples; Figure S2: Volcano plots showing DEMs between TZI × DUOII and DUOII × TZI; Figure  S3: Correlation analysis of significant metabolic pathways of DEMs; Table S1: Identification of differentially expressed metabolites in positive ion mode; Table S2: Identification of differentially expressed metabolites in negative ion mode; Table S3: Identification and classification of known differentially expressed metabolites; Table S4: The KEGG pathways of differentially expressed metabolites.
Author Contributions: J.G., G.Z. and H.Z. conceived and designed the study. J.G. performed the experiments, analyzed the data, prepared the figures and tables, and wrote the manuscript. Y.L. and Y.S. contributed reagents, materials, and analytical tools. G.Z. revised the manuscript. All authors have read and agreed to the published version of the manuscript.