Association analysis for detecting significant single nucleotide polymorphisms for phosphorus-deficiency tolerance at the seedling stage in soybean [Glycine max (L) Merr].

Tolerance to low-phosphorus soil is a desirable trait in soybean cultivars. Previous quantitative trait locus (QTL) studies for phosphorus-deficiency tolerance were mainly derived from bi-parental segregating populations and few reports from natural population. The objective of this study was to detect QTLs that regulate phosphorus-deficiency tolerance in soybean using association mapping approach. Phosphorus-deficiency tolerance was evaluated according to five traits (plant shoot height, shoot dry weight, phosphorus concentration, phosphorus acquisition efficiency and use efficiency) comprising a conditional phenotype at the seedling stage. Association mapping of the conditional phenotype detected 19 SNPs including 13 SNPs that were significantly associated with the five traits across two years. A novel cluster of SNPs, including three SNPs that consistently showed significant effects over two years, that associated with more than one trait was detected on chromosome 3. All favorable alleles, which were determined based on the mean of conditional phenotypic values of each trait over the two years, could be pyramided into one cultivar through parental cross combination. The best three cross combinations were predicted with the aim of simultaneously improving phosphorus acquisition efficiency and use efficiency. These results will provide a thorough understanding of the genetic basis of phosphorus deficiency tolerance in soybean.

Introduction

Phosphorus (P) is an essential macronutrient required for plant growth and development. However, P availability is limited in most soils worldwide due to P fixation by Fe and Al oxides and the formation of organic P complexes, especially in tropical and subtropical areas (Vance , Wang ). Soybean, one of the most important food and oil crops in the world, is generally considered to have a high P requirement for optimal seed yield and quality. Thus, a lack of available phosphate in soil is one of the primary stresses that reduce soybean productivity. Application of P fertilizers can alleviate this problem, but excess P application leads to water source contamination and reduction of limited P resources (Cordell ). Therefore, breeding soybean cultivars with improved ability to access soil-bound P and/or efficiently utilize endogenous P is a sustainable and economical approach to solve this problem.

Phosphorus deficiency affects soybean development and physiology throughout its ontogeny, e.g., by reducing the biomass and P concentration of seedlings, decreasing the leaf area (Fredeen ), and changing the balance of diurnal starch accumulation and utilization (Qiu ). In general terms, early developmental stages are more vulnerable to phosphorus deficiency than mature stages. Phosphorus deficiency during early developmental stages results in a large reduction in the number of soybean nodes and kernels, and recovery cannot occur even if sufficient P is supplied at later stages (Wang 1963, Wu ). A similar phenomenon has been reported for wheat and barley (Brenchley 1929, Gericke 1924): P deficiency during the first 4–6 weeks of growth decreases the development of the tiller and the formation of the head in barley and wheat. These results indicate that rapid and vigorous plant seedling growth under P deficient conditions may be an important prerequisite for optimal soybean yield. Therefore, variation in certain traits at the seedling stage might be of importance in breeding for phosphorus-deficiency tolerance. Further, determination of the genetic control of phosphorus-deficiency tolerance during early vegetative development is necessary for developing P-efficient varieties.

However, the genetic basis of phosphorus-deficiency tolerance is complex and controlled by multiple genes. Quantitative trait locus (QTL) analysis is often used to dissect phosphorus-deficiency tolerance-related complex traits. In recent years, several QTL mapping studies have been conducted to better understand P efficiency in rice (Wissuwa ), wheat (Su ), maize (Zhu ), Brassica napus (Yang ), and common bean (Beebe ). Nevertheless, to our knowledge, the genetic analysis of phosphorus efficiency in soybean is only in the early stages compared with that in other crops. Currently, some QTLs underlying phosphorus-deficiency tolerance at the seedling stage in soybean have been detected using bi-parental segregating populations. Three QTLs for fresh shoot weight at the seedling stage were mapped on the F linkage group by Li . Zhang identified 34 additive QTLs for five parameters related to phosphorus-deficiency tolerance at the seedling stage by using a RIL population derived from Nannong 94–156 and Bogao. Liang detected 31 QTLs for twelve phosphorus-deficiency tolerance-related traits, which were distributed on five linkage groups. However, classical marker development methodologies for trait dissection can be inefficient and expensive because bi-parental design suffers from two main limitations, namely, limited allelic diversity and poor mapping solution, due to limited recombination events during population creation.

Fortunately, association analysis studies can overcome these limitations. Association analysis is another analytical approach for QTL detection that exploits linkage disequilibrium (LD) between markers and closely linked QTLs in irregular population structures. This powerful technique enables researchers to exploit natural diversity and locate candidate genes in the genome (Korte and Farlow 2013). With the development of assay-based high-throughput marker systems, association mapping for QTL discovery has gained popularity in a range of plants (Atwell , Huang , Kump , Sukumaran ). Association mapping has been successfully used for soybean to identify favorable alleles for various traits, including yield-related traits (Sonah ), quality-related traits (Hwang ), and resistance traits (Mamidi ). Recently, a genome-wide association study was conducted to accelerate molecular breeding for the improvement of phosphorus efficiency related traits during the reproductive stage in a natural population. GmACP1, which is an important gene related to P efficiency in soybean, was identified by combining genome-wide association analysis, linkage analysis, and other analyses (Zhang ). GmACP1 encodes an acid phosphatase which was assessed by expression in Escherichia coli (BL21). When GmACP1 was over-expressed in soybean roots, the hairy roots of transgenic plants secreted more acid phosphatase to catalyze the hydrolysis of ρ-nitrophenyl phosphate and release inorganic phosphate for soybean absorption (Zhang ). These results demonstrate that association mapping can also be used to identify favorable alleles for phosphorus efficiency in soybean.

