Mapping QTL underlying plant height at three development stages and its response to density in soybean [Glycine max (L.) Merri.]

Abstract

To keep plant height (PH) stable in various densities, mapping QTL underlying te response of PH to density is of importance for soybean breeding. In this study, two soybean RIL populations were planted under two densities in two environments and PH at the pod-bearing stage (PBS), seed-filling stage (SFS) and maturity stage (MS) were measured to identify QTL controlling PH and its response to density increment. Expression of QTL controlling PH varied due to the change of seeding density, and the response to density (RD) of PH at SFS and MS surpassed that at PBS in both populations. Total 25 PH QTL and 8 RD QTL were identified which functioned on PH at three developmental stages. Of 25 PH QTL, 14, 8 and 8 were detected at PBS, SFS and MS, respectively. Three QTL could be expressed at two stages, and thirteen, seven and six QTL played roles separately at PDS, SFS and MS, respectively. Four QTL were expressed in the two seeding densities, and twelve and ten QTL showed effects specifically at low and high seeding density, respectively. Among eight RD QTL, five, two and one acted on PH at PBS, SFS and MS, and five and three QTL functioned in two environments, respectively. The information found in the present research lays the foundation for molecular design breeding on improving PH with stability in multiple densities in soybean.

Keywords:

• Soybean
• plant height
• planting density
• growth stages
• QTL

Introduction

