Geographic differentiation and phylogeographic relationships among world soybean populations

Xueqin Liu,Jianbo He, Yufeng Wang, Guangnan Xing, Yan Li, Shouping Yang,Tuanjie Zhao, Junyi Gai*

Soybean Research Institute,Nanjing Agricultural University,Nanjing 210095,Jiangsu,China

MARA National Center for Soybean Improvement, Nanjing Agricultural University,Nanjing 210095,Jiangsu,China

MARA Key Laboratory of Biology and Genetic Improvement of Soybean(General),Nanjing Agricultural University, Nanjing 210095,Jiangsu,China

State Key Laboratory for Crop Genetics and Germplasm Enhancement,Nanjing Agricultural University,Nanjing 210095,Jiangsu,China Jiangsu Collaborative Innovation Center for Modern Crop Production,Nanjing Agricultural University,Nanjing 210095,Jiangsu,China

A B S T R A C T

1.Introduction

Originating in ancient China, soybean (Glycine max (L.) Merr.)is now cultivated worldwide. It has been recognized as a miracle crop in North and South America over the past six decades for its nutritional value(the seed contains about 40%protein and 20% oil) and high economic profitability [1].Worldwide soybean production is increasing year by year and amounted to 337 Mt. in 2017/2018, reaching a historical peak(Table S1).

In China, the major soybean-producing areas are the northern single-cropping area(including northeast and northwest China),the Huang-Huai(Yellow River-Huai River)double cropping area,and the southern multiple-cropping area[2-4].For the center of origin of soybean, the Huang-Huai [4-6] and Changjiang valleys and the region to their south (southern China) [7-9] have been suggested. Over its long history,farmers in ancient China developed many soybean landraces for their own production in each geographic region.Scientific breeding of soybean cultivars began in the early 20th century and about 2500 cultivars have been released.These landraces and released cultivars constitute the soybean germplasm collection.

Singh and Hymowitz [10] roughly mentioned on the basis of historical records and ancient legends that soybean was disseminated from China to Asia, Europe, and the USA. The literature [11-13] records its dissemination to Canada, South America, and Africa. The putative paths of soybean dissemination from China to other regions of the world are shown in Fig. 1. Soybean reached Asian countries around the first century A.D. and other continents in the last three centuries[10,14].

With the dissemination and adaptation of soybean to other continents, soybean landraces were developed by local farmers before scientific breeding began, and commercial soybean cultivars were developed by breeders after that point.Given that introductions to different ecological-regions may have come from different and even multiple sources, genetic backgrounds may differ among geographic- and ecologicalregions. In particular, landraces and improved varieties adapted to individual regions were developed as a consequence of adaptation to local environments and diverse breeding objectives, resulting in differentiation among regional populations and differing agronomic performance[15-17].

As a short-day plant, soybean is greatly affected by daylength and temperature variation among growing regions,and depends particularly on local latitude and altitude as well as the sowing seasons in multiple cropping systems. During dissemination and adaptation, sharp differentiation is expressed in traits influenced by daylength and temperature,such as flowering date, maturity date, and main stem node number[18-20].Apart from these geographic-ecological traits,the genetic diversity and differentiation of soybeans can be detected as genotypic variation.Abe et al.[21]used 20 SSRs to genotype 131 soybean accessions introduced from 14 Asian countries and found two germplasm pools, namely Japanese and parts of Korean peninsula populations clustered in one group, and most populations in Southeast and South and Central Asia with Chinese and parts of Korean peninsula populations clustered in the other group. Ude et al. [22] used five AFLPs to investigate the genetic diversity and relationships among North American soybean ancestors (NASA),North American soybean cultivars (NASC), and modern Chinese and Japanese cultivars, and found that Japanese and Chinese cultivars were distinct from NASA and NASC.However, knowledge of worldwide soybean variation and differentiation is limited. In particular, little is known about the phylogeographic relationships among soybean populations worldwide. This is because of difficulties in collecting worldwide soybean germplasm,given that governments have diverse germplasm exchange policies and it is difficult to test multiple maturity groups (MGs) in the same natural environment(location and season).

Fig.1- The paths of soybean dissemination from the center of origin to other world regions.Arrows indicate directions of soybean dissemination.Thirteen regional populations are shown in different colors.

