Phenotypic and genetic variation in phosphorus-deficiency-tolerance traits in Chinese wheat landraces

Phosphorus deficiency is a major limiting factors for affecting crop production globally. To understand the genetic variation of phosphorus-deficiency-tolerance, a total of 15 seedling traits were evaluated among 707 Chinese wheat landraces under application of phosphorus (AP) and non-application of phosphorus (NP). A total of 18,594 single-nucleotide polymorphisms and 38,678 diversity arrays technology sequencing markers were used to detect marker-trait associations under AP and NP. Top ten genotypes with extremely tolerance and bottommost ten genotypes with extremely sensitivity were selected from 707 Chinese wheat landraces for future breeding and genetic analysis. A total of 55 significant markers (81 marker-trait associations) for 13 traits by both CMLM and SUPER method. These were distributed on chromosomes 1A, 1B, 2A, 2B, 2D, 3A, 4B, 5A, 5B, 6A, 6B, 6D, 7A and 7B. Considering the linkage disequilibrium decay distance, 25 and 12 quantitative trait loci (QTL) were detected under AP and NP, respectively (9 QTL were specific to NP). The extremely tolerant landraces could be used for breeding phosphorus-deficiency-tolerant cultivars. The QTL could be useful in wheat breeding through marker-assisted selection. Our findings provide new insight into the genetic analysis of P-deficiency-tolerance, and will be helpful for breeding P-deficiency-tolerant cultivars.


Background
Phosphorus (P) is an essential macronutrient for plant growth, and yet it is a major limiting factor for crop production globally [1,2]. Although total P in soil is abundant, it is largely unavailable for uptake by plants. Phosphorus readily forms complexes with metal cations, and microorganisms in the soil convert phosphate P into organic compounds [3][4][5], thus making P unavailable in soils. Phosphorus fertilizer is commonly used to overcome this problem. However, applying large quantities of fertilizer causes environmental concerns about fertilizers leaching into water systems, and accelerates the exhaustion of nonrenewable phosphate resources [6,7]. It has been estimated that the global P resources will be exhausted by the end of this century [2,8]. Thus, it is important to breed P-deficiency-tolerant varieties for sustainable agriculture.
Wheat (Triticum aestivum L.) is among the earliest domesticated crop plants, being cultivated 10,000 years ago in the pre-pottery Neolithic Near East Fertile Crescent [9,10]. Currently, wheat is the most widely cultivated food crop, contributing about a 20% of the calories consumed by humans. P is an important limiting nutrient for wheat growth and development. Wheat yield has been severely limited by P deficiency globally [11,12]. Hence, development of P-deficiency-tolerance wheat is critical. Landrace genotypes are an important resource for wheat improvement. Chinese wheat landraces have been shown to be enriched with genes and alleles that are tolerant or resistant to abiotic and biotic stress [13][14][15][16]. Therefore, screening tolerant genotypes and understanding the genetic basis of P-deficiency tolerance in Chinese wheat landraces could provide important insights for the breeding of tolerant wheat cultivars.
The genetic control of P-deficiency-tolerance traits has been investigated extensively using linkage mapping of biparental wheat populations. Quantitative trait loci (QTL) for P-deficiency-tolerance traits have been successfully identified using this approach [17][18][19][20]. An alternative method, which may have greater potential for improving P-deficiency tolerance, is to identify the allelic variations due to divergent selection pressures, within a large genetic diversity panel. Through genome-wide association studies (GWAS), such natural allelic variations for P-deficiencytolerance have been identified in Arabidopsis [21], Aegilops tauschii [22], soybean [23], but not yet in wheat. Thus, our objective was to identify such allelic variations in wheat through GWAS. These findings will provide new insights for the improvement of P-deficiency tolerance.
Phosphorus nutrition during the early growing stage is critical for wheat final yield. P deficiency during the early growing stage causes large reductions in tiller development and head formation, and plants cannot recover from this, even if sufficient P is later supplied [24,25]. Therefore, it is important to understand the genetic control of Pdeficiency-tolerance during the seedling stage, and to develop P-deficiency-tolerant cultivars. In this study, we first examined phenotypic variation in 15 P-deficiency-tolerant traits in a core of 707 wheat landraces at seedling stage. We then evaluated the P-deficiency-tolerance of these landraces to identify suitable wheat germplasm for future breeding of tolerant wheat cultivars. Finally, we performed a GWAS using 57,272 polymorphism markers to identify markertrait associations (MTAs).

