Identification of Non-Pleiotropic Loci in Flowering and Maturity Control in Soybean

Pleiotropy is considered to have a significant impact on multi-trait evolution, but its roles in the evolution of domestication-related traits in crop species have been unclear. In soybean, several known quantitative trait loci (QTL) controlling maturity, called the maturity loci, are known to have major effects on both flowering and maturity in a highly correlated pleiotropic manner. Aiming at the identification of non-pleiotropic QTLs that independently control flowering and maturity and dissecting the effects of pleiotropy in these important agronomic traits, we conducted a QTL mapping experiment by creating a population from a cross between domesticated soybean G. max and its wild ancestor G. soja that underwent stringent selection for non-pleiotropy in flowering and maturity. Our QTL mapping analyses using the experimental population revealed novel loci that acted in a non-pleiotropic manner: R1-1 controlled primarily flowering and R8-1 and R8-2 controlled maturity, while R1-1 overlapped with QTL, affecting other agronomic traits. Our results suggest that pleiotropy in flowering and maturity can be genetically separated, while artificial selection during soybean domestication and diversification may have favored pleiotropic loci such as E loci that control both flowering and maturity. The non-pleiotropic loci identified in this study will help to identify valuable novel genes to optimize soybean’s life history traits and to improve soybean’s yield potential under diverse environments and cultivation schemes.


Introduction
Genetic pleiotropy, correlation among traits due to genes that possess pleiotropic roles or the linkage of multiple genes, is known to have a major impact on the multi-trait evolution of organisms. Neutral theories and evolutionary models suggest that pleiotropy is a cause of evolutionary constraints, preventing organisms from evolutionary progress [1,2]. In contrast, recent genome-wide observations indicate that pleiotropy promotes the evolution of complexity [3,4]. Crop domestication provides a unique opportunity to dissect the effects of pleiotropy on multi-trait evolution. Intense selection by humans during crop domestication and diversification processes causes rapid and directional changes in a set of traits, typically including non-seed shattering, reduced branching, apical dominance, large fruit or grain size, and synchronized flowering, which is known as the domestication syndrome [5]. However, much remains unknown about the extent of pleiotropy among domestication-related traits and its consequences in the evolution of crop species.
Cultivated soybean (Glycine max) is considered to be domesticated from its wild ancestor Glycine soja approximately 6000 to 9000 years ago [6,7]. Preceding the domestication event, recent discoveries suggest a prolonged period of low intensity management or semi-domestication of wild soybeans at multiple locations in East Asia [8]. Compared with its ancestor, G. max shows reduced pod dehiscence, branching and lodging, larger seeds, and more synchronized flowering and maturation. A handful of genes that are implicated in soybean domestication, diversification, and improvement processes have been identified. Among them, three genes are thought to have played central roles at the early stage of soybean domestication: SHAT1-5 controlling pod dehiscence [9,10], GmHs1-1 controlling seed hardness [11], and GmFT2c controlling photoperiodic flowering [12].
An emphasis of soybean domestication and diversification was selection for plants adapted to regional photoperiods that confines soybean varieties carrying differing photoperiod sensitivities for flowering and maturity to specific latitudinal boundaries [13,14]. Several QTL controlling diversity in flowering and maturity among cultivated soybean varieties have been identified: the maturity loci E1-E10 and J [15,16]. The E1 locus encodes a transcription factor carrying a B3 domain [17], E2 is an ortholog of the Arabidopsis flowering regulator GIGANTEA [18], and E3 and E4 are homologs of the photoreceptor PHYTOCHROME A (PHYA) [19,20]. Recently, additional maturity loci E9 and J are identified as the FLOWERING LOCUS T homolog GmFT2a [21] and the circadian clock gene EARLY FLOWERING 3 [22], respectively. Major maturity loci are observed to affect both flowering and maturity in a highly correlated manner [23,24]: early flowering is associated with early maturity, while late flowering is associated with late maturity, suggesting that the pleiotropy of these loci controlling both flowering and maturity may have been a limiting factor for the diversification of soybean's life cycle in breeding programs.
In order to dissect the effects of pleiotropy in the evolution of flowering and maturity control in soybean, in this study, we conducted a selection experiment by creating a population from a cross between G. max and G. soja that underwent stringent selection for non-pleiotropy in flowering and maturity, and carried out QTL mapping analysis to identify novel loci that independently control flowering or maturity in a non-pleiotropic manner.

