For linkage mapping, two populations were used in this study. A recombinant inbred line (RIL) population was used for the initial QTL detection. The RIL population (NJMN-RIL) consisted of 104 F_2:8 lines developed from a cross between M8108 (a landrace material, seed oil content exceeding 20%) and NN1138-2 (an elite cultivar characterized by high yield, seed oil content ∼ 19.0%). NJMN-RIL and its two parents were grown in five environments: Jiangpu Experimental Station, Nanjing, Jiangsu Province, in 2012 and 2014 (12JP and 14JP); Fengyang Experimental Station, Chuzhou, Anhui Province, in 2012 (12FY); Yancheng Experimental Station, Yangcheng, Jiangsu Province, in 2014 (14YC); and Huaian Experimental Station, Huaian, Jiangsu Province, in 2014 (14HA).

A secondary F₂ population was used to confirm the most stable QTL from the initial mapping results. The F₂ population consisted of 203 individuals developed from a backcross using MN-5 with NN1138-2. MN-5 was a line from the NJMN-RIL population with the M8108 type–homozygous qOil-5-1 region, while the NN1138-2 type–homozygous was the other major QTL region identified in this study. The F₂ population was grown at Jiangpu Experimental Station in 2015.

A total of 279 soybean accessions from the Yangtze-Huai soybean breeding germplasm population were used to perform the association study. The genotype of these soybean accessions can be obtained from our previous study (Li et al., 2016). The soybean accessions were grown in two environments: 13JP and 14JP.

The NJMN-RIL and germplasm populations were grown in a randomized complete block design, with three replications from approximately June to October. A 1 m single row plot of 10 plants per RIL/accession was grown, with a distance of 10 cm between plants and a row spacing of 50 cm. The secondary F₂ population was grown in rows with a length of 2 m, with 50 cm between rows and 20 cm between plants. Field management was performed under normal conditions.

Seed oil content was measured by an Infratec^TM1241 near infrared analysis (NIR) Grain Analyzer (Foss, Hillerød, Denmark) on a 10% moisture basis using approximately 15–20 g samples. Seed oil content for each RIL and soybean accession was an average of three replications within a planting environment.

The frequency distribution of phenotypic data was calculated using the SPSS Statistics 20.0 software (SPSS Inc., Chicago, IL, United States). The broad-sense heritability (h²) and analysis of variance (ANOVA) tests were conducted using the SAS PROC generalized linear model (GLM) program. The h² of seed oil content in RIL and germplasm population was estimated using the following equation:

where σ²_g is the genotypic variance, σ²_ge is the variance of the genotype-by-environment interaction, σ²_e is the error variance, n is the number of environments, and r is the number of replications within an environment (Nyquist and Baker, 1991).

For the NJMN-RIL population, approximately 1 g of fresh leaves obtained from two parents and individuals of each RIL were used to extract the genomic DNA using the cetyltrimethylammonium bromide method (Doyle, 1990). All RIL individuals and the two parents were genotyped using SLAF-seq technology to generate genome-wide SLAF markers. SLAF library construction, high-throughput sequencing, and high-quality SLAF marker acquisition were performed as described by Sun et al. (2013). The soybean reference genome (G. max Wm82.a1) was used to determine the physical position of the high-quality SLAF markers.

For the secondary F₂ population, a total of 31 pairs of simple sequence repeat (SSR) markers based on the soybean SSR database (Song et al., 2010) and four pairs of InDel markers based on the results of Lam et al. (2010) and distributed across the qOil-5-1 region (38–40.5 Mb on the Chr. 05), were used to screen for polymorphisms between MN-5 and NN1138-2. Finally, 10 polymorphic pairs (one InDel marker and nine SSR markers) were used to genotype the F₂ population and then construct a genetic map covering the qOil-5-1 region.

A set of 279 soybean germplasm lines was assembled for GWAS. All lines were genotyped for SNP markers using high-throughput genotyping platforms, as described by Li et al. (2016). Finally, 59,845 SNPs, with minor allele frequencies (MAF) > 0.05, were available for the GWAS.

After genotyping the 104 RILs, all polymorphic SLAF markers were filtered four times and further quality assessed as described by Sun et al. (2013). A SLAF with less than three SNPs and an average depth for each sample above three was considered a high-quality SLAF marker. Then, all high-quality SLAF markers were used for construction of the linkage map. The linkage map was constructed using High Map software (Liu et al., 2014). The markers were divided into linkage groups (LGs) using the single-linkage clustering algorithm at a logarithm of odds (LOD) threshold ≥3.0. The Kosambi mapping function (Kosambi, 1943) was used to calculate the map distances in cM from the recombination frequencies.