In this study, 219 soybean accessions were used to conduct an association analysis by combining information on 1142 single nucleotide polymorphism (SNP) markers. Five traits of soybean seedlings, including plant shoot height (PH), shoot dry weight (SDW), phosphorus concentration (PC), phosphorus acquisition efficiency (PAE) and use efficiency (PUE), were phenotyped in 2012 and 2013 under two different P conditions. The aims of this study were to determine the genetics of phosphorus-deficiency tolerance at the soybean seedling stage and to detect SNPs associated with phosphorus-deficiency tolerance with the aim of improving soybean stress resistance using marker-assisted selection (MAS).

Materials and Methods

Plant materials and culture

A set of 219 soybean landraces with diverse geographic origins and variant phenotypes were phenotyped to perform association mapping. For all accessions, the germplasm name, origin and the probability of being classified into each subpopulation are listed in Supplemental Table 1. This experimental population, which included landraces, cultivars and breeding lines, originated from 26 different provinces and six ecological regions of soybean cultivation in China (latitude 53 to 24ºN and longitude 134 to 97ºE). The experimental population was used to identify quantitative trait nucleotides (QTN) associated with yield (Hao ), SC7 resistance (Yan ), and protein content (Zhang ). These accessions were selected because they display diverse genetic variation in their phosphorus efficiency.

Two independent pot trials at the soybean seedling stage were carried out in 2012 and 2013. The soil used in the pot trials was collected from the Jiangpu Experimental Station of Nanjing Agricultural University (32º12′ N, 118º37′ E), Nanjing, China. According to previous studies, soil is considered phosphorus deficient if the concentration of soil available P is less than 10 mg kg−1, whereas an ideal soil for plant growth has an available P concentration higher than 20 mg kg−1 (Bray and Kurtz 1945, Colwell and Esdaile 1968). The total and available P in the soil used in this study were 942 and 8.51 mg kg−1, respectively; therefore, this soil was considered phosphorus-deficient. The experiment was performed using a completely randomized design with a split-plot restriction according to Furlani . The main plots were defined by their P treatment (plus P and minus P), while the subplots were the 219 genotypes. The genotypes were sown in a pot (two plants per pot) containing 4 kg of dry soil. There were three replicates with six plants in each trial. The soil for the plus P treatment received 60 mg kg−1 of N (as urea) and 20 mg kg−1 of K (as potassium chloride) and 20 mg kg−1 of P (as potassium dihydrogen phosphate). The soil for the minus P treatment contained the same amounts of N and K but no P. Chemical fertilizer was applied as a basal dressing before sowing using a broadcast application. Tap water was used for irrigation when the growth conditions were consistent with drought.

Evaluation of plant performance and phenotypic analysis

Plant shoot height (PH) was measured before harvest as the length of the cotyledon node to the top of the growing point. The aerial parts of the plant from the cotyledon node were harvested after 35 days of growth. Shoot dry weight (SDW) was determined by incubating the sample at 105ºC for 60 min and then oven-drying at 65ºC for 72 hours. The phosphorus concentration (PC) was determined using nitric acid digestion and ICP emission spectroscopy as described by Li . The phosphorus acquisition efficiency (PAE) was calculated as the total amount of phosphorus in each shoot, and shoot P use efficiency (PUE) was calculated as the shoot dry weight produced per unit of phosphorus (g dry weight/mg phosphorus) (Gourley ). We analyzed the phenotypes of five traits measured under the plus P and minus P conditions, as well as the modified phenotype calculated from the phenotypes under the two P conditions. The modified phenotype was obtained as the residual in regressing the phenotypic value under the minus P condition on that under the plus P condition. This modified phenotype is referred to as the conditional phenotype and is considered to be a suitable index for phosphorus-deficiency tolerance. The conditional phenotype was calculated using QGAStation 2.0 software (Zhu 1995).

Statistical analysis

Statistical analysis for all traits was conducted using SAS 8.1 (SAS Institute, Cary, NC, USA). Analysis of variance (ANOVA) for all phenotypic data was performed using PROC GLM. Correlation analyses were conducted using PROC CORR for each trait based on the mean values measured over two years. Trait measurements were averaged based on the replications prior to association analysis. The broad-sense heritability estimate (h2) was calculated for each trait across different P environments according to the approach described by Shi .

Association analysis

