Introduction {#s1}
============

Soybean \[Glycine max (L.) Merr.\] is an important crop around the world, and worldwide soybean consumption has increased rapidly in recent years. Increasing production per unit area thus remains an important breeding target. Yield is influenced by a variety of traits, such as number of pods, number of seeds, and number of nodes. Plant height is an important factor in the formation of yield traits, and at the same time promotes or inhibits other yield components. Plant height mainly affects yield by influencing lodging resistance and the number of pods per plant ([@B1]). With the development of molecular genetics and molecular marker technology ([@B28]), the traditional breeding process has been greatly improved, producing cultivars with excellent traits.

Most agronomic traits in soybean, such as plant height and pod number, are controlled by quantitative trait loci (QTLs) ([@B42]), making it difficult to understand their genetic basis and molecular mechanism due to insufficient precision. A large number of plant height-related QTLs have been reported. [@B52] used a recombinant inbred line (RIL) population containing 236 individuals derived from a cross between "Minsoy" and "Noir 1" to identify nine soybean plant height QTLs distributed on eight chromosomes. [@B54] used a RIL population containing 143 individuals derived from a cross between "Charleston" and "Dongnong 59" to uncover 17 plant height QTLs distributed on 10 chromosomes. [@B14] identified three soybean plant height QTLs using three backcross-derived populations. [@B61] used an F2-derivative group constructed from "Jiyu50" and "Jnong18" and containing 236 individuals to identify nine plant height QTLs distributed across linkage groups A1, C2, G, M, F, and K.

Genome-wide association studies (GWAS) exploit the abundant mutations and recombinations existing in diverse populations. They have a relatively high localization accuracy compared to traditional QTL mapping methods, enabling prediction and identification of genes. With the development of genotyping and sequencing technologies, GWAS have been widely used to study the genetic basis of traits in major crops such as maize ([@B47]; [@B39]; [@B11]), rice ([@B20]; [@B2]), wheat ([@B21]; [@B53]), *Arabidopsis thaliana* ([@B31]), and tomato ([@B65]). GWAS are also widely used for genetic analysis of major traits in soybean such as flowering time ([@B64]), pod number ([@B10]), seed oil ([@B70]), seed protein ([@B3]), salt tolerance ([@B27]), and 100-seed weight ([@B63]). However, there are few reports on GWAS of high-density single-nucleotide polymorphisms (SNPs) in soybean ([@B22]). Therefore, the application of GWAS in soybean, especially for agronomic traits, remains to be explored. Previous studies have shown that soybean plant height is a typical quantitative trait controlled by multiple genes, and its genetic effect is relatively small. However, more effective methods are required for detecting QTLs. Multi-locus GWAS is a suitable method for detecting significant quantitative trait nucleotides (QTNs) and has been used in several studies ([@B17]; [@B69]).

In addition to natural populations, GWAS are widely used for genetic analysis of other mapping populations, such as nested association mapping (NAM) and multi-parent advanced generation inter-cross (MAGIC) populations ([@B30]). [@B57] used a maize NAM population for GWAS to determine the genetic basis of important leaf structural features. [@B8] constructed a NAM population of 25 RILs and identified beneficial alleles for improving maize kernel starch, protein, and oil. [@B67] used a four-way RIL to perform GWAS on the genetic mechanism of soybean protein. [@B16] conducted GWAS using a hop population constructed from the parents "Teamaker" and USDA 21422M, revealing the genetic basis of fiber quality traits and yield components in 231 lines. [@B41] used a F~6:8~ RIL derived from cotton cultivars "Lumianyan28" and "Xinluzao24" to detect 134 QTLs for fiber quality traits and 122 QTLs for yield components, with 35 common candidate genes. Therefore, analyzing the genetic basis of complex traits in single bi-parental populations is feasible using GWAS. At the same time, [@B46] consider that linkage analysis and association analysis together are more accurate and effective than single analysis methods in revealing the genetic mechanisms of traits.

To investigate the genetic basis of soybean plant height variation, we chose two parents with significant genetic differences in plant height, "Dongnong L13" and "Henong 60", to construct a RIL population with 156 lines. We used simple-sequence repeat (SSR) markers for genotyping and obtained QTLs related to plant height from linkage analysis between genotypic data and phenotypic data obtained in different environments. We then performed SNP genotyping of the RIL population using gene chip technology and conducted association analysis between SNP marker genotypes and phenotypic traits. Association analysis results not only further verified the QTLs obtained by linkage analysis, but also narrowed the search for candidate genes related to plant height, which will accelerate molecular breeding to select and improve agronomic traits in soybean.

Materials and Methods {#s2}
=====================

Plant Materials {#s2_1}
---------------

Two soybean cultivars, "Dongnong L13" (P1) and "Henong 60" (P2), with wide genetic differences in plant height were used as parents for constructing RILs. "Dongnong L13" was bred by crossing "Heinong 40" and "Jiujiao 5640", and "Henong 60" was bred by crossing "Beifeng 11" and "Hobbit". The cross P1 × P2 was carried out in Harbin (E 126.63°, N 45.75°), Heilongjiang Province, China, and F~1~ seeds were harvested in 2008. In the same winter, F~1~ seeds were sowed in Yacheng City, Hainan Province (E 109.00°, N 17.50°). From 2010 to 2014, F~1~ plants were self-crossed in Harbin and Yacheng following the single-seed descent method, that is, single seeds were selected from individual plants at each generation until all individuals showed homozygous genotypes. RIL6013 containing 156 families was obtained and used for construction of genetic maps and QTL mapping.

Field Experiment and Trait Evaluations {#s2_2}
--------------------------------------

Experiments were carried out at four locations: Keshan (E 125.64°, N 48.25°), Harbin, Acheng (E 127.63°, N 45.82°), and Shuangcheng Cities (E 126.92°, N 45.75°). Field tests of nine environments across 4 years and four locations were carried out: Keshan (E1) in 2013, Harbin (E2) in 2014, Keshan (E3) and Harbin (E4) in 2015, Acheng (E5), Shuangcheng (E6), and Harbin (E7) in 2016, Shuangcheng (E8) in 2017, and Acheng (E9) in 2018. All plant materials in each environment were grown in three-row 3 m × 0.7 m plots in a completely randomized block design with three replications. Experimental plots were managed identically to local soybean crops.

Plant height was investigated in the field after maturity and defined as length from the cotyledon mark to the top of the main stem. Plant height of the 156 RILs was measured with a meter ruler. Plants with marginal utility were avoided during the survey, and 10 plants were randomly selected from each plot to conduct field investigations to determine plant height. The mean value of the 10 observations was taken as the observation value for the plot, and the average of the observation values of the three block replicates was used as the phenotypic data for QTL mapping.

Statistical Analyses of Phenotypic Data {#s2_3}
---------------------------------------

Minimum, maximum, mean, standard deviation, kurtosis, and skewness of the sample observations were calculated using the software SAS V9.21 (<https://www.sas.com/en_us/home.html>) and histograms of plant height frequency distribution of the RIL6013 population in the nine environments were drawn. The significance of the genotypic variance in each group of observations was then calculated using the general linear model (GLM) for analysis of variance (ANOVA). Finally, the VARCOMP command was used to estimate the variance component of the mixed linear model (MLM), and the generalized heritability of plant height in a single environment was calculated according to the following equation:

h

2

=

σ

G

2

σ

G

2

+σ

ϵ

2

(where $\sigma_{\text{G}}^{2}$ is the genotypic variance, $\sigma_{\text{ϵ}}^{2}$ is the error variance).

SSR Marker Map and QTL Analysis {#s2_4}
-------------------------------

Construction of the RIL6013 map was based on the method of [@B45]. The genetic map contained 20 soybean linkage groups and 137 SSR markers covering a total genome length of 1886.80 cM. The genetic distance of each linkage group was 19.68 cM (H) -- 163.67 cM (F), and the average genetic length was 94.34 cM. Each linkage group contained four to 11 markers, and the average genetic distance between the two markers was 16.13 cM.

Interval mapping (IM-ADD) and inclusive composite interval mapping (ICIM-ADD) under the .bip (QTL mapping in bi-parental populations) function built into the software QTL Icimapping V 4.1 ([@B59]) were used to detect additive QTLs. The Scan Step was set to 1.00 cM and the LOD threshold was set to 2.50. In addition, for ICIM-ADD, the PIN value was set to 0.001. After the QTL positioning results were obtained, they were named according to the method of [@B43].

SNP Genotyping {#s2_5}
--------------

DNA samples extracted were used for SNP genotyping using the SoySNP660K BeadChip at Beijing Boao Biotechnology Co., Ltd. A total of 54836 SNPs across 20 chromosomes remained after quality filtering; SNP markers were excluded by minor allele frequency (MAF \< 0.05), and the maximum missing sites per SNP was \<10% ([@B4]). These SNPs were used for analysis of population structure and GWAS.

Analysis of Population Structure {#s2_6}
--------------------------------

The analysis of population structure was performed using STRUCTURE V2.3.4 ([@B48]). For each run, the number of burn-in iterations was set to 100000, and the number of Marko Chain Monte Carlo (MCMC) was set to 100,000, while the mixed model and allele frequency correlation model were considered in the analysis. Set the K number of the subpopulations in the population from 1 to 10, and number of iterations was set to 5.To determine the best K value using STRUCTURE HARVESTER ([@B9]) (<http://taylor0.biology.ucla.edu/structureHarvester>).

Linkage Disequilibrium Analysis {#s2_7}
-------------------------------

TASSEL 5.0 ([@B6]) was used to analyze linkage disequilibrium (LD) by analyzing r^2^ values of all pairs. The LD decay trend was analyzed using negative natural logarithm, and the physical distance of LD decay was estimated to where r^2^ dropped to 0.2.

Genome-Wide Association Studies {#s2_8}
-------------------------------

mrMLM.GUI software was used to perform GWAS with the following the five methods for identifying significant QTNs: mrMLM ([@B58]), FASTmrMLM ([@B56]), FASTmrEMMA ([@B60]), pLARmEB ([@B66]), and ISIS EM-BLASSO ([@B56]). The critical P value of FASTmrEMMA was set to 0.005 while the critical P value parameter of the other methods was set to 0.01, and the critical LOD value of significant QTNs was set to 3 ([@B58]). QTNs located in at least two environments or detected using two different methods were considered to be significant. The kinship matrix used in the analysis process was also calculated by the software.

Identification of Potential Candidate Genes {#s2_9}
-------------------------------------------

QTNs were used to predict potential candidate genes based on GWAS. Genes highly expressed in stems according to the Phytozome website (<https://phytozome.jgi.doe.gov>) were searched for in the interval of LD decay distance when r^2^ dropped to 0.2 on both sides of the QTN. These genes were used for pathway analysis by combining gene annotation information and protein conserved domain functions from the NCBI database (<https://www.ncbi.nlm.nih.gov/>) and previous studies. Potential candidate genes were identified from the pathway analysis and used to identify homologous genes on the ensemble plant website (<http://plants.ensembl.org/index.html>) and speculate on their potential functions based on their gene ontology (GO) number (<https://www.ebi.ac.uk/QuickGO/>).

Results and Analysis {#s3}
====================

Phenotypic Variation Analysis {#s3_1}
-----------------------------

The RIL6013 population showed a bimodal distribution in the nine environments, except for E8 and E6, and plant height observations generally showed a unimodal distribution. We speculate that the plant height of this population may be regulated by multiple major and minor genes ([**Figure 1**](#f1){ref-type="fig"}). The plant height range of the population was generally higher than that of the parents, indicating a large separation within the population. Genotypic variances between lines were extremely significant (P \< 0.01), and the kurtosis and skewness of statistical observations ranged from \[−1, 1\]. We therefore considered that this population was suitable for plant height research and QTL identification ([**Table 1**](#T1){ref-type="table"}).

![Frequency histograms of plant height in RIL6013 in population in nine environments.](fpls-11-00009-g001){#f1}

###### 

Statistical analysis of soybean plant height for parents and the RIL6013 population in nine environments.

  Env.[^a^](#fnT1_1){ref-type="table-fn"}   Parents[^b^](#fnT1_2){ref-type="table-fn"}   RIL6013 population                                                                          
  ----------------------------------------- -------------------------------------------- -------------------- ------- -------- -------- ---------------- ------- ------- ----------- -------
  E1                                        --                                           --                   77.00   144.00   67.00    117.71 ± 15.12   −0.61   −0.19   9.69\*\*    74.30
  E2                                        --                                           --                   67.40   151.40   84.00    120.57 ± 18.55   −0.72   −0.21   13.12\*\*   80.05
  E3                                        131.00                                       102.00               72.80   145.80   73.00    116.98 ± 16.32   −0.78   −0.04   11.75\*\*   78.07
  E4                                        \-                                           \-                   72.00   157.00   85.00    120.83 ± 19.14   −0.52   −0.46   13.90\*\*   81.02
  E5                                        114.00                                       88.67                59.00   135.33   76.33    104.69 ± 17.87   −0.76   −0.21   12.76\*\*   79.67
  E6                                        135.67                                       109.84               65.33   149.00   83.67    122.41 ± 18.12   −0.67   −0.42   12.80\*\*   79.82
  E7                                        137.50                                       107.00               63.50   152.00   88.50    122.88 ± 17.77   −0.95   0.54    11.47\*\*   77.59
  E8                                        91.33                                        84.00                66.50   172.50   106.00   118.57 ± 24.73   0.08    −0.74   23.81\*\*   88.61
  E9                                        116.67                                       95.17                51.33   156.00   104.67   116.08 ± 23.21   −0.58   −0.10   27.30\*\*   89.78

Env, environment: E1, Keshan City, 2013; E2, Harbin City, 2014; E3, Keshan City, 2015; E4, Harbin City, 2015; E5, Acheng City, 2016; E6, Shuangcheng City, 2016; E7, Harbin City, 2016; E8, Shuangcheng City, 2017; E9, Acheng City, 2018.

Parents: P~1~, female cultivar "Dongnong L13"; P~2~, male cultivar "Henong 60".

Min, minimum observed value among samples.

Max, maximum observed value among samples.

Range, difference between maximum and minimum value.

Mean ± standard deviation of the observed values.

F-value for genotypic variation; \*\*extremely significant difference at level of P \< 0.01.

h^2^, Broad-sense heritability of plant height in a single environment.

The coefficient of variation (CV) for plant height was between 12.84% and 20.86%. The range of plant heights in RIL6013 demonstrated that plant height variation is large. Mean heights of the population under E4, E6, and E7 were higher than those under the other six environments, and the CVs in E8 and E9 were larger than those in the other seven environments ([**Table 1**](#T1){ref-type="table"}). This suggested that different QTLs might be detected in different environments.

From a breeding point of view, except in E8, the plant height of the parents fell within the range of \[Mean−STD, Mean+STD\] for the population, indicating that the population has strong transgressive heterosis and plant height breeding goals will have a larger choice within the population.

The heritability of RIL6013 was between 74.30% and 89.78% in a single environment ([**Table 1**](#T1){ref-type="table"}). This high heritability indicates that the genetic variance is superior to the error variance, and RIL6013 is suitable for high-quality selection of plant height in soybean.

Correlation analysis revealed a strong correlation between the various environments ([**Table 2**](#T2){ref-type="table"}). The ANOVA results showed extremely significant differences (P\<0.01) between the genotypes of this population, and also significant differences between the various environments. There were significant differences between genotypes, environment, and the genotype-by-environment interaction ([**Table 3**](#T3){ref-type="table"}), indicating that soybean plant height is highly influenced by all of these factors.

###### 

Correlation analysis of nine environments.

       E1             E2       E3       E4       E5            E6            E7       E8            E9
  ---- -------------- -------- -------- -------- ------------- ------------- -------- ------------- ----
  E1   1                                                                                            
  E2   0.731^\*\*^    1                                                                             
  E3   0.055          0.048    1                                                                    
  E4   0.062          0.044    0.091    1                                                           
  E5   −0.053         -0.046   0.002    −0.040   1                                                  
  E6   0.012          −0.007   0.023    −0.046   0.681^\*\*^   1                                    
  E7   −0.066         −0.105   0.035    0.015    0.616^\*\*^   0.669^\*\*^   1                      
  E8   −0.231^\*\*^   −0.112   −0.115   −0.060   0.070         −0.053        0.005    1             
  E9   −0.217^\*\*^   −0.121   −0.093   −0.089   −0.147        −0.114        −0.125   0.688^\*\*^   1

\*\*Extremely significant difference at level of P \< 0.01.

###### 

Joint ANOVA of plant height in RIL6013 in multiple environments.

                          DF     SS           MS         F              Pr \> F
  ----------------------- ------ ------------ ---------- -------------- ----------
  Environment             8      106218.84    13277.36   170.54^\*\*^   \<0.0001
  Genotype                155    240233.10    1549.89    19.91^\*\*^    \<0.0001
  Environment\*Genotype   1219   1343337.47   1102.00    14.15^\*\*^    \<0.0001
  Error                   2765   215274.53    77.86                     

\*\*Extremely significant difference at level of P \< 0.01.

QTL Mapping of the Plant Height {#s3_2}
-------------------------------

We detected 48 plant height QTLs ([**Table S1**](#ST1){ref-type="supplementary-material"}), distributed across the 20 soybean linkage groups. One (I) to five (M) QTLs mapped to each linkage group with LOD values between 2.54 and 13.46, and each QTL explained between 0.55% and 19.55% of the phenotypic variation. Four QTLs (*qPH-C1-1*, *qPH-M-1*, *qPH-F-1*, *qPH-L-1*) explaining \>10% phenotypic variation can be considered the major QTLs controlling plant height.

Twelve QTLs were detected by both IM and ICIM methods and were distributed on 10 linkage groups A2, B1, C1, C2, D1a, D1b, F, H, L, and N ([**Table S2**](#ST2){ref-type="supplementary-material"}). LOD values were between 2.54 and 13.46, and the phenotypic variation explained (PVE) ranged from 1.45% to 19.55%. The additive effects of *qPH-F-1* detected in the E8 environment were negative, indicating that the allele derived from the male parent "Henong 60" increased plant height. The additive effects of the remaining QTLs were positive, indicating that the allele derived from the female parent "Dongnong L13" increase plant height.

We identified 21 QTLs in two or more environments, distributed over 14 linkage groups except D1b, A1, C2, B2, E, and I ([**Table S3**](#ST3){ref-type="supplementary-material"}). LOD values were 2.58 to 13.46, and the PVE ranged from 0.55% to 14.83%. Additive effects of *qPH-K-2* detected in E1 and E2, *qPH-F-1* detected in E8, and *qPH-J-1* detected in E1 and E2 were negative, indicating the genes causing greatest effect were from the male parent; for the remaining QTLs, the genes making the greatest contribution were from the female parent.

We detected nine QTLs across multiple environments and using multiple methods ([**Table 4**](#T4){ref-type="table"}): E2 and E7, one QTL (*qPH-A2-1*); E1 and E3, one QTL (*qPH-N-2*); E1 and E2, five QTLs (*qPH-B1-1*, *qPH-C1-1*, *qPH-D1a-3*, *qPH-H-1*, and *qPH-L-1*); E1, E2, E5, and E7, one QTL (*qPH-D1a-2*); E2, E5, E6, E7, and E8, one QTL (*qPH-F-1*) ([**Figure 2**](#f2){ref-type="fig"}). *qPH-A2-1* was detected in the A2 linkage group, with LOD values from 2.97 to 3.13, PVE ranging from 1.13% to 6.93%, and a positive additive effect indicating that the dominant gene was derived from the female parent. *qPH-N-2* was detected in the N linkage group with LOD value of 2.92 to 3.14 and PVE ranging from 1.12% to 1.75%, in which the gene with greatest contribution to the QTL under E1 was derived from the male parent while that with the greatest contribution under E3 was derived from the female parent. *qPH-B1-1* was detected in the B1 linkage group, with LOD value of 4.86 to 7.27, PVE ranging from 2.12% to 9.79%, and the dominant gene derived from the female parent. *qPH-C1-1* was detected in the C1 linkage group, with LOD value of 3.79 to 13.46, PVE ranging from 1.91% to 19.55%, and the dominant gene derived from the female parent. *qPH-D1a-3* was detected in the D1a linkage group, with LOD value of 2.64 to 8.35, PVE ranging from 1.60% to 9.38%, and the dominant gene derived from the female parent. *qPH-H-1* was detected in the H linkage group, with LOD value of 4.84 to 7.00, PVE from 1.02% to 7.04%, and the dominant gene derived from the female parent. *qPH-L-1* was detected in the L linkage group, with LOD value of 4.58 to 8.20, PVE ranging from 3.02% to 14.00%, and the dominant gene derived from the female parent. *qPH-D1a-2* was detected in the D1a linkage group, with LOD value of 2.83 to 5.27, PVE from 2.52% to 5.71%, and the dominant gene derived from the female parent. *qPH-F-1* was detected in the F linkage group, with LOD value of 3.29 to 6.97, and PVE ranging from 1.14% to 10.72%; the additive effect in E2 and E8 was negative, indicating that the dominant gene was derived from the male parent.

![Total numbers of QTLs detected across multiple environments and by multiple methods. Pink bars indicated the number of QTLs detected in each environment.](fpls-11-00009-g002){#f2}

###### 

Stably expressed QTLs identified by multiple methods and in multiple environments.

  QTL           Environment                  Chromosome   Left marker   Right marker   LOD[^a^](#fnT4_1){ref-type="table-fn"}         PVE(%)[^b^](#fnT4_2){ref-type="table-fn"}       Add[^c^](#fnT4_3){ref-type="table-fn"}                  Method
  ------------- ---------------------------- ------------ ------------- -------------- ---------------------------------------------- ----------------------------------------------- ------------------------------------------------------- ------------------------------------
  *qPH-A2-1*    E2,E7,E7                     8            Satt119       Sat_406        2.97,3.13,3.13                                 1.13,4.31,6.93                                  13.12,16.38,16.38                                       IM,IM,ICIM
  *qPH-B1-1*    E1,E2,E2                     11           Satt197       Sat_123        5.15,7.27,4.86                                 3.99,2.12,9.79                                  12.07,15.44,10.64                                       IM,IM,ICIM
  *qPH-C1-1*    E1,E2,E2                     4            Satt713       Satt565        3.79,13.46,10.84                               1.91,2.74,19.55                                 6.95,15.47,12.80                                        IM,IM,ICIM
  *qPH-D1a-2*   E1,E2,E5,E7,E5               1            Sat_413       Sat_160        2.83,5.27,4.75,3.65,4.75                       4.54,2.52,2.64,5.71,2.64                        12.18,16.28,17.27,17.56,17.27                           IM,IM,IM,IM,ICIM
  *qPH-D1a-3*   E1,E2,E1,E2                  1            Satt515       AZ302047       7.91,8.35,4.90,2.64                            3.18,1.60,9.38,3.22                             9.25,11.92,6.33,5.14                                    IM,IM,ICIM,ICIM
  *qPH-F-1*     E2,E5,E6,E7,E8,E5,E6,E7,E8   13           Sat_039       Sat_417        6.97,3.29,5.16,4.31,3.33,3.29,5.16,3.58,3.33   1.14,3.02,2.28,5.22,3.49,3.02,2.28,7.69,10.72   −9.55,15.15,16.14,14.20,−9.59,15.15,16.14,13.27,−9.59   IM,IM,IM,IM,IM,ICIM,ICIM,ICIM,ICIM
  *qPH-H-1*     E1,E2,E1                     12           Satt181       Satt434        7.00,6.50,4.84                                 2.33,1.02,7.04                                  8.41,10.05,5.82                                         IM,IM,ICIM
  *qPH-L-1*     E1,E1,E2                     19           Sat_195       Satt448        4.67,4.58,8.20                                 3.02,14.83,14.00                                13.19,11.54,16.03                                       IM,ICIM,ICIM
  *qPH-N-2*     E1,E3,E3                     3            Satt584       Satt660        3.14,2.92,2.92                                 1.12,1.45,1.75                                  −8.28,9.67,9.67                                         IM,IM,ICIM

LOD, logarithm of odds.

PVE, phenotypic variation explained by QTL.

Add, contribution of parents to the additive effect.

Population Structure {#s3_3}
--------------------

Selected 5000 of the 54836 with more polymorphic SNP markers ([@B49]), which were randomly distributed on 20 soybean chromosomes. Calculating ΔK using STRUCTURE 2.3.4 ([**Figure 3A**](#f3){ref-type="fig"}, K = 1--10), revealing two subgroups (K = 2) based on the ΔK value ([**Figure 3B**](#f3){ref-type="fig"}). These two subgroups contain 64 (46.04%) and 75 (53.96%) lines, respectively.

![Population structure based on 5,000 SNPs distributed across 20 chromosomes. **(A)** Plot of ΔK calculated for K = 1--10. **(B)** Population structure (K = 2); the areas of the two colors (green and red) illustrate the proportion of each subgroup.](fpls-11-00009-g003){#f3}

Linkage Disequilibrium Analysis {#s3_4}
-------------------------------

Because we used a bi-parental population, the LD decay distance was much longer than that in the natural population. As shown in [**Figure S1**](#SF1){ref-type="supplementary-material"}, r^2^ decreased gradually with increased distance, and the LD decay distance was estimated at 950 kb when r^2^ was 0.2. As this distance is too large, we determined the range for potential candidate genes according to the region showing the fastest LD decay. LD decayed fastest within 200 kb, and then gradually slowed down, so we searched for potential candidate genes in the interval of 100 kb on either side of each QTN.

QTN Location by Multi-Locus GWAS Methods {#s3_5}
----------------------------------------

We used mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, and ISIS EM-BLASSO to detect 10 QTNs ([**Table 5**](#T5){ref-type="table"}), and the five methods detected three, one, one, two, and seven QTNs, respectively ([**Figure 4A**](#f4){ref-type="fig"}). We detected two, one, one, two, one, and three QTNs in environments E1, E3, E5, E6, E8, and E9, respectively ([**Figure 4B**](#f4){ref-type="fig"}). No significant QTN was identified in environments E2, E4, or E7.

###### 

QTNs identified by multiple methods.

  Environment   Method[^a^](#fnT5_1){ref-type="table-fn"}   Marker         Chromosome   Position (bp)   QTN effect           LOD score        r2 (%)[^b^](#fnT5_2){ref-type="table-fn"}
  ------------- ------------------------------------------- -------------- ------------ --------------- -------------------- ---------------- -------------------------------------------
  E6            5                                           AX-157411267   2            14266435        −9.02                3.87             9.76
  E3            1                                           AX-157176525   4            10948200        6.82                 3.25             16.47
  E6            5                                           AX-157526374   4            2455146         5.94                 3.76             8.86
  E1            4,5                                         AX-157388309   7            35573836        6.56,6.56            4.93,5.77        14.25,16.42
  E9            5                                           AX-157136514   8            46855033        −6.50                3.38             6.11
  E9            5                                           AX-157144941   9            38573690        −8.21                4.19             10.65
  E5            1,2,3                                       AX-157278476   13           11184314        −8.61,−7.25,−14.48   4.49,3.46,4.44   16.93,12.08,11.98
  E9            1                                           AX-157296110   13           42895901        10.15                3.55             16.30
  E8            5                                           AX-157143002   16           36971410        9.90                 3.42             9.44
  E1            4,5                                         AX-117466184   20           46210693        7.04,7.04            3.68,4.05        9.16,10.55

mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, and ISIS EM-BLASSO are indicated by 1 to 5, respectively.

r2 (%), proportion of total phenotypic variation explained by each QTN.

![**(A)** Total number of significant QTNs detected by five methods. **(B)** Total number of significant QTNs detected in nine environments using five methods.](fpls-11-00009-g004){#f4}

We detected three QTNs using multiple methods ([**Table 6**](#T6){ref-type="table"}, [**Figure 4B**](#f4){ref-type="fig"}). The pLARmEB and ISIS EM-BLASSO methods detected a single QTN (AX-157388309) on chromosome 7, with LOD values ranging from 4.93 to 5.77 and PVE of 14.25% to 16.42%, and a QTN (AX-117466184) located on chromosome 20, with LOD value of 3.68 to 4.05 and PVE ranging from 9.16% to 10.55%. Both QTN effects were positive. The mrMLM, FASTmrMLM, and FASTmrEMMA methods detected a single QTN (AX-157278476) located on chromosome 13, with LOD value ranging from 3.46 to 4.49 and PVE ranging from 11.98% to 16.93%; the QTN effect was negative. The ISIS EM-BLASSO method detected the largest number of QTNs located in multiple environments, and the QTN effects detected by various methods were consistent (both positive or negative).

###### 

Stably expressed QTNs identified by multiple methods.

  Environment   Method[^a^](#fnT6_1){ref-type="table-fn"}   Marker         Chromosome   Position (bp)   QTN effect           LOD score        r2 (%)[^b^](#fnT6_2){ref-type="table-fn"}
  ------------- ------------------------------------------- -------------- ------------ --------------- -------------------- ---------------- -------------------------------------------
  E1            4,5                                         AX-157388309   7            35573836        6.56,6.56            4.93,5.77        14.25,16.42
  E5            1,2,3                                       AX-157278476   13           11184314        −8.61,−7.25,−14.48   4.49,3.46,4.44   16.93,12.08,11.98
  E1            4,5                                         AX-117466184   20           46210693        7.04,7.04            3.68,4.05        9.16,10.55

mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, and ISIS EM-BLASSO were indicated by 1 to 5, respectively.

r2 (%), proportion of total phenotypic variation explained by each QTN.

Potential Candidate Gene Determination {#s3_6}
--------------------------------------

We searched for potential candidate genes within 100-kb intervals on both sides of four QTNs, AX-157278476, AX-157296110, AX-157176525, and AX-157526374, which were repeatedly identified by linkage analysis and GWAS. Five genes were in the interval 11.08--12.08 Mb and 25 genes in the interval 42.80--43.00 Mb on chromosome 13, and three genes were in the interval 10.84--11.04 Mb and 29 in the interval 2.36--2.56 Mb on chromosome 4. We identified a total of 62 genes, of which 56 were expressed in stems and 24 were highly expressed in stems according to the Phytozome website. We used these 56 genes for pathway analysis, of which 19 (33.9%) were annotated in 16 pathways in the KEGG database (<http://www.kegg.jp/>) ([**Table 7**](#T7){ref-type="table"}, [**Figure 5**](#f5){ref-type="fig"}). Through the ensemble plant website, we identified homologs of the 56 candidate genes in *Arabidopsis thaliana*, and predicted their possible functions ([**Table 8**](#T8){ref-type="table"}).

###### 

Details of 19 candidate genes annotated in the KEGG database.

  QTN name           Gene name                 Chromosome   Position                 KO number    Annotation
  ------------------ ------------------------- ------------ ------------------------ ------------ ------------------------------------------------------------------------
  AX-157296110       *Glyma.13G334200.1*       chr13        42794878..42801645       K12881       THOC4; THO complex subunit 4
  **AX-157296110**   ***Glyma.13G334300.1***   **chr13**    **42807703..42808635**   **K11251**   **H2A; histone H2A**
  **AX-157296110**   ***Glyma.13G334500.1***   **chr13**    **42815054..42819704**   **K08678**   **UXS1; UDP-glucuronate decarboxylase \[EC:4.1.1.35\]**
  **AX-157296110**   ***Glyma.13G334800.1***   **chr13**    **42840520..42845919**   **K10575**   **UBE2G1; ubiquitin-conjugating enzyme E2 G1 \[EC:2.3.2.23\]**
  AX-157296110       *Glyma.13G335000.1*       chr13        42851863..42856533       K14833       NOC2; nucleolar complex protein 2
  AX-157296110       *Glyma.13G335600.1*       chr13        42916206..42918785       K10355       ACTF; actin, other eukaryote
  AX-157296110       *Glyma.13G335800.1*       chr13        42925815..42927692       K12471       EPN; epsin
  AX-157296110       *Glyma.13G336200.1*       chr13        42957851..42960820       K01803       TPI; triosephosphate isomerase (TIM) \[EC:5.3.1.1\]
  AX-157526374       *Glyma.04G029000.1*       chr04        2354851..2358499         K07056       rsmI; 16S rRNA (cytidine1402-2\'-O)-methyltransferase \[EC:2.1.1.198\]
  AX-157526374       *Glyma.04G029400.1*       chr04        2378504..2382290         K03064       PSMC6; 26S proteasome regulatory subunit T4
  AX-157526374       *Glyma.04G029700.1*       chr04        2405521..2408263         K20535       MPK1_2; mitogen-activated protein kinase 1/2 \[EC:2.7.11.24\]
  AX-157526374       *Glyma.04G029800.1*       chr04        2411114..2416765         K12402       AP4M1; AP-4 complex subunit mu-1
  **AX-157526374**   ***Glyma.04G030000.1***   **chr04**    **2429647..2437011**     **K10572**   **IPPK; inositol-pentakisphosphate 2-kinase \[EC:2.7.1.158\]**
  AX-157526374       *Glyma.04G030300.1*       chr04        2461743..2467168         K00975       glgC; glucose-1-phosphate adenylyltransferase \[EC:2.7.7.27\]
  **AX-157526374**   ***Glyma.04G030400.1***   **chr04**    **2466301..2470762**     **K19355**   **MAN; mannan endo-1,4-beta-mannosidase \[EC:3.2.1.78\]**
  AX-157526374       *Glyma.04G030500.1*       chr04        2481992..2488564         K00611       OTC; ornithine carbamoyltransferase \[EC:2.1.3.3\]
  AX-157526374       *Glyma.04G030700.1*       chr04        2493579..2495558         K01726       GAMMACA; gamma-carbonic anhydrase \[EC:4.2.1.-\]
  AX-157526374       *Glyma.04G030800.1*       chr04        2497398..2501640         K01726       GAMMACA; gamma-carbonic anhydrase \[EC:4.2.1.-\]
  AX-157526374       *Glyma.04G031600.1*       chr04        2525754..2530508         K01247       alkA; DNA-3-methyladenine glycosylase II \[EC:3.2.2.21\]

Bold font indicates genes that we propose correlate with plant height in soybean.

![KEGG pathway annotation for 19 candidate genes.](fpls-11-00009-g005){#f5}

###### 

Annotation information of homologous genes in *Arabidopsis thaliana*.

  Gene name               GO pathway annotation information   Homologous gene   Homologous gene description
  ----------------------- ----------------------------------- ----------------- ----------------------------------------------------------------------------------
  ***Glyma.13G035900***   **GO:0004672**                      ***AT1G21250***   **Wall-associated receptor kinase 1**
  *Glyma.13G036000*       GO:0006886                          *AT1G21380*       TOM1-like protein 3
  *Glyma.13G334200*       GO:0003676                          *AT5G59950*       RNA-binding (RRM/RBD/RNP motifs) family protein
  *Glyma.13G334500*       GO:0003854                          *AT2G28760*       UDP-XYL synthase 6
  *Glyma.13G334700*       --                                  *AT3G46450*       SEC14 cytosolic factor family protein / phosphoglyceride transfer family protein
  *Glyma.13G334800*       --                                  *AT3G46460*       UBC13
  *Glyma.13G335100*       --                                  *AT5G59350*       Transmembrane protein
  *Glyma.13G335300*       --                                  *AT2G28840*       Putative E3 ubiquitin-protein ligase XBAT31
  *Glyma.13G336000*       --                                  *AT1G07590*       Pentatricopeptide repeat-containing protein At1g07590, mitochondrial
  *Glyma.13G336500*       GO:0031047                          *AT1G13790*       Factor of DNA methylation 4
  *Glyma.04G029400*       GO:0006281                          *AT1G45000*       26S proteasome regulatory subunit S10B homolog B
  *Glyma.04G029500*       --                                  *AT1G75400*       RING/U-box superfamily protein
  *Glyma.04G029700*       GO:0005524                          *AT1G10210*       MPK1
  *Glyma.04G029700*       GO:0005524                          *AT1G59580*       Mitogen-activated protein kinase
  *Glyma.04G029900*       --                                  *AT1G19650*       Phosphatidylinositol/phosphatidylcholine transfer protein SFH4
  ***Glyma.04G030000***   **GO:0005524**                      ***AT1G22100***   **Inositol-pentakisphosphate 2-kinase**
  *Glyma.04G030300*       GO:0009058                          *AT1G27680*       Glucose-1-phosphate adenylyltransferase large subunit 2, chloroplastic
  *Glyma.04G031600*       GO:0006284                          *AT1G75230*       DNA glycosylase superfamily protein
  *Glyma.04G031800*       --                                  *AT1G75210*       Cytosolic IMP-GMP specific 5-nucleotidase, putative

Bold font indicates genes that we propose correlate with plant height in soybean.

Discussion {#s4}
==========

There are two parents, Dongnong L13 (110 cm) and Heongong 60 (70 cm), with significant differences in plant height, were distributed among 156 families. At the same time, the plant height observations of each family in the population showed a unimodal distribution, indicating that the plant height is a typical quantitative trait and controlled by multiple genes. The variation of plant height among the families in the population was large, and the genetic basis was rich. The range of RIL6013 is large than population conducted by "Minsoy" and "Noir 1", indicated this population is more suitable for QTL mapping. At the same time, the results of the *F* test showed that the differences between the families in each environment were extremely significant (P \< 0.01), also indicated that the population is suitable for analysis of variance and QTL mapping.

There is increasing research on QTLs related to soybean yield traits employing different methods. Most previous studies were single-environment analyses, affected by loss of related QTLs and data loss when using high LOD values or false positive detection of QTLs when using low LOD values. There are many detection methods for QTL mapping of soybean traits, including single marker analysis (SMA), interval mapping (IM), composite interval mapping (CIM), and inclusive complete interval mapping (ICIM). SMA is generally used in the case of few molecular markers, and genetic map construction is difficult. There are limitations on the development of QTL research using SMA ([@B15]), and genetic effects of QTLs and markers between the linkage distance cannot be distinguished. IM was first proposed by [@B32]. The possible location of QTLs can be inferred by marker interval, and the number of individuals required for the method is small. However, linkage between QTLs easily creates false positives, resulting in low detection efficiency of QTLs. [@B62] proposed a CIM method, which introduces molecular markers closely linked to other QTLs as background genetic control, thereby reducing the residual variance, eliminating the influence of other QTLs, and improving the accuracy of QTL mapping. However, this method cannot be used in complex situations such as those involving epistasis or environment interactions. Compared with CIM, [@B59] proposed that ICIM simplifies the process of controlling genetic background variation and improves the detection efficiency of QTLs ([@B38]). We used both IM and ICIM to effectively detect QTL intervals accurately in this study.

GWAS are effective methods for analyzing the genetic basis of complex traits in natural plant populations. Meanwhile, linkage analysis of segregating populations can effectively eliminate false positives, a defect of GWAS, but the results of linkage analysis usually have large intervals making it harder to find target genes. The four multi-locus GWAS methods can accurately map QTNs, and using two methods to verify each other can improve the accuracy of mapping. Analysis of multi-year and multi-environment data improves the accuracy of QTL position and effect estimation ([@B24]) and is conducive to searching for stable QTLs. We also used multi-locus GWAS to detect QTNs related to soybean plant height. Linkage analysis based on multiple methods and environments detected nine QTLs ([**Table 4**](#T4){ref-type="table"}), while GWAS analysis detected 10 QTNs ([**Table 5**](#T5){ref-type="table"}). These included *qPH-F-1* located on chromosome 13 at 1.91--43.31 Mb; AX-157278476 and AX-157296110 located at 11.11 Mb and 42.90 Mb on chromosome 13, respectively; *qPH-C1-1* located at 0.55--46.13 Mb on chromosome 4; and AX-157526374 and AX-157176525 at 2.46 Mb and 10.95 Mb on chromosome 4, respectively ([**Figure 6**](#f6){ref-type="fig"}, [**Table S4**](#ST4){ref-type="supplementary-material"}). Use of two methods not only narrowed the search range for candidate genes, but also supported the reliability of our results.

![Details of overlapping loci on two chromosomes.](fpls-11-00009-g006){#f6}

The population structure is an important factor leading to false positives in the results of association analysis, and the population structure is evaluated to determine the statistical model used. In statistics, confounding factors refer to irrelevant variable that are related to both dependent variables and independent variables ([@B44]). The population structure variables are confounding factors and are related to genotypes and phenotypes. At the same time, many traits are directly related to the subsets. There are differences in the phenotype and genotype between the families of the RIL population, so it is necessary to analyze the population structure of the RIL population. We performed multi-locus GWAS by a model removed Q, and the results were significantly different from the model Q+K ([**Table S5**](#ST5){ref-type="supplementary-material"}), indicating that the population structure has a great influence on the GWAS, so we adopted the model Q+K.

There is a close relationship between quantitative traits and environment. Plant height is a typical quantitative trait, and also affects soybean yield. The interaction between plant height and environmental conditions is an important factor affecting soybean yield. [@B19] reported that QTLs are sensitive to environmental performance and there are large differences in QTL detection between different environments, but QTLs for different traits have different stability, and those with high heritability are more easily detected across different environments. [@B12] believe that QTLs that detect higher effect values in a single environment are not as efficient as those detected in multiple environments. Among the QTNs mapped by GWAS, three (AX-157388309, AX-157278476, and AX-117466184) that could be located by multiple methods explained phenotypic variation ranging from 9.16% to 16.93%. These stable main QTNs can be used for fine mapping and related research on marker-assisted selection.

[@B51] proposed the soybean genetic map GmComposite2003 based on the SoyBase database (<https://www.soybase.org/>) and results of recent studies, covering 20 chromosomes of soybean, with a total length of 2562.28 cM, and including 3334 markers. We used this as a reference map for integrating the results of this experiment. We detected nine soybean plant height QTLs using multiple methods and in multiple environments, and eight of these QTLs---*qPH-D1a-2* ([@B33]; [@B7]; [@B18]), *qPH-D1a-3* ([@B18]), *qPH-N-2* ([@B26]), *qPH-C1-1* ([@B34]; [@B29]), *qPH-A2-1* ([@B34]), *qPH-B1-1* ([@B34]; [@B54]; [@B13]; [@B55]), *qPH-F-1* ([@B34]; [@B26]; [@B50]; [@B13]; [@B25]; [@B37]; [@B35]; [@B61]; [@B61]; [@B68]), and *qPH-L-1* ([@B33]; [@B34]; [@B36]))---showed inclusion within or overlap with chromosomal regions of QTLs detected in previous studies ([**Table 9**](#T9){ref-type="table"}). *qPH-H-1* detected in E1 and E2 on chromosome 12 is not reported on the genetic map. We consider this a newly discovered QTL with PVE ranging from 1.02% to 7.04%. The 10 QTNs located in this study are not overlapping with or close to previously described regions, and represent newly located QTNs, which can explain phenotypic variation of 6.11% to 16.93%.

###### 

QTLs located in previous studies.

  QTL         Chr   Marker interval    Reference
  ----------- ----- ------------------ --------------------------------------------------------------------------------
  qPH-D1a-2   1     Sat_413-Sat_160    [@B7]; [@B18]; [@B33]
  qPH-D1a-3   1     satt515-AZ302047   [@B18]
  qPH-N-2     3     Satt584-Satt660    [@B26]
  qPH-C1-1    4     Satt713-satt565    [@B29]; [@B34]
  qPH-A2-1    8     satt119-sat_406    [@B34]
  qPH-B1-1    11    satt197-sat_123    [@B13]; [@B54]; [@B34]; [@B55]
  qPH-F-1     13    sat_039-Sat_417    [@B26]; [@B13]; [@B50]; [@B61]; [@B35]; [@B25]; [@B37]; [@B68]; [@B34]; [@B61]
  qPH-L-1     19    sat_195-satt448    [@B36]; [@B34]; [@B33]

Based on the analysis of 19 significant genes in the KEGG pathway, we conclude that five genes may regulate the growth and development of soybean cells ([**Table 7**](#T7){ref-type="table"}, bold text), which in turn affects soybean plant height and yield. *Glyma.13G334500* is involved in the steroid biosynthesis process and regulates oxidoreductase activity, participates in the redox process, and is closely related to plant growth. *Glyma.04G030000* is involved in the binding of ATP and specifically binds to 5-adenosine triphosphate ([@B23]). *Glyma.04G030400* is involved in all carbohydrate pathways in the cell. *Glyma.13G334300* is involved in the biosynthesis of histone H2A, which plays an important role in regulating plant flowering, growth and development, immune response, mitosis, and DNA repair ([@B5]). *Glyma.13G334800* regulates the production of ubiquitin ligase, which plays an important role in plant growth and development and stress regulation, including plant height and seed weight ([@B40]).These genes are all thought to be closely related to the stem growth of soybean plants.

We identified 19 genes in Arabidopsis homologous to the 56 candidate soybean genes, two of which may be associated with soybean plant height ([**Table 8**](#T8){ref-type="table"}). *Glyma.13G035900* appears in the GO: 0004672 pathway. The homologous gene in Arabidopsis is *AT1G21250*, which regulates the synthesis of cell wall receptor kinase and controls cell wall elongation, which may be related to plant stem growth. *Glyma.04G030000* appears in the GO: 0005524 pathway; the homologous gene in Arabidopsis is *AT1G22100*, which is involved in the synthesis of ATP and is involved in cell growth and development.

Summary {#s5}
=======

We detected nine significant QTLs by linkage analysis and 10 QTNs using five multi-locus GWAS methods. Combining the two analysis methods, we obtained four significant QTNs and predicted five candidate genes closely related to soybean plant height around these four QTNs. We identified 19 homologous genes in Arabidopsis, two of which may regulate plant height and development.

Data Availability Statement {#s6}
===========================

All datasets generated for this study are included in the article/[**Supplementary Material**](#SF1){ref-type="supplementary-material"}.

Author Contributions {#s7}
====================

W-XL and HN conceived and designed the experiments. YF, QD, XL, ZQ, YW, XT, JS, JW, CY, and SJ performed the field experiments. SL and ZT performed the genome sequencing. YF, KZ, and HN analyzed and interpreted the results. YF and HN drafted the manuscript, and all the authors contributed to manuscript revision. WL provided laboratory conditions.

Funding {#s8}
=======

The authors gratefully acknowledge the financial support for this study provided by grants from the National Key Research and Development Program of China (2017YFD0101303-6) to W-XL and Project of Research and Development on Applied Technology of Harbin in Heilongjiang Province (2017RAXXJ019) to HN.

Conflict of Interest {#s9}
====================

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material {#s10}
======================

The Supplementary Material for this article can be found online at: <https://www.frontiersin.org/articles/10.3389/fpls.2020.00009/full#supplementary-material>

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

[^1]: Edited by: Yuan-Ming Zhang, Huazhong Agricultural University, China

[^2]: Reviewed by: Haiming Xu, Zhejiang University, China; Zhenyu Jia, University of California, Riverside, United States

[^3]: This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science
