Honglei Ren, Jianan Han, Xingrong Wang, Bo Zhang, Lili Yu, Huawei Gaoa,,

Huilong Honga,b, Rujian Suna,b, Yu Tianb, Xusheng Qic, Zhangxiong Liub,*, Xiaoxia Wua,*,Li-Juan Qiub,*

aCollege of Agriculture, Northeast Agricultural University, Harbin 150030, Heilongjiang, China

bNational Key Facility for Crop Gene Resources and Genetic Improvement, Institute of Crop Sciences, Chinese Academy of Agricultural Sciences,Beijing 100081, China

cInstitute of Crop Sciences, Gansu Academy of Agricultural Sciences, Lanzhou 730070, Gansu, China

dSchool of Plant and Environmental Sciences, Virginia Polytechnic Institute and State University, Blacksburg, VA 24060, USA

eHeilongjiang Academy of Agricultural Sciences, Harbin 150086, Heilongjiang, China

Keywords:Soybean drought tolerance Simplified genome sequencing Quantitative trait loci Plant height Seed weight per plant

A B S T R A C T Drought stress is an important factor affecting soybean yield.Improving drought tolerance of soybean varieties can increase yield and yield stability when the stress occurs.Identifying QTL related to drought tolerance using molecular marker-assisted selection is able to facilitate the development of drought-tolerant soybean varieties. In this study, we used a high-yielding and drought-sensitive cultivar‘Zhonghuang 35’and a drought-tolerant cultivar ‘Jindou 21’ to establish F6:9 recombinant inbred lines. We constructed a highdensity genetic map using specific locus amplified fragment sequencing (SLAF-Seq)technology. The genetic map contained 8078 SLAF markers distributing across 20 soybean chromosomes with a total genetic distance of 3780.98 cM and an average genetic distance of 0.59 cM between adjacent markers. Two treatments (irrigation and drought) were used in the field tests, the Additive-Inclusive Composite Interval Mapping (ICIM-ADD) was used to call QTL,and plant height and seed weight per plant were used as the indicators of drought tolerance.We identified a total of 23 QTL related to drought tolerance.Among them,seven QTL (qPH2, qPH6, qPH7, qPH17, qPH19-1, qPH19-2, and qPH19-3) on chromosomes 2, 6, 7, 17,and 19 were related to plant height,and five QTL(qSWPP2,qSWPP6,qSWPP13,qSWPP17,and qSWPP19) on chromosomes 2, 6, 13, 17, and 19 were related to seed weight and could be considered as the major QTL. In addition, three common QTL (qPH6/qSWPP6, qPH17/qSWPP17, and qPH19-3/qSWPP19) for both plant height and seed weight per plant were located in the same genomic regions on the same chromosomes. Three (qPH2, qPH17, and qPH19-2)and four novel QTL(qSWPP2,qSWPP13,qSWPP17,and qSWPP19)were identified for plant height and seed weight per plant, respectively. Two pairs of QTL (qPH2/qSWPP2 and qPH17/qSWPP17)were also common for both plant height and seed weight per plant.These QTL and closely linked SLAF markers could be used to accelerate breeding for drought tolerant cultivars via MAS.

1. Introduction

Soybean [Glycine max (L.) Merr.] is one of the most important crops providing protein and oil for human and animal consumption [1]. However, drought is a serious abiotic stress for any crops [2], and soybean is highly sensitive to water deficit especially in arid and semi-arid areas [3]. More than 50% of water consumption is used for global agriculture, of which 64.61%of water consumption is used for agriculture in China and more than 90% in India and Africa [4]. In arid regions, natural precipitation and supplementary irrigation cannot meet the demand of soybean growth and development,which can reduce soybean yield up to 40%[5].The worst scenario is that there is no production at all when water deficit occurs during the reproductive stage, preventing seed formation of crops [6]. Therefore, it is very important to improve drought tolerance in soybean to increase yield in water deficient areas.

