ZHANG Hua,WU Hai-yan,TlAN Rui,KONG You-bin,CHU Jia-hao,XlNG Xin-zhu,DU Hui,JlN Yuan,Ll Xi-huan,ZHANG Cai-ying

State Key Laboratory of North China Crop Improvement and Regulation/Hebei Agricultural University,Baoding 071001,P.R.China

Abstract Insufficient available phosphorus in soil has become an important limiting factor for the improvement of yield and quality in soybean.The mining of QTLs and candidate genes controlling soybean phosphorus utilization related traits is a necessary strategy to solve this problem.In this study,11 phosphorus utilization related traits of a natural population of 281 typical soybean germplasms and a recombinant inbred line (RIL) population of 270 lines were evaluated under different phosphorus conditions at two critical stages:the four-leaf stage as the seedling critical stage was designated as the T1 stage,and the six-leaf stage as the flowering critical stage was designated as the T2 stage.In total,200 single nucleotide polymorphism (SNP) loci associated with phosphorus utilization related traits were identified in the natural population,including 91 detected at the T1 stage,and 109 detected at the T2 stage.Among these SNP loci,one SNP cluster (s715611375,ss715611377,ss715611379 and ss715611380) on Gm12 was shown to be significantly associated with plant height under the low phosphorus condition at the T1 stage,and the elite haplotype showed significantly greater plant height than the others.Meanwhile,one pleiotropic SNP cluster (ss715606501,ss715606506 and ss715606543)on Gm10 was found to be significantly associated with the ratio of root/shoot,root and total dry weights under the low phosphorus condition at the T2 stage,and the elite haplotype also presented significantly higher values for related characteristics under the phosphorus starvation condition.Furthermore,four co-associated SNP loci (ss715597964,ss715607012,ss715622173 and ss715602331) were identified under the low phosphorus condition at both the T1 and T2 stages,and 12 QTLs were found to be consistent with these genetic loci in the RIL population.More importantly,14 candidate genes,including MYB transcription factor,purple acid phosphatase,sugar transporter and HSP20-like chaperones superfamily genes,etc.,showed differential expression levels after low phosphorus treatment,and three of them were further verified by qRT-PCR.Thus,these genetic loci and candidate genes could be applied in markerassisted selection or map-based gene cloning for the genetic improvement of soybean phosphorus utilization.

Keywords:soybean,phosphorus utilization,GWAS,linkage mapping,genetic loci,candidate genes

1.lntroduction

Soybean is an important crop worldwide,which provides high quality food protein and oil for human use (Phansaket al.2016;Leamyet al.2017;Liuet al.2021).However,in recent years,insufficient available phosphorus in the soil has become an important limiting factor for soybean yield and quality improvements in many countries (Bilyeuet al.2008;George and Richardson 2008;Konget al.2018;Liaoet al.2020).Therefore,the mining of genetic loci and candidate genes controlling soybean phosphorus utilization and the breeding of new soybean varieties with high phosphorus absorption and utilization efficiencies has become the most effective and indispensable strategy to solve this problem (Konget al.2018).

Some recent studies have focused on mining the genetic loci that control soybean phosphorus utilization related traits such as shoot fresh and dry weights,root fresh and dry weights,plant phosphorus content and utilization efficiency,etc.For example,Liet al.(2005) identified seven QTLs with phenotypic variation explanation (PVE) 8.8–25.6% controlling three phosphorus utilization related traits.Cuiet al.(2007) detected seven QTLs on five linkage groups (LGs) controlling phosphorus utilization related traits with PVE 4.8–17.0%.Genget al.(2007) identified 15 QTLs controlling five related traits with PVE 4.0–13.8%.Suet al.(2009) obtained three QTLs controlling plant height,root dry weight and leaf phosphorus on LG N and LG O with PVE 4.84–18.20%.Zhanget al.(2009) reported 34 QTLs linked with five related traits on nine linkage groups with PVE 6.6–19.3%.Lianget al.(2010) obtained 31 QTLs controlling 12 related traits on five linkage groups with PVE 9.1–31.1%.

Meanwhile,Liet al.(2016) detected 172 QTLs on 29 genomic regions conferring photosynthesis and phosphorus utilization related traits.Zhanget al.(2016)identified 20 QTLs on 17 chromosomes with PVE 5.14–12.59%,and eight of them were novel loci for phosphorus utilization.Caiet al.(2018) identified 14 major QTLs linked with root morphology related traits with PVE 4.23–13.98%,and the major QTL (qP10-2) was found to contain an acid phosphatase geneGmHAD1,andGmHAD1was ultimately demonstrated to have functions in improving phosphorus utilization efficiency in transgenic plants.Recently,Yanget al.(2019) fine-mapped the phosphorus efficiency locus 1(PE1) to a relatively narrow region (～200 kb) on Gm11,where 16 candidate genes were predicted.

