Introduction

Pod dehiscence is necessary for spreading seeds of wild plant species that bear seeds in pods and allows wild plants to introduce their progeny into a broad range of environments (Fuller 2007). However, in cultivated crops, seeds from plants that undergo pod dehiscence cannot be harvested, resulting in serious yield losses. During planting and harvesting, resistance to pods dehiscence prior to harvest is critical (Hancock 2004; Harlan et al. 1973). Therefore, during soybean domestication, pod dehiscence resistance was artificially selected as an important trait (Hideyuki et al. 2012). Although researchers intensively screen for and avoid pod dehiscence during crop domestication, pod dehiscence before harvest is still a problem in soybean breeding (Christiansen et al. 2002). Soybean breeding has made great progress in incorporating resistance to pod dehiscence, but little is known about the genetics of pod dehiscence in soybean (Funatsuki et al. 2014). The main limiting factor for genes conferring resistance to pod dehiscence is low marker coverage. This leads to large intervals for detecting and mapping QTLs related to pod dehiscence; a large number of predicted candidate genes are contained within these intervals (Zhang et al. 2016).

The rapid development of soybean genomics research has greatly enhanced the speed and the accuracy of QTL mapping for many important agronomic traits. The completion of the soybean reference genome, Williams 82, has accelerated the construction of a high-density genetic linkage map, allowing for more accurate QTL mapping and gene mining in soybean (Schmutz et al. 2010; Ju et al. 2017; Zhang et al. 2016). In soybean, researchers have constructed many genetic linkage maps using markers: restriction fragment length polymorphisms (RFLPs), simple sequence repeats (SSRs) and single-nucleotide polymorphisms (SNPs) (Gutierrez-Gonzalez et al. 2011; Hyten et al. 2010; Jeong 2002; Keim et al. 1990). To date, multiple studies have identified major QTLs for pod dehiscence in soybean using recombinant inbred lines (RILs) with RFLP (Bailey et al. 1997) and SSR markers (Funatsuki et al. 2006). These studies have shown that pod dehiscence is dominated by several major and minor QTLs. Due to the various genetic backgrounds in different soybean populations, pod dehiscence was found to be controlled by a single recessive major gene or multiple genes (Tsuchiya 1986, 1987; Kang et al. 2005; Funatsuki et al. 2008; Yamada et al. 2009). A major pod dehiscence QTL was discovered on chromosome 16 using RILs, and other minor QTLs were identified on chromosomes 2, 15 and 19 (Bailey et al. 1997). After fine mapping, the candidate gene Glyma16g25580 was confirmed to control pod dehiscence and was designated as Pod dehiscence 1 (Pdh1) (Funatsuki et al. 2014). Meanwhile, a NAC gene named SHAT1-5, Glyma.16g019400, was found to regulate pod dehiscence in soybean (Dong et al. 2014). Using SSR markers, three novel minor QTLs for pod dehiscence were confirmed on chromosomes 5, 10 and 14 from two sets of RILs (Kang et al. 2009).

Next-generation sequencing (NGS) has provided many sequencing technologies that have been used to construct high-density genetic maps on many populations. NGS techniques include restriction-site-associated DNA tagging sequencing (RAD-seq), genotyping-by-sequencing (GBS) and specific-locus amplified fragment sequencing (SLAF-seq) (Baird et al. 2008; Elshire et al. 2011; Sun et al. 2013). These techniques have become effective tools for QTL mapping and allow researchers to identify, sequence and genotype individuals in a population. Since NGS can directly identify DNA sequence differences with high accuracy, NGS techniques have been widely used for plant and animal genomics (Bhatia et al. 2013). SLAF-seq is a rapid, high-resolution technique based on NGS technology for broad-scale SNP genotyping, which has been extensively used in genomic analyses due to its characteristics of high molecular marker density, large marker number, even distribution, high accuracy and lack of repetitive markers (Sun et al. 2013). This technique can improve the efficiency of QTL mapping associated with important agronomic traits (Cao et al. 2017).

