Genetic basis and adaptation trajectory of soybean from its temperate origin to tropics

Resequencing of soybean accessions from low latitudes

To investigate the genomic basis for the natural variation in soybean adaptation to low latitudes, we conducted whole-genome resequencing of a panel of 329 soybean accessions collected from 15 countries and covering all soybean subgroups in which 165 accessions are from in low-latitude regions (Supplementary Fig. 1a and Supplementary Data 1). Using the whole-genome single-nucleotide polymorphism (SNP) marker set, we performed phylogenetic analysis and principal component analysis (PCA) to quantify the population structure of these 329 soybean accessions. These analyses clearly classified the accessions into three groups: wild soybeans, landraces, and improved cultivars, and four main regions: China, Southeast Asia, South Asia, and South America (Supplementary Fig. 1b, c, Supplementary Fig. 2 and Supplementary Data 1). Consistent with previous findings^14,15, the decay of linkage disequilibrium (r²) with the physical distance between SNPs occurred in all three groups (Supplementary Fig. 1d).

Identification of the Tof16 locus

Using a linear mixed model for GWAS of the panel of 329 accessions, we identified one consistent significantly associated locus on chromosome 16 (P < 1.13 × 10⁻⁸) controlling flowering time under natural SD conditions in both 2018 and 2019 in Guangzhou. (Fig. 1a, b). This GWAS peak is consistent with the SD flowering QTL on chromosome 16 that we previously identified using two F₂ segregation populations, PI591429 × PI628930 and PI240664 × BR121¹⁶. These analyses across multiple biparental and natural populations indicate that variation at this locus (hereafter referred to as Time of Flowering 16, Tof16) is widespread and substantially contributes to the control of flowering time under SD conditions.

a–b GWAS scan for flowering time (R1 stage) using data from the 329-accession panel grown over the 2018 (a) and 2019 (b) field seasons in Guangzhou, China. c Phenotypes of tof16 mutation and wild-type Harosoy under SD (12 h light/12 h dark) conditions. Scale bar, 10 cm. d Flowering time. e Time to maturity. f Grain yield per plant of tof16^CR mutation and Harosoy. All data were given as mean ± s.e.m. (n = 10 plants), the value of each plant was represented by a dot. One-tailed Student’s t-test was used to generate the P values. g Diurnal variation in transcript levels of in wild-type Harosoy and tof6^CR mutant under SD. All data were given as mean ± s.e.m. (n = 5 plants). h Location of the AATATC-motif in E1 gene promoter. i ChIP-qPCR results demonstrated that Tof16 directly bound to the promoter of E1 gene. Values are mean ± s.e.m. (n = 3 biologically independent samples), the value of each replication was represented by a dot. Source data underlying Fig. 1d–g and  i are provided as a Source Data file.

Tof16 encodes LHY1a

We generated a large (n = 2418) inbred F₆ population from the PI591429 × PI628930 cross for fine-mapping of the Tof16 locus by recurrent selection for heterozygosity at Tof16 from the F₂ to F₅ generations. Analysis of this population located Tof16 within a 120-kilobase (kb) region harboring 12 genes based on the Williams 82 reference genome¹⁷ (Supplementary Fig. 3a, b and Supplementary Table 1). We cloned and sequenced all 12 of the predicted genes in the two parents; of these, the sequence of the circadian clock gene LHY1a (Glyma.16G017400) differed between the two parents (PI591429 and PI628930) (Supplementary Fig. 3b). The late flowering parent PI628930 harbored two SNPs predicted to cause a gain of a stop codon, resulting in premature termination of translation after 159 amino acids in the 750-amino-acid LHY1a protein (Supplementary Figs. 3b and 4). We also sequenced the coding region of LHY1a in PI240664 and BR121, finding that in PI240664, this gene harbored one SNP predicted to convert a serine into a cysteine, which is a conserved site in leguminous plants (Supplementary Figs. 4 and 5). The presence of different mutations in these two lines suggested that the LHY1a gene was a strong candidate for the Tof16 locus.