The QTLs were detected using WinQTLCart 2.5 software. The window size, the working speed, and the control marker number were set at 10 cM, 1 cM and 5, respectively. Model 6 (standard model) for composite interval mapping was used to identify QTLs in each environment, with a log likelihood of 2.5 set as the threshold of an existing QTL (Wang et al., 2007).

The Q and the kinship matrix (K) for these 279 soybean accessions were calculated as previously described (Li et al., 2016). The 279 soybean materials could be divided into three major subpopulations. The LD was approximately 480 kb, whereas r² dropped to half its maximum value (Li et al., 2016). The LD blocks were within the interval of a QTL based on an LD analysis using the r² correlation between each marker using the HAPLOVIEW 4.2 software (Barrett et al., 2005). The parameters in the program included MAF (≥0.05) and the integrity of each SNP (≥50%).

In this study, the GWAS was performed using the TASSEL 5.0 software (Bradbury et al., 2007). Two statistical models were considered: (1) the GLM model without the Q and K; and (2) the mixed linear model (MLM) with the Q and K, which regarded Q and K as fixed and random effects. Compared with the GLM method, the MLM model can reduce the number of false positive results and was considered more appropriate for this study. To obtain more reliable results, the Bonferroni threshold (P < 1/59845, -log₁₀P > 4.78) was used as the threshold to identify significant association SNPs. The association region was also estimated by extending 480 kb (LD distance) upstream and downstream of the most significantly associated SNP position.

The phenotypic performance is presented in Table 1 and Supplementary Figure S1. For the NJMN-RIL population, in the female parent, M8108 had an average seed oil content of ∼20.5%, and NN1138-2 of ∼19%, with these differences being significant in all environments (P < 0.001). A continuous distribution and transgressive segregation were observed for the trait in the RIL population in all environments. Due to the kurtosis and skewness (absolute value) being less than 1, the phenotypic frequency distributions had a relatively normal distribution. This indicates that the phenotypic data was suitable for QTL mapping. For the germplasm population, the seed oil content ranged from 15.63 to 22.90% and from 14.67 to 22.63% in 2013 and 2014, respectively. The phenotypic frequency distributions of the germplasm population also had a normal distribution.

In addition, the ANOVA results indicated significant differences in genotype, environment, and genotype-by-environment (Supplementary Table S1). The h² of seed oil content in the NJMN-RIL population and germplasm population was high (0.93 and 0.88, respectively), which was consistent with previous studies (Chung et al., 2003; Panthee et al., 2005).

The RIL population was genotyped using SLAF-seq technology. After being filtered and quality assessed, a total of 2,086 SLAF markers were available to construct genetic maps for the NJMN-RIL population. Finally, 2,062 SLAF markers were grouped into 20 LGs. The total genetic distance of this map was 2054.50 cM. The average distance between adjacent markers was approximately 1 cM. The genetic length of 20 LGs ranged from 40.21 cM (Chr. 17) to 199.77 cM (Chr. 12). The largest LG was Chr. 18, with 338 SLAF markers. The LG with the minimum number of markers was Chr. 05, with 28 SLAF markers. The mean LG length was 102.73 cM and each chromosome contained an average of 103 markers. Detailed map information is provided in Figure 1 and Supplementary Table S2.

Based on the genetic map, a total of eight QTLs for seed oil content, distributed on seven chromosomes, were identified in the NJMN-RIL population across the different environments by WinQTLCart 2.5 software, with LOD scores ranging from 2.7 to 9.1 and 6.3 to 26.2% of the phenotypic variation explained by individual QTLs (Table 2). Among them, qOil-5-1 was the most stable QTL detected in all environments and could explain 7.5–26.2% of the phenotypic variation. The qOil-5-1 QTL was mapped at an approximate 6.8 cM (43.4–50.2 cM) interval in all environments (Figure 2A and Table 2). This QTL interval has a corresponding physical distance of 2.5 Mb, from 38.04 to 40.51 Mb on Chr. 05 (Figure 2B). The qOil-10-1 and qOil-14-1 QTLs could also be considered major QTLs because they were detected in four environments. The other five QTLs were detected in a single environment. In addition, the positive alleles of qOil-3-1 and qOil-10-1 came from NN1138-2, while the others came from M8108.

