Introduction

Drought is a major adverse climatic factor that severely restricts plant growth and sustainable agricultural development1. Breeding drought-tolerant cultivars has been a major objective of many crop improvement programs and is becoming an important strategy of adapting crops to climate changes2, 3.

Barley (Hordeum vulgare L.) is the fourth most grown cereal crop worldwide, and is characterized by having relatively high drought tolerance4. Consequently, barley is frequently used as a model crop for investigating the genetic basis of drought tolerance5, 6. Drought-adaptive traits derived from wild relatives and landraces of barley may contribute to increased yields and yield stability in improved cultivars under drought conditions7. The ability of different barley genotypes to adapt to drought stress resides in their genetic diversity8, 9. The exploration of genetic variation in genes related to drought response in natural populations may facilitate a better understanding of the molecular mechanisms of drought tolerance in barley.

Previous research identified the P5CS gene, encoding delta-1-pyrroline-5-carboxylate synthase (P5CS), as the main drought-responsive gene in barley10,11,12. In plants, P5CS is a rate-limiting enzyme involved in the biosynthesis of proline from glutamate, and proline is considered one of the most accumulated osmolytes for adaptation to drought stress13, 14. The P5CS protein is encoded by two genes, P5CS1 and P5CS2, that have different expression patterns and perform distinct functions during crop development. Studies have shown that the expression of P5CS1 is significantly induced by several abiotic stresses such as drought and salt14,15,16. However, P5CS2 is apparently a housekeeping gene that mainly functions in basic proline metabolism in Arabidopsis, Sorghum bicolor, and Zea mays 15,16,17. Transgenic plants with overexpressed P5CS1 significantly increased biomass under abiotic stress conditions compared with non-transgenic controls18,19,20. Arabidopsis P5CS1 T-DNA insertion mutants reduced P5CS1 transcription level, and proline synthesis, hypersensitivity to low water potential, and accumulation of reactive oxygen species16. Kesari et al.21 found that intron sequence variations in Arabidopsis P5CS1 promoted alternative splicing, produced a non-functional P5CS1 transcript variant, and demonstrated that P5CS1 was a gene under selection for environmental adaptation. Barley genotypes with drought tolerance manifested higher proline contents and higher expression levels of HvP5CS1 compared to drought-sensitive genotypes under drought stress22, 23. However, little is known about the underlying effect of genetic variations in HvP5CS1 on traits related to drought tolerance.

Single nucleotide polymorphisms (SNPs) and small insertions and deletions (indels) are the most common forms of nucleotide variations in natural populations, and may be significantly associated with phenotypic variations and plant adaptation to diverse environments24. SNPs are ideal molecular markers for genome mapping, marker-assisted selection and association analyses in plant breeding programs25,26,27,28. Compared to the analysis of SNPs in arbitrary sequences, careful analysis of SNPs in genes of specific interest presents a much greater value for association analysis between allelic variants with phenotypic differences and the subsequent marker-assisted selection of the associated traits29, 30. To date, allelic variations of HvP5CS1 in barley have not been systematically examined. Exploration of genetic variation in HvP5CS1 may provide a better understanding of the functions of HvP5CS1 and useful information to improve drought tolerance in barley.

In this study, EcoTILLING (a variant of Targeting Induced Local Lesions in Genomes) technology31 was used to detect allelic variation in the targeted region of HvP5CS1 in 287 barley accessions collected from diverse geographical areas, and to evaluate nucleotide diversity and the neutral test of detected regions. In addition, association analysis between the allelic variations and drought tolerance related traits was performed to validate gene function and provide a potential source of beneficial alleles for marker-assisted selection to improve barley growth under drought stress.

Results

Nucleotide polymorphisms in HvP5CS1

Across the amplified region of HvP5CS1, 41 natural variation sites were identified among the 287 accessions by EcoTILLING analysis and sequencing. The polymorphism information content (PIC) values of those sites ranged from 0.014 to 0.483, with an average of 0.243 per site (Table 1). These variations contained 16 SNPs and 25 indels, with a mean of one SNP per 67 bp and one indel per 43 bp in the genomic sequence. Of the 16 SNPs, only one was located in the exon region of HvP5CS1, which was a silent synonymous change, while eight were in the intron region and seven were in the 3′ downstream of the non-coding region. Of the 25 indels, eight were in the exon regions, 11 in the intron regions and six in the 3′ downstream of the non-coding region (Table 1 and Fig. 1).

Table 1 Polymorphism information content (PIC) values and positions in the gene for nucleotide polymorphisms in HvP5CS1.
Full size table
Figure 1
figure 1