Due to the high protein and oil content, soybean [Glycine max (L.) Merr.] has become a major source of plant oil and edible protein for humans and animals. Plant height (PH) is an important agronomic trait which determines soybean yield [1–5], for PH influences adaptive traits [6], the number of pods on the main stem [7],and the harvest index [3, 8], which affect the yield [2, 3, 6, 8–11]. PH is a quantitative trait that is controlled by multiple genes and is easily influenced by environmental conditions [2–4, 11–16]. With the development of molecular marker techniques, quantitative trait loci (QTL) have been gradually mapped to plant chromosomes [11, 15, 17, 18]. In the last 20 years, linkage maps constructed from multiple populations have been used for QTL mapping studies on plant height [1, 4, 9, 11, 13, 19, 20]. Lee et al. [3] identified a major effector locus for plant height, Dt1, on linkage group (LG) L, which could explain 67.7% of the total variation for plant height. Mansur et al. [9] identified a locus that could explain 31.6% of the total variation in plant height. Zhang et al. [20] identified three loci on LG C2 that could explain 20% of the total variation. By June 2019, 256 plant height QTL had been mapped to 20 soybean chromosomes (https://www.soybase.org/). Furthermore, Diers et al. [21] identified 15 QTL related to PH by a Nested Association Mapping (NAM) Recombinant Inbred Lines (RIL). Xue et al. [22] identified 19 novel QTL for PH in a four-way RIL population. Su et al. [23] mapped 9 QTL controlling PH. In previous studies, QTL for soybean PH were identified predominantly based on the phenotypes at maturation, which cannot explain the genetic effects at specific development stages [14]. In contrast, dynamic development trait mapping has been performed in some studies, such as in rice and cotton [24, 25]. Sun et al. [14] mapped dynamic development changes in plant height in soybean using 10-day intervals. Seeding density is one of the important cultivation techniques that affects PH in soybean, and generally PH extends with the increase of seeding density [26–30]. Additionally, different genotypes have different responses to the change in seeding density [31]. However, the influence of seeding density on the developmental process of PH and its genetic basis has not been explored.

In the present study, two soybean recombinant inbred line (RIL) populations (RIL3613 and RIL6013) with a common female parent were used to identify QTL associated with PH and its response to densities at the pod-bearing stage, seed-filling stage and maturity stage. The objective is to analyze the QTL basis of PH development under two seeding densities in soybean.

Materials and methods

Design and construction of the genetic populations

In 2005, three soybean cultivars with large plant height differences were selected, including Dongnong L13 (PH 135 cm) as a common female parent and Heihe 36 (PH 78 cm) and Henong 60 (PH 62 cm) as male parents to produce the crosses (Dongnong L13 × Heihe 36 and Dongnong L13 × Henong 60). Single seed descent was used in each generation to obtain the progeny, and the seperate individuals were derived from single seeds from the F₂ to F₅ generations. In the F₆ generation, individuals were planted to form plant lines, and two recombinant inbred populations, RIL3613 from Dongnong L13 ×Heihe 36 and RIL6013 from Dongnong L13 × Henong 60, were constructed with 134 and 156 lines [32].

Plant density and field experiment design

The entire field experiment was carried out in Keshan (E125°N48°) and Harbin (E 128°N45°) in 2015. In Keshan, the planting density was 2.6 × 105 plants/ha (termed 2015KD1) and 3.5 × 105 plants/ha (termed 2015KD2), and 2.2 × 105 plants/ha (termed 2015HD1) and 3.1 × 105 plants/ha (termed 2015HD2) in Harbin. The main planting factors are listed in Supplemental Table S1. A complete randomized block design with three replicates was adopted for the field experiments. Three-line plots were used with 5 m length and 0.67 m width. The plant distances for 2015KD1, 2015KD2, 2015HD1 and 2015HD2 were 6, 4.4, 7, and 10 cm, respectively.

Plant height measurements

Five plants were marked and their PH was measured at the pod-bearing period (R4 stage), seed-filling period (R6 stage) and maturity period (R8 stage), which were specified according Fehr et al. [33] for 2015KD1, 2015KD2, 2015HD1 and 2015HD2. The specific measurement method included randomly selecting five individuals from every line and measuring the distance from the cotyledonary nod to the top of the growing point. The mean value of five individuals was taken as the PH at different periods.

Data statistics and analysis

Descriptive statistics were used to evaluate the mean plant height values for each stage of each density for each environment. Analysis of variance (ANOVA) including the effects of environment, genotypes, densities and genotype × density interaction effects was conducted for the mean of PH, and the variance components were estimated by a mixed model method. The following equation was used to estimate heritability: $$h^2=\frac{{\sigma }_G^2}{{\sigma }_G^2+{\sigma }_{GD}^2/d+{\sigma }^2/de}$$ where, $h^2$ is heritability, ${\sigma }_G^2$ is the genetic variance component, ${\sigma }_{GD}^2$ is the variance of genotype × density interaction component, ${\sigma }^2$ is the error variance component, d is the level of density treatment, $e$ is the number of environments.

Microsoft Excel 2010 (Microsoft Corporation, Albuquerque, NM, US) was used for the descriptive statistics and analysis of the data. ANOVA was carried out using the generalized liner model (GLM) procedure in SAS 9.2 software, and ${\sigma }_G^2,$ ${\sigma }_{GD}^2$ and ${\sigma }^2$ were estimated via the MIXED procedure in SAS9.2.

Estimation of response to density

Response to density (RD) could be measured by the net effect of density on PH, which was carried out according to the condition variable method proposed by Zhu [34]. The specific calculation method for RD was to subtract the background of normal seeding density (the first density) from the PH phenotype of the higher density (second density). The specific equation is shown as follows: $$RD={PH}_2-C_{(1,2)}\cdot (P{H}_1-\overline{P{H}_1})/V_1$$ where PH₂ is the plant height at the second planting density, C_(1,2) is the covariance of plant height under the two planting densities, $\overline{P{H}_1}$ is the mean plant height under the first planting density, and V₁ is the variance of plant height under the first planting density.

RD was calculated using the formulas available in Excel 2010.

Mapping QTL

PH QTL mapping was carried out using the data of PH under various seeding density in each location, and RD QTL mapping was conducted via the estimated RD values.

Based on the constructed genetic map in a previous study [32], inclusive composite interval mapping (ICIM) [35], composite interval mapping (CIM) [36] and single marker analysis (SMA) was applied to detect GDE and DE QTL. The log of odd (LOD) scores for ICIM and CIM was 2.5, and a probability level of p < 0.05 of F-value from SMA was used as a significance threshold to detect QTL. QTL that were simultaneously detected by at least two methods and with a phenotypic variance explained (PVE) over 2% or coefficient of determination (R²) greater than 0.1 were regarded as QTL. ICIM-ADD was implemented using IciMapping 4.1 (http://www.isbreeding.net/software/). CIM and SMA were implemented using WinQTLCart 2.5 (https://brcwebportal.cos.ncsu.edu/ qtlcart/WQTLCart.htm).

Naming conventions for QTL

We used q to represent QTL, PH to represent plant height traits, and RD to represent response of PH to density. QTL were named according to qPH/qRDPH + chromosome name + sequence number. QTL that were mapped to the same marker region were given the same sequence number. Unified numbering was used for the two populations (RIL3613 and RIL6013).

Results and discussion

Phenotypic analysis

The genotype variation of PH for each population reached significance at 0.01 level under each density treatment in the studied environments, which showed that difference in gene expression existed between the two populations. Comparing the two populations, there were no distinct differences in the mean values, standard deviations and ranges, which showed that there was maybe a similar power to detect QTL in these two populations. Furthermore, the mean of variance from genotype × density interaction effects reached significance, which showed that the expression of QTL controlling PH varied according to the change of seeding density. The RD on PH varied depending on developmental stages, for both populations. The mean of RD at SFS and MS exceeded that at PBS, which showed that the influence of density on PH at SFS and MS surpassed that at PBS. The standard deviation and range of RD in Harbin was higher than that in Keshan, which indicates that the extension of the response of genotypes to density in Harbin was larger than in Keshan. There were significant variations between the two environments, and the means, standard deviations and ranges of PH showed difference. The significant genetic variance, genotype × density interaction effects variance and environmental variance of PH also indicated the possibility to identify different QTL under the variation of seeding density and plant locations. The heritabilities of PH ranged from 0.57 at PBS to 0.71 at MB for RIL3613, and from 0.42 at MB to 0.70 at SFS for RIL6013, respectively (Supplemental Table S2).

Kurtosis and skewness showed that the data of PH and RD of two the populations in all environments followed continuous normal distribution with an extensive phonemic variation, which indicate it is feasible to conduct genetic analysis based on these data (Supplemental Table S3, Supplemental Table S4 and Supplemental Figure S1, Supplemental Figure S2).

QTL detected by combination of multiple statistical methods

QTL mapping methods include SMA, CIM, ICIM, etc., and each mapping method uses different statistical principles with differing parameter estimation and significance testing approaches [31]. A combination of multiple statistical methods could decrease false positive QTL [37]. SMA is the simplest analysis method and only tests one marker at a time without background control. This may result in false positive QTL [38]. But the region where QTL are located in IM is easily disturbed by background variables, resulting in biased estimation. CIM and ICIM can improve QTL detection power and decrease deviations from estimations of QTL location and effects by background controlling [35]. While differences existed in selection of molecular markers as background which led to different QTL mapping results. Therefore, the joint use of multiple mapping methods can increase the reliability of the QTL identification results [39]. To decrease the false positive ratio, we used at least two methods simultaneously in QTL selection. Furthermore, there were differences in the estimation of genetic contribution of QTL to Phenotype according to the statistical methods. To select QTL with bigger contribution to target traits, the second standard was set of a Phenotypic variance explained (PVE) over 2% or coefficient of determination (R²) greater than 0.1. Combining the above two standards, 25 PH and eight RD QTL for PH were screened from 58, 35, and 570 QTL for PH and RD detected in the RIL3613 and RIL6013 populations using the CIM, ICIM, and SMA methods, respectively. The results supported each other and some false positive QTL that were present from the single mapping methods were filtered out, ensuring the authenticity of the detected QTL loci. Selecting QTL that are authentic and have high contribution rates according to the plant height target has significance in selective breeding.

Mapping of plant height QTL at different developmental stages

In the case of dynamic traits, the phenotype changes over time, which indicates that the genetic role on dynamic traits differs in the process of development. Based on the previous research on PH of rice [25], triticale [40], soybean [14], biomass and developmental stage of triticale [41], QTL associated with dynamic traits functioned at one stage specifically or multiple periods. In the present study, a total of 13, 8 and 6 QTL (Tables 1, 2 and 3) were detected by over two methods in the RIL3613 and RIL6013 populations at PBS, SFS and MS with PVE (ICIM) of 2.0% −12.63 or R² (CIM) of 0.08 - 0.42, PVE (ICIM) of 7.21%-11.8% or R² (CIM) of 0.06 – 0.32, PVE (ICIM) of 5.31%-11.69 or R² (CIM) of 0.10-0.30 for single QTL, respectively, which were located at the A2, B1, B2, C1, D1a, E, F, G, H, I, K, L, M and N linkage groups. Among the whole PH QTL, three QTL were identified in two developmental stages, i.e. qPH-C1-1 and qPH-M-4 at SFS and MS, and qPH-M-2 at PBS and SFS, whereas the remaining 21 QTL were specific for PBS, SFS or MS. In line with previous reports [14, 22, 24, 25, 43], in this experiment, QTL for PH are generally expressed dynamically, and PH is controlled by different combinations of QTL throughout PH development.

Table 1. QTL for plant height at pod bearing stage.

QTL	Left Marker	Right Marker	Genomic position † cM	Environment ‡	Population §	LOD ¶	a ǁ	PVE or R^2 lll
qPH-A2-1	Satt424	Sat_294	60.59-131.97	HD2	RIL3613	2.8/*/1.4/*	6.30/*/3.33/*	10.21/*/0.04/*
qPH-B2-1	satt577	Satt070	6.05-72.80	HD1	RIL6013	8.4/*/1.2/*	19.11/*/4.19/*	5.19/*/0.02/*
qPH-C1-2	Satt396	Sat_140	24.11-41.43	KD1	RIL3613	*/3.0/1.0/*	*/5.17/3.20/*	*/0.08/0.04/*
qPH-D1a-1	Sat_413	Sat_160	5.93-104.27	KD1	RIL6013	3.6/2.6/*/*	13.98/12.93/*/*	9.99/ 0.42/*/*
qPH-E-1	Satt720	Satt651	20.79-32.09	HD1	RIL3613	*/6.0/2.0/1.4	*/-8.59/-3.74/-3.36	*/0.38/0.07/0.05
qPH-I-1	satt367	satt330	27.98-77.83	HD1	RIL6013	3.2/3.1/*/2.2	15.37/13.78/*/5.92	2.00/0.25/*/0.00
qPH-I-2	Satt239	Satt354	36.93-46.22	HD1	RIL3613	*/3.3/1.5/0.9	*/6.22/3.89/2.93	*/0.13/0.05/0.01
qPH-K-1	Sat_243	satt673	50.79-86.77	HD1	RIL6013	5.3/*/1.0/*	−18.41/*/-4.21/*	4.98/*/0.03/*
qPH-L-2	Satt497	Sat_099	33.70-78.23	HD2	RIL3613	3.3/2.6/1.8/2.5	7.53/6.50/3.62/7.10	9.84/0.15/0.05/0.08
qPH-M-1	Satt567	Satt346	33.47-112.79	HD2	RIL3613	3.8/*/2.4/3.2	−5.87/*/-3.88/-4.46	10.65/*/0.08/0.09
qPH-M-2	Satt567	Satt697	33.47-85.34	KD1	RIL3613	3.0/3.2/2.5/*	−6.86/-7.04/-4.20/*	9.13/0.24/0.09
	Satt567	Satt697		KD2	RIL3613	3.9/2.9/3.5/*	−7.82/-5.33/-4.98/*	12.63/0.15/0.12/*
qPH-M-3	Satt346	Sat_121	103.98-112.79	HD2	RIL3613	*/3.2/*/1.0	*/-4.15/*/-2.48	*/0.12/*/0.03
qPH-N-1	Satt584	Sat_166	37.98-38.59	HD1	RIL6013	*/2.5/*/1.0	*/13.91/*/10.24	*/0.10/*/0.03
qRDPH-A2-1	Satt424	Sat_294	60.59-131.97	HRB	RIL3613	2.8/*/2.0/*	6.19/*/3.60/*	10.99/*/0.07/*
qRDPH-F-1	satt343	satt030	3.04-3.95	HRB	RIL6013	*/3.6/*/1.3	*/6.20/*/2.92	*/0.18/*/0.02
qRDPH-G-1	satt570	AW734137	12.74-15.63	KS	RIL6013	*/3.5/1.1/1.5	*/-3.96/-2.30/-2.65	*/0.09/0.02/0.01
qRDPH-G-2	sat_210	AW734137	3.71-5.63	KS	RIL6013	*/6.8/2.2/*	*/-7.21/-3.18/*	*/0.25/*/*
qRDPH-M-1	Satt567	Satt346	33.47-112.79	HRB	RIL3613	3.1/4.4/1.5/2.2	−5.28/-5.89/-2.91/-3.48	9.30/0.13/0.06/0.08

†: Genomic position means the genetic area in soybean public genetic map [42].

‡: HD1 means planting density of 2.2 × 10⁵ plants/hm² in Harbin; HD2 means planting density of 3.1 × 10⁵ plants/hm² in Harbin; KD1 means planting density of 2.6 × 10⁵ plants/hm² in Keshan; KD2 means planting density of 3.5 × 10⁵ plants/hm² in Keshan. HRB represent Location Harbin and KS represent Keshan.

§: RIL3613 refers to recombinant inbred lines population derived from Dongnong L13 × Heihe36; RIL6013 refers to recombinant inbred lines population derived from Dongnong L13 × Henong60.

¶: In LOD₁/ LOD ₂/ LOD ₃/ LOD ₄, LOD ₁ refer to log of odd (LOD) of inclusive compositive interval mapping methods (ICIM); LOD ₂ refer to log of odd (LOD) of compositive interval mapping methods (CIM); LOD ₃ refer to log of odd (LOD) of single marker analysis methods (SMA) for left marker; LOD ₄ refer to log of odd (LOD) of single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

ǁ: In a₁/ a₂/ a₃/a ₄, a₁ refer to additive effects derived from inclusive composite interval mapping methods (ICIM); a₂ refer to additive effects derived from composite interval mapping methods (CIM); a₃ refer to additive effects derived from single marker analysis methods (SMA) for left marker; a₄ refer to additive effects derived from single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

lll:In PVE/ R²₁/ R²₂/ R²₃, PVE refer to PHenotypic variation ratio of explanation of the QTL computed by inclusive composite interval mapping methods (ICIM); R²₁ refer to determination coefficient of the QTL computed by composite interval mapping methods (CIM); R²₂ refer to determination coefficient of the QTL computed by single marker analysis methods (SMA) for left marker; R²₃ refer to determination coefficient of the QTL computed by single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

Table 2. QTL for plant height at seed filling period.

QTL	Left Marker	Right Marker	Genomic position † cM	Environment ‡	Population §	LOD ¶	a ǁ	PVE or R^2 lll
qPH-B1-1	Sat_270	Satt359	21.99-102.55	HD1	3613	*/3.3/2.0/*	*/9.00/7.69/*	*/0.14/0.07/*
qPH-C1-1	satt565	Satt713	0.00-88.94	KD2	6013	2.9/4.1/1.2/1.8	6.39/7.89/3.02/3.62	7.21/0.24/0.00/0.02
qPH-F-1	sat_039	Sat_417	27.87-135.94	HD1	6013	*/3.0/1.9/*	*/12.18/5.75/*	*/0.32/0.00/*
qPH-G-1	AZ254740	Satt688	8.23-12.54	HD2	3613	*/3.9/2.2/1.1	*/-6.63/-4.79/-3.63	*/0.12/0.07/0.04
qPH-L-1	Satt182	Satt497	14.03-33.7	HD2	3613	*/3.2/* /1.1/	*/6.03/*/3.59	*/0.13/*/0.03
qPH-M-2	Satt567	Satt697	33.47-85.34	HD2	3613	2.6/3.2/2.1/*	−6.56/-7.36/-4.56/*	11.80/0.19/0.07/*
qPH-M-4	Satt536	Satt655	62.13-76.41	KD1	3613	3.4/*/1.7/2.6	5.11/*/3.50/4.33	9.12/*/0.06/0.09
qPH-N-2	Satt675	Satt080	34.65-45.13	HD2	3613	2.9/3.9/ 0.9/2.7	5.17/5.30/3.00/5.06	7.29/0.17/0.03/0.08
qRDPH-D1b-1	AI856451	Satt041	50.11-84.04	HRB	3613	*/3.2/*/2.8	*/-5.05/*/-5.03	*/0.12/*/0.08
qRDPH-M-2	Satt536	Satt655	62.13-76.41	HRB	3613	*/3.3/1.8/*	*/4.78/3.48/*	*/0.14/0.07/*

†: Genomic position means the genetic area in soybean public genetic map [42].

‡: HD1 means planting density of 2.2 × 10⁵ plants/hm² in Harbin; HD2 means planting density of 3.1 × 10⁵ plants/hm² in Harbin; KD1 means planting density of 2.6 × 10⁵ plants/hm² in Keshan; KD2 means planting density of 3.5 × 10⁵ plants/hm² in Keshan. HRB represent Location Harbin and KS represent Keshan.

§: RIL3613 refers to recombinant inbred lines population derived from Dongnong L13 × Heihe36; RIL6013 refers to recombinant inbred lines population derived from Dongnong L13 × Henong60.

¶: In LOD₁/ LOD₂/ LOD₃/ LOD₄, LOD₁ refer to log of odd (LOD) of inclusive compositive interval mapping methods (ICIM); LOD ₂ refer to log of odd (LOD) of compositive interval mapping methods (CIM); LOD ₃ refer to log of odd (LOD) of single marker analysis methods (SMA) for left marker; LOD ₄ refer to log of odd (LOD) of single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

ǁ: In a₁/ a₂/ a₃/a ₄, a₁ refer to additive effects derived from inclusive composite interval mapping methods (ICIM); a₂ refer to additive effects derived from composite interval mapping methods (CIM); a₃ refer to additive effects derived from single marker analysis methods (SMA) for left marker; a₄ refer to additive effects derived from single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

lll: In PVE/ R²₁/ R²₂/ R²₃, PVE refer to PHenotypic variation ratio of explanation of the QTL computed by inclusive composite interval mapping methods (ICIM); R²₁ refer to determination coefficient of the QTL computed by composite interval mapping methods (CIM); R²₂ refer to determination coefficient of the QTL computed by single marker analysis methods (SMA) for left marker; R²₃ refer to determination coefficient of the QTL computed by single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

Table 3. QTL for plant height at maturity period.

QTL	Left Marker	Right Marker	Genomic position † cM	Environment ‡	Population §	LOD ¶	a ǁ	PVE or R^2 lll
qPH-B1-2	Satt197	Satt359	46.38-102.55	KD1	3613	2.9/2.9/1.9/1.0	−5.68/-5.98/ −3.84/-2.94	5.31/0.13/ 0.07/0.04
				KD2	3613	2.6/*/1.8/*	−5.22/*/-3.82/*	8.62/*/0.06/*
qPH-C1-1	satt565	Satt713	0-88.94	KD2	6013	*/3.0/*/1.8	*/5.65/*/4.02	*/0.10/*/0.00
qPH-C1-3	Satt396	Sat_367	24.11-28.04	KD1	3613	*/3.0/1.7/1.2	*/5.47/4.25/3.27	*/0.10/0.05/0.04
qPH-G-2	Satt503	Satt288	68.76-76.76	KD1	3613	2.7/*/0.9/*	−8.16/*/-2.62/*	11.69/*/0.03/*
	Satt503	Satt288		KD2	3613	*/2.7/0.9/*	*/-8.53/-2.61/*	*/0.30/0.03/*
qPH-H-1	satt293	Satt181	89.08-91.12	HD1	6013	*/4.6/2.2/*	*/7.92/5.37/*	*/0.11/0.00/*
qPH-M-4	Satt536	Satt655	62.13-76.41	KD2	3613	2.7/4.4/1.8/2.5	5.07/5.67/3.75/4.35	8.50/0.15/0.06/0.09
qRDPH-B1-1	Sat_247	Satt484	49.73-118.52	KS	3613	*/4.6/3.1/*	*/5.33/3.79/*	*/0.23/0.04/*

†: Genomic position means the genetic area in soybean public genetic map [42].

‡: HD1 means planting density of 2.2 × 10⁵ plants/hm² in Harbin; HD2 means planting density of 3.1 × 10⁵ plants/hm² in Harbin; KD1 means planting density of 2.6 × 10⁵ plants/hm² in Keshan; KD2 means planting density of 3.5 × 10⁵ plants/hm² in Keshan. HRB represent Location Harbin and KS represent Keshan.

§:RIL3613 refers to recombinant inbred lines population derived from Dongnong L13 × Heihe36; RIL6013 refers to recombinant inbred lines population derived from Dongnong L13 × Henong60.

¶: In LOD₁/ LOD ₂/ LOD ₃/ LOD ₄, LOD ₁ refer to log of odd (LOD) of inclusive compositive interval mapping methods (ICIM); LOD ₂ refer to log of odd (LOD) of compositive interval mapping methods (CIM); LOD ₃ refer to log of odd (LOD) of single marker analysis methods (SMA) for left marker; LOD ₄ refer to log of odd (LOD) of single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

ǁ: In a₁/ a₂/ a₃/a ₄, a₁ refer to additive effects derived from inclusive composite interval mapping methods (ICIM); a₂ refer to additive effects derived from composite interval mapping methods (CIM); a₃ refer to additive effects derived from single marker analysis methods (SMA) for left marker; a₄ refer to additive effects derived from single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method.

lll: In PVE/ R²₁/ R²₂/ R²₃, PVE refer to PHenotypic variation ratio of explanation of the QTL computed by inclusive composite interval mapping methods (ICIM); R²₁ refer to determination coefficient of the QTL computed by composite interval mapping methods (CIM); R²₂ refer to determination coefficient of the QTL computed by single marker analysis methods (SMA) for left marker; R²₃ refer to determination coefficient of the QTL computed by single marker analysis methods (SMA) for right marker. * represent that the QTL could not be detected the corresponding method..

Stable and specific QTL at different density

Soybean plant height is a quantitative trait that is controlled by both genetic and environmental factors [11, 14, 16]. Therefore, QTL can be classified as stable QTL and environmentally specific QTL. The former refer to QTL that are expressed in multiple environments, whereas the latter represents QTL functioning in single environments. QTL mapping of these two types of loci is usually achieved through multi-year and multi-site experiments. Lee et al. [44] used single factor ANOVA to map six QTL that were stably expressed in multiple environments. In this study, we detected QTL controlling PH at two plant densities. At low density (2015HD1 and 2015KD1), 20 QTL associated with PH were detected only in one stage, whereas at high density (2015HD2 and 2015KD2), four QTL identified in two stages, i.e. qPH-M-2 at PBS and SFS, qPH-L-2 at PBS and MS, and qPH-G-1 and qPH-C1-1 at SFS and MS. There were three QTL detected simultaneously at the two seeding densities, i.e. qPH-M-2 at PBS, and qPH-B1-2 and qPH-G-2 at MS, and 13 QTL and 11 QTL were only discovered at low density (2015HD1 and 2015KD1) and high density (2015HD2 and 2015KD2) (Figure 1). By comparing the number of two kinds of QTL, it could be shown that selection of a variety suitable for a special environment is easier than breeding a line with a stable PH in multiple environments (Figure 2).

Figure 1. Temporal contribution of QTL.

Temporal development of PVE or R² of the QTL detected at PBS, SFS and MS for PH under low seeding density (Left) and high density (Right)

Figure 2. Venn diagram for QTL related to PH detected under two densities at three developmental stages.

PBS: Pod bearing stage; SFS: Seeds filling stage; MS: Maturity stage

D1: Planting density of 2.2 × 10⁵ plants/hm² in Harbin OR 2.6 × 10⁵ plants/hm² in Keshan; KD2 means planting density of 3.1 × 10⁵ plants/hm² in Harbin OR 3.5 × 10⁵ plants/hm² in Keshan

Genetic basis for the response to seeding density of PH

Density is an important environmental factor that affects plant height in soybeans: plant height increases generally as seeding density increases [26, 28, 29]. Different varieties have different responses to increase in seeding density [45]. Generally, the existence of the integral variation of this response could be tested by the significance of the genotype-by-environment interaction, but little attention has been paid to the response of single lines to density increase. In this research, we used RD to indicate the lines’ response to density increase, which showed wide variation for the two populations under two environments. Therefore, genetic factors play an important role in altering plant height when density increases. However, no studies on the genetic basis in soybeans in response to changes in density have been conducted previously. Here, based on the phenotypic data of RD, we identified five RD QTL (qRDPH-A2-1, qRDPH-F-1, qRDPH-G-1, qRDPH-G-2 and qRDPH-M-1), two RD QTL (qRDPH-D1b-1 and qRDPH-M-2) and one RD QTL (qRDPH-B1-1) underlying the response to the change in the density identified at PBS, SFS and MS with PVE (ICIM) of 9.30% −10.99% or R² (CIM) of 0.01 - 0.25, R² (CIM) of 0.07 – 0.08, R² (CIM) of 0.04-0.23 for single QTL, respectively, which were located at the A2, B1, D1b, F, G and M linkage groups. These QTL adjusted the phenotype of PH according to planting density through additive effects, i.e. the allelic genes from Dongnong L13 at qRDPH-A2-1, qRDPH-M-2, qRDPH-B1-1 and qRDPH-F-1, the allelic genes from Heihe 36 at qRDPH-M-1 and qRDPH-D1b-1, and the allelic genes from Henong 60 at qRDPH-G-1 and qRDPH-G-2 would increase PH as the planting density was increased. These were the molecular bases of the influence of density on plant height changes. Regarding the developmental stages, the number and effects of RD QTL in PBS (five QTL) were larger than those in SFS (two QTL) and MS (one QTL) (Tables 1, 2 and 3), which showed the most active response to density appeared in PBS. Under the effects of these QTL, soybean plant height exhibited large differences under different densities.

Co-mapping in multiple genetic backgrounds

By anchoring the QTL for PH found in the present and previous researches on a public linkage map [42], a comparison of the genomic region could be conducted (Supplemental Table S5). The corresponding locus of the qPH-M-4 QTL that was mapped in this study lies at 62.13–76.41 cM (Satt536–Satt655) in the public map. There are reports on five similar QTL in the region where this QTL is located [9, 15, 43, 46–48]. The corresponding locus of the qPH-B1-2 QTL that was mapped in this study lies at 46.38–102.55 (Satt197–Satt359) in the public map. Four similar QTL in the region where this QTL is located have been described [2, 14, 49]. The corresponding locus of the qPH-C1-1 QTL that was mapped in this study lies at 0.00–88.94 cM (satt565–Satt713) in the public map. Three similar QTL in the region where this QTL is located have been reported [2, 50]. The corresponding locus of the qPH-D1a-1 QTL that was mapped in this study lies at 5.93–104.27 cM (Sat_413–Sat_160) in the public map, and four similar QTL in the region where this QTL is located have been discussed [14, 43, 51, 52]. The corresponding locus of the qPH-G-2-2 QTL that was mapped in this study lies at 68.76–76.76 cM (Satt503–Satt288) in the public map. There is a report on one similar QTL in the region where this QTL is located [14].

The corresponding locus of the qPH-M-4 QTL that was mapped in this study lies at 33.47–85.34 cM (Satt567–Satt697) in the public map; and one similar QTL in the region where this QTL is located has been reported [46]. Although qRDPH-A2-1, qPH-B1-1, qRDPH-B1-1, qPH-B2-1, qPH-E-1, qPH-F-1, qRDPH-F-1, qRDPH-G-2, qPH-I-1, qPH-I-2, qPH-L-1, qPH-L-2, qPH-M-1, qRDPH-M-1 and qRDPH-M-2 were not mapped to multiple environments, these QTL have been mapped in some studies (Supplemental Table S5). The QTL qPH-C1-3, qPH-C1-2, qPH-G-1, qRDPH-G-1, qPH-H-1 and qPH-N-2 that were mapped in this study were not mapped in other studies and may constitute new loci, so these six QTL are novel and need further analysis.

Molecular marker-assisted selection protocol for plant height

The QTL for PH and RD identified in a specific combination of density and location could be used to design a molecular breeding scheme to achieve a breeding goal of PH. For QTL associated with PH, allele genotype of QTL carried by a specific parent that showed an additive effect on PH could be estimated by the posterior probability [35]. Then the breeding goal could be predicted in each developmental stage separately by the sum of additive effects from alleles in the combination of whole QTL (Figure 3). For example, of the combination of QTL detected in RIL3613 in MS under the planting density 2.6 × 105 plants/ha in Keshan, three QTL, i.e. qPH-C1-3 (Satt396- Sat_367), qPH-B1-2 (Satt197- Satt359) and qPH-G-2 (Satt503- Satt288), were detected with positive, negative and negative additive effect, respectively. There may be eight combinations of alleles which showed that the sum of additive effects ranged from -19.61 to 19.61 (Figure 4). There was significant difference of the sample means among the eight combinations of alleles (p < 0.01) (Supplemental Table S6, Figure 5). For each combination of alleles, the difference between the theoretical mean and the sample mean was non-significant (Supplemental Table S7). Depending on the analysis above, the allelic combination could be selected based on the breeding objective of PH, i.e. if the breeding objective of PH is 126.23 ∼ 137.17 cm, the allelic combination should be q1q2q3 and Q1q2q3, and if we hope to breed varieties with PH under 97.95 cm, the optimal allelic combination should be q1Q2Q3. In a similar way, we could choose an optimal allelic combination of QTL for a special environment.

Figure 3. Additive effects of QTL for PH and AED detected in two populations under two kinds of densities over two location.

In the name of QTL, 36- or 60- represent QTL detected in RIL3613 or RIL6013 population; D1 or D2 indicate QTL detected under the first or second kind of density; qph and qphae represent QTL of PH and AED; B1, I, …, represent linkage groups; the number, i.e., 1,2, 3,4 represent the sequence of QTL in one linkage group.

Figure 4. Additive effects of allelic combination of QTL qph-C1-3, qph-B1-2, qph-G-2 associated with PH at maturity stages.

Q1 and q1 represent allele of qph-C1-3 carried by Dongnong L13 and Heihe36, Q2 and q2 represent allele of qph- B1-2 carried by Dongnong L13 and Heihe36, Q3 and q3 represent allele of qph- G -2 carried by Dongnong L13 and Heihe36.

Figure 5. Theory and sample means of allelic combination of QTL qph-C1-3, qph-B1-2, qph-G-2 associated with PH at maturity stages.

Q1 and q1 represent allele of qph-C1-3 carried by Dongnong L13 and Heihe36, Q2 and q2 represent allele of qph-B1-2 carried by Dongnong L13 and Heihe36, Q3 and q3 represent allele of qph-G-2 carried by Dongnong L13 and Heihe36.

^{a, b, c} indicate the significance of difference between the phenotype means of allelic combination, i.e., no same letter means difference between two allelic combination was significant at 0.05 level, the reverse is also true.

Conclusions

Fourteen, eight and eight QTL controlling PH, and five, two and one QTL associated with response of PH to density were detected at PBS, SFS and MS. The individual and cumulative actions of these QTL lead to difference in development and final PH. The QTL identified in the present research are of importance for the molecular breeding design of PH.

Disclosure statement

No potential conflict of interest was reported by the authors.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This study was supported by the grants from Project of Research and Development on Applied Technology of Harbin in Heilongjiang Province (2017RAXXJ019) to H.N.