Construction of a high-density genetic map and mapping of QTLs for soybean (Glycine max) agronomic and seed quality traits by specific length amplified fragment sequencing

Background Soybean is not only an important oil crop, but also an important source of edible protein and industrial raw material. Yield-traits and quality-traits are increasingly attracting the attention of breeders. Therefore, fine mapping the QTLs associated with yield-traits and quality-traits of soybean would be helpful for soybean breeders. In the present study, a high-density linkage map was constructed to identify the QTLs for the yield-traits and quality-traits, using specific length amplified fragment sequencing (SLAF-seq). Results SLAF-seq was performed to screen SLAF markers with 149 F8:11 individuals from a cross between a semi wild soybean, ‘Huapidou’, and a cultivated soybean, ‘Qihuang26’, which generated 400.91 M paired-end reads. In total, 53,132 polymorphic SLAF markers were obtained. The genetic linkage map was constructed by 5111 SLAF markers with segregation type of aa×bb. The final map, containing 20 linkage groups (LGs), was 2909.46 cM in length with an average distance of 0.57 cM between adjacent markers. The average coverage for each SLAF marker on the map was 81.26-fold in the male parent, 45.79-fold in the female parent, and 19.84-fold average in each F8:11 individual. According to the high-density map, 35 QTLs for plant height (PH), 100-seeds weight (SW), oil content in seeds (Oil) and protein content in seeds (Protein) were found to be distributed on 17 chromosomes, and 14 novel QTLs were identified for the first time. The physical distance of 11 QTLs was shorter than 100 Kb, suggesting a direct opportunity to find candidate genes. Furthermore, three pairs of epistatic QTLs associated with Protein involving 6 loci on 5 chromosomes were identified. Moreover, 13, 14, 7 and 9 genes, which showed tissue-specific expression patterns, might be associated with PH, SW, Oil and Protein, respectively. Conclusions With SLAF-sequencing, some novel QTLs and important QTLs for both yield-related and quality traits were identified based on a new, high-density linkage map. Moreover, 43 genes with tissue-specific expression patterns were regarded as potential genes in further study. Our findings might be beneficial to molecular marker-assisted breeding, and could provide detailed information for accurate QTL localization. Electronic supplementary material The online version of this article (10.1186/s12864-018-5035-9) contains supplementary material, which is available to authorized users.


Background
Soybean is not only an important oil crop, but also an important source of edible protein and industrial raw material [1]. Agronomic traits, such as yield, plant height (PH), lodging and seed weight (SW), have been the primary focus of breeders for many years. As the major factors of the market price of soybean, seed quality traits are increasingly attracting the attention of breeders. However, the negative correlation between yield and quality of crops makes it much difficult to select for these traits [2]. Therefore, simultaneous improvement of yield and quality has become a major problem for soybean breeders.
Molecular marker-assisted selection (MAS) might be an alternative to fit the increasing global demand for soybean products [3]. A number of QTLs underlying important agronomic traits and seed quality traits have been reported over the past decades. So far, at least 196, 265, 297 and 221 QTLs controlling PH, SW, Oil and Protein have been identified respectively (www.soybase.org), based on the different genetic backgrounds, environments and statistical methods. Furthermore, large confidence intervals around QTLs make the causative gene identification difficult.
With the development of next generation sequencing technology, several methods for single nucleotide polymorphisms (SNP) discovery, such as restriction-site associated DNA sequencing (RADseq) [4,5], genotyping-by-sequenci ng (GBS) [6], specific length amplified fragment sequencing (SLAF-seq) [7] have been produced, which make it possible to obtain thousands of SNPs suitable for high-density genetic map throughout the genome. SLAF markers, which have the properties of being present in large amount, being evenly distributed and avoiding repeated sequences, has been used for genetic analysis in plants, such as sesame [8], walnut [9], rice [10], sorghum [11], wax gourd [12], grape [13] and soybean [14][15][16][17]. Since the first high-density map was constructed by SLAF-seq [7], there have been several maps reported so far. Qi et al. constructed a map, including 5308 markers with 2655.68 cM in length, using a RIL population derived from a cross between 'Charleston' and 'Dong-nong594' [16]. Li et al. constructed a high-density map, using a F5:8 population of 110 RILs from a cross between 'Luheidou2' and 'Nanhuizao' , which was used to identified QTLs associated the isoflavone content and fatty acid composition in soybean [15,18]. Zhang et al. reported 20 QTLs associated with phosphorus efficiency-related traits based on a high-density map constructed by SLAF-seq [17]. Cao et al. mapped QTL associated with plant height and flowering time according the map constructed by SLAF-seq using a population of 236 RILs derived from a cross between two summer planting varieties, 'ZXD' and 'NN1138-2' [14]. Nevertheless, based on the high-density map, little QTLs related to seed weight and/or protein have been reported. Therefore, we reported a high-density genetic linkage map using the SLAF-seq approach, which was based on an F 8:11 RIL population with 149 individuals. Moreover, the QTLs associated with plant height, seed weight, oil and protein content were located and analyzed. The results presented here will aid molecular marker-assisted breeding and provide detailed information for accurate QTL localization.