Diagram of PCR amplification and distribution of SNPs in HvP5CS1. White boxes are the coding sequence of exons, black boxes are the 5′-UTR and 3′-UTR, and black lines are introns. The relative positions of the PCR products amplified for EcoTILLING are indicated as arrows and black bars. Black up arrows indicate changes in coding regions of DNA and do not affect the amino acid sequence of the protein product. Black down arrows indicate changes in the non-coding regions of DNA. Black bars refer to either insertions or deletions.

Full size image

Nucleotide diversity (π) of the targeted region of HvP5CS1 was 0.00917 across 287 barley accessions. π values were different among geographic regions, ranging from 0.00683 for European accessions (nine accessions) to 0.00997 for Arabian Peninsula accessions (14 accessions) (Table 2). To test whether the identified nucleotide polymorphisms were selectively neutral, D*, F* and Tajima’s D statistic were calculated. For 287 barley accessions, Tajima’s D value was 1.45687, but not significant (P < 0.05). However, the values of D* and F* were highly significant (P < 0.02). Apart from the Australian and Europe accessions, the other subgroups, including those from Africa, Arabian Peninsula, Middle East Asia and North East Asia, showed evidence of non-neutral variation (Table 2). In addition, analysis of the HvP5CS1 sequence in a sliding window plot for Tajima’s D value identified three regions with significant or marginally significant positive Tajima’s D values in the sliding windows between nucleotides 505–525, 705–745 and 825–865, respectively (Supplementary Fig. S1).

Table 2 Barley HvP5CS1 nucleotide diversity (π), haplotype diversity and selection (D* and F*, and Tajima’s D) statistics for each geographic region.
Full size table

The barley population was grouped into putative haplotype categories based on the cleaved banding pattern in evaluated gel-frames. For the 41 sequence-validated nucleotide variations across the amplified region of HvP5CS1, 13 distinct haplotypes were identified in 287 accessions (Supplementary Table S1), and eight of them (60%) were rare haplotypes with frequencies below 0.05. The HvP5CS1_H1 was the most-common haplotype with 78 accessions. The reference haplotype (HvP5CS1_H2) was the second most-common in amplicon, which included 61 accessions. The non-reference haplotype frequencies (NHF) ranged from 0.007 to 0.272, which differed significantly among the geographical regions of the tested accessions. Accessions from North East Asia had the highest NHFs while those from the Arabian Peninsula had the lowest (Supplementary Table S2).

Population structure and genetic diversity

Association analysis requires the population structure to be taken into account in order to avoid false positive associations. Analysis of genetic distance and population structure confirmed a significant population structure in the barley population. The Bayesian approach implemented in the Structure software revealed a steadily increasing curve for LnP(D) (Supplementary Fig. S2), therefore, it is difficult to deduce the optimal number of subgroups based on the values of LnP(D). In contrast, the ΔK value attained a clear maximum at K = 3 (Supplementary Fig. S2) and the three groups (or clusters) revealed relatively low levels of admixture with G1, G2 and G3 containing 71, 113 and 103 accessions, respectively (Fig. 2a).

Figure 2
figure 2

Genetic diversity of the 287 barley germplasm. (a) Population structure of the population. G1–G3 stand for three groups. The y-axis is the subgroup membership, and the x-axis is the accessions in three groups (G1, G2, and G3). (b) Cluster tree based on UPGMA. C1–C7 stand for cluster 1 to 7. The UPGMA tree is color-coded based on the results from population structure analysis (red G1, green G2, and blue G3); (c) Number of accessions in each of the three groups from different sources.

Full size image

The results from the unweighted pair subgroup method with arithmetic (UPGMA) tree analysis on SSR marker data were generally consistent with the results from the STRUCTURE analysis with a few exceptions (Fig. 2b). UPGMA cluster C1, C2 and C3 contained 108 accessions with 85 from the G3, 20 from the G2 and three from the G1. The Cluster C4 contained 69 accessions with 66 from the G1 and three from the G3. Clusters C5, C6 and C7 contained 110 accessions with 93 from the G2, 15 from the G3 and two from the G1. Based on this information, we selected K = 3 as the optimal subgroup number. Of the ten runs for K = 3, the run with the highest likelihood value was selected to assign the posterior membership coefficients (Q-matrix) to each accession.

Among the 287 accessions tested, 34 (89.47%) wild accessions were classified into group G1, 52 cultivars (63.41%) into G2, 167 landraces into the three groups with 33 (19.76%) in G1, 58 (34.73%) in G2, and 76 (45.51%) in G3 (Fig. 2c).

Association analysis between SNPs/haplotype and traits

Association analysis was performed to find possible associations between nucleotide variations in HvP5CS1 and drought tolerance related agronomic traits. Since 28 of the nucleotide variations showed either linkage disequilibrium (LD) within groups or rare alleles (frequency <5%), only 13 distinct variations were used for association analysis. Of which, five were significantly (P < 0.001) associated with at least one of six phenotypic traits, and four were highly significant (P < 0.0001) (Table 3). The SNP at position 305 bp in the amplicon of HvP5CS1 was significantly associated with days to heading (DH), total above ground biomass dry yield (BY) and grain yield (GY) at Tel Hadya location and plant height (PH) at Breda location, and explained 4.97%, 6.11%, 6.31% and 6.66% of the phenotypic variations, respectively. Another indel at position 735 bp exhibited significant associations with growth vigour (GV), GY and BY, and explained 6.83%, 6.65% and 6.62% of the phenotypic variations, respectively, at Tel Hadya, and 6.01% and 7.90% for GY and BY at Breda, respectively. Similarly, an indel at position 846 bp explained 6.98%, 7.71% and 7.71% of the phenotypic variations for GV, GY and BY at Tel Hadya, and 9.08% and 9.38% for GY and BY at Breda, respectively. The SNP at position 768 bp was significantly associated with spike length (SL), BY and GY, and explained 7.36%, 4.48% and 9.81% of the phenotypic variation at Tel Hadya, respectively, and 4.43% and 5.95% for SL and GY at Breda, respectively. The SNP at position 920 bp was significantly associated with DH, BY and GY at Tel Hadya and with PH at Breda, explaining 4.95%, 4.97%, 5.13% and 5.35% of the phenotypic variations, respectively. Therefore, the five SNPs significantly associated with drought tolerance related traits in this study may have potential as DNA markers to be used in marker-assisted selection for improving drought tolerance in barley after further validation.

Table 3 Significant associations between SNPs of HvP5CS1 and agronomic traits of barley.
Full size table

Haplotype-trait associations were restricted to haplotypes with a higher than 5% frequency of minor alleles. Two haplotypes were significantly (P < 0.001) associated with the phenotypic traits (Table 4), and contained desired alleles at four variation sites. The haplotype HvP5CS1_H1 showed highly significant (P < 0.0001) associations with GY and SL at Tel Hadya, and explained 11.89% and 7.55% of the phenotypic variations, respectively. Similarly, the haplotype HvP5CS1_H1 explained 5.23% and 5.30% of the phenotypic variations for GV and BY, respectively, at Tel Hadya, and 7.05% for GY at Breda. The haplotype HvP5CS1_H4 exhibited significant associations with BY at Tel Hadya and PH at Breda, and explained 5.00% and 5.92% of the phenotypic variations, respectively.

Table 4 Association results for significant (P < 0.001) haplotypes identified by the GLM for certain agronomic traits in 227 barley accessions.
Full size table

Discussion

Using transcriptome analysis, the P5CS gene has been cloned in various plant species17, 32, 33, but its allelic variation has only been characterized in common bean13. In this study, 41 natural variation sites were identified from 1065 bp DNA sequence of HvP5CS1 across 287 barley accessions. The average frequency of polymorphisms was 67 bp per SNP and 43 bp per indel of the genomic sequence. The SNPs were distributed unevenly along the DNA sequence with nine variation sites in the exon regions, 19 in the intron regions and 13 in the 3′ downstream of the non-coding region. Several studies have shown that a low level of polymorphisms in exon regions but a high level in non-coding regions of barley34,35,36,37. Genetic diversity analysis in a new P5CS gene from common bean identified that the frequency of SNPs in intron sequences was 1.5-fold higher than in exon sequences. All indels were exclusively found in the intron region13. A C-terminal glutamic-γ-semialdehyde dehydrogenase (GSA-DH) domain in the HvP5CS1 protein, which is from 217–270 and 394– 444 nucleotide acids in the gene region, was an extremely conserved domain without any nucleotide variation, which agrees with Chen et al.13. Compared to previous reports on barley38,39,40, nucleotide diversity (π = 0.00917) and haplotype diversity (0.823) of HvP5CS1 was relatively high in this study. As predicted, accessions from the Arabian Peninsula and Middle East Asia had higher nucleotide diversity (π) and haplotype diversity than those from the other five regions, which is consistent with previous reports41, 42. One plausible reason is that these geographic regions were the main distribution areas of wild barley43.

In this study, Tajima’s D value was not significant despites a high positive estimate of 1.45687, but the value of D* and F* was significant. These results may be attributed to the intermediate frequency of polymorphisms observed, which may result from balancing selection, diversifying selection or population subdivision13. In addition, analysis of the HvP5CS1 sequence in a sliding window plot for Tajima’s D value illustrated that three regions had significantly or marginally significantly positive Tajima’s D values in the sliding windows between nucleotide midpoint numbers of 505–525, 705–745 and 825–865. The three domains may be the targets of naturally-acquired protective responses of barley under natural environments.

Association analysis emerged as a powerful approach to investigate the role of sequence polymorphisms in phenotype variation in response to environmental stresses44. In the present study, five variations in HvP5CS1 were significantly (P < 0.001) associated with up to four drought tolerance related traits in barley. Of these, one at 305 bp of HvP5CS1 was in an intron region, while the other four at positions 735 bp, 768 bp, 846 bp and 920 bp were at 3′ downstream of the non-coding region. These variations may affect the expression of HvP5CS1. In our previous studies, the transcription level of HvP5CS1 in Tadmor (drought tolerant) was 1.22 fold that of WI2291 (drought sensitive) under drought stress22. Tadmor had positive alleles at the positions 305 bp, 735 bp and 846 bp in contrast with those in WI2291. Except for position 920, the other four sequence variations were significantly (P < 0.0001) associated with BY and GY at Tel Hadya (a favorable environment). Meanwhile two variations (at the positions 735 bp and 846 bp) were associated with BY and GY at Breda (a drought stress environment), which have marginally significantly positive Tajima’s D values in the sliding window. Haplotype-trait associations showed that two haplotypes, HvP5CS1_H1 and HvP5CS1_H4 carrying positive alleles of the two variations at positions 735 bp and 846 bp were significantly (P < 0.001) associated to drought tolerance related traits, and explained 5.00~11.89% of the phenotypic variations. Our previous research showed that Tadmor carrying HvP5CS1_H1 had high osmotic adjustments capacity and high yield stability, while Er/Apm and Express, both carrying HvP5CS1_H3, showed a low and constitutive osmotic capacity45. The variations at the positions 735 bp and 846 bp may provide useful markers as a selection tool to improve barley yield under stress conditions.

Several studies showed that GY was positively correlated with proline content in barley46 and wheat47, 48, and proline content had weak positive but non-significant correlations with DH, PH, SL and TKW in wheat48. Our study showed that GY had a significantly positive correlation with BY and TKW, but a significantly negative correlation with PH and SL in barley (Table 5), which agrees with several earlier reports49,50,51.

Table 5 Phenotypic scores from field trials for 227 barley accessions.
Full size table

This study identified 13 unique haplotypes based on 41 variations in HvP5CS1 among 287 barley accessions collected from 35 countries using EcoTILLING technology. Five polymorphisms and two haplotypes were significantly (P < 0.001) associated with drought tolerance related traits in barley. These polymorphisms can serve as DNA markers in breeding to improve drought tolerance in barley after further validation.

Materials and Methods

Plant materials

A set of 287 barley (Hordeum vulgare L.) accessions consisting of 167 landraces, 82 cultivars and 38 wild accessions were obtained from the International Center for Agricultural Research in the Dry Areas (ICARDA). Among them, 241 were originally collected from 35 countries in six adjacent geographic regions, including Middle East Asia, North East Asia, Arabian Peninsula, Africa, Europe and Australia, while the origin of the other 46 accessions is unknown (Table 6 and Supplementary Table S3).

Table 6 The geographic origins of the barley accessions used in the study.
Full size table

Genomic DNA was isolated from young leaf tissue of each accession using a modified cetyltrimethyl ammonium bromide (CTAB) protocol52, quantified using a spectrophotometer, and normalized to a concentration of 20 ng/μl.

Evaluation of drought tolerance related traits

As described by Varshney et al.53, 227 barley accessions were grown in the 2004–2005 grown season at two locations in Syria—Tel Hadya (36°01′N; 37°20′E, elevation 300 m asl, long term average annual rainfall = 338 mm) and Breda (35°56′N; 37°10′E, elevation 354 m asl, long term average annual rainfall = 269 mm). Tel Hadya represented a favorable condition and Breda represented a drought stressed environment. Field plots were arranged in a row and column grid (10 rows, 28 columns). Trials were sown on 3 m2 plots with 1.2 m wide (4 rows 30 cm apart) and 2.5 m long under rainfed conditions. Rainfall at the growing season was 303 mm at Tel Hadya and 268 mm at Breda. Phenotypic data were recorded for seven developmental and yield-related traits: days to heading (DH in days), growth vigour (GV 1 = good to 5 = bad), plant height (PH in cm), spike length (SL in cm), 1000 kernel weight (TKW in g), grain yield (GY in kg/ha) and total above ground biomass dry yield (BY in kg/ha) as described in Varshney et al.53. These phenotypic traits were tested for correlations using SPSS 16.0 software (SPSS Inc., Chicago, IL, USA). Phenotypic scores of measured traits are listed in Table 5.

PCR amplification and EcoTILLING assays

Natural variations in barley HvP5CS1 were screened using EcoTILLING techniques. The primers were designed according to the mRNA sequence of HvP5CS1 (GenBank:AK251855.1) with melting temperatures around 60 °C by using Primer 5.0 software (Premier Biosoft International, Palo Alto, CA, USA). The DNA region amplified by the primer pair is around 1065 bp, which includes a 3′-UTR and small coding sequence interrupted by three introns (Fig. 1). The forward primer labeled with IRDye700 at 5′-end and the reverse primer labeled with IRDye800 at 5′-end were synthesized by LI-COR Inc. (Lincoln, NE, USA).

The accession ICARDA IG 138264 (drought sensitive) was selected as a reference accession. PCR amplification of the target region was performed in a volume of 10 μl containing a mix of 10 ng each of reference and sample genomic DNAs. Heteroduplex formation, heteroduplex digestion, and denaturing polyacrylamide gel electrophoresis were performed as described by Xia et al.54. During electrophoresis, the LI-COR DNA Analyzer captured two images in IRD700 and IRD800 channels, respectively. Tiff images were manually scored using GelBuddy program55. Big dark bands of different sizes in both IRD700 and IRD800 channels were considered a polymorphic site (Supplementary Fig. S3). Total length PCR products from both channels should be equivalent to the fragment size of the undigested PCR product. Samples were grouped into putative haplotype categories based on the cleaved banding pattern in evaluated gel-frames.

DNA sequencing and statistical analysis

Once a polymorphism in a barley accession was identified, the corresponding DNA sample was amplified with the primers used in EcoTILLING. The resulting PCR fragment was directly sequenced to confirm the polymorphisms. Each polymorphic site was sequenced from more than one accession that showed the polymorphism in the same site to confirm that only two alleles segregated at the specific site. Multiple sequence alignment was conducted using ClustalW software56. Nucleotide diversity (π), haplotype diversity, Tajima’s D57 and D* and F* values58 were calculated using DnaSP v5.059. PIC was calculated as described in Kota et al.60.

Population structure analysis

The barley germplasm collection was structured into geographical groups in which some individuals are possibly related, so background genetic effects must be controlled to avoid spurious associations. Incorporating structure components as covariates in association analysis helps to reduce the false associations61.

All of the accessions were genotyped using 21 SSR markers assigned to seven chromosomes in barley (Supplementary Table S4). Amplification of the SSRs was carried out as described by Hayden et al.62. Amplified PCR products were analyzed in the LI-COR 4300 DNA analyzer (LI-COR, Lincoln NE, USA) and scored by comparing sizes between PCR products and a molecular weight ladder. Data from the 21 SSR markers were used to determine the population structure using the STRUCTURE software version 2.3.363, 64. The program was run for K ranging from 1 to 12 with ten replications per K, using the admixture model with 10,000 burn-in period and 100,000 iterations. The optimal number of subpopulations was determined by the ad hoc statistic ΔK based on the rate of change of the likelihood value65.

To validate the genetic structure and test for different models, genetic similarities for each pair of lines were calculated using the NTSYS-pc 2.11 software package (Biostatistics Inc., USA). Subsequently, a phylogenetic tree using the unweighted pair subgroup method with arithmetic means (UPGMA) based on the genetic similarities matrix was constructed using the MEGA 5.0 program66.

Association analysis

Association between SNP markers in HvP5CS1 and drought tolerance related traits was evaluated using a general linear model (GLM_Q) in the TASSEL v3.0 (http://www.maizegenetics.net/tassel), where the SNP was considered as a fixed effect, and the factor and matrix of subpopulation membership (Q matrix) were used as cofactors to account for population structure. The significance of associations between markers and traits was tested using an F-test. The association between a marker and a trait is represented by its R2 value, an estimate of the percentage of phenotypic variation explained by the marker.

As described above for single marker association analysis, haplotype-trait association test was also performed using GLM_Q in the TASSEL v3.0. Haplotypes with a frequency <5% were discarded to avoid biased associations. Only five haplotypes were identified for trait-haplotype association analysis. Accessions carrying the same haplotype were coded as “1”, whereas these accessions without the target haplotype were coded as “0” as described by Li et al.67. All other analyses were the same as described for single marker analysis.

Data availability

All data generated or analysed during this study are included in this published article and its Supporting Information files.