Efforts have been made [10,14] to collect a representative sample of soybean accessions worldwide as with the aim of investigating phylogeographic relationships among soybeans in diverse geographic regions. The present study aimed to characterize the geographic-ecological differentiation among soybeans worldwide, phylogeographic relationships among geographic-regional populations,and ultimately the historical dissemination and evolutionary development of cultivated soybean. The phylogeographic relationships among geographic-regional populations involved two axes, namely time (history) and area (geographic region). In the present study, the area axis was real, but no time axis could be constructed because the accessions collected for this study were modern rather than historical (already extinct) accessions. Accordingly, the historical phylogeographic relationship among geographic-regional populations could be inferred only from present-day geographic-regional populations (13 geographic-populations composed of 371 accessions). For this purpose, two indicators were used to describe geographic patterns.One was geographically and ecologically characterized traits (geographic-ecological traits) sensitive to geographically associated daylength and temperature(flowering date, maturity date, main stem node number,geographic-ecological-trait-based cluster). The other was gene/allele-derived or marker/haplotype-derived traits,which are genetic traits passed in regional populations:allele richness, allele frequency dispersion, specific-present allele(SPA),specific-deficient allele(SDA),and allele constitution,as well as gene/allele-based or marker/haplotype-based cluster.

2. Materials and methods

2.1. Plant materials

Based on available information about soybean dissemination,world soybean production, and collected soybean accessions,the 13 geographic regions were grouped as shown in Fig. 1. A total of 371 accessions from the 13 geographic regions were sampled from the germplasm collection of the National Center for Soybean Improvement(NCSI),Nanjing Agricultural University, Nanjing, China (Table S2). The 13 regional populations were labeled with the assumed stages of their dissemination for convenience of description. Here, O represents the centers of origin in the Huang-Huai Valleys (OHCHN)and Changjiang Valley and south of Changjiang Valley(O-SCHN) in China because they were both nominated in the literature [4-7,9,23] as centers of origin. A, B, C, and D represent four secondary centers extending outward from the center of origin.“A”represents northeast China(A-NCHN,including a few accessions from northwest China and in fact introduced mainly from northeast China in recent decades.These two parts comprise the northern single-cropping soybean area, which is conventionally abbreviated as northeast China).“B”represents the Korean peninsula(B-KORP)and Japanese islands (B-JPAN). “C” represents Southeast Asia (CSEAS).“D”represents northern North America(D-NNAM)and southern North America (D-SNAM). A1, A2, C1, C2, and D1 represent tertiary centers derived from A, B, C, and D,including the far east of Russia (A1-RUFE), southern Sweden(A2-SSWE,which ought to be a sample collected from Europe,but only accessions from southern Sweden could be used because of few European accessions in the germplasm collection),South Asia(C1-SASI),Africa(C2-AFRI),and Central and South America (D1-CSAM). No accessions from Australia were included in the study, owing to the short history of soybean cultivation in Australia and the absence of Australian soybeans from the NCSI collection. The 371 accessions originated in China, the USA, Brazil, and 24 other countries,spanning the belt between roughly 58°N and 34°S(Table S2).

2.2. Field experiments

From 2011 to 2012, the 371 accessions were tested in one-row plots length 1 m, an interval of 0.4 m, and two replications at Jiangpu Experimental Station (32°07′N, 118°62′E) of Nanjing Agricultural University. At this latitude, all accessions could mature or nearly mature under spring sowing conditions,contributing to the comparability of their geographic-ecologicaltraits. Sowing dates were April 28 in 2011 and April 25 in 2012.This was a full-season test (rather than local double cropping with sowing after the harvest of wheat). Geographic-ecologicaltraits were observed, including phenology stages of emergence(VE), beginning bloom (R1) and full maturity (R8) following Fehr and Caviness[24].Flowering and maturity dates were calculated as the periods from sowing to R1 and to R8,respectively.If some accessions were unable to mature normally, the maturity date was estimated by adding the maturity date of MG VIII to the difference in flowering date between MG VIII and immature MG IX and X.Main stem node number was recorded at maturity with the cotyledonary node numbered zero,with counting to the top of the main stem,in 2012 only.

2.3. Statistical and clustering analysis of geographicecological-traits

PROC MEANS in SAS 9.1(SAS Institute Inc.,Cary,NC,USA)was used to calculate the descriptive statistics (mean, minimum,and maximum) of flowering date, maturity date, and main stem node number. PROC GLM was used to fit an analysis of variance (ANOVA) from which the genetic coefficient of variation (GCV) was obtained from the equation. GCV(%) =σg/μ × 100%, where σgis the genotypic standard deviation calculated according to the expected mean squares in ANOVA and μ is the population mean. PROC CLUSTER was used to cluster the 371 accessions and 13 regional populations to obtain phenotypic data using the method of average linkage.

2.4. Genotyping

