W931A and W931B were cultivated under consistent conditions in the test field at the Crop Research Institute, Anhui Academy of Agricultural Sciences (Hefei, Anhui, China). At the UM stage (when flower buds were 1.0-1.5 mm) and the BP stage (1.5-2.0 mm flower buds) (Li et al., 2005), anthers were collected and frozen in liquid nitrogen, then stored at -80°C prior to further experiments. There were three biological replicates per sample. W931B at the UM stage (W931B_UM) and at the BP stage (W931B_BP) were treated as the controls, whereas W931A at the UM stage (W931A_UM) and at the BP stage (W931A_BP) were considered the experimental samples.

Anthers from the BP and mature pollen (MP) stages were fixed in 2.5% (w/v) glutaraldehyde in 0.1 mol/L phosphate-buffered saline (PBS) at pH 7.4. After washing twice in PBS, samples were post-fixed in 1.0% OsO4, dehydrated in an ethanol series, and embedded in Epon812 resin as previously described (Fang et al., 2021). Embedded samples were cut into 2-µm-thick sections with an ultramicrotome (LEICA, Germany). Specimens were stained with uranyl acetate followed by lead citrate, then photographed under an H-7650 microscope (Hitachi, Japan) (Ning et al., 2018).

Total RNA was isolated from three biological replicates of the W931A and W931B samples using RNAiso Plus reagent (TaKaRa, Japan), then treated with RNase-free DNase I (Promega, Madison, WI, USA). Subsequently, 0.5µg RNA per sample was used for first-strand cDNA synthesis using a HiScript II 1st Strand cDNA Synthesis kit (Vazyme, Nanjing, China). qRT−PCR was performed in technical triplicate on a LightCycler 96 (Roche, Switzerland) using the SYBR Green Premix Ex Taq™ II quantitative PCR system (TaKaRa, Japan). The amplification program consisted of a denaturation step at 95°C for 15 min; 40 cycles of 95°C for 10 s and 60°C for 20 s; and 72°C for 30 s. Gene expression was normalized using the 2^−ΔΔCT method (Livak and Schmittgen, 2001) using TUBULIN1 as the internal reference gene. All primer sequences are listed in Supplementary Table S1 .

RNA samples were sequenced on an Illumina HiSeq 2000 platform, generating 100-bp paired-end reads. The quality of each sample file was measured using FastQC v0.11.9. Illumina adapters, reads containing Ns (indicating that the corresponding bases could not be determined), and low-quality reads (those with Qphred values ≤ 20 for more than 50% of the total read length) were removed using Trim-galore v0.6.6, resulting in a set of clean reads. The reference genome file (Gmax_275_v2.0.fa) and the reference annotation file (Gmax_275_Wm82.a2.v1.gene.gtf) were downloaded from Phytozome v12. Clean reads were mapped to the soybean reference genome using Hisat2 v2.2.1. After mapping, reads were counted with HTSeq v0.13.5. Differential expression analysis was performed in the ‘DESeq2’ (v1.32.0) R package using the following screening criteria: |log2(Fold Change)| ≥ 1; p < 0.05; and false discovery rate (FDR) < 0.05. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) biochemical pathway enrichment analyses were performed on the DEGs using the ‘clusterProfiler’ (v4.0.5) package in R.

Total protein was extracted as previously described (Li et al., 2016) from three biological replicates of soybean samples at the UM and BP stages. Precipitated proteins were washed with ice-cold 90% ethanol containing 10 mM dithiothreitol. After drying at room temperature, precipitates were lysed in sodium dodecyl sulfate (SDS). Samples were centrifuged at 12,000 × g at room temperature for 10 min, then the supernatant was collected. The sample extraction procedure was repeated, and the resulting supernatant for each sample was combined with the supernatant from the first extraction to generate the total protein solution. Bicinchoninic acid reagent (Thermo Fisher Scientific, San Jose, CA, USA) was used to measure the protein content.