To validate whether LHY1a is the causative gene of the Tof16 locus, we generated loss-of-function mutations of LHY1a (named tof16^CR) in the Harosoy background using CRISPR/Cas9-mediated gene editing and evaluated the phenotypes of the mutants vs. wild-type Harosoy (Supplementary Fig. 6). The tof16^CR plants showed significantly delayed flowering time and maturity, altered yield-related traits, and improved overall grain yield relative to Harosoy (Fig. 1c–f, Supplementary Fig. 7a–d). These results confirm the notion that LHY1a is the causative gene of the Tof16 locus and that two natural mutations arose in PI628930 and PI240664 (hereafter referred to as tof16–1 and tof16–2, respectively).

To examine the specific effects of the Tof16 locus, we also compared the phenotypes of two F₇ near-isogenic lines (NILs) carrying either the functional Tof16 allele (NIL-Tof16) or the non-functional tof16–1 alleles (NIL-tof16-1 or NIL-tof16-2) in different genetic backgrounds. NIL-tof16–1 and NIL-tof16–2 showed significantly delayed flowering time and maturity (Supplementary Fig. 8a–c, e–g), along with increased plant height, node number, pod number, branch number (Supplementary Fig. 7e–l), and grain yield compared to the functional NILs (Supplementary Fig. 8d, h), confirming the notion that the tof16–1 and tof16–2 alleles delay flowering and greatly enhance grain yields in soybean under SD conditions. Therefore, like J, Tof16 also functions as a flowering enhancer; the loss of function of both genes might have contributed to the adaptation of soybean to the tropics.

Tof16 is genetically dependent on the legume-specific flowering repressor E1

E1 plays a central role in photoperiodic flowering by repressing the expression of two key FT homologs, FT2a and FT5a^14,18. There are three major natural alleles of E1 in soybean: E1, e1^as, and e1^fs; e1^as is a weak mutant allele and e1^fs is a null functional allele¹⁸. To examine the genetic interaction of Tof16 and E1, we developed a NIL set for E1/Tof16, E1/tof16^CR, e1^as/Tof16, e1^as/tof16^CR, e1^fs/Tof16, and e1^fs/tof16^CR in the Harosoy background and performed phenotypic analysis. The tof16^CR allele delayed flowering and maturity in both the E1 and e1^as genetic backgrounds, but the effect was weaker in the e1^as background. By contrast, the effect of Tof16 on flowering was completely eliminated in the e1^fs null functional background (Supplementary Fig. 9a–c), implying that the full effect of Tof16 on flowering depends on E1.

We evaluated the effect of Tof16 on the transcriptional regulation of E1, FT2a, and FT5a under SD (12 h light/12 h dark) conditions using tof16^CR and Harosoy or independent NIL pairs for each locus. As expected, E1 was expressed at higher levels in tof16^CR (Fig. 1g), and FT2a and FT5a were expressed at lower levels in this line compared to Harosoy (Supplementary Fig. 10a, b). A similar result was obtained for the NILs. These results indicate that functional alleles of Tof16 repress E1 expression and increase FT2a and FT5a expression relative to the mutant alleles (Supplementary Figs. 10c, d and 11).

To further explore the molecular nature of the relationship between Tof16 and E1, we determined whether Tof16 directly binds to the promoter of E1 in vivo. We generated transformants overexpressing Tof16-6HA in the Williams 82 background and subjected them to chromatin immunoprecipitation (ChIP)-qPCR assays (Supplementary Fig. 12). Tof16 directly associated with the E1 promoter regions that contained AATATC motif (a part of the EE motif, Fig. 1h, i). These results consist with our previous finding that Tof16 protein could bind to AATATC motif in the E1 promoter in vitro¹⁴. Taken together, these results indicate that Tof16 enhances early flowering and maturity by direct associating with the E1 promoter to suppress E1 expression, thus releasing FT2a and FT5a transcription.

Four LHY homologs redundantly control soybean flowering and grain yield