Analysis of SLAF-seq and SLAF markers
DNA sequencing generated about 400.91 M pair-end reads. The Q30 (indicating a 0.1% chance of error) was 90.69% and guanine-cytosine (GC) content was 40.3%. The numbers of SLAFs in the female and male parents were 312,740 and 275,046, respectively. The numbers of SLAFs in each individual ranged from 167,933 to 237,666, with an average of 207,105. Among the 391,476 SLAF markers detected, 53,132 markers were polymorphic. All polymorphic SLAFs were then genotyped separately for all individuals. After discarding the SLAF markers lacking parent information, 30,415 markers were genotyped successfully and were classified into eight segregation types (Fig. 1). As the population was derived from a cross between two fully homozygous parents, only 27,472 markers with aa×bb type might be suitable for map construction. After filtering low-quality SLAF markers, segregation distortion markers and makers with the MLOD value ≤3, 5111 markers were used for the map construction (Additional files 1, 2, and 3). The average depth of the markers was 81.26-fold in the female parent, 45.79-fold in the male parent, and 19.84-fold in the offspring.

The basic characteristics of the genetic map
The length of final map was 2909.46 cM, with an average distance of 0.57 cM between adjacent markers (Table 1; Fig. 2). There were 8597 SNP loci among the 5111 markers on the map. For each chromosome, the average distance ranged from 0.24 cM to 2.55 cM ( Table 1). The largest linkage group was LG18 (chr18) with 480 markers, a length of 202.52 cM, and an average distance of only 0.42 cM between adjacent markers. The smallest linkage group was LG7 (chr7) with 63 markers, a length of 54.57 cM, and an average distance of 0.87 cM between adjacent markers.

Phenotypic evaluation
'Huapidou' and 'Qihuang26' showed a significant difference in PH, SW and Protein, but did not differ from each other in Oil significantly (Table 2; Fig. 3). However, the phenotypic values were all in a condition of continuous distribution approximately (Fig. 3). The coefficients of variation for four traits were about 20%. The heritabilities of four traits ranged from 49.58 to 82.73%. However, the heritability of Protein was only 49.58%, indicating that other factors affected Protein should be considered.

Analyses of additive QTLs
In total, 35 additive QTLs for PH, SW, Oil and Protein were identified on 17 chromosomes by ICIM method (Table 3; Figs. 2, 4). A single QTL explained 2.66% (qPH16-1) to 37.61% (qPH8-1) of phenotypic variance. Among the QTLs, 14 QTLs of them were observed for the first time (Table 3). A total of 21 QTLs were related to the region of the QTLs reported previously, and 19 of them were co-located in the regions with shorter intervals than previously reported, which might provide more detailed information for gene identification.
CIM method was also used to identify QTLs separately for 2 years. In total, 18 QTLs were observed on 11 chromosomes in 2013, and 21 QTLs were found on 15 chromosomes in 2014 (Additional file 4; Fig. 2). Of these QTLs, 31 QTLs were both identified by two methodsMoreover,  q2013PH19-1 and q2014PH19-2 with shorter intervals than previously reported [19], which were located in the same confidence intervals, might be stable across both years. It was noteworthy that q2013Oil1-1 was placed in the same confidence intervals as q2013Protein1-1, and q2014Oil10-2 was placed in the same confidence intervals as q2013protein10-1, which might be useful in the coordinated improvement of seed quality for soybean breeding.

