Genome-Wide Association Studies for Dynamic Plant Height and Number of Nodes on the Main Stem in Summer Sowing Soybeans

Plant height (PH) and the number of nodes on the main stem (NN) serve as major plant architecture traits affecting soybean seed yield. Although many quantitative trait loci for the two traits have been reported, their genetic controls at different developmental stages in soybeans remain unclear. Here, 368 soybean breeding lines were genotyped using 62,423 single nucleotide polymorphism (SNP) markers and phenotyped for the two traits at three different developmental stages over two locations in order to identify their quantitative trait nucleotides (QTNs) using compressed mixed linear model (CMLM) and multi-locus random-SNP-effect mixed linear model (mrMLM) approaches. As a result, 11 and 13 QTNs were found by CMLM to be associated with PH and NN, respectively. Among these QTNs, 8, 3, and 4 for PH and 6, 6, and 8 for NN were found at the three stages, and 3 and 6 were repeatedly detected for PH and NN. In addition, 34 and 30 QTNs were found by mrMLM to be associated with PH and NN, respectively. Among these QTNs, 11, 13, and 16 for PH and 11, 15, and 8 for NN were found at the three stages. A majority of these QTNs overlapped with the previously reported loci. Moreover, one QTN within the known E2 locus for flowering time was detected for the two traits at all three stages, and another that overlapped with the Dt1 locus for stem growth habit was also identified for the two traits at the mature stage. This may explain the highly significant correlation between the two traits. Our findings provide evidence for mixed major plus polygenes inheritance for dynamic traits and an extended understanding of their genetic architecture for molecular dissection and breeding utilization in soybeans.