The qOil-5-1 QTL was detected in all environments and so could be considered the most stable QTL in the NJMN-RIL population. To confirm and refine the position of this QTL, the secondary F₂ population was used. The seed oil content ranged from 17.50 to 21.90% in this F₂ population (Table 1 and Supplementary Figure S1). We also constructed a genetic map with 10 polymorphic pairs for the F₂ population (Figure 2C). This genetic map covered the qOil-5-1 region (from ∼37.8 to 40.8 Mb on Chr. 05). Based on the phenotype and genetic map, a QTL was identified. The highest peak of the QTL was found at the marker InDel-5-7 [primer sequence: F′ (5′–3′): CAAAAATACATTGAATGGAAGGCAA; R′ (5′–3′): TGACAGAGAGATATTATCCCCAGC], with an LOD score of 6.87, which strongly exceeded the threshold level. The InDel-5-7 marker corresponded to position 38494541-38494952 bp on the reference genome (G. max Wm82.a1). The 1-LOD confidence interval of the QTL spanned approximately 1 Mb, from position 37.85 Mb (BARCSOYSSR_05_1217) to 38.90 Mb (BARCSOYSSR_05_1261) of the physical map, which was consistent with the RIL population mapping results (Figures 2B,C).

In this study, 279 soybean gremplasm lines with 59,845 (MAF > 0.05) SNPs were used to perform the association analysis. The quantile–quantile (QQ) plot showed that the MLM method was appropriate for this study and could reduce the number of false positive results (Supplementary Figure S2). Therefore, the MLM model was used for the association analysis.

There were two association peaks identified on Chr. 05 and on Chr. 20 in 2013; and one association peak on Chr.05 and two association peaks on Chr. 20 were identified in 2014. The phenotypic variation explained by each peak ranged from 8.5 to 12.9% (Table 3 and Figure 3A).

The loci for seed oil content on Chr. 05 could be detected in both the linkage and association population. Therefore, qOil-5-1 can be considered an important locus for seed oil content in the summer planting soybeans. To reduce the qOil-5-1 region, we focused on the results from the association population (Figure 3B). The peaks on Chr. 05 were located at 38506373 and 38473956 bp in 2013 and 2014, respectively. Based on the LD distance, we extended the region of 480 kb upstream and downstream of the significantly associated SNP position. The LD was estimated from the r² correlation between each marker in this region, with r² values less than 0.2 considered unlinked. This region possessed two large LD blocks, one from 38.13 to 38.57 Mb and the other from 38.59 to 38.95 Mb (Figure 3C). However, all the SNPs showing a significant association with seed oil content in the 2 years were located within the first LD block region. This suggested that the major QTL for seed oil content was located in a small region from 38.13 to 38.57 Mb on Chr. 05.

Quantitative trait loci mapping is an effective approach to analyze quantitative traits in plants. Linkage mapping and GWAS are commonly used to map the QTLs of complex quantitative traits. However, the linkage mapping method has a relatively low genome resolution and the GWAS results are affected by the Q and individual relationships. A combination of GWAS and linkage mapping could provide mutual authentication and enable more reliable results to be obtained. This strategy has been successfully applied in many studies (Konishi et al., 2006; Zhang D. et al., 2014; Zhang Y. et al., 2016; Zhao et al., 2015; Sun et al., 2016). Because the soybean cultivars in the Jiang-Huai River Valley usually have a relatively low seed oil content (Zhang J. W. et al., 2014), it is difficult to breed high oil content varieties by traditional breeding in this region. Molecular MAS is a possible alternative. However, little is known about the genetics of oil content in these soybean. In this study, we used 104 RILs and 279 accessions to dissect the genetic basis of seed oil content in soybean, which would provide useful information for MAS to breed high oil content varieties in this region.