Restriction-site-associated DNA sequencing (RAD-seq) was performed by BGI Tech,Shenzhen,China.Genomic DNA of all 371 accessions was extracted from young leaves using the CTAB method [25]. The sequence depth (120.17 Gb of sequence) of the 371 accessions was approximately 4.08×with 4.64% coverage. Sequence was generated with an Illumina HiSeq 2000 instrument using multiplexed shotgun genotyping [26]. All sequences were aligned against Williams 82 [27] using SOAP 2 [28], using sequence similarity, pairedend relationships, and sequence quality. Using RealSFS software [29], 98,482 single nucleotide polymorphism (SNPs)were identified for the 371 accessions with a proportion of missing and heterozygous alleles ≤30% and minor-allele frequencies ≥1%. FastPHASE software [30] was used for SNP genotyping imputation after heterozygous loci were transformed into missing alleles. After that, SNPs within a linkage disequilibrium (LD) block (D' ＞0.7) were organized into an SNPLDB (SNP linkage disequilibrium block) marker with its haplotypes using Haploview software[31,32].A total of 20,701 SNPLDBs and their 55,404 haplotypes were identified.Liu et al.[33]used the same genotype dataset in a previous study.

2.5.Analysis of the genetic marker-haplotypes(gene-alleles)of populations

Because marker-haplotype data were used for the genetic analysis of the populations a marker with its haplotype was analogous to a gene or QTL (quantitative trait locus) with its alleles. In the following description, gene-allele and markerhaplotype, which are analogous, are used alternately for simplicity. Genetic or gene diversity was analyzed in its broad sense, including allele richness (the first-rank statistic of gene diversity)and allele frequency dispersion(also known as gene diversity but the second-rank statistic of gene diversity in a narrow sense). These statistics were calculated with PowerMarker 3.25 [34]. Allele richness (A) refers to the total number of alleles at all marker loci in a population and is calculated as A = ∑Ai, where Airepresents the number of alleles at the ith locus in the population. Allele frequency dispersion (D) means the variation of allele frequency in the population and is calculated as＝1-, whererepresents the mean allele frequency dispersion at the lth locus, ~plurepresents the frequency of the uth allele at the lth locus,and k is the total number of alleles at the lth locus.

Wen et al.[35]defined the SPA of a population as an allele present in that population and in no other population and SDA as an allele present in all other populations but absent from that population. These two indicators were used to evaluate the specificity of the 13 geographic populations.

2.6. Genetic differentiation and phylogeographic analysis

Arlequin 3.11 [36] was used to calculate differentiation coefficient FSTstatistics as where niis the allele number of the ith population, P is the total number of populations,and

where SSD(AP)is the sum of the squares of populations,SSD(WPi) is the sum of the squares of the ith population, N is the total number of genotypes, and σT2is the sum of expected mean squares.

Genetic distance (dij) was calculated for all pairs of accessions: dij= 1 - sij, where sijis the coefficient of genetic similarity,calculated as

where nijkis the common allele number of the ith and jth individual at the kth locus,and m is the total number of loci.

Clustering was performed using the genetic distance matrix of dijfor both the 371 accessions and the 13 regional populations with the neighbor-joining procedure of PowerMarker 3.25 and displayed by MEGA 4 [37]. On this basis, the results of geographic-ecological traits were combined to infer phylogeographic relationships among the geographic populations.

3. Results

3.1.Phylogeographic relationships among geographic-regional populations in terms of geographic-ecological traits

As shown in Table 1,world soybean accessions with their MGs from 000 to X covering all kinds of MGs had flowering dates ranging from 30 to 123 days, maturity dates from 75 to 200 days, and main stem node number from 4 to 35 nodes under uniform full-season environmental conditions in Nanjing, China. This result revealed an extremely wide range of variation in these traits. Mean flowering dates, maturity dates, and main stem node numbers in centers of origin (O)were respectively 51-62 days (ranging from 37 to 104 days),128-141 days (ranging from 95 to 178 days) and 15-16 nodes(ranging from 8 to 28 nodes).Those in northern centers(A,A1,and A2) were shortened; those in the Korean peninsula and Japanese centers(B)were similar to those in centers of origin;those in southern centers(C,C1,C2)were markedly increased;those in new American centers (D) showed both lower and higher values. Table 1 shows the differences in means and ranges of these three traits among the 13 regional populations and even within each population,indicating the adaptation of soybeans to different geographic environments (daylength and temperature conditions) with the dissemination of soybeans to different areas around the world.

The 13 geographic ecological-populations were phenotypically clustered based on the three daylength-and temperaturesensitive traits(Fig.2-A).Generally,A-NCHN,A1-RUFE,A2-SSWE,and D-NNAM were in a daylength and temperature mostinsensitive group; O-HCHN, O-SCHN, B-KORP, and B-JPAN were in a less sensitive group; and C-SEAS, C1-SASI and C2-AFRI, DSNAM,and D1-CSAM were in a sensitive group,with“C”groups more sensitive than “D” groups. For individual traits, the clustering results for maturity date and main stem node numberwere very similar to those for the three traits, whereas the clustering results for flowering date were somewhat different from those of the other two,with D-SNAM and D1-CSAM joining the O-SCHN group(Fig.S1,Fig.2-A).

Table 1-Phenotypic diversity and differentiation among regional populations of three daylength- and temperaturesensitive traits of soybeans grown in Nanjing(32°07′N,118°62′E),China.

3.2.Phylogeographic relationships among geographic-regional populations in terms of genetic marker-haplotype (gene-allele)traits

The analysis of molecular variance showed highly significant genetic differences among and within regional populations(Table S3). In the centers of origin (primary centers), O-HCHN and O-SCHN harbored 42,519 and 43,349 alleles accounting for 76.7% and 78.2% of the alleles (55,404) in all the populations(Table 2), which meaning that 23.3% and 21.8% of the alleles arose in geographic regions other than the center-of-origin regions. In the secondary centers, A-NCHN, B-KORP, and BJPAN harbored 41,439 (74.8%), 41,200 (74.4%), and 43,942(79.3%) alleles, respectively, showing proportions similar to those in the center of origin, although their areas were much smaller. However, in other secondary centers, C-SEAS and DNNAM harbored 44,141 (79.7%) and 44,986 (81.2%) alleles,respectively,showing a marked increases.Allele richness in the newly formed tertiary centers A2-SSWE, C2-AFRI, and D1-CSAM was lower than that in the major centers O,A,B,C,and D.

Fig.2-Unrooted trees showing the relationships among 13 regional populations.(A)Phenotypic clustering by three daylengthand temperature-sensitive traits using the method of average linkage.(B) Genotypic clustering using neighbor joining.

A total of 247 SPAs and 4112 SDAs were detected. In the older centers O-HCHN, O-SCHN, B-KORP, and C-SEAS, SPAs were only 8, 11, 2, and 4, respectively, whereas in the newly formed center D-NNAM and in the tertiary centers A1-RUFE,A2-SSWE, and C1-SASI, SPAs were as many as 28, 41, 41, and 67, respectively (Table 2, Fig. 3). To be specific, 24 SPAs were found in the older center B-JPAN. No SPAs were found in the younger tertiary centers C2-AFRI and D1-CSAM. SDAs followed similar trends.

Gene diversity, D in a narrow sense (allele frequency dispersion) indicates the evenness of allele frequencies. In the process of evolution, new mutants led to a low D-value,which increased to equilibrium after a long period of adaptation. In the present study, the center of origin “O” and secondary centers “B” and “C” showed relatively stable Dvalues of 0.204-0.235. However, new secondary centers DNNAM and D-SNAM showed relatively unbalanced D-values of 0.157-0.175, as did the secondary center A-NCHN and the tertiary centers A2-SSWE and D1-CSAM. This finding indicates the emergence or introduction of more new alleles in these secondary and tertiary centers(Table 2).

The above measurements of genetic diversity showed differences of allele frequency dispersion among geographic populations, while differentiation coefficient (Fst) and genetic similarity coefficient (sij) measured the degree of differentiation between two populations (Table 3). The overall mean Fstbetween populations was 0.188, with individual Fstvalues ranging from 0.036 to 0.422. Very small Fstwas observed between O-HCHN and O-SCHN(0.054),indicating high consistency between the two center-of-origin populations. Similarities were observed between the two “B”s (0.053) but not the two“D”s(0.153),indicating a larger differentiation between DNNAM and D-SNAM.The differentiation of the center of origin“O”from the secondary center“B”was not large(0.110-0.136),but was larger from the secondary center A-NCHN(0.158-0.196) and still larger from D-NNAM and D-SNAM(0.196-0.233). A2-SSWE was quite different from all other populations, with Fstranging from 0.279-0.422 and a mean of 0.354. Moreover, “D” centers differed greatly from Asian and African populations(0.200-0.212).

Values on the upper diagonal are similarity coefficients sijand those on the lower diagonal are differentiation coefficients Fst.

Thus, allele richness and allele frequency dispersion were usually higher,SPA and Fstwere not higher in centers of origin(“O”s) and older secondary centers (such as “B”s), but these genetic diversity values in younger secondary centers(such as“C”and“D”)and tertiary centers may be higher or lower than those in older centers. This finding means that alleles emerged or vanished after dissemination and adaptation,giving rise to changes in allele frequency and thus differentiation among 13 regional populations.

Geographic and genetic clustering revealed the relationship among all the populations on the basis of the genetic structure of each population. Fig. 2-B shows the results of clustering using genetic distances,dividing the 13 populations into five groups.(1)From the top,one group is composed only of the A2-SSWE population, which was genetically quite distinct from other groups. (2) The second group consists of all secondary and tertiary centers in the Americas, including D-NNAM,D-SNAM,and D1-CSAM,among which D-SNAM and D1-CSAM show a closer genetic relationship because the latter was developed based on the former. (3) The third group contains only the A-NCHN population included in the upper cluster with three American populations rather than centers of origin in China. This composition matches the results of geographic-ecological clustering and suggests a close relationship between major sources of American populations and A-NCHN. This explains why A-NCHN was not included in center of origin but nominated as a secondary center. (4) The fourth group was largest, consisting of six regional populations all from Asia. In this group, the center-of-origin populations O-HCHN and O-SCHN were included showing a closer relationship with C-SEAS, C1-SASI, and C2-AFRI than with A-NCHN. (5) The fifth group contains the secondary centers B-KORP and B-JPAN, which were genetically more distant from the centers of origin O-HCHN and O-SCHN than from the derived centers C-SEAS,C1-SASI,and C2-AFRI.

A-NCHN, A1-RUFE, and A2-SSWE were clustered in three different groups with A1-RUFE closer to“O”s than to A-NCHN.The reason may be that its source was mainly from the earlier accessions of northeast China, closer to O-HCHN, whereas some new introductions from other populations were included in the present northeast China population, whereas newly developed early-maturing accessions had not yet been introduced to A1-RUFE.

Table 2-Genetic diversity and differentiation of 13 regional populations.

3.3. Phylogeographic relationships among world soybean accessions

Clustering individual accessions may provide further information on the genetic structure of an ecological-population.In the clustering analysis of 371 world soybean accessions based on three geographic-ecological-traits, these accessions were grouped into six clusters characterized by different magnitudes of the three traits(Table 4,Fig.S2).(1)The largest group (group IV) was composed of 160 accessions insensitive to daylength and temperature and from all regional populations except the earliest-maturing population A2-SSWE and the latest-maturing populations C2-AFRI and D1-CSAM. (2)The second largest cluster(group VI)comprised 89 accessions sensitive to daylength and temperature and mainly from OSCHN, B-JPAN, D-SNAM, and D1-CSAM. (3) The third largest cluster(group V)consisted of 76 accessions very insensitive to daylength and temperature and mainly from A-NCHN, A1-RUFE, A2-SSWE, B-JPAN, and D-NNAM. (4) The fourth cluster(group III)comprised 34 accessions very sensitive to daylength and temperature and mainly from O-SCHN, C-SEAS, C1-SASI,C2-AFRI, and D1-CSAM. (5) The other two clusters (groups I and II) contained respectively only two and 10 accessions,which were sensitive to daylength and temperature.

In Table 4, O-SCHN, B-JPAN, C-SEAS, and C1-SASI are scattered across respectively 5, 4, 4, and 4 clusters among geographic-regional populations, indicating a wide range of daylength- and temperature-response features in these regions. However, O-HCHN, A-NCHN, A1-RUFE, A2-SSWE, DSNAM, and D1-CSAM were scattered across respectively only 2,2,2,1,2,and 2 clusters and concentrated on daylength-and temperature-insensitive or -sensitive clusters. By sensitivity to daylength and temperature, clusters ranked iii, ii, vi, i, iv,and v in descending order,while geographic regions ranked OSCHN (covering five clusters), B-JPAN, C-SEAS, B1-SASI (covering four clusters),B-KORP,C2-AFRI,D-NNAM(covering three clusters),O-HCHN,A-NCHN,A1-RUFE,D-SNAM,and D1-CSAM(covering two clusters), and A2-SSWE (covering only one cluster) in descending order of sensitivity. Thus, there was higher geographic-ecological variation in the O-SCHN and“B”groups and lower variation in the O-HCHN and “A” and “D”groups, possibly owing to the difference between multiple planting seasons and mainly summer-planting season.

Fig.3-The distribution of specific-present alleles(SPAs)in regional populations other than C2-AFRI and D1-CSAM.The vertical axis,labeled“Allele”,shows the specific-present alleles and the horizontal axis,labeled“Accession”,shows accessions of the 11 regional populations.Different colors represent the sources of regional populations with 247 SPAs in total.

Table 3-Genetic differentiation among 13 regional populations.

Table 4-The frequency distribution of geographic accessions in six phenotypic clusters of world soybeans based on three daylength- and temperature-sensitive traits.

In genotypic clustering based on individual accessions, the 371 accessions were classified into six genetically different groups(Table 5,Fig.4-A).(1)Cluster IV was the largest cluster,comprising 132 accessions from all geographic regions by A2-SSWE,including almost all accessions from the centers of origin O-HCHN and OSCHN.This cluster contained major parts of C-SEAS,C1-SASI,and C2-AFRI, including parts from A-NCHN, B-KORP, and B-JPAN,which maintained closer genetic relationships with accessions of the center of origin. (2) Cluster V was the second largest,containing 117 accessions from nine geographic regions,including a major part from A-NCHN, D-NNAM, and A2-SSWE. This result indicated that both American and European earlymaturing accessions could trace their sources to northeast China. (3) Cluster III comprised 47 accessions from seven geographic regions, including a major part from B-KORP, BJPAN, and A1-RUFE. B-KORP and B-JPAN germplasm exerted no obvious genetic influence on American accessions,given that not many“D”accessions were included in this cluster.(4)Cluster VI comprised 61 accessions from six geographic regions, mainly from D-SNAM and D1-CSAM and partly from D-NNAM,indicating that the accessions in D-SNAM had their germplasm from DNNAM but that the genetic structure of D-SNAM diverged from that of D-NNAM and formed a separate cluster along with D1-CSAM after intensive breeding.(5)Cluster I was very specific with only one accession from a specific type in Africa (there were others in cluster IV). (6) Cluster II contained 13 accessions including 11 from A-NCHN and two from A1-RUFE,all of which were early-maturing, but genetically different from earlymaturing ones in clusters III,IV,and V.

It is interesting that the accessions in the centers of origin O-HCHN and O-SCHN were concentrated in cluster IV,indicating the similar genetic bases of the two populations.However, accessions in other 11 geographic regions were distributed in several clusters, indicating that multiple sources were introduced from outside regions and that genetic differentiation occurred because of multiple introductions. In particular, O-SCHN covered only two clusters genetically, but five clusters geographic-ecologically. Accessions in A-NCHN and A1-RUFE were all early-maturing, but covered four clusters, suggesting that early-maturing accessions might have quite different genetic structures. Accessions in secondary centers B-KORP and B-JPAN were distributed in four or three clusters genetically and geographic-ecologically, indicating a large variation of their maturity groups and mixed genetic sources(Table 5).

Table 5-The frequency distribution of geographic accessions in six genotypic clusters of world soybeans.χ2= 727.3,χ20.01,60 =88.4, P ＜0.01.

These six genotypic clusters were further clustered. Fig. 4-B shows a close genetic relationship between clusters III and IV as well as clusters II and V.Most accessions of cluster IV were from O-HCHN, O-SCHN, C-SEAS, and C1-SASI, whereas cluster III comprised accessions mainly from B-KORP and B-JPAN,suggesting that soybeans were introduced from China into the Korean peninsula, Japan, Southeast Asia, and South Asia and that soybeans in the Korean peninsula and Japanese islands formed different germplasm pools later. The formation of new germplasm pools might be attributable to the long history of breeding efforts of historical farmers and modern breeders in the Korean peninsula and the Japanese islands. Another reason might be genomic introgression through interspecific hybridization,given the abundance of local wild soybeans naturally present in the Korean peninsula and the Japanese islands,as described by Xu et al.[38]and Wang et al.[39].Cluster V comprised 117 accessions,most of which were from high latitude regions including A1-RUFE,A2-SSWE,A-NCHN,and D-NNAM.Cluster VI comprised 61 accessions,of which 51(84%)came from D-SNAM and D1-CSAM.The differentiation between clusters V and VI suggests that soybeans in the south of North America and Central and South America formed a specific genetic structure after their introductions from northern North America.

4. Discussion

4.1.Features of the genotypic and phenotypic clustering of the world soybeans

Integrating genetic marker-haplotype (gene-allele) and geographic-ecological-trait clustering analysis, the present study focused on the same world soybean population. Both genotypic and phenotypic clustering put O-HCHN and OSCHN into the same cluster and also did so for A-NCHN and DNNAM; C-SEAS, C1-SASI, and C2-AFRI; B-KORP and B-JPAN;and D-SNAM and D1-CSAM (Fig. 2). Thus, 11 of the 13 geographic-populations were grouped into the same five sets of sensitivity-similar geographic-populations. Such consistency must reflect common factors, of which the major one might be the common genetic constitutions of daylength-and temperature-sensitive traits. In addition to consistent sets,however, some differences were also found, and may reflect different genetic constitutions of other traits.

Both approaches showed the largest clusters covering most of the regional populations (iv vs. IV), namely a cluster covering a major part of D-SNAM and D1-CSAM(vi vs.VI)and a cluster covering O-SCHN, A-NCHN, A1-RUFE, A2-SSWE, BKORP, B-JPAN, D-NNAM (v vs. V) in the phenotypic and genotypic clustering among the accessions. These three pairs of major clusters displayed similarities among the three traits(Table 4 and Table 5).However,the two clustering results also showed some differences. For example, three minor clusters in both clustering results had quite different components; O-SCHN covered five clusters in phenotypic clustering but only two clusters in genetic clustering,and the same disparity was observed for A-NCHN and A1-RUFE(2 vs.4 clusters) and C1-SASI (4 vs. 2 clusters). This finding suggests that different genetic constitutions may show similar daylength and temperature responses and vice versa.

The results of geographic-ecological and genetic clustering showed some consistency at a lower cluster level and some inconsistency at a higher cluster level. This finding may be explained as follows.First,the genotypic clustering reflected the relationship among accessions based on the whole genome,whereas geographic-ecological clustering involved only the part of the genomic information that conditioned the geographicecological traits. Second, multiple differing sources of introductions and germplasm exchanges might result in differential genetic constitutions and even similar phenotypes of geographicecological traits. Third, in addition to geographic conditions a variety of breeding goals are likely to have caused extra genetic variation.In summary,both geographic-ecological clustering and genotypic clustering shed light on the relationships among regional populations,and their results can be integrated to reveal the phylogeographic features of soybeans.

4.2. Historical dissemination paths of soybeans summarized from present phylogeographic analysis

Historical dissemination paths of soybeans may be summarized by integration of the phylogeographic results obtained from the analysis of geographic-ecological-traits and molecular markerhaplotype(gene-allele)traits with information in the literature.

Fig.4- Unrooted trees showing genetic relationships among 371 world soybean accessions. (A)Unrooted tree showing the genetic relationships among 371 world soybean accessions. Colors represent the sources of regional populations.Roman numerals beside vertical bars indicate six clusters assembled from 371 accessions;(B)Second-order genotypic clusters based on the first-order clusters.

The first path of dissemination was from the center of origin to northeast China, although previous researchers [10]considered northeast China as part of the center of origin.However,though it is part of China,it is not necessarily part of the center of origin of soybean.Population A-NCHN clustered apart from O-HCHN and O-SCHN with the earliest-maturing groups and SPAs more than those in the centers of origin.For this reason, we assigned it as the first secondary center ANCHN.A-NCHN,A1-RUFE,and A2-SSWE were not in the same cluster in Fig. 2-B, but shared the same genetic cluster V in Table 5, while A1-RUFE was very similar to A-NCHN in clustering distribution (Table 5, cluster II, III, IV, and V).Accordingly,A-NCHN,A1-RUFE,and A2-SSWE were placed in Group“A”,with the latter two as populations derived from ANCHN. Here, it was understood that the early soybeans in ANCHN were introduced mainly from the center of origin,where very early-maturing varieties, such as MG 00 and MG 000, were developed intensively and independently along with introductions from other regions used in breeding programs, which caused the present A-NCHN, A1-RUFE, and A2-SSWE grouped into different cluster branches.This grouping is in accord with a previous [14] grouping of Russia and northeast China as a soybean pool, Swedish and Russian varieties maintain a close relationship [11] and several Swedish accessions originated in northern Manchuria(northern northeast China)[40].

The second path of dissemination was from the center of origin eastward to B-KORP and then to B-JPAN from about the first century A.D.Multiple paths were taken,one through the border area between northeast China and the Korean peninsula and others via ocean voyages, which possibly brought soybeans directly to B-KORP or B-JPAN. These two secondary“B”s span the same range of maturity groups 0-VII even in a small territory, with a wide range of flowering and maturity dates,they were genetically clustered in a separate cluster apart from that of the center of origin (Fig. 2-B) and contained wider genetic clusters(II,IV,V,and even VI).These secondary“B”s adapted to maritime conditions and formed a new geographic group. In the literature [3,21,41], Korean peninsula and Japanese soybeans are genetically similar,but different from Chinese soybeans, in agreement with the present finding that B-KORP and B-JPAN shared a common origin but are genetically grouped into a cluster other than that containing the“O”population.

The third path of dissemination was from the center of origin southward to Southeast Asia via commerce and migration, afterwards to South Asia, and then to Africa (CSEAS, C1-SASI, and C2-AFRI). Soybeans were first cultivated in Africa in 1896 [12]. They were grown first in Algeria and then introduced to South Africa,Mauritius,and Tanzania.We supposed that this was a joint path with B-KORP or B-JPAN.However, geographic-ecological data showed that their maturity groups were mainly later ones (MG III-IV and MG VII-X in C-SEAS, MG VI-X in C1-SASI and C2-AFRI) and adapted to short daylength and warm temperature,with the greatest number of new SPAs in C1-SASI. These derived populations were genetically clustered mainly in the same cluster(Group IV)with centers of origin O-HCHN and O-SCHN(Fig. 2-B, Table 5), suggesting that these three populations showed small genetic differences from “O”s. Thus, C-SEAS was a secondary center and the other two were tertiary centers.

The fourth or the most recent dissemination path was mainly from China and later Japan westward to Europe with missionaries and sailors,with soybeans then grown in France and England in the late 18th century [42]. Soybeans were introduced by Samuel Bowen from China into the USA via London in 1765 [43]. In 1893, Zavitz introduced soybeans into Canada, and bred forage soybean cultivars [11]. Afterwards,soybean germplasm was reintroduced from China several times, especially from northeast China to the USA. The first soybean accessions were introduced from the southern USA into Brazil around the end of the 19th century[13],followed by the introduction of some accessions from Philippines, Thailand and Pakistan.Soybeans were brought from Brazil directly to Argentina and other countries. With suitable climate and soil conditions, Brazil and Argentina became large soybean producers. In the meantime, breeding progress was made,contributing to the formation of a potential new gene center of soybeans there.The present study genetically clustered the secondary centers D-NNAM and D-SNAM and the tertiary center D1-CSAM in the same group with A-NCHN (Fig. 2-B)grouped with the major part(55 of 64 accessions)of D-NNAM in the same major group V,meaning that the germplasm of DNNAM, D-SNAM, and D1-CSAM was mainly from northeast China. Thus, the fourth path of dissemination was actually from the center of origin to A-NCHN via Europe or directly to North America, in general agreement with the literature [44],which indicates that many accessions in USDA soybean collection were from China and the Korean peninsula, with these introductions forming the basis of North American germplasm. These introductions were also used in the breeding programs of southern North America.After adapting to southern North America, D-SNAM was different from DNNAM.Its major cluster was accordingly Group VI in Table 5.Soybeans in Central and South America were originally introduced from D-SNAM and formed the tertiary center D1-CSAM after a few decades.

Thus, different paths of dissemination resulted in adaptation to different geographic-ecological environments, improvements in meeting local requirements and the formation of different geographic-populations. Among the paths, the dissemination to the north and then to the west (Americas)was the most important because the soybean production of the western world began with the soybean germplasm of northeast China. However, the southward path may represent the starting point for breeding later-maturing soybeans in tropical and subtropical areas.

5. Conclusions

World soybean accessions numbering 371 from 13 geographic ecological-regions, their phylogeographic features, and the relationships among their geographic-populations were studied using two sets of trait data: molecular gene/allele or marker/haplotype traits and geography-related photo- and temperature-sensitive traits (flowering date, maturity date,and main stem node number). The latitude (daylength and temperature) of the testing site Nanjing in China allows all materials to blossom and mature or approximately mature naturally,favoring the comparability of geographic-ecological characteristics. The results confirmed the historical differentiation forming geographic populations,including two centerof-origin populations,four secondary centers,and five tertiary centers with different genetic and geo-ecological components.The results of both genetic and phenotypic clustering showed that 11 of 13 geographic-populations were grouped into the same five sets of sensitivity-similar geographic-populations,indicating that geographic-related photo- and temperaturesensitive traits were influential in determining the genetic differences among geographic-populations. Both clustering results grouped the population of northeast China into a secondary center rather than the center-of-origin populations together with the population of northern North America,correcting the conventional model considering it part of the center of origin.The phylogeographic analysis integrated with historical information also provided an explanation for the model of four dissemination paths: from centers of origin to the north (northeast China, Russia, and Sweden/Europe in turn), the east (Korean peninsula and Japan), the south(Southeast Asia, South Asia and Africa) and from northeast China to the west (Europe, northern USA, southern USA, and South America). Among these, the most recent path of dissemination to the Americas via northeast China was the most influential, given that the soybean production of the Western world has grown rapidly in the past half century,accounting for around 90% of current world soybean production. Thus, the soybean geographic-population of northeast China was the one most profoundly influencing the development of soybean production worldwide in modern society.

Declaration of competing interest

None declared.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (2017YFD0101500,2016YFD0100304), the National Natural Science Foundation of China (31701447, 31671718, 31571695), the Ministry of Education 111 Project (B08025), the Ministry of Education Program for Changjiang Scholars and Innovative Research Team in University (PCSIRT_17R55), the China Agriculture Research System (CARS-04), the Priority Academic Program Development of Jiangsu Higher Education Institutions, the Fundamental Research Funds for the Central Universities and the Jiangsu Collaborative Innovation Center for Modern Crop Production, Natural Science Foundation of Jiangsu Higher Education Institutions(19KJB210010).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2019.09.010.