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.

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.

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).

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:

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 x̄D1 and V_D1 are the average and variance of MSNN under the density of D1, respectively.

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:

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, σG2 is the genetic variance, σE2 is the variance of density effect, σGD2 is the variance of the QTL × density interaction effect, σp2 is the phenotypic variance, and g, d, and r are the numbers of allelic genotype, density, and replication. On the basis of estimated σG2, σD2, σGD2, and σp2, the phenotypic variation explanation ratio (%) of additive (PVE_A) and additive × density interaction (PVE_AD) effect were estimated by the following formula:

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.

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/).

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).

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.

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).

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%.

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).

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.

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%.

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.

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).

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.

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.

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.

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.

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.