*3.3. Blocked Carbohydrate Metabolism and Sugar Transport Leads to Abnormal Pollen Development in SXCMS5A*

Carbohydrate metabolism and sugar transport were the most basic metabolic processes in plant, providing energy and carbon for anther development, and starch and sucrose serve as energy reserves for pollen maturation [34,35]. In this study, starch and sucrose metabolism, pentose and glucuronate interconversions, and glycolysis/gluconeogenesis were three enriched pathways for carbohydrate metabolism. These 3 key pathways contained 58 DEGs of which 55 were downregulated in SXCMS5A. Male sterility in the soybean CMS line NJCMS1A was associated with these three pathways [21]. Similarly, we found that many sugar transport-related DEGs, such as *STP11* and *sucrose transport protein SUC4-like* (*SUC4*), were downregulated in SXCMS5A. In cucumber, Sun et al. [16] found that downregulation of the sugar transporter *CsSUT1* inhibited pollen germination and caused male sterility. In soybean CMS-based F1, Ding et al. [36] had found downregulation of sugar transport-related DEGs and reduction of sugar accompanied by the decrease of male fertility under heat stress. Furthermore, substance amount analysis also showed that energy supply was decreased in SXCMS5A compared to SXCMS5B. Thus, the downregulated expression of carbohydrate metabolism and sugar transport related genes might lead to insufficient energy supply, which was related to pollen abortion and male sterility in SXCMS5A.

#### *3.4. Abnormal ROS Metabolism Leads to Pollen Abortion in SXCMS5A*

PCD was a common phenomenon in the development of animals and plants, regulated by genes under specific physiological or pathological conditions [37]. Tapetum provides nutrients for pollen development, and its abnormal PCD process is one of the direct causes of plant male sterility [14,15,38,39]. In this study, cytological analysis showed that the tapetum was vacuolated and abnormally developed in SXCMS5A compared to SXCMS5B, which had typical morphological characteristics of abnormal PCD. Although the male-sterile lines formed microspores through meiosis, they could not provide related materials for the development of pollen and could not ultimately form functional pollen. Most importantly, there was a close relationship between PCD and ROS metabolism [40,41]. Studies had shown that abnormal ROS metabolism during anther or spikelet development was related to male sterility [42–44]. Ascorbic acid and glutathione had important physiological functions in plants, with special roles in maintaining the redox balance of cells in the plant antioxidant system [45,46]. In this study, *L-AO* and *GST* genes (components of ascorbic acid and glutathione metabolism, respectively) were downregulated in flower buds of the soybean CMS lines compared to their maintainer lines. Furthermore, enzyme activity analysis also showed that ROS scavenging were decreased in SXCMS5A compared to SXCMS5B. Thus, the downregulated expression of ROS metabolism genes might lead to abnormal PCD, which was related to pollen abortion and male sterility in SXCMS5A.

#### *3.5. Abnormal Expression of TF Related DEGs Causes Pollen Abortion in SXCMS5A*

TFs are important regulators of gene expression, and their expression changes may have an important impact on plant growth and development [47,48]. For example, MYB, bHLH and WRKY participated in regulation of rates of gene transcription and regulation of meiosis, which was very important for stamen development and maturation [49–52]. Most of these TFs played a key role in the process of tapetum PCD and pollen formation, and their abnormal functioning often caused male sterility [13,53,54]. In this study, 20 coding TF DEGs were found between SXCMS5A and SXCMS5B among which 15 downregulated and 5 upregulated in SXCMS5A compared to SXCMS5B. Among these DEGs, 6 were related to MYB TFs, 2 were related to bHLH TFs, and 1 was related to WRKY TFs. In addition, effective activation of the ethylene signaling pathway was required for plant responses to growth and environmental signals, but continuous over activation of the ethylene signaling pathway had obvious inhibitory and toxic effects on plant growth and reproduction [55]. Previous research had shown that increasing ethylene concentration could lead to male sterility in wheat [55]. Up-regulation of the *Glyma.16G164800* encoding ethylene response

transcription factor was also consistent with this result. Thus, the abnormal expression of these TFs might cause disturbance in expression of genes related to tapetum and pollen development, which was related to pollen abortion and male sterility in SXCMS5A.

#### *3.6. Proposed Model for the Mechanism of Male Sterility in Soybean CMS Line SXCMS5A*