Phenotypic variation
To evaluate variation in the phenotypic response to P deficiency, 707 wheat landraces were grown under both application of phosphorus (AP) and non-application of phosphorus (NP). Fifteen traits were evaluated to determine the effect of P deficiency and their genetic variation at the seedling stage. There was significant variation among genotypes for all traits (p < 0.001; ANOVA; Table 1). The P treatment had highly significant effects on all traits (p < 0.001; Table 1). Phosphorus deficiency had negative effects on almost all traits, except for the dry root-shoot ratio (DRS), fresh rootshoot ratio (FRS), and root diameter (RD). The coefficients of variation for the traits were 1.22 to 49.09% under AP, and 1.21 to 48.40% under NP (Table 2).  (Table 2). RD showed low heritability under both conditions. The Shannon-Weaver diversity index (H′), which shows the diversity of each trait, was 0.28 to 0.86 under AP and 0.28 to 0.87 under NP; it reflected moderate to high diversity for the traits, except for RT (Table 2).

Correlation analysis
Under AP, correlation coefficients ranged from 0.018 to 0.976 (Table 3; correlations under AP are shown below the diagonal). With the exception of DRS, FRS, and RD, the correlations between all pairs of the 15 traits were  Using '1' as diagonal line, the correlation of all tested traits under AP and NP are below and above the diagonal line, respectively;* and ** represent significance level of P < 0.05 and P < 0.01, respectively significantly positive. DRS was significantly negatively correlated with shoot dry weight (SDW), shoot fresh weight (SFW), shoot length (SL), and total length of shoot and root (TL). DRS was significantly positively correlated with the other traits, and was not significantly correlated with RD, RF, and RT. FRS was significantly negatively correlated with SL, and significantly positively correlated with the other traits. RD was significant positively correlated with FRS, and significantly negatively correlated with the other traits, except for DRS and root volume (RV). Under NP, the correlation coefficients ranged from 0.033 to 0.965 (Table 3; correlations under NP are shown above the diagonal). With the exception of DRS, FRS, and RD, all pairs traits were significantly positively correlated. DRS was significantly negatively correlated with RD, SDW, and SL, and significantly positively correlated with the others except for SFW. FRS was significantly negatively correlated with RD, SFW, and SL, and significantly positively correlated with the others. RD was significantly negatively correlated with all traits except root dry weight (RDW), RV, SDW, and SFW.
Under both AP and NP, root surface area (RSA) and total root length (TRL) had the highest correlations. With the exception of DRS, FRS, and RD, all pairs of correlations were significantly positive. Among the six root morphological traits (RF, root length (RL), RSA, RT, RV, and TRL), the correlations were moderately to highly positive. The other root morphological trait, RD, was negatively correlated with RF, RL, RSA, RT, and TRL under both conditions.

Principal component (PC) analysis
The PC analysis were performed using the relative trait values. The cumulative amount of phenotypic variation explained (PVE) by the first three PCs was 82.98% (Table 4). PC1 explained 57.84% of the phenotypic variation. With the exception of DRS, FRS, and RD, all of the other traits were important factors within the characteristic vector of PC1. PC1 represented plant biomass and root architecture, and can thus be defined as the biomass and root factor. PC2 (eigenvalue 2.53) explained 16.84% of the phenotypic variation. DRS and FRS, which influence P absorption in shoots and roots under Pdeficiency, were important factors within the characteristic vector of PC2. PC2 can thus be defined as the rootshoot ratio factor. The relative values of TL, SL, SFW, SDW, and RT for PC2 were negative, indicating that an increase in the root-shoot ratio will reduce TL, SL, SFW, SDW, and RT. PC3 (eigenvalue 1.24), for which RD was the only important factor, explained 8.29% of the phenotypic variation. The relative values of TL, RL, SL, FRS, DRS, TRL, and RT were negative for PC3, indicating that an increase in RD will correspond to a reduction in TL, RL, SL, FRS, DRS, TRL, and RT.

Screening for wheat tolerant genotypes
Using a weighting method [26], the synthesis value (S value) was calculated to evaluate wheat tolerance to P deficiency, and the P-deficiency tolerance index (PDTI) for each trait was calculated. The accessions with extremely high or low S values are listed in Additional file 1: Table S1. A high S value indicates high tolerance. The wheat landraces were classified into three groups, ranging from − 2.16 (accession AS661384) to 2.52 (accession AS661809) (Fig. 1). There were 173 accessions in the first group (S ≥ 0.5; classified as high tolerance). Group 2 (− 0.5 ≤ S < 0.5; intermediate tolerance) included 353 accessions. The remaining accessions (S < 0.5) were classified as sensitive. Accessions with higher S values also had higher PDTI (Additional file 1: Table  S1). This indicates that both of these indicators are effective for screening wheat landraces under P-deficiency.