The soybean genome contains four LHY/CCA1 homologs (LHY1a, LHY1b, LHY2a, and LHY2b)^14,19 and its amino acid sequences were high homology (Supplementary Fig. 5). Since we demonstrated that Tof16/LHY1a controls flowering time, maturity, and grain yield in soybean, we asked whether the other LHY family members also control these soybean traits and whether these homologs could potentially be used for agricultural applications. We crossed the lhy1a/1b/2a/2b quadruple mutants with wild-type Harosoy and obtained all 15 homozygous mutational combinations of LHY (Fig. 2a). We examined the phenotypic differences of the mutants grown in fields in Guangzhou under natural SD conditions. Among the single mutants, lhy1a and lhy1b, but not lhy2a or lhy2b, showed significantly delayed flowering time and maturity and improved overall grain yield relative to Harosoy (Fig. 2a–e, Supplementary Fig. 13a–d). All of the multiple mutants except the lhy2c/2d double mutants showed significantly delayed flowering time and increased grain yield in a quantitative manner compared to wild-type Harosoy (Fig. 2a–e, Supplementary Fig. 13a–d). Strikingly, the lhy1a/1b/2b triple mutants exhibited the best architecture and highest grain yield but shorter flowering time and earlier maturity compared to the lhy1a/1b/2a/2b quadruple mutants under natural SD conditions (Fig. 2a–e, Supplementary Fig. 13a–d). Furthermore, all mutants possessing the lhy1a (tof16) mutation had higher grain yields than those without this mutation, indicating that lhy1a plays a crucial role in controlling flowering time, maturity, and grain yield under SD conditions (Fig. 2a–e, Supplementary Fig. 13a–d).

a Phenotypes of lhy mutants. Scale bar, 10 cm. b Flowering time. c Time to maturity. d Total pods per plant height. e Grain yield per plant. All data were given as mean ± s.e.m. (n = 10 plants), the value of each plant was represented by a dot. The presence of the same lowercase letter above the histogram bars in (b–e), denotes nonsignificant differences across the two panels (P > 0.05). One-way ANOVA was used to generate the P values. Source data underlying Fig. 2b–e are provided as a Source Data file.

To gain further insight into how LHY homologs control flowering time and grain yields under long-day (LD) conditions, we also evaluated the mutants in the field in Changchun under natural LD conditions. All mutants except lhy2a, lhy2b, and lhy2a/2b showed delayed flowering time and maturity compared to the wild type (Supplementary Fig. 14a–i). Due to the larger effect of LHY1a on flowering time, multiple mutants carrying the lhy1a mutation failed to mature, and it was difficult to harvest the seeds naturally until the end of the growing season (Supplementary Fig. 14a–i). However, unlike their performance under SD conditions, the lhy1a and lhy1b2a2b mutants exhibited the best architecture and improved grain productivity under natural LD conditions. These results suggest that LHY homologs have redundant but divergent functions in controlling flowering time, maturity, and grain yield in soybean under both SD and LD conditions. We further examined the expression of four LHY genes under LD and SD conditions, and found that the expression of four LHY genes were no significant difference under LD and SD conditions (Supplementary Fig. 15). We also tested the expression levels of E1 in lhy multiple mutants, the result showed that the amount of E1 correlated with flowering and maturity under LD and SD conditions, but not yield (Supplementary Figs. 16–19). Therefore, manipulating the combinations of these alleles and modulating the genetic complexity of the LHY homologs could help create the appropriate genotypes to maximize the adaptation and yield potential of soybean at different latitudes.

Tof16 and J additively control flowering time and grain yield