INTRODUCTION
Soybean (Glycine max L. Merr.) is an important source of plant protein and oil for human consumption. Improving seed yield is the major target for soybean breeders. Plant architecture can strongly affect the suitability and productivity of seed yield in agricultural crops (Li R. et al., 2014). Plant height (PH) and the number of nodes on the main stem (NN) as key plant type traits have obvious effects on the seed yield of soybean because they are related to some important characteristics such as lodging and adaptability (Chapman et al., 2003;Liu et al., 2011). PH and NN are highly correlated with some other soybean agronomic traits such as days to flowering (DF) and days to maturity (DM), which are thought to be mainly adaptive traits in response to the photoperiod, allowing each cultivar to adapt to limited geographic regions (Zhang et al., 2004Panthee et al., 2007).
Plant height and NN are complex traits governed by many quantitative trait loci (QTL) in soybeans (Lee et al., 1996;Zhang et al., 2004;Liu et al., 2011Liu et al., , 2013Yao et al., 2015;Cao et al., 2017). To date, more than 200 and 30 QTLs for PH and NN have been reported on SoyBase 1 via linkage mapping. Recently, a study showed that highly significant correlations were observed among yield-related traits such as PH and NN. In addition, 23 novel QTLs and 8 QTL hotspots were identified for yield and qualityrelated traits by QTL analysis in a soybean RILs population. In particular, most loci associated with these traits were co-located in the same genomic region on three chromosomes (Chr04, Chr06, and Chr19), which was consistent with the results of phenotypic correlation analysis . Fang et al. (2017) identified 245 loci for 84 agronomic traits via genome-wide association studies (GWAS) in 809 soybean accessions and further dissected the genetic networks underlying the phenotypic correlations of traits. Of these traits, PH and NN exhibited a significant positive correlation. Some major genes were also cloned to reveal the molecular mechanism of PH and NN. Two known loci, Dt1 for stem growth habit and E2 for DF, were involved in regulating PH and NN and other agronomic traits in soybeans (Kato et al., 2015;Zhang et al., 2015). Dt1 plays a primary role in determinate stem varieties and has an epistasis effect on the Dt2 locus, another stem growth habit locus involved in the development of PH in soybeans (Bernard, 1972;Liu et al., 2016). The E2 locus encodes a homolog of GIGANTEA, which regulates the expression of CO and FT in Arabidopsis and controls soybean flowering through regulating GmFT2a (Watanabe et al., 2011). On the other hand, a target trait such as PH or NN performs dynamically when plants grow gradually. However, the phenotypes of PH and NN were mostly investigated at the mature stage. Sun et al. (2006) reported that different QTL architectures have been found for PH at the different developmental stages through linkage mapping. Although several studies of the developmental behavior of quantitative traits have been reported in soybeans (Vodkin et al., 2004;Li W. et al., 2007;Xin et al., 2008;Teng et al., 2009), the genetic architecture of dynamic development behavior of complex traits remains to be further explored.
With the wide application of next-generation sequencing techniques, high-throughput single nucleotide polymorphism (SNPs) have been discovered and utilized to construct highresolution genetic maps and to conduct GWAS (Hyten et al., 2008;Michael and VanBuren, 2015;Song et al., 2016). GWAS is a powerful approach because it takes full advantage of all recombination events that occur in the evolutionary process of a natural population. It has been successfully used to explore the genetic basis for a broad range of complex traits in many plant 1 www.soybase.org species such as Arabidopsis (Atwell et al., 2010;Horton et al., 2012), rice Yang et al., 2014), maize (Kump et al., 2011;, and soybean (Hwang et al., 2014;Sonah et al., 2015;Zhou et al., 2015;Zhang et al., 2016;Chang and Hartman, 2017).
The mixed linear models (MLMs) have been widely used for GWAS. The compressed MLM (CMLM) was also utilized to reduce computing time and to improve statistical power for quantitative trait nucleotide (QTN) detection (Zhang et al., 2010). Nevertheless, the current GWAS methods such as MLM and CMLM are mainly based on the single-locus genome-wide scan, which often requires correction for multiple tests. The typical Bonferroni correction is so conservative that some small-effect loci may not reach the significance threshold. With the rapid development of statistical methods, several multi-locus GWAS approaches have been developed to improve the power of QTN detection (Cho et al., 2010;Segura et al., 2012;Moser et al., 2015). The obvious advantage of these methods is no Bonferroni correction due to the nature of multi-locus methods. Recently, Wang et al. (2016) proposed a new multi-locus random-SNPeffect mixed linear model (mrMLM) method to improve the power and accuracy of GWAS. Differing from the other multilocus methods, the mrMLM is a two-stage method. At the first stage, the SNP effect is viewed as being random, and all the potentially associated markers are selected by a random-SNP-effect MLM with a modified Bonferroni correction for significance test. At the second stage, all the selected markers are placed into one model and all the non-zero effects are further detected by a likelihood ratio test for QTN identification.
Summer-planting soybean is a major soybean crop grown in the region between the Yangtze River and the Huai River in the southern region of middle China, an important soybean production area. Although the genetic architecture of some agronomic traits such as PH was reported in our previous study (Cao et al., 2017) in the summer planting soybean, the genetic bases of dynamic PH and NN for them remain largely unknown. The aim of this study was to dissect the genetic basis of PH and NN at three different developmental stages in 368 summer planting soybean genotypes using the GWAS strategy. Our findings will provide useful genetic information for soybean molecular breeding.

Plant Materials, Field Trials and Phenotypic Evaluation
A soybean breeding line (SBL) population containing 368 accessions was established to service the local soybean breeding. All these pure lines were obtained from the National Center for soybean improvement, Nanjing Agricultural University, Nanjing, China. All experimental materials were planted at Jiangpu (JP) (32 • 12 N and 118 • 37 E) and Fengyang (FY) (32 • 47 N and 117 • 19 E) Station in Jiangsu and Anhui province, respectively, on June 20, 2011. At each location, the experimental design was a randomized complete block with 50 cm × 50 cm hill plots and three replications. The phenotypes for PH and NN were measured at the three different developmental stages over two locations: Stage 1 (35 days after the emergence of seedlings), Stage 2 (50 days after the emergence of seedlings) and Stage 3 (90 days after the emergence of seedlings). All the phenotypes were named PH1, PH2 and PH3 for PH, and NN1, NN2 and NN3 for NN. PH and NN were the averages of three measurements per plot.

Statistical Analysis
The analysis of variance (ANOVA) was performed for all traits using the PROC GLM procedure of SAS version 9.3 (SAS Institute, Inc., Cary, NC, United States). The model for the phenotype of a trait was y where µ is the total mean, G i is the effect of the i th genotype, E j is the effect of the j th environment, GE ij is the interaction effect between the i th genotype and the j th environment, R k(j) is the effect of the k th block within the j th environment, and e ijk is a random error following N(0, σ 2 e ). Descriptions of all traits were determined by the mean of each trait over two locations. The broad-sense heritability (h 2 ) was calculated as: h 2 = σ 2 g / σ 2 g + σ 2 ge /n + σ 2 e /nr for combined environments and h 2 = σ 2 g / σ 2 g + σ 2 e for an individual environment, where σ 2 g is the genotypic variance, σ 2 ge is the genotype by environment interaction variance, σ 2 e is the error variance, n is the number of environments, and r is the number of replications. Variance components and correlation coefficients were estimated by the PROC VARCOMP and CORR procedure of SAS, respectively. To minimize the effects of environmental variation, the best linear unbiased predictors (BLUPs) of individual lines for each trait were calculated using the R package lme4 (Bates et al., 2015).

Genotyping, SNPs Polymorphism and Haplotype Block Estimation
High-throughput SNPs were generated by RAD-seq. The quality control of sequencing data and methods of calling variations are described in our previous study . A total of 62423 SNPs with a minor allele frequency (MAF) ≥ 5% were used for further analysis in the present study.
The MAF of the SNPs was calculated using VCFtools software (Danecek et al., 2011). Haplotype blocks were estimated using pLINK V1.90 software (Purcell et al., 2007) with the command option -blocks, following the default algorithm as described by Gabriel et al. (2002). The visualization of haplotype blocks was carried out with the R package LDheatmap (Shin et al., 2006). The estimated parameters for SNPs polymorphism were displayed using circos (Krzywinski et al., 2009).

Linkage Disequilibrium Estimation
Linkage disequilibrium (LD) between pairwise SNPs was calculated as the squared correlation coefficient (r 2 ) of alleles using the linkage disequilibrium tools option of RTM-GWAS V1.1 software . The r 2 value was calculated for all pairwise SNPs with a 100 kb summary bin setting within the 5 Mb distance and then averaged across the whole genome. Because of the substantial difference in recombination rate between euchromatic and heterochromatic regions, the r 2 value was calculated separately for the two chromosomal regions. The physical length of the euchromatic and heterochromatic regions for each chromosome was defined as in the G. max 1.0 reference genome. The LD decay rate was measured as the chromosomal distance at which the average pairwise r 2 dropped to half its maximum value . Only r 2 for SNPs with pairwise distances less than 5 Mb in either the euchromatic or heterochromatic region was used to draw the average LD decay figure by R script.

Population Structure and Principal Component Analysis
Filtering SNPs used the -indep-pairwise command option of pLINK. The pruned data were then used to estimate population structure using ADMIXTURE V1.3.0 software (Alexander et al., 2009). In the ADMIXTURE setting, the number of clusters (K) was set from 1 to 10 initially; then, each Q and the relevant P-value was estimated. The most likely number of subpopulations was determined using the method described in Evanno et al. (2005). A principal component analysis (PCA) of whole-genome SNPs was performed using EIGENSOFT V5.0.2 software (Price et al., 2006) smartpca program, and the first two eigenvectors were plotted in two dimensions. The neighbor-joining tree was constructed using TASSEL V5.2 software (Bradbury et al., 2007).

Genome-Wide Association Studies
After excluding SNPs with an MAF < 5%, 62423 SNPs were retained for 368 soybean accessions. To minimize false positives and increase statistical power, the population structure (Q) and kinship (K) matrix were estimated for the population. For the MLM, both the regular MLM and compressed CMLM involve the Q and K matrices as a fixed effect and random effect, where the Q matrix was replaced by the principal components (PCs) in CMLM. CMLM was implemented by the R package GAPIT (Genome Association Prediction Tool) V2 (Tang et al., 2016). Another R package, mrMLM V2.1, representing the mrMLM method was adopted . Thus, GWAS was conducted by combining the CMLM and mrMLM methods in this study. The critical threshold of significance for SNP-trait association was set at a P-value = 1.0 × 10 −4 in CMLM according to the empirical value and at a LOD value of 3 in mrMLM. A QTN was defined as a haplotype block possessing SNPs identified as significantly associated with a trait (Schneider et al., 2016). The QTNs were named following the nomenclature described by McCouch et al. (1997). In addition, the abbreviation was used for the loci associated with the traits at the different stages. Thus, qPH(NN)(1,2,3)-10-1 indicated a locus located on chromosome 10 associated with both PH and NN at all three stages.

Prediction of Candidate Genes
Genes annotated in G. max Williams 82 reference gene model 1.0 were used as the source of candidate genes. The prediction of candidate genes mainly referred to the genes with a known function in soybeans related to the trait or the orthologs in Arabidopsis. 1 | Descriptive statistics for plant height (cm) and number of nodes on the main stem at three stages over two environments in the SBL population.

Trait Performance of the Tested Population
The phenotypic characteristics of PH and NN for the 368 soybean lines are shown in Table 1. Averaged over two environments, PH and NN showed a large variation at the three different stages with range values of 27.97-109.06 (cm) and 7.67-20.39, respectively. The absolute values of kurtosis and skewness were approximately 1 for both PH and NN (Supplementary Figure S1). Significant positive correlations were observed for PH and NN among the three stages and between PH and NN (Supplementary Figure S2) (r > 0.60, P < 0.0001). PH1 and PH3 were moderately correlated with NN3 and NN1 at r = 0.49 and 0.48, P < 0.0001, respectively. Analysis of variance (ANOVA) indicated that there were significant differences in the effects of genotypes, environments, and their interactions for the traits at all stages. Additionally, a relatively high heritability (≥60%) was estimated for PH and NN at all stages, indicating that the genetic effects play a primary role in the performance of PH and NN.

Characterization of the SNPs, Population Structure, LD, and LD Haplotype Block Estimation
A total of 62423 SNPs with an MAF ≥ 0.05 were used for further analyses, with an average marker density of 1 SNP every 16.42 kb genome-wide, varying across chromosomes from 29.39 kb per SNP on chromosome 5 to 10.51 kb per SNP on chromosome 15. The MAF and haplotype block (>50 kb) for the population characteristics are presented in Figure 1. The most likely K-value was K = 3 based on the analysis of population structure (Figure 2A), which suggested that the overall population could be divided into three subpopulations. This result was also supported by the phylogenetic analysis (Figures 2B,C) and PCA ( Figure 2D). The LD decay rate of the population was estimated at 400 kb in euchromatin, where r 2 dropped to half of its maximum value (r 2 = 0.23) (Figure 3). In heterochromatin, however, r 2 did not drop to half of its maximum value until 3.5 Mb. The haplotype analysis showed that 62423 SNPs were grouped into 5697 haplotype blocks. The size of the blocks ranged from 6 bp to 200 kb across the whole genome. The distribution of haplotype blocks is shown in Supplementary  Figure S3.

GWAS for the Traits via CMLM
Genome-wide association studies was conducted using the BLUPs of individual performance over two environments. A total of 11 loci for PH and 13 loci for NN were identified in the CMLM association panel at the suggestive significance level (P = 1 × 10 −4 ). There were 8, 3, and 4 loci for PH and 6, 6, and 8 loci for NN at the three developmental stages. Three and six loci were detected for PH and NN at more than two stages, respectively (Table 2 and Figure 4). Notably, two major loci associated with PH and NN were identified under the Bonferroni correction for multiple tests (0.05/N) (bold in Table 2). qPH (NN)(1,2,3)-10-1 was identified on chromosome 10 at all stages, while qPH(NN)3-19-1 was identified on chromosome 19 only at the last stage. The haplotype analysis showed that the peak SNPs for these two major loci were located within the two haplotype blocks named H2842 and H5441. Three main haplotypes (total frequency > 80%), H2842-1, H2842-2, and H2842-3, were identified for H2842, whose frequencies were 19.02, 58.42, and 10.87%, respectively, and two main haplotypes, H5441-1 and H5441-2, were identified for H5441, with frequencies of 46.47% and 33.70%, respectively. We further analyzed the effect of these haplotypes for PH and NN. The results showed that there were significant differences for both PH and NN among the main haplotypes of H2842 or H5441 (Supplementary Figure S4). The LD and haplotype blocks for these two major loci are presented in Figure 5.

GWAS for Traits via mrMLM
To validate the reliability of the loci determined by the CMLM method and identify more loci associated with PH and NN, a multi-locus random effect MLM method was used to conduct GWAS.

Prediction of Candidate Genes
To predict candidate genes for loci significantly associated with both PH and NN, we selected the putative genes tagged by the most significant SNPs. Two major loci coassociated with PH and NN pinpointed the known genes. E2(Glyma10g36600), encoding a homolog of GIGANTEA (GI) protein, is one of the key genes regulating soybean flowering and maturity by regulating GmFT2a (Watanabe et al., 2011). It was found 370 kb upstream of the peak SNP (Gm10_44346474, MAF = 0.41) of qPH1-10-1 (qNN1-10-1) on Gm10, which was associated with both PH and NN at all three stages (Tables 3, 4). Dt1 (Glyma19g37890) is a homolog of Arabidopsis TERMINAL FLOWER1 (TFL1) and plays a primary role in the soybean stem growth habit (Bernard, 1972). It was found 41 kb upstream of the peak SNP (Gm19_44938780, MAF = 0.35) of qPH3-19-1 (qNN3-19-3) on Gm19, which was associated with both PH and NN at the last stage. In addition to these two major loci, we predicted the candidate genes for two other loci associated with both PH and NN. A putative gene, Glyma04g40640, was found 35 kb from the peak SNP (Gm04_46671367, MAF = 0.11) of qPH1-4-1 (qNN1-4-1) on Gm04, which was associated with both PH and NN at the first stage. It was homologous to Arabidopsis APRR5, which is involved in various circadianassociated biological events such as flowering time in long-day photoperiod conditions and red light sensitivity of seedlings during early photomorphogenesis (Kamioka et al., 2016). The putative gene Glyma13g06811 was identified 125 kb from the peak SNP (Gm13_6859748, MAF = 0.46) of qPH2-13-1(qNN3-13-2) on Gm13, which was associated with both PH and NN at the second and third stages, respectively. It is homologous to Arabidopsis AGAMOUS-LIKE 8 (AGL8), which is a MADs-box transcription factor involved in various biological events such as flower and fruit development and the maintenance of inflorescence meristem characteristics (Ma et al., 1991).

Linkage Disequilibrium and the Statistical Method Are of Great Significance in GWAS
Genome-wide association studies are a powerful tool to elucidate the genetic architecture for complex quantitative traits in crops (Morris et al., 2013;Su et al., 2016;Sun et al., 2017). The mapping resolution and statistical power are the main considerations in GWAS. LD is one of the factors limiting the mapping resolution of GWAS. The recombination rate is one of the major factors affecting LD extension and is different between euchromatic and heterochromatic regions . Previous studies suggested that there was a very diverse range of LD values in different crops and different chromosomal regions in a specific crop (Li Y.H. et al., 2014;Sonah et al., 2015;Zhang et al., 2015). Similar to previous studies, a large difference in LD decay rate was also observed between the euchromatin (400 kb) and heterochromatin regions (3.5 Mb) of soybeans in this study. A longer LD has been observed in self-pollinated crops such as soybean compared to cross-pollinated crops such as maize . A previous study also reported a longer LD for the chromosomes involved in the domestication process. The QTLs for domestication traits such as seed weight and flowering were mapped in these regions (Li Y.H. et al., 2013).
Another concern with GWAS is the statistical method. In the present study, the CMLM method was used for the singlelocus GWAS. Consistent with the previous study, we also found that the CMLM method could take less computing time than regular MLM and reduce false positive results simultaneously . However, the single-locus GWAS methods often need correction for multiple tests. For instance, the typical Bonferroni correction corrects an α = 0.05 to α = 0.05/m, where m is the number of statistical tests performed. For a GWAS with 500,000 markers, the statistical significance threshold for an association would be corrected to 1e −7 , such that no or a few loci could reach the significance threshold  after the correction. Such a situation is not always suited to the nature of complex traits. A previous study showed that no significantly associated locus for soybean seed weight was detected, possibly for this reason (Fang et al., 2017). Wang et al. (2016) also reported that some small-effect loci were not significantly associated with the traits in the single-locus approach under the Bonferroni correction but significantly associated with that in the mrMLM method. Actually, the smalleffect loci should not be neglected in the genetic system of complex traits. Fortunately, several multi-locus GWAS methods have already been reported in previous studies (Segura et al., 2012;Wang et al., 2016), where no Bonferroni correction is needed because of the multi-locus nature. Thus, the multilocus GWAS method may play an increasingly important role in dissecting the genetic architecture of complex traits in the post-association time.

Combination of CMLM and mrMLM GWAS Methods to Identify the Major and Minor Loci for PH and NN in the SBL Population
Plant height and NN are quantitative traits controlled by numerous loci in soybeans 2 . In the present study, 11 loci for PH and 13 loci for NN were detected via the CMLM method at the suggestive threshold level, only two of which were identified after the Bonferroni correction. To confirm the loci determined by the CMLM method and identify additional loci for PH and NN, mrMLM as a multi-locus method was used for GWAS in this study. As expected, we found 32 and 28 loci, except for two major loci, for PH and NN via the mrMLM method. Some novel loci in comparison with that obtained from the CMLM were located around the previously reported QTLs. For PH, the loci qPH1-18-1 and qPH1-18-2 detected by the mrMLM method have been previously identified by Sun et al. (2006), which were also located in the same genomic regions of the known QTLs. However, only six loci in addition to the two major loci, qPH1-4-1, qPH1-6-1, and qPH1-14-1 for PH and qNN1-4-1, qNN1-14-1, and qNN2-13-1 for NN, identified at a suggestive threshold level (P = 10 −4 ) in CMLM were validated in mrMLM association panel. One potential reason was that the mrMLM method improved the power and accuracy for QTN detection due to the nature of the statistical model. Thus, many novel loci were detected by the mrMLM method. Another possible reason was that the relatively stringent threshold still applied to the CMLM method. More loci might be commonly detected in two association panels if a lower threshold was adopted in the CMLM method, but the false positive results might increase under such conditions. Undoubtedly, mrMLM can identify not only the major loci but also the minor loci for quantitative traits compared with the CMLM method.

Genetic System of Dynamic PH and NN in Summer Planting Soybeans
Many agronomic and yield-related traits such as PH and NN were highly correlated in soybeans. Previously identified loci for these traits usually co-located in the same genomic regions (Malik et al., 2007;Liu et al., 2011Liu et al., , 2017Fang et al., 2017). Similarly, significant positive correlations were observed between PH and NN at different stages in this study. Furthermore, six loci, including major and minor loci, were shared for PH and NN, suggesting that both PH and NN have a similar genetic system controlled by major and minor loci. Previous studies showed that the E2 and Dt1 loci have an effect on many agronomic and yield-related traits in soybean (Liu et al., 2011;Kato et al., 2015;Zhang et al., 2015). We also found that PH and NN shared these two loci and further confirmed that the genetic pattern of the E2 locus was different from that of the Dt1 locus. The former was detected for both PH and NN at all stages while the latter was only detected at the last stage. Plant height and NN in soybeans are dynamic traits, as the phenotype changes constantly during the plant lifecycle. However, studies of the genetic basis of PH and NN have mainly measured the final PH and number of nodes on the main stem, especially when the phenotype was investigated at the mature stage. We accessed PH and NN in the summerplanting accessions at three different stages, which included one vegetative stage and two reproductive stages when plants were in late vegetative growth, as well as the flowering and mature stages, to reveal the genetic control underlying the dynamic PH and NN based on the GWAS strategy in this study. A previous report showed that the haplotype was more appropriate than the single SNP to uncover the genetic variation and improve the efficiency of breeding for target traits due to the existence of multiple alleles (Hao et al., 2012). We identified the E2 and Dt1 genes in two haplotype blocks. The analysis of haplotypes revealed that the main haplotypes of these two haplotype blocks were related to PH and NN over stages. Our results suggested that E2 and Dt1, as the major loci, play different roles in regulating the development of PH and NN at different stages. Nine and 20 novel loci were identified for PH and NN at different stages via a new multi-locus GWAS method, respectively. Moreover, the differential loci were identified for both PH and NN among the different stages. These common and specific loci for PH and NN at different stages unveil the genetic architecture underlying the dynamic PH and NN. Although most of them have not been confirmed yet, candidate genes were predicted for several loci  Table. associated with both PH and NN, some of which are located in the same regions as known QTLs for PH and NN in SoyBase. Further studies should confirm these loci and identify candidate genes for them.

CONCLUSION
More loci, including 34 loci for PH and 30 loci for NN, were identified by the mrMLM method than by the single-locus CMLM method. A few loci were commonly identified for PH and NN via the two methods at the different developmental stages. One stable locus that overlapped with the E2 gene was identified for PH and NN at all three stages, while another major locus, referred to as the Dt1 gene, was determined at the last stage by both methods. Most loci were mainly detected at only one or two of the examined developmental stages. The dynamic PH and NN was controlled by a set of specific loci and a few common loci in summer planting soybeans.  Table. FIGURE 6 | Venn diagrams for loci associated with PH (A) and NN (B) over three stages.
Frontiers in Plant Science | www.frontiersin.org