Studies have demonstrated that drought tolerance is a complex trait controlled by multiple genes[7].Identifying QTL associated with drought tolerance and understanding its mechanism can accelerate soybean genetic researches and varietal improvement. So far, 14 QTL related to drought tolerance have been identified in independent studies. Using 184 recombinant inbred lines (RILs) derived from a cross between ‘Kefeng 1’ and ‘Nannong 1138-2’, ten QTL were identified that were associated with drought tolerance using yield and drought susceptibility indexes in different environments including water-stressed and well-watered conditions in the field and greenhouse tests. These QTL are located on chromosomes 1, 5, 6, 7, 12, 16, 17, and 20, respectively [8].Taking the phenotypic response to a silver inhibitor as an indicator of drought tolerance, four QTL related to transpiration rate were localized on chromosomes 3,5,10,and 12 using a RIL population derived from the ‘PI 416937’ × ‘Benning’cross, which explained 17.7%–24.7% of the phenotypic variations [9]. Using a novel GWAS (RTM-GWAS, restricted twostage multi-locus genome-wide association analysis) technique and relative root length and relative shoot length as the indicators, 111 drought tolerance-related QTL and 262 alleles were localized on all the soybean chromosomes except for chromosomes 5,16,and 17,explaining 88.55%–95.92%phenotypic variations in a nested association mapping population[10].Twenty-one genes were involved in the drought response pathway in the SoyBase database(https://soybase.org/),which belong to a drought-response protein family and a droughtinduced protein family according to the gene annotation in Arabidopsis thaliana L.

Since the soybean reference genome of Williams 82 was released in 2010[11],many novel genes have been discovered.Furthermore, numerous single nucleotide polymorphisms(SNPs) and insertion-deletions (InDels) have been identified by comparing the sequences of different soybean cultivars[12–14]. Several SNPs and InDels are associated with important traits to adapt to abiotic and biotic stresses [15]. In soybean, SNPs, simple sequence repeats (SSRs), and InDels have been used to construct genetic maps and a large number of QTL have been found and determined to be related to drought tolerance. However, the genetic research of drought tolerance in soybean is difficult and impeded due to the complexity of the traits.The main limiting factors in drought tolerance breeding are technical limitations and low density of markers in genetic maps, limiting the efficiency and accuracy of QTL mapping. With the rapid development of massively parallel sequencing techniques, whole genome sequencing using next generation sequencing technologies,such as specific locus amplified fragment sequencing (SLAFseq), have provided powerful methods to discover SNPs and genotype different cultivars in large populations.

Using the SLAF-seq technology, a number of QTL have been identified in soybean.Li et al.[16]identified a major QTL related to isoflavone content, qIF20-2, which explained 19.6%of the phenotypic variation across different environments. A RIL population derived from the ‘Charleston’ × ‘Dongnong 594’cross was used to identify twelve QTL associated with oil content, of which four pairs of epistatic QTL explained approximately 70% of the phenotypic variation and environmental interaction [17]. Cao et al. [18] used 104 RILs to construct a genetic map with 2062 SLAF markers and found eight QTL related to oil content. Eight out of twenty QTL associated with phosphorous use efficiency were identified across multiple years and treatments[19].Nine novel QTL for fatty acids were detected and explained 0.4%–37.0% of the phenotypic variations [20]. Recently, twenty-four stable QTL were identified for isoflavone components, explaining 4.2%–21.2%of phenotypic variations[21].

The objective of this study was to identify QTL associated with soybean drought tolerance using a high-density genetic map developed by SLAF-Seq.The QTL identified in this study can be utilized as a resource for future identification and genetic determination of drought tolerance genes in soybean.

2. Materials and methods

2.1.Soybean populations and drought resistance assessments

The soybean RILs were derived from the female parent, a high-yielding but drought-sensitive cultivar ‘Zhonghuang 35’(developed by the Institute of Crop Sciences, Chinese Academy of Agricultural Sciences),and the male parent,a droughttolerant cultivar ‘Jindou 21’ (developed by the Institute of Economic Crops, Shanxi Academy of Agricultural Sciences).The F2population containing 234 individuals was developed in 2008.The F3to F6populations were generated by single seed descent and the F6:7to F6:9RILs were generated by the multiple seed descent method. In 2013, 2014, and 2015, the RIL individuals were planted in the Dunhuang Experimental Station (N40°17′, E94°65′) of the Gansu Academy of Agricultural Sciences,which has an annual precipitation of<40 mm.Following the methods of Qi et al. [22], the flat planting land was divided into two parts by an isolation strip 3.0 m wide and 1.0 m high to separate the water stress (WS; drought)treatments and the sufficient water(SW;control)treatments.The two parents and the RILs had three replicates of each treatment in a randomized block design with row length of 3.0 m, row spacing of 50 cm, plant spacing of 10 cm, and 3 rows per plot. At April 11th each year, seeds were sown by hand(with each hole containing two seeds).If both seedlings emerged, one plant was pulled after producing the second trifoliate leaves. The entire field was watered in a rate of 120 m3once before sowing.After seedling emergence,the control(SW) treatment was watered five times (60 m3each) and the WS treatment was not watered during the entire growing period. The field managements were practiced following the local production recommendations.The soil water contents of three specific growth periods are shown in Table 1. Ten random plants were harvested from each plot and used to measure plant height and seed weight of each plant. Soil water contents showed no significant difference among the three years under the WS and SW conditions(Table 1).

2.2. Phenotypic analysis

Plant height and seed weight per plant, as well as the drought resistance coefficient (DC), were calculated for the parental cultivars and the RILs following the “Description Specification and Data Standard for Soybean Germplasm Resources”[23].Plant height(cm)of each plant was measured from the cotyledon node to the growing point of the main stem.Seed weight per plant(g)was an average value of grain weight from 10 plants.The drought resistance coefficients of plant height and seed weight per plant were calculated using the Eq.(1),respectively:

Descriptive statistics, analysis of variance (ANOVA) for phenotypic data, and correlation analysis were performed by OriginPro 8(OriginLab)and R 3.1.3(https://www.r-project.org/).Broad-sense heritability was calculated following the established method[24]:

where, σ2G, σ2GE, and σ2Eare genetic variance, variance of interaction between genotype and environment, and environmental variance, respectively. r and n are the number of repeats and environments,respectively.

2.3. SLAF library construction and sequencing

Fig.1– Frequency distributions of plant height,seed weight per plant and drought resistance coefficient in the RILs under different environments in three years.Arrows indicate‘Zhonghuang 35’and‘Jidou 21’.(a–c)Plant height of the RILs under water stressed treatment(WS)and sufficient water treatment(SW)in 2013,2014,and 2015.(d–f)Seed weight per plant of the RILs under WS and SW in 2013,2014 and 2015.(g–i)Drought resistance coefficient of the RILs in 2013,2014 and 2015.PH,plant height;SWPP,seed weight per plant;DC-PH,drought resistance coefficient by plant height;DC-SWPP, drought resistance coefficient by seed weight per plant.

Genomic DNA was extracted from leaves using the modified CTAB method [25]. Construction and sequencing of the SLAF libraries were performed on 234 RILs and the two parents using the SLAF-seq technique (Biomarker Technologies,Beijing, China) [26]. Briefly, Rsa I and Hae III restriction enzymes were used to digest genomic DNA according to the soybean reference genome (http://phytozome.jgi.doe.gov/pz/portal.html).Then,an adenine was added to the 3′end of the digested fragments and ligated with the Dual-index that was used to distinguish the raw sequencing data [27]. Using PCR amplification, purification, and mixture PCR products, these fragments were obtained and used to construct a sequencing library. Sequencing was performed by the Illumina HiSeqTM 2500 platform after quality certification, including sequence quantity, base distribution, and sequencing quality distribution.To evaluate the accuracy of the SLAF libraries,rice(Oryza sativa L.ssp.japonica cv.Nipponbare)(http://rice.plantbiology.msu.edu/) was used as a control with the same process of library construction and sequencing. The sequences of each sample were clustered according to the sequence similarities and were defined as the SLAF markers.The polymorphisms of SLAF markers were identified by comparing them on different sample sequences.

2.4. SLAF-seq polymorphism analysis and genotyping

To ensure the quality of the bioinformatic analysis, the original sequence data were filtered, and grouping and genotyping of the SLAF-seq data were performed following a standard protocol of data processing.The original sequencing read length of the SLAF-seq library was 125 bp. The reads were filtered according to the following criteria:(1)containing adapter sequence; (2) N content >10%; and (3) restriction enzyme fragments residues.Therefore,the sequences ranging from 4 bp to 103 bp were selected for further analysis. The BWA software (http://bio-bwa.sourceforge.net/) was used to compare the filtered sequencing reads with the reference genome. The same pair-end reads were considered as the same SLAF markers. Polymorphic analysis was performed based on the number of alleles and the differences of gene sequences, and the SLAF makers were classified into polymorphic, non-polymorphic, and repetitive categories. The SLAF markers were mapped to the reference genome (http://phytozome.jgi.doe.gov/pz/portal.html).

For the subsequent genetic analysis,the polymorphic SLAF markers were genotyped. According to the general genetic coding rules, the model aa×bb was used to genotype SLAF markers in the RIL populations. The parental genotypes were aa (male parent) and bb (female parent), and the offspring genotype was aa, bb, and ab. To ensure the quality of the genetic map,the polymorphic SLAF markers were filtered out if the parental sequencing depth was below 10×, the number of SNPs was >3, integrity filtration to screen markers covered<70% of all offspring genotypes, or the parental information on the filtered sites was missing.The final polymorphic SLAF markers were used to construct the high-density genetic map.

2.5. Construction of the high-density genetic map

The selected polymorphic SLAF markers were assigned to different chromosomes by alignment against the reference genome (http://phytozome.jgi.doe.gov/pz/portal.html). The markers with the modified logarithm of odds (MLOD) values<5 were filtered out. Each chromosome was considered as a linkage group and the linear arrangement of markers was obtained using HighMap composition software (http://highmap.biomarker.com.cn/). The genetic distances between adjacent markers were estimated to obtain the high-density genetic map.

2.6. QTL mapping for drought resistance traits in soybean

ICIM-ADD [28] in the QTL Icimapping V4.1 software (http://www.isbreeding.net/) was used to map the drought tolerance related QTL. The likelihood of odd (LOD) threshold was determined as the 95% confidence interval by 1000 permutations and the LOD = 3.9 was considered as a putative QTL related to a certain trait in a genomic region.

2.7. Gene annotations in the QTL intervals

Based on the soybean reference genome sequence(Wm82.a2.v1, https://phytozome.jgi.doe.gov/ pz/portal.html), the sequence of each QTL interval was determined by the BLASTx program (https://blast.ncbi.nlm.nih.gov). The Nr (https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/), Swiss-Prot (http://www.gpmaw.com/html/ swiss-prot.html), and KOG/COG (http://www.ncbi.nlm.nih.gov/COG/) databases in BLASTx were used for gene annotation.

3. Results

3.1. Phenotype analysis of parents and populations

The plant height and seed weight per plant of the two parents and the RIL population were recorded under the SW and WS conditions across three years. Compared to ‘Zhonghuang 35’,‘Jindou 21’ showed significantly higher plant height and greater seed weight per plant, and had a higher drought resistance coefficient (P < 0.05). Under the WS condition, the RIL population showed decreased plant height and seed weight per plant compared to the SW(Fig.1),which indicated that drought stress reduced the growth and development of soybean plants. The heritability of plant height (71.36%–88.44%) was greater than that of seed weight per plant(49.09%–85.46%) in the drought stress group compared to the SW control group.For the RIL population,both genotypes and environmental effects significantly affected plant height and seed weight per plant (Table 2), indicating that there is a variation in drought tolerance in the RIL population.

3.2. SLAF library construction and genotyping of the RIL population

The RsaI and HaeIII enzymes were used to digest genomic DNA of the parents and the RIL population prior to library construction and sequencing. The sequencing depths for‘Zhonghuang 35’ and ‘Jindou 21’ were 32.46× and 34.64×,respectively, and the average sequencing depth of the offspring was 3.87×. The enzyme digestion efficiency was 93.17%. The sequencing Q30 was 84.43% and the GC content was 38.29%. The pair-end alignment efficiency of rice as the control was 84.01%. A total of 297,054 SLAF markers were obtained(Fig.S1A),of which 54,021(14.5%)were polymorphic(Fig.S1B).A total of 43,074 markers were used for genetic map construction(Fig.S2).

3.3. Construction of the soybean genetic map

After filtration and quality assessment, 8078 out of 43,074 SLAF markers were used to construct a high-density genetic map across 20 chromosomes. The map had a total genetic distance of 3480.98 cM and an average genetic distance of 0.59 cM between two adjacent markers(Fig.2).The number of markers ranged from 117 to 722, the genetic length ranged from 155 to 189 cM, and the average distance between twomarkers ranged from 0.24 to 1.37 cM per chromosome (Table 3). The quality of the soybean high-density genetic map was evaluated by collinearity analysis of the position and the genetic map of the SLAF markers on the reference genome(Fig.3),which showed that the SLAF markers were consistent with the soybean genome. The average Spearman coefficientwas 94.8%, indicating that the map had a high collinearity with the consensus map(Table 4, Fig.3).

Table 2–Analysis of variance of agronomic traits of RILs(F-test value).

Fig.2– High density genetic map.Distribution of markers on 20 chromosomes. Black bars on each chromosome represent mapped SLAF-seq markers.The chromosome number is shown on the x-axis,and genetic distance(cM)is shown on the yaxis.

Table 3–Characteristics of the high density genetic map.

3.4. QTL mapping of the RIL population

A total of 23 QTL were identified for plant height and seed weight per plant. Fourteen QTL for plant height were located on chromosomes 2,4,6,7,10,14,17,and 19,respectively.The LOD values of these QTL were 3.50–20.58, and the phenotypic variations explained by these QTL ranged from 2.30% to 24.18% (Table S1). Four QTL (qPH6, qPH7, qPH19-1, and qPH19-2) were detected in different treatments and years (Table 5).Eleven QTL for seed weight per plant were located on chromosomes 2, 6, 8, 13, 17, 18, and 19, respectively (Table S1).The LOD values of these QTL ranged from 3.70 to 15.87 and the phenotypic variation explained were from 4.97% to 21.33%. Two QTL (qSWPP6 and qSWPP13) were detected in different treatments and years (Table 5). The QTL that were detected in different treatments and years were considered as the stable QTL related to the target treats.Furthermore,three pairs of QTL (qPH6/qSWPP6, qPH17/qSWPP17, and qPH19–3/qSWPP19) for both plant height and seed weight per plant were located on the same chromosomal regions,so they may be related to drought tolerance. qPH6 (LOD = 5.47, PVE =10.29%)/qSWPP6 (LOD = 6.01, PVE = 9.99%) was detected not only in different treatments and years but also in different traits,which indicated that it might be a stable and major QTL controlling drought tolerance.In order to identify the additive effects of several QTL, we compared the phenotypes of the RILs that had one QTL with the phenotypes of RILs that had multiple QTL using drought DC-PH and DC-SWPP (Fig. 4). By combining qPH6 and qPH19-2, the RILs possessed qPH6 and qPH19-2 (AA genotype) had significantly higher DC-PH compared to the other combinations (Fig. 4a). By combining qSWPP19 and qSWPP17, the RILs possessed qSWPP19 (AA, AB genotypes)had higher DC-SWPP compared to the RILs that did not possess qSWPP19 (BA,BB genotypes)(Fig.4b).

Fig.3–Collinearity analysis of the genetic map and genome.The x-axis is the genetic distance of each linkage group.The yaxis is the physical length of each linkage group,with the collinearity of genomic markers and genetic maps represented in scatter.Different colors represent different chromosomes or linkage groups.

Table 4–Spearman coefficients of the chromosomes.

3.5. Candidate genes prediction

To further select candidate genes based on the Williams 82 soybean reference genome, we compared all sequence variations within nine stable and major QTL intervals located on chromosomes 2, 6, 7, 13, 17, and 19 in 800 genes. Within these intervals, we screened out 1528 SNPs and 268 InDels in 74 genes that were different between the two parents (Tables S2 and S3). Among these genes, 65 genes possessed annotated functions and located on chromosomes 2, 6, 13, 17, and 19 (Table S4). Within qPH2/qSWPP2, 22 annotated genes were identified and according to its’ functions, Glyma.02G049400, a peptide methionine sulfoxide-reductase family protein, was involved in detoxification, stress responses, and secondary metabolism and screened out for further analysis. Within qDT6/qSWPP6, 31 annotated genes were identified, 5 genes(Glyma.06G220600, Glyma.06G223800, Glyma.06G213100,Glyma.06G223500, and Glyma.06G230400) maybe be involved in stress response, which are related to drought tolerance. Glyma.06G220600 has a ubiquitin-binding domain and is involved in plant stress responses.Glyma.06G223800 is involved in potassium ion transmembrane transport. Glyma.06G213100, Glyma.06G223500, and Glyma.06G230400 are involved in both biotic and abiotic stress responses and related to abscisic acid (ABA) and auxin signaling. There were 4 genes within qSWPP13, but we didn’t find drought-related genes. Within qPH19-3/qSWPP19, Glyma.19G152200 is a leucine-rich repeat receptor-like protein kinase and responses to drought.

DC, drought resistance coefficient;WS, water stress;SW, sufficient water.

4. Discussion

Fig.4–Additive effects of several QTL.(a)Additive effects of qPH6 and qPH19-2 using DC-PH.(b)Additive effects of qSWPP19 and qSWPP17 using DC-SWPP.A and B indicate different genotypes,respectively.n,the number of RILs.Asterisks indicate significant differences as determined by ANOVA(****P < 0.001,**P <0.01,*P <0.05).

The SLAF-seq technology has been widely used for genetic map construction in the genetic researches[29–33].Compared to other sequencing technology, such as Genotyping-bysequencing (GBS) and Restriction-site associated DNA sequencing (RAD-seq), SLAF-seq markers are able to capture a large number of SNPs distributed evenly across a genome at low sequencing costs. Considering the depth of sequencing,SLAF-seq gave the best balance between sequencing cost and sequencing depth [16,26]. Using SLAF-seq, Li et al. [16]constructed a high-density genetic map in soybean containing 3541 SLAF markers with an average genetic distance of 0.72 cM. Cao et al. [18] constructed a soybean genetic map containing 3255 SLAF markers with an average genetic distance of 0.66 cM. Increase of the marker density can improve the resolution of genetic maps as well as the accuracy and quality of QTL mapping [18,34–37]. Using the SLAF-seq technology, we constructed a high-density soybean genetic map containing 8078 SLAF markers with a total genetic distance of 3780.98 cM and an average genetic distance of 0.59 cM (Table 3). Compared to previously constructed soybean genetic maps, the number of SLAF markers was greatly increased and the map distance between markers was reduced in our map, providing a high-density reference genetic map for soybean researches and could be used for gene mining easily. Based on this high-density genetic map,we investigated three drought tolerance related traits including plant height,seed weight per plant,and drought tolerance coefficient, and identified QTL for drought tolerance in different environments. Six QTL intervals (qPH2, qPH6/qSWPP6, qSWPP13, qPH17/qSWPP17, qPH19-2, and qPH19-3/qSWPP19) were also related to other traits, such as soybean branching, pod number, and seed protein and oil contents(Table S5). Some QTL for plant height and seed weight were detected on the same or close regions of the same chromosomes in different environments, indicating that the highdensity genetic map and the measurement of phenotypes were accurate and reliable.Phenotypic data showed that plant height, seed weight per plant, and drought tolerance coefficient in the three-year tests were stable under either SW or WS conditions. The drought tolerance coefficients calculated using plant height and seed weight per plant ranged from 72%–76% and 89%–94% across the three years,respectively.

Drought tolerance has become one of the major focuses of research in plant abiotic stress.SoyBase had recorded 65 QTL related to drought tolerance. Du et al. [8] used 184 F2:7:11populations from‘Kefeng 1’and‘Nannong 1138‐2’to calculate the index of yield and drought tolerance and identified 10 QTL related to drought tolerance under irrigation and drought conditions. Specht et al. [38] used 236 RILs and different indexes including carbon isotope discrimination, lodging,plant height,pod maturity,seed oil,seed protein,seed weight and seed yield, and identified a total of 55 QTL related to drought tolerance.Previous studies[39,40]have indicated that it was inaccurate to value and identify the effect of QTL in a single environment due to the instability of environment.Therefore, a QTL is considered as a stable QTL if it can be detected in different environments. A total of 23 intervals were identified and seven QTL (qPH2/qSWPP2, qPH6/qSWPP6,qPH7, qSWPP13, qPH19-1, qPH19-2, and qPH19-3/qSWPP19)conferred both plant height and seed weight per plant under different environments. Among the stable QTL, four for plant height (qPH6, qPH7, qPH19-1, and qPH19-2) were detected in different environments, and three of them (qPH6, qPH19-2,and qPH19-3) were also reported to be associated with plant height in previous studies [41–53]. Two QTL for seed weight per plant (qSWPP6 and qSWPP13-1) were detected in different environments and qSWPP6 was reported to control seed weight [54–58]. More importantly, qPH6/qSWPP6 was the common QTL conferring plant height and seed weight per plant in different environments.This common QTL controlled unrelated traits might be due to pleiotropy of the locus. qPH6 on chromosome 6 overlapped with a previously reported QTL.In our study, the position of qPH6 on chromosome 6 was 19,480,011–38,332,813. Specht et al. [38] and Du et al. [8] had identified the same locus as qPH6.Its physical positions were 18,320,494–23,848,501 and 35,541,410–38,332,813, respectively.Interestingly, three novel QTL (qPH2, qPH17/qSWPP17, and qPH19-2) were related to plant height and four novel QTL(qSWPP2, qSWPP13, qSWPP17, and qSWPP19) were related to seed weight. Therefore, these major and stable QTL were reliable and could provide a reference for the molecular marker-assisted selection and could be used for further fine mapping genes for drought tolerance. Furthermore, the additive effects between qPH6 and qPH19-2 showed that the RILs possessed qPH6 and qPH19-2 (AA genotype) simultaneously had significantly higher DC-PH compared to the other combinations(Fig.4a).The additive effects between qSWPP19 and qSWPP17 showed that the RILs with qSWPP19 (AA, AB genotypes)had higher DC-SWPP than those without qSWPP19(BA, BB genotypes) (Fig. 4b). Therefore, we speculated qSWPP19 might be a major QTL that controls drought tolerance.

Drought tolerance in plants is very complex controlling by complicated networks of multiple genes. Based on these major and stable QTL detected in different environments, by annotating the genes within the 10 QTL intervals against the publically available databases, 74 protein-encoding genes differed in sequences between the two parents, of which 65 genes have GO annotations (Table S4). Although some SNPs were in untranslated regions (UTRs), many studies have reported that the UTRs could play regulatory roles in RNA translation, RNA stability, and RNA transcription. Mutations in the UTRs could affect gene expression levels and control the phenotype in mutants. Besides, within the QTL interval,many loci are tightly linked. Therefore, these genes possessing SNPs in the UTRs or causing nonsynonymous mutations might cause critical mutations that are related to our target traits. Because SLAF-seq is a selective sequence technique and couldn't cover the entire sites,the locations of the whole mutations are not yet shown for further experimental verification.

Several functional genes are involved in osmotic stress,salt stress, and water deprivation, which are related to drought tolerance. Using a restricted two-stage, multi-loci,genome-wide association study (RTM-GWAS) of relative root length (RRL) and relative shoot length (RSL), Khan et al. [10]identified 73 and 38 QTL with 174 and 88 alleles,respectively.The QTL explained 40.43%and 26.11%of the variation in these traits, and the QTL-environment interaction (QEI) explained 24.64%and 10.35%of the phenotypic variance of RRL and RSL,respectively. Khan et al. [10] also identified 134 candidate annotated genes, including 88 genes harboring SNPs associated with RRL and 46 genes harboring SNPs associated with RSL in the NAM population. Among the 74 candidate genes identified, Glyma.02G049400 was also detected and there was a nonsynonymous mutation (G/T) between the two parents.Glyma.02G049400 was a peptide methionine sulfoxidereductase family protein that was involved in the cellular protein modification process, oxidation-reduction process and protein metabolic process. Oxides, such as hydrogen peroxide and oxylipins, can regulate plant development,stress adaptation, and programmed cell death. Low oxide can respond to environmental and developmental stress,and high oxide triggers cell death[59].In Arabidopsis,genes related to OPDA REDUCTASE3 (OPR3) were involved in detoxification,stress responses,and secondary metabolism[60].

Glyma.06G220600 is a ubiquitin-binding domain involved in plant stress responses. Ubiquitin-conjugating enzyme accumulation could improve salt tolerance by activating the BR signaling pathway in Arabidopsis[61–64].Transforming E2 proteins of mung bean (Vigna radiata (L.) Wilczek.) and soybean into Arabidopsis can enhance the resistance to salt and drought stress[63,64].Although we didn’t detect any SNP/InDels located at the gene coding region in Glyma.06G220600,there were three SNPs in the 3′UTR. It was found that the ionotropic glutamate receptor was a type of cation channel protein and responded to a variety of biotic and abiotic stresses by different pathways,including plant light signaling,apical meristem activity,pollen tube growth,and cytoplasmic Ca2+flow [65]. Glyma.06G223800 was annotated as an ion transport protein and K+channel related protein involved in potassium ion transmembrane transport.A base InDel(G)was detected in the 5′UTR of this gene between the two parents.Srivastava et al. [66] found that mRNA variation in the UTRs increased the coding capacity of the genome by producing multiple types of transcripts to respond to stimuli from plant growth and environmental stress [66]. Therefore, genes with mutations in the UTR regions might control drought tolerance directly through gene expression or indirectly regulate drought tolerance in plants.

Glyma.06G213100, Glyma.06G223500, and Glyma.06G230400 are involved in specific aspects of abscisic acid (ABA) and auxin signaling. Increased anion channel activation and ABA treatment were reported to reduce water loss during drought stress [67]. Johnson et al. [68] identified geranylgeranyl transferase I to be involved in specific aspects of abscisic acid(ABA)and auxin signaling in Arabidopsis.ABA is involved in both biotic and abiotic stress responses and regarded as a stress hormone. Under stress conditions, the phytohormone ABA regulates many pathways including stomatal closure,and also induces expression of stress-related genes[68].

It was found that the ER genes (SbER1-1 and SbER2-1)from sorghum (Sorghum bicolor L.), a drought-tolerant model plant,could increase drought tolerance in Arabidopsis and maize(Zea mays L.) by increasing water-use efficiency and photosynthetic rate under drought stress. These ER genes encode leucine-rich repeat-receptor-like kinases [69]. In our study,Glyma.19G152200 was also annotated as leucine-rich repeat receptor-like protein kinase. More importantly,Glyma.19G152200 has a nonsynonymous SNP (T/G) between the two parents.

The candidate genes identified in this study may be related to plant stress, so they can be used for further research on soybean drought tolerance.However,the QTL should be finely mapped,and the functions of these genes need to be verified.The identified QTL related to drought tolerance are valuable resources for fine mapping, gene cloning, and functional analysis of drought tolerance in soybean, which provides excellent genetic resources for breeders to develop drought tolerant soybean varieties.

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Acknowledgments

This research was supported by the National Key Research and Development Program of China(2016YFD0100201 and 2016YFD0100304) and the National Science and Technological Innovation Program of China.

Author contributions

Honglei Ren and Xingrong Wang conducted the RIL populations and phenotyping. Honglei Ren and Jianan Han conducted the data analysis,QTL mapping,genomic comparative analysis and wrote the manuscript. Lili Yu, Huawei Gao,Huilong Hong, Rujian Sun, Yu Tian and Xusheng Qi provided advice on experimental implementation. Lijuan Qiu, Xiaoxia Wu and Zhangxiong Liu conceived and supervised the project,and reviewed the manuscript. Lijuan Qiu and Bo Zhang revised the manuscript.

Appendix A. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.04.004.