The adaptation of soybean to low latitude or tropical regions largely depends on the natural loss of function of the flowering enhancer J¹². Therefore, it is critical to explore the genetic relationship between Tof16 and J. We tested the reciprocal transcriptional regulation between Tof16 and J using mutants or NIL sets of tof16 or j. No mutual transcriptional regulation was observed between Tof16 and J (Supplementary Fig. 20). To further explore the genetic interaction of Tof16 and J, we developed two NIL sets for the four different homozygous allelic combinations at two loci from a cross between tof16^CR/E1 and NIL-j/E1 or a cross between tof16^CR/e1^as and NIL-j/e1^as in the Harosoy background and subjected them to phenotypic evaluation. The presence of a recessive allele of either Tof16 or J delayed flowering and maturity and enhanced grain yield in both the E1 and e1^as genetic backgrounds. However, the double recessive mutant tof16 j showed significantly later flowering, maturity, and higher grain yield than either the j or tof16 single mutant in both the E1 and e1^as genetic backgrounds (Fig. 3a–d, Supplementary Fig. 21a–l), suggesting that Tof16 and J additively control flowering and grain yield in a genetically independent manner.

a Phenotypes of NILs possessing different allelic combinations at Tof16 and J in E1 background under SD (12 h light/12 h dark) conditions. Scale bar, 10 cm. b Flowering time. c Time to maturity. d Grain weight per plant. All data were given as mean ± s.e.m. (n = 10 plants), the value of each plant was represented by a dot. The presence of the same lowercase letter above the histogram bars in (b–d) denoted nonsignificant differences across the two panels (P > 0.05). One-way ANOVA was used to generate the P values. e–g Diurnal variation in transcript levels of E1 (e), FT2a (f), FT5a (g) in possessing different allelic combinations at Tof16 and J in E1 background under SD conditions. All data were given as mean ± s.e.m. (n = 5 plants). h Combining of various CRISPR/Cas9 generated mutants of LHY allows improve soybean adaption to tropic regions and yield. The value represents the average grain weight per plant in fields of Guangzhou under natural SD conditions. i Combining of natural or gene-edited of Tof16 and J alleles in back ground E1 or e1^as enhance soybean yield in tropic regions. The value represents the average grain weight per plant in fields of Guangzhou under natural SD conditions. Source data underlying Fig. 3b–g are provided as a Source Data file.

Consistent with this genetic effect, E1 transcript levels were highest in NIL-tof16^CR/j, followed by NIL-Tof16/j or NIL-tof16^CR/J, and were lowest in NIL-Tof16/J. As expected, FT2a and FT5a showed the opposite expression pattern (Fig. 3e–g). These results, together with our previous findings¹², suggest that the positive regulators of flowering Tof16 and J both promote flowering, which depends on the function of E1. However, Tof16 and J independently but additively regulate flowering time, maturity, and grain yield in soybean.

Based on these findings, we propose two possible methods for the quantitative improvement of flowering time and grain yield in soybean in tropical regions (Fig. 3h, i): (1) Due to the genetically redundant and divergent effects of LHY homologs, genotypes with a quantitative series of flowering time, maturity, and yield traits could be created by combining various CRISPR/Cas9-generated mutants (Fig. 3h). (2) Due to the genetic effects of Tof16, J, and E1, their various natural or artificial alleles could be combined, allowing another set of genotypes conditioning quantitative traits to be produced (Fig. 3i). How these genotypes could be selected or utilized remains to be explored and depends on the photoperiod or latitudinal environment.

Stepwise selection of Tof16 and J during soybean adaptation

Our findings indicate that the positive regulators of flowering Tof16 and J play essential roles in the adaptation of soybean to SDs and yield development. To gain insight into the evolutionary trajectory of soybean adaptation from its temperate origin to the tropics, we examined the genomic variations in the Tof16 and J coding sequences in 1624 resequenced soybean accessions, including 1295 previously described accessions^14,20 and the 329-accession panel used in the current study. We identified 15 unique high-confidence haplotypes in Tof16. In addition to haplotypes H10 (tof16-1) and H11 (tof16-2), we identified two novel loss-of-function haplotypes: H1 (named tof16-3) and H8 (named tof16-4) (Supplementary Fig. 22a, Supplemental Data 2). Analysis of haplotype origin indicated that the tof16-2 allele (H11, Tof16-SNP^A1276T) first occurred in wild soybeans originating in the Huanghui region and were subsequently domesticated into landraces in areas where soybean domestication occurred (Fig. 4a, Supplementary Fig. 22b). Following domestication and during the adaptation of soybean to low latitudes, H11 (tof16-2) was intensively selected in the accessions that adapted to low-latitude regions, suggesting that this haplotype plays critical roles in soybean adaptation to the tropics (Fig. 4a). Interestingly, H10 (tof16-1) and H8 (tof16-4) originated from H11 (tof16-2), suggesting that loss-of-function alleles of tof16-2, tof16-1, and tof16-4 were under stepwise selection during adaptation to low latitudes. However, H8 (tof16-3) originated from H1 and only occurred in India and Nepal, while tof16-1 and tof16-4 mainly occurred in Brazil, indicating that unique loss-of-function alleles of tof16 were selected in different regions to enhance the adaptation of soybean (Fig. 4a, Supplementary Fig. 22b).