According to previous reports and the data presented in this study, we made the following speculation on the mechanism of male sterility in soybean CMS line SXCMS5A (Figure 9). First, the rearrangement of soybean CMS line SXCMS5A mitochondrial genome generated CMS genes, including *ORF178* and *ORF103c*. The production of *ORF178* or *ORF103c* leads to mitochondrial dysfunction, such as blocked energy synthesis and massive production of ROS. Subsequently, mitochondrial defects directly/indirectly lead to the down-regulation of genes related to carbohydrate metabolism, sugar transport and pollen wall development in the nucleus, leading to further energy shortage and abnormal pollen development during pollen development. In addition, the downregulation of enzymatic ROS scavenging related genes in the nucleus leads to the dysfunction of the enzymatic ROS scavenging system, resulting in the inability of effective ROS clearance and the accumulation of ROS, which affects the normal PCD process of anther tapetum. The combination of these processes eventually leads to male sterility in soybean CMS line SXCMS5A. Further studies are needed to validate this proposed model. *Int. J. Mol. Sci.* **2022**, *23*, x FOR PEER REVIEW 14 of 20

**Figure 9.** A proposed model for the mechanism of male sterility in soybean CMS line SXCMS5A. The upregulated and downregulated genes or metabolite contents are in red and green **Figure 9.** A proposed model for the mechanism of male sterility in soybean CMS line SXCMS5A. The upregulated and downregulated genes or metabolite contents are in red and green backgrounds, respectively.

#### backgrounds, respectively. **4. Materials and Methods**

#### **4. Materials and Methods**  *4.1. Plant Materials and Sample Collection*

*4.1. Plant Materials and Sample Collection*  Three soybean CMS lines, SXCMS5A, SXCMS6A, and SXCMS7A with their respective maintainer lines were used in the present study. SXCMS5A was developed by continuous backcross with the CMS line H3A as the donor parent and the variety JY20 (designated as SXCMS5B afterwards) as the recurrent parent; SXCMS6A was developed Three soybean CMS lines, SXCMS5A, SXCMS6A, and SXCMS7A with their respective maintainer lines were used in the present study. SXCMS5A was developed by continuous backcross with the CMS line H3A as the donor parent and the variety JY20 (designated as SXCMS5B afterwards) as the recurrent parent; SXCMS6A was developed by continuous backcross with the CMS line H3A as the donor parent and the strain LX11 (designated as

by continuous backcross with the CMS line H3A as the donor parent and the strain LX11 (designated as SXCMS6B afterwards) as the recurrent parent; and SXCMS7A was

line H3A was developed by continuous backcross with the CMS line JLCMS1A as the donor parent and the strain H3 as the recurrent parent, whereas JLCMS1A was

SXCMS5A and SXCMS5B were planted in the summer of 2017 at Dangtu Experimental Station of Nanjing Agricultural University. Because it was difficult to judge the precise pollen development stage from flower bud appearance in soybean, the mixture of flower buds with different sizes were collected from three individual plants in the afternoon as three biological replicates for SXCMS5A and SXCMS5B. The samples were immediately placed in liquid nitrogen and then stored at −80 °C for RNA-Seq.

SXCMS5A and SXCMS5B, SXCMS6A and SXCMS6B, and SXCMS7A and SXCMS7B were planted in the summer of 2019 at Dangtu Experimental Station. The mixture of flower buds with different sizes were collected in the afternoon during flowering period, immediately placed in liquid nitrogen, then stored at −80 °C for qRT-PCR. All qRT-PCR

SXCMS5A and SXCMS5B were planted in the spring of 2019 at Dongyang Experimental Station of Shanxi Agricultural University. Flower buds with different sizes

introduced from Jilin Academy of Agricultural Sciences.

reactions were performed with three biological replicates.

SXCMS6B afterwards) as the recurrent parent; and SXCMS7A was developed by continuous backcross with the CMS line H3A as the donor parent and the strain JDX (designated as SXCMS7B afterwards) as the recurrent parent. Here, the CMS line H3A was developed by continuous backcross with the CMS line JLCMS1A as the donor parent and the strain H3 as the recurrent parent, whereas JLCMS1A was introduced from Jilin Academy of Agricultural Sciences.