In this study, we selected SLAF markers based on the whole soybean genome through SLAF-seq. Using a set of 260 RILs derived from the cultivars Heihe 43 (ZDD24325) and Heihe 18 (ZDD23620) as parents, a high-density genetic linkage map covering the entire soybean genome was constructed to map QTLs for pod dehiscence. We obtained several QTLs associated with pod dehiscence. The high-density genetic map we generated will facilitate exploration of effective QTLs, and the novel QTLs detected in our population will be helpful for further research on pod dehiscence in soybean.

Materials and methods

Plant materials and phenotypic evaluation

We used 260 F2:7 recombinant inbred lines (RILs) that were developed by a single-seed descendent (SSD) method from an F2 population of the cross between the cultivars Heihe 43 (ZDD24325) and Heihe 18 (ZDD23620) as parents. The two parents and their derived RIL population were grown at the same location of Heihe Experiment Station, Heilongjiang Academy of Agricultural Sciences, from 2015 to 2017. The parents and RIL population were planted four rows per plot (5.0-cm plant spacing, 66.7-cm row spacing and 4.0-m row length). Heihe 43 is resistant to pod dehiscence and is the largest precocious cultivated variety in Heilongjiang province, developed from Heihe 18 and Heihe 23 (ZDD23625) by systematic breeding programs. Heihe 18 is sensitive to pod dehiscence. Although Heihe 18 and Heihe 43 have minimal differences among numerous agronomic traits, the pod dehiscence in Heihe 18 is extremely obvious (Han et al. 2015a).

We used a modified oven-dry method (Kang et al. 2005) to determine the phenotype of pod dehiscence. A total of 30 pods per line per year were carefully collected using scissors at the R8 stage and stored in sealed bags to prevent water evaporation. The pods were then incubated at 80 °C for 5 h in an oven, and the number of dehiscent pods was calculated. The PD was estimated using the following equation: PD = (dehiscent pods number/total pods number) × 100% (Peng et al. 1991).

Statistical analysis

One-way analysis of variance (ANOVA) was performed by the statistical package SAS version 9.1, including the frequency distribution, the mean of the RIL population, the coefficient of variation (CV), the broad-sense heritability (h2). The h2 for PD was estimated using the following equation: h2 = \({\sigma }_{\mathrm{G}}^{2}\)/\({\sigma }_{\mathrm{P}}^{2}\) (Nyquist and Baker 1991), where \({\sigma }_{\mathrm{G}}^{2}\) and \({\sigma }_{\mathrm{P}}^{2}\) are the genotypic and the phenotypic variance, respectively.

DNA extraction

Young, healthy and fresh leaves from both parent lines and all 260 RILs were collected in centrifuge tubes, frozen in liquid nitrogen, ground in a tissue grinder and then stored at − 80 °C. Total genomic DNA was extracted from every leaf sample following the modified CTAB protocol (Saghaimaroof et al. 1984). The quality and concentration of the extracted DNA were assessed by electrophoresis on 1% agarose gels and using a spectrophotometer (UV–Vis Spectrophotometer Q5000).

SLAF library construction and sequencing

