Quantitative Trait Locus Analysis of Protein and Oil Content in Response to Planting Density in Soybean (Glycine max [L.] Merri.) Seeds Based on SNP Linkage Mapping

Soybean varieties suitable for high planting density allow greater yields. However, the seed protein and oil contents, which determine the value of this crop, can be influenced by planting density. Thus, it is important to understand the genetic basis of the responses of different soybean genotypes to planting density. In this study, we quantified the protein and oil contents in a four-way recombinant inbred line (FW-RIL) soybean population under two planting densities and the response to density. We performed quantitative trait locus (QTL) mapping using a single nucleotide polymorphism (SNP) linkage map generated by inclusive composite interval mapping. We identified 14 QTLs for protein content and 17 for oil content at a planting density of 2.15 × 105 plant/ha (D1) and 14 QTLs for protein content and 20 for oil content at a planting density 3.0 × 105 plant/ha (D2). Among the QTLs detected, two oil-content QTLs was detected at both plant densities. In addition, we identified 38 QTLs for the responses of protein and oil contents to planting density. Of the QTLs detected, 70 were identified in previous studies, while 33 were newly identified. Fourty-five QTLs accounted for over 10% of the phenotypic variation of the corresponding trait, based on 23 QTLs at a marker interval distance of ~600 kb detected under different densities and with the responses to density difference. Pathway analysis revealed four candidate genes involved in protein and oil biosynthesis/metabolism. These results improve our understanding of the genetic underpinnings of protein and oil biosynthesis in soybean, laying the foundation for enhancing protein and oil contents and increasing yields in soybean.