SXCMS5A and SXCMS5B were planted in the summer of 2017 at Dangtu Experimental Station of Nanjing Agricultural University. Because it was difficult to judge the precise pollen development stage from flower bud appearance in soybean, the mixture of flower buds with different sizes were collected from three individual plants in the afternoon as three biological replicates for SXCMS5A and SXCMS5B. The samples were immediately placed in liquid nitrogen and then stored at −80 ◦C for RNA-Seq.

SXCMS5A and SXCMS5B, SXCMS6A and SXCMS6B, and SXCMS7A and SXCMS7B were planted in the summer of 2019 at Dangtu Experimental Station. The mixture of flower buds with different sizes were collected in the afternoon during flowering period, immediately placed in liquid nitrogen, then stored at −80 ◦C for qRT-PCR. All qRT-PCR reactions were performed with three biological replicates.

SXCMS5A and SXCMS5B were planted in the spring of 2019 at Dongyang Experimental Station of Shanxi Agricultural University. Flower buds with different sizes were collected in the afternoon at the flowering stage and fixed in formaldehyte-alcohol-acetic acid (FAA) for cytological examination. The mixture of flower buds of different sizes was collected in the afternoon at the flowering stage, immediately placed in liquid nitrogen, then stored in at −80 ◦C for enzyme activity assay and substance content analysis. All enzyme activity assay and substance content analysis were performed with three biological replicates.

### *4.2. Cytological Examination*

To observe the cyto-morphological characteristics of pollen development of SXCMS5A and SXCMS5B, flower buds with different sizes were fixed, dehydrated, embedded, sectioned and stained according to a previous report [56]. To observe the pollen fertility of SXCMS5A and SXCMS5B, the anthers of unopened flowers (the flowers that will open in the morning of next day) in the afternoon were taken and stained with a 1% I2-KI solution [57]. All samples were observed using a light microscope (Nikon Eclipse CI, Tokyo, Japan), and photographed under the imaging system (Nikon DS-U3, Tokyo, Japan).

#### *4.3. Total RNA Extraction, Library Construction, and Sequencing*