Analyses of epistatic effects
A total of 3 pairs of epistatic QTLs involving 6 loci on 5 chromosomes were identified for Protein ( Table 4). The epistatic effect explained 5.49%, 4.49% and 4.06% of the PV, respectively. Pair one was composed of 2 QTLs, qProtein2-1 and qProtein3-2, with the PVE of 5.49%. qProtein5-2 was observed to have the epistatic effect with qProtein12-1. Meanwhile, qProtein3-3 showed epistatic interaction with qProtein17-1. However, no epistatic effect was observed for PH, SW and Oil.

Prediction of candidate genes
After filtering QTLs by the PVE and physical distance, 18 QTLs were used to mine candidate genes. According to the physical map, a total of 89, 144, 16 and 64 genes were screened in the interval of the filtered QTLs associated with PH, SW, Oil and Protein, respectively. Based on the expression data of candidate genes from phytozome and soybase ( Fig. 5), 43 genes were considered to be potential candidates ( Table 5). All genes from specific QTLs intervals were evaluated based on their expression pattern in different organs. In the case of genes from the QTLs associated with PH, 13 genes showed higher expression in stem and shoot apical meristem, indicating it might be considered as candidate genes related to PH (Fig. 5a). A total of 14 genes in the interval of the QTLs  associated with SW, expressed in seed development stages (10 to 42 days after flowering), might participate in the pathways affecting SW ( Fig. 5b and 5c). As the accumulation of oil and protein in soybean seed was throughout the seed development stage [20,21], gene expressed sustainably in seed development stage might affect the biological process associated with oil and protein. In the present study, there were 7 genes found in the region of QTLs related to Oil, which expressed stably throughout the seed development stage, suggesting it might be associated with Oil (Fig. 5d). Meanwhile, 9 genes in the interval of QTLs for Protein, showed expression in the seed development stage, might be associated with Protein ( Fig. 5e).

Construction of a high-density genetic map based on SLAF markers
QTL mapping has been used as an efficient approach to analyze quantitative traits in plants. Parental genetic diversity and marker density are the major factors affecting the efficiency and accuracy of QTL mapping. In this study, the female parent, 'Huapidou' , was a semi-wild soybean gerplasm, which showed high resistance to whitefly [22]. 'Qihuang26' , with more than 46% of protein conten in seeds, was a main variety in Huang-Huai-Hai region of China. In the present study, four traits of the RIL population derived from Huapidou and Qihuang26 showed to be continuous with normal or skew normal distributions. Increasing marker density could improve the resolution of genetic map for a given mapping population [23]. SLAF-seq is an effective sequencing-based method for large-scale marker discovery and genotyping, which has been used for genetic analysis in many species [8][9][10][11][12][13]24]. In the present study, we used 5111 high-quality SLAF markers to construct a high-density map, and a total of 8597 SNP loci were integrated into 20 LGs ultimately. This high-density genetic map, making QTL mapping more accurate and reliable, would be beneficial to MAS breeding.