Molecular markers and population structure
After excluding markers with missing data > 20% and minor allele frequency (MAF) < 0.05, 57,272 polymorphism markers were retained. Based on the 'Chinese Spring' physical map v1.0, 21,503 were mapped on A genome, 25,365 were mapped on B genome, and 10,404 markers were mapped D genomes.
A GWAS was also performed using settlement of mixed linear model under progressively exclusive relationship (SUPER) method in genome association and prediction integrated tool (GAPIT) to verify the results of TASSEL. Compared the results from TASSEL, a total of 55 significant markers (81 MTAs) identified by TAS-SEL were confirmed by SUPER method in GAPIT. These 55 significant markers were used for further analysis. Based on the linkage disequilibrium decay distance, we considered significant markers within a 5.98 Mb region to constitute a single QTL [28]. According to this, 25 QTL were detected under AP, 12 under NP, and 3 under both. Thus, 9 QTL were specific to NP.
Three QTL, QTL-2D-3, QTL-4B-2 and QTL-7A-1 were identified under both conditions (Additional file 3: Table S2). These three QTL, located on chromosomes 2D, 4B and 7A, were stably expressed under both conditions. Thus, nine QTL were specific to NP conditions. This indicates that these QTL occurred exclusively under P-deficiency stress.
Pleiotropy was revealed by GWAS. Seven QTL were identified for multiple traits. QTL-2D-3 was associated with six traits (RDW, RF, RFW, RSA, RV, and TRL) under NP, and with RF under AP. QTL-4B-2 was associated with six traits, including RF, RFW, RSA, RT, SDW, TRL under AP, and with RT under NP. This finding is supported by Pearson's correlation analysis, with the correlation coefficients among these six traits ranging from 0.578 (between SDW and TRL) to 0.976 (between RSA and TRL).

Discussion
As an important plant organ for P uptake and utilization, the root system has been used in screening P-deficiency-tolerant wheat genotypes [29]. The seedling stage often determines the traits of the mature stage, and the traits of the two stages are closely related. The QTL associated with N-absorption rate detected in a field experiment, and those associated with seedling character in greenhouse hydroponic culture, were colocated [30]. Further, the QTL of root hair on Chr. 2A and 6A were also associated with yield-related traits [31]. These results indicate that nutrient absorption and utilization during the seedling stage can affect the phenotype in the mature stage. The traits showed significant variation in their responses to P deficiency. Most showed moderate to high heritability (Table 2), revealing their potential for analysis in our next GWAS. Most of the traits were significantly correlated, and the six root traits were moderately to highly correlated. High correlations between wheat-root seedling traits have also been reported in previous studies [32][33][34][35]. Those studies revealed that seedling root traits are inherited together, and that it is difficult to independently select for one of these traits.
Here, the first three PCs explained more than 80% phenotypic variation and included all tested traits ( Table  4). The first PC mostly reflected the contribution of 12 traits (with the exception of DRS, FRS and RD). The second and third PCs reflected the contributions of rootshoot ratio and RD, respectively. We found that increasing root-shoot ratio led to reduction in aboveground biomass. DRS, FRS, and RD were higher under NP than under AP, whereas the other traits were lower under NP. Increases in the root-shoot ratio in response to stress have been reported previously [28,36]. Plants can increase their root-shoot ratio to promote P uptake and utilization in response to low-P stress [37]. 'Chinese Spring' was identified as sensitive to P deficiency (S = − 1.16), consistent with previous studies [19,26,38]. It indicated that using S value was a reliable method to screen wheat tolerant genotypes. Finally, ten extremely sensitive and ten extremely tolerant landraces were selected for further genetic analysis and for breeding tolerant cultivars (Additional file 1: Table S1).
The SUPER method is a power GWAS for identifying QTL as it extracts a small subset of single-nucleotide polymorphisms (SNPs) and use them in FaST-LMM [39]. This method can retain the computational advantages of FaST-LMM and increase statistical power [39]. A total of 81 of 116 (69.83%) MTAs by CMLM were confirmed by SUPER method. Comparing with CMLM, the number of MTAs identified using SUPER were greatly increased (Additional file 2: Fig. S1). However, quantile-quantile (Q-Q) plots from different algorithms revealed that CMLM was fitted better than SUPER (Additional file 2: Fig. S1). Previous studies revealed the similar results [40,41].
Pleiotropy was identified in seven QTL. QTL for different highly correlated traits may be located in the same chromosomal regions [42]. We found that all traits except DRS, FRS, and RD, were significantly positively correlated under both conditions. The six root morphological traits were moderately to highly correlated (Table 2). Pleiotropy has been observed in previous studies [19,32,33,43]. QTL-2D-3, located on Chr. 2D at 186.65 Mb, was associated with RDW, RF, RFW, RT, TRL, RSA, and RV under NP, and with RF under AP. In a previous study, QTL associated with P-use efficiency were identified on Chr. 2D [44]. We found that QTL-7A-1, located on Chr. 7A at 52.03 Mb, was associated with RT and RFW under AP, and with RT under NP. This indicates that QTL-7A-1 may regulate and control RT under both AP and NP, and root development under AP. Previous studies have found that Chr. 7A may play important roles in response to P deficiency [44,45].