Trizol (Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from the flower buds of SXCMS5A and SXCMS5B. The construction of cDNA library referred to prokaryote, considering that the plant mitochondrial genomes were similar to its ring genome. So, after total RNA was extracted, sample mRNA was enriched by removing rRNA by a Ribo-ZeroTM Magnetic Kit (Epicenter, Madison, WI, USA). Next, the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse-transcribed into cDNA with random primers. Second-strand cDNA was synthesized with DNA polymerase I, RNase H, dNTPs, and buffer. The cDNA fragments were then purified with a QiaQuick PCR extraction kit, end-repaired, poly(A) was added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeqTM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

#### *4.4. Raw Sequencing Data Analysis and Bioinformatics Analysis*

The raw data from the sequencing machines were initially filtered to get clean data. The short-reads alignment tool Bowtie2 [58] was used to compare and remove reads containing rRNA. Tophat2 software [59] was used to compare the reads of the filtered rRNA to the reference genome (Nucleus, Wm82.a2.v1; mitochondria, JX463295.1; chloroplast, DQ317523.1). Next, the transcripts of a group of different repeats were fused into

comprehensive transcripts with Cufflinks software [60]; transcripts of multiple groups were then merged into a group of final comprehensive transcripts for further analysis of downstream differential expression. Transmembrane domain of ORF was predicted using DeepTMHMM (V.1.0.12, https://dtu.biolib.com/DeepTMHMM, accessed on 15 September 2022) [61]. The FPKM values of DEGs between SXCMS5A and SX-CMS5B in soybean root, stem, leaf and flower tissues were obtained from the RNA-seq data in Phytozome v12.0 (https://phytozome.jgi.doe.gov/pz/portal.html#, accessed on 15 August 2020), and the heat map was conducted using MeV 4.9 software.

#### *4.5. Quantification of Gene Abundance and DEG-Analysis*

Gene abundance was quantified with RSEM software [62]. The gene expression level was normalized with fragments per kilobase of transcript per million mapped reads (FPKM) method. To identify differentially expressed genes, the edgeR package was used. A standard of "FDR < 0.05 and |log2FC| > 1" was used as the threshold to screen for significant DEGs. DEGs were then subjected to enrichment analysis of GO functions and KEGG pathways.

#### *4.6. GO and KEGG Pathway Enrichment Analysis*

GO enrichment analysis identified all GO terms that were significantly enriched in DEGs comparing to the genome background, and filtered the DEGs that correspond to biological functions. All DEGs were mapped to GO terms in the gene ontology database (http://www.geneontology.org/, accessed on 15 December 2017); gene numbers were calculated for every term, and significantly enriched GO terms in DEGs compared to the genome background were defined by hypergeometric test. The calculated *p* values were run through FDR correction, taking a Q value ≤ 0.05 as a threshold. GO terms meeting this condition were defined as significantly enriched GO terms in DEGs. This analysis enabled identification of the main biological functions correlated with the DEGs in question.

KEGG was the major public pathway-related database [63]. Pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways in DEGs compared to the whole genome background. The calculated *p* value was run through FDR correction, taking a Q value ≤ 0.05 as a threshold value. Pathways meeting this condition were defined as significantly enriched pathways in DEGs.

#### *4.7. qRT-PCR Analysis*

qRT-PCR was used to validate the gene expression levels of DEGs detected by RNA-Seq. All primers (Table S6) were designed based on the mRNA sequences, and synthesized commercially (General Biosystems, Chuzhou, China). Total RNAs from the flower buds of SXCMS5A and SXCMS5B, SXCMS6A and SXCMS6B, SXCMS7A and SXCMS7B were used for the validation of RNA-Seq. Using the protocol provided in the HiScript Q RT SuperMix for qPCR kit (+gDNA wiper, Vazyme, Nanjing, China), 1 µg of total RNA was reversetranscribed using oligo (dT) primers. qPCR analysis was carried out using the AceQ qPCR SYBR Green Master Mix (Vazyme, Nanjing, China) on a Bio-Rad CFX96 instrument (CFX96 Touch, BIO-RAD, USA). *GmTubulin* (accession number: *NM\_001252709.2*) was used as the internal control gene [36]. The maintainer lines were used as the control of their male sterile lines. Relative expression levels of the genes were quantified using the 2 <sup>−</sup>∆∆Ct method [64].

#### *4.8. Substance Contents and Enzyme Activity Assay*

Soluble sugar and starch contents were measured by visible light spectrophotometer according to the operation procedure of soluble sugar content detection kit (Solarbio, Beijing, China) and starch content detection kit (Solarbio, Beijing, China), respectively. ATP content was measured by UV spectrophotometer according to the operation procedures of ATP content detection kit (Solarbio, Beijing, China).

CAT and AAO activities were measured by UV spectrophotometer according to the operation procedure of CAT activity test kit (Solarbio, Beijing, China) and AAO activity test kit (Solarbio, Beijing, China), respectively. GPX activity was measured by visible light spectrophotometer according to the operation procedure of GPX activity test kit (Solarbio, Beijing, China).

#### **5. Conclusions**

In this study, two ORFs in mitochondria, including *ORF178* and *ORF103c*, were upregulated in sterile lines and had transmembrane domains, which were identified as two candidate CMS genes of soybean CMS line SXCMS5A as well as its two half-sib sister lines with a same cytoplasm source (SXCMS6A and SXCMS7A). Our study showed that pollen wall development, carbohydrate metabolism, sugar transport, ROS metabolism related genes and TFs were involved in the process of pollen abortion and male sterility. The male sterility mechanism of SXCMS5A might be the rearrangement of soybean mitochondrial genome to produce CMS gene, which directly or indirectly affected a series of biological processes, such as the decrease of energy supply and the outbreak of ROS, leading to the abnormal development of anther tapetum and finally pollen abortion. Future research will focus on cloning CMS related candidate genes in soybean.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms232012227/s1.

**Author Contributions:** J.G. and S.Y. designed the experiments. Z.B., X.D., R.Z., Y.Y. and B.W. performed the experiments. Z.B. and X.D. wrote the manuscript. Z.B., X.D., S.Y. and J.G. revised the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by grants from the National Key R&D Program of China (2016YFD0101500, 2016YFD0101504).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The datasets generated by this study can be found in the NCBI using accession number PRJNA887481.

**Acknowledgments:** We thank Huan Sun (Jilin Academy of Agricultural Sciences, China) for kindly providing seeds of the soybean CMS line JLCMS1A, and Baoguo Wei (Center for Agricultural Genetic Resources Research, Shanxi Agricultural University, China) for providing seeds of the soybean CMS line H3A derived from JLCMS1A for this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