We constructed a SLAF library and sequenced individuals from 260 RILs and their parents. The Williams 82 soybean reference genome sequence (Wm82.a2.v1, https://phytozome.jgi.doe.gov) was used to predict digestion, and we chose a combination of RsaI and HaeIII restriction enzymes to digest the genomic DNA. Amplified fragments ranging from 364 to 414 bp in size were defined as SLAF markers; 132,516 SLAF markers were predicted. A single-nucleotide (A) overhang was added to the obtained digested fragments. Dual-index sequencing adapters were then ligated to the A-tailed fragments. These fragments were generated using PCR and were then purified and mixed. The target fragments were isolated to generate a sequencing library. After purification and dilution, paired-end (each end 125 bp) sequencing was performed on an Illumina HiSeq platform (Illumina, Inc, San Diego, CA, USA). To evaluate the accuracy of the SLAF library, Oryza sativa L. Japonica was selected as a control for the same treatment and was used in library construction and sequencing. In comparison with the control data, the efficiency of enzyme digestion was evaluated to assess accuracy and effectiveness. Using reads clustering, SLAF markers were developed in parents and their offspring. Polymorphic SLAF markers were screened and used for further genetic map construction.

SLAF-seq data grouping and genotyping

The procedure of Sun et al. (2013) was adapted for SLAF marker grouping and genotyping of the 260 RILs as follows. The original sequencing read length of the SLAF-seq library was 125 bp. In order to ensure analysis quality, the original sequencing reads were filtered by removing reads containing: (1) barcode adapters; (2) > 10% base content (N); and (3) residue of restriction enzyme fragments. We selected reads ranging from 4 to 103 bp to analyze the data. Filtered reads were aligned to the reference genome using the Burrows–Wheeler Aligner (BWA) (Os et al. 2005), and reads with the same paired end were identified as the same SLAF marker. SLAF markers were analyzed for polymorphism and classification based on the number of alleles and the differences between sequences. SLAF markers were then mapped to the reference genome. The distributions of SLAF markers on chromosomes were used to map SLAF markers and polymorphic SLAF markers on the chromosomes.

Then, all polymorphic SLAF markers for genetic mapping were filtered using the following criteria: (1) average sequence depths of more than tenfold; (2) filtering out SLAF markers with more than five SNPs; and (3) integrity filtering. Markers were required to cover at least 70% of the genotypes in all progeny. Markers with missing data were filtered based on parental genotypes. The final polymorphic SLAF markers were used to construct a high-density genetic map.

Construction of a high-density genetic map

Polymorphic SLAF markers were aligned with the reference genome and mapped onto 20 chromosomes. We then calculated the MLOD (the modified logarithm of odds) scores between markers and filtered the markers with MLOD values lower than 5. Each chromosome was a linkage group (LG) used as a unit. HighMap software (Li et al. 2008) was used to analyze the linear arrangement of markers in the linkage group and to calculate the map distances. Finally, we obtained a high-density genetic map.

QTL mapping for pod dehiscence

According to the phenotypes in three environments, QTLs for pod dehiscence were detected using inclusive composite interval mapping (ICIM) in the R/qtl package (Broman et al. 2003). A total of 1000 permutation tests at a 95% confidence level were used to set the LOD threshold. Based on 1000 permutations, LOD = 2.5 was used to determine the presence of a putative QTL associated with a target trait in a certain genomic region.

Gene expression analysis by real-time qPCR (qRT-PCR)

Total RNAs of pods (R6 stage) were extracted from Heihe 43 and Heihe 18, respectively, using TransZol Up Plus RNA Kit (TransGen Biotech). For reverse transcription, the first-strand cDNA synthesis was performed using the TransScript First-Strand cDNA Synthesis SuperMix Kit (TransGen, China). For qRT-PCR, gene expressions were examined using cDNA templates on an Applied Biosystems 7300 Real-Time PCR System. Gene-specific primers for the candidate genes were designed using Primer3.0. The relative gene expression levels followed the \({2}^{-\Delta\Delta\mathrm{Ct}}\) method (Pfaffl 2001). The mRNA level of GmActin (Glyma18g52780) gene as a reference for normalization and three biological replicates were used for each gene.

Candidate gene prediction within QTL intervals

Sequences within QTLs were defined according to the Williams 82 soybean reference genome sequence (Wm82.a2.v1, https://phytozome.jgi.doe.gov). The functions of candidate genes were annotated using Blastx program (https://www.geneontology.org/) in Nr (nonredundant), Swiss-Prot and KOG/COG (clusters of orthologous groups). All genes were categorized by Gene Ontology (GO) annotation (https://www.geneontology.org/).

Development and application of molecular markers within QTLs

According to the different loci within QTLs, we developed CAPS (cleaved amplified polymorphic sequence) and dCAPS (derived cleaved amplified polymorphic sequence) markers and used these molecular markers to detect genotypes of the RILs. The online software dCAPS Finder 2.0 (https://helix.wustl.edu/dcaps/dcaps.html) was used to confirm all enzymes. Within qPD05-1, we chose a restriction site and developed a CAPS marker, in which enzyme Mnl I was used as a marker. The qPD01 QTL has no restriction enzyme sites, so we used the reverse complementary sequence and artificially introduced a mismatch base. With this method, we found a restriction enzyme site to develop a dCAPS marker with HpaII as a candidate marker enzyme. Indel markers were developed based on the Indel sites within qPD08-1. PCR primers were designed on both sides of the candidate CAPS/dCAPS markers and the Indel markers using Primer3 software (https://primer3.ut.ee/).

Genomic DNA from parents and 260 RILs were used as templates. The 20-μL PCR contained 100 ng of template DNA, 10 × PCR buffer, 2 mmol/L dNTPs, 2 mmol/L primers and 1 unit of Taq polymerase (Takara). The PCR amplification program was as follows: 95 °C for 4 min, 34 cycles of 95 °C denaturation for 30 s, 58 °C annealing for 40 s, 72 °C extension for 1 min and a final 72 °C extension for 10 min. PCR products were detected using 0.8% agarose gel electrophoresis. In accordance with the restriction enzyme digestion method (New England BioLabs, NEB), the 10-μL system consisted of 5 μL of PCR product, 0.2 μL 10 U/μL of enzyme, 1.5 μL of buffer (NEB, www.neb.com/) and 3.3 μL of ddH2O and was digested at 37 °C in a water bath for 40 min. The digested products of CAPS/dCAPS were detected using 7% polyacrylamide gel electrophoresis, and PCR products of Indels were denatured and then separated using 7% denaturing polyacrylamide gel electrophoresis (PAGE).

Results

Phenotypic analysis of pod dehiscence

We identified the PD of parents and the RILs in Heihe, Heilongjiang province, in 2015, 2016 and 2017 (Fig. 1 and Table 1). Heihe 18 had higher PD than Heihe 43 at 37% and 3% in 2015, 47.03% and 11.36% in 2016 and 75% and 3.33% in 2017. RIL population showed more resistance of pod dehiscence over 3 years, showing that the female Heihe 43 contributed to resistance of pod dehiscence and the positive skewness indicated a transgressive segregation toward lower PD. The broad-sense heritability of pod dehiscence was 81.3% in our population.

Fig. 1
figure 1

Frequency distributions of pod dehiscence in RILs from 2015 to 2017. The cultivars Heihe 43 and Heihe 18 indicated by the arrows designate the degree of pod dehiscence in parental lines

Full size image
Table 1 Statistical analysis of pod dehiscence in the parents and the RIL population in 3 years in Heihe
Full size table

SLAF-seq and genotyping of RILs

The parent lines and RILs were sequenced and genotyped using SLAF-seq. The sequence depths in Heihe 43 and Heihe 18 were 39.74 × and 34.87 × , respectively, and the number of SLAF markers was 195,279 and 205,644, respectively. The average sequence depth was 12.35 × . The enzyme digestion efficiency was 89.82%. There were 398,386 controlled sequencing reads that were used to evaluate the accuracy of library construction; these reads were compared with the reference soybean genome. The alignment efficiency was 86.29%. The average Q30 of sequencing was 80.40%, and the average GC content was 38.12%. In total, we obtained 364,461 SLAF markers, of which 24,249 were polymorphic and the proportion of polymorphisms was 6.65%. After filtering and quality assessment, 4593 SLAF markers were used to construct a high-density genetic map. Chromosomal distribution of all SLAF markers and polymorphic SLAF markers was plotted according to the distribution of SLAFs on the chromosomes (Fig. 2). Compared with the soybean reference genome, we identified a total of 6834 SNPs, of which 4576 were transition type (Tri) and 2258 were the transversion type (Trv) (Table 2).

Fig. 2
figure 2

Distribution of SLAF markers and polymorphic SLAF markers on chromosomes. The abscissa is the length of the chromosome, and each yellow stripe represents a chromosome. The genomes were divided in units of 1 M. Different colors indicate the number of SLAFs, the color is deeper and the SLAFs are more. Darker regions in the figure are concentrated distributions of SLAF markers. The upper panel of the figure represents the distribution of the SLAF tag and the bottom panel of the figure represents the distribution of polymorphism of SLAF markers (with M as the unit)

Full size image
Table 2 The SNP information of the mapped SLAF markers
Full size table

Construction of a high-density genetic map in soybean

A total of 4593 polymorphic markers were mapped on 20 chromosomes based on the reference genome, and a high-density genetic map was constructed, with a total length of 1478.86 cM and an average distance between markers of 0.53 cM (Fig. 3). The mean chromosome length was 73.94 cM. Seven hundred and six markers were assigned on chromosome 3 with a genetic length of 101.35 cM, which was the longest of all chromosomes. Seventeen markers were assigned on chromosome 12 with a genetic length of 17.27 cM, which was the shortest of all chromosomes. The minimum number of markers was consistent with its physical length (Schmutz et al. 2010). The proportion of gaps < 5 cM between two markers was 95.72% (Table 3). Collinearity analysis (Fig. 4) of the position of markers and the genetic map on the genome showed that the order of markers on the 20 chromosomes was consistent with the genome, indicating high collinearity, indicating that the gene annotation within QTL intervals was reliable.

Fig. 3
figure 3

High-density linkage map. SLAF markers are distributed on 20 chromosomes. The black bars in each linkage group represent mapped SLAF-seq markers. The linkage group number is shown on the x axis, and genetic distance is shown on the y axis (cM is the unit)

Full size image
Table 3 Characteristics of the high-density genetic map
Full size table
Fig. 4
figure 4

Collinearity analyses of the genetic map and genome. The abscissa is the genetic distance of each linkage group; the ordinate is the physical length of each linkage group, which scatters the form of markers in the genome and genetic map collinearity. Different colors represent different chromosomes or linkage groups

Full size image

QTL for pod dehiscence

The R/qtl package was used to identify QTL associated with pod dehiscence in 3 years. In total, six novel QTLs were detected and were located on chromosomes 1, 5, 8 and 14 (Fig. 5). All of QTLs were not found in previous studies; thus, these QTLs could be considered as novel QTLs for pod dehiscence. The LOD score of these QTLs ranged from 2.61 to 11.00, and they could explain 7.22% to 24.44% of the phenotypic variation (Table 4). Of these, four QTLs (qPD01, qPD05-1, qPD05-2 and qPD08-1) explained high phenotypic variations ( > 10%). In addition, qPD01, qPD05-1 and qPD08-1, distributed on chromosomes 1, 5 and 8, respectively, could be detected in different years. Thus, these three QTLs might be the major and stable QTLs in our population. The additive effect of these QTLs was negative, which showed that the positive alleles came from the female Heihe 43 and the positive alleles improved the resistance of pod dehiscence.

Fig. 5
figure 5

Mapping of QTLs for pod dehiscence on Chr01, Chr05, Chr08 and Chr14. The curves indicate the physical position of markers against LOD score of QTL detected on chromosomes. Different lines represent different years

Full size image
Table 4 QTL for pod dehiscence in soybean across three environments
Full size table

Candidate genes prediction and expression analysis

Based on the Williams 82 soybean reference genome, six QTL intervals contained 639 genes located on chromosomes 1, 5, 8 and 14. Currently, only 217 genes functions have been annotated by GO annotation (Fig. 6). To further select candidate genes, we compared all these sequence variations within the QTL intervals and screened out 191 SNPs and 23 Indels in 34 genes that were different between the two parents (Table 5). Of 34 genes, 32 genes were annotated their functions (Table 6). Except three genes without suitable primers to be tested, using gene-specific primers for the candidate genes (Table 7), the expression levels of the 31 candidate genes (Fig. 7) were investigated in the pods (R6 stage) collected from Heihe 43 and Heihe 18 by qRT-PCR and different levels of transcription abundance were observed among these genes. Among 31 genes, nine genes had relatively high expression levels, including Glyma.01g045700, Glyma.01g045800, Glyma.01g046000, Glyma.05g005600, Glyma.05g225900, Glyma.05g227400, Glyma.08g271900, Glyma.08g274500 and Glyma.14g195200. Seven of nine genes (Glyma.01g045800, Glyma.01g046000, Glyma.05g005600, Glyma.05g225900, Glyma.05g227400, Glyma.08g271900 and Glyma.08g274500) had significant differences between two parents.

Fig. 6
figure 6figure 6

Gene ontology (GO) annotation of genes within the QTLs. The y axis on the right shows the number of genes in each category, and the y axis on the left shows the percentage of a specific category of genes in that main category

Full size image
Table 5 SNP and Indel information of the QTL interval
Full size table
Table 6 Information of candidate genes
Full size table
Table 7 Primer sequences for qRT-PCR
Full size table
Fig. 7
figure 7

Expression level of 31 candidate genes. Y axes indicate the ratios of relative fold expression levels. Asterisks indicate significant differences as determined by ANOVA (**P < 0.01). Relative expression was calculated based on the expression level of the GmActin gene

Full size image

Development of molecular markers

According to different sites within the QTL interval, PCR amplifications were conducted using the designed primers (Table 8), and we obtained single PCR products with the same size as the target fragments. Based on alignment with Williams 82 reference sequence, there are differences between the two parents, which was consistent with the prediction.

Table 8 Information of molecular markers developed within QTLs
Full size table

After digestion with Mnl I, three digested products were produced theoretically by dCAPs marker in the qPD01 interval. However, one fragment of 25 bp was too small to detect; therefore, only two kinds of digestion products were detected. The DNA fragment size without digestion was 157 bp, and the allelic variation was T. Material carrying the T allele showed the pod dehiscence. After digestion, there were two fragments, 132 bp and 25 bp. The allelic variation was C, and material carrying the C allele showed the resistance of pod dehiscence. Within the qPD05-1 interval, the PCR products could be digested by HpaII and produced three digestion products theoretically. Only two kinds of digestion products could be detected due to the similar length of the two fragments after digestion. The DNA fragment size without digestion was 540 bp and the allelic variation was A that provided the resistance of pod dehiscence; the digested DNA fragment sizes were 294 bp and 246 bp and the allelic variation was G that provided the pod dehiscence. Within qPD08-1, we detected two bands: One DNA fragment size was 422 bp and the allelic variation was GTTT; the other DNA fragment size was 419 bp and the allelic variation was G that existed in the resistance of pod dehiscence. These markers were validated by using 260 RILs derived from two parents. The genotypes and phenotypes of five RILs are listed in Fig. 8 and Table 9. We found eight allele combinations in 260 RILs, and three combinations had more than ten individuals (Table 10). The identification efficiency to resistance of pod dehiscence for combination 1 with all three resistant alleles was the highest (97.4%), and both combination 2 and combination 4 with two resistant alleles were 84.2% and 82.4% respectively. The results indicated that the developed CAPS/dCAPS and Indel markers could be used to identify genotypes and predict the phenotypes of soybean varieties.

Fig. 8
figure 8

Validation of the molecular markers to QTLs (qPD01, qPD05-1 and qPD08-1). Electrophoresis results of the polymerase chain reaction (PCR) of molecular markers for Heihe 43, Heihe 18 and five recombinant inbred lines (RILs)

Full size image
Table 9 The PD of the RILs used to identify the molecular markers
Full size table
Table 10 Identification efficiency of three molecular markers different allelic combination on resistant materials on 260 RILs
Full size table

Discussion

QTL mapping is an effective method for analyzing quantitative traits in soybean. The quality of a genetic map has a significant impact on the resolution of QTL mapping. For a specific population, increasing marker density improves the accuracy and efficiency of a genetic map (Gutierrez-Gonzalez et al. 2011; Yu et al. 2011; Zou et al. 2012). Nevertheless, both wild and cultivated soybeans exhibit high linkage disequilibrium (LD) (cultivated soybean: ~ 150 kb; wild soybean: ~ 75 kb), and the average genetic distance is greater than in other crops (Lam et al. 2010). Therefore, a high-density genetic map is needed to improve the efficiency and accuracy of QTL mapping in soybean and can be useful for molecular marker-assisted selection (Cao et al. 2017; Staff 2014).

SLAF-seq is an efficient method for identifying QTLs and genotyping based on NGS. This technique is advantageous for large numbers of individuals due to its convenient library preparation procedure (Elshire et al. 2011). Compared with SLAF-seq, genome resequencing provides a more accurate view of the genome and captures large variations as well as small ones that may be missed by SLAF-seq. However, SLAF-seq is relatively low cost and can be used for mapping in large-scale populations (Li et al. 2014; Sun et al. 2013). SLAF-seq provides the best balance between sequencing cost and sequencing depth. SLAF-seq also obtains genome-wide variation sites at a lower cost and provides for researchers to select functional polymorphic markers, construct high-density genetic maps and further identify functional genes of important agronomic traits (Wei et al. 2017). This method has been extensively used to construct high-density genetic maps and identify major QTLs related to complex traits in several plants (Xia et al. 2015; Xu et al. 2015; Han et al. 2015a, b; Qin et al. 2015; Su et al. 2016). Pod dehiscence is one of the most important agronomic traits in crops and is directly related to yield. Therefore, in many crops, including soybean, it is vital to identify the QTLs related to pod dehiscence and the markers tightly linked to pod dehiscence. In soybean, several studies have found QTLs for pod dehiscence located on chromosome 16 and have performed fine mapping using different plant materials (Bailey et al. 1997; Dong et al. 2014; Funatsuki et al. 2014; Gao and Zhu 2013; Hideyuki et al. 2012; Kang et al. 2009), but only two genes for these pod dehiscence QTLs have been cloned. Therefore, we detected novel QTLs associated with pod dehiscence in soybean using SLAF-seq and predicted candidate genes within QTLs based on the high-density genetic maps, which could be applied for identifying new QTLs for pod dehiscence and other valuable agronomic traits. The QTLs for pod dehiscence also provided promising candidate genes for further characterization.

Constructing a high-density genetic map

Because of multiple generations of self-fertilization during the process of development of RIL populations, the observed recombination of RIL populations is greater than that of the corresponding F2 population. The increased recombination for a reference population that consists of RILs, when compared to an F2 reference population improves the accuracy and efficiency of the map. Moreover, in map-based cloning, secondary mapping populations are usually mapped after primary populations; however, doing so is laborious and time-consuming. In this study, we constructed a high-density genetic map and sequenced the whole genome in soybean RIL7 using SLAF-seq. High-quality markers were evenly distributed on 20 chromosomes. Recently, Cao et al. (2017) and Li et al. (2014) used SLAF-seq to construct high-density soybean genetic maps using 236 RILs and 200 RILs, respectively. These two genetic maps contained 3255 and 3541 SLAF markers, respectively, and the average genetic distances between two markers were 0.66 cM and 0.72 cM, respectively. Researchers have shown that increasing population number and marker density can significantly improve the efficiency of QTL mapping (Li et al. 2014). In our study, we used 260 RILs to construct a high-density genetic map containing 4593 SLAF markers. Compared with the genetic maps of Cao et al. (2017) and Li et al. (2014), the genetic distance between two adjacent SLAF markers was reduced to 0.53 cM in our map. Each chromosome contained 230 markers on average, and the shortest genetic distance between adjacent markers was 0.11 cM. Our high-density soybean genetic map may be due to the large size of the population, the increased number of SLAF markers and narrowed gaps among markers, which showed that our genetic map is an effective reference for soybean and can be used for further genetic mapping and prediction of pod dehiscence-related genes.

QTLs for pod dehiscence

The trait of pod dehiscence has considerable genetic diversity in soybean cultivars (Caviness 1965; Tsuchiya 1986; Helms 1994; Romkaew and Umezaki 2006; Yamada et al. 2009). Tsuchiya (1986) predicted that three genes control pod dehiscence. Several studies have identified two major genes controlling pod dehiscence, named qPDH1 and SHAT1-5, which are located on chromosome 16 (Bailey et al. 1997; Dong et al. 2014; Funatsuki et al. 2014; Gao and Zhu 2013; Hideyuki et al. 2012; Kang et al. 2009). Other minor QTLs related to pod dehiscence are located on chromosomes 2, 5, 10, 14, 15 and 19 (Bailey et al. 1997; Kang et al. 2009). In China, Luo et al. (2012) used 112 RILs and the composite interval mapping method and reported a QTL controlling pod dehiscence located in chromosome 6 named qPDH6-1.

Kang et al. (2009) identified seven minor QTLs for pod dehiscence using two RIL populations. Compared to this genetic map, our genetic map was constructed using a specific RIL population, which allowed us to construct a more efficient QTL map for a given phenotypic trait. We identified six novel QTLs for pod dehiscence on chromosomes 1, 5, 8 and 14 that were not detected in previous studies (Table 2). Of these six QTLs, qPD05-1 and qPD05-2 were located on chromosome 5, the same chromosome as a previously reported minor QTL (Satt 385) (Kang et al. 2009). In the USDA (United States Department of Agriculture) map, qPD05-1 is located between Sat_374 and BARCSOYSSR_05_1304, and qPD05-2 is located between BARCSOYSSR_05_0375 and BARCSOYSSR_05_0374. Interestingly, qPD05-1 and qPD05-2 were not at the same position as QTL Satt 385. Meanwhile, qPD05-1 could be found across 2 years (2015, 2016), of which the LOD and the phenotypic variation reached maximum values of 11.19% and 15.14%, respectively. Therefore, this QTL could be identified as a major, stable QTL for pod dehiscence. The other QTL, qPD14, was located on chromosome 14, the same chromosome as a previously reported minor QTL (Satt 577) (Kang et al. 2009). Similar in the previous case, qPD14 was located between Satt 063 and Sat_424. These two QTLs were not at the same position at chromosome 14.

As given in Table 2, the QTLs qPD01 and qPD08-1 could also be detected in different environments with phenotypic variation of 24.44% and 19.61%, respectively. Therefore, in our population, qPD01, qPD05-1 and qPD08-1 might be the major QTLs for pod dehiscence. qPD05-2, qPD08-2 and qPD14 were detected in a single environment; the phenotypic variation explained by qPD05-2 was 17.57%, qPD08-2 was 8.37%, and qPD14 was 7.22%. QTLs detected in the progeny could be considered stable in different environments. Therefore, these QTLs and closely linked SLAF markers can be used for molecular marker-assisted selection in soybean.

It is interesting to note that the previously reported major QTLs on chromosome 16 were not identified in our study in all environments; this result might be explained by the different genetic background between the parents of the mapping populations; the same phenomenon was found in the other trait (Li et al. 2014). The QTLs that have been identified for pod dehiscence QTL such as Pdh1 (Funatsuki et al. 2014) may be present in elite populations developed from crossing parents that are used for commercial soybean varieties including the parent we used Heihe 43, but we did not map this locus in our population. However, the QTLs that we have identified by crossing two parents that also used for developing commercial soybean varieties may have different mechanism of pod dehiscence resistance. Heihe 43 is the largest planted variety in Heilongjiang province of the People’s Republic of China at present; therefore, the breeding value of the QTLs we identified for pod dehiscence resistance can be contributed to genetic research of pod dehiscence in soybean.

Gene annotation and putative candidate genes for pod dehiscence

We predicted candidate genes within a narrow region between two adjacent markers based on high-density genetic maps and the high-quality genome sequences of Williams 82. Gene prediction depends on high accuracy and collinearity between the genetic map and the reference genome (Dan et al. 2016). The high collinearity of the markers and genetic map in our map allowed us to identify gene candidates related to pod dehiscence within QTLs (Fig. 4). Gene annotation within six QTLs was performed using Nr, Swiss-Prot and KOG/COG databases. These QTLs contained 34 protein-encoding genes, of which 32 have GO annotations (Table 6). In order to obtain more information of candidate genes in pod dehiscence, the expression profiles of 31 genes (Fig. 7) between Heihe 43 and Heihe 18 were analyzed by qRT-PCR. Among 31 genes, nine genes had relatively high expression levels, and seven of nine genes had significant differences between two parents, including Glyma.01g045800, Glyma.01g046000, Glyma.05g005600, Glyma.05g225900, Glyma.05g227400, Glyma.08g271900 and Glyma.08g274500. Previous studies have suggested that some enzymes control cell wall degradation (Christiansen et al. 2002) and synthesis (Romkaew et al. 2008). In our study, four genes involved some enzymes including fatty acid/sphingolipid desaturase (Glyma.01G046000), DNA methyltransferase (Glyma.05G005600), NADP-dependent malic enzyme (Glyma.05G227400) and oxidoreductase (Glyma.08G274500), but the function of these genes needs to be verified. Our results provided information for future study.

In this study, we constructed a high-density genetic map containing 4593 SLAF markers using 260 RILs. Our finding showed that this map was accurate and efficient for QTL mapping. We identified three stable QTLs associated with pod dehiscence resistance in three environments. These novel and stable QTLs not only contribute to study of the genetic mechanism of pod dehiscence but also facilitate molecular marker-assisted selection to increase production in soybean.

Author contribution statement

JH conducted the data analysis, QTL mapping, genomic comparative analysis and wrote the manuscript. DH and HY conducted the RIL populations and phenotyping. YG, ZW and YT provided advice on experimental implementation. LQ conceived and supervised the project, reviewed and revised the manuscript.