Based on the calculated protein content, 100 μg of total protein was removed from each sample and digested with Trypsin Gold (Promega, Madison, WI, USA) using a protein:trypsin ratio of 30:1. Samples were incubated at 37°C for 16 h. After trypsin digestion, peptides were labeled with iTRAQ tags following the manufacturer’s protocol for the 8-plex iTRAQ reagent (Applied Biosystems, Foster City, CA, USA). The labeled peptide mixtures were then pooled and pre-separated with strong cation exchange chromatography using the LC-20AB high-performance liquid chromatography (HPLC) pump system (Shimadzu, Kyoto, Japan). iTRAQ analysis was performed on a TripleTOF 5600 system (AB SCIEX, Concord, ON, Canada) combined with a Famos autosampler (LC Packings) and an LC20-AD Nano HPLC instrument (Shimadzu) as previously reported (Du et al., 2014). Proteome Discoverer v1.2.0.339 (Thermo Fisher Scientific) was used to transform the raw data into MGF files. Protein identification and quantitation were performed in Mascot v2.3.0 (Matrix Science, London, UK) to compare the iTRAQ data against data in Soybase (https://soybase.org/). Differentially expressed proteins (DEPs) were classified as those with p < 0.05 in a least significant difference test, FDR < 0.05, and > 1.5-fold change between W931A and W931B.

DEP interaction analysis was conducted using the STRING database (https://string-db.org/). Both known and predicted protein-protein interactions were included. The option “Multiple sequences” was selected to analyze interactions among 53 differentially expressed soybean proteins of interest. Protein sequences were downloaded from Soybase.

W931A is an excellent male sterile line with a strong stem that produces more than twice as many floral buds as regular fertile varieties (Ren, 2005). To determine the cytological basis of functional deficiencies at the BP stage, we compared the ultramicroscopic structures of developing anthers and pollen grains between W931A and its maintainer line, W931B, using TEM. The results showed that the anthers of W931A and W931B developed similarly; for example, the tapetal layers of both lines were degraded at the BP stage ( Figures 1A, D ).

At the BP stage, W931B pollen grains were filled with condensed cytoplasm. Starch granules and lipid materials were extensively accumulated, and the protectum was smooth ( Figures 1A–C ). In contrast, W931A pollen grains showed an irregular shape and abnormal development (e.g., few starch granules, large vacuoles, irregularly thickened intine, and uneven protectum) ( Figures 1D–F ). At the MP stage, W931B pollen grains were round and exhibited a concentrated cytoplasm, in which a number of elaioplasts and abundant lipid compounds could be observed ( Figures 1G–I ). Comparatively, W931A pollen grains showed severe degeneration and vacuolization of the cytoplasm ( Figures 1J–L ).

To reveal the molecular mechanism underlying CMS, we analyzed the transcriptome and the proteome of developing W931A and W931B pollen grains at the UM and BP stages. In total, 12 cDNA libraries were constructed and sequenced (three biological replicates of two genotypes at two stages). Clean reads were mapped to the Wm82 reference genome, and approximately 468 million reads were aligned, with an average mapping rate of 86.48% per sample. The average GC content was 44.18% ( Supplementary Table S2 ).

The transcriptomes were then compared between W931A and W931B pollen grains at the UM and BP stages. A total of 2,259 differentially expressed genes (DEGs) were identified at the UM stage, including 598 genes that were up-regulated and 1,661 that were down-regulated in W931A ( Supplementary Figure S1A , Table S3 ). In contrast, a total of 1,962 DEGs were identified at the BP stage, including 1,063 up-regulated and 899 down-regulated genes ( Supplementary Figure S1B , Table S4 ). The expression levels of 15 randomly-selected genes were assessed using qRT-PCR to validate the sequencing results. The coefficient of determination (R²) for the two data types was near 1 (R² = 0.9304), indicating that the RNA-seq data were robust and suitable for further analysis ( Supplementary Figure S2 ).

To gain insight into the functions of the DEGs, GO term enrichment (p < 0.05) was assessed, including GO biological process (GBP), cellular component (CC), and molecular function (MF) terms. At the UM stage, GO analysis of the 2,259 DEGs revealed significant enrichment of the GBP terms “pectin catabolic process”, “photosynthesis, light harvesting in photosystem I”, “cell wall modification”, “regulation of pH”, and “response to light stimulus” ( Supplementary Table S5 ). At the BP stage, the 1,962 DEGs were significantly enriched in the GBP terms “pectin catabolic process”, “photosynthesis, light harvesting in photosystem I”, “cell wall modification”, “response to light stimulus”, and “protein-chromophore linkage” ( Supplementary Table S6 ). KEGG biochemical pathway enrichment (p < 0.05) was also investigated, and revealed enrichment of pathways including “pentose and glucuronate interconversions”, “photosynthesis-antenna proteins”, “metabolic pathways”, “plant-pathogen interaction”, and “ribosome biogenesis in eukaryotes” at the UM stage ( Supplementary Table S8 ), and “metabolic pathways”, “biosynthesis of secondary metabolites”, “flavonoid biosynthesis”, “photosynthesis-antenna proteins”, and “tyrosine metabolism” at the BP stage ( Supplementary Table S9 ).

Further analysis showed that a total of 865 DEGs were shared between the two developmental stages ( Figure 2A ). We also found that most of the DEGs were down-regulated; at the UM stage, there were 10-fold more down-regulated than up-regulated genes; at the BP stage, there were approximately five-fold more down-regulated than up-regulated genes ( Figure 2B ). To determine the putative functions of the shared DEGs, we conducted a GO enrichment analysis. The 865 DEGs were primarily associated with the GBP terms “pectin catabolic process”, “cell wall modification”, “photosynthesis, light harvesting in photosystem I”, “response to light stimulus”, “protein-chromophore linkage”, and, notably, “pollen development” ( Supplementary Table S7 ). The most significantly enriched MF terms were “pectinesterase activity” and “pectinesterase inhibitor activity” ( Figure 2C , Supplementary Table S7 ). Many of the shared DEGs that had those MF annotations encoded pectinesterases (PMEs)/PME inhibitors and pectate lyases. For example, Rfk1 encodes a PME/PME inhibitor protein that is involved in pectin metabolism and regulates fertility restoration in wheat (Chen et al., 2021). Glyma.09G042200 encoded a PME/PME inhibitor enzyme that was strongly down-regulated in W931A at both the UM and BP stages. qRT-PCR results also showed marked down-regulation of three PME/PME inhibitor genes in W931A compared to W931B ( Figure 3 ).

In W931A, 1,394 and 1,097 DEGs were uniquely expressed at the UM and BP stages, respectively ( Supplementary Tables S3, S4 ). GO term analysis of the 1,394 unique DEGs at the UM stage revealed enrichment of the GBP terms “mitochondrion organization”, “hydrogen peroxide catabolic process”, “positive regulation of mitochondrial translation”, “response to oxidative stress”, and “cell wall modification” ( Supplementary Table S5 ), whereas the 1,097 DEGs unique to the BP stage were significantly enriched in the GBP terms “auxin-activated signaling pathway”, “defense response to other organism”, “flavonoid biosynthetic process”, “fatty acid biosynthetic process”, and “brassinosteroid metabolic process” ( Supplementary Table S6 ). These data indicated that many genes expressed in the pollen had unique spatiotemporal expression patterns.

We next carried out a proteomic analysis to identify factors involved in soybean CMS at the protein level. To analyze DEPs between W931A and W931B, proteins expressed in the anthers of W931A and W931B plants at the UM and BP stage were analyzed using iTRAQ technology. A total of 630 DEPs were identified in W931A compared to W931B at the UM stage, with 305 up-regulated and 325 down-regulated ( Figure 4A ; Supplementary Table S11 ). At the BP stage, 242 proteins were up-regulated in W931A compared to W931B, whereas 384 were down-regulated ( Figure 4A ; Supplementary Table S12 ).

GO term enrichment analysis revealed that the GBP terms “glycolytic process”, “one-carbon metabolic process”, “carbohydrate metabolic process”, “response to toxic substance”, and “cell wall modification” were significantly enriched in W931A at the UM stage ( Supplementary Table S14 ). At the BP stage, markedly enriched GBP terms included “glycolytic process”, “glucose metabolic process”, “one-carbon metabolic process”, “pectin catabolic process”, and “cell wall modification” ( Supplementary Table S15 ). However, several key molecular functions were also enriched at the UM stage, such as “nuclear import signal receptor activity”, “glutathione transferase activity” and “NAD binding”, indicating that enzymes with these functions tended to be more active during UM development ( Supplementary Table S14 ). In contrast, GO terms that were specifically enriched at the BP stage included “electron carrier activity”, “structural constituent of ribosome”, “mRNA binding”, and “phosphatidic acid binding” ( Supplementary Table S15 ), indicating that proteins associated with these MF terms played important roles at the BP stage.

Further analysis revealed 343 DEPs in common between the two plant developmental stages ( Figure 4B ; Supplementary Table S13 ). Enriched GO terms in the shared DEP set included “glycolytic process”, “cell wall modification”, “pectin catabolic process”, “glucose metabolic process”, and “carbohydrate metabolic process” ( Supplementary Table S16 ). The largest category of enriched MF GO terms was “pectinesterase activity” ( Supplementary Table S16 ). The shared DEPs were mapped to a total of 31 KEGG pathways ( Supplementary Table S18 ). The most significantly enriched pathway was “metabolic pathways”, followed by “carbon metabolism” and “glycolysis/gluconeogenesis”. The primary physiological function of carbohydrate metabolism in an organism is to provide carbon sources. Some proteins related to monosaccharide and polysaccharide metabolic pathways were also evidently expressed; proteins with the pathway annotations “galactose metabolism”, “fructose and mannose metabolism”, and “starch and sucrose metabolism” were enriched in W931A ( Supplementary Table S18 ).

To identify whether there was general concordance between mRNA and associated protein levels during pollen development, the transcriptomic and proteomic data were compared by conducting a correlation analysis between the DEGs and DEPs. One hundred and forty-four common-DEGs and DEPs were expressed at the UM stage ( Supplementary Tables S3, S11 ). Seventy-eight common-DEGs and DEPs were expressed at the BP stage ( Supplementary Tables S4, S12 ). The expression profiles of 53 sets of DEGs and associated DEPs found at both the UM and BP stages are shown in Supplementary Figure S3 and Table 1 . Further analysis of the 53 shared DEGs with similar corresponding protein expression showed that they primarily encoded PMEs/PME inhibitors, polygalacturonases (PGs), peroxidases, ATPases, and sugar transport proteins.

A total set of 53 DEGs/DEPs that were differentially expressed at both the UM and BP stages were identified as genes/proteins that affect biological processes during soybean pollen development. GO term analysis of this set showed marked enrichment of the GBP terms “cell wall modification”, “pectin catabolic process”, “carbohydrate metabolic process”, “hydrogen peroxide catabolic process”, and “response to oxidative stress” ( Supplementary Table S19 ). The MF terms “pectinesterase activity”, “carboxylic ester hydrolase activity”, “polygalacturonase activity”, “carbohydrate binding”, and “peroxidase activity” were also enriched ( Supplementary Table S19 ). Analysis of the transcriptomic and proteomic data demonstrated consistent activation of PMEs/PME inhibitors at the mRNA and protein levels, indicating the importance of these proteins in regulating soybean fertility. In addition, the KEGG pathways “pentose and glucuronate interconversions”, “metabolic pathways”, “biosynthesis of amino acids”, “biosynthesis of secondary metabolites”, and “glycolysis/gluconeogenesis” were significantly enriched ( Figure 5A ; Supplementary Table S20 ).

Pollen development is an intensely energy-consuming process (Hanson and Bentolila, 2004; Linke and Börner, 2005; Geng et al., 2018). The process of carbon fixation is primarily responsible for generating usable carbon for plants (Saxena et al., 2020). Here, transcriptomic analysis revealed 30 DEGs in pollen that participated in pentose and glucuronate interconversions ( Supplementary Table S10 ) proteomic analysis showed that there were 26 DEPs associated with carbon metabolism, 19 related to glycolysis/gluconeogenesis, 10 related to pentose and glucuronate interconversions, six associated with galactose metabolism, six related to fructose and mannose metabolism, eight associated with starch and sucrose metabolism, and four related to the pentose phosphate pathway ( Supplementary Table S18 ).

We also found that five DEPs were involved in the tricarboxylic acid (TCA) cycle, with three of these proteins up-regulated in W931A and two down-regulated ( Supplementary Table S17 ). The two down-regulated DEPs associated with the TCA cycle were phosphoenolpyruvate carboxykinase (ATP) proteins. We also found that seven DEPs were associated with oxidative phosphorylation, with five up-regulated and two down-regulated in W931A compared to W931B ( Supplementary Table S17 ). The two down-regulated proteins were components of ATPase 9.

MADS-box genes encode TFs that have important biological functions in the formation of male flower organs (Linke et al., 2003). In Arabidopsis thaliana, AGL104 (At1g22130) is expressed specifically in the pollen (Adamczyk and Fernandez, 2009). Two agamous-like MADS-box (AGL) genes with the GBP “pollen development” were markedly expressed at both the UM and BP stages ( Supplementary Table S7 ). Glyma.13G086400 (log₂(FC) = -5.83) and Glyma.14G168700 (log₂(FC) = -5.99) were clearly down-regulated in W931A at the UM stage ( Supplementary Table S3 ) and similarly down-regulated at the BP stage (log₂(FC) = -2.09 and -2.48, respectively) ( Supplementary Table S4 ). AGL104 expression levels changed during anther abortion, suggesting that AGL104 may have an important molecular function in pollen production. Moreover, the Arabidopsis gene tyrosine-phosphatase 1 (PTEN1) encodes a pollen-specific phosphatase that is essential for pollen development (Gupta et al., 2002). Phospholipase A2 (PLA2) plays critical roles in Arabidopsis pollen development, germination, and tube growth (Kim et al., 2011). In the present study, two genes with the GBP “pollen development”, Glyma.08G241500 (PTEN1) and Glyma.07G127900 (PLA2-like), were also down-regulated at both the UM and BP stages ( Supplementary Tables S3, S4 ). They were significantly expressed at the transcriptional level, but not at the protein level, indicating that a complex post-transcriptional regulatory network may exist in W931A to regulate male sterility.

Pectin is one of the most important polysaccharides for plant cell wall growth (Zhang et al., 2020). We found that pectin catabolic processes occurred at both the transcriptional and protein levels. In W931A at the UM stage, 39 DEGs (two up-regulated and 37 down-regulated) had annotations associated with pectin catabolic processes ( Supplementary Table S5 ). Similarly, there were 34 DEGs (seven up-regulated and 27 down-regulated) in W931A at the BP stage encoded pectate lyases or PMEs ( Supplementary Table S6 ), and seven DEPs (all down-regulated) at the UM stage ( Supplementary Table S14 ), and 10 DEPs (all down-regulated) at the BP stage that were pectate lyases or PMEs ( Supplementary Table S15 ). Five PMEs were also identified in the set of 53 shared DEGs and DEPs; these were down-regulated at both the mRNA and protein levels ( Figure 5B ). Furthermore, among the 53 candidate DEGs were three DEGs encoding PGs (Glyma.06G207300, Glyma.14G006100, and Glyma.04G143200), three DEGs encoding peroxidases (Glyma.02G037800, Glyma.01G027000, and Glyma.10G198100), one DEG encoding tubulin alpha-3 (Glyma.16G154000), and one DEG encoding ATPase 9 (Glyma.15G004300) ( Table 1 ).

Interactions between the candidate proteins corresponding to the 53 shared DEGs and DEPs were next predicted using STRING 10.0. The results indicated that these proteins may form a complex network regulating pollen development in W931 CMS ( Figure 6 ). The putative protein regulatory network included PME, β-galactosidase, actin-depolymerizing factor, and S-adenosylmethionine synthetases as core proteins involved in pollen development ( Figure 6 ). There was a potential interaction mechanism between four PME proteins (PME28, PME58, and two PME1s) and β-galactosidase 15 (BGAL 15). The PMEs and BGAL 15 may be involved in the cell wall loosening associated with pollen expansion during microspore development. In okra (Abelmoschus esculentus), AeADF1 has been shown to promote pollen germination and pollen tube growth (Dong et al., 2021). Here, actin depolymerizing factor 7 (GmADF7) and its interacting proteins were significantly down-regulated in W931A, indicating that they may inhibit pollen germination and thus cause soybean sterility.