Identification of genes for drought resistance and prediction of gene candidates in soybean seedlings based on linkage and association mapping

Ynjun Zhng,Zhngxiong Liu,Xingrong Wng,Yue Li,Yongsheng Li,Zuowng Gou,Xingzhen Zho,Huilong Hong,Honglei Ren,Xusheng Qi,*,Lijun Qiu,*

a Institute of Crop,Gansu Academy of Agricultural Sciences,Lanzhou 730070,Gansu,China

b National Key Facility for Gene Resources and Genetic Improvement/Key Laboratory of Crop Germplasm Utilization,Ministry of Agriculture/Institute of Crop Sciences,Chinese Academy of Agricultural Sciences,Beijing 100081,China

c Maize Research Institute,Heilongjiang Academy of Agricultural Sciences,Harbin 150086,Heilongjiang,China

Keywords:Candidate genes Drought tolerance Glycine max QTL SNP

ABSTRACT Drought is one of the primary abiotic stress factors affecting the yield,growth,and development of soybeans.In extreme cases,drought can reduce yield by more than 50%.The seedling stage is an important determinant of soybean growth:the number and vigor of seedlings will affect growth and yield at harvest.Therefore,it is important to study the drought resistance of soybean seedlings.In this study,a recombinant inbred line (RIL) population comprising 234 F6:10 lines (derived from Zhonghuang 35 Jindou 21) and a panel of 259 soybean accessions was subjected to drought conditions to identify the effects on phenotypic traits under these conditions.Using a genetic map constructed by single nucleotide polymorphism (SNPs) markers,18 quantitative trait loci (QTL) on 7 soybean chromosomes were identified in two environments.This included 9 QTL clusters identified in the RIL population.Fifty-three QTL were identified in 19 soybean chromosomes by genome-wide association analysis(GWAS) in the panel of accessions,including 69 significant SNPs (-log10 (P) ≥3.97).A combination of the two populations revealed that two SNPs (-log10 (P) ≥3.0) fell within two of the QTL (qPH7-4 and qPH7-6) confidence intervals.We not only re-located several previously reported drought-resistance genes in soybean and other crops but also identified several non-synonymous stress-related mutation site differences between the two parents,involving Glyma.07g093000,Glyma.07g093200,Glyma.07g094100 and Glyma.07g094200.One previously unreported new gene related to drought stress,Glyma.07g094200,was found by regional association analysis.The significant SNP CHR7-17619(G/T)was within an exon of the Glyma.07g094200 gene.In the RIL population,the DSP value of the ‘‘T” allele of CHR7-17619 was significantly (P ＜0.05) larger than the ‘‘G” allele in different environments.The results of our study lay the groundwork for cloning and molecular marker-assisted selection of droughtresistance genes in soybeans at the seedling stage.2021 Crop Science Society of China and Institute of Crop Science,CAAS.Production and hosting by Elsevier B.V.on behalf of KeAi Communications Co.,Ltd.This is an open access article under the CC BY-NCND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1.Introduction

Drought is the most important abiotic constraint to crop production worldwide.With increasing effects of climate change the frequency and duration of droughts are increasing at an alarming rate,causing huge agricultural losses every year[1].Soybeans(Glycine max) are a globally important crop due to their value in producing oil,plant protein,and food [2].Growing soybeans requires significant amounts of water.Since the root system is not well-developed,soybeans are sensitive to drought conditions.As such,drought can affect the yield and quality of soybeans across the world [3–5].

Drought continues to be a threat to soybean production around the world [6].The seedling stage is the basis of soybean growth.The intensity of drought stress determines the number and health of seedlings after emergence.In addition,the identification of drought resistance at the seedling stage is cheap,easy to perform and requires a shorter experimental cycle.Therefore,it is useful to study drought resistance at the soybean seedling stage [7,8].Since 2007,we have screened many accessions and analyzed several morphological and physiological indices related to drought resistance in soybeans and identified a number of droughtresistance germplasms.After years of multi-point identification,we found that drought stress seriously affected soybean yield and that the yield of many accessions under drought stress was decreased by 40%–70%[9–12].This could have a devastating effect on soybean production.Therefore,it is important to breed drought-resistance varieties and study the mechanisms of drought resistance of soybean crops as a contribution to global food security.