Mapping Population
The initial cross between the G. max reference variety Williams 82 (Maturity group III) and the G. soja accession PI 549046 (Maturity group III) was made in 1996 ( Figure 1). The F1 plants were backcrossed to Williams 82 in 1997, and the BC1F1 plants were grown in 1998. From the 3000 BC1F2 plants, approximately 60 plants were selected based on agronomic appearance similar to Williams 82. From the 60 BC1F3 lines, 87 plants out of 17 BC1F3 lines with the best agronomic appearance similar to Williams 82 were selected in 2000. In the BC1F4 generation, nine lines that showed agronomic appearance similar to Williams 82 were selected and bulk harvested. These lines were grown in the field, and phenotypes were observed in 2001 and 2002. In order to identify non-pleiotropic loci that control either flowering or maturity independently unlike the known major maturity loci that control both flowering and maturity in a highly pleiotropic manner, our selection scheme in this study targeted a type of non-pleiotropy affecting flowering but not maturity. We selected the lines LG01-7909 and LG01-7919 that flowered 10 to 12 days later and matured 7 to 13 days earlier than Williams 82 in 2001, and flowered 9 to 11 days later and matured 11 to 12 days earlier than Williams 82 in 2002. In 2003, these two BC1F5 lines were backcrossed again to Williams 82, and the BC2F1 plants were grown in the following winter in Puerto Rico. These populations were advanced through single-seed descent from the F2 to the F4 generation. In the BC2F4 generation, 109 plants were harvested from the population