a Loss-of-function alleles of Tof16 frequency is highly correlated with low latitude regions. Data from 1624 diversity panels. b–c Flowering time(R1) variations in 329 accessions possess Tof16 and tof16 alleles in Guangzhou 2018 (b), and Guangzhou 2019 (c). Proportions of j alleles in upper pie chart. Dark gray represents J alleles and light gray represents j alleles. d Loss-of-function alleles of J frequency is highly correlated with low latitude regions. Data from 1624 diversity panels. j⁹²⁰ represents j-11 allele including J-SNP^G920T. j^920bm represents SNP920-based mutations including j-1, j-3, j-6 and j-10 mutational alleles. j-others represents j mutations other than j^920bm. e–f Flowering time(R1) variations in 329 accessions possess J and j alleles in Guangzhou 2018 (e), and Guangzhou 2019 (f). Proportions of tof16 alleles in upper pie chart. Dark gray represents Tof16 alleles and light gray represents tof16 alleles. g Loss-of-function Tof16 and J improve soybean adaptation to low latitudes. Data from 1624 diversity panels. h–i Flowering time variations of eight allelic combinations of Tof16 and J in Guangzhou 2018 (h), and Guangzhou 2019 (i). The presence of the same lowercase letter above the histogram bars in (b, c, e, f, h, i) denotes nonsignificant differences across the two panels (P > 0.05). One-way ANOVA was used to generate the P values. All data were given as mean ± s.e.m. The value of each plant was represented by a dot. Source data underlying Fig. 4b, c, e, f, h, and i are provided as a Source Data file.

To further validate the functional significance of tof16-1, tof16-2, tof16-3, and tof16-4, we performed a transient transfection assay in Arabidopsis thaliana protoplasts. tof16-2 partial impaired the ability of Tof16 to repress the expression of E1, but tof16-1, tof16-4, and tof16-3 completely impaired this ability (Supplementary Fig. 23). These results imply that partial loss of function of tof16-2 (standing variations from soybean in the central area of origin) was first selected during soybean adaptation to the tropics but was not sufficient for full adaptation. In this genetic background, null alleles of tof16-1 and tof16-4 occurred and were further selected for better adaptation to the tropics. Population genetic association analysis of flowering time in the 329-accession panel showed that tof16-2 flowered the latest, followed by tof16–3, tof16-1,4, and Tof16 in both 2018 and 2019 (Fig. 4b, c). The unexpected finding that the weak functional allele tof16-2 flowered later than the loss-of-function alleles tof16-1, tof16-4, and tof16-3 could be explained by the interference of the loss of function of j (Fig. 4b, c). Taken together, these data suggest that selection at tof16-1or tof16-4 and tof16-2 arose in a stepwise manner. Tof16 loss-of-function alleles independently originated and were selected in two important soybean planting areas (Brazil and India) in the tropics.