The marker density should be considered during QTL mapping, regardless of which method is used. Because of the technology limitations, the genetic maps used in most of the studies mentioned above were constructed for mapping QTLs associated with seed oil content in soybean with only hundreds of restriction fragment length polymorphism, amplified fragment length polymorphism, and/or SSR markers. In addition, the interval distance between adjacent markers was large, which limits the efficiency and accuracy of QTL mapping. In our study, SLAF-seq was used to genotype a total of 104 individuals and the two parents. We constructed a genetic map with 2,062 SLAF markers for linkage mapping, and the average distance between adjacent markers was only 1 cM. For GWAS, 59,845 polymorphism SNPs with MAF >0.05 were used. Previous studies have shown that the LD of cultivated soybean is approximately 150 kb (Lam et al., 2010; Hao et al., 2012; Wen et al., 2014). Hwang et al. (2014) reported that the mean LD (r²) rapidly declined to 0.2 within 360 kb in euchromatic regions and declined to 0.2 at 9600 kb in heterochromatic regions. However, there is a high LD in the soybean genome, implying a limitation on the effectiveness of increasing marker density to improve the resolution. This genetic map and high-number of SNPs can meet the demands of linkage and a GWAS of the complex traits in soybean.

Soybean is a primary source of plant oil for humans (Wilcox, 2004). The breeding of high oil content varieties is an important goal of oil crop breeding programs. The QTLs controlling seed oil content have been reported many times in bi-parental populations (Mansur et al., 1996; Brummer et al., 1997; Kabelka et al., 2004; Qi et al., 2011; Wang et al., 2015). Many loci associated with soybean seed oil content have been detected by GWAS methods (Hwang et al., 2014; Bandillo et al., 2015; Wen et al., 2015). However, most QTLs are environmentally sensitive. Brummer et al. (1997) considered that in molecular-assisted breeding programs, breeders should use QTLs that are stable in a range of different environments. Thus, there is a need to identify stable QTLs that can be conveniently used for MAS and map-based cloning to determine the mechanisms of seed oil content.

In this study, the distribution of seed oil content followed a normal distribution in both the natural soybean population and in RILs, indicating that the genetic control of seed oil content is complex (Supplementary Figure S1). Although seed oil content had a high h², the environmental instability of this trait made it difficult to detect stable QTLs in different environments. However, there were some major loci detected in different environments using one method in our study. In linkage mapping, qOil-10-1 and qOil-14-1 were detected in four environments, had adjacent or physically overlapping QTLs, as reported by previous studies, and might actually be the same QTL (Tajuddin et al., 2003; Wang et al., 2014). In the GWAS, the significant SNPs on Chr. 20 were found to be associated over the 2 years investigated. The peak marker, Gm20_31164168, was similar to the one identified by Vaughn et al. (2014) at 31150279 bp on Chr. 20. Therefore, our results were in good agreement with previous studies, and the QTLs identified can be treated as important targets for the identification of candidate genes involved in the modification of oil in future studies.

Although the number of QTLs controlling seed oil content have been reported by linkage mapping and/or GWASs in soybean, the confirmation of such QTLs is rare. Here, we confirmed and narrowed down the interval of a stable QTL on Chr. 05 by a combination of linkage and association mapping. In the NJMN-RIL population, qOil-5-1 was mapped at an approximate 6.8 cM interval in all environments. Based on the physical position of the SLAF markers, this genetic map interval has a corresponding physical distance of 2.5 Mb, from 38.04 to 40.51 Mb on Chr. 05, which has adjacent or physically overlapping QTLs, as reported by Mansur et al. (1996), Orf et al. (1999), and Pathan et al. (2013). In the secondary F₂ population, the position of this QTL could be refined at a 3 cM interval between the BARCSOYSSR_05_1226 and BARCSOYSSR_05_1250 markers, and the peak marker corresponded to position 38494541–38494952 bp. These results indicated that the region of qOil-5-1 is from ∼37.85 to 38.90 Mb. In the germplasm population, based on the GWAS results and the partial LD of Chr. 05, we found that all the SNPs that had a significant association with seed oil content over the 2 years investigated were located in one LD block. Taken together, this evidence strongly suggests that this stable QTL for seed oil content was located in a small region from 38.13 to 38.57 Mb on Chr. 05.

Plant oil is typically biosynthesized and accumulated in the cotyledon or endosperm tissues of seeds. The pathways for lipid biosynthesis and oil accumulation have been studied and are known to be controlled by multiple genes. Within the LD block region, there were more than 60 functional genes². Thirteen of those genes were highly expressed in seed (information from Phytozome database³) and may be related to oil biosynthesis or oil accumulation (Table 4). However, the functions of these genes need to be verified in future studies. Finally, such fine-mapping of QTLs and close markers (including SNP, SSR and InDel markers) must be used for MAS and map-based cloning to determine the mechanisms of seed oil content.