Introduction Phosphorus (P) is essential for all living cells and organisms. P occurs in complex DNA and RNA structures that contain and translate genetic information and control all living processes in plants and animals. Animals need to obtain an adequate supply of P from their food [1]. As for plants, P is one of the most essential mineral nutrients required for growth and development [2]. However, most unmanured soils do not contain sufficient readily available P to meet the high demands of crops, particularly during certain periods of the growing cycle. Therefore, fertilizers containing P must be supplied. In recent years, the annual global consumption of phosphates was more than 50 million tons [3]. However, excess P application is problematic because the excess dissolution of phosphates contaminates water sources [4], [5]. Furthermore, most of the annually fertilized phosphates are fixed in the soil in organic forms by adsorption, sedimentation and transformation. Consequently, 50–80% of the total P found in soil exists as organic phosphates, which are unavailable to plants in the absence of mineralization [6]. Therefore, agricultural soil is a large ‘potential P pool’ that must be developed and used. Furthermore, the global reserve of rock phosphate is a non-renewable resource. A recent study conducted by the International Fertilizer Development Center (IFDC) concluded that at the current extraction rates, the global commercial phosphate reserves will be depleted in approximately 300–400 years [7]. Therefore, the development of plants that can efficiently utilize endogenous and added P is a sustainable, economic approach for plant production. In response to persistent P deficiency, plants have developed many adaptive mechanisms to enhance the availability and uptake of P, including modifications in root architecture, increased activities of internal and extracellular acid phosphatases, greater exudation of small molecular organic acids, and symbiosis with mycorrhizal fungi [2]. Hundreds of genes associated with the plant P metabolic pathway have been identified [2], but few of these have been applied in crop breeding programs, perhaps because most of the identified genes were cloned by reverse genetics methods from a few model plants [6], [8]–[10], while few were identified by forward genetics methods [11]. For example, tremendous efforts have been made to dissect the genetic basis of soybean P efficiency by assessing factors such as biomass [12], root architecture [13], P concentration, acid phosphatase activity [14] and flower/pod abscission rates [15]. These studies, which were based on genetic analyses in segregating populations, identified several QTLs that likely control P efficiency. However, no QTLs associated with soybean P efficiency have been cloned. The reasons for these limited results may include the following: (a) the QTLs associated with P efficiency in soybean are generally localized to large chromosomal regions (10–20 centimorgans), (b) the effects of the QTLs are minor and often represent only a small fraction of the phenotypically relevant variation, and (c) QTLs are complex, and their evaluation requires information regarding germplasm diversity. Therefore, many challenging questions remain unanswered. For example, how can we identify the candidate P efficiency genes in these QTLs? Is P efficiency modulated by genetic variation of the candidate genes? What are the most favorable alleles and haplotypes useful for the breeding of high P efficiency? These questions reveal the highly challenging nature of the genetic dissection of P efficiency in soybean and other plants. Although QTL linkage mapping provides useful information on genetic loci, it is typically difficult to isolate candidate genes based on a single QTL mapping experiment. Furthermore, the genes identified by this method are restricted to those that segregate in the considered cross [16]. Genome-wide associations (GWAs) overcome this limitation and have recently been successfully applied in studies of Arabidopsis [17], rice [18], [19], maize [20] and other plants [21], [22]. This method can provide increased power for the localization of QTLs because of the higher recombination rates between markers and QTL alleles in randomly mating populations [23]–[25]. However, GWAs also have limitations because they could generate false positives as a result of population structure. Although population structure can be controlled by statistical methods [13], [17], [18], the complementary use of family-based linkage analyses in controlled crosses is also an option [26], [27]. The complementarity of GWAs and classical linkage analyses has been well demonstrated by studies of Arabidopsis flowering time [16] and rice Al tolerance [28]. In addition, because P efficiency is a typical quantitative trait controlled by multiple genes [29] and the molecular mechanism underlying soybean P efficiency is poorly understood, a candidate-gene association analysis could be an effective means to functionally characterize target genes. This strategy to study complex quantitative trait genes has been successfully applied for many plant traits such as maize flowering time [23], carotenoid content [30], architecture of floral branch systems [31] and rice seed shattering-loss [32]. Soybean is a highly important crop. Low P availability is the most significant soybean production constraint and is more problematic than other nutrient deficiencies, toxicities or diseases [33]. Because soybean is a model legume plant, the characterization of P efficiency related genes and the mechanisms responding to low-P stress in soybean would eventually facilitate P efficiency studies in other legumes and plants. We aimed to clone the soybean candidate gene for P efficiency, elucidate its effects, determine its favorable haplotypes and develop valuable functional markers. Therefore, we performed a series of experiments including linkage mapping, genome-wide association mapping, candidate-gene association mapping, gene expression and plant transformation. These studies revealed that a soybean acid phosphatase-encoding gene, GmACP1, located within the major QTL qPE8, is associated with P efficiency, and its genetic variation modulates P efficiency related traits in soybeans. Results Significant variation in P efficiency among soybean germplasms To determine the genetic variation of P efficiency in soybean plants, four P-efficiency-related traits were determined using 152 soybean recombination inbred lines (RILs) (Table 1) and 192 soybean accessions (Table 2). These traits included plant height (PHt), acid phosphatase activity (APA), leaf P concentration (PC) and 100-seed weight (100-SW) under low-P (−P, soil available P 20 mg kg−1) conditions; the relative values of the traits −P/+P are denoted as RPH, RAPA, RPC and 100-RSW, respectively. As shown in Table 1, the phenotypic RIL values ranged from 0.12–0.62 mg g−1 for PC, 0.13–0.58 mg g−1 for the pod P concentration (PPC) and 0.70–3.01 µmol ρ-NP min−1 mg protein−1 for APA in the low-P condition. The transgressive segregation of P-efficiency-related traits was obvious, and the phenotypic variation was significantly affected by the genotypes and treatments. In addition, analysis of variance (ANOVA) indicated that the phenotypic APA variation between the two parents (Bogao and Nannong 94–156) was significant (P = 0.01) in different P conditions (Table 1). The mean APA values for the individual accessions in the natural population ranged from 1.34–2.12 µmol ρ-NP min−1 mg protein−1, and the mean PC values ranged from 0.07–1.84 mg g−1. Among the lines of diverse soybean accessions, the PC levels reached 1.84 mg g−1; however, one soybean accession had a PC of only 0.07 mg g−1 (Table 2). Overall, the soybean plants clearly exhibited considerable natural variation in the traits related to P efficiency and displayed very high genetic diversity. 10.1371/journal.pgen.1004061.t001 Table 1 Descriptive statistical results for traits related to phosphorus (P) efficiency in soybean recombinant inbred lines (RILs) and their parents in experiments conducted in 2006 and 2007. Population Parents Trait Treatment Year Mean±SD Range Skew h2 (%)a Rb Gc Td Ye Bogao Nannong 94–156 PC −P 2006 0.30±0.08 0.12–0.53 0.94 57.1 ns s s ns 0.31 0.41 2007 0.31±0.09 0.12–0.62 0.12 0.31 0.45 +P 2006 0.40±0.09 0.20–0.77 0.69 41.3 0.41 0.45 2007 0.40±0.08 0.20–0.72 0.46 0.43 0.57 PPC −P 2006 0.33±0.07 0.15–0.58 0.54 57.6 ns s s ns 0.38 0.41 2007 0.30±0.08 0.13–0.56 0.33 0.39 0.46 +P 2006 0.41±0.08 0.26–0.66 0.72 50.6 0.49 0.55 2007 0.38±0.08 0.22–0.65 0.57 0.51 0.54 APA −P 2006 1.41±0.32 0.70–3.01 0.81 31.6 ns s s s 1.23 1.52 2007 1.45±0.40 0.75–2.80 0.98 1.31 1.42 +P 2006 1.30±0.29 0.53–2.08 −0.02 38.9 1.01 1.22 2007 1.02±0.37 0.45–2.27 0.84 1.12 1.13 PHt −P 2006 51.62±14.00 30.00–88.00 0.25 73.4 ns s s s 60.38 33.5 2007 57.1±14.08 25.83–97.33 −0.24 61.33 35.2 +P 2006 56.79±14.24 23.75–83.75 0.24 84.5 65.75 37.25 2007 64.62±13.24 24.33–89.05 −0.26 69.56 40.2 100-SW −P 2006 17.1±2.55 7.00–23.00 0.02 45.5 ns s s ns 13.36 14.51 2007 15.99±2.65 8.72–22.45 −0.02 12.34 13.5 +P 2006 18.46±2.31 12.98–27.77 0.74 62.5 17.28 19.12 2007 17.67±2.64 11.76–26.60 0.53 16.45 17.56 s: significant difference at P = 0.01; ns: not significant; PC: P concentration in the plant (mg g−1); PPC: P concentration in the pod at the seed filling stage (mg g−1); APA: acid phosphatase activity (µmol min−1 mg protein−1); PHt: plant height (cm); 100-SW: 100-seed weight (g). a heritability; b replication; c genotype; d treatment; e year. 10.1371/journal.pgen.1004061.t002 Table 2 Descriptive statistical results for traits related to phosphorus (P) efficiency in 192 soybean accessions in experiments conducted in 2008 and 2009. Trait Treatment Year/Site Mean±SD Range Skew Ra Gb Tc Yd PC −P 2008HN 0.28±0.08 0.14–0.52 0.14 ns s s s 2008NJ 0.28±0.10 0.14–1.00 1.5 2009NJ 0.50±0.31 0.07–1.53 0.97 +P 2008HN 0.35±0.07 0.16–0.67 0.37 2008NJ 0.40±0.12 0.13–1.00 0.74 2009NJ 0.82±0.45 0.15–1.84 0.37 APA −P 2008HN 2.12±0.66 0.12–6.47 0.99 ns s s ns 2008NJ 1.80±0.69 0.49–7.25 0.86 2009NJ 1.96±0.60 0.55–6.15 1.2 +P 2008HN 1.61±0.62 0.10–2.62 −0.35 2008NJ 1.34±0.46 0.39–3.80 0.44 2009NJ 1.54±0.38 0.18–3.77 0.73 s: significant difference at P = 0.01; ns: not significant; PC: P concentration in the plant (mg g−1); APA: acid phosphatase activity (µmol min−1 mg protein−1); HN and NJ: field experiments in Henan and Nanjing, China, respectively. a replication; b genotype; c treatment; d year. P efficiency related QTL linkage mapping in the RIL population Several P-efficiency-related QTLs have been identified in our previous works [14], [15], most of which focused on P efficiency at the soybean seedling stage. However, 90% of soybean P absorption occurs during the reproductive stage [33], [34]. Therefore, studying the P-efficiency-related traits during the reproductive stage is expected to be better suited for identifying the QTL associated with P efficiency. In this study, we initially aimed to identify the QTLs related to P efficiency during the soybean reproductive stage using 152 RILs in various environments. We used the relative values RPH, RPC, RPPC, RAPA and 100-RSW as indices to assess P efficiency in soybeans. Three primary QTL were found to be associated with P efficiency which were located on three soybean chromosomes (8, 14 and 18) (Table 3). Of the three loci, a stable QTL (for RAPA, RPPC and RPC) termed qPE8 (QTL for P efficiency related traits on chromosome 8) was identified; this QTL encompasses a 6.3-Mb region and explained up to 41% of the phenotypic variation (Table 3, Figure S1). This locus was consistent with previously identified QTL associated with P efficiency in soybean [14], [15]. 10.1371/journal.pgen.1004061.t003 Table 3 QTL analysis for traits related to phosphorus (P) efficiency in experiments conducted in 2006 and 2007 using 152 soybean recombinant inbred lines (RILs). Trait Year Chr. Marker interval Confidence Interval LOD ADD. R2 (%) RPC 2006 Chr.8 Satt089-Sat-310 155.8–161.7 9.8 0.17 41.05 Chr.13 Sat_262-Satt030 30.0–53.3 2.8 −0.08 7.72 2007 Chr.8 Satt089-Sat-310 156.0–163.4 13.4 0.17 37.74 RPPC 2006 Chr.8 Satt089-Sat-310 154.9–162.1 11.9 0.11 31.07 2007 Chr.5 Satt276-Sat_368 27.8–35.1 2.7 −0.05 5.67 Chr.8 Satt089-Sat-310 150.7–154.3 8.5 0.11 22.41 RAPA 2006 Chr.8 Satt089-Sat-310 150.9–160.2 6.8 0.2 18.05 2007 Chr.8 Satt089-Sat-310 141.0–160.2 3.5 0.14 19.02 RPHt 2006 Chr.14 Sat_424-Satt063 37.5–40.5 9.1 −0.03 21.28 Chr.1 Satt436-Sat_343 79.7–91.0 3.1 −0.02 8.2 2007 Chr.14 Sat_424-Satt063 37.6–40.6 10.2 −0.04 23.2 Chr.1 Satt436-Sat_346 88.0–92.7 3.2 −0.03 6.75 Chr.19 Sat_191-Satt313 201.1–214.3 3.2 0.02 7.05 100-RSW 2006 Chr.18 Sat_185-Satt012 63.8–65.8 13.5 0.08 26.06 Chr.18 Sat_352-Satt138 106.60–112.2 4.6 0.04 7.54 2007 Chr.18 Sat_185-Satt012 63.6–67.4 4.9 0.13 14.18 ADD: additive effect; R2 : the contribution ratio of QTL effect; RPC: relative P concentration in the plant at two P conditions; RPPC: relative P concentration in the pod at the filling stage at two P conditions; RAPA: relative acid phosphatase activity at two P conditions; RPHt: relative plant height at two P conditions; 100-RSW: relative value of 100-seed weight at two P conditions. GWAs for P efficiency related traits in a natural population To overcome the limitations of linkage analysis, we conducted GWAs to identify the loci associated with P efficiency using 1,536 single nucleotide polymorphisms (SNPs) [35] and the relative phenotypic values at several sites and over several years as P efficiency phenotypes generated in 192 soybean accessions. GWAs were conducted across all 192 genotypes using SNPs with minor allele frequencies (MAF)>0.05. In addition, the GWAs were separately associated with the P efficiency related phenotype, APA and PC across sites and years (Figure 1). A mixed linear model (MLM), which controlled for the complex population structure and pedigree relationships, was used in each analysis to correct for the confounding effects of the subpopulation structure and relatedness between individuals. The relative performance of the MLM in all traits was evaluated to control false positives in the studied population, as shown in the quantile-quantile (QQ) plots (Figure S2). The smaller deviation of the observed P value from the expected value indicated that the MLM was suitable for the control of type I errors. Seventy-four significant SNPs associated with P efficiency were identified and organized by the GWAs into six major clusters on chromosomes 8, 9, 11, 13, 18 and 19 across various environments (Figure 1). 10.1371/journal.pgen.1004061.g001 Figure 1 Genome wide associations (GWAs) of phosphorus (P) efficiency across traits, sites and years in 192 soybean accessions. GWAs across traits, sites and years (RAPA-HN/NJ2008 denotes the relative acid phosphatase activity under −P and +P obtained in Henan/Nanjing in 2008, RPC-HN/NJ2008 denotes the relative P concentration under −P and +P obtained in Henan/Nanjing in 2008). Chromosomes listed in the red font denote the six major significant SNP clusters associated with P efficiency identified by the GWAs. The color bands indicate the 10 major bi-parental QTL clusters identified in previous reports (green) or in the current study (yellow) by linkage mapping. Next, the SNPs identified by the GWAs were compared to the positions of QTL regions identified by linkage mapping in this study and previous reports (Figure 1). First, the results indicated that the GWAs identified more phenotype-genotype associations and provided higher resolution; however, linkage QTL mapping identified rare specific alleles that were undetectable by the GWAs. Second, most of the highly significant associated regions identified by the GWAs were also identified and validated by linkage QTL mapping (Figure 1). Finally, a highly significant SNP cluster (including four SNPs, P = 7.5×10−8) on chromosome 8, co-localized to qPE8, was found to be flanked by Satt089 and Sat_310. In addition, previous studies have shown that this locus is associated with P efficiency at the soybean seedling stage [9], [10]. Therefore, we analyzed qPE8 to identify the P-efficiency-related genes located within this region. At this stage of the study, however, it was not possible to discriminate among the candidate genes for P efficiency in this region. Therefore, we decided to investigate the region between Satt089 and Sat_310 (Figure 2A). 10.1371/journal.pgen.1004061.g002 Figure 2 Fine mapping and positional cloning of GmACP1. (A) A phosphorus efficiency quantitative trait locus (QTL) qPE8 was mapped to the interval between the markers Satt089 and Sat_310 on soybean chromosome 8 using 152 recombinant inbred lines (RILs). (B) This QTL was further delimited to an approximately 250 kb region on chromosome 8 using a natural population composed of the 192 accessions. (C) The arrow indicates the predicted gene between the markers Sat_233 and BARC-039899-07603, including the five candidate genes Glyma08g20700, Calcineurin B; Glyma08g20710, Phospholipase D; Glyma08g20800 and Glyma08g20820, Putative Phosphatase; and Glyma08g20830, Protein Phosphatase. (D) The GmACP1 gene model showing the allelic variation (frequency >5%) of the GmACP1 sequence. (E) Linkage disequilibrium (LD) plot between GmACP1 SNPs. The physical position of each SNP is shown above the plot. The magnitude of LD indexed by the D′ statistic is also shown. Red squares without numbers indicate complete LD (D′ = 1, P 5%) were tested, and the P-value for each site was estimated based on 1,000 permutations of the dataset under a mixed linear model (MLM). Markers were defined as being significantly associated with traits on the basis of a significant association threshold (−LogP≥2.00, P≤0.01). Associations for the relative acid phosphatase activity (RAPA) and relative P concentration (RPC) were performed on the 192 soybean accessions; only the significant association sites were reported. We performed a multiple linear regression (MLR) using SAS version 9.0 (GLM) (SAS Institute Inc., Cary, NC, USA) to detect the effect of GmACP1 polymorphisms and haplotypes in this natural soybean population. Furthermore, we used a 5-fold cross validation to estimate accuracy of the algorithm (MLR). We divided the phenotypic data into five segments, four of which were used for training, whereas one segment was omitted and used for testing. This was performed 1000 times; in the first instance, the first segment was used for testing and the remainder was used for training; then, the second segment was used for testing and the remaining segments were used for training, and so on. For the regional association mapping studies, the 192 soybean accessions were genotyped with 18 SSR markers and 22 SNP markers based on the soybean physical map (http://soybase.org/gbrowse/cgi-bin/gbrowse/gmax1.01/#search). The regional association analysis was performed as described above. Expression and purification of GmACP1 in E. coli The GmACP1 coding sequences (including the sequences of high/low P efficiency haplotypes) were was cloned into the pET28a expression vector and expressed in E. coli (BL21) cells, as described by Baldwin et al. [39]. The recombinant clones contained His tags at both the N- and C-terminal ends of the GmACP1 peptide. The E. coli cells were induced with IPTG to produce the recombinant protein and then lysed by sonication and run through a Ni-affinity column, according to the manufacturer's recommendations (GenScript, TM0217). The purity of the affinity-purified GmACP1 protein was confirmed on an SDS-PAGE gel. The GmACP1 protein concentration was measured using the Bradford Protein Assay Kit (Bio-Rad, Hercules, CA). The reactions contained 150 ng of GmACP1 protein and 1 µg of synthetic substrate ρ-nitrophenyl phosphate in 500 µL of NaAc-HAc solution (100 mM). Control assays were performed using an extract from E. coli containing only pET-28a. All reactions were performed in triplicate. The GmACP1 phosphatase activity was measured at 37°C for 10 min with pH values ranging from 3.0 to 9.0. Soybean hairy root transformation The GmACP1 overexpression vector was constructed using the Gateway technology with the Clonase II Kit (Invitrogen, Carlsbad, CA). The GmACP1 open reading frame (ORF) was amplified from the cDNA of soybean accession Kefeng No. 1, which is a variety with high P efficiency, using gene-specific primers (forward: 5′-GGGGACAAGTTTGTACAAAAAAGCAGGCTTCCAACATGTCTGGAACCGTGAT-3′, reverse: 5′-GGGGACCACTTTGTACAAGAAAGCTGGGTCTGGTCTA CTGGGAGGACT-3′). Two recombination reactions (BP and LR reactions) constitute the basis of the Gateway cloning technology. Finally, the amplified fragment was introduced into the pMDC83 vector according to the manufacturer's instructions and confirmed by sequencing. An empty vector was constructed using a single enzyme digestion method. The plasmid pMDC83 was digested with the restriction endonuclease KpnI and re-annealed following the removal of the fragment of the recombination site between attR1 and attR2. This process was confirmed by PCR (forward: 5′-GGTTGGCCATGGAACAGGTA-3′, reverse: 5′-GAGGACCTCGA CTCTAGAAC-3′) and sequencing analyses. Next, the expression vector containing the GmACP1 gene and empty vector were introduced into Agrobacterium rhizogenes K599 (kindly provided by Prof. Peter Gresshoff) through the freeze-thaw method. Soybean hairy root transformation was performed using the accession Kefeng No. 1, as previously described by Kereszt et al. [43]. The seeds were surface-sterilized with chlorine gas for 4 h prior to germination in vermiculite. Supporting Information Figure S1 Phosphorus (P) efficiency related QTL mapped on chromosome 8 using 152 RILs in 2006 and 2007. The black arrow indicates the P efficiency related QTL (qPE8) mapped to the location of GmACP1 on chromosome 8. (TIF) Click here for additional data file. Figure S2 Quantile-quantile plots of estimated −log10 (p) for phosphorus efficiency related traits from association analysis based on MLM with Q and K. HNAPA2008 and NJAPA2008 denote the phenotypic data obtained in Henan and Nanjing in 2008. The black triangles represent the P values expected under the null distribution. (TIF) Click here for additional data file. Figure S3 qRT-PCR of five candidate genes in two representative accessions with different phosphorus (P) efficiency values (the Y-axis denotes the gene expression levels). Nau3 (L) is an accession with low P efficiency, and Nau4 (H) is an accession with high P efficiency (gene annotation: Glyma08g20700, Calcineurin B; Glyma08g20710, Phospholipase D; Glyma08g20800, Putative Phosphatase; Glyma08g20820, Putative Phosphatase; and Glyma08g20830, Protein Phosphatase). (TIF) Click here for additional data file. Figure S4 Comparison of GmACP1 with other related proteins. Invariant residues are shown in bold, and other conserved residues are highlighted. Alignment of the amino acid sequences revealed two peptide motifs that are conserved in the active site of the HAD and DDDD superfamilies of hydrolytic phosphotransferases. The Asp residue predicted to be transiently phosphorylated during the phosphate transfer reaction is indicated with an arrow. The abbreviations and GenBank accession numbers for the acid phosphatase sequences analyzed are as follows: Phaseolus vulgaris putative phosphatase (ABP52095); Ricinus communis putative phosphatase (XP_002527425); Arabidopsis thaliana putative phosphatase (NP 173213); LePS2, Lycopersicon esculentum (AAG40473); and Zea mays putative phosphatase (NP_001151156). (TIF) Click here for additional data file. Figure S5 Linkage disequilibrium (LD) across GmACP1 in 192 soybean accessions. The bp positions of the polymorphisms in the alignment are shown on the left. Lower left triangle: P-value derived from Fisher's exact test. Upper right triangle: D′ values. (TIF) Click here for additional data file. Figure S6 Relationship between the acid phosphatase activity and the GmACP1 polymorphisms/haplotype in the two parents. The parent with low P efficiency (Bogao) had the haplotype “11TTGA” and the parent with high P efficiency (Nannong94–156) had the favorable haplotype “40TTAT.” ** indicates significance at P = 0.01. (TIF) Click here for additional data file. Figure S7 R2 of the relative acid phosphatase activity (RAPA) and relative phosphorus concentration (RPC) across three environments based on a 5-fold cross validation. The red arrows indicate the means calculated by multiple linear regression (MLR) analysis. (TIF) Click here for additional data file. Table S1 Twenty-eight annotated genes between Sat_233 and BARC-039899-07603, five of which were considered putative candidate genes related to phosphorus efficiency (bold). (DOCX) Click here for additional data file. Table S2 Fourteen previously reported soybean accessions with diverse phosphorus efficiencies. (DOCX) Click here for additional data file. Table S3 Summary of the 192 soybean accessions used. (DOCX) Click here for additional data file.