In addition to the related QTLs detected from a biparental segregating population through the linkage mapping strategy above,Ninget al.(2016) identified 19 associated single nucleotide polymorphisms (SNPs)on nine soybean chromosomes (Gm01,Gm02,Gm03,Gm04,Gm08,Gm13,Gm15,Gm17 and Gm19) from a natural population (219 accessions) through a genomewide association study (GWAS) approach.The SNPs were associated with five phosphorus utilization related traits,such as plant height,shoot dry weight,plant phosphorus concentration,plant phosphorus acquisition and utilization efficiencies.Furthermore,Zhanget al.(2014) identified one major QTL (qPE8) with PVE 41%that was associated with soybean phosphorus utilization through a combination of the GWAS (192 accessions)and linkage mapping (152 lines) methods,then the functional geneGmACP1in QTLqPE8was shown to improve soybean hairy root phosphorus efficiency,and the gene functional markers were developed based on the polymorphisms ofGmACP1genomic sequences in different soybean germplasms (including the parents).

Although some studies have already identified several QTLs associated with phosphorus utilization in soybean,the genetic basis and molecular mechanism of phosphorus utilization still remain largely unknown.In addition,the genetic improvement of high phosphorus efficiency varieties has become increasingly necessary in soybean breeding programs.In view of these considerations,the objectives of this study are (1)to identify the associated QTLs controlling soybean phosphorus utilization related traits using two types of populationsviathe combination of the GWAS and linkage mapping strategies at two crucial phosphorus requirement stages;and (2) to identify and verify the candidate genes associated with phosphorus utilization related traits in the consistent QTL regions.

2.Materials and methods

2.1.Soybean materials

A natural population consisting of 281 germplasms(55 landraces and 226 cultivars) derived from different provinces,including Hebei,Shandong,Shanxi and Henan,etc.,was used for the GWAS in this study.In addition,a recombinant inbred line (RIL) population named ZN-RIL(high phosphorus utilization efficiency variety Zhonghuang 15×low phosphorus utilization efficiency germplasm Niumaohuang,F6:7,270 lines) was used for linkage mapping and verification of the genetic loci obtained from GWAS.

2.2.Methods

Soybean planting and experimental designThe soybean natural population was planted in circular plastic pots (diameter～11 cm) filled with vermiculite.The plants were grown in an artificial climate chamber with 12 h of supplemental light (～10 000 lx) per day,and the room temperature was set at (28±1)°C.Two different phosphorus treatments were conducted one week after sowing.The low phosphorus treatment used phytate (the major form of organic phosphorus in soil,which cannot be absorbed and utilized by plants directly) as the sole phosphorus supply (1.0 mmol L–1phosphorus) and the normal phosphorus treatment used KH2PO4as the sole phosphorus supply (1.0 mmol L–1phosphorus) according to previous studies (Wu 2016;Konget al.2018;Wuet al.2020).The experiment was conducted by a split-block design (phosphorus treatment as main plot and soybean germplasm as subplot) with four replications (two plants per replication).In addition,the ZN-RIL population was planted,managed,and treated nearly the same as the natural population above.

Phosphorus efficiency-related trait analysisAt 21 d(T1stage,four-leaf stage as the seedling critical stage)and 36 d (T2stage,six-leaf stage as the flowering critical stage) after planting,11 phosphorus utilization related traits (PH,plant height;RL,root length;R/S,ratio of root/shoot;SFW,shoot fresh weight;SDW,shoot dry weight;RFW,root fresh weight;RDW,root dry weight;TFW,total fresh weight;TDW,total dry weight;PC,phosphorus concentration;PUE,phosphorus utilization efficiency) of the natural population and the ZN-RIL population were analyzed using the methods described in Table 1.

Table 1 Measurement methods of phosphorus utilization related traits in soybean