Field Evaluation
The BC2F6 RILs, the parent line Williams 82, and maturity checks (IA2023 and LD00-3309) were evaluated at Urbana and Stonington, Illinois in 2009 and at Urbana, Villa Grove, Bellflower, and Stonington, Illinois in 2010. All field plots were four rows wide with 0.76 m spacing and 3.1 m long and were planted at a rate of 30 seeds per meter. At all locations and in both years, the lines were planted in a randomized complete block design with two replications. The phenotype data collected include flowering date (R1), which was recorded when approximately 50% of the plants had at least one flower [25], plant maturity date (R8), which was recorded when approximately 95% of the pods had reached mature pod color [25], plant lodging scored at maturity based on a 1 to 10 scale (1 = all plants are erect, 10 = all plants are prostrate), plant height (cm) measured from the soil surface to the top node of the main stem and seed yield (kg ha −1 ) recorded at 13% moisture, and a stem vininess score between 1 (typical

Field Evaluation
The BC2F6 RILs, the parent line Williams 82, and maturity checks (IA2023 and LD00-3309) were evaluated at Urbana and Stonington, Illinois in 2009 and at Urbana, Villa Grove, Bellflower, and Stonington, Illinois in 2010. All field plots were four rows wide with 0.76 m spacing and 3.1 m long and were planted at a rate of 30 seeds per meter. At all locations and in both years, the lines were planted in a randomized complete block design with two replications. The phenotype data collected include flowering date (R1), which was recorded when approximately 50% of the plants had at least one flower [25], plant maturity date (R8), which was recorded when approximately 95% of the pods had reached mature pod color [25], plant lodging scored at maturity based on a 1 to 10 scale (1 = all plants are erect, 10 = all plants are prostrate), plant height (cm) measured from the soil surface to the top node of the main stem and seed yield (kg ha −1 ) recorded at 13% moisture, and a stem vininess score between 1 (typical indeterminate cultivar) and 5 (viney stem termination). The 115 BC2F6 RIL population was divided by maturity based on the data collected in 2008. Lines that matured between 18 and 24 Agronomy 2020, 10, 1204 4 of 15 September were blocked together (Maturity set 1) and those that matured between 25 and 30 September were blocked together (Maturity set 2) for further analysis.

Statistical Analysis of Phenotypic Data
Statistical analysis of phenotypic distribution and multiple factor analysis (MFA) were carried out using R scripts. Box plots were created using the Psych package [26] and multiple factor analysis (MFA) was conducted using the FactominR package [27]. Correlation analysis was conducted in R.

Genotyping
A total of 1536 single nucleotide polymorphisms (SNPs) from the 1536 Universal Soy Linkage Panel (USLP 1.0) [28] was used for genotyping assay. A single leaf from a single plant was collected for each line between 30 and 40 days after sowing. DNA samples from 115 BC2F6 RILs and Williams 82 were obtained using the DNeasy Plant Mini kit (QIAGEN), and genotyping was carried out using the Illumina GoldenGate ® Assay [29] by Perry Cregan at USDA-ARS, Maryland.
The E1 alleles were genotyped as described previously [30], and the SNP marker BARC-056323-14257 was genotyped via the direct sequencing of fragments amplified by polymerase chain reaction using the primers 5 -ACCATCATTGTTGTGAACCCTAC-3 and 5 -AACGTCATCCACTTGAAACTTGG-3 .

SNP Selection
SNPs of interest were selected using Illumina's Genome Studio ® 1.0. SNPs with low call numbers (below 5%) or no segregation in our mapping population were excluded. In addition, SNPs for which both parents shared the same allele type were removed. Of the 1536 SNPs genotyped, 545 were selected for quantitative trait loci (QTL) mapping (Table S1).

QTL Mapping
Marker distances and linkages were calculated using the Kosambi function in JoinMap ® 4.0 with a minimum logarithm of the odds (LOD) threshold of 3.0. QTL analysis was conducted using composite interval mapping (CIM) in the software Windows QTL Cartographer 2.5 [31]. Stepwise selection was used to identify background co-factors in the model. For CIM, a 1000-iteration permutation test was performed for all traits to generate a genome-wide critical alpha of 0.05. To control the background marker effect, a window size of 1cM was employed.

Candidate Gene Identification
The physical locations of significant flowering and maturity SNPs were compared with a list of flowering-related genes and previously identified flowering and maturity QTL (http://soybase.org/). Candidate genes were called if they were located within 500 kb of QTL-associated SNPs.

Plant Growth Condition
The flowering time and maturation time measurement of the BC3F2 population was conducted in the greenhouse during November through February under a moderate short-day photoperiod (12 h light) with supplementary white fluorescence light. Highest temperatures during the day ranged between 24.0 and 26.5 • C, and lowest temperatures in the evening ranged between 20.0 and 22.5 • C.

Selection for Non-Pleiotropy
The G. max reference variety Williams 82 and the G. soja accession PI 549046 were used as the parents to create the mapping population evaluated in this study ( Figure 1). PI 549046 was one of the lines selected from four Chinese provinces that were genetically the most distinct from cultivated soybean [32]. Mimicking the domestication and diversification processes of soybean, two backcrosses to Williams 82 were performed after the initial cross between Williams 82 and PI 549046, and plants that showed agronomic appearance similar to Williams 82 in traits such as lodging, height, and stem vininess were selected through multiple generations. As an attempt to identify non-pleiotropic loci, we implemented a selection scheme for a type of non-pleiotropy that affects flowering but not maturity; in the BC1F5 generation, we selected the lines LG01-7909 and LG01-7919 that flowered significantly later but matured moderately earlier than Williams 82 in both 2001 and 2002. These lines were backcrossed again to Williams 82 (BC2) and advanced to the F5 generation. Then, we further selected lines that showed a maturity time similar to Williams 82, resulting in a total of 115 BC2F5 RILs with a narrow maturity range of 12 days. Strong positive correlation was observed in all comparisons between different locations for each trait, Maturity set, and year ( Figure S1), except for height in Urbana and Stonington in Maturity set 2 in 2009, which appeared marginal (p = 0.088).

Phenotypic Distribution
The 115 BC2F6 RIL population was used for phenotyping at multiple field locations: Urbana and Stonington, Illinois in 2009 and Urbana, Villa Grove, Bellflower, and Stonington, Illinois in 2010. All traits displayed a strong yearly and/or location effect in the distribution (Figure 2a Aiming at identifying loci that affect flowering in a non-pleiotropic manner, the 115 RILs were divided into two sets for further analysis based on their maturity: early maturity (Maturity set 1; n = 64) and late maturity (Maturity set 2; n = 51). The phenotypic distribution differed significantly between maturity sets 1 and 2 in flowering, maturity, and height ( Figure 2c); a difference of 4.6 days in flowering (p = 2.12 × 10 −7 ), 4.8 days in maturity (p = 8.19 × 10 −27 ), and 9.2 cm in height (p = 1.48 × 10 −3 ). In contrast, lodging, stem vining, or yield did not differ between maturity sets 1 and 2. No significant difference in phenotypic distribution was observed between the LG01-7909 progenies (Subpopulation 1; n = 68) and the LG01-7919 progenies (Subpopulation 2; n = 47) in any traits.

Multiple Factor Analysis (MFA)
Multiple factor analysis (MFA) further demonstrated the different relationship among the traits ( Figure 3). Dimensions 1 and 2 explain 45.3% and 23.1% of the total variation, respectively, and Dimension 3 makes up 16.7%. The traits are observed to cluster into three groups. Height, stem vining, and lodging cluster together and make up major components in Dimension 1 (0.84, 0.80, and 0.86, respectively) ( Figure 3; Table 1), while flowering and maturity times are represented together in Dimension 2 (0.82 and 0.56, respectively), consistently with the previous implication regarding the highly correlated behavior of maturity loci in flowering and maturity [23,24,33,34]. Finally, yield is the primary component in Dimension 3 (0.83). Aiming at identifying loci that affect flowering in a non-pleiotropic manner, the 115 RILs were divided into two sets for further analysis based on their maturity: early maturity (Maturity set 1; n = 64) and late maturity (Maturity set 2; n = 51). The phenotypic distribution differed significantly between maturity sets 1 and 2 in flowering, maturity, and height ( Figure 2c); a difference of 4.6 days vining, and lodging cluster together and make up major components in Dimension 1 (0.84, 0.80, and 0.86, respectively) ( Figure 3; Table 1), while flowering and maturity times are represented together in Dimension 2 (0.82 and 0.56, respectively), consistently with the previous implication regarding the highly correlated behavior of maturity loci in flowering and maturity [23,24,33,34]. Finally, yield is the primary component in Dimension 3 (0.83).

Construction of Linkage Map
Genetic linkage maps were constructed for this mapping population, resulting in 30 linkage groups with 335 SNP markers spanning 846 cM ( Figure S2). Marker locations were generally consistent with the Universal Soy Linkage Panel 1.0 [29]; however, we found that chromosomes 1, 5, 6, 7, 12, 13, 14, and 16 were divided into two or three linkage groups in our maps. This is likely caused by the two backcrosses to Williams 82 and intense selection, limiting the number of SNP markers segregating in the population. On average, each linkage group contained 12 markers with a mean distance of 3.6 cM between the nearest markers. Extensive segregation distortion was observed, resulting in G. soja alleles being substantially underrepresented for many markers in our population. Each RIL possessed approximately 3-20% of SNP markers homozygous for G. soja alleles (Table S1). Linkage group 7a contained the highest percentage of markers homozygous for G. soja alleles at 64%, and linkage group 12a contained the least at 1.3%.  Table 2). One QTL was found for flowering, two were found for maturity, two were found for height, two were found for yield, seven were found for lodging, and five were found for stem vining. Most of these QTL were found in specific locations or subpopulations.

Construction of Linkage Map
Genetic linkage maps were constructed for this mapping population, resulting in 30 linkage groups with 335 SNP markers spanning 846 cM ( Figure S2). Marker locations were generally consistent with the Universal Soy Linkage Panel 1.0 [29]; however, we found that chromosomes 1, 5, 6, 7, 12, 13, 14, and 16 were divided into two or three linkage groups in our maps. This is likely caused by the two backcrosses to Williams 82 and intense selection, limiting the number of SNP markers segregating in the population. On average, each linkage group contained 12 markers with a mean distance of 3.6 cM between the nearest markers. Extensive segregation distortion was observed, resulting in G. soja alleles being substantially underrepresented for many markers in our population. Each RIL possessed approximately 3-20% of SNP markers homozygous for G. soja alleles (Table S1). Linkage group 7a contained the highest percentage of markers homozygous for G. soja alleles at 64%, and linkage group 12a contained the least at 1.3%.

QTL Mapping
Composite interval mapping (CIM) analysis was carried out on both the data in its entirety and  Table 2). One QTL was found for flowering, two were found for maturity, two were found for height, two were found for yield, seven were found for lodging, and five were found for stem vining. Most of these QTL were found in specific locations or subpopulations.  The QTL R1-1 on chromosome 9 was significantly associated with flowering (R1). This QTL explains 26% of the flowering variation observed and is centered on the SNP marker BARC-056323-14257 with an LOD score of 4.10 and an additive effect of 6.41. Despite the late flowering parent lines BC1F5 LG01-7909 and LG01-7919 that were selected for the second backcrossing, this result indicates that the G. soja allele confers earlier flowering than the G. max allele, which is potentially due to epistatic interactions with other loci in the genome. Two QTL, R8-1 and R8-2, were associated with maturity (R8) and located near the markers BARC-041649-08056 and BARC-015067-02556 on chromosomes 13 and 18, respectively. R8-1 explains 32% of the variation observed with an LOD score of 5.88. R8-2 explains 17% of the variation observed with an LOD score of 3.14. The additive effects of R8-1 and R8-2 are −1.31 and −1.18, respectively, with the G. soja allele conferring faintly later maturation than the G. max allele.
Two QTL, HT-1 and HT-2, were associated with height and located near the markers BARC-056323-14257 and BARC-040851-07854 on chromosomes 9 and 11, respectively. HT-1 explains 10% of the variation observed with an LOD score of 3.15. HT-2 explains between 15% and 26% of the variation observed with a LOD score between 4.92 and 5.48. The additive effects of HT-1 and HT-2 are above 8.29 with both G. soja alleles conferring a large reduction in height.
Two QTL, YD-1 and YD-2, were associated with yield and located near BARC-014395-01348 and BARC-060195-16470, respectively, both on chromosome 18. YD-1 explains 15% of the variation observed with an LOD score of 4.71. YD-2 explains 13% of the variation observed with an LOD score of 3.77. The additive effects for YD-1 and YD-2 are 3.83 and 3.60, respectively, with the G. soja allele conferring a reduction in yield.
Seven QTL were associated with lodging, locating on chromosomes 9, 11, 13, 14, and 16. R 2 values ranged from 13% to 20%, with LD-7 explaining the highest percentage of the variation at 20%. LOD scores spanned between 3.18 and 5.36, with LD-3 showing the strongest peak. The additive effects of QTL associated with lodging range between −0.52 and 1.39, with the majority of G. soja allele conferring a decline in lodging. LD-4 and LD-6 are the exceptions to this trend.
Five QTL were associated with stem vining, locating on chromosomes 2, 6, 9, 11, and 14. R 2 values ranged from 9% to 24% variation explained, with SV-1 showing the highest value. LOD scores were spread between 2.90 and 4.96, with SV-5 showing the strongest peak. The additive effects of QTL affecting stem vining ranged from −0.42 to 1.36. With the exception of SV-5, all G. soja alleles confer a reduction in stem vining.

Pleiotropic QTL and Verification
Unlike the pleiotropic E loci, the flowering QTL R1-1 identified in this study did not overlap with the maturity QTL R8-1 and R8-2, nor with previously reported maturity QTL (Figure 4; Table 2). However, the R1-1 locus overlapped with the QTL controlling other agronomic traits: HT-1 and LD-2, co-localizing on chromosome 9. We found another cluster of QTL controlling multiple traits: HT-2, LD-3, and SV-4 on chromosome 11.
To verify the QTL identified in this study, previously reported QTL listed in SoyBase (http://soybase. org/) were compared with the locations of our QTL. QTL for flowering (R1-1), height (HT-1 and HT-2), and yield (YD-1 and YD-2) identified in this study, as well as one of maturity QTL (R8-2) and 3 of lodging QTL (LD-1, LD-3 and LD-6), overlapped with previously reported QTL (Figure 4; Table 3). The stem vining QTL SV-1 overlaps with the QTL TH-C for twinning habit [35], but none of the other stem vining QTL has been previously reported. To verify our QTL controlling flowering or maturity, we compared their locations with those of known flowering gene homologs [8]. Two genes homologous to Phytochrome E (PHYE) and TFL1/FT that are known to play a role in flowering control in Arabidopsis are located within the close range of the flowering QTL R1-1 ( Figure 4; Table 4). In addition, one of the maturity QTL, R8-1, appears to be in the vicinity of homologs of the flowering-related genes SRF6 and SPA3.
To provide further verification for the flowering QTL R1-1, the RIL LG08-7379 carrying the G. soja allele of the peak marker BARC-056323-14257 was backcrossed with Williams 82. In the F2 population, the plants carrying the homozygous G. soja allele at the marker flowered earlier than the plants carrying the homozygous G. max allele by 5.8 days (p < 0.01), while the maturation time of these plants was statistically indifferent ( Figure 5).    To provide further verification for the flowering QTL R1-1, the RIL LG08-7379 carrying the G. soja allele of the peak marker BARC-056323-14257 was backcrossed with Williams 82. In the F2 population, the plants carrying the homozygous G. soja allele at the marker flowered earlier than the plants carrying the homozygous G. max allele by 5.8 days (p < 0.01), while the maturation time of these plants was statistically indifferent ( Figure 5).

Discussion
Pleiotropy is considered to have significant implications for multi-trait evolution. In this study, we focused on the key domestication-related traits flowering and maturity that were known to behave in a highly correlated manner in domesticated soybean. Our selection experiment for non-pleiotropy resulted in a different evolutionary trajectory from that of soybean domestication. Our QTL mapping analysis using the experimental populations did not identify the E loci that were selected during soybean domestication and improvement, but it did identify novel loci affecting either flowering or maturity in a non-pleiotropic manner. This result suggests that artificial selection during the soybean domestication and diversification processes favored pleiotropic loci over non-pleiotropic loci in the evolution of flowering and maturity control, despite the fact that genetic variations in the wild progenitor that allow diverse patterns of life history are available.
In addition to flowering and maturity, pleiotropy likely played an important role in the evolution of multiple agronomic traits. A previous study reported a QTL on chromosome 11 in a mapping population derived from G. max parents that affected flowering, maturity, height, and lodging, and as well as another QTL on chromosome 6 that affected flowering, height, and lodging [36]. The latter locus maps near the known major maturity locus E1. The E1 locus has also been identified in multiple mapping populations derived from a cross between G. max and G. soja parents, and it affects several agronomic traits [35,37]. Our QTL mapping did not identify this locus, which was likely due to the selection for non-pleiotropy between flowering and maturity. Indeed, we found zero lines in Subpopulation 2 and only 2 homozygous and 1 heterozygous lines in Subpopulation 1 carrying the G. soja E1 allele (Table S2), suggesting that our selection scheme effectively purged pleiotropic loci such as E1 from the population. In addition, segregation distortion may have played a role in removing certain G. soja alleles, including E1, from the population. Multiple studies have reported segregation distortion in mapping populations derived from a cross between G. max and G. soja that results in disproportionally low percentages of G. soja alleles at heterologous segregating loci [35]. The causes of segregation distortion are yet to be clarified, but in our selection scheme, it is possible that selection for agronomic appearance while advancing generations through single seed descent, as well as selection for a narrow maturity range at the BC2F5 stage, may have skewed the proportion of G. soja alleles. Despite the strict selection scheme for agronomic appearance and segregation distortion, our QTL mapping found two genomic locations controlling multiple agronomic traits (Figure 4), suggesting the prevalence of pleiotropy among loci controlling agronomic traits. It has been reported that artificial selection during domestication causes rapid evolutionary changes [38]. Selection for prevalent pleiotropic loci may have accelerated the evolution of cultivated soybean during the domestication and diversification processes.
While our selection experiment successfully identified non-pleiotropic novel loci, it allowed only a small fraction of genetic variation subjected to the QTL mapping analysis due to our selection scheme and the biparental nature of our mapping population. In addition, the relatively small size of our mapping population may have limited the power of QTL mapping. Future experiments to follow-up this work include parallel evolutionary experiments under different selection regimes using expanded mapping populations created from diverse parental varieties in G. max and G. soja. The use of a larger mapping population, as well as a gentler selection scheme for agronomic appearance would also help maintain the level of G. soja alleles in a mapping population.
Another limitation of our experiment is that it does not assess the magnitude of factors other than pleiotropy that account for phenotypic evolution, nor their interaction with pleiotropy. For example, artificial selection toward pleiotropic loci may be in part a result of widespread epistatic interactions in the genome, which are shown to have a prominent effect on narrowing the evolutionary path [39]. Moreover, the possibility in which pleiotropy was favored as an indirect result of selection for gene-by-environment effects that eliminates alleles sensitive to environmental variants cannot be ruled out, as genetic loci that exhibit phenotypic robustness to diverse environments are shown to have a tendency to affect multiple traits [38]. Indeed, several QTL identified in this study appeared under specific years or locations.

Conclusions
This study represents one of the first attempts to experimentally dissect the roles of pleiotropy in the multi-trait evolution of crops. Our work demonstrates that pleiotropy in flowering and maturity can be genetically separated. Non-pleiotropic loci identified in this work provide novel genetic resources for soybean breeding programs and ultimately allow the diversification of soybean's growth habit and life history optimized for regional environments and different cultivation strategies.