Conclusions
We evaluated 15 traits in 707 Chinese wheat landraces, under application of P or non-application of P. Ten extremely tolerant and ten extremely sensitive accessions were selected as germplasm materials for further study. In total, 25 and 12 QTL were identified under AP and NP, respectively, while nine of 12 were specific to NP. In total, seven QTL showed pleiotropy, and several QTL had been previously identified. Our findings provide new insight into the genetic analysis of P-deficiencytolerance, and will be helpful for breeding P-deficiencytolerant cultivars.

Plant Germplasm
The 707 accessions used in this study were from a core collection of wheat landraces [46] originating from ten agroecological zones in China (Additional file 4: Table S3).

Glasshouse experiments and phenotypic data collection
All landraces were grown hydroponically in a greenhouse at the Triticeae Research Institute, Sichuan Agricultural University. A completely randomized design, each with three replications, was used in this study. The greenhouse environment, hydroponic system, and phenotypic date collection were as previously described [26,47,48]. Briefly, the AP and NP treatments contained the modified from Hoagland's nutrient solution [26,[47][48][49] with and without NH 4 H 2 PO 4 (1 mmol/L), respectively. Two sets of seedlings were grown for 3 day under AP, and were then grown for 12 day under AP and NP, respectively. Solutions were replaced every 4 day. After 12 day of growth, phenotypic data of all traits were gathered as our previous described method [28].

Phenotypic data analysis
To eliminate environmental effects, the best linear unbiased prediction (BLUP) values across three repetitions were conducted using the MIXED procedure in SAS [50]. The BLUP values for each trait were used to determine descriptive statistics, for ANOVA testing, and to obtain the H′ and Pearson correlation coefficients, using IBM SPSS Statistics for Windows 20.0 (IBM Corp., Chicago, IL, USA). The H 2 was calculated using the formula H 2 = V G /(V G + V E /r), where V G is the genotypic variance, V E is the environment variance, and r is the number of replications [51]. To screen for P-deficiency-tolerance, we used a weighting method to acquire the S value of each landrace genotype. The S value was calculated using the following S ¼ P k i¼1 riYi= P k i¼1 ri [26].

Genotyping and population structure analyses
Genomic DNA was extracted using the CTAB method. Genotyping-by-sequencing libraries (96-plex) were constructed via the two-enzyme method [52], and sequenced on the Illumina HiSeq 2500 system. SNP calling was done using the Tassel pipeline [53]. The physical distances of SNP were based on the Chinese Spring reference sequence v1.0 [54]. SNPs without physical distances were removed. The linkage disequilibrium K-number neighbor imputation method was used for imputation accuracies [55]. Finally, 18,594 SNPs were retained with missing data ≤20% and MAF ≥ 0.05. Besides, 38,678 DArT-seq markers from our previous study were also used for GWAS [46]. Population structure was evaluated using STRUCTUR E 2.3.4, implementing model-based Bayesian cluster analysis [56]. Based on the admixture model, population genetic clusters of K = 1 to K = 10 were estimated with 10,000 replicates for burn-in and 10,000 replicates for MCMC. Five runs were set for each K. The optimal K value was determined using STRUCTURE HARVESTER [57] implementing the Evanno method [27]. The optimal alignment of the five repeated runs was determined using CLUMPP [58].