Statistical analysis of population phenotypesThe descriptive statistics of the related traits were obtainedviaMicrosoft Excel 2007.The variance analysis of the split-block design was conducted using the general linear model in Microsoft Excel 2007,with the block factor treated as a random effect,the phosphorus treatment as a fixed effect,and genotype as a random effect (Gai 2000).The broad-sense heritability was estimatedviathe formulah2=Vg/VP.The correlation coefficients of related traits were calculated using IBM SPSS Statistic 21 Software,and the correlation hot-plot was obtained through the Pheatmap Software in R Statistical Program(https://www.r-project.org/) (Wanget al.2015).

Genotype and association analysis of the natural populationThe natural population was genotyped with the SoySNP6K array as described in previous studies(Akondet al.2013;Songet al.2013;Huoet al.2019),and STRUCTURE 2.3.4 Software was used to analyze the soybean population structure (Evannoet al.2005).Additionally,linkage disequilibrium (LD) decay was analyzed with TASSEL 5.0 Software (Bradburyet al.2007).The LD decay distances were chosen using the method of decreasingr2to half of the maximum value (Huoet al.2019).The association analysis of related traits and SNP markers was carried out with GAPIT Software in the R Statistical Program,and theP-value was set to 0.001 as the Type I error significance threshold (Hwanget al.2014;Tanget al.2016).The nearby associated SNPs with a physical distance of less than 270 kb were considered the same locus (LD decay distance in the present study)according to the description in Huoet al.(2019).

Genotype analysis and linkage mapping of the RlL populationThe ZN-RIL population was genotyped and the linkage mapping was analyzed using similar methods as our previous study (Liet al.2018).Briefly,ZN-RIL wasgenotyped with 1 015 SSR markers,and JoinMap 4.0 Software was used to construct the linkage map.The relatively wide genetic variations and significant differences had already been demonstrated in our previous study (Wu 2016;Wuet al.2020).Therefore,the QTL mapping was conductedviathe IciMapping 4.2 Software with the combination of the linkage map and the phenotypes of phosphorus utilization related traits according to Menget al.(2015).

Analysis of candidate genes in associated regionsTo further study the candidate genes in associated regions,an analysis of the expression levels of candidate genes was performed using the transcriptome data of Zhonghuang 15(a high phosphorus utilization efficiency variety) under low and normal phosphorus treatments (1.0 mmol L–1phytate phosphorus as the low treatment,and 1.0 mmol L–1KH2PO4as the normal treatment;with root samples harvested at 0,7,28,49 and 70 days after treatment).The qRT-PCR method was also used to further demonstrate the expression patterns of selected candidate genes as described in our previous study (Konget al.2018).

3.Results

3.1.Genetic variations and correlation coefficients of phosphorus utilization related traits in the soybean natural population

Significant differences in 10 of the 11 phosphorus utilization related characteristics (except for the root fresh weight) between the two phosphorus treatments(low and normal phosphorus) were found in the soybean natural population (Table 2),which indicated that the low phosphorus treatment influenced the growth and development of soybean plants.Specifically,some of the root related traits (such as the root dry weight and the ratio of root/shoot) were much higher under the low phosphorus condition than under the normal phosphorus condition,while the root fresh weights were equal,which indicated that the plants could stimulate many low-phosphorus molecular mechanisms in the root to adapt the phosphorus stressed environment.Meanwhile,the broad genetic variations of these related traits were found through the CV% values (11.67–31.64%),which indicated that many genetic loci and candidate genes existed in the natural population resulting in these wide variations of phenotype.The normal distributions were also indicated by the absolute kurtosis and skewness values (0.01–0.95),which implied that these related characteristics were attributed to quantitative traits and were controlled by multiple genes.Additionally,many significant correlations of these phosphorus utilization related traits were detected under the different phosphorus conditions at the two stages,such as shoot fresh weight,shoot dry weight,root fresh weight,root dry weight,total fresh weight,and total dry weight under the low phosphorus condition,etc.(Fig.1),which indicated that some pleiotropic genetic loci might exist which confer these significantly correlated traits.

Fig.1 Correlation coefficients of phosphorus utilization related traits at four-leaf (T1) and six-leaf (T2) stages.A and B,top right,T1 and T2 stages under low phosphorus,respectively;bottom left,T1 and T2 stages under normal phosphorus,respectively;the diagonal indicates the same trait between low and normal phosphorus at T1 and T2 stages.LP,low phosphorus treatment;NP,normal phosphorus treatment.PH,plant height;RL,root length;R/S,ratio of root/shoot;SFW,shoot fresh weight;SDW,shoot dry weight;RFW,root fresh weight;RDW,root dry weight;TFW,total fresh weight;TDW,total dry weight;PC,phosphorus concentration;PUE,phosphorus utilization efficiency.＊ and ＊＊,significant differences at 0.05 and 0.01 levels,respectively.

3.2.Genetic structure and optimal model analysis of the natural population

The soybean natural population was divided into two subgroups (SG1 and SG2) according to the population structure analysisviadelta K and Q value methods,and SG1 consisted of 220 germplasms,while SG2 consisted of 61 genotypes.Furthermore,three different models(MLM,GLM and FarmCPU) were tested to select the optimal model in this study,and the results showed that the MLM Model was optimal because it effectively reduced the statistical hypothesis errors in GWAS (Fig.2).

Fig.2 Optimal model analysis of soybean natural population.

3.3.SNP loci associated with phosphorus utilization at the T1 and T2 stages in the natural population

SNP loci at the T1 stageIn total,91 significant SNP loci were found to be associated with the 11 phosphorus utilization characteristics described above,65 of which were identified under the low phosphorus condition on 20 chromosomes,while the other 26 were identified under the normal phosphorus condition on 12 chromosomes(Table 3;Appendices A–D).Among these SNP loci,29 pleiotropic loci were identified on 15 chromosomes,of which,18 pleiotropic loci were detected under the low phosphorus condition (Table 3).

The first SNP locus (Locus1),which contained two SNP markers (ss715597964 and ss715597976) on Gm07,was associated with shoot fresh and dry weights,and root and total dry weights (Table 3).Locus2 on Gm08,which included two SNP markers (ss715601736 and ss715601748),was associated with soybean root fresh and dry weights,and total fresh and dry weights.Locus3 on Gm12,which included four SNP markers(ss715611375,ss715611377,ss715611379 and ss715611380),was found to be associated with soybean plant height,and Locus4 (ss715604183,ss715604186 and ss715604222) on Gm09 was associated with shoot,root and total dry weights.Meanwhile,Locus5(ss715608757) on Gm11 was found to be associated with shoot dry weight,root fresh weight,and total fresh and dry weights,while Locus6 (ss715614266) on Gm13 was associated with total dry weight,and shoot fresh and dry weights.Since the soybean plant height and weight were seriously influenced by soil phosphorus starvation,and the reductions of plant height and weight have been regarded as the typical symptoms for soybean growth under low phosphorus stress conditions,these genetic loci controlling soybean plant height and weight could be applied in soybean breeding for phosphorus improvement.

Table 2 Phenotypic variations of phosphorus utilization related traits under the different treatments at two stages in the soybean natural population

Table 3 Associated single nucleotide polymorphism (SNP) loci under the low phosphorus (LP) condition at the four-leaf (T1) stage in the natural population

Table 3 (Continued from preceding page)

In addition to the SNP loci above,another 59 SNP loci were found to be associated with phosphorus utilization related traits under the low phosphorus condition (Table 3;Appendix A).Meanwhile,26 SNP loci were found to be associated with 11 related traits under the normal phosphorus condition (Appendix B),and they were located on 12 chromosomes (Gm01,Gm02,Gm03,Gm05,Gm06,Gm07,Gm08,Gm09,Gm11,Gm14,Gm18 and Gm20).Among these loci,11 pleiotropic loci were found to be associated with nine related traits (except for the plant height and root length) on eight chromosomes.

SNP-loci at T2stageIn total,109 significant SNP loci were found to be associated with 11 phosphorus utilization characteristics,including 53 loci on 18 chromosomes identified under the low phosphorus condition,and 56 loci on 19 chromosomes identified under the normal phosphorus condition (Table 4;Appendices E–H).Among these SNP loci,30 pleiotropic loci were identified on 15 chromosomes,of which,17 pleiotropic loci were detected under the low phosphorus condition (Table 4).

Table 4 Associated single nucleotide polymorphism (SNP) loci under the low phosphorus (LP) condition at the six-leaf (T2) stage in the natural population

Table 4 (Continued from preceding page)

The first SNP locus (Locus1,containing three SNP markers,ss715606501,ss715606506 and ss715606543)on Gm10 was associated with the ratio of root/shoot,root dry weight and total dry weight (Table 4).Locus2 on Gm10,which included two SNP markers (ss715607012 and ss715607022),was associated with plant height,the ratio of root/shoot,root dry weight and total dry weight.Since the soybean plant could increase the root growth and development under the low phosphorus condition,there might be some candidate genes controlling soybean root growth in this region on Gm10.

Meanwhile,Locus3 on Gm04 (ss715588661),was found to be associated with the ratio of root/shoot,shoot fresh weight,and root and total dry weights.Locus4 on Gm05 (ss715590328) and Locus5 on Gm16(ss715623346) were found to be associated with the ratio of root/shoot,and root and total dry weights.More importantly,Locus6 on Gm06 (ss715594989) and Locus7 on Gm18 (ss715628776) were found to be associated with root phosphorus concentration and phosphorus utilization efficiency under the low phosphorus condition.

In addition to these SNP loci,another 46 SNP loci were found to be associated with phosphorus utilization related traits under the low phosphorus condition (Table 4;Appendix E),while there were 56 loci on 19 chromosomes(except for Gm02) associated with 11 related traits under the normal phosphorus condition (Appendix F).Among these loci,13 pleiotropic loci on 10 chromosomes were associated with 10 related traits (except for the plant height).

SNPs co-associated with phosphorus utilization related traits at the two stagesTo identify the genetic loci co-associated with phosphorus utilization related traits at the T1and T2stages,the significantly associated SNPs were re-analyzed and the four co-associated SNPs were found on Gm07,Gm08,Gm10 and Gm15 (Table 5).For these loci,the co-associated traits included the plant height,shoot fresh and dry weights,and root and total dry weights.

Table 5 Co-associated SNPs for phosphorus utilization related traits at two stages in soybean

Among these co-association SNPs,the pleiotropic genetic locus on Gm07 (ss715597964) was associated with root dry weight,total dry weight,and shoot fresh and dry weights at the T1stage,and also associated with root dry weight and total dry weight at the T2stage under the low phosphorus condition.Meanwhile,the SNPs ss715607012 on Gm10 and ss715622173 on Gm15 were associated with plant height under the low phosphorus condition at both stages.The SNP ss715602331 on Gm08 was associated with shoot fresh weight at the T1stage and with plant height and shoot dry weight under the low phosphorus condition at the T2stage.These associated SNPs at the T1and T2stages not only demonstrate the genetic loci controlling soybean phosphorus utilization,but also offer some important loci for soybean phosphorus improvements from vegetative growth to reproductive growth.

Elite genotypes in the significant associated SNP lociTo further validate the association results and select elite genotypes,the haplotypes of the associated loci were studied (Figs.3 and 4).The results showed that the SNP loci on Gm12 (ss715611375,ss715611377,ss715611379 and ss715 611380;physical position 11.01–11.16 Mb) generated two haplotypes (haplotype-I and haplotype-II),and the average plant height of haplotype-I(G–T–C–A,19 germplasms) was 32.13 cm,which was much greater than the height of haplotype-II (T–C–T–G,261 germplasms) under the low phosphorus condition.Meanwhile,the SNP loci on Gm10 (ss715606501,ss715606506 and ss715 606543;physical position 37.55–37.70 Mb) also had two haplotypes,and the G–A haplotype (haplotype-I,33 germplasms) showed much higher root dry weight (0.29 g),total dry weight(0.87 g) and ratio of root/shoot (0.36) than the A–G haplotype (haplotype-II,245 germplasms) under the low phosphorus condition.

Fig.3 Analysis of elite haplotype associated with plant height (PH) under low phosphorus (LP) condition at fourleaf (T1) stage.A,physical position of single nucleotide polymorphisms (SNPs) and linkage disequilibrium (LD) based on pairwise D′ (standardized disequilibrium coefficient)values between SNPs (values in boxes indicate the LD between SNP pairs,and the colors were applied via the standard Haploview color scheme;red indicates the LOD＞2 and D′=1;white indicates the LOD＜2 and D′＜1).11.01 and 11.40 Mb indicated the physical positions on chromosome 12.B,haplotype analysis of associated SNPs.C,phenotypic differences of PH among haplotypes.Bars mean SE.Different letters mean significant at 0.05 level.

3.4.Verification of the natural population associated SNP loci with phosphorus utilization in another soybean RlL population

To further verify the above associated SNP loci in the natural population,one RIL population named ZN-RIL was used to identify the linkage markers of the same phosphorus utilization related traits under the phosphorus treatments.The relative broad genetic variations and significant correlations of related traits demonstrated in our previous studies (Wu 2016;Wuet al.2020) also suggested that QTLs and candidate genes existed which confer the differences of the phenotypes.Based on this consideration,the linkage mapping analysis was conducted and the results showed that 12 QTLs were consistent with the genetic loci in the natural population(Appendix I).The SNP ss715622173 on Gm15 was associated with plant height in the GWAS population at the T1and T2stages under the low phosphorus condition(Table 5),and it was also located in the QTL region(qPH-E-1,Sat_380-BE347343,PVE 5.45%) for the plant height at the same stages under the low phosphorus condition in the ZN-RIL population.Meanwhile,the SNPs ss715622756 and ss715622776 on Gm15,which were associated with the ratio of root/shoot at the T1stage in the GWAS population,also overlapped with the QTLqPH-E-1in the ZN-RIL population under the low phosphorus condition.

Fig.4 Analysis of elite haplotype associated with root/shoot ratio (R/S),root dry weight (RDW) and total dry weight (TDW)under low phosphorus (LP) condition at six-leaf (T2) stage.A,physical position of single nucleotide polymorphisms (SNPs)and linkage disequilibrium (LD) based on pairwise D′ (standardized disequilibrium coefficient) values between SNPs (values in boxes indicate the LD between SNP pairs,and the colors were applied via the standard Haploview color scheme;red indicates the LOD＞2 and D′=1;shades of pink/red indicated the LOD＞2 and D′＜1;white indicated the LOD＜2 and D′＜1;blue indicates the LOD＜2 and D′=1).37.29 and 37.92 Mb indicated the physical positions on chromosome 10.B,haplotype analysis of associated SNPs (ss715606506 was not analyzed since its heterozygous genotypes in soybean varieties).C,phenotypic differences of R/S,RDW and TDW among haplotypes.Bars mean SE.Different letters mean significant at 0.05 level.

In addition to the co-located QTLs on Gm15,there were other co-located genetic loci on Gm08 between these two populations.Specifically,SNPs ss715601736 and ss715601748 (Locus2 at T1stage) were associated with root fresh and dry weights,and total fresh and dry weights under the low phosphorus condition (Table 3),and these two markers were located in the QTLqPH-A2-1region (Satt209-Sat_138) with PVE 22.39%of the ZN-RIL population.In addition,there were also co-located genetic loci on Gm04 between these two populations.The SNPs ss715589670 and ss715589694(Locus9 at T1stage) were found to associate with root phosphorus concentration,and shoot and total fresh weights under the low phosphorus condition (Table 3),and these markers were located in the QTLqRL-C1-1region (Satt646-Sat_416,PVE 5.44%) of the ZN-RIL population.

3.5.Candidate genes associated with phosphorus utilization related traits in soybean

To further analyze the candidate genes associated with phosphorus utilization related traits,the co-located SNP loci (the co-associated SNPs at the two stages in Table 5,the co-located SNP loci across the two populations in Appendix I,and the co-located SNP loci with strongly pleiotropic effects in Tables 3 and 4) were used to blast against Soybase data (www.soybase.org).The results showed that 14 candidate genes,including MYB transcription factor gene,purple acid phosphatase gene,sugar transporter and HSP20-like chaperones superfamily genes,etc.,adjacent to the associated SNP loci (～270 kb)were screened out according to the gene annotations and expression patterns (Table 6).The expression levels of these candidate genes were relatively higher after low phosphorus stress (Fig.5-A),especially for the candidate genesGlyma.10g176400,Glyma.12g114200andGlyma.04g101900,as determined through the analysis of root transcriptome data of high phosphorus utilization efficiency variety Zhonghuang 15 under the low and normal phosphorus treatments.

Table 6 Candidate genes with induced expression under the low phosphorus condition in soybean root

To demonstrate the expression patterns of these candidate genes,we selected three genes (Glyma.10g176400,Glyma.12g114200andGlyma.04g101900) for further analysis through qRT-PCR (Appendix J).The results showed that the expression levels of these genes under the low phosphorus condition were significantly higher than under the normal condition at 28 d (Fig.5-B),which further demonstrated the expression levels of these candidate genes in the root transcriptome sequencing analysis.Thus,these candidate genes and the associated SNP markers could be used for gene functional analysis and marker-assisted selection in the future.

Fig.5 Expression analysis of candidate genes association with soybean phosphorus utilization related traits in Zhonghuang 15(expression levels were indicated by the ratio of low phosphorus/normal phosphorus,and the root sample at 0 d was treated as the control).A,expression analysis of candidate genes by RNA-seq,data were shown by the log2Fold change.B,expressions of the candidate genes Glyma.04g101900,Glyma.10g176400 and Glyma.12g114200 by qRT-PCR.Bars mean SE (n=3).

4.Discussion

4.1.Co-located genetic loci controlling soybean phosphorus utilization related traits

Phosphorus is the major macro-element for soybean growth and development,and it participates in many important physiological and biochemical processes in soybean (Konget al.2018).However,most of the phosphorus in soil exists as the organic phosphate (such as phytate phosphate) which cannot be absorbed and utilized directly.Meanwhile,phosphorus utilization has been demonstrated as a quantitative trait that is controlled by minor-effect genes (Liet al.2005,2016;Cuiet al.2007;Suet al.2009).Thus,the discovery and application of the genetic loci and candidate genes controlling phosphorus utilization have become increasingly important in soybean breeding.

Association analysis and linkage mapping have become the two major strategies for QTL mining (Zhanget al.2014,2019;Juet al.2017;Liuet al.2019).Through the combination of these two strategies,some functional markers and genes have been identified and applied in plant breeding.Zhanget al.(2014) had found the soybean phosphorus efficiency controlling major QTLqPE8through the combination of linkage mapping and association analysis,and the functional geneGmACP1was identified.Meanwhile,the functional marker was developed based on the polymorphisms ofGmACP1in different germplasms with different phosphorus efficiencies,and the marker was finally used to identify soybean varieties in phosphorus utilization breeding.

In view of this success,the combination of association analysis and linkage mapping strategy was used in our study.When the association and segregation populations were treated through the same low phosphorus starvation treatment method at the same growth stages,the same phenotypes of these two populations were evaluated.Some co-located genetic loci were detected in these two populations (Appendix I),which indicated that some similar response mechanisms might exist in different soybean populations and varieties,and this also demonstrated the reliability of the genetic loci identified in our study.

Among these co-located genetic loci,the SNP ss715622173 on Gm15 associated with plant height in the association population was located in the plant height QTL(qPH-E-1) region of the mapping population (ZN-RIL) at the T1and T2stages under the low phosphorus condition in this study.Meanwhile,two SNPs (ss715622756 and ss715622776) on the same chromosome associated with the ratio of root/shoot at the T1stage in the association population also overlapped with the plant height QTLqPH-E-1under the low phosphorus condition.More importantly,Caiet al.(2018) also found an acid phosphatase activity QTLqRAPA-LPE1(41.66–42.00 Mb)and a root weight QTLqRN-LPE2(49.57–49.64 Mb) in this region,and some candidate genes were found to control phosphorus utilization related traits.

Besides the co-located QTL on Gm15,there was also another co-located QTL on Gm08 in this study.In the natural population,the SNPs ss715601736 and ss715601748 (40.64–40.70 Mb) were associated with root dry weight and total fresh weight under the low phosphorus condition at the T1stage.These two markers are located exactly in the QTLqPH-A2-1region with PVE above 10% of the ZN-RIL population.Meanwhile,Zhanget al.(2009) detected a ratio of root/shoot QTL(q2,15.07–39.91 Mb) in this region under low phosphorus conditions,and Caiet al.(2018) identified two QTLs(qTRL-LPE2,43.31–43.39 Mb;qP8-3,43.31–43.70 Mb)controlling root length in this region under phosphorus starvation conditions.

Meanwhile,there were other co-located genetic loci on Gm04 between these two populations.SNPs ss715589670 and ss715589694 were associated with root phosphorus concentration,and shoot and total fresh weights at the T1stage under the low phosphorus condition in the natural population,and they were also located in the QTLqRL-C1-1region in the ZN-RIL population.Meanwhile,Zhanget al.(2016) found a plant dry weight QTL (q4-1,11.44–18.52 Mb) in this region under low phosphorus conditions,and some reported QTLs conferring plant dry weight were also located in this region.The SPX family candidate geneSPX3(Glyma.04g147600) was identified near this region,and it was found to increase the root phosphorus contents of soybean transgenic hairy roots (Yaoet al.2014).More importantly,when thisSPX3was further analyzed through the RNA-seq results of Zhonghuang 15 in our study,relatively higher expression levels were found after phosphorus starvation in the high efficiency variety,which indicated thatSPX3might be one important functional gene for phosphorus utilization in this region on Gm04.

In addition,Liet al.(2016) detected a chlorophyll content related QTL with PVE 9.5% (q10,～36.60 Mb)on Gm10 under low phosphorus conditions;and Caiet al.(2018) identified root dry weight related QTL with PVE 6.19% (qRDW-LPE2,40.37–40.48 Mb) on Gm10 under low phosphorus conditions.In our study,some SNPs were found near these QTLs.For example,SNPs ss715606501,ss715606506 and ss715606543 on Gm10 were near QTLq10and were associated with the ratio of root/shoot,root dry weight and total dry weight in our study.Furthermore,SNP ss715607012 on the same chromosome near QTLqRDW-LPE2was associated with plant height at the T1and T2stages under the low phosphorus condition in our study.When the haplotypes of ss715606501,ss715606506 and ss715606543 were analyzed,the favored haplotype showed much higher values for the ratio of root/shoot,and root and total dry weights than the others,which indicated that there were some candidate genes in this region which control the plant weights under phosphorus starvation environments.

Other than these co-located QTLs above,there were also some novel and newly discovered genetic loci in our study.For example,SNP ss715597964 on Gm07 was associated with root and total dry weights,and shoot fresh and dry weights at the T1stage,and with root dry weight and total dry weight at the T2stage under the low phosphorus condition.In addition,SNPs ss715611375,ss715611377,ss715611379,and ss715611380 on Gm12 were associated with soybean plant height at the T1stage under the low phosphorus condition,and the favored haplotype showed much greater plant height than the unfavored haplotype.These novel genetic loci could also be applied in soybean phosphorus utilization molecular breeding in the future.

4.2.Comparison of genetic loci under different low phosphorus stress conditions in the same population

In China,different soil phosphorus contents exist in different soybean ecological regions.Usually,a soil available phosphorus concentration above 20 mg kg–1is regarded as an adequate supply for plants,and a level of ＜10 mg kg–1is considered a low phosphorus concentration in soil (Zhanget al.2010;Hwanget al.2014).Thus,some references have selected ＜10 mg kg–1as the low phosphorus stress concentration to identify the phenotypes of soybean germplasms.For example,Zhanget al.(2010) evaluated a RIL population using 6.51 mg kg–1as the low phosphorus treatment,and 20 mg kg–1as the normal treatment.Zhanget al.(2014)analyzed RIL and natural populations using ＜0.5 mg kg–1as low phosphorus,and ＞20 mg kg–1as the normal treatment.Ninget al.(2016) evaluated 219 germplasms using 8.51 mg kg–1as the low phosphorus treatment,and＞20 mg kg–1as the normal treatment.Yanget al.(2019)studied various germplasms using 11.24 mg kg–1as low phosphorus,and 25.33 mg kg–1as the normal treatment.Lianget al.(2010) evaluated a group of germplasms using 11.2 mg kg–1as the low phosphorus treatment.

To identify the genetic loci under different low phosphorus stress conditions,and to demonstrate the genetic loci detected in this study,another low phosphorus concentration experiment (EXP–0.5,0.5 mmol L–1phytate)using the same soybean natural population was designed by our research team (Zhanget al.2020),and the results of the 0.5 mmol L–1phytate (EXP–0.5) and 1.0 mmol L–1phytate (EXP–1.0,this study) conditions were compared.This comparison revealed many co-associated genetic loci (Appendix K),which indicated that similar response mechanisms existed under the different low phosphorus starvation conditions in soybean,and also demonstrated the reliability of the associated genetic loci in this study.

Among the co-associated genetic loci under these two different phosphorus starvation conditions (0.5 and 1.0 mmol L–1phytate),the SNPs ss715611375,ss715611377 and ss715611379 (i.e.,Locus3,T1stage)were associated with plant height under both of these low phosphorus stress conditions (Appendix K).Meanwhile,the markers ss715604183 and ss715604186 (i.e.,Locus4,T1stage) were found to be associated with root dry weight under these two low phosphorus stress conditions,and ss715612154 (i.e.,Locus12,T1stage) was associated with both the shoot fresh and dry weights.Thus,these genetic loci identified under the different low phosphorus stress concentrations could be used in different ecological regions with various soil phosphorus contents,and could also be used for marker-assisted selection or gene cloning in soybean phosphorus utilization genetic improvement efforts in the future.

4.3.Candidate genes controlling phosphorus utilization related traits in soybean

Plants have evolved many mechanisms to adapt to the low phosphorus stress environment in soil,including the molecular mechanism through the high efficiency controlling genes.Some related genes have been reported to control plant phosphorus absorption and utilization under phosphorus starvation conditions,such as the MYB transcription factor geneAtPHR1(Nilssonet al.2007) andOsPHR2(Wu and Wang 2008),the purple acid phosphatase geneAtPAP15(Wanget al.2009) andGmPAP14(Konget al.2018),the phosphorus transporter geneAtPT1(Muchhalet al.1996),AtPT2(Mukatiraet al.2001) andOsPht1;1(Zhanget al.2021),the aluminum activated malate transporter geneGmALMT5(Penget al.2018) andAtALMT1(Balzergueet al.2017;Mora-Macíaset al.2017).To discover new functional genes related to phosphorus utilization,the candidate genes near the associated SNPs were investigated,and some related genes that control phosphorus utilization were screened out.

Among these candidate genes,there were two MYB transcription factor genes namedGlyma.04g101900(MYB93) andGlyma.08g293300(MYB149) on Gm04 and Gm08,respectively,which were found to be induced by low phosphorus stress in the high efficiency variety in our study.Glyma.04g101900is a homologous gene toAtMYB93(Gibbset al.2014) that was induced after 28 d of phosphorus starvation in our study;while08g293300is a homologous gene toIAR4(Quintet al.2009) that was induced after 7 and 28 d of phosphorus starvation in our study.InArabidopsis,AtMYB93was reported to be specifically expressed in lateral root primordia cells,and showed functions in promoting the growth and development of plant lateral roots (Gibbset al.2014).IAR4inArabidopsiswas also reported to have functions in regulating the development of lateral roots and root hairs (Quintet al.2009).Since the roots play an irreplaceable function in allowing plants to adapt to the low phosphorus stress environment,and soybean usually flowers at～30 d after sowing,we deduced that these two MYB transcription factor genes on Gm04 and Gm08 might participate in the absorption and utilization of phosphorus at the soybean flowering stages.

In addition to the MYB transcription factor genes,a purple acid phosphatase geneGlyma.08g291600(GmPAP17) on Gm08 was identified in our study.Purple acid phosphatase is an important acid phosphatase that is usually secreted into the plant rhizosphere to resolve the organic phosphorus in soil for plant growth (Konget al.2018).Some purple acid phosphatase genes have been identified and a few of them are used in plant breeding.One report indicated that there were 35 purple acid phosphatase genes in the soybean genome,andGmPAP17was reported to have higher expression levels in soybean roots after phosphorus starvation (Liet al.2012).In our study,SNP ss715601736 on Gm08 was associated with phosphorus utilization related traits in the natural population,and was also located in the QTLqPH-A2-1region with PVE 22.39% of the ZN-RIL population.Furthermore,the candidate geneGmPAP17was very near to SNP ss715601736 (physical distance～93 kb),thus the purple acid phosphataseGmPAP17was considered as one of the important candidate genes for this genetic locus on Gm08 for adapting to the low phosphorus stress environment.

5.Conclusion

In total,200 SNP loci associated with phosphorus utilization related traits were identified,with 91 detected at the T1stage,and 109 detected at the T2stage.Among these loci,one SNP cluster (ss715611375,ss715611377,ss715611379 and ss715611380) located on Gm12 was found to be significantly associated with plant height under the low phosphorus condition at the T1stage,and the elite haplotype showed a much greater plant height than the others.Meanwhile,the pleiotropic SNP cluster(ss715606501,ss715606506 and ss715606543) on Gm10 was found to be significantly associated with the ratio of root/shoot,root dry weight and total dry weight under the low phosphorus condition at the T2stage,and the elite haplotype had much higher values for related traits than the others.Furthermore,four co-associated SNP loci were identified at the two stages under the low phosphorus condition,and 14 candidate genes which were close to the associated pleiotropic loci,such asGlyma.10g176400,Glyma.12g114200andGlyma.04g101900,were selected through the transcriptomic data.Three of them were further demonstrated by qRT-PCR and showed strongly induced expression levels under the low phosphorus condition.These SNP loci and candidate genes could be applied in marker-assisted selection or map-based gene cloning for genetic improvements in soybean phosphorus utilization.

Acknowledgements

This research was funded by the Project of Hebei Province Science and Technology Support Program,China (17927670H and 16227516D-1).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm