Genome-Wide Association Analysis of Aluminum Tolerance in Cultivated and Tibetan Wild Barley

Tibetan wild barley (Hordeum vulgare L. ssp. spontaneum), originated and grown in harsh enviroment in Tibet, is well-known for its rich germpalsm with high tolerance to abiotic stresses. However, the genetic variation and genes involved in Al tolerance are not totally known for the wild barley. In this study, a genome-wide association analysis (GWAS) was performed by using four root parameters related with Al tolerance and 469 DArT markers on 7 chromosomes within or across 110 Tibetan wild accessions and 56 cultivated cultivars. Population structure and cluster analysis revealed that a wide genetic diversity was present in Tibetan wild barley. Linkage disequilibrium (LD) decayed more rapidly in Tibetan wild barley (9.30 cM) than cultivated barley (11.52 cM), indicating that GWAS may provide higher resolution in the Tibetan group. Two novel Tibetan group-specific loci, bpb-9458 and bpb-8524 were identified, which were associated with relative longest root growth (RLRG), located at 2H and 7H on barely genome, and could explain 12.9% and 9.7% of the phenotypic variation, respectively. Moreover, a common locus bpb-6949, localized 0.8 cM away from a candidate gene HvMATE, was detected in both wild and cultivated barleys, and showed significant association with total root growth (TRG). The present study highlights that Tibetan wild barley could provide elite germplasm novel genes for barley Al-tolerant improvement.


Introduction
Aluminium (Al 3+ ) toxicity is considered as a major factor limiting crop production on acid soils [1,2]. Al 3+ rapidly inhibits plant root growth, thus preventing uptake of water and nutrients [2]. The initial and visual symptom of Al toxicity is the inhibition of root elongation and development of lateral roots.
There has been no commonly-recoginized standard criterion for evaluating plant Al tolerance. Many researches have used relative root growth (RRG) to estimate Al tolerance [3][4][5][6][7]. However, a weak correlation (R 2 = 0.17) was found between rice Al tolerance based on RRG of the longest root and RRG of the total root system [8]. In addition, an absolute root elongation was also used for Al tolerance evaluation in short-term Al stress experiments [9][10][11][12]. It was reported that utilization of different root index in QTL or genome-wide association analysis could be helpful to identify novel loci [13]. To make the association analysis more efficient, accurate determination of root parameters is crucial and essential. Famoso et al [8] developed a novel and high-throughput Al tolerance phenotyping platform, which consists of a root imaging system and root computer program. So far, this technology has been successfully applied in many research projects [14][15].
The Al tolerance machanism of cereals, such as wheat, barley, rice, maize and sorghum, has been extensively studied [4,[16][17][18][19][20][21]. The most common mechanism is secretion of organic acids, which may chelate Al in the rizosphere, thus preventing Al being absorbed by roots [22]. However, it has been reported that Al resistance in maize cannot be fully explained by root organic acid release [23]. Thus, multiple mechanisms for Al tolerance may exist in plants. Degenhardt et al. [24] reported that a decrease of Al 3+ activity is caused by root-mediated increase in rizosphere pH in an Arabidopsis mutant. Furthermore, a recent research reported that a T-DNA insertional mutant of XTH31, a gene controlling xyloglucan endohydrolase (XEH) and xyloglucan endotransglucosylase (XET) activities, gains the function of Al resistance by modulating cell wall xyloglucan content and Al binding capacity in Arabidopsis [25]. In short, the underlying mechanism of Al tolerance is complex and still not totally understood.
Altough barley is considered as one of the most sensitive species to Al toxicity among cereal crops, there is still a wide variation in Al tolerance among cultivars or genotypes [7]. So it is meaningful to investigate the underlying mechanism of barley Al tolerance. It was widely reported that the aluminum tolerance of barley was controlled by a dominant loci on chromosome 4H. Furukawa et al [4] and Wang et al [26] identified the same candidate gene HvMATE (HvAACT1) in cv. Murasakimochi and Dayton respectively, which is responsible for Al-induced citrate secretion and located on the long arm of chromosome 4H. However, Gruber et al [27,28] reported that over-expression of HvALMT1, which encodes an anion channel protein to facilitate organic anion transport, confers enhanced Al tolerance. Moreover, multigenic control of barley aluminum tolerance was reported in a few QTL studies [29,30]. Hence, barley Al tolerance is a complex trait, which may be quantitatively inherited. On the other hand, QTL analysis has the limitation that only the alleles segregated between the parents of DH or RIL population could be identified. Hence the genetic diversity of tolerant genotypes used in QTL analysis was limited. Compared with QTL mapping, GWA can always detect more alleles in one experiment. Because natural populations, especially Tibetan wild barley accessions used in the present study, possess wider genetic variation and may contain richer tolerant alleles in comparison with a DH population.
During the extended period of domestication, especially the modern breeding and intensified cultivation, genetic diversity of cultivated barley declined gradually and rare alleles disappeared. Barren genetic background became the bottleneck of breakthrough for breeding. The Tibet plateau is considered as one of the centers of cultivated barley domestication [31][32][33]. The Tibetan wild barley has been proved to be rich in genetic diversity of stress adaptation or tolerance [34][35]. Due to harsh environments in Tibet, special resistant mechanisms may be developed and some excellent alleles may be formed in Tibetan wild barley. In fact, there is clear evidence that Tibetan wild barley contains elite alleles of HKT1, HKT2 [36] and HvCBF4 [37], conferring enhanced salinity tolerance.
In the present study, four root parameters were precisely determined by using a root morphology analysis system, which is modified according to Famoso et al [8], consisting of a root scanning system WinRHIZO and root image analysis software RootReader2D. The relationships between these root parameters were analyzed. Subsequently, population structure was analyzed in whole accessions, and linkage disequilibrium (LD) decays of cultivated and Tibetan wild barleys were determined and compared. Finally, genome-wide association analysis (GWAS) was conducted by using the four root parameters and 469 DArT markers within or across the cultivated barley and Tibetan wild accessions.

Plant materials and growth conditions
A total of 110 Tibetan wild barley accessions and 56 cultivated cultivars were used for Al tolerance identification and association mapping. All wild barleys were kindly provided by professor Sun of Huazhong Agricultural University, China [36,37], and 56 cultivated barleys are widely planted in the different areas of China, but no one was from Tibet. Barley seeds were soaked in deionized water for 6 h at 20uC in the dark and then germinated on moist filter paper in an incubator (20/14uC, day/night) for 2 days. Then, nine seedlings (about 2 cm root length) were selected for both control and treatment, scanned by WinRHIZO [38], a root-measuring system, and images were acquired for precise analysis of root length. After scanning, seedlings were placed on the net cups floating on a 0.2 mM CaCl 2 solution, pH 4.3, in 50 ml centrifugal tubes as control. For Al stress treatment, seedlings were exposed to a 0.2 mM CaCl 2 solution, pH 4.3, containing 5 mM AlCl 3 . The roots were imaged for both control and Al 3+ -treated plants after 6 days treatment.

Measurement and calculation of root parameters
The images were analyzed by RootReader2D (http://www. plantmineralnutrition.net/rootreader. htm), and longest root length (LR) and total root length (TR) were recorded. To evaluate  Al 3+ tolerance for each genotype, root growth parameters were calculated as follows: Longest root growth (LRG) = LR T6d -LR T0d *, which represents absolute longest root elongation under Al stress; Relative total root growth (RTRG)~(TR T6d {TR T0d )= (TR CK6d {TR CK0d ); * T and CK represent Al 3+ treatment and control respectively, and 0d and 6d represent plants with initial and 6 days treatments, respectively.
The correlation analysis among these root parameters was conducted by SPSS (v19.0).

Population structure, cluster and linkage disequilibrium (LD) analysis
The samples of the whole-genome profiling were analyzed using the Barley PstI (BstNI) version 1.7 array with Diversity Arrays Technology (DArT P/L) in Australia (http://www.triticarte.com. au/content/barley_diversity_analysis.html). A total of 469 barley DArT markers (MAF.0.03) distributed over the whole genome were used to detect population structure using the STRUCTURE software v2.3.3 [39,40], performing ten independent runs, setting the number of clusters (k) from 1 to 10, with 100,000 MCMC (Markov Chain Monte Carlo) iterations in an admixture model. The largest value of statistic index Dk, which was the change rate of LnP(D) (from STRUCTURE output), was used as an indicator of true number of clusters (k) [41]. We calculated genetic distance with NTSYSpc (version 2.10e) and developed a neighbor-joining tree according to Nei [42].
Linkage disequilibrium (LD) between every two linear DArT markers was estimated by TASSEL (v3.0) [43] for Tibetan wild and cultivated barleys individually. The r 2 values, the squared allele frequency correlations, which is a measurement of the correlation between a pair of variables, and genetic distances were selected to display LD decay curve and fitted equation using Origin Pro (v8.0). LD decay distance in the whole genome was calculated as the genetic distance, when r 2 = 0.1 using the fitted equation.

Genome-wide association analysis (GWAS)
Association analysis between 469 DArT markers and four root parameters was performed using three approaches in all accessions (166) using TASSEL (v3.0) [43,44]. The first approach was Q method, in which Q matrix was included as a cofactor in the regression model to correct population structure. The second approach was k method, and a kinship, which represented the pair-wise relationship matrix, was estimated using software SPAGeDi [44,45], and considered as a cofactor in the regression model. The third approach was Q+k method, considering both population structure and kinship as cofactors. According to the Q2Q plot from the output of TASSEL, both k and Q+k methods were appropriate for the present study. Then, we conducted GWAS using k method in Tibetan wild and cultivated barleys, respectively. The Benjamini-Hochberg false discovery rate (BH-FDR) [46] of q-value = 0.01 was applied in association significance test. Manhattan plots were displayed by R software v2.14.2.

Root index distribution analysis
According to the polymorphism (two type: 1 and 2) of detected markers and subpopulation groups (Wild barley, w; cultivated barley, c) and the root index of accessions were classified into four groups in a box plot. Five detected markers (bpb-6949, bpb-9458, bpb-8524, bpb-0631 and bpb-8021) were analyzed independently. ANOVA analysis was conducted among four groups (w1, w2, c1 and c2) using LSD test.

Al tolerance evaluation based on root indexes and their relationships
A total of 110 Tibetan wild accessions and 56 cultivated barleys (Table 1) were evaluated for Al tolerance in response to 5 mM Al 3+ . A wide range of variations in LRG, TRG, RLRG and RTRG were observed ( Fig. 1a-d). Under Al 3+ stress, LRG and TRG were distributed around a mean of 1.4060.28 (SD) and 6.6861.70, ranging from 0.64-2.14 and 1.84-10.92, respectively. RLRG and RTRG were present with a mean of 0.37360.094 and 0.39760.107, ranging from 0.183-0.767 and 0.116-0.696, respectively. Weak correlations (less than 0.15) were found between root absolute elongation (LRG, TRG) and relative growth (RLRG, RTRG) ( Table 1), suggesting that absolute and relative growth could be different in evaluation of barley Al tolerance.

Genetic division between Tibetan wild and cultivated barleys
Population structure was analyzed using 469 DArT markers. The k2Q1 group, consisting of forty-eight Tibetan accessions with Q value .90%, could be distinctly detected when k was 2 or more ( Fig. 2a, b). The result was consistent with the data from the cluster analysis (Fig. 3), in which Tw1 consisting of 48 Tibetan wild accessions showed significant genetic division from the remaining accessions. It is suggested that wide genetic diversity is present in Tibetan wild barley.
To determine the number of subpopulations suitable for association analysis, Dk criterion was applied, and the most evident level of differentiation was observed with k = 3 (Fig. 2c). The k3Q1 consisted of 48 Tibetan accessions, and the same was for k2Q1. Interestingly, k3Q3 (Q.60%) consisted of 32 cultivated barleys (Fig. 2b), revealing the existence of considerable genetic diversity between Tibetan and cultivated barleys. Cluster analysis showed a similar result that 35 cultivated barleys belonged to C1 group (Fig. 3).
Linkage disequilibrium (LD) decay of genetic distance in the genome of Tibetan and cultivated barleys was 9.30 cM and 11.52 cM (r 2 = 0.1) (Fig. 4), respectively. Thus 469 DArT markers used in the present study could cover the whole genomic region, and are sufficient for genome wide association analysis.

Identification of root parameters mediated-Al tolerance loci through GWAS
To identify Al tolerance loci, a genome wide association study was conducted using 469 DArT markers and four root growth  parameters (LRG, TRG, RLRG and RTRG) (Fig. 5). When Q method was applied in association analysis across all accessions, 97 marker-trait associations, most of which belonged to TRG and LRG groups, were detected (p,0.01). However, most of these associations were identified as false positive results by Q2Q test (result by TASSEL), because the P values from TRG and LRG groups deviated from the expected value (Fig. 6a). P values from k model and Q+k model were close to expected values, indicating both models are suitable for association analysis (Fig. 6b, c). Eighteen and 16 associations were detected (p,0.01) using the k model and Q+k model, respectively, and 15 associations detected using both models.
Association analysis was also performed in subgroups. Twentyseven and 24 associations were detected in the cultivated and wild groups, respectively (p,0.01). When the threshold of significant test was set at p,0.001, only three regions were detected. In the cultivated group, bpb-6949 on Chr.4H and bpb-0631 on Chr.1H explained 25.6% and 23.1% of the phenotypic variation, respectively. In the wild group, bpb-9458 and bpb-0590 on Chr.2H both explained 12.9% of the phenotypic variation (p,0.001). It was worth mentioning that three root parameters (LRG, RLRG and TRG) were associated with bpb-8524 and bpb-6707 on Chr.7H, suggesting these markers were powerfully correlated with Al tolerance, and explaining 9.7% and 9.5% of phenotypic variation. Bpb-8021 on Chr.3H was detected as the most significant marker across all accessions using both k and Q+k methods, but certified as false positive association by following distribution analysis.

Root index distribution of the detected five markers
Based on marker polymorphisms, the distribution of root index were examined within cultivated and Tibetan barleys based on markers bpb-6949, bpb-9458, bpb-8524, bpb-0631 and bpb-8021 (Fig. 7a-e). W1-1 (95 accessions) showed higher TRG than w1-2 (15 accessions), which was similar to the result in the cultivated group (Fig. 7a). However, bpb-6949 was detected in all accessions and cultivated group by GWA, but not in the wild group. It may be explained that this Al 3+ -tolerant locus commonly existed in barley accessions, but its effect may be reduced or covered by other loci in the Tibetan group due to abundant genetic diversity.
W2-1 (10 accessions) and w3-1 (7 accessions) showed higher RLRG than that the other three groups, whereas no significant difference in RLRG was found between groups in the cultivated group (Fig. 7b, c). Moreover, only one accession was presented in haplotype c3-1 with low RLRG. These results confirmed that bpb-9458 and bpb-8524 associated with Al 3+ -tolerant genes were Tibetan group-specific.
C4-2 showed lower TRG than the other three groups, indicating bpb-0631 was cultivated group-specific, which was consistent with the results obtained by GWA analysis (Fig. 7d). Although bpb-8021 showed the most significance in GWA analysis across all accessions, no significant difference was found among groups within this marker, suggesting that it is a false positive association (Fig. 7e).

GWAS results are affected by phenotyping methods
Many researchers usually consider relative longest root growth as an indicator for Al tolerance evaluation [3][4][5][6][7]. However, it was reported that the application of different phenotypic indexes directly impacted the significance of QTLs detected in rice Al tolerance [13], which was confirmed in the current study, that different loci were detected as associated with Al tolerance by using different root parameters in GWAS. Although there were significant corrections between LRG and TRG, and between RLRG and RTRG, most of loci detected are different.
Therefore, the question may be raised as to which index could be recognized as reasonable indicators for Al tolerance evaluation. In the present study, different subpopulations needed different root parameters to identify Al tolerance loci. When absolute root growth (LRG, TRG) indexes and relative root growth (RLRG, Figure 5. GWA analysis of Al tolerance within and across 56 cultivated and 110 Tibetan wild subgroups. Four root parameters were applied for evaluating Al tolerance: TRG (6), RTRG (%), LRG (+) and RLRG (%;). GWA analysis was firstly conducted using three different methods: Q, k and Q+k methods. Then, k method was selected to perform GWA analysis in cultivated and Tibetan wild subgroups individually. Significant associations were identified using criterion of 2log10(P) .2 or 3. doi:10.1371/journal.pone.0069776.g005 RTRG) indexes were used, 22 and 5 loci were detected respectively in the cultivated group (Table 2, Fig. 5). The different results were obtained in the Tibetan group, with 6 loci being detected using absolute root growth (LRG, TRG) indexes and 18 loci using relative root growth (RLRG, RTRG) indexes. The number of detected loci showed that absolute root growth (LRG, TRG) indexes are more suitable for the cultivated group, and relative root growth (RLRG, RTRG) indexes for the Tibetan group. However, the underlying reason is not known. Furthermore, the most significant locus bpb-6949 in the cultivated group Figure 6. Quantile-quantile (Q-Q) plots of estimated 2log 10 (P). Q2Q plots were displayed in four marker-trait association analysis using three models: a Q method, b k method, c Q+k method. The black line is the expected line under the null distribution. Red symbol represents the observed P values for LRG, yellow for TRG, blue for RLRG and green for RTRG. doi:10.1371/journal.pone.0069776.g006 was associated with TRG index, while bpb-9458 was detected using RLRG index in the Tibetan group. This result suggests that Al tolerance strategies differ among barley subpopulations and is mediated by developmental types of roots.

Tibetan wild barley provides elite germplasm for barley improvement in stress tolerance
Genetic differentiation is present between Tibetan wild and cultivated barleys. This is well supported by population structure partition (k = 2) and the genetic division of Tw1 group in cluster analysis in the current study ( Fig. 2a and 3). Therefore, the Tibetan wild barley could provide elite germplasm for barley improvement in abiotic and biotic stress tolerance to fight against unpredictable climate changes in the future.
Structure analysis and Dk showed that three estimated populations were present in all 166 accessions (Fig. 2b, c). A subpopulation k3Q3 group consisting of 32 cultivated cultivars, was segregated from k2Q2 group. It may be suggested that cultivated barley gradually developed and engendered its own novel physiological or molecular traits or Al tolerant strategies during its expansion into specific environments. A good example is that an Al-resistant cultivar Murasakimochi acquires enhanced Al tolerance by a 1 kb insertion in the upstream region of HvAACT1 [47]. This insertion happened during the expansion of barley cultivation onto acid soil area and was identified to be originated from Eastern Asian region [47]. However, there is no such insertion in all wild barley accessions used in the present study (Data not shown).
LD decay analysis showed that DArT markers used in this study could cover four times as the size of the barley genome, indicating the density of DArT markers sufficient for GWAS. LD decayed more rapidly in Tibetan barley (9.30 cM) than cultivated barley (11.52 cM), thus the results of GWAS in Tibetan wild group can provide higher resolution of fine mapping and gene discovery than those in cultivated barley.
In short, it may be concluded that Tibetan wild barley is appropriate and efficient for genome-wide association by which we may detect rare elite alleles contributing to elevated Al tolerance in barley.
Novel Al-tolerant loci are detected by GWAS GWA mapping by using 110 Tibetan wild and 56 cultivated barleys genotyped with 469 DArT markers identified four significant regions, all of which are subpopulation-specific.
In the cultivated group, bpb-6949 was localized at 0.8 cM away from a major QTL in Chr.4H and a candidate gene HvMATE, which is considered as a major gene conferring barley Al tolerance [26,48]. This marker was also detected (p,0.01) when all barley genotypes were examined, suggesting that it is a common locus related to Al tolerance in barley.
However, the genetic diversity of tolerant genotypes used in a previous QTL study was limited [49][50][51][52][53]. Only the alleles segregating between parents of DH population could be identified by the QTL analysis. HvMATE (HvAACT1) was identified by using Dayton and Murasakimochi [4,26]. The two genotypes, as well as WB229 and FM404 were used as tolerant parents in the previous QTL studies. All are cultivated barleys. In the present study, GWAS was also performed in the Tibetan group and novel loci were identified on Chr.2H and Chr.7H. This result strongly suggests that different Al tolerant mechanism exists in Tibetan  Table 2. Lists of DArT markers with a significant marker-trait association using different model (k and Q+k methods) and within subgroups.