Phenotypic characterization of ms6 mutant
The offspring of heterozygous ms6 plants (T295H) were planted in a greenhouse condition. During vegetative stage, all plants grew well just like wild types (WT), while in reproductive stage about a quarter of plants (ms6) were male-sterile and unable to develop pods after blooming (Fig. 1a). Compared to the anthers in WT plants, the anthers of ms6 plants were more whitish and shrinking (Fig. 1b). Pollen grains released from WT anthers were round and turned dark blue after stained by I2-KI solution (Fig. 1c), but no pollen grains were produced in ms6 sterile anthers (Fig. 1d). These results are consistent with the previous report (Skorupska and Palmer 1989).
WT and ms6 flowers at late tetrad and late microspore stages were cross-sectioned and observed under light microscopy. At late tetrad stage, WT anther wall was composed of 5-layers, including epidermis, endothecium, middle layer, parietal layer, and tapetum from outside to inside, in which the cytoplasm of tapetum cells is highly condensed; meanwhile, callose surrounding the tetrads in locule appeared to start degeneration (Fig. 1e). Anthers at same stage in Arabidopsis and rice show similar cytological features except that their anther walls lack of the parietal layer (Sanders et al. 1999; Zhang et al. 2011). On the other hand, ms6 anther had radically enlarged and highly vacuolated parietal and tapetal layers; in the locule, callose encompassed partially or non-separated irregular microspores with multiple nuclei, indicating that cytokinesis II of meiocytes was abnormal (Fig. 1f). This is similar to the phenotype of ms6 (Ames2) (T354H), but tapetal layer of ms6 (Ames2) was degenerated more rapidly and was almost completely degraded at this stage (Ilarslan et al. 1999). By the time of early pollen stage, enlarged pollens with thick walls were observed in the locule of WT anther, and the anther wall was composed by an epidermis, an enlarged endothecium, and a narrow parietal layer with attachment of remnant tissue from degraded tapetum cells (Fig. 1g). Comparatively, in ms6 anther, tapetum cells were completely dissolved and pollens were crushed, while parietal layers were abnormally vacuolated and swollen, similar to the situation reported in ms6 (Ames2) (T354H) (Fig. 1h; Ilarslan et al. 1999).
A SNP mutation of Glyma.13g066600, an R2R3-MYB transcription factor encoded gene, is likely responsible for the male sterility in ms6
The ms6 loci was mapped previously within a 3.72 Mb region on Chr13 between two SSR markers Satt149 (Chr13:13,134,055 bp) and Satt030 (Chr13:16,855,019 bp) (Yang et al. 2014). To narrow down the region, WT and ms6 bulks were constructed from a BC5F2 mapping population derived from T295H and ‘JiuB’, and subject to BSA-seq analysis. Plotting the ∆SNP-index values between two bulks against their genomic positions showed that ms6 was associated with a 1.5 Mb region (15,853,267 − 17,349,424 bp) on Chr13 (Fig. 2a), consistent to the previously reported interval (Yang et al. 2014). In addition, 249 variations between two bulks were identified in this region, including 214 SNP and 35 InDel.
Subsequently, a fine mapping was conducted with 328 individual plants in the BC5F2 mapping population with 9 polymorphic SSR markers identified in the 1.5 Mb interval, including BARCSOYSSR-13-0243, BARCSOYSSR-13-0244, BARCSOYSSR-13-0245, BARCSOYSSR-13-0249, BARCSOYSSR-13-0257, BARCSOYSSR-13-0259, BARCSOYSSR-13-0275, BARCSOYSSR-13-0277, and BARCSOYSSR-13-0283 (Fig. 2b; Fig. S1). Finally, ms6 was restricted to a 255 kb region (Chr13: 16,428,596 − 16,683,664 bp) between SSR markers BARCSOYSSR-13-0259 and BARCSOYSSR-13-0275 (Fig. 2b), which harbored 23 annotated genes. According to the BSA-Seq results, amongst these genes, only Glyma.13g066600 had an SNP (T-to-A, Chr13: 16,641,429 bp) in the CDS region (Fig. 2c). The SNP in ms6 destroyed an MseI restriction site (TTAA) that was present in WT sequence (Fig. 2c).
To verify the mutation, a 126-bp region covering the SNP site was amplified by PCR from homozygous WT (+/+), heterozygotes (+/ms6), and homozygous ms6 (ms6/ms6), and subsequently digested with MseI. As expected, the amplicon of WT was cleaved to a 99-bp band and a 27-bp band that was invisible on the agarose gel. Comparatively, the amplicon from ms6 plants could not be digested at all, while about half of the amplicons from heterozygotes could be cut (Fig. 2e), evidencing the mutation of T to A at Chr13:16,641,429 bp in ms6. Moreover, the 126-bp segment could be used as a CAPS marker to differentiate the plant genotype at ms6 locus at any stage.
Glyma.13G066600 encodes a typical R2R3-MYB transcription factor with two MYB motifs close to the N terminus served as DNA binding domain (Fig. 2d). This protein showed a strong homology to TDF1, amino acid sequence exhibiting 48% identity and 59% similarity to AtTDF1 in Arabidopsis (Fig. S2). TDF1 is known as a key transcription factor in regulating tapetum development and function. Null mutant of TDF1 in Arabidopsis and rice both promoted vacuolization in tapetal cells and suppressed the degradation of the callose surrounding the tetrads, resulting in squeezed microspores and no pollen formation (Zhu et al. 2008; Cai et al. 2015), which was similar to ms6. The SNP in Glyma.13g066600 CDS in ms6 leads to the amino acid substitution of leucine to histidine at the residue 76 (L76H) (Fig. 2d), which is a well-conserved residue in the R2 MYB motif (Fig. 2f). These results suggested that mutation at Glyma.13g066600 was responsible for the male sterility of soybean ms6, and therefore its wild-type allele was termed as GmMS6.
GmMS6 is a homolog of TDF1 that is only present in angiosperm
Blast search showed that GmMS6 had a homolog GmMS6L with 92% amino acid sequence identity, encoded by Glyma.19g017900 (Supplementary Fig. S2). Similar to AtTDF1, GmMS6L doesn’t have the N30 extension in the sequence (Supplementary Fig. S2). To further confirm the relationship between GmMS6 and TDF1, we conducted a phylogenetic analysis to GmMS6, GmMS6L, and the homologs of AtTDF1 and OsTDF1 in different land plant evolutionary lineages. AtMYB80 and OsMYB80 were used as outgroup sequences for MYB80 is the closest homolog of TDF1, also known as MYB35 (Dubos et al. 2010).
The result showed that TDF1 homologs in angiosperm were clustered into a monophyletic group with two branches. All the dicots were grouped in a branch with the basal angiosperm species A. trichopoda, while all the monocots were grouped in the other branch with the basal or near basal monocot species A. comosus and Z. marina (Fig. 3). Comparatively, the TDF1 homologs with highest bit-score in lycophyte S. moellendorffii and moss P. patens were clustered with MYB80 and those in gymnosperm species, G. biloba and P. taeda, had even a further evolutionary relationship. These data suggested that TDF1 is only present in angiosperm and has diverged at a very early stage in angiosperm evolution before monocots and dicots were differentiated.
In the phylogenetic tree, GmMS6 and GmMS6L were positioned together in the TDF1 clade, evidencing that they were homologs of TDF1 and evolved from a recent duplication (Fig. 3). It is interesting to know whether L76H in GmMS6 (GmMS6L76H) could lead to male sterility, and if so, why GmMS6L fails to compensate the GmMS6’s function in ms6 mutant.
GmMS6 but not GmMS6L76H could complement the attdf1 male sterile phenotype in Arabidopsis
The function TDF1 was likely conserved in all the angiosperm species. Firstly, phylogenetic results showed that TDF1 in angiosperm were clustered into a monophyletic group (Fig. 3). Secondly, depletion of TDF1 in Arabidopsis and rice caused similar detrimental effects on tapetum development as lack of GmMS6 did (Zhu et al. 2008; Cai et al. 2015). Therefore, we speculated that GmMS6 could substitute the role of AtTDF1 in Arabidopsis although GmMS6 is 30 amino acids longer at the N terminus, while GmMS6L76H could not.
To verify this, we firstly transformed heterozygous attdf1 (+/attdf1) plants with the CDSs of AtTDF1, GmMS6, GmMS6L, and GmMS6L76H driven by native promoter AtTDF1p. T1 transgenic plants in homozygous attdf1 background were subsequently assessed for anther fertility. As it turned out, most attdf1 transformed with AtTDF1 (14/20, rescued/transformants), GmMS6 (19/26, rescued/transformants) and GmMS6L (8/10, rescued/transformants) were fully complemented, producing functional pollens and elongated siliques as wild type did (Fig. 6a-e, 6h-l, 6o-p), and demonstrating that both GmMS6 and GmMS6L are functional TDF1. Therefore, we termed these two proteins as GmTDF1a and GmTDF1b, correspondingly. By contrast, AtTDF1p:GmMS6L76H (0/25, rescued/ transformants) failed to complement the attdf1’s sterility (Fig. 6a-b, 6f, 6h-i, 6m, 6r), confirming that GmMS6L76H was a malfunctioned protein and responsible for the aberrant male development of ms6 mutant.
Additionally, expression of AtTDF1p-driven AtTDF1 L46H, the AtTDF1 mutant corresponding to GmMS6L76H, in attdf1 displayed consistent result. The fertilities of transgenic AtTDF1 L46H/attdf1 plants could not be restored because all the transformants (0/16, rescued/ transformants) developed just like attdf1(Fig. 6a-b, 6g, 6h-i, 6n, 6s), indicating that leucine at this position (46th in AtTDF1 and 76th in GmTDF1) is crucial for the function of TDF1.
GmMS6/GmTDF1a is the major functional TDF1 in soybean
Complementation assay above exposed that GmMS6 and GmMS6L were both functional proteins, which rose up the question why ms6 would exhibit male sterility when there is another TDF1 coding gene in the genome. To answer that, we analyzed the expression patterns of GmMS6 and GmMS6L by qRT-PCR in roots, stems, leaves, young flowers, siliques, and immature seeds. GmMS6 displayed a typical tissue-specific expression pattern, with a much higher expression level in young flowers, which is about 6-fold higher than the 2nd highest level shown in leaves (Fig. 5). We further analyzed the expression level of GmMS6 in petals, sepals, and pistils, and found that GmMS6 is barely expressed in these floral parts (Fig. 5). Therefore, the high expression level detected in young flowers should be contributed by the gene expression in anthers, demonstrating that GmMS6 is an anther-specific gene, similar to the TDF1 proteins in Arabidopsis and rice (Zhu et al. 2008; Cai et al. 2015). On the contrary, GmMS6L was expressed at a super low level in all examined tissues, indicating it is likely in the process of pseudolization (Fig. 5). Different expression patterns of GmMS6 and GmMS6L illustrated that GmMS6 was the major functional TDF1 in soybean anther development, and explained why mutation at ms6 locus would lead to male sterility.
L76H does not alter the subcellular localization or transactivation activity of GmMS6
The subcellular localization of GmMS6, GmMS6L76H and GmMS6L were analyzed by transiently expressing their GFP fusion proteins driven by 35S promoter in N. benthamiana leaves. Free GFP was also expressed as a control, which showed signals all over the cells (Fig. 4a). Comparatively, the fluorescence from GmMS6-GFP and GmMS6L-GFP were restricted in nucleus, the general subcellular distribution of transcription factors (Fig. 4a). Additionally, GmMS6L76H-GFP was also present in nucleus, showing that L76H would not vary the subcellular localization of GmMS6.
TDF1 is a transcriptional activator in Arabidopsis and rice, therefore, GmMS6 that could compensate AtTDF1’s function should possess transcription activation activity as well. We then performed a transactivation activity test in yeast AH109 strain (Fig. 4b). Yeast clones expressing GAL4 DNA binding domain (BD) could only grow on SD medium lack of Trp (SD/-Trp) but not the selective medium (SD/-Trp-His) due to no transactivation activity in BD region. Similar phenomenon was observed for yeast strain expressing BD fused GmMS6DBD (BD-GmMS6DBD), which was a truncated form of GmMS6 only containing the N-terminal 191 residues (the DNA binding domain, 39–146 aa). In contrast, yeast clones expressing BD-GmMS6, BD-GmMS6L and BD-GmMS6L76H could grow well on both SD/-Trp and SD/-Trp-His medium, showing that these three proteins all possessed transactivation activity. Therefore, L76H doesn’t affect the transactivation activity of GmMS6. As L76 is a conserved residue in R2 motif of the DNA binding domain (Fig. 2d and 2f) and L76H has no effects on protein’s subcellular location or the transactivation activity, we suspected that L76H mutation likely disrupted GmMS6’s function by altering its DNA binding capacity.
The DYT1-TDF1-AMS-MYB80/MYB103/MS188-MS1 genetic pathway is present in soybean
In Arabidopsis, TDF1 is the critical component in a well-characterized genetic pathway regulating the tapetal development and pollen wall formation, which is composed of five transcription factors, including two basic helix-loop-helix (bHLH) factors DYSFUNCTIONAL TAPETUM 1 (DYT1) and ABORTED MICROSPORES (AMS), two MYB factors TDF1 and MYB80/MYB103/MS188, and one PHD-finger protein MALE STERILITY 1 (MS1) (Zhu et al. 2011; Lu et al. 2020). In this pathway, DYT1 directly activates the expression of TDF1, and TDF1 subsequently promotes the expression of AMS. Then, AMS is required for the expression of the gene encoding MYB80/MYB103/MS188, which is an activator critical for expressing MS1 (Fig. 7a). Depletion of any member in this cascade would lead to distorted tapetum and aborted pollens (Wilson et al. 2001; Sorensen et al. 2003; Zhang et al. 2006; Zhang et al. 2007). The same regulatory cascade known as UDT1-TDF1-TDR-OsMS188-PTC1 pathway was also identified in rice as well (Cai et al. 2015).
Soybean ms6 exhibited similar cytological abnormities to the null mutants attdf1 and ostdf1, such as vacuolated tapetum cells, undissolved callose, and crushed microspores (Fig. 1f-h; Ilarslan et al. 1999; Zhu et al. 2008; Cai et al. 2015). AtTDF1p-driven GmMS6/GmTDF1a was able to recover the fertility of attdf1 mutant like AtTDF1p-driven OsTDF1 did (Fig. 4d, 4k; Cai et al. 2015). These showed that TDF1’s function was conserved in Arabidopsis, rice, and soybean, and implied that DYT1-TDF1-AMS-MYB80/MYB103/MS188-MS1 genetic pathway was likely present in soybean as well. Therefore, we performed a homology search of these transcription factors in soybean genome, and found that the whole pathway was indeed present in soybean. Moreover, all the members in this cascade had multiple paralogs, particularly, four for DYT1 and two for the others (Fig. 7a). Their expression levels in WT and ms6 young flowers were assessed by qRT-PCR analysis. As a result, the expression level of GmDYT1(Glyma.13g250200) increased ~ 40% in ms6 (Fig. 7b), which is consistent to that TDF1 negatively feedback-regulates the expression of DYT1 in Arabidopsis (Cai et al. 2015). However, the expression levels of the other three GmDYT1s were not changed significantly in ms6, indicating the functions of DYT1s are diverged. Similarity, the expressions of two GmTDF1 were not affected in ms6 (Fig. 7b), suggesting that GmTDF1a is not involved in its own expression. On the other hand, the expressions of genes encoding the downstream transcription factors, like AMS, MYB80/MYB103/MS188, and MS1, were mostly downregulated in ms6 compared to WT (Fig. 7b), which is similar to the reports in attdf1 and ostdf1 (Zhu et al. 2008; Cai et al. 2015). One exception was Glyma.01g047400, encoding a MS1 homolog, which expressed similarly in WT and ms6 (Fig. 7b). This suggests that the expression of Glyma.01g047400 is not controlled through the DYT1-TDF1-AMS-MYB80/MYB103/MS188-MS1 genetic pathway and implies that Glyma.01g047400 is functionally diverged from its paralog Glyma.02g107600 or even might lost the function in anther development. Noticeably, Glyma.02g107600 is located at Gm02:10270908...10267934 in the region downstream of Satt157 (Gm02:9240028), which is rich of fertility controlling candidate genes, including ms3, msMOS, and female-partial sterile-1 (Fsp1) (Cervantes-Martinez et al. 2009).
Additionally, we examined the expressions of two enzyme coding genes that are regulated by this transcription factor cascade and related to microspores releasing and pollen wall formation, that is, A6 and MALE STERILE 2 (MS2) (Zhang et al. 2007; Zhu et al. 2008; Wang et al. 2018). A6 protein is proposed as a member of callase complex required to degrade the callose encompassing tetrads, as its protein sequence is highly similar to β-1,3-glucanase and its gene expression is correlated to callase synthesis temporally and spatially (Hird et al. 1993). The results showed that A6 gene only had one copy in soybean genome (Fig. 7a) and was barely expressed in ms6 (Fig. 7b), explaining why the callose surrounding tetrads was not dissolved in ms6 mutant. MS2 encodes a fatty acyl reductase catalyzing the palmitoyl acyl-carrier protein into a fatty alcohol, a precursor of sporopollenin (Aarts et al. 1997; Chen et al. 2011). These precursors are subsequently transported from tapetum cells into locule and deposited onto the microspore surface to produce the sexine layer (Wang et al. 2018). MS2 has two paralogs in soybean genome (Fig. 7a), both downregulated drastically in ms6 (Fig. 7b). These results showed that the DYT1-TDF1-AMS-MYB80/MYB103/MS188-MS1 genetic pathway is present and functions conservatively in soybean.