In the present experimental population, 192 of 219 soybean landraces were genotyped for 1536 SNPs described by Hao . To reduce the risk of false-positive marker-trait associations, SNPs showing residual allelic heterogeneity in the accession were considered missing data as well as SNPs with minor allele frequencies of 10% or less. Thus, a total of 1142 SNPs were utilized for the association analysis study. Association analysis based on conditional phenotypes is termed conditional association analysis (Cui , Zhang and Wang 2015, Zhu 1995), whereas analysis using the original phenotype is referred to as unconditional association analysis. TASSEL 4.0 (Bradbury ) software was used to conduct association mapping. Association mapping was performed using the mixed linear model (Q+K), which greatly reduces false positive rates (Yu ). The population structure (Q), inferred using the software program STRUCTURE 2.2 (Pritchard ), was previously described by Hao . The kinship matrix (K) was estimated using TASSEL 4.0 software. Only those markers that met the criterion of −log (P) > 2.0 for their P values in both years were considered to be significant SNPs. A Manhattan plot of the P values for all traits examined in the association analyses was constructed using Origin 8.0 (Origin. ver.8.0 for Windows 7, OriginLab., 2008, USA).

Prediction of novel parental combinations

To identify novel elite combinations of parents from the experimental population, we evaluated the novel parental combinations according to the predicted phenotypic values of the hypothetical offspring from the crosses between parental cultivars. First, we selected the elite accessions, which were with relatively higher mean values of phosphorus efficiency-related traits in both years, and contained more than one favorable allele of stable SNPs determined by t-test, for crossing to generate hypothetical progenies. To predict the phenotypic values of the hypothetical offspring, we assumed that the markers that were significantly associated with traits co-segregated with the loci controlling the traits, and that the alleles at different loci were in accordance with the law of independent assortment. Then, the phenotypic values of CPAE and CPUE were predicted for the hypothetical crossed offspring as the sums of the QTL effects that were estimated based on their QTL genotypes and the mean values of the experimental population over two years for the conditional phenotypes. The hypothetical offspring of a parental combination is termed the elite offspring, which was with highest phenotypic values of CPAE and CPUE in all the hypothetical progenies. Finally, the best three cross combinations were selected for breeding in future according to relatively higher predicted phenotypic values of CPAE and CPUE of those elite offspring.

BLAST analysis

