QTL for Main Stem Node Number and Its Response to Plant Densities in 144 Soybean FW-RILs

Li, Wen-Xia; Wang, Ping; Zhao, Hengxing; Sun, Xu; Yang, Tao; Li, Haoran; Hou, Yongqin; Liu, Cuiqiao; Siyal, Mahfishan; Raja veesar, Rameez; Hu, Bo; Ning, Hailong

doi:10.3389/fpls.2021.666796

ORIGINAL RESEARCH article

Front. Plant Sci., 20 August 2021

Sec. Plant Breeding

Volume 12 - 2021 | https://doi.org/10.3389/fpls.2021.666796

This article is part of the Research Topic Soybean Molecular Breeding and Genetics View all 25 articles

QTL for Main Stem Node Number and Its Response to Plant Densities in 144 Soybean FW-RILs

$\nWen-Xia Li&#x;$ Wen-Xia Li¹^†

Ping Wang^1,2^†

Hengxing Zhao¹

Xu Sun¹

Tao Yang¹

Haoran Li¹

Yongqin Hou¹

Cuiqiao Liu¹

Mahfishan Siyal¹

Rameez Raja veesar¹

Bo Hu¹^*

Hailong Ning¹^*

¹Key Laboratory of Soybean Biology, Ministry of Education, Key Laboratory of Soybean Biology and Breeding/Genetics, Ministry of Agriculture, Northeast Agricultural University, Harbin, China
²High Education Institute, Huaiyin Institute of Technology, Huai'an, China

Although the main stem node number of soybean [Glycine max (L.) Merr. ] is an important yield-related trait, there have been limited studies on the effect of plant density on the identification of quantitative trait loci (QTL) for main stem node number (MSNN). To address this issue, here, 144 four-way recombinant inbred lines (FW-RILs) derived from Kenfeng 14, Kenfeng 15, Heinong 48, and Kenfeng 19 were used to identify QTL for MSNN with densities of 2.2 × 10⁵ (D1) and 3 × 10⁵ (D2) plants/ha in five environments by linkage and association studies. As a result, the linkage and association studies identified 40 and 28 QTL in D1 and D2, respectively, indicating the difference in QTL in various densities. Among these QTL, five were common in the two densities; 36 were singly identified for response to density; 12 were repeatedly identified by both response to density and phenotype of two densities. Thirty-one were repeatedly detected across various methods, densities, and environments in the linkage and association studies. Among the 24 common QTL in the linkage and association studies, 15 explained a phenotypic variation of more than 10%. Finally, Glyma.06G094400, Glyma.06G147600, Glyma.19G160800.1, and Glyma.19G161100 were predicted to be associated with MSNN. These findings will help to elucidate the genetic basis of MSNN and improve molecular assistant selection in high-yield soybean breeding.

Introduction