Traditional research on drought resistance has focused on selection using drought resistance indices,evaluation methods,and identification of drought-resistance germplasm [13].In recent years,molecular biological methods such as the QTL mapping,genome-wide association analyses (GWAS),and structural and functional genomics have been widely used to analyze the mechanisms underlying drought resistance [1].The inheritance of quantitative traits has been described as a ‘moving target’ since these traits are largely influenced not only by the actions of multiple QTL/genes,but also by interactions between QTL/genes (epistasis)and between QTL/genes and environmental factors [14].Linkage and association analyses are the most common methods used to dissect complex traits.Linkage analysis in hybrid populations involves two specific parents,meaning that there are only two alleles at any locus.This decreases the potential value of QTL mapping.Association analysis is based on linkage disequilibrium (LD) and can identify relationships between groups of target traits with complex backgrounds,genetic markers,or candidate genes.In addition,when the population used in association analysis is small linkage analysis can be used to corroborate the findings.Therefore,linkage analysis and association analysis were used in this study to reduce analysis errors.

QTL for drought tolerance in soybeans were analyzed using the linkage analysis [15–24] of traits such as water use efficiency[25],yield [16] and root traits [26].Kaler et al.[17] found 61 significant SNPs for soybean canopy wilting tolerance.Khan et al.[15] used restricted two-stage multi-locus genome-wide association studies (RTM-GWAS) to identify 111 drought-tolerant QTL and 262 alleles with a high proportion of QEI (QTL-environment interaction) and genetic variation.These accounted for 88.55%–95.92% of the phenotypic variance in nested association mapping(NAM),from which QTL-allele matrices were established and candidate genes annotated.However,no study has yet combined two mapping (linkage and association) strategies to perform a genetic analysis of quantitative traits in relation to drought resistance in soybean.

Many QTL or SNPs related to agronomic and physiological traits of soybeans (www.soybase.org) under abiotic and biotic stress have been identified;however,few related to soybean drought resistance have been identified [27].Even fewer are related to drought resistance in soybean at the seedling stage.In this study,we investigated the drought survival proportion (DSP),plant height (PH),root length (RL),seedling fresh weight (SFW),root fresh weight (RFW),seedling dry weight (SDW),and root dry weight (RDW) of soybean seedlings subjected to drought stress conditions in two environments.We used 8078 SNPs and 4616 SNPs to perform linkage and association analyses and mapped the drought resistance QTL/SNPs identified in soybean seedlings.We used two methods to identify stable QTL/SNPs for drought resistance and to identify candidate genes related to seedling drought resistance.

2.Materials and methods

2.1.Population analysis

We used 234 F6:10recombinant inbred lines (RIL) derived from Zhonghuang 35 Jindou 21 in linkage analysis.Zhonghuang 35 is a drought-sensitive,high-yielding cultivar,whereas Jindou 21 is strongly drought resistant at the adult stage.Comparison of the two cultivars under drought at the seedling stage revealed significant phenotypic differences.Additionally,there were significant differences in other agronomic traits,such as plant height and the growth period between two parents at the adult stage.The association panel consisted of 259 diverse cultivars from 11 provinces in China[28].According to seven phenotype values at seedling stage,70 cultivars (out of 259) were used to construct the regional association analysis group for subsequent regional associations,because these can represent all cultivars well,and their resequencing has been completed (Table S1).

2.2.Phenotype identification and statistical analysis

The drought resistance assays were performed in Dunhuang and Zhangye,Gansu province,from 2018 to 2019.In a droughtresistance shed (Fig.1),a meshed plastic turnover box(45 30 25 cm)was placed into a water tank that allowed drainage.The box was filled with 20 cm of thick,arable soil of medium uniform fertility,and the water was drained once after saturation of the soil [29].The RIL population and association panel were planted in randomized complete block design in three replications.Each row contained one cultivar/accession(20 seedlings)and each box contained 3 rows.

For the first drought stress rehydration treatment we stopped supplying water at the V2 stage when drought stress began.When the seedlings appeared permanently wilted(absolute soil moisture content,5.9% ± 0.5%),the water content of the soil was measured,and rehydration began[29].On the third day after rehydration,we examined the number of surviving seedlings and recorded them as having survived if the leaves of the seedlings turned green.In the second drought stress rehydration treatment,we stopped the water supply after the first rehydration,at which point the second drought stress began.When the seedlings were permanently wilted (absolute soil moisture content 3.2% ± 0.5%),we measured the soil water content and began rehydration.We then examined the number of surviving seedlings on the third day after rehydration,recording them as having survived if the leaves of the seedlings turned green.

The drought survival proportion of seedlings under repeated drought conditions [7]:

where,DSP represents the surviving proportion of seedlings under repeated drought,DSP1 represents the surviving proportion from the first drought treatment,DSP2 represents the surviving proportion from the second drought treatment,represents the average of the total number of seedlings in the three replications before the first drought,represents the average number of seedlings in the three replications after the first rehydration,andrepresents the average number of seedlings in the three replications after the second rehydration.

After the second drought stress rehydration test,five plants with uniform growth potential were selected from each cultivar/accession to measure the following traits.Plant height (PH) and root length (RL) were measured using the straightedge method.Seedling fresh weight (SFW) and root fresh weight (RFW) were measured using an electronic balance.After the fresh weights were measured,the dry weights (SDW and RDW) were measured after being placed in Kraft paper bags and held at 121 C for 2 h and followed by 80 C for 24 h.Average values were used in subsequent data analysis.

Fig.1.Conditions for assaying drought resistance in seedling stage soybeans.

SPSS 17.0 and R 3.6.1 were used to analyze the phenotypic data.Broad-sense heritability of each trait was calculated using the formula [30].

The best linear unbiased prediction(BLUP)under the mixed linear model was performed for each trait value by the following formula using the R package ‘‘lme4”;the resulting BLUP value was used for QTL mapping and GWAS analysis [31,32].

where Pikrefers to the observed phenotype of the kthaccession in the ithyear;μ is the mean value of the trait;Yiis the fixed effect of the ithyear,gkis the random effect of the k genotype;and eikis the residuals of the model.

2.3.Genotyping

In previous research,we used SLAF-seq (Specific-Locus Amplified Fragment Sequencing) technology,independently developed by Beijing Biomarker Technologies Co.,Ltd.,to develop highdensity molecular tags for RILs in a segregating population(2 parents and 232 lines)and a panel of accessions (259 cultivars).After an analysis of polymorphism we identified 8078 polymorphic SNP markers in the RIL population [33] and 4616 polymorphic SNP markers in the accession panel [28].Total genomic DNA of the regional association analysis groups were extracted from young leaf tissue using a standard cetyltrimethylammonium bromide(CTAB) method [34].We performed tunable genotyping by sequencing,with one base pair of selectivity,using the Life Technology Ion Proton Systems by Data2Bio LLC,according to methods previously described [35].Each individual sequence read was scanned and trimmed for regions of low-quality sequence(defined as having a PHRED quality score ＜15).Trimming was conducted in two steps:(i) low-quality nucleotides of each read end were removed;(ii) remaining nucleotides were then scanned using overlapping windows of 10 bp,and sequences beyond the last window with average quality values less than PHRED 15 were truncated.Trimmed reads were aligned to the Glycine max Wm82.a2.v1 reference genome using GSNAP[36,37].Each SNP was required to have a minor allele frequency(MAF)of ≥0.05,a missing rate ≤70%,and a heterozygosity ≤20%.From this,we obtained a highquality set of 2473 SNPs.The sequences obtained were used in subsequent analysis (Table S2).

2.4.QTL mapping