A number of candidate genes and pathways that may be involved in the evolution of phosphorus-deficiency tolerance have been obtained in Arabidopsis thaliana (Fang , Plaxton , Zhang ). These genes were used as a query in a BLASTP analysis against soybean proteins in Phytozome v10.2 (http://www.phytozome.net). Genes were selected based on the top 40 hits with an E-value cut-off of 10−20 (Mamidi ).

Results

Significant variation among soybean accessions with respect to phosphorus-efficiency traits

The mean values, standard deviations, ranges, and broad-sense heritability of the five traits related to seedling growth under the two P conditions are shown in Table 1. The phenotypic values of these five traits showed large variation among the association panel in the two years of experiments, with coefficients of variation (CVs = SD/Mean × 100%) ranging from 17.43% for PC in the minus P condition in 2013 to 57.14% for PUE in the minus P condition in 2012. The CVs of most traits exceeded 20%. Phosphorus was significantly correlated with soybean growth in all the accessions in the field. Under the minus P condition, the mean phenotypic values for PH, SDW, PC and PAE significantly decreased, but PUE significantly increased in both years. Analysis of variance showed that genotype, P level, and P level × genotype had significant effects on the five traits in the experimental population, suggesting that the treatments and genetic background were very important in explaining the overall phenotypic variation. The distribution of the trait values among all individual lines fell into the expected normal continuous distributions (Fig. 1), indicating polygenic control.

Table 1

Means, range, variation, ANOVA, and broad-sense heritability results of five phosphorus efficiency traits from the association population

Trait	Year	Treatment	Mean ± SD	CV%	Range	ANOVA			h2
						Pa	Gb	P × G
PH (cm)	2012	Minus-P	31.72 ± 6.74	21.25	14.00–48.67	**	**	**	0.63
		Plus-P	38.20 ± 9.39	24.58	19.67–65.33
	2013	Minus-P	37.17 ± 8.14	21.9	18.50–58.00	*	**	**
		Plus-P	43.97 ± 11.95	27.18	21.33–76.5
SDW (g)	2012	Minus-P	2.74 ± 0.83	30.29	1.16–5.79	**	**	**	0.47
		Plus-P	3.55 ± 1.19	33.52	1.39–8.00
	2013	Minus-P	2.84 ± 0.99	34.86	1.15–6.62	**	**	**
		Plus-P	3.95 ± 1.53	38.73	1.36–10.21
PC (μg/g)	2012	Minus-P	3968 ± 861.21	21.7	1720–6055	*	**	**	0.53
		Plus-P	6165 ± 1398	22.68	2547–10644
	2013	Minus-P	3417 ± 595.65	17.43	1512–4971	**	**	**
		Plus-P	6476 ± 1245	19.22	4107–12159
PAE	2012	Minus-P	10.82 ± 3.91	36.14	2.95–24.41	**	**	**	0.39
		Plus-P	21.63 ± 8.14	37.63	6.28–50.30
	2013	Minus-P	9.67 ± 3.63	37.54	3.44–22.70	*	**	**
		Plus-P	25.48 ± 11.24	44.11	9.47–79.20
PUE	2012	Minus-P	0.14 ± 0.08	57.14	0.08–0.35	**	**	**	0.51
		Plus-P	0.09 ± 0.05	55.56	0.04–0.27
	2013	Minus-P	0.15 ± 0.06	40	0.10–0.33	**	**	**
		Plus-P	0.08 ± 0.03	37.5	0.04–0.12

Phosphorus level;

Genotype.

PH: plant shoot height; SDW: shoot dry weight; PC: phosphorus concentration; PAE: phosphorus acquisition efficiency; PUE: phosphorus use efficiency.

Fig. 1

Frequency distribution of the five phosphorus efficiency traits under minus P and plus P conditions based on the 2-year data in the association population. PH: plant shoot height; SDW: shoot dry weight; PC: phosphorus concentration; PAE: phosphorus acquisition efficiency; PUE: phosphorus use efficiency.

Relationships among the five traits

Simple correlations among the five traits were calculated based on the means of two years of continuous data from 219 accessions (Table 2). With the exception of SDW and PUE, which were not significantly correlated, both positive and negative correlations were detected between the other traits under both the minus and plus P conditions. PAE was strongly positively correlated with SDW and moderately positively correlated with PH and PC at both P levels, which implies that PAE and other traits may share related genetic determinants or mechanisms. However, PAE was moderately negatively correlated with PUE, which is consistent with observations from previous studies (Cui , Geng ). The highest significant negative correlation (r = −0.94) was found between PUE and PC for the minus P condition. PUE and PAE were significantly negatively correlated under both P levels. PC showed a low but significant correlation with SDW under both P conditions, which suggests that shoot biomass may have exerted a “dilution effect” on the shoot phosphorus concentration in the experimental population.

Table 2

Sample correlation coefficients between different phosphorus traits in the association population under minus P (below diagonal) and plus P (above diagonal) conditions

	PH	SDW	PC	PAE	PUE
PH		0.72***	−0.20**	0.56***	0.16*
SDW	0.65***		−0.15*	0.84***	0.23
PC	−0.23***	−0.11*		0.38***	−0.92***
PAE	0.43***	0.80***	0.48***		−0.38***
PUE	0.24***	0.18	−0.94***	−0.47***

indicate significance at the 0.05, 0.01, and 0.001 levels, respectively.

PH: plant shoot height; SDW: shoot dry weight; PC: phosphorus concentration; PAE: phosphorus acquisition efficiency; PUE: phosphorus use efficiency.

The conditional and unconditional phenotypes of the five traits from the experimental population were used for association analysis. We only selected those markers that were significant in both years and met the criteria of having a significant association threshold (P < 0.01). The SNPs associated with the conditional phenotypes were referred to as conditional SNPs, whereas the SNPs associated with the original phenotypes were termed unconditional SNPs. For unconditional association analysis, a total of 39 significant unconditional SNPs-trait associations involving 23 SNPs were mapped to nine chromosomes for the five traits under the plus and minus P conditions (Supplemental Table 2). Of these 39 unconditional SNP-trait associations, 19 were identified based on the phenotype data from the plus P condition, and 20 were detected based on the minus P condition values. Chromosome 3 had the largest number of SNPs (a total of 11 SNPs), whereas only one SNP was observed on chromosome 13.

Nineteen conditional SNP-trait associations involving 13 SNPs met the significance criteria for conditional association analysis (Fig. 2, Table 3). These SNPs were distributed on nine chromosomes, including Chr1–Chr4, Chr8, Chr13, Chr15, Chr17, and Chr19. Interestingly, almost half of the SNPs were mapped to chromosome 3 (9 SNPs).

Fig. 2

Manhattan plot of the P values in conditional association analysis for the five phosphorus efficiency traits. The values of −log (P) were determined by association analysis of data over the two years and the average data. The blue dots represent the observed data in 2012, and the dark yellow dots represent the data in 2013, the dark cyan dots represent the observed combined data. The gray line represents the genome-wide line (−log (P) > 2.0). CPH: conditional phenotype of plant shoot height; CSDW: conditional phenotype of shoot dry weight; CPC: conditional phenotype of phosphorus concentration; CPAE: conditional phenotype of phosphorus acquisition efficiency; CPUE: conditional phenotype of phosphorus use efficiency.

Table 3

Significant SNPs detected using the conditional trait values for the five phosphorus efficiency traits in two years

Trait	Marker	CHRa	Position (bp)	2012		2013
				−log (P)	R2 (%)	−log (P)	R2 (%)
CPH	BARC-028709-05992	3	40654334	2.01	3.61	2.31	4.28
CPH	BARC-029599-06229	15	5897144	2.12	5.30	2.08	5.18
CSDW	BARC-030665-06919	2	10233524	2.11	3.87	2.95	5.86
CSDW	BARC-028205-05792	3	40462124	3.02	7.67	3.28	8.31
CSDW	BARC-028709-05992	3	40654334	2.48	4.70	2.50	4.73
CSDW	BARC-028671-05985	4	48121001	2.17	5.46	2.14	5.35
CPC	BARC-024477-04900	1	2699256	3.70	10.12	3.58	9.01
CPC	BARC-030965-06980	3	2785339	2.78	6.90	3.49	8.69
CPC	BARC-018809-03022	3	41074807	5.23	13.37	2.38	5.85
CPC	BARC-045047-08867	8	8109865	2.67	5.04	2.11	3.76
CPC	BARC-013257-00462	13	32207135	2.48	4.60	2.42	4.46
CPAE	BARC-028205-05792	3	40462124	2.21	5.50	2.55	6.32
CPAE	BARC-028709-05992	3	40654334	2.09	4.57	2.40	4.54
CPAE	BARC-018809-03022	3	41074807	5.00	12.88	2.03	5.01
CPAE	BARC-031827-07220	17	1785893	2.80	7.03	3.30	8.29
CPAE	BARC-013481-00498	19	42334572	2.97	7.47	2.96	7.38
CPAE	BARC-014509-01564	19	48637997	2.22	5.52	3.06	7.65
CPUE	BARC-024477-04900	1	2699256	3.70	10.12	2.78	9.16
CPUE	BARC-030965-06980	3	2785339	4.59	12.68	3.17	10.48

Chromosome.

CPH: conditional phenotype of plant shoot height; CSDW: conditional phenotype of shoot dry weight; CPC: conditional phenotype of phosphorus concentration; CPAE: conditional phenotype of phosphorus acquisition efficiency; CPUE: conditional phenotype of phosphorus use efficiency.

BARC-028709-05992 on chromosome 3 and BARC-029599-06229 on chromosome 15, which were detected using the data from 2012 and 2013, were significantly associated with the conditional phenotype of plant shoot height (CPH). Four significant SNPs (BARC-030665-06919, BARC-028205-05792, BARC-028709-05992 and BARC-028671-05985) for the conditional phenotype of shoot dry weight (CSDW) accounted for 7.38–10.48% of the phenotypic variance. BARC-028205-05792 and BARC-028709-05992 were both located on chromosome 3. BARC-030665-06919 and BARC-028671-05985 were distributed on chromosome 2 and chromosome 4, respectively (Table 3). A total of five significant SNPs with a phenotypic variance explained (PVE) of 3.76–9.01% were found to control the conditional phenotype of phosphorus concentration (CPC) over the two years. Six SNPs that were significantly associated with the conditional phenotype of phosphorus acquisition efficiency (CPAE), and had PVE values ranging from 4.28% to 8.31%, were repeatedly found in two years. BARC-024477-04900 and BARC-030965-06980, which were located on chromosomes 1 and 3, respectively, were two coincident SNPs associated with the conditional phenotype of phosphorus use efficiency (CPUE) across different years.

It is worth noting that a novel SNP cluster with three significantly associated SNPs (BARC-028205-05792, BARC-028709-05992, and BARC-018809-03022) was observed using the conditional phenotype values. The three SNPs that showed significant conditional SNP-trait associations were located within the same genomic region of chromosome 3 and were significantly co-associated with more than one trait. Among the three significant SNPs, BARC-028205-05792 (SNP position: 40462434 bp) was associated with CSDW and CPAE and was located near SNP BARC-018809-03022 (SNP position: 41075023 bp), which was co-associated with CPC and CPAE. In addition, the distance between BARC-028205-05792 and BARC-028709-05992 (SNP position: 40654716 bp), which was associated with controlling CPH, CSDW, and CPAE, was only 0.19 Mb. The coincidence of the SNPs in one chromosomal region across separate years not only indicated the reliability of the SNPs reported above but also implied that these markers may exhibit high efficiency in MAS.

Mining favorable alleles

We determined the favorable alleles of the SNP associated markers for phosphorus efficiency based on the means of conditional phenotypic values of each trait over the two years. A t-test revealed that out of 19 significantly associated SNPs, the phenotypic effects of the alleles of 14 SNPs were significantly different for the corresponding associated traits (Fig. 3, Table 4). Interestingly, the 14 distinct associated markers were included in the components of the novel SNP cluster. For BARC-028709-05992, which was co-associated with CPH, CSDW and CPAE, a comparison of 80 A-type and 111 G-type soybean accessions demonstrated that allele A of BARC-028709-05992 was the favorable allele for the three traits. Two significant SNPs, BARC-024477-04900 and BARC-030965-06980, which were associated with CPC and CPUE, showed different alleles as favorable for these two traits. It was found that allele G of the two SNPs was the favorable allele for CPUE, whereas allele A was the favorable allele for CPC. Each of the remaining seven SNPs have their own separate favorable alleles for each trait. For example, the significant SNP BARC-028205-05792, associated with CSDW across the two years, was identified in 79 A-type and 112 C-type soybean accessions. The CSDW of A-type soybean accessions was significantly increased by 0.07 relative to that of the C-type soybean accessions (t-test, P < 0.001). For the significant SNP BARC-014509-01564, which was associated with CPAE, 155 C-type soybean accessions showed significantly higher CPAE than the 36 A-type accessions, with an average increase of 0.08 per C allele (t-test, P < 0.001). The results of the t-test confirmed that these loci might be responsible for the genetic variation in phosphorus efficiency between the genotypes.

Fig. 3

Association of marker allele polymorphisms with the five phosphorus efficiency traits. The sandybrown plot represents the association of marker allele polymorphisms with the conditional phenotype of plant shoot height (CPH), the blue plot represents the association of marker allele polymorphisms with the conditional phenotype of phosphorus use efficiency (CPUE), the yellow plot represents the association of marker allele polymorphisms with the conditional phenotype of phosphorus concentration (CPC), the purple plot represents the association of marker allele polymorphisms with the conditional phenotype of shoot dry weight (CSDW), and the green plot represents the association of marker allele polymorphisms with the conditional phenotype of phosphorus acquisition efficiency (CPAE).

Table 4

T-test of five phosphorus efficiency traits among accessions grouped by the alleles of association SNPs across two years

Trait	Marker	CHRa	Position (bp)	Alleles	Number	Mean	T-test
CPH	BARC-028709-05992	3	40654334	A	80	0.87	P < 0.01
				G	111	0.82
CPH	BARC-029599-06229	15	5897144	A	73	0.88	P < 0.001
				G	118	0.82
CSDW	BARC-030665-06919	2	10233524	A	147	0.77	P < 0.001
				T	44	0.7
CSDW	BARC-028205-05792	3	40462124	A	79	0.8	P < 0.001
				C	112	0.73
CSDW	BARC-028709-05992	3	40654334	A	80	0.79	P < 0.01
				G	111	0.73
CPC	BARC-024477-04900	1	2699256	A	137	0.63	P < 0.001
				G	54	0.54
CPC	BARC-030965-06980	3	2785339	A	170	0.61	P < 0.001
				G	21	0.49
CPC	BARC-018809-03022	3	41074807	A	46	0.58	P < 0.05
				C	145	0.61
CPC	BARC-013257-00462	13	32207135	A	148	0.62	P < 0.001
				G	43	0.55
CPAE	BARC-028709-05992	3	40654334	A	80	0.47	P < 0.05
				G	111	0.44
CPAE	BARC-013481-00498	19	42334572	A	77	0.5	P < 0.001
				T	114	0.42
CPAE	BARC-014509-01564	19	48637997	A	36	0.39	P < 0.001
				C	155	0.47
CPUE	BARC-024477-04900	1	2699256	A	137	1.63	P < 0.01
				G	54	1.81
CPUE	BARC-030965-06980	3	2785339	A	170	1.63	P < 0.001
				G	21	2.05

Chromosome.

CPH: conditional phenotype of plant shoot height; CSDW: conditional phenotype of shoot dry weight; CPC: conditional phenotype of phosphorus concentration; CPAE: conditional phenotype of phosphorus acquisition efficiency; CPUE: conditional phenotype of phosphorus use efficiency.

The typical carrier variety containing favorable alleles with the best phenotypic value might be used to design parental combinations for crop breeding (Niu ). PAE and PUE are two components of phosphorus efficiency. To simultaneously improve P uptake and P utilization efficiency in soybean, we might aggregate all the favorable alleles for the traits into one cultivar to the greatest degree possible in the following way. First, the favorable alleles of stable SNPs were selected using a t-test, and the elite accessions containing more than one favorable allele and with a relatively high mean value of phosphorus efficiency-related traits in both years were selected from the experimental population. Then the phenotypic values of CPAE and CPUE were predicted for the hypothetical crossed offspring of elite accessions as the sums of the QTL effects that were estimated based on their QTL genotypes and the mean values of the experimental population over two years for the conditional phenotypes. Finally, we selected the elite parental combinations for crossing to generate high phosphorus-efficiency progenies according to the elite offspring with higher phenotypic values of CPAE and CPUE. To simultaneously improve P uptake and P utilization efficiency five favorable alleles could be pyramided by the parental combinations of Y8D6002 × Y8D6051, Y8D6002 × Y8D6096, and Y8D6096 × Y8D6179 (Tables 5, 6).

Table 5

Significant SNPs with favorable alleles among the selected accessions

Locus	Selected accessions
	Y8D6002	Y8D6051	Y8D6096	Y8D6179
BARC-028709-05992	A	/	/	A
BARC-024477-04900	/	A	A	/
BARC-030965-06980	A	A	A	A
BARC-013481-00498	A	/	A	/
BARC-014509-01564	C	C	C	C

Table 6

Top three superior crosses for phosphorus efficiency improvement

P1					P2					RILs
Taxa	CPAE	CPUE	Pa-CPAE	P-CPUE	Taxa	CPAE	CPUE	P-CPAE	P-CPUE	Mb-CPAE	M-CPUE
Y8D6002	0.62	1.52	0.49	1.96	Y8D6051	0.48	1.57	0.44	2.13	0.49	2.13
Y8D6002	0.62	1.52	0.49	1.96	Y8D6096	0.38	2.00	0.47	2.13	0.49	2.13
Y8D6096	0.38	2.00	0.47	2.13	Y8D6179	0.49	1.86	0.46	1.96	0.49	2.13

P1 and P2: the two parental lines of the cross; RILs: the recombinant inbred lines from the parental cross; CPAE: conditional phenotype of phosphorus acquisition efficiency; CPUE: conditional phenotype of phosphorus use efficiency;

the phenotypic values of CPAE and CPUE were predicted as the sums of the QTL effects that were estimated based on their QTL genotypes and the mean values of the experimental population over two years for the conditional phenotypes;

the maximum predicted phenotypic value.

Discussion

Effect of phosphorus deficiency on seedling growth

In the present study, five morphological and physiological traits at the seedling stage were used as indicators to identify the genomic regions that control phosphorus efficiency. A strong influence of the traits was observed under the minus P condition (Table 1). PAE was more severely reduced by phosphorus deficiency than PH and SDW under the minus P condition. Among other components of phosphorus efficiency, PUE was increased, which is in agreement with a previous study (Su ). Therefore, this result implied that PAE may be the limiting factor for shoot development at the seedling stage under phosphorus deficient conditions. This was verified by correlation analysis (Table 2). Higher positive correlations were observed between PAE and agronomic traits than those observed with PUE, PH and SDW. These results indicate the importance and complexity of understanding the genetic basis of phosphorus-deficiency tolerance in soybean breeding.

Conditional association analysis for phosphorus-deficiency tolerance

Phosphorus-deficiency tolerance in plants is a complex trait represented by different phenotypic traits in different growth periods. The aim of the QTL analysis of conventional phosphorus-deficiency tolerance is to analyze the expression of phenotypic traits at different P levels (Li , Liang , Zhang ). However, it is extremely difficult to discover the genetic relationship between these phenotypic traits and phosphorus-deficiency tolerance and to identify the loci concerned with phosphorus-deficiency tolerance because of the difference in trait expression under different P levels. Most unconditional QTLs for traits of soybean seedling under phosphorus-deficiency stress did not contribute directly to phosphorus-deficiency tolerance but only to the development of those traits. Conditional genetic analysis can effectively eliminate genetic inter-relationships among individual QTLs at different P levels and reveal QTLs that cannot be identified under unconditional analysis. This is because this method can exclude the influence of trait expression itself and isolate the variations that were only caused by different phosphorus treatments (Cui , Zhang and Wang 2015, Zhu 1995). In the present study, association analysis based on conditional phenotypes was defined as conditional association analysis, whereas that based on the original phenotypes was defined as unconditional association analysis. Less than one-third of the significantly associated conditional SNPs (five SNPs) were identified by unconditional analysis, after comparing of the results of conditional and unconditional analysis (Fig. 2, Table 3, Supplemental Table 2). The inconsistency in SNPs detected using conditional and unconditional analyses is in agreement with earlier reports (Cui , Zhang and Wang 2015, Zhang ). These findings imply that some significant SNPs that are associated specifically with phosphorus-deficiency tolerance could be identified only by eliminating the influence of the trait expressed under plus P conditions. In addition, these findings also imply that trait expression itself may suppress the expression of most conditional SNPs associated with phosphorus-deficiency tolerance. Hence, conditional association analysis is more viable and effective for phosphorus-deficiency-tolerance related traits.

Quantitative trait nucleotides (QTNs) for phosphorus efficiency and important SNP clusters

Previous reports have descried the determination of the genetic basis of phosphorus efficiency in soybean using bi-parental segregating populations and/or association panels (Li , Liang , Zhang , 2014b). It is very difficult to compare the results of different studies due to the use of different genetic maps, especially when marker density is very low. Therefore, we determined whether the conditional SNPs detected in our study were located at or near the same locations as those identified in previous reports by comparing the chromosomal positions of these QTLs (http://www.soybase.org/). Zhang performed QTL mapping and detected a total of 34 significant QTLs for five phosphorus efficiency-related traits at the soybean seedling stage in a RIL population. In our study, only BARC-029599-06229 for CPH and BARC-013481-00498 for CPAE were located near or in the region of the confidence interval for the two QTLs detected by Zhang . This result confirmed that a certain causal gene or genes for P homeostasis might exist in this genomic region. Li identified seven QTLs for three traits: fresh shoot weight, leaf phosphorus content, and root phosphorus content. Liang mapped three clusters of QTLs associated with root-related traits and P efficiency traits under a low P condition and one cluster under a high P condition. However, no SNPs were detected at or near these locations in the present study. The different mapping populations and experimental environments that were used in the studies mentioned above may lead to some inconsistency in the QTLs mapped in different studies. In a recent study, a number of QTLs/QTNs for phosphorus efficiency-related traits during the reproductive stage in a soybean RIL population and natural population were mapped (Zhang ). However, only a small number of our QTNs (BARC-013481-00498 and BARC-014509-01564) were near the QTLs/QTNs detected by Zhang because different soybean growth stages were used in the two studies. This difference is not surprising given that plant growth and development during different growth stages is controlled by different genes.

According to genomic sequencing, most identified genes are distributed in the euchromatic regions, whereas a few genes are found in the heterochromatic regions (Schmutz ). Hence, it is common that QTLs for different traits are clustered in particular genetic regions, and this phenomenon has been reported in different species (Hu , Obara ). In our study, three SNPs that were significantly associated with more than one trait were located in close proximity to each other on chromosome 3. An interesting SNP, BARC-028709-05992, which was detected for CPH, CSDW, CPAE, PlusPH (plant shoot height under the plus P condition), PlusSDW (shoot dry weight under the plus P condition), and MinusSDW (shoot dry weight under the minus P condition) in our study, was coincident with a QTL for phosphorus efficiency in soybean (Hang ). This QTL region has also been mapped in several different studies for 100-seed weight, seed yield, seed size, seed phosphorus concentration and chlorophyll fluorescence parameters (Hao , 2012b, Hu , Ning ). Another significant SNP, BARC-028205-05792, which is co-associated with PlusSDW, MinusPAE (phosphorus acquisition efficiency under the minus P condition), CSDW and CPAE, is also associated with chlorophyll and chlorophyll fluorescence parameters (Hao ). Moreover, GmPT2 (Glyma.03g162800), a Pht1 family gene, is located near BARC-028205-05792 (Supplemental Table 3). The expression of GmPT2 can complement the yeast P-uptake mutant (MB192) under P-limiting conditions (Qin ). The region of BARC-028709-05992 was located 168124 bp away from Glyma.03g194200, which is orthologous to the purple acid phosphatase gene (PAPs) (Supplemental Table 3). PAPs play a vital role in the adaptive strategies of leguminous plants to phosphorus deficiency (Liang , Zhang ). These findings suggest that the genetic regions of the SNP cluster mentioned above are important for increasing genotype selection intensity by MAS and should be given more attention in pyramiding breeding.

Predicted novel parental combination for molecular design breeding

Peleman and van der Voort (2003) first introduced the concept of “Breeding by Design” based on the extensive study of the technology and potential of marker-assisted selection. Breeders can design their crosses such that superior alleles are combined in one individual using a “menu” listing the superior alleles of breeding materials. PAE and PUE are two components of phosphorus efficiency. However, it is difficult to improve PAE and PUE simultaneously because of their negative correlation. Pyramiding more positive favorable alleles for PAE and PUE into one cultivar, as far as possible, may solve this problem. In our study, three parental combinations (Y8D6002 × Y8D6051, Y8D6002 × Y8D6096, and Y8D6179 × Y8D6096) using four selected accessions with high a P value in the t-test and more than one favorable allele could be used in soybean phosphorus-deficiency tolerance breeding programs (Tables 5, 6). For the parental cross combination of Y8D6179 × Y8D6096, the maximum predicted CPAE phenotypic value of the progeny (0.49) was higher than the predicted CPAE values of its two parents (0.47 and 0.46) because the favorable alleles that controlled CPAE were pyramided into the best progeny. Therefore, the best RIL of Y8D6096 × Y8D6179 with maximum CPAE and CPUE values may be selected as a more phosphorus efficient cultivar by breeders. The results of predicted parental cross combinations can be used for future molecular design breeding. According to previous reports (Wan 2006, Wang ), we know that the predicted results can be confirmed by further crop breeding. Wang selected three chromosome segment substitution lines (CSSL5, 16, and 19) as parental cross combinations from 65 CSSLs on the basis of the effects of the detected loci for rice grain length and width. Wang reported that the predicted length and width of the best offspring line were 6.21 mm and 3.07 mm, respectively. Recently, a japonica cultivar carrying a target genotype of indica genome fragments (>8.5 mm length × 3.2 mm width), derived from the designed target genotypes for large grain, was bred using molecular design breeding (Niu , Wang ). Hence, the predicted combinations in the present study will be valuable for guiding breeders in future soybean phosphorus-deficiency tolerance breeding programs.

Function of putative phosphorus efficiency-related genes

QTNs should be sufficiently robust for further fine mapping or for use in a phosphorus deficiency stress-breeding scheme. It is important to identify the key genes responsible for these phenotypic traits. A previous study determined that soybean and Arabidopsis diverged from a common ancestor approximately 90 million years ago (Grant ). The presence of homologs of partial genomic sequences and the conserved nature of certain residues of most genes suggest that many genes are likely to be functionally conserved across soybean and Arabidopsis (Grant ). These factors make comparative genome mapping/homology cloning possible between the two species. Thus far, hundreds of genes involved in P acquisition, transport, translocation and homeostasis have been reported in Arabidopsis (Fang , Plaxton , Zhang ). Therefore, we predicted specific candidate genes involved in phosphate metabolism by performing a BLASTP analysis and synteny comparisons between soybean and Arabidopsis according to Mamidi . Briefly, gene prediction in the target genomic region was performed based on soybean genome information and bioinformatics. Genes that were not only related to phosphorus efficiency (annotated as such or orthologous to genes of Arabidopsis) but also located in the decay distance of LD were considered candidate genes. The decay distance of LD is approximately 500 kb according to Hao . In total, 17 genes involved in phosphate homeostasis in Arabidopsis were mapped to the soybean LD regions of QTNs, including transporter genes, PAP family members and others (Supplemental Table 3). The approach of associating QTNs with functional genes in Arabidopsis by synteny comparisons will not only provide more information for dissecting the genetics of the mechanisms of phosphorus-deficiency tolerance in soybean but will also facilitate the utilization of phosphorus-deficiency stress or responsive genes in soybean breeding.