QTL mapping in soybean using a high-density map
Soybean is a primary source of plant oil and protein for humans due to its high nutritional value. PH and SW were main yield-related traits in soybean. So far, markers associated with the QTL underlying PH, SW, Oil and Protein have been mapped onto all linkage groups. In total, there were 35 QTLs for PH, SW, Oil and Protein observed using a high-density map based on an F8:11 RIL population with 149 individuals from the cross between 'Huapidou' and 'Qihuang26'. Furthermore, there were 14 novel QTLs related to PH, SW, Oil and Protein, indicating the distinct genetic architecture in the population derived from cultivated soybean and semi-wild soybean. Among the novel QTLs, qPH8-1 had the highest PVE value and the highest LOD value might be the major QTL related to PH. It was notable that qSW13-1 explained the hightest PV in the QTLs identified for SW. More remarkably, four novel QTLs for Oil, inculding qOil1-1, qOil1-2, qOil10-1 and qOil10-2 explained up to 72.73% of the PV for Oil, which suggested it might be potential loci to Oil. qProtein1-1, which explained 17.68% of the PV, might be an major QTL for further fine mapping. So many novel QTLs observed in the present study indicated that more germplasms need to be used for revealing the complex genetic basis of soybean.
The stability of QTL is essential for the use in a breeding programme. In the study, 31 QTLs were identified by both ICIM and CIM methods. Furthermore, one QTL for PH was identified by CIM in both experiments    [14,19,25,26]. Two major QTLs associated with SW, qSW13-1 and qSW15-1, both with the physical distance of approximate 11 Kb, explained 18.55% and 12.21 of the PV for SW, respectively. qSW13-1 had been reported as being associated with L050-14 [27], Satt144 [28,29] and Sat_103 [30]. qSW15-1 had been detected in two soybean populations, derived from 'Young' and 'PI416937' (Pop1), 'PI97100' and 'Coker 237' (Pop2) [31]. Han et al. also identified the similar QTL on chr15 in the population from a cross between 'Hefeng25' and 'Conrad' [32]. Therefore, qSW13-1 and qSW15-1 might be considered as major and stable QTLs for further fine mapping and map-based cloning to elucidate the mechanisms of SW.
In the present study, four QTLs related to Oil had been reported [28,[33][34][35], inculding qOil6-1, qOil13-1, qOil19-1 and qOil19-2, but none of them explained more than 10% of the PV. Lee et al. [36] reported cr274_1 associated with Protein on chr15 using a population derived from 'Young' and 'PI416937'. The QTL for Protein between Satt173 and Satt581 on chr10 had been identified [37], similar with the result of Liu et al. [38]. Our study detected two QTLs related to Protein, qProtein10-1 and qProtein15-1, with 16.83% and 14.36% of the PVE, respectively, mapped on the same area as previous studies [36][37][38], might be good for MAS breeding and accurate QTL localization. Several QTLs of various traits can map to the same locus [14,39]. In this study, two pairs of QTLs, q2013Oil1-1 and q2013Protein1-1 as well as q2014Oil10-2 and q2013pro-tein10-1, with inverse additive effect for Oil and Protein, were located in the same marker interval ( Fig. 4; Additional file 4), which implies that q2013Oil1-1 and q2014Oil10-2 not only control oil content in seeds but also affect protein content in seeds. It is consistent with previous reports that an negative correlation is in agreement between protein and oil concentration in soybean seeds [40,41].
Knowledge of epistasis effect, which is defined as interactions between alleles of two or more genetic loci, is essential to understand the genetic mechanism and the  gene networks underlying complex traits. In this study, 3 pairs of epistatic QTLs for Protein were identified by ICIMapping-EPI. However, these epistatic QTLs did not display additive effect alone. It might be considered modifying genes that have no significant effects alone but might affect the expression of Protein related genes through epistatic interactions. Nevertheless, epistatic interaction could not be detected in some map populations [42]. It might be the reason that no epistatic effect observed for PH, SW and Oil in the present study.

Gene mining based on precise QTLs
As the average ratio of gene to physical distance is about 1 gene per 20 Kb in soybean genome [43], the accuracy of QTL mapping is of great benifit to gene localization and identification. The physical distance of 11 QTLs in the current study was shorter than 100 Kb, which might lead to a direct opportunity to find candidate genes by bioinformatics tools. For example, the minimum confidence interval of qPH8-2 was 8.3 Kb, which was much shorter than 0.09 Mb detected previously [26]. Furthermore, Glyma.08 g337400, encoding a transducin/WD40 repeat-like superfamily protein, was predicted in the interval of qPH8-2, which might be a promising target to engineer transgenic plants with higher biomass and improved growth development for plant-based bioenergy production [44]. In the interval of qPH17-1, Glyma.17 g169100, encoding a 2OG-Fe(II) oxygenase superfamily protein, was one of the important gibberellin oxidase genes [45], which might affect plant height directly. Glyma.17 g167700, encoding a growth regulator CYCLIN D3-2, expressed in growing shoot apices preferentially [46][47][48]. In the interval of qSW5-1, Glyma.05 g055700, encoding beta vacuolar processing enzyme, was involved in seed coat formation at the early stage of seed development [49]. In the interval of qSW12-1, ribosomal protein L23AB encoded by Glyma.12 g066700, was required for normal development [50]. Glyma.01 g106000 and Glyma.01 g106100 in the interval of qOil1-2, encoding Glutathione  S-transferase TAU 8, might influence Oil by suppressing lipid peroxidation [51]. There were 3 genes found in the confidence interval of qOil10-2, which explained 20.17% of the PV for oil, inculding Glyma.10 g181800, Glyma.10 g181900 and Glyma.10 g182000. Moreover, Glyma.10 g181900, encoding a trigalactosyldiacylglycerol 1 protein (TGD1), affected the metabolic flux of chloroplast lipid synthesis and photosynthetic capacity, which resulted in the change of fatty acid in leaf and seed [52][53][54][55]. It was noteworthy that qOil10-2 was placed in the same confidence intervals as q2013protein10-1. The inverse relationship between oil and protein in soybean seed is well documented in the previous reports [40,41]. However, little study of TGD1 had been reported on the function of protein accumulation in seeds. In the interval of qProtein10-1, Glyma.10 g183900, encoding peptide transporter 3, contributed to nitrogen allocation and grain Yield [56], Glyma.10 g184900, encoding a ureidoglycolate amidohydrolase, played a key role in nitrogen transport and storage [57][58][59]. In a word, on the basis of the physical position of these precise QTLs detected using a high-density map in the present study, it would be easy to find candidate gene.

Conclusions
In this study, we genotyped a RIL population (Huapidou × Qihuang26) by SLAF-seq. A high-density genetic map for soybean was constructed and used to identify QTLs associated with four traits, including plant height, seed weight, oil content seed and protein content in seed. A total of 35 QTLs related to four traits were identified. Of these QTLs, 21 QTLs were coincident with previous research. Furthermore, three pairs epistatic QTLs involving 6 loci on 5 chromosomes were identified for Protein. In addition, 43 genes with tissue-specific expression patterns were considered to be potential genes in further study. Our findings might be of great useful for MAS breeding, and could provide detailed information for accurate QTL localization.

Methods
Plant material and phenotyping F8:9 and F8:10 populations of 149 RILs derived from a cross between 'Huapidou (ZDD09982)' and 'Qihuang26 (ZDD23189)' were planted in the experiment field of Shandong Academy of Agricultural Sciences in Jinan, Shandong Province, China, in 2013 and 2014, respectively. Each individual was planted in one row using single seed sowing; each row was 3 m, with 50 cm row spacing and 10 cm plant spacing, with three replicates. Five plants in each replicate were selected randomly to calculate the plant height (PH). The weight of 100 random filled seeds was measured as seed weight (SW). Oil and protein in soybean seed were detected by DA 7200 NIR food analyzer (Perten, Switzerland). SW, Oil and Protein were repeat 3 times in each replicate. Frequency distribution, descriptive statistics, the broad-sense heritability (h 2 ) and the analysis of variance for RIL population and parents were analyzed with the SPSS statistics 17.0 and Microsoft Excel 2010. The h 2 was estimated as described by previous study [60].

DNA extraction and genotyping
Seedlings of the F8:11 population of 149 RILs and parents were planted in 2016. Young healthy leaves from the two parents and RIL individuals were collected and genomic DNA was extracted by the CTAB method [61]. DNA was quantified with NanoDrop and by electrophoresis in 1% agarose gels with a λ DNA standard. SLAF-seq was used to genotype a total of 151 samples (149 individuals and two parents) as described by previous study [7]. All polymorphic SLAFs were genotyped with consistency in the offspring and parental SNP loci. All SLAF markers should be filtered in quality assessment. A SLAF marker with parental homozygous, which had less than three SNPs, average depths of each sample above 3, was used as a high quality SLAF marker.

Linkage map construction
Before map construction, SLAF marker should be filtered by linkage analysis, markers with the MLOD value > 3 were used to construct genetic linkage map. SLAF markers with high quality were located into 20 LGs. HighMap Strategy was used to order SLAF markers and correct genotyping errors [24]. All LGs should be undergone these procedures: first, markers were arranged by their locations on choromosome; second, genotyping errors or deletions were corrected by SMOOTH [62], according to the relationship between ordered markers; then MSTmap was used to order the map [63]; after that, SMOOTH was used again to corrected the new ordered genotype. High-quality map would be obtained after 4 or more cycles. Map distance was estimated using the Kosambi mapping function.

QTL mapping
Based on the high-density genetic map, the QTLs underlying PH, SW, Oil and Protein were identified by QTL ICIMapping V3.3 software [64]. Inclusive Composite Interval Mapping (ICIM) and Composition interval mapping (CIM) methods were used to identify the QTLs. The threshold of logarithm of odds (LOD) score for evaluating the statistical significance of QTL effect was determined using 1000 permutations at the 5% level of significance. The location of a QTL was described according to its LOD peak location and thesurrounding region with 95% confidence interval [65]. As a result,