INTRODUCTION
Soybean (Glycine max [L.] Merri.) is a major source of vegetable oil and feed protein. Dry soybean seeds are composed of approximately 40% protein and 20% oil (Rajcan et al., 2005). Increasing the protein content (PC) and oil content (OC) of soybean seeds is an important breeding objective.
The protein and oil contents of seeds are inherited as quantitative traits controlled by multiple genes, leading to a low efficiency of soybean improvement based on phenotypic selection. Therefore, much research has focused on quantitative trait locus (QTL) mapping for protein and oil content. To date, 248 QTLs for protein content and 327 for oil content have been deposited in SoyBase (http://www.soybase.org). Most of this research was carried out using conventional molecular marker techniques (Powell et al., 1996), such as restriction fragment length polymorphism (RFLP), amplified fragment length polymorphism (AFLP), and simple sequence repeat (SSR) mapping, which has reduced the accuracy of QTL mapping and led to reduced genetic stability due to the low densities of these markers and their uneven distributions in the genome. In recent years, single nucleotide polymorphism (SNP) markers have become powerful tools for exploring plant genomes due to their high density, good genetic stability, and suitability for accurate, high-throughput genotyping (Lee et al., 2015;Zhang et al., 2018). SNPs markers have been successfully used in soybean research. For example, Liu et al. (2017) used parental Zhong huang × Hua xia3 hybridization RIL populations obtained using SNP markers and detected four QTLs for protein content and nine QTLs for oil content in two environments in 2 years. Akond et al. (2014) detected one QTL for protein content and 11 QTLs for oil content using an F 5 : 8 RIL population derived from MD 96-5722 × Spencer using 5376 SNP markers. Wang et al. (2015) used an RIL population from V97-1346 × R05-4256 and detected 13 SNPs for protein content. In these studies, QTL mapping for protein and oil content in soybean plants was conducted using bi-parental populations. This technique has several limitations, as it results in the detection a small number of QTLs and limits the richness of the alleles and phenotypic variation. However, more recently, four-way recombinant inbred line (FW-RIL) populations have been used for QTL mapping, which overcomes the limitation posed by populations with a narrow genetic basis. In addition, these lines contain four alleles in a single locus, which greatly improves the QTL detection capacity. Moreover, high-density linkage maps can be constructed using FW-RIL populations, and the genetic markers are highly polymorphic (Ning et al., 2016;Liu et al., 2019).
In recent decades, the planting density of soybean has gradually been increasing, which has greatly affected the quality of soybean seeds. Bellaloui et al. (2014) investigated the effects of planting density on the composition of soybean seeds and observed that the highest protein and the highest oil contents occurred at different planting densities. Akond et al. (2012) detected two QTLs for protein and 6 QTLs for oil content at a higher planting density (25 cm row spacing) and three QTLs for oil content at a lower plant density (50 cm row spacing). By comparing the phenotypes and QTLs detected, many different genetic bases for protein and oil content at different planting densities have been uncovered. It is therefore important to breed plant varieties that are suitable for various planting densities and to analyze the genetic basis of the effects of planting densities on protein and oil content.
In the current study, we analyzed the genetic basis of the responses of protein and oil contents in soybean to different planting densities based on phenotypic data from a 144 member FW-RIL population and high-density SNP maps. We identified 45 QTLs with phenotypic variance explained (PVE) values of >10% under different planting densities and with the responses to density difference. Finally, we identified four candidate genes that might control protein and oil content under the influence of planting density that could be useful to improve molecular breeding and increase the protein and oil contents of soybean.  (E5). The field experiments were arranged in a split-plot design with three replications. The main treatment comprised a normal planting density of 2.15 × 10 5 plants/ha (D1) and an increased planting density of 3.0 × 10 5 plants/ha (D2), which were implemented using a plant spacing of 0.07 and 0.05 m, respectively, in each row of a three-ridge block 3 meters long and 2.1 meters wide. The four parents and FW-RILs were planted in a subplot with a random permutation. The environmental conditions, fertility, year, planting date, planting location, annual rainfall, and annual accumulated temperature are summarized in Table S2. Field management was performed using local field cultivation conditions for soybean.

Plant Materials
Five mature plants obtained from the four parental lines and the FW-RIL population were randomly selected from the middle row of each plot. Protein and oil content were measured using an INFRAC TM1214 near-infrared grain quality analyzer (Foss). The detailed process is described in Zhang et al. (2018). The average values of three samples were used as the phenotypic values of the parents and FW-RILs.

Statistical Analysis
Estimating the Effect of Density The response to density (RD) value refers to the change in the quality of a trait due to an increase in the planting density. RD was calculated as the remainder of the trait value at high density (D2) after subtracting the trait value at normal density (D1). Specifically, RD was estimated using the conditional variable method (Zhu, 1995) with the following formula: where RD is the response to density difference, y D1 is the phenotypic value of the first (normal) density, y D2 is the phenotype value of the second (high) density, C D1D2 is the covariance between phenotypes of quality traits under the two densities, and y D1 and V D1 are the average and variance, respectively, of quality traits under the first (normal) density.
The maximum, minimum, standard deviation, range, skewness, and kurtosis of the protein and oil content for each density in each environment were also calculated. Analysis of variance (ANOVA) was conducted to compare the phenotypic values of protein and oil contents at each planting density in each environment or jointly in multiple environments and in response to density difference.
Similar to the method used to estimate average heritability in multiple macro-environments, the average broad-sense heritability (h 2 ) at multiple planting densities was estimated using the following equation: where h 2 is broad-sense heritability, σ 2 G is the variance of the genotype, σ 2 GD is the genotype × density interaction, σ 2 is the variance of error, e is the number of environments, d is the number of planting densities, and r is the number of replications. The data for protein and oil content were analyzed using Proc MIXED in SAS 9.2 statistical software (SAS Institute, Cary, NC, US).

SNP Genotyping
DNA was extracted from the fresh leaves of plants from the FW-RIL population from the cross (Kenfeng 14 × Kenfeng 15) × (Heinong 48 × Kenfeng 19) using the CTAB method (Doyle and Doyle, 1990). The DNA was used for SNP genotyping analysis with the SoySNP660K BeadChip at Beijing Boao Biotechnology Co., Ltd. The SNP markers were filtered for minor allele frequency (MAF > 0.05), with a maximum of <10% missing sites per SNP (Belamkar et al., 2016). A linkage map of soybean containing 2332 SNP markers (https://figshare.

QTL Analysis
Based on the SNP linkage map, the average protein and oil contents under each planting density and RD in every environment were used to map the QTLs with the inclusive composite interval mapping (ICIM) method (Zhang S. et al., 2017) using GAPL V1.2. Firstly, The LOD (likelihood of odds) score for putative QTLs was determined after 1,000 permutations at a significant level of P < 0.05 with objective to find major QTL. Then mapping QTL was re-analyzed by setting LOD score of 3 to screen minor QTL. The QTLs were named as follows: q-trait name-chromosome name-sequence number, where q represents QTL, PC, and OC represent protein content and oil content, respectively, and RD represents the response to density. The QTLs that were mapped to the same marker region were given the same sequence number.

Identification of Candidate Genes for Protein and Oil Content
Twelve QTLs and eleven QTLs with PVE > 10% were detected within a 600 kb genomic region under different planting densities and with the responses to density difference, respectively. To further explore whether these QTLs are related to protein and oil content in soybean, we attempted to identify the candidate genes associated with the QTLs. We used the Glyma.Wm82.a2.v1 gene model in SoyBase (https://soybase.org/SequenceIntro. php) to identify all gene sequences based on the intervals of each QTLs with PVE > 10%. As a result, various highly expressed genes that controlled protein and oil content were identified based on the BAR website (http://www.bar.utoronto.ca). Finally, we used the KEGG website (http://www.kegg.jp/blastkoala/) to identify candidate protein and oil content-related genes based on KEGG pathway analysis results.

Phenotypic Variation
We constructed bar chats of the protein and oil contents of 144 FW-RILs (Figure 1). Protein and oil contents at different planting densities in the same environment exhibited a continuous normal distribution. Analysis of the protein and oil contents of the four parents at different densities located in the interval between the minimum and maximum values of the FW-RIL population (Table 1) indicated that the protein and oil contents of the FW-RILs was significantly higher than those of the parents, suggesting that transgressive inheritance was involved in protein and oil contents in the FW-RILs. Most of the kurtosis and skewness values of the data were between −1 and 1, suggesting that the population was suitable for ANOVA of related traits. An F test revealed significant (P < 0.01) differences in quality traits among the densities in a single environment. Therefore, the FW-RILs were constructed from two high-oil (Kenfeng14, Kenfeng15) and two high-protein varieties (Kenfeng19, Heinong48), providing an ideal basis for QTL analysis. ANOVA of the values of the parents and FW-RILs indicated that the genotype and genotype × density interaction effect was significant (P ≤ 0.01; Table 2). Compared to direct density effects, the genotype × density interaction effect was significant, indicating that density affects protein and oil content indirectly by altering the genetic basis of quality formation and that different genotypes have different responses to density increase.
The estimated broad-based heritability varied for protein and oil content at different planting densities in each environment, indicating that density influences the genetic basis of protein and FIGURE 3 | Distribution of QTLs for protein content on linkage groups detected in the soybean FW-RIL population. Red and blue represent QTLs detected for protein content under different planting densities and the response of density difference, respectively. oil content in soybean seeds ( Table 2). The effect of environment on the response of protein and oil content to density was significant (P ≤ 0.01; Table 3), indicating that the response of protein and oil content to planting density differed depending on the environment. In addition, the effect of increased planting density on protein and oil content varied markedly in terms of both direction and magnitude according to genotype (Figure 2). Therefore, we analyzed the effects of QTLs on protein and oil contents at different planting densities.

QTLs for PC and OC at Two Planting Densities
In this study, 65 QTLs associated with PC and OC were detected on 15 and 14 of the 20 soybean chromosomes under LOD value of 3 and under permutation method, respectively (Figures 3, 4), which 5 PC QTLs and 6 OC QTLs were located in two method. Of these, 28 PC QTLs (Table 4) and 37 OC QTLs (Table 5) were detected at different planting densities (D1 and D2). Among the QTLs detected, 14 for PC and 17 for OC were detected at D1. The LOD values ranged from 3.01 to 6.43. A single QTL accounted for 5.22% (qPC-17-2) to 14.58% (qOC-1-1) of phenotypic variance. The 14 remaining QTLs for PC and 20 QTLs for OC were detected at D2. The LOD values ranged from 3 to 8.93. The PVE values of the QTLs ranged from 4.77% (qOC-7-3) to 25.05% (qPC-6-1). Finally, 26 QTLs accounted for over 10% of the phenotypic variation. These findings indicate that protein and oil contents are controlled by both major-effect (PVE ≥ 10%) and minor-effect (PVE < 10%) QTLs (Figures 5, 6). Among these QTLs, only two QTLs (qOC-7-3, qOC-15-1) was simultaneously identified at both planting densities (Figure 7), indicating that the genetic basis for protein and oil content differed at different densities.
Of the QTLs detected, one protein-content QTLs (qPC-18-2) and three oil-content QTLs (qOC-18-1, qOC-10-1, qOC-5-1) FIGURE 4 | Distribution of QTLs for oil content on linkage groups detected in the soybean FW-RIL population. Red and blue represent QTLs detected for oil content under different planting densities and the responses to density difference, respectively.
were detected in more than two environments with PVE of 6.87-9.16%, while the 27 remaining protein-content QTLs and 34 oilcontent QTLs were detected in specific environments (D1: five PC QTLs and two OC in E1, three PC QTLs and two OC QTLs in E2, two PC QTLs and two OC QTLs in E3, two PC QTLs and five OC QTLs in E4, one PC QTLs and four OC QTLs in E5; D2: two PC QTLs and one OC QTLs in E1, one OC QTLs in E2, six PC QTLs and four OC QTLs in E3, two PC QTLs and six OC QTLs in E4, four PC QTLs and seven OC QTLs in E5). These results indicate that the genetic basis for protein and oil content differed in various environments (Figure 8).

QTLs for Response to Density
We identified 38 QTLs for the response of protein and oil content to density (RD) (Figures 3, 4) under LOD value of 3 and under permutation method, which 5 RDPC QTLs were located in two method. Briefly, 20 RD QTLs for protein content were identified on various chromosomes, including Chr01, Chr03, Chr05, Chr06, Chr07, Chr09, Chr12, Chr13, Chr14, Chr15, Chr17, Chr18, and Chr19, with LOD values ranging from 3.05 to 13.68 and accounting for a phenotypic variance of 3.64% (qRDPC-3-2)-38.02% (qRDPC-12-1). Moreover, 18 QTLs for oil content  Diers et al., 1992;Mansur et al., 1996;Orf et al., 1999;Chapman et al., 2003;Tajuddin et al., 2003;Eskandari et al., 2013;Lu et al., 2013;Mao et al., 2013;Asekova et    were located on 10 different chromosomes. The PVE values of single QTLs for oil content ranged from 5.88% (qRDOC-7-2) to 24.68% (qRDOC-6-2), with a LOD value ranging from 3.04 to 4.85 (Table 6). There were 19 RD QTLs with a PVE >10%, indicating that the RD trait is controlled by both major-effect and minor-effect QTLs (Figure 9). Among the RD QTLs, one QTL for protein content were detected in two environments, and the 19 remaining RD QTLs for protein content and 18 for oil content were environment specific (Figure 10). The finding that not all RD QTLs were detected in all environments might explain the differences in protein and oil contents in response to planting density. Among the RD QTLs, 11, 8, 12 and 9 positive alleles from Kenfeng14, Kenfeng15, Heinong48 and Kenfeng19, respectively increased the protein content at higher planting density. Similarly, 8, 13, 6 and 7 positive alleles from Kenfeng14, Kenfeng15, Heinong48 and Kenfeng19, respectively increased the oil content when the density increased from D1 to D2 ( Figure S1).

Analysis of Potential Candidate Genes
We identified 484 genes based on 23 QTLs under different planting densities and with the responses to density difference. Seventy-six genes were highly expressed within these regions in seeds, with 26 annotated genes in 20 pathways and 3 protein families identified in the KEGG database (Figure 11). Among these, four genes were selected as potential candidate genes affecting protein and oil content due to their annotations and functions in various metabolic pathways ( Table 7; see bold text).

QTLs With Effects on Protein and Oil Content in Soybean at Different Planting Densities
Protein and oil contents are quantitative traits that are affected by environmental conditions. Competition for nutritional resources occurs when the planting density changes, and different plant   varieties can respond differently. In the current study, the protein content of four soybean varieties and the oil content of five soybean varieties increased when the planting density increased from D1 to D2 in five environments (Figure 2). However, the responses of the remaining varieties to density were inconsistent, indicating that the environment affects the expression of genes that control protein and oil content. Bellaloui et al. (2014) used four soybean varieties to study the effects of planting density on protein and oil content. The protein contents of soybean varieties P93M90, AG 3906, P94B73, and V52N3 reached their maximum levels at a planting density of 518,700, 518,700, 180,000, and 150,000 plants ha −1 , respectively. Similarly, the oil contents of the four varieties were highest at a planting density of 247,000, 444,600, 150,000, and 60,000 plants ha −1 , respectively. These results indicate that the responses of individual protein and oil contents vary depending on the planting density. Therefore, the results of the present study are in agreement with those reported by Bellaloui et al. (2014).
Similar to the physiological responses to abiotic stresses, such as water deficiency, waterlogging, low phosphorus levels, cold temperatures, and light and nitrogen deficiency, a specific molecular mechanism also controls the responses of protein and oil content to increasing planting density (Osman et al., 2013). There are two aspects of the genetic basis of the effects of QTLs on the variation of protein and oil content at increasing planting density. First, the cumulative effects of QTLs could be detected directly based on protein and oil contents at a specific density in a particular environment. Here, these phenotypic values reflect the cumulative effects of genotype, macro-environment (location and years), interactions between genotype and environment, and planting density. The results of QTL mapping of the same trait can vary in different environments. Ku et al. (2015) evaluated the effects of two different planting densities (60,000 and 120,000 plant/hm 2 ) on three plant height-related traits (plant height, ear height, and ear height-to-plant height ratio) in maize. Nine QTLs were detected at the low planting density, and seven QTLs were detected at the high planting density. Akond et al. (2012) used a RIL population derived from a cross between soybean lines PI 438489B and Hamilton to detect QTLs for PC and OC under two planting densities. Three QTLs for protein content were detected  Tajuddin et al., 2003;Chen et al., 2007;Liang et al., 2010;Qi et al., 2011b;Eskandari et al., 2013;Mao et al., 2013 qRDOC at a low planting density (50 cm row space), while 2 QTLs for protein content and 6 QTLs for oil content were detected at a high planting density (25 cm row space). In the current study, we used an FW-RIL population to detect density-specific QTLs for protein and oil content. We identified 28 protein-content QTLs and 35 oil-content QTLs at planting densities D1 and D2, respectively. Among these QTLs, two QTLs was identified at both planting densities, whereas the remaining QTLs were detected at different densities. These results indicate that the genetic basis for the cumulative effects of QTLs on protein and oil content under two densities strongly differed, as indicated by significant genotype by density interaction effects (P ≤ 0.01; Table 2). Second, the net effects of QTLs on planting density could be identified based on increases in protein and oil content in response to density. To analyze the net effect of changes in planting density, the effects of all factors except planting density must be removed. Zhu (1995) proposed a conditional variable method to exclude the background from the covariance of related traits, allowing the net effects of some factors to be identified. This method has been used to estimate the net effects of various developmental stages  and correlated traits (Li et al., 2020). Here, we estimated the responses of protein and oil content to planting density using this method. We performed QTL mapping for the effect of increased planting density from 2.15 × 10 5 plants/ha (D1) to 3 × 10 5 plants/ha (D2) on protein and oil content. Using linkage analysis, we identified 20 and 18 RD QTLs controlling the responses of soybean protein and oil content to density, respectively, indicating that a specific molecular mechanism regulates the response to increasing planting density from D1 to D2. Therefore, when high-quality soybean varieties suitable for high-density planting were bred, alleles with positive effects should be pyramided. Conversely, when high-quality soybean varieties suitable for lowdensity planting was selected, alleles with negative effects should be combined.

Potential Candidate Genes Associated With Protein and Oil Contents
Based on the QTLs detected under different planting densities and the responses to density difference, as well as their pathway annotations, we identified four genes that might be related to differences in protein and oil contents under different plant densities ( Table 6; see bold text).
Energy produced in plants via photosynthesis is stored in the form of proteins, lipids, and other organic compounds. Glyma.10G215400.1 encodes pyruvate dehydrogenase E2 component, is mainly involved in carbon metabolism, citrate cycle (TCA cycle), and glycolysis/gluconeogenesis. This enzyme catalyzes the formation of pyruvate, which is the main substrate for the Calvin cycle, so we believe that the gene is affected by planting density, and closely related to protein and oil content. Glyma.19G190100.1 encodes an enzyme that regulates the formation of pyruvate kinase (PK) and plays an important role in carbon dioxide fixation and glycolysis in chloroplasts under light stimulation (Grodzinski et al., 1999). The enzyme participates in the glycolysis pathway, and metabolites can provide a premise for acetyl-CoA to form fatty acids and provide a carbon skeleton for fatty acid synthesis (Ambasht and Kayastha, 2002;Sébastien et al., 2007). Sucrose and starch produced by the glycolysis pathway can promote mitochondrial respiration and the TCA cycle, The TCA cycle is a key metabolic pathway for the oxidation of sugars, fats, and proteins (Robinson et al., 1991;Horchani et al., 2010). Therefore, we believe that this enzyme affects the protein and oil contents of soybean. Glyma.03G014600.1 encodes malate dehydrogenase (MDH), an important enzyme that catalyzes malic acid formation and is significantly affected by light stimulation, it is photo regulatory enzyme (Li et al., 1999). Malate is an important intermediate metabolite in plant cells, which are many biological functions in metabolic pathways (glyoxylate cycle, tricarboxylic acid cycle, glucose synthesis, amino acid synthesis, and redox stability), we know that these metabolic pathways are related to protein and oil synthesis (Minarik et al., 2002;Matsuda et al., 2010). Moreover, this enzyme is a key enzyme in the C4 pathway in Wheat (Hata and Matsuoka, 1987) and Soybean (Hao et al., 1991), which maintain a high photosynthetic carbon assimilation capacity when the light capacity and carbon dioxide content decrease. The activity of this enzyme is affected by oxygen concentration. Changes in planting density affect the oxygen concentration around plants, and changes in the expression of Glyma.03G014600.1 affect photosynthesis and the synthesis of proteins and oils. Therefore, changes in planting density cause competition between plants, which affect the ability of plants to absorb light, as well as oxygen. Therefore, we conclude that the gene responsible for this enzyme plays a role in regulating protein and oil contents in soybean. Glyma.19G190900.1 encodes an enzyme that catalyzes the enolase, which is involved in the changes in protein and oil content resulting from change of temperature. The activity of enolase continues to increase after temperature change, which strengthen the glycolysis process, thus promoting plant growth (Thomashow, 1999). Changes in planting density result in temperature in different plant varieties. Increasing the activity of enzymes that promote the glycolysis pathway process will enhance the protein and oil contents of plants. These findings could be used to enhance the protein and oil contents of soybean in the future.

Comparison of QTL Mapping Using LOD Value Is 3 and Permutation Method
LOD (log of odd) threshold is usually used to assure the existence of QTL. There are two methods to choose LOD threshold according to research goal for detect QTL. The first method is by permutation tests (Doerge and Churchill, 1996) which will generate a higher LOD threshold to decrease false positive QTL in each detection procedure. By this method a small amount of major effect QTL could be detected (Panthee et al., 2005;Bachlava et al., 2009;Tucker et al., 2010). The shortage of this method is consumption of time in the larger amount of calculation, and even failure in a large genome data. Sun et al. (2013) propose a method to certain the LOD threshold in QTL mapping, based on the genome-wide significance level, the population type, marker density and genome size, and Zhang P. et al. (2017) specified the threshold LOD value 3.766 for ICIM and IM of four-way cross recombinant inbred lines population to control the genome-wide typed one error at an equal level of 0.05. The second way is to set a lower LOD threshold (LOD = 2.0, equivalently P = 0.002 and significant at P < 0.01) with objective to avoid missing of QTLs due to slightly lower significance and the putative QTL were assured by repeated detected in multiple environments or multiple genetic background (Cornelious et al., 2005;Yesudas et al., 2013;Wang et al., 2014a,b). In this research, we firstly detected QTL depending on LOD threshold generated by permutation (1,000 repeats and 0.05 type one error probability), and 5, 6, 5 QTLs for PC, OC, RDPC with higher PVE (7.19%-38.02%) were detected. Among the 16 QTLs, only qOC-1-1 was repeatedly identified in E1 and E2.The major QTL could be detected directly by permutation which will miss minor QTL. The detection of minor QTL is necessary the inheritance and breeding of quality traits under multiple environments. The comparison on minor QTL detected under multiple environments and plant densities could explain the difference of genetic basis (molecular ecotypes). Furthermore, the accumulation of minor QTL could increase prediction of protein and oil content under some single planting density. Aimed to screen common and specific QTL under various densities, considering our linkage map characters (there are 1,031 marker intervals in 20 chromosome of 3539.66 cM length) ( Table S3), we choose a lower LOD threshold 3 (equivalently genome-wise P = 0.005) to screen QTL in the whole genome by ICIM. A total of 23, 31, 15 and 18 QTLs controlling PC, OC, RDPC, and RDOC with lower PVE (3.64-24.68%) were added. Then, five (qOC-1-1, qOC-5-1, qOC-10-1, qPC-18-2, and qRDPC-14-1) and two (qOC-15-1 and qOC-7-3) of all 103 QTLs were detected repeatedly in in two environment and two planting densities, respectively. Comparing two kinds of LOD threshold, five QTLs (qOC-10-1, qOC-15-1, qOC-5-1, qOC-7-3, qRDPC-14-1) could be verified the facticity by repeated identification under lower LOD threshold. In summarization, the permutation is suitable to detect the major QTL for higher PVE, and the repeated identification in multiple environments under lower threshold is feasible to discover QTL with low PVE.

SUMMARY
In this study, we identified 65 QTLs for protein and oil content under different planting densities and 38 QTLs with density responses based on SNP mapping. Based on these QTLs, four candidate genes were identified: these genes are affected by planting density and control protein and oil content. The molecular mechanism for the formation of protein and oil content under multiple planting densities involves the cumulative effects of QTLs and the response to increases in planting density. These findings lay the foundation for enhancing protein and oil contents and increasing yields in soybean at specific planting densities.

DATA AVAILABILITY STATEMENT
The datasets GENERATED for this study can be found in Figshare

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2020.00563/full#supplementary-material Table S1 | Comparison of parental source and quality traits.