We used JoinMap Version 3.0 (http://www.joinmap.nl) mapping software [38] and Kosambi mapping function to construct a local molecular linkage map and calculate genetic distances.We performed the QTL mapping analysis using ICIM software Icimapping v4.1 (http://www.isbreeding.net/software/?Type=detail&id=18).QTL mapping was performed separately with the average value of replications in a single environment.To co-localize genetic loci in different genetic backgrounds,the LOD threshold was set to 2.0[39,40].If the LOD value was ＞2.0 the LOD value corresponded to a drought-resistance QTL.The detected QTL were named according to the English abbreviation of the corresponding chromosome number and the rank of QTL on the chromosome.For example,qDSP1-1 represents the first QTL on the first chromosome.

2.5.Linkage disequilibrium and kinship analysis

Tassel 5.2.13 software[41]was used to analyze the linkage disequilibrium,and the LD attenuation window size was set to 100,meaning we calculated the linkage imbalance between 100 pairs of markers before and after each marker.When P ＜0.001,there was significant LD [42].We used the kinship option in Tassel 5.2.13 software to analyze the kinship.

2.6.Genome-wide association and regional association analysis

We used the general linear model (GLM) and the mixed linear model (MLM) in Tassel 5.2.13 software for genome-wide association analysis[43].We then determined the optimal model according to the Quantile-Quantile plot (Q-Q plot) of the BLUP values of seven traits,after which the Manhattan map was drawn.The Q value of the group structure was obtained using Structure 2.3.4 software analysis,and the K value and PCA value of kinship were obtained using Tassel 5.2.13 software [28].Significant association threshold values were set at -log10(0.5/n) (n,total filtered SNP number).This was used to determine whether there were SNP sites related to drought resistance in soybean seedlings.The general linear model (GLM) in Tassel 5.2.13 software for the regional association analysis used the Q model.Significant association threshold values were set at -log10(0.5/n) (n,total filtered SNP number).Associated QTL were defined across regions of SNPs with LD value(r2＞0.2) between the significant association and surrounding SNPs,and adjacent SNPs were within a physical distance of 3 Mb[44–45].

2.7.Candidate gene prediction

Along with QTL co-located by the linkage analysis and association analysis,an LD analysis of SNPs in the co-located QTL was performed with Haploview v.4.2 software [46] to find existing LD blocks.Using the genome sequence of Wm82.a2.v1 (https://www.soybase.org/),we found the drought-resistance genes in the marker regions on both sides of the blocks and annotated the candidate genes.

2.8.Isolation and bioinformatic analysis of protein sequences

The protein sequences of non-synonymous genes between the two parents were downloaded from SoyBase (https://www.soybase.org/).Two sequence alignments of the non-synonymous mutant genes and the reference genes were conducted using DNAMAN V6 software (https://www.lynnon.com/).Protein domains were predicted using the CD search tool on the National Center for Biotechnology Information database (NCBI,https://www.ncbi.nlm.nih.gov).Protein secondary structure was predicted using the SOMPA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?-page=npsa_sopma.html) online prediction program [47].Protein tertiary structure homology modeling was constructed using SWISS-MODEL (https://swissmodel.expasy.org/interactive) online software [48] and Swiss-Pdb Viewer software [49].

2.9.Specific PCR marker development and candidate gene validation

Total genomic DNA the RIL were extracted from young leaf tissues using the CTAB method [34].Cleaved amplified polymorphic sequence (dCAPS) markers used for PCR amplification were designed through dCAPS Finder 2.0 (http://helix.wustl.edu/dcaps/dcaps.html),the sequences for which were:F:TTACAACAAGGGACTTGAAGCTACC,R:TCTTGGGGCAAAAGTTCTCGT.PCR amplification products of the RIL population were cut with enzyme AciI(NEW ENGLAND Biolabs (NEB)) and analyzed by polyacrylamide gel electrophoresis.The RIL population was classified into three groups based on the selected SNP.A t-test was applied to detect the significance of the phenotypic difference between the groups.

3.Results

3.1.Changes in soil water content after drought stress

Before drought stress treatment,the soil water content was kept within a range required for normal soybean growth to ensure that it was not affected by drought stress.In 2018,the soil water contents at Dunhuang after each of two drought stresses were 5.88% and 3.28%.In 2019,the soil water contents at Zhangye after two drought stresses were 5.95%and 3.14%(Fig.S1).The soil water contents did not differ significantly between the two environments following both treatments of drought stress.This indicated that trait differences in seedlings following drought stress were mainly related to genotype rather than environment.

3.2.Statistical analysis of phenotypic differences after drought stress

In 2018–2019,the drought resistance of soybean seedlings was studied by evaluating seven traits.Statistical analysis showed that difference between the seven traits were significant (Table S3).Some cultivars died and others survived following drought stress.The coefficient of variation for each trait in different environments ranged from 13% to 103%,indicating that the genetic variation in the two populations was relatively large.

According to normality tests the drought stress traits showed continuous variation;most were normally distributed slightly skewed to the right.Thus,following drought stress,seven traits(such as plant height) were quantitatively controlled by multiple genes(Fig.S2)with significant difference a between the genotypes in each population(Table S4).Broad-sense heritability ranged from 0.07 to 0.94(Table S3).Correlation analysis showed highly significant positive correlations between the seven traits (P ＜ 0.01,Table S5).This demonstrated that the seven traits reflected the drought response of soybean seedlings and indicated that the investigated traits were reliable for the genetic analysis of seedling stage drought resistance.Therefore,the DSP,PH,RL,SFW,RFW,SDW,RDW,and BLUP values were used in the subsequent genome-wide association analysis.

3.3.QTL mapping of drought resistance in soybean seedlings

Based on the linkage map using 8078 SNPs we detected 46 QTL(Table S6) using interval mapping (IM) and inclusive complete interval mapping (ICIM) in the two environments.These were located on soybean chromosomes 2,6,7,9,11,12,and 17.The LOD scores ranged from 2.02 to 3.16,and the QTL explained 2.82% to 6.40% of the corresponding phenotypic variance,respectively.Two different mapping methods were used to detect the twelve primary QTL that included qSDW2-1,qRDW6-1,qRL6-3,qSDW6-3,qSFW6-3,qDSP6-4,qDSP6-5,qRFW7-2,qRL7-2,qRL11-1,qDSP12-2,and qRDW17-1.Combining QTL located at the same loci using IM and ICIM we mapped 18 QTL in soybean seedlings following drought stress.The additive effects of the 18 QTL ranged from-0.30 to 6.16.Twelve positive values indicated that the alleles for increasing drought-resistance phenotypic values came from resistant parent Zhonghuang 35,indicating that QTL mapped at the seedling stage could improve the drought resistance of soybeans.There were nine QTL clusters on chromosomes 2,6,7,9,11,12,and 17 (Fig.2).There were two clusters on each of chromosomes 6 and 7.Six QTL were in the same cluster on chromosome 6 and four and five QTL were in separate clusters on chromosome 7.These results indicated a high possibility of major QTL or candidate genes in these clusters.

3.4.Genome-wide association analysis

Using Tassel 5.2.13 software,we analyzed the relationships between association population and calculated the relationship between any two cultivars,and then graphed these relationships(Fig.S3A).About 55% of cultivars had a relationship value of 0 and 9%had a relationship value of 0–0.05(Fig.S3B).This indicated that the relationships members were distant,thus meeting the requirements of GWAS.

At r2=0.2,the attenuation distance of the soybean genome is 3 Mb and its deceleration rate is slow (Fig.S4).This could be related to the span of evolutionary time and artificial selection.However,the status of soybean as a self-pollinating crop has resulted in high homozygosity and large LD.

Fig.2.Chromosome locations of QTL for drought resistance in the RIL population.Color intensity of the bar graph represents the genetic distance between markers:red indicates a small distance,and blue indicates a large distance.Red italics represent co-detected QTL with association analysis.C1 to C20,chromosome 1 to chromosome 20.

Using the BLUP values of the seven traits,we performed an association analysis using the general GLM,Q,and PCA linear models and K,Q+K,and PCA+K models for mixed linear modules using Tassel 5.2.13 software.By comparing the distribution of the Q-Q map under different models,the K,K+Q,and K+PCA models under the MLM model can better control false positives than the GLM,Q,and PCA models under the GLM model.The P-value detected by the K+PCA model was closer to the expected value,so the K+PCA model was selected as the best method for identifying drought-resistance loci in soybean seedlings.In the MLM analysis model,the PCA matrix was used as a fixed effect and the K matrix of kinship was used as a random effect to control false positives during our GWAS and the drawing of our Manhattan diagram and Q-Q plot (Fig.S5).

At -log10(P) ≥3.97 we detected 121 significant SNPs in different environments;71 in a single environment in 2018–2019.Fifteen SNPs were related to RDW,53 to RL,and 3 to RFW.They were distributed on all chromosomes except chromosome 9 and explained 7.39% to 30.05% of the phenotypic variances.Fifty SNPs were detected by the BLUP values for different traits,including 46 related to RL and 4 related to RDW.These were mapped on 17 chromosomes (except for 2,9 and 12) and accounted for 4.64% to 25.93% of the phenotypic variance.Excluding SNPs at the same locus for different traits we detected 69 significant SNPs.These were delineated into 53 QTL for three measured traits(Table S7).QTL qRL3-3,qRDW6-8,qRL7-8,qRL13-1,qRL14-1,qRL16-1,qRL16-2,qRFW20-1,and qRL20-2 contained 3,4,2,2,3,2,2,2,and 5 SNPs,respectively.The other QTL contained only one associated SNP.

3.5.Identification of candidate genes

Generally,genetic loci co-localized in different genetic backgrounds are thought to have stable effects on phenotype.Therefore,we also focused on loci that were detected in both populations.According to a previous study [50],we lowered the threshold of LOD to 2.0 and -log10(P) ≥3.0 to identify stable loci across the two populations.By comparison,the significant SNP sites ss246200757 and ss246238768 detected by RDW-2019 and RDW-BLUP were in the QTL (qPH7-4 and qPH7-6) in chromosome 07(Table 1).The physical intervals of QTL qPH7-4 and qPH7-6 were 489.69 kb and 264.26 kb and accounted for 4.69%and 6.40%of the corresponding phenotypic variances,respectively. SNPs ss246200757(RDW-2019)and ss246238768(RDW-BLUP)explained 6.49% and 6.01% of the respective phenotypic variances,respectively.This could be the primary segment of the gene conferring to drought resistance in soybean seedlings.

We used Haploview v.4.2 software to perform an LD analysis of the significant SNPs ss246200757 and ss246238768 extending for 300 kb on the left and right sides of the SNPs.Three LD blocks were found (Fig.3) and we chose blocks 1 (6,219,969–6,339,218) and 2(8,669,670–8,864,954) in the physical intervals of qPH7-4 and qPH7-6 to predict candidate genes.

The reference genome sequence Wm82.a2.v1 (https://www.soybase.org/) has 35 genes in the two blocks (Table S8).For the 35 genes we identified 48 SNP and 2 InDels differences between the two parents(Table S9).Among the 35 candidate genes there were eight non-synonymous differences in Glyma.07g093000,Glyma.07g093200,Glyma.07g094100,and Glyma.07g094200.

Table 1 Summary of association and linkage mapping results related to plant height and root dry weight.

3.6.Sequence analysis of genes with non-synonymous mutations

The amino acid sequence in the Glyma.07g093200 at 8,699,308 bp(T/C)was not altered by the non-synonymous mutation.Four non-synonymous mutations in Glyma.07g093000 caused four amino acid changes from proline,serine,threonine and alanine to leucine,leucine,arginine,and serine,respectively(Fig.S6A).One non-synonymous mutation (8,788,120 bp) in Glyma.07g094100 caused an amino acid change of asparagine to histidine.Another,at 8,787,829 bp caused no amino acid change(Fig.S6B).A non-synonymous mutation at 8,797,367 bp in Glyma.07g094200 caused a change from threonine to alanine(Fig.S6C).The protein domain was unaltered after the amino acid changes in three genes that belonged to the protein kinase superfamily (Fig.S7).

Fig.3.LD analysis of co-detected SNPs on chromosome 07.(A) The physical intervals and LOD scores of QTL qPH7-4 and qPH7-6.(B) LD heat map of surrounding the significant SNPs.Red dashed lines indicate the QTL intervals.

The secondary structures of the proteins encoded by the three genes were mainly random coils and alpha helixes,and the proportions of extended strand and beta region structure were low(Table S10).The number of amino acids in the beta turn and random coil did not change after the mutations in Glyma.07g094100(Fig.S6E).The secondary structures of Glyma.07g093000 (Fig.S6D)and Glyma.07g094200 (Fig.S6F) both changed.The amino acid sequences of the three genes were not matched to highly homologous sequences,and the tertiary structures had a large proportion of areas that were not covered.Therefore,the amino acid at positions 359 and 242 does not show Glyma.07g093000(Fig.S8A) and Glyma.07g094100 (Fig.S8B),respectively.We constructed the protein tertiary structure of the mutant sequence according to the sequence prior to the mutation.After the mutation,the tertiary structures of the Glyma.07g093000 and Glyma.07g094200 proteins (Fig.S8C) changed to varying degrees.These changes might be caused by the van der Waal and helix folding structures.

3.7.Further screening of candidate genes

In order to further narrow the range of candidate genes,2473 SNPs obtained by re-sequencing the QTL between qPH7-4 and qPH7-6 were used for regional association analysis.By comparing the GLM models,we found that the P-value was closer to the expected value,so GLM was selected as the best model for regional correlation analysis.Using a general linear model,the strict significant P-value was set at 2.0 10-4(0.5/2473).Five associations in drought-resistance traits were identified based on the resequencing of subsets of accessions (Table 2).We also identified one SNP (CHR7-17619 at 8796932) located within exons of Glyma.07g094200 (8792240–8797837) that was associated with thevariation in DSP-2018(Fig.4)and accounted for 25.62%of the phenotypic variance.This could be new drought-resistance gene.

Table 2 Significant SNPs for all drought-resistance traits identified by regional association analysis.

Fig.4.Identification of the significant SNP CHR7-17619 within exons of Glyma.07g094200.(A) Gene model of Glyma.07g094200. (B) Manhattan plot of DSP-2018 surrounding CHR7-17619.Area between the black dashed lines indicates the physical interval of Glyma.07g094200.Red dots indicate the significant SNP CHR7-17619.

3.8.Verification of candidate genes

Through sequencing and comparative analysis we found that the SNP CHR7-17619 was a ‘‘T” and ‘‘G” change in Zhonghuang 35 and Jindou 21,respectively,and ‘‘T” and ‘‘G” showed positive and negative effects on DSP,respectively.PCR amplification and polyacrylamide gel electrophoresis showed that 105 lines were ‘‘TT”,94 lines were ‘‘GG”,and 29 lines were ‘‘GT” in the RIL population.A t-test showed that the DSP value of the lines containing the‘‘TT” allele was significantly (P ＜0.05) larger than that of the lines containing the‘‘GG”allele in each environment,and the DSP value of ‘‘GT” was between ‘‘GG” and ‘‘TT” (Fig.5).

4.Discussion

Soybeans are sensitive to lack of water and drought resistance is related to the period,intensity,and duration of drought,as well as cultivar genotype.Soybeans are significantly affected by water during seed germination and early growth,which are important for crop establishment and therefore for final yield potential.Drought stress of seedlings leads to fewer and shorter roots and low survival rates.Symbiotic nitrogen fixation in legumes is also sensitive to water deficit and salinity stress [51,52].

Drought resistance of crops is a quantitative trait controlled by multiple genes and is greatly affected by environments.This makes phenotyping extremely important when studying drought resistance in crops.Identifying drought resistance at the seedling stage is typically performed in the laboratory using a hypertonic solution treatment and the repeated drought method.However,only partially simulates conditions in the field.This study improved the repeated drought method [7,53] by identifying new techniques(Fig.1)that can be performed in a simple drought-resistance shed.The soil was from the field,and the drought-resistance shed was closed on rainy days to prevent rainfall from affecting the experiment.The conditions in the testing environment were similar to those in the field,providing novel and accurate phenotypic identification.

Seedling height,root length,seedling dry weight,and root dry weight were the parameters considered closely related drought resistance in soybean seedlings [54,55].In this study,there were significantly different genotypic response in all seven traits after drought stress with coefficients of variation ranging from 13% to 103%.DSP directly reflects the drought resistance of soybean seedlings.Drought-resistant genotypes have high DSP in contrast to drought-sensitive genotypes.Correlation analyses demonstrated that the DSP,PH,RL,SFW,RFW,SDW,and RDW were significantly positively correlated (P ＜0.01),indicating that all seven traits are affected by drought stress.

Drought resistance is a complex quantitative trait controlled by multiple genes/QTL and is significantly affected by the environment.Molecular marker-assisted selection using QTL for drought resistance is a particularly effective method.Linkage mapping and association mapping are common procedures for dissecting the genetic architecture of complex traits.The two mapping strategies are complementary in accuracy and breadth of mapping,amount of information obtained,and methods of statistical analysis.Combining the two strategies greatly improves the value of genetic analysis of complex quantitative traits [56–58].

To identify stable genetic loci,we focused on QTL located by both procedures.We found that the significant SNP loci ss246200757 and ss246238768 related to RDW were in the QTL qPH7-4 (5,949,323–6,439,015 bp) and qPH7-6 (8,648,912–8,913,1 72 bp),respectively.There are currently few reports on the location of drought resistance sites in soybean seedlings.Khan et al.[15]identified 73 and 38 QTL with 174 and 88 SNPs,which accounted for 40.43% and 26.11% of phenotypic variance (PV),respectively,and QTL-environment interaction (QEI) effects accounted for 24.64%and 10.35%of PV for relative root length(RRL)and relative shoot length (RSL),respectively.One hundred and thirty-four annotated candidate genes were involved in nine biological processes.QTL RSL2.1 (position:4,253,109–4,444,638) and RSL7.1(3,637,584–3,637,602) overlapped with qSFW2-1 (4,440,872–4,878,598) and qRL7-2(3,250,143–3,908,906).However,the two QTL co-localized in the present study were not reported by Khan et al.[15].These major and stable QTL could provide targets for marker-assisted selection.

Fig.5.Mean effects of SNP CHR7-17619 genotype on DSP for the RIL population compared by t-test.Violin plots of the DSP-2018(A),DSP-2019(B)and DSP-BLUP(C).*and**indicate significance at P ＜0.05 and P ＜0.01,respectively.DSP,proportions of drought survival;2018,year 2018;2019,year 2018;BLUP,best linear unbiased prediction.

There were 35 genes in two blocks with 17 having at least one mutant site difference between the two parents.Of them,four genes (Glyma.07g093000,Glyma.07g093200,Glyma.07g094100,and Glyma.07g094200) had eight non-synonymous mutations.Seven genes(Glyma.07g069000,Glyma.07g069300,Glyma.07g069500,Glyma.07g069700,Glyma.07g093200,Glyma.07g094000,and Glyma.07g094800) were related to drought resistance.Block 1 contained four drought-resistance genes.We found that Glyma.07g069000 was located in the 6,232,846–6,235,874 bp interval and encoded a receptor-like kinase (RLK) that improves plant performance under drought stress [59],primarily through regulation of the ABA signal transduction pathway [60].Glyma.07g069000 is involved in drought stress response in Arabidopsis [61],maize[62],wheat[63]and other crops,but has not been reported in soybean.Glyma.07g069300,identified in the 6,262,886–6,265,779 bp interval,encodes a ubiquitin regulating enzyme (UBC).After six hours of treatment with 20% polyethylene glycol (PEG),GmUBC2 transcripts were significantly increased compared with the control.The GmUBC2 gene is expressed in soybean plant tissues and its overexpression improves drought resistance in soybean and Arabidopsis [64].Glyma.07g069500,located in the 6,304,097–6,306,91 0 bp interval,encodes a cytochrome P450 (CYP).There are 10 CYP707A gene families in soybean.Expression of CYP family genes can improve drought resistance of crop plants [65,66].Ten genes from the CYP707A gene family were expressed in different soybean tissues following drought stress,indicating that Glyma.07g069500 might participate in the response of soybean seedlings to drought stress [67].Glyma.07g069700 (6,316,520–6,321,126 bp) encodes glycerol-3-phosphate acyltransferase (GPAT),which is mainly involved in formation of the cell cuticle.Knockout of GPAT4 and GPAT8 caused structural changes to stomata in Arabidopsis thaliana with loss of cuticles in the guard cells.Following overexpression,the C16 and C18 contents in cutin monomers in leaves and stems increased by 80%,and transpiration of water was significantly decreased [68].

Block 2 contained six drought-resistance genes.Glyma.07g093000 (8,681,254–8,683,612 bp) was annotated as an Slocus lectin protein kinase family protein.In our study,Glyma.07g093000 had four non-synonymous mutations (T/C,C/T,C/G,and G/T) between two parents.Glyma.07g093200,in the 8,694,869–8,700,465 bp interval,encodes an armadillo repeat protein involved in the abscisic acid response.This protein interacts with transcription factor ABF2,which controls ABA-dependent gene expression via G-box-type ABA-responsive elements.The ABF subfamily of bZIP factors consists of four members,ABF1-ABF4 [69].Analyses of ABF overexpression and knockout lines showed that ABF2,ABF3,and ABF4 each play overlapping but distinct roles in ABA function and abiotic stress response [70].ABF2 plays an essential role in regulation of seedling growth and glucose response,and overexpression results in altered ABA/stress responses.ABF2 is unique in that,unlike ABF3 and ABF4,its overexpression phenotypes are dependent on the developmental stage[71].Additionally,there was a non-synonymous mutation (T/C)difference in Glyma.07g093200 between the parents of the RIL population.Annotation of Glyma.07g094000 shows that it is related to ATP-binding cassette A2;members of the Glyma.07g094000 ABCA gene subfamily showed differential expression after drought stress in wild soybean,indicating that the present gene was related to drought resistance [72].Glyma.07g094100 and Glyma.07g094200 encode protein kinase superfamily proteins,and there were two and one non-synonymous mutations,respectively,between the two parents.Glyma.07g094800 is a drought resistance-related gene in block 2,and encodes a protein kinase superfamily member with similar function to Glyma.07g069000.

The amino acid sequences of genes Glyma.07g093000,Glyma.07g094100,and Glyma.07g094200 were changed by nonsynonymous mutations,and the secondary and tertiary structures also changed accordingly.Regional association analysis showed that SNP CHR7-17619 was significantly associated with Glyma.07g094200,explaining 25.62% of the phenotypic variance.A non-synonymous mutation was also present within the resequenced subset.The proteins of the tertiary structure predicts that gene Glyma.07g094200 encodes brassinosteroid insensitive 1-associated receptor kinase 1 (BAK1).BAK1 functions by binding to transmembrane ligand receptors activated by transphosphorylation.BAK1 is a common co-receptor in plants where it plays important roles in growth and development and in biotic and abiotic stress.Several studies have demonstrated that due to its multiple functions BAK1 can be used as a growth regulator to increase crop yield and to improve adaptation to adverse conditions [73–76].A quadruple brassinosteroid (BR) receptor mutant (bri1,brl1,brl3,bak1) was likely to be the determining factor for promoting drought resistance without impairing yield in Arabidopsis thaliana[77,78].However,no study on drought resistance of this mutant has been performed in soybean.Phenotype verifications show the DSP value of the lines containing the ‘‘TT” allele in Glyma.07g094200 was significantly (P ＜0.05) larger than that of the lines containing the ‘‘GG” allele.Thus markers based on SNP CHR7-17619 in Glyma.07g094200 can be used in selection for drought resistance in soybeans.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement

Lijuan Qiu and Xusheng Qi:Conceptualization and Writing–review &editing.Zhangxiong Liu and Honglei Ren:Resources.Yue Li,Yongsheng Li and Zuowang Gou:Investigation.Xingzhen Zhao and Huilong Hong:Data curation.Yanjun Zhang,Zhangxiong Liu,and Xiongrong Wang:Writing– original draft.

Acknowledgments

This work was supported by National Key Research and Development Program of China (2016YFD0100201),Scientific Research Conditions Construction and Achievement Transformation Project of Gansu Academy of Agricultural Sciences (Modern Biological Breeding) (2019GAAS07),Science and Technology Major Project of Gansu Province(18ZD2NA008),and Crop Germplasm Resources Protection (2017NWB036-5).

Appendix A.Supplementary data

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