Since soybean was one of the important crops worldwide, it has been an ongoing aim for soybean breeders to breed high yield cultivars in order to meet increasing global demand. As a major plant architecture trait, main stem node number (MSNN) affects soybean seed yield (Yao et al., 2015; Chang et al., 2018), for it is related with seed yield characters, such as logging, number of pods per plant, and days to flowering (Chapman et al., 2003; Zhang et al., 2004; Egli, 2013). MSNN is a typical quantitative trait, and the interaction of the genotype and the environment complicates the study on genetic basis. Therefore, molecular markers are widely used to locate quantitative trait loci (QTL) to reveal the molecular mechanism of MSNN in soybean yield. To date, Soybase (https://www.soybase.org/) has listed 37 QTL for MSNN by genetic linkage analysis (Zhang et al., 2004; Chen et al., 2007; Gai et al., 2007; Li et al., 2010; Liu et al., 2011; Moongkanna et al., 2011; Yao et al., 2015) and 11 quantitative trait nucleotides (QTNs) by genome-wide association study (GWAS) (Fang et al., 2017).

Genetic linkage analysis is an effective and traditional method to identify genetic intervals that associated plant phenotypes of traits (Tanksley et al., 1992). With the further development of DNA chip technology, single nucleotide polymorphism (SNP) has been widely used in high-density genetic linkage map construction and mapping of QTL (Hyten et al., 2010; Kim et al., 2010; Akond et al., 2013; Jun et al., 2014; Lee et al., 2015). GWAS can identify QTNs in genome regions based on the density of SNP, with advantages of high detection accuracy, high throughput, low cost, and time-saving. However, high false positive ratio is its inevitable defect. More researchers supported the viewpoint that the combination of linkage and association analysis was more accurate and effective than single methods (Ott et al., 2011; Liu et al., 2018; Fang et al., 2020; Zhang Y. C. et al., 2020). However, identifying the location of MSNN QTL for soybean using the combination of linkage analysis and GWAS analysis has not yet been discussed.

Population selection is the most important foundation for the linkage and GWAS analyses to map QTL. Linkage analysis typically is based on populations derived from two parents, and its detection power is usually relatively lower because of its less genetic diversity (Zhang S. et al., 2017). GWAS analysis generally uses natural populations or germplasm resources, and the population structure problem reduces the accuracy of the results. In order to solve the problems, scientists suggested constructing a special population, such as multi-parent advanced generation inter-cross (MAGIC) (Kover et al., 2009). The great opportunity for recombination in multiple parent populations increases mapping accuracy. Abundant genetic variation improves the efficiency of the detection of QTL, and clear kinship in progenies solves the population structure problem (Cavanagh et al., 2008). Kover et al. (2009) first constructed a MAGIC population with 19 Arabidopsis thaliana parents, and proved that the MAGIC population had great advantages in the location of QTL by mapping several known QTL with high precision. Huang et al. (2012) created a wheat MAGIC population from four excellent Australian varieties and identified QTL for plant height and hectoliter weight successfully. Butrón et al. (2019) also identified QTL for resistance to Fusarium ear rot in a MAGIC maize population.

In this study, in order to identify more accurate QTL and further perform gene mining precisely, an FW-RIL derived from a four-way cross was used to identify QTL for the MSNN of two densities in five environments by the combination of linkage analysis and GWAS. This research will enrich MSNN QTL and improve the precision of gene mining, as well as reveal the molecular mechanisms of MSNN in response to density, which will subsequently lay the foundation for marker-assisted selection breeding to increase soybean yield.

Materials and Methods

Plant Materials

To construct a four-way recombinant inbred line population, four soybean varieties with different node numbers in the main stem, Kenfeng 14, Kenfeng 15, Heinong 48, and Kenfeng 19, were used as parents. In 2008, two single crosses of Kenfeng 14 × Kenfeng 15 and Heinong 48 × Kenfeng 19 were obtained in Harbin, Heilongjiang province, China, and the F₁ was crossed as (Kenfeng 14 × Kenfeng 15) × (Heinong 48 × Kenfeng 19) in 2009. From 2010 to 2014, the progeny was self-crossed following the single seed descent method in Harbin and Yacheng, Hainan province, China. Finally, an FW-RIL population with 144 homozygous individuals was obtained and used for genetic map construction and mapping of QTL.

Field Experiment and Trait Measurement

The field experiment was conducted in Harbin (E126.63°, N45.75°) in 2015 (E1), Keshan (E125.64°, N48.25°) in 2015 (E2), Acheng (E127.63°, N45.82°) in 2016 (E3), Shuangcheng (E126.92°, E45.75°) in 2016 (E4), and Harbin in 2016 (E5). The parents and FW-RILs were planted in a three-row 5 × 0.7 m plot in a split block design of three replications. The main block arranged the plant densities, namely, 2.2 × 10⁵ plants/ha (D1) and 3 × 10⁵ plants/ha (D2). The sub-blocks were planted lines. The management procedures followed the normal production practices.

Five mature plants of the four parents and 144 four-way recombinant inbred lines (FW-RILs) were selected randomly in the middle of each row to measure MSNN before the harvest in the field for each replication. MSNN indicated the number of nodes from the cotyledonary node to the top of the main stem. The average of the three replications was used for phenotypic data analysis.

Genotyping and SNP Map Construction

Juvenile leaves were frozen in liquid nitrogen from the parents and FW-RIL plants, and then were ground into powder. Total genomic DNA was extracted with the CTAB method (Doyle et al., 1990) and eluted in 50-μl deionized water. SNP genotyping was conducted with SoySNP660K BeadChip at Beijing Boao Biotechnology Co. Ltd. A total of 109,676 SNPs were selected from 600,010 across 20 chromosomes, with minor allele frequency (MAF) > 0.05 and maximum SNP deletion locus <20% as criteria for the screening of SNP quality, and heterozygous loci were marked as missing to better estimate marker effect. Then, the locus was selected at each 100 kb interval along each chromosome from 3′-bottom to 5′-bottom. 2,292 high-quality SNPs on 20 chromosomes following Mendelian segregating ratio was applied to construct linkage map by the software GAPL V1.0 (Zhang S. et al., 2017). The length of the 20 linkage groups ranged from 76.4 to 329.7 cm, and the total length was 3,539.7 cm. The markers in each linkage group ranged from 16 to 316, with an average interval distance of 4.09 cm (ranging from 1.92 to 10.93 cm).

Statistical Analysis

Phenotypic Variation Analysis

The maximum, minimum, and standard deviations, skewness, and kurtosis of MSNN were calculated for each density in each environment. ANOVA was conducted with SAS V 9.2. ANOVA for single environment was carried out according to the following equation:

\begin{array}{l} x_{i j r} = μ + R_{r} + D_{j} + R D_{j r} + G_{i} + G D_{i j} + ε_{i j r} \end{array}

where x_ijris the rth observation of the ith genotype under the jth density in an environment; μ is the grand mean; R_r is the effect of main block r; D_j is the effect of density j; RD_jr is error of main block; G_i is the effect of genotype i; GD_ij is the interaction effect of genotype i by density j; and ε_ijr is the residual error, ε_ijr ~ N(0, σ²).

For multiple environments, joint ANOVA was conducted according to the following equation:

\begin{array}{l} x_{e i j r} = μ + E_{e} + E_{e} (R_{r}) + D_{j} + E D_{e j} + E_{e} (R D_{r j}) + G_{i} + G D_{i j} \\ + G E_{e i} + G D E_{e i j} + ε_{e i j r} \end{array}

where x_eijr is the rth observation of the ith genotype under the jth density in eth environment; μ is the grand mean; E_e is effect of eth environment; E_e (R_r) is the effect of rth main block in eth environment; D_j is the effect of density j; ED_ej is the interaction effect of density j by environment e; E_e (RD_rj) is error of main block in eth environment; G_i is the effect of genotype i; GD_ij is the interaction effect of genotype i by density j; GE_ei is the interaction effect of genotype i by environment e; GDE_eij is the interaction effect of genotype i by density j by environment e; and ε_eijr is the residual error, εe_ijr ~ N(0, σ²).

Genotype variance, genotype × density interaction variance, and error variance were estimated via a mixed model. The heritability (h²) for single environment was calculated with the following equation:

\begin{array}{l} h^{2} = \frac{σ_{G}^{2}}{σ_{G}^{2} + σ_{G D}^{2} / d + σ^{2} / d r} \end{array}

The heritability (h²) for multiple environments was calculated with the following equation:

\begin{array}{l} h^{2} = \frac{σ_{G}^{2}}{σ_{G}^{2} + σ_{G D}^{2} / d + σ_{G D E}^{2} / d e + σ^{2} / e d r} \end{array}

where h² is heritability; $σ_{G}^{2}$ is the variance of genotype; $σ_{G D}^{2}$ is the variance of genotype × density interaction; $σ_{G D E}^{2}$ is the variance of genotype × density × environment interaction; σ² is the variance of error; e is the number of environments; d is the number of planting density; and r is the number of repetitions.

Response to Density Estimation

Response to density refers to the difference in node number in response to change in density (from D1 to D2). Response to density (RD) could be evaluated according to the conditional variable method (Zhu, 1995) with the following equation:

\begin{array}{l} R_{D} = x_{D_{2}} - C_{D_{1} D_{2}} (x_{D_{1}} - {\bar{x}}_{D_{1}}) / V_{D_{1}} \end{array}

where R_D is the response to density; x_D1 is the phenotype value under the density of D1;x_D2 is the phenotype value under the density of D2; C_D1D2 is the covariance between phenotypes of MSNN under the two densities; and ${\bar{x}}_{D 1}$ and V_D1 are the average and variance of MSNN under the density of D1, respectively.

Linkage Analysis

Based on the SNP linkage map constructed above, interval mapping (IM) and inclusive composite interval mapping (ICIM) methods were used to map the QTL for MSNN in every density and environment through the PLQ function of GAPL software V1.0 (Zhang S. et al., 2017). In order to determine the existence of QTL, the scanning step was set to 1 cm, and the likelihood of odds (LOD) threshold was set to 3. The QTL were named qlNN-chromosome-sequence number or qlRDNN-chromosome-sequence number. The QTL mapped to the same marker region were given the same sequence number. QTL mapping results were mapped on chromosomes with MapChart2.1 (https://www.wur.nl/en). For QTL for MSNN detected in one interval, QTL by density effect in each environment, i.e., the additive effect over two densities, and additive × density interaction effect, were estimated. The formulas are shown as follows:

\begin{array}{l} μ_{. .} = \frac{1}{g d r} \sum_{i, j, k} y_{i j k} \\ μ_{i j} = \frac{1}{r} \sum_{i k} y_{i j k} \\ G_{i} = \frac{1}{d r} \sum_{j, k} y_{i j k} - μ_{. .} \\ D_{j} = \frac{1}{g r} \sum_{i, k} y_{i j k} - μ_{. .} \\ {G D}_{i j} = μ_{i j} - μ_{. .} - G_{i} - D_{j} \\ σ_{G}^{2} = \sum_{i} f_{i} G_{i}^{2} \\ σ_{D}^{2} = \frac{1}{d} \sum_{j} D_{j}^{2} \\ σ_{G D}^{2} = \sum_{i, j} f_{i j} {G D}_{i j}^{2} \\ σ_{p}^{2} = \frac{1}{g d r} \sum_{i, j, k} {(y_{i j k} - μ_{. .})}^{2} \end{array}

where y_ijk is the kth phenotype of ith allelic genotype in jth environment, μ_.. is the grand mean of all observation, μ_ij is the mean of ith allelic genotype in jth environment, G_i is the ith allele effect genotype of putative QTL, D_j is the jth density effect, GD_ij is the QTL × density interaction effect of ith allele genotype under jth density, $σ_{G}^{2}$ is the genetic variance, $σ_{E}^{2}$ is the variance of density effect, $σ_{G D}^{2}$ is the variance of the QTL × density interaction effect, $σ_{p}^{2}$ is the phenotypic variance, and g, d, and r are the numbers of allelic genotype, density, and replication. On the basis of estimated $σ_{G}^{2}$ , $σ_{D}^{2}$ , $σ_{G D}^{2}$ , and $σ_{p}^{2}$ , the phenotypic variation explanation ratio (%) of additive (PVE_A) and additive × density interaction (PVE_AD) effect were estimated by the following formula:

\begin{array}{l} {P V E}_{A} = σ_{G}^{2} \times 100 / σ_{p}^{2} \\ {P V E}_{A D} = σ_{G E}^{2} \times 100 / σ_{p}^{2} \end{array}

Genome-Wide Association Studies

The analysis of population structure was performed with the software STRUCTURE V 2.3.4. The number of subpopulations value (K) was determined with STRUCTURE HARVESTER (http://taylor0.biology.ucla.edu/structureHarvester). Linkage disequilibrium (LD) was analyzed with TASSEL 5.0. The K value was 2, and the LD was 1.63 Mb. The procedure is described in detail in a previous study (Zhang et al., 2018). Then, the GWAS was conducted with the software mrMLM.GUI V3.0 (Zhang Y. W. et al., 2020). Five multiple locus GWAS methods, mrMLM (Wang et al., 2016), FASTmrMLM (Tamba et al., 2017), FASTmrEMMA (Wen et al., 2018), pLARmEB (Zhang J. et al., 2017), and ISIS EM-BLASSO (Tamba et al., 2017), were used to identify significant QTL that control MSNN and its response to density. The probability P in the first step was set at 0.01 for mrMLM, FASTmrMLM, pLARmEB, ISIS EM-BLASSO, and 0.005 for FASTmrEMMA. The critical LOD score was set at 3 to determine significant QTL. The QTL were named qnNN-chromosome-sequence number or qnRDNN-chromosome-sequence number.

Candidate Gene Prediction

The QTL used to search candidate genes should satisfy the following conditions: (1) for QTL detected by linkage: should be detected in different densities, methods, or environments; explain the phenotypic variation more than 10%, and the interval length should be <600 kb; (2) for QTL detected by GWAS: should be detected in different densities, with more than two multiple locus GWAS methods, in multiple environments, or by co-location with QTL; and explain the phenotypic variation more than 10%. The Glyma.Wm82.a2.v1 gene model in Soybase (https://soybase.org/) was used to identify genes at the interval of each of the QTL (at the interval of 100 kb on either side, determined by the rate of LD decay). According to the Phytozome website (https://phytozome.jgi.doe.gov), genes highly expressed in the stem or shoot tip were selected among them. Then, the selected genes were put together to conduct pathway analysis on the Kyoto Encyclopedia of Genes and Genomes (KEGG) website (http://www.kegg.jp). Finally, the candidate genes were predicted through the results of pathway analysis combined with their homologous genes information on other crops and potential functions in GO number (https://www.ebi.ac.uk/QuickGO/) and the NCBI database (http://www.ncbi.nlm.nih.gov/).

Results and Analysis

Phenotypic Analysis

The summary of MSNN phenotype is presented in Table 1. The data showed that the node number of the parents and FW-RIL varied with density and environment. The range of FW-RIL covered the parents, which indicated strong bilateral transgressive segregation. The skewness and kurtosis values of the FW-RIL ranged from −1 to 1, and the phenotypic data displayed a typical normal distribution (Figure 1). All the characters of phenotypic variation indicated that MSNN was controlled by large- and small-effect QTL. The significant genotypic and genotype × environment interaction variance indicated a substantial genetic variation of MSNN existed among the FW-RILs and response of genotypes to environment varied among different environments. On the basis of the significant genotype × density interaction variance and genotype × density × environment interaction variance, it was implied that the MSNN response differed to the densities and that the response varied among the environments (Table 2).

TABLE 1

Table 1. Summarization of phenotype of nod number in main stem.

FIGURE 1

Figure 1. Frequency distribution of main stem node number (MSNN) in the four-way recombinant inbred line (FW-RIL) population in different treatments. (A–E) show the distribution in E1, E2, E3, E4, and E5, respectively; E1, Harbin in 2015; E2, Keshan in 2015; E3, Acheng in 2016; E4, Shuangcheng in 2016; E5, Harbin in 2016; D1, normal density (2.2 × 10⁵ plants/ha); D2, high density (3 × 10⁵ plants/ha).

TABLE 2

Table 2. Variance and heritability of nod number in main stem of four-way recombinant inbred lines.

Comparing the two densities, the mean of MSNN in D2 was higher than that in D1 (Figure 2), indicating the existence of MSNN response to density. The difference in heritability among the environments suggested the genetic basis for the formation of MSNN change according to environment. The extreme difference in heritability between joint and single environments showed that the response of genotype to change in density varied among environments. The whole variation in density under environments showed it was possible to detect different QTL in various environments.

FIGURE 2

Figure 2. Variation in main stem node number (MSNN) under two densities in five environments for four-way recombinant inbred line (FW-RIL). E1, Harbin in 2015; E2, Keshan in 2015; E3, Acheng in 2016; E4 Shuangcheng in 2016; E5, Harbin in 2016; D1, normal density (2.2 × 10⁵ plants/ha); D2, high density (3 × 10⁵ plants/ha).

Mapping of QTL for MSNN

In this study, 38 QTL for MSNN were detected on 18 chromosomes (except chromosome 7 and 15) with LOD value of over 3, which explained 3.44–14.93% of phenotypic variance (Figure 3, Supplementary Table 1). Among these QTL, 13 were identified in D1, seven in D2, and five in both D1 and D2. For 14 QTL underlying the response of MSNN to density, three and three were associated singly with MSNN in D1 and D2, and one of which was associated simultaneously with MSNN in D1 and D2 (Figure 4).

FIGURE 3

Figure 3. Distribution of quantitative trait loci (QTL) underlying main stem node number (MSNN) on 20 chromosomes. The red and blue colors represent QTL controlling NN and response to density change identified in present research, respectively, and the green color represents QTL underlying NN identified in previous research listed in Soybase (www.soybase.org).

FIGURE 4

Figure 4. Venn figure for frequency of linkage (left) and genome-wide association study (GWAS) (right) quantitative trail loci (QTL) underlying main stem node number (MSNN) and response to density increment. D1, normal density (2.2 × 10⁵ plants/ha); D2, high density (3 × 10⁵ plants/ha); RD, response to density.

In the five different environments, eight, three, three, six, and 12 quantitative trait loci (QTL) were detected in E1, E2, E3, E4, and E5, respectively. Six QTL could be found in more than two environments: qlNN-6-2 in E3 and E4; qlNN-8-2 in E3 and E5; qlNN-9-2 (qlRDNN-9-2) in E1, E2, and E3; qlNN-10-2 (qlRDNN-10-1) in E2 and E3; qlNN-13-1 (qlRDNN-13-1) in E4 and E5; and qlNN-17-1 (qlRDNN-17-1) in E1 and E5, respectively (Figure 5). Besides, 13 QTL could be detected with both the IM and ICIM methods, and in two densities and more than two environments, they could be considered as stable (Table 3). Among them, 12 QTL explained the phenotypic variation over 10%. qlNN-1-3 was detected in E2D1 with PVE of 11.01%; qlRDNN-3-1 was detected in E2RD with PVE of 7.57–12.07%; qlNN-6-1 (qlRDNN-6-1) was detected in E1D1, E1D2, and E1RD with PVE of 5.4–12.5%; qlNN-6-2 was detected in E3D1, E4D1, and E4D2 with PVE of 6.85–11.68%; qlNN-8-2 was detected in E3D1, E5D1, and E5D2 with PVE of 6.82–11.56%; qlNN-9-2 (qlRDNN-9-2) was detected in E1D2, E2D2, E2RD, and E3D2 with PVE of 6.63–10.50%; qlNN-9-3 was detected in E4D1 and E4D2 with PVE of 5.81–11.75%; qlNN-10-2 (qlRDNN-10-1) was detected in E2RD and E3D1with PVE of 6.12–11.02%; qlNN-13-1 (qlRDNN-13-1) was detected in E4RD and E5D2 with PVE of 5.16–10.07%; qlNN-17-1 (qlRDNN-17-1) was detected in E1D1 and E5RD with PVE of 3.44–12.96%; qlNN-18-1 was detected in E5D1 and E5D2 with PVE of 6.72–10.02%; and qlNN-19-1 was detected in E4D1 with PVE of 10.73–14.93%.

FIGURE 5

Figure 5. Venn figure for frequency of linkage (left) and genome-wide association study (GWAS) (right) quantitative trail loci (QTL) underlying main stem node number (MSNN) in five environments. E1, Harbin in 2015; E2, Keshan in 2015; E3, Acheng in 2016; E4, Shuangcheng in 2016; E5, Harbin in 2016.

TABLE 3

Table 3. Stable quantitative trait loci (QTL) for main stem node number (MSNN) identified under different densities in different environments with different methods.

Of all the QTL, the genome length of 15 QTL was <600 kb, which included six of the 13 ones: qlRDNN-3-1, qlRDNN-3-2, qlNN-6-1 (qlRDNN-6-1), qlNN-6-2, qlNN-17-1 (qlNN-17-1), and qlNN-19-1. The other stable QTL were repeatedly identified at a wide interval. In addition, all of the stable QTL with genome length of <600 kb could explain the phenotypic variation more than 10% except qlRDNN-3-2. Consequently, these intervals might play a critical role in mining genes to regulate MSNN.

Among all the alleles from the 13 stable QTL for MSNN, the parent Kenfeng 14 carried the positive additive effect alleles for 7 QTL, Kenfeng 15 for 11 QTL, Heinong 48 for 7 QTL, and Kenfeng 19 for 8 QTL. Four, four, one and one QTL from Kenfeng14, Kenfeng 15, Heinong 48, and Kenfeng 19 could obviously increase MSNN (additive effect > 1). Oppositely, Kenfeng 14 carried the negative additive effect alleles for 0 QTL, Kenfeng 15 for 6 QTL, and Heinong 48 and Kenfeng 19 for 10 QTL. Three, two, seven and four QTL from Kenfeng 14, Kenfeng 15, Heinong 48, and Kenfeng 19 could obviously decrease MSNN (additive effect < −1) (Table 3).

In the five environments, the total of PVE_A and PVE_AE varied extremely, ranging from 1.97 (in E2) to 45.98% (in E5). It was shown that the genetic basis of MSNN response to density varied in different environments (Table 4).

TABLE 4

Table 4. Additive and additive by density effect of quantitative trait loci (QTL) under two densities in each environment.

By comparison of PVE_A and PVE_AD, 15 QTL (qlNN-20-1, qlNN-12-1, qlNN-3-1, qlNN-3-2, qlNN-13-1, qlNN-19-1, qlNN-8-2, qlNN-8-1, qlNN-6-4, qlNN-17-2, qlNN-9-3, qlNN-18-1, qlNN-9-1, qlNN-18-2, and qlNN-12-2) expressed stably in two densities. Of these, qlNN-3-2 in E1, qlNN-9-2 and qlNN-14-1 in E4, and qlNN-3-1, qlNN-6-4, qlNN-8-1, qlNN-8-2, qlNN-9-1, qlNN-12-1, qlNN-12-2, qlNN-13-1, and qlNN-18-1 in E5 showed consistency with PVE_A over 2%.

Twelve QTL (qlNN-8-3, qlNN-2-1, qlNN-11-1, qlNN-14-1, qlNN-1-2, qlNN-17-1, qlNN-1-3, qlNN-10-2, qlNN-6-3, qlNN-10-1, qlNN-18-3, and qlNN-1-1) showed larger inconformity in various densities. Among these, qlNN-1-2 and qlNN-11-1 in E1, and qlNN-10-2 in E3 were expressed differently in specific density with PVE_AD more than 2%. Three QTL (qlNN-5-1, qlNN-9-2, qlNN-6-2) responded differently to density change in various environments.

QTL by GWAS

By GWAS analysis, QTL associated with MSNN were detected all over the genome on 18 chromosomes except chromosomes 16 and 17 (Figure 3). Thirty-four QTL were found in D1, 18 in D2, and 34 in RD, in which five were simultaneously found in D2 and RD (Figure 4). In other words, a total of 81 QTL were found, 47 of which could explain 10.1–38.38% phenotypic variation (Supplementary Table 2). From the different environments, 10, 26, 25, 9, and 10 QTL were identified specifically in E1, E2, E3, E4, and E5, respectively, and one QTL (qnRDNN-13-3) was repeatedly identified in E1 and E5 (Figure 5). Twenty-four stable QTL could be found with multiple methods or in the environments, 19 of which could explain phenotypic variation more than 10% (Table 5). qnNN-4-1 (qnRDNN-4-1) was identified in E4D2 and E4RD with PVE of 11.684–28.71%; qnNN-4-2 was identified in E5D1 with PVE of 15.1315–15.3806%; qnRDNN-5-3 was identified in E2RD with PVE of 8.3101–18.7664%; qnNN-6-2 was identified in E4D2 with PVE of 10.7297–17.4975%; qnRDNN-7-1 was identified in E5RD with PVE of 9.9725–24.8721%; qnNN-7-2 was identified in E2D1 with the PVE of 20.9359–26.3968%; qnNN-7-4 was identified in E2D2 with PVE of 12.1529–13.9883%; qnNN-9-1 was identified in E2D1 with PVE of 5.1691–10.674%; qnRDNN-9-2 was identified in E3RD with PVE of 11.9235–13.3229%; qnNN-10-2 was identified in E4D1 with PVE of 9.1925–10.6283%; qnNN-11-1 was identified in E5D1 with PVE of 26.5631–26.5932%; qnNN-12-1 was identified in E1D1 with PVE of 15.6919–22.8492%; qnRDNN-13-1 was identified in E2RD with PVE of 10.5682–20.7672%; qnNN-13-2 (qnRDNN-13-2) was identified in E3D2 and E3RD with PVE of 4.6444–16.7223%; qnNN-13-3 (qnRDNN-13-3) was identified in E1D2, E1RD, and E5RD with PVE of 14.1841–18.659%; qnNN-14-1 (qnRDNN-14-1) was identified in E1D2 and E1RD with PVE of 24.0645–38.3834%; qnNN-15-1 was identified in E4D1 with PVE of 13.4194–24.9554%; qnNN-18-3 was identified in E2D1 with PVE of 11.7346–20.1023%; and qnNN-19-2 was identified in E2D1 and could explain phenotypic variation of 4.3512–10.302%.

TABLE 5

Table 5. Genome-wide association studies for main stem node number (MSNN) detected by multiple methods under different densities in different environments.

Among the whole QTL identified by GWAS, 24 were co-located in the interval of QTL detected by linkage analysis (Figure 3), 15 of which could explain phenotypic variation more than 10%: qnRDNN-5-3 in qlNN-5-1, qnRDNN-9-2 in qlRDNN-9-1, qnNN-10-2 (qlRDNN-10-1) in qlNN-10-2, qnRDNN-13-1, qnNN-13-2 (qnRDNN-13-2) and qnNN-13-3 (qnRDNN-13-3) in qlNN-13-1 (qlRDNN-13-1), qnNN-18-3 in qlNN-18-1, qnNN-1-1 in qlNN-1-1, qnNN-1-3 in qlNN-1-3, qnNN-4-3 in qlRDNN-4-1, qnRDNN-5-1 in qlNN-5-1, qnRDNN-9-1 in qlNN-9-1, qnNN-9-4 in qlNN-9-2 (qlRDNN-9-2), qnNN-13-4 in qlNN-13-1 (qlRDNN-13-1), and qnNN-18-2 in qlNN-18-1. These QTL also could be considered stable because of the detection by linkage analysis and GWAS.

Candidate Gene Prediction

In this study, genes were screened based on the physical position of the five stable QTL (genome length <600 kb and PVE > 10%) and the 27 stable one (PVE > 10%) mentioned above. In total, 549 genes were found, among which 265 were highly expressed in the stem or shoot tip. Then these genes were used to conduct pathway analysis in the KEGG database (http://www.kegg.jp).

A total of 106 genes (which accounted for 40%) were annotated and divided into 36 catalogs and three protein families (Figure 6). Among these genes, four (Glyma.06G094400, Glyma.06G147600, Glyma.19G160800, and Glyma.19G161100) were speculated as potential candidate genes to regulate MSNN (Table 6).

FIGURE 6

Figure 6. Information on pathways (left) and orthologous protein families (right) of main stem node number (MSNN)-annotated candidate genes.

TABLE 6

Table 6. Details of four potential candidate genes for main stem node number (MSNN).

Discussion

Superiority of Using FW-RIL Population

In this study, a FW-RIL population was used for mapping, a simple mode of MAGIC (Kover et al., 2009) with four parents. It kept the advantage of the MAGIC population in abundant genetic variation. When mapping QTL, it could be applied to analyze allelic additive effects value from four parents. As long as there were differences between any two parents, QTL could be detected. For example, qlNN-1-1, qlNN-8-2, and qlNN-11-1 could not be detected in a bi-parent population derived from Kenfeng 14 × Kenfeng 15, because the allelic additive effects from the parents were approximately equivalent (Supplementary Table 1); so, the more allelic additive effect differences in the FW-RIL population, the greater the improvement made in QTL detection. Moreover, FW-RIL was an artificial population without the population structure problem, which was suitable for GWAS analysis.

Although the population size was relatively small, the combination of the linkage and GWAS analyses could improve mapping power. Furthermore, the experiment was conducted in three environments, which could compensate for the shortage in lower power in single environments. In summary, the statistical methods and multiple environment design could increase QTL detection power.

Combination of Genetic Linkage Analysis and GWAS Analysis

Both the genetic linkage and GWAS analyses were the main methods to identify genome regions related to quantitative traits. As mentioned above, research studies have already combined the two methods to conduct target trait location analysis. In soybean, the combination of the two methods was used in several traits, such as seed size and shape (Hu et al., 2013), seed protein and oil content (Zhang et al., 2019), number of pods (Song et al., 2020), and plant height (Fang et al., 2020). Similar to the MAGIC population, four parents carried multiple allelic genotypes in FW-RIL, so it could conduct linkage and association analysis (Zhang et al., 2018; Li et al., 2019, 2020, 2021; Liu et al., 2019; Qi et al., 2020; Song et al., 2020; Tian et al., 2020; Wang et al., 2021). First, the principles of linkage and association were different: the former associated an interval (region) with a target trait, and the latter associated a position (SNP) with a target trait. Second, the genotype data used in linkage and association analysis were different: the former was based on a small number of markers, and the latter depended on the large amount of markers over the whole genome, so the combination of linkage and GWAS could increase the identification of genomic regions associated with target traits in FW-RIL. In this research, the IM and inclusive composite interval mapping methods were used in the linkage analysis, and mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, and ISIS EM-BLASSO were used in the GWASs. A total of 38 QTL were identified by linkage analysis and 81 QTL were identified by GWAS, and 24 QTL were co-located in the interval of the identified QTL. The results indicated that the difference of the two methods in statistical principles and the genetic basis could complement each other and facilitate the detection of QTL.

Density Response of MSNN

Plant density is considered to be an important factor affecting soybean yield and yield components, such as MSNN. Ikeda et al. (1994) reported that soybean yield increased as density increased because of increase in total node number, especially branch node number. In this study, the MSNN of most of the lines increased as the density increased. Some lines showed the opposite response to density increase, indicating that the expression of gene for MSNN was probably affected by the change in density. By combining linkage and GWAS analyses, 55 QTL were identified in D1, and 33 QTL were identified in D2, respectively. Only five of them were identified in both of the densities, while the rest were detected in the single density. The results showed that the genetic basis of the QTL for MSNN was significantly different in the two densities, and for the genotype, environments, densities and their interaction were all at work. Inspired by the conditional genetic effects (Zhu, 1995) based on a net-effect analysis, the effect of MSNN response to density was estimated by the removal of other factors except planting density increment. In total, 48 QTL for MSNN response to density were identified when planting density increased from 2.2 × 10⁵ (D1) to 3 × 10⁵ plants/ha (D2), which is more valuable for molecular assistance selection (MAS) on MSNN in a specific planting density. Besides, in terms of the additive effects of all QTL for MSNN response to density, Heinong 48 and Kenfeng 14 were relatively suitable parents for increasing and decreasing MSNN in the MAS of soybean breeding, respectively.

Comparison of QTL Identified in Various Genetic Backgrounds

There were 25 QTL identified by linkage and 11 identified by GWAS listed in Soybase (https://www.soybase.org/search/qtllist_by_symbol.php). Among the 119 QTL identified this study, 10 had genome intervals that overlapped with published QTL node number (Figure 3). qlNN-2-1 was identified on chromosome 2 in genome intervals of 29,959,409–41,608,316 bp, overlapping with Node number 4-1 (38,221,027–40,699,300 bp) (Liu et al., 2011). qlNN-5-1was identified on chromosome 5 in genome intervals of 22,088,622–41,360,809 bp, overlapping with Node number 3-1 (35,971,621–38,939,759 bp) (Chen et al., 2007). qlNN-6-2 was identified on chromosome 6 in genome intervals of 11,860,267–12,150,538 bp, overlapping with Node number 5-1 (10,251,126–12,336,492 bp) (Moongkanna et al., 2011). qlNN-13-1 (qlRDNN-13-1) was identified on chromosome 13 in genome intervals of 444,838–43,052,819 bp, overlapping with Node numbers 1-5, 1-6, 1-7, 1-8 (Gai et al., 2007), and 2-3 (Zhang et al., 2004). qlNN-17-2 was identified on chromosome 17 in genome intervals of 7,296,590−9,660,500 bp, overlapping with Node number 7-1 (5,788,551–9,576,644 bp) (Li et al., 2010). qnNN-5-4 (37,951,491 bp) and qnRDNN-5-4 (38,349,709 bp) were identified in chromosome 5 and fell in the interval of Node number 3-1 (35,971,621–38,939,759 bp) (Chen et al., 2007); qnNN-6-3 (19,386,897 bp) was identified on chromosome 6 and fell in the interval of Node number 2-2 (19,370,872-20,218,893 bp) (Zhang et al., 2004). qnRDNN-13-1 (12,074,020 bp) and qnNN-13-1 (14,139,382 bp) were identified on chromosome 13 and fell in the interval of Node number 1-5 (10,199,530–15,306,234 bp) (Gai et al., 2007). The rest of 33 QTL identified by linkage and 76 QTL identified by GWAS were newly discovered, among which 37 with PVE > 10% were repeatedly identified with multiple density, environments, or methods (Supplementary Tables 1, 2). Consequently, this study probably would provide a great number of available genome regions and some potential high-confident candidate genes for MSNN.

Candidate Gene Related With MSNN

It is known that only few genes were directly related to MSNN in different crops. A novel ricMT gene was highly expressed in stem nodes (Yu et al., 1998). ZmMADS3 was expressed in the stem nodes of maize, and the transgenic maize reduced the number of nodes (Heuer et al., 2001). Dt1 controlled the number of nodes in soybean by regulating stem growth habit (Bernard, 1972). Therefore, it is of great significance to explore potential candidate genes for MSNN. In this study, four among 106 genes were predicted for MSNN.

Brassinosteroids are essential plant hormones with significant effect on cell proliferation and elongation. Glyma.06G147600 was annotated as protein brassinosteroid insensitive 1 (BRI1). It has been demonstrated that BRI1 is a receptor kinase that transduces steroid signals across the plasma membrane, which is likely to be the primary brassinosteroids (BR) receptor in Arabidopsis (Wang et al., 2001). Glyma.06G094400 was annotated as 14-3-3 protein epsilon. 14-3-3 proteins were highly conserved regulatory proteins, which interact with diverse target proteins in a sequence-specific and phosphorylation-dependent manner (Bridges and Moorhead, 2005). They have been proved to be involved in many processes of metabolism, hormone signaling introduction, cell division, and responses to abiotic and biotic stress in plants (Chen et al., 2006; Takahashi et al., 2007; Swatek et al., 2011). 14-3-3 proteins participate in BR signal transduction by regulating the subcellular localization and activity of both BZR1 and BZR2/BES1, which are the key transcription factor of BR signal transduction (Gampala et al., 2007). Chae et al. (2016) found that 14-3-3 proteins bound to BRI1, a kind of BR-receptor kinase, and phosphorylated in a BR-dependent manner, demonstrating that 14-3-3 proteins play an important role in the BR signaling of A. thaliana. Therefore, Glyma.06G094400 and Glyma.06G147600 could play an important role in MSNN because they probably would regulate stem growth through BR signaling pathway.

Glyma.19G160800 is annotated as tryptophan synthase alpha chain. Tryptophan synthase is an enzyme that catalyzes the final two steps in the biosynthesis of tryptophan, which could be converted to indole acetic acid (IAA) via the indole acetaldehyde or indole acetonitrile pathway. Glyma.19G161100 is annotated as auxin-responsive protein IAA. In other words, the two genes function in plant growth by IAA indirectly or directly. IAA is well-known for its strong effect on stimulating elongation in isolated stem segments (Yang et al., 1996), which has previously demonstrated that stem elongation strongly responded to exogenous IAA in light-grown pea (Murayama and Ueda, 1973; Yang et al., 1993). Recent research studies further showed that an auxin gradient was involved in cell proliferation in Arabidopsis and rice (Wang et al., 2018), and that auxin could convert to other forms to keep homeostasis to regulate soybean stem growth and development through various pathways (Jiang et al., 2020). Obviously, the two genes have a function in soybean stem growth and might have a certain relationship with MSNN. The four genes are all related with plant hormone signal transduction. It is necessary to verify the function of these genes in the future.

Conclusion

In this study, by combining linkage analysis and GWAS analysis, a total of 119 QTL associated with MSNN were identified in the FW-RIL population. Among them, 24 were simultaneously identified by the two methods. On the basis of the five QTL repeatedly detected in D1 and D2 and the 36 QTL for MSNN response to density, it was implied that a specific molecular mechanism controlled the MSNN response with the increase in plant density. In addition, 109 QTL were newly found, and four candidate genes were predicted to be closely related to MSNN. These genes could be of great value for MAS of soybean breeding.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

HN and W-XL: conceptualization. HZ, MS, RR, TY, HL, YH, CL, and XS: investigation. HN: resources, data curation, supervision, and funding acquisition. PW: writing—original draft preparation. BH and HN: writing—review and editing. W-XL: project administration. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Hundred-Thousand and Million Project of Heilongjiang Province for Engineering and Technology Science's Soybean Breeding Technology Innovation and New Cultivar Breeding (2019ZX16B01) to HN, the National Natural Science Foundation of China (32072015), and the Science Foundation of Heilongjiang Province, China (LH2021C029) to W-XL.

Conflict of Interest

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.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank Professor Zhixi Tian (State Key Laboratory of Plant Cell and Chromosome Engineering, Institute of Genetics and Developmental Biology, Chinese Academy of Sciences, Beijing, China) for conducting genotype analysis on the lines of FW-RIL in previous research studies.

Supplementary Material

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

References

Akond, M., Liu, S. M., Schoener, L., Anderson, J. A., Kantartzi, S. K., Meksem, K., et al. (2013). A SNP-based genetic linkage map of soybean using the SoyS-NP6K illumina infinium beadchip genotyping array. J. Plant Genome Sci. 1, 80–89. doi: 10.5147/jpgs.2013.0090