We screened for natural variation of the J coding sequence in the same 1624-accession panel. In total, 28 haplotypes were identified, including seven distinct loss-of-function alleles and two weak loss-of-function alleles. These haplotypes included the previously reported alleles j-1, j-2, j-3, j-4, j-6, j-8 (including j-8-1 and j-8-2), and j-9 (e6)^12,13 and the newly discovered alleles haplotype H18 (named j-10) and haplotype H21 (named j-11) (Supplementary Fig. 24a, Supplementary Data 2). Like tof16-2, haplotype origin network analysis indicated that haplotype H21 (j-11, a SNP^G920T resulting in an amino acid substitution) also first occurred in wild soybean originating from Huanghuai and was later was domesticated into landraces, but this allele was substantially selected in accessions grown in low-latitude regions (Fig. 4d, Supplementary Fig. 24b). Furthermore, loss-of function alleles j-1, j-3, j-6, and j-10 all occurred in the tropics and originated from haplotype H21 (j-11), indicating that, similar to tof16, stepwise selection of the weak allele of j-11 and other loss-of-function alleles of j occurred during the adaptation of soybean to the tropics (Fig. 4d, Supplementary Fig. 24b).

Transient transfection assays also demonstrated that j-11 partially impaired the ability of J to repress the expression of E1 (Supplementary Fig. 25), indicating that H21 (j-11, SNP j⁹²⁰) is a weak loss-of-function allele that might contribute to adaptation to low latitudes. We further classified all haplotypes of J into four groups, J, j⁹²⁰, j^920bm (SNP920-based mutations, including the j-1, j-3, j-6, and j-10 mutant alleles), and j^other (j mutations other than j^920bm) and evaluated their associations with flowering in the 329-accession panel. j^920bm flowered the latest, followed by j^other, j⁹²⁰, and J in both 2018 and 2019 (Fig. 4e, f). Taken together, our data suggest that like the selection at Tof16, the selection of j^920bm and j⁹²⁰ also arose in a stepwise manner.

Selection of natural mutations of Tof16 and J allowed soybean to move into the tropics

We firstly investigated that the distribution of the various natural alleles of E1, Tof16, and J in tropical soybean accessions. We found that all of the tropical accessions harbor the dominant E1 allele (Supplementary Data 3), and the frequency of the tof16-2 allele (26.75%) and j-11 allele (47.42%) were highly variable in tropical accessions (Supplementary Tables 2 and 3). To further explore the contributions of mutations of Tof16 and J to the adaptation of soybean to low latitudes, we grouped eight allelic combinations (Tof16/J, Tof16/j⁹²⁰, Tof16/j^others, tof16-2/J, tof16-2/j⁹²⁰, tof16-2/j^920bm, tof16-1,4/j⁹²⁰, tof16-3/j⁹²⁰) and examined the geographic distributions of the eight Tof16/J allelic combinations in the 1624 accessions covering all latitudes. Accessions carrying loss-function-of Tof16 alleles were enriched in Brazil, but accessions carrying loss-function-of J alleles were enriched in Southern China, Southeast Asia, and Brazil (Fig. 4g), suggesting that the selection of mutations of tof16 and j might have occurred independently. Interestingly, by extracting the 165 accessions from low latitude, we found that more than 80% accessions harbor the loss of functions of Tof16 and J suggest the two genes are the major genetic forces to drive soybean adaptation from temperate into tropics (Supplementary Data 4, Supplementary Table 4). By contrast, most accessions adapted to high-latitude regions such as Northern China, the United States, and Canada, where early flowering is required, carried the full functional alleles of Tof16 and J (Fig. 4 a, d, g), suggesting that the mutations of these two genes helped soybean move from its temperate origins into the tropics.

Finally, to further explore the functional significance of tof16 and j, we examined their association with flowering time in the 329-accession panel in two environments in Guangzhou (2018 and 2019). Accessions carrying two recessive alleles flowered significantly later than accessions carrying single recessive alleles tof16 or j, and accessions carrying both functional alleles (Tof16/J) flowered earlier than the other three genotypic groups in all environments (Fig. 4h, i). These results further confirm the genetic additive effects of Tof16 and J revealed from their interactions in NILs (Fig. 3b–d). Therefore, the selection of natural mutations of both genes has played substantial roles in expanding soybean cultivation from its temperate origins into the tropics and has facilitated the improvement of soybean adaptation and yield in the tropics.