Marker-Trait Associations of Agronomic Traits in Barley (Hordeum vulgare L.) under Semi-Arid Conditions

Abstract

Heading date, plant height, and spike-related phenotypes are key traits that affect the yield potential and adaptation of barley to Mediterranean environments. These traits are controlled by a complex network of genes as well as environmental effects. Marker-trait associations (MTAs) were performed using a set of 361 barley genotypes, evaluated phenotypically for three years under semi-arid conditions. In total, 24 MTAs were detected using the recently developed GWAS method (BLINK) at –log p ≥ 5. These included one MTA with awn length (Awl) on chromosome 2H; twelve MTAs with heading date (HD) distributed over all chromosomes except 4H; three MTAs with grains per spike (GPS) on chromosomes 1H, 2H, and 6H; six MTAs were detected for plant height (PH) on 2H (2), 3H, 4H, 6H, and 7H; and two MTAs with spike length (SL) on 3H and 7H. The results showed novel and known associations between specific SNP markers and heading date, plant height, and spike-related traits, suggesting that these markers can be used in breeding programs to improve these traits in barley. This study provides valuable information for breeders and geneticists working to develop new barley varieties that are better adapted to semi-arid climates and have improved yield and quality characteristics.

1. Introduction

Barley (Hordeum vulgare L.) is a cereal grain that is widely cultivated around the world. It is the fourth most important cereal crop after rice, wheat, and corn. Barley is a diploid species with a genome size of approximately 5.1 Gb [1], which is larger than that of rice. As the world’s population is expected to reach 9.7 billion in 2050, the total global food demand is projected to grow by between 35% and 56% from 2010 to 2050. When considering the impact of climate change, the figures for total food demand change slightly to a range of 30% to 62% [2]. However, abiotic stress, particularly drought, poses a significant challenge to achieving this target. Terminal heat stress during grain filling is a major constraint for barley production in Tunisia. To mitigate this problem, breeding efforts have focused on developing early heading cultivars with reduced plant stature, such as dwarf barley, which can escape the high temperatures during grain filling [3,4].

Understanding the genetic factors that influence agronomic characteristics is a significant scientific challenge in enhancing crop improvement. Many of the traits that are crucial for crop yield are complex and are affected by both genetic and environmental factors [3].

As genome sequence data becomes more readily accessible and technology for identifying and analyzing genetic markers improves, the use of linkage disequilibrium (LD) in genome-wide association studies (GWAS) has become more popular in crops, particularly in wheat [5,6,7], barley [8,9], and rice [10,11], for determining the genetic factors that contribute to complex trait variation. In barley, association mapping studies have been used for fast identification and scoring of genetic markers for a wide range of agronomic traits, including heading date [12,13,14], plant height [15,16,17] and spike-related traits [18,19]. The heading date determines the time of flowering and grain filling, which is a crucial factor in yield formation. Early heading enables the crop to escape harsh weather conditions, pests, and diseases, while late heading can lead to yield losses due to unfavorable environmental conditions. Heading date is regulated by various environmental and genetic factors, including photoperiod, vernalization, and developmental stage. One of the central regulators of heading date in barley is the Ppd-H1 gene, which promotes flowering under long days [20], while the Ppd-H2 gene is responsible for acceleration of flowering in response to short photoperiods [21]. HvCO1 is one of the barley homologs of the Arabidopsis CONSTANS gene and acts in parallel to Ppd-H1 to promote flowering [22]. Overexpression of HvCO1 up-regulates HvFT1 (an ortholog of the Arabidopsis FLOWERING LOCUS T gene) and accelerates heading. HvCEN/eps2 [23] is another paralog of HvFT1 that has a large effect on time to flowering, particularly in Mediterranean regions. The vernalization pathway is regulated by a different set of genes. Vrn-H1 [24] is closely related to the Arabidopsis gene APETALA1 and regulates the transition to the reproductive stage [25]. Allelic variation at Vrn-H2 [26] arises from the presence or deletion of a cluster of three genes (ZCCT-H), which are characterized by a putative zinc finger and a CCT-domain. Vrn-H2 is considered to play the role of repressor of flowering time by inhibiting expression of Vrn-H3 [27]. Vrn-H3 was shown to correspond to HvFT1, the ortholog of the Arabidopsis FT gene [21,28]. This gene plays a central role in promoting flowering as an integrator of the vernalization and photoperiod pathways in barley.

The decrease in the height of plants has been a crucial factor in enhancing crop productivity. The effectiveness of a breeding initiative depends on the identification of genes that cause semi-dwarfism. The semi-dwarfing sdw1/denso gene is among the most significant genes regulating plant height in barley and has been integrated into various high-yield modern cultivars [29,30]. This gene is located on chromosome 3H and linked to SNP (11_10754) [31]. Barley spike traits are also determined by the presence of awns. Awned varieties of barley produce higher yields compared to awnless varieties, especially in warm and dry conditions [32]. The awns inside the spike play a crucial role in photosynthetic activity [33]. Previous studies have identified QTLs for awn length located on 3H [34], 6H [35], and 7H [36].

The aim of this study was to identify the marker-trait associations of single nucleotide polymorphism (SNP) markers with heading date, plant height, and spike-related traits using 361 barley accessions evaluated for three years under a semi-arid climate in Tunisia. Therefore, improving these traits through breeding programs can lead to increased barley yield, making it an essential consideration for barley farmers and researchers.

2. Materials and Methods

2.1. Plant Materials and Phenotyping

A collection of 361 barley genotypes, including 16 Tunisian old landraces and exotic accessions (from North Africa, the Middle East, and Ethiopia) derived from the U.S. National Plant Germplasm System (NPGS), were used in this study, including two-rowed (27%) and six-rowed (73%) row types with spring, winter, and facultative growth habits. During three cropping seasons: 2017/2018 (E1), 2019/2020 (E2), and 2021/2022 (E3), field trials were conducted in El Kef, Tunisia, characterized by a semi-arid climate with an average annual rainfall of 380 mm. The experiment was arranged in a randomized complete block design (RSBD). Each accession was planted in two replicate plots, each comprising two 2.5 m long rows planted with 65 seeds with an inter-row spacing of 25 cm. The traits considered in this study are as follows: heading date (HD): recorded as the number of days when 50% of the ears in a plot emerged; plant height (PH): measured in cm from soil surface to tip of the spike (excluding awns); spike length (SL): measured in cm from the base of the spike to the tip of the terminal spikelet; awn length (Awl): the length of the awn in the central spikelets; and number of grains per spike (GPS).

2.2. Statistical Analysis

Descriptive statistics, analysis of variance, and heritability were checked using the variability package in R software. Phenotypic correlations matrix between traits were analyzed using the ggpair function from the GGally package in R (an extension package to ggplot2) [37]. The effects of genotype, replication, location (environment), and genotype by location interaction were tested via analysis of variance (ANOVA) using the anova_joint function of the metan package [38]. The estimate of the response variable for the ith genotype in the jth environment using the additive main effect and multiplicative interaction (AMMI) model, is given as follows:

yij = μ + αi + τj + ∑k = 1pλkaiktjk + ρij + εij

where λk is the singular value for the kth interaction principal component axis (IPCA); aik is the ith element of the kth eigenvector; and tjk is the jth element of the kth eigenvector. A residual ρij remains, if not all p IPCA are used, where p ≤ min (g⁻¹; e⁻¹).

2.3. Phylogenetic Tree, Kinship, and Principal Component Analysis (PCA)

The population structure was estimated using the similarity matrix (kinship matrix) and the principal component analysis (PCA). Both analyses were constructed using the GAPIT package [39].

2.4. Association Mapping (AM)

The genotyping data of barley accessions were from the sequencing experiment, BarleyNB_9K (Platform: Infinium 9K) including 6913 SNP markers, obtained from the T3/barley database. Markers with a minor allele frequency (MAF) of less than 5% were removed, resulting in 6118 SNPs used for genome-wide association studies (GWAS). GWAS were conducted using Bayesian-information and linkage-disequilibrium iteratively nested keyway (BLINK), a statistical method implemented in the GAPIT3 R package [40]. The kinship model (K) and the first three principal components (PCs) derived from GAPIT were considered for the association mapping. Manhattan plots were created to display the distribution of SNPs across the chromosome. The significance levels used in the GWAS analysis on these plots are represented by both Bonferroni and FDR thresholds (p = 0.01) and (p ≤ 0.05), respectively.

3. Results

3.1. Descriptive Statistics

The research took place in El Kef, Tunisia, which experienced a low and inconsistent distribution of rainfall during three cropping seasons: 2017/2018, 2020/2021, and 2021/2022. Due to this irregular temperature and rainfall pattern, particularly during the grain filling stage, each year was considered as a distinct environment. Combined descriptive statistics over three years are summarized in

Table 1. Descriptive statistics of the studied traits combined over three years.

Traits	Mean	Min.	Max.	σ2g	σ2p	σ2e	GCV	PCV	ECV	H2	GA
HD	130.5	110	155	69	74.2	5.2	6.3	6.6	1.7	0.92	16.5
PH	74	40	106	67	99.8	32.7	11	13.4	7.7	0.67	13.8
SL	7.5	4	16.2	2.9	3.6	0.7	22.8	25.4	11.2	0.80	3.1
GPS	53.7	18	93	524	540	16.3	42.6	43.2	7.5	0.96	46
Awl	11.9	0	21	5	6.5	1.4	19	21.4	10	0.78	4.1

σ²g: genotypic; σ²p: phenotypic variance; σ²e: environmental variance; GCV: genotypic coefficient of variation; PCV: phenotypic coefficient of variation; ECV: environmental coefficient of variation; H²: broad-sense heritability; GA: genetic advance.

. A wide range of variation was observed for most of the studied phenotypic traits. Phenotypic variances (σ²p) of all studied traits were greater than those of genotypic variances (σ²g). All traits had a PCV coefficient that was higher than the GCV coefficient, but the values were similar to each other, indicating that there is a small range of genetic variation and little impact from environmental factors. The estimated heritability for different traits ranged between 67% and more than 96%. In PH, the broad sense heritability was medium (67%) and genetic advance was moderate (13.8%), indicating the existence of both additive and non-additive gene effects. SL and Awl showed a high heritability of 80% and 78%, respectively, but a very low genetic advance ranging from 3.1% and 4.1%, respectively, suggesting the non-additive gene effect of these traits. High heritability (92%–96%) accompanied by moderate to high genetic advance (16.5%–46%) in the case of HD and GPS would suggest the existence of additive gene actions.

3.2. Analysis of Variance Using AMMI Model

The ANOVA analysis combined with the AMMI model for the 361 barley accessions evaluated over three years showed highly significant (p < 0.001) differences among environments, genotypes, and GE interactions for all traits (

Table 2. Analysis of variance for studied traits using an additive main effect and multiplicative interaction (AMMI) approach.

Source		ENV	GEN	GEN: ENV	Residuals	Total	%IPC1	%IPC2	Variation Explained (%)
Df		2	360	720	1080	2885	361	359	ENV	G:E	GEN
Mean Square	HD	6019.2 ***	212.3 ***	11.7 ***	0.9	36.8	59.5	40.5	12.3	8.6	78
	PH	71976.6 ***	234.0 ***	70.8 ***	6.2	117.2	51	49	50.1	17.8	29.3
	SL	211.4 *	9.5 ***	1.2 ***	0.3	2.1	58.1	41.9	8.3	17	66.8
	Awl	757.2 *	16.7 ***	2.2 ***	0.8	4	57.8	42.2	15	15.7	60
	GPS	7207.2 ***	1586.5 ***	35.2 ***	6.9	223.1	60.2	39.8	2.3	4.1	92.3

All genotype and genotype-by-environment interactions are significant at *** p < 0.001 and * p < 0.05.

). The G:E interaction was divided into two interaction principal component axes (IPCA). For all traits, IPCA1 explained more than 57% of the total G:E interaction sums of squares percentage, except for plant height, for which it explained 51%. The total variation explained by the main effects of genotypes ranged from 29.3% for PH to 92.3% for GPS. Environment had the most important effect on plant height, accounting for 50.1% of the total treatment variation.

3.3. Population Structure Analysis

Based on SNP marker data, a phylogenetic tree was constructed using 361 barley accessions. Two subpopulations associated with row type could be clearly segregated (Figure 1A). The relationship between accessions is shown in the heat map of the kinship matrix (K-matrix) (Figure 1B). The kinship matrix, which is based on SNP markers, can also be referred to as co-ancestry or half-relatedness. This matrix helps to illustrate the coefficients of co-ancestry among individuals, with a stronger red color indicating a higher degree of relatedness. Principal components analysis (PCA) is the most commonly utilized approach to address population structure in genetic association studies conducted on unrelated individuals [41]. The projection of genotypes against the first two PCs shows a clearly structured population according to row-type phenotype (2 row and 6 row) (Figure 1C).

3.4. GWAS Analysis for Markers Trait Association

A set of 6118 mapped SNP markers were used for GWAS, which were evenly distributed over all seven chromosomes. The number of markers varied among chromosomes, with a minimum of 600 SNPs on chromosome 1H and a maximum of 1221 SNPs on chromosome 5H. The highest number of loci detected were on 2H (6) and 7H (6), and the lowest was found on 4H (1). The expected and observed p-values were displayed using a Q-Q plot matched for all traits, then diverged to indicate a positive association. Thus, the association analysis is reliable and not likely to produce false negatives (Figure 2). Marker trait association (MTA) analysis using 361 barley accessions identified 23 markers associated with the five studied traits studied over three years (Figure 2,

Table 3. List of identified SNP markers associated with evaluated traits.

Traits	SNP	Chr	Position (bp)	P Value	Effect
Awl	12_20326	2H	54584401	1.11 × 10−8	−0.77
HD	11_21140	1H	539810979	2.27 × 10−7	−1.5
	12_30872	2H	29124597	3.76 × 10−6	1.32
	SCRI_RS_174077	2H	746743803	7.92 × 10−6	1
	SCRI_RS_159598	3H	23309866	6.81 × 10−6	1.27
	SCRI_RS_214153	5H	572518499	9.79 × 10−8	−1.11
	11_11489	5H	648004337	5.49 × 10−7	−1.11
	SCRI_RS_209824	6H	30634644	2.42 × 10−7	1.1
	SCRI_RS_236959	6H	146371783	6.88 × 10−7	−1
	SCRI_RS_169904	7H	55242186	2.55 × 10−5	−1.37
	12_30026	7H	560756345	1.84 × 10−6	1.64
	SCRI_RS_222863	7H	599463445	1.30 × 10−8	−1.44
	SCRI_RS_204483	7H	624191837	3.38 × 10−7	1.22
GPS	SCRI_RS_198546	1H	383207519	1.00 × 10−7	−4.65
	12_30896	2H	652032221	1.40 × 10−20	8.55
	SCRI_RS_170058	6H	123873661	3.01 × 10−10	−3.48
PH	SCRI_RS_142400	2H	132659897	1.89 × 10−6	−1.89
	SCRI_RS_135355	2H	687098008	6.65 × 10−6	−1.24
	12_30921	3H	673188759	1.57 × 10−6	−1.32
	11_10639	4H	544748028	2.05 × 10−6	1.48
	SCRI_RS_154574	6H	17007008	3.78 × 10−9	−1.7
	SCRI_RS_168494	7H	60197487	1.62 × 10−7	−1.8
SL	SCRI_RS_826	3H	590142670	5.66 × 10−7	0.32
	SCRI_RS_198005	7H	647499514	1.70 × 10−6	0.43

). For HD, 11 MTA were identified on all chromosomes except chromosome 4H, with an effect ranging from −1.5 to 1.64 days. For PH, six associated markers were located on chromosomes 2H, 3H, 4H, 6H, and 7H. Only one marker located on chromosome 4H showed a positive effect on plant height, and increased height by 1.48 cm. Only one SNP, located on the short arm of chromosome 2H, was found to be associated with the variation in awn length with a negative effect. Three significant markers on chromosomes 1H, 2H, and 6H were detected for GPS. Those located on 1H and 6H reduced the grain number by 4.65 and 3.48 grains, respectively. However, SNPs located on 2H increased seed number by 8.55. Two significant SNP markers were located on 3H and 7H, controlling SL with a positive effect.

4. Discussion

4.1. Genetic Variation of the Studied Traits

The analysis of variance for each trait showed the existence of significant effects of the genotypes over the environment, indicating important genetic variability and differences among barley accessions evaluated under semi-arid conditions (Table 1). PCV and GCV values ranged from 6.6% to 43.2% and 6.3% to 42.6%, respectively. The lowest values were estimated for the heading date, and the highest was attributed to GPS. ECV variance was quite high for SL and low for other traits, suggesting that environmental effects may bring favorable changes in the expression of SL. These observations are in accordance with earlier results from [42,43]. Heritability values were medium to high, ranging from 0.67 for PH to 0.96 for GPS. Traits with a very high heritability indicate a low effect of environment on the phenotypic expression and the effectiveness of selection based on their phenotypic selection [43].

Assessing the genotypic responses in various environments is crucial for identifying suitable genotypes and evaluating environments that can enhance breeding and cultivar development [44,45,46]. The variance analysis conducted in this study revealed significant differences (p ≤ 0.001) between genotypes, test environments, and the impact of genotype–environment interaction. Analysis of variance for studied traits revealed that the genotypic effect was considerably more important than that of the environment, ranging from 60% to 90.3%, except for plant height, for which the environmental effect was much stronger (50.1%), suggesting that PH is modulated by a combination of genetic and non-genetic factors [47].

4.2. Marker Trait Association of the Studied Traits

Genome-wide association analysis has become a highly effective method for studying quantitative traits. It has the potential to offer more detailed mapping of QTLs compared to bi-parental studies [48]. In the current study, GWAS analysis was conducted using a panel of diverse barley accessions for several traits.

4.2.1. Heading Date

Heading date is an important trait that determines the crop’s ability to mature before the onset of unfavorable conditions such as heat stress, water stress, or disease [49]. This is particularly critical in regions with changing climates, as earlier flowering may help crops to avoid environmental stress and increase yields. In this study, 11 MTAs were associated with HD and distributed over all chromosomes except chromosome 4. The SNP (11_21140) on the long arm of chromosome 1H coincides with the position of the earliness per se gene eam8/eps1. In barley, the eam8 mutation causes an early-flowering-day phenotype with rapid flowering under both short days (SDs) and long days (LDs) [50]. The QTL on the short arm of chromosome 2H (12_30872) could correspond to the Ppd-H1 gene. Under LDs, the Ppd-H1 dominant allele causes early heading and is prevalent in Mediterranean conditions [20]. In this study, this gene caused a delay in the heading date, suggesting a mutation in the conserved CCT domain or the sixth exon of the Ppd-H1 gene [20,51]. In the long arm of chromosome 2H, 1 SNP (SCRI_RS_174077) exceeded the Bonferroni level. This region was previously identified in several studies as a region that included the HvAP2 gene [52,53,54], conferring a late heading date and co-segregating with the Flt-2L gene in wheat [55]. The detected locus in 3HS in this study is approximately 20 cM from the known heading date genes HvFT2 [56] and HvGI [57]. However, this region co-located with an associated SNP marker (SCRI_RS_189757) during the time from tipping to heading [58].

Another significant region is co-localized with the positions of the Vrn-H1 and HvPHYC genes on 5H. Both loci have a negative effect, causing an early heading date. It is known that deletion in the first exon of Vrn-H1 is associated with spring growth in barley [59]. Barley accessions carrying spring alleles of the Vrn-H1 gene showed early heading under drought conditions [60]. HvPHYC, a gene in barley that is similar to the PHYTOCHROME C locus, plays a key role in determining early flowering and regulates the process under long-day conditions (Nishida et al., 2013). Through its interactions with other genes such as Vrn-H1, sdw1, and Ppd-H1, HvPHYC can trigger early heading specifically in response to long daylight hours (Pankin et al., 2014).

On chromosome 6H, two associations were found. A locus on chromosome 6H is co-located with the flowering gene (HvNAM-1) [61]. The other locus is close to the putative heading and flowering genes: HvCO7 [17,62], HvCYP734A7 [17], HvCO5 [63], HvCry1a [64], and HvCry2 [65].

On the long arm of chromosome 7H, three MTAs were identified. Two SNPs (12_30026, SCRI_RS_204483) were associated with a delay in the heading date, and one (SCRI_RS_222863) had an early heading date. This region harbors several candidate genes, including HvCO6 and HvCMF7. In the short arm of 7H, one SNP has passed the FDR cutoff and is located near the vrn-H3 gene. According to our results, this locus promotes early heading, which coincides with previous research suggesting the central role of vrn-H3 in the flowering pathways, leading to the transition from vegetative to reproductive state under long-day conditions [28,66].

4.2.2. Plant Height

Six MTAs were identified as showing a negative effect on plant height, except the SNP (11_10639) located on 4H, which increases plant height. The genomic region on the short arm of 2H is co-localized distal to the HvCEN gene. The wild alleles of the HvCEN locus were found to reduce the phenotypic value of plant height [67]. The locus on the long arm of 2H (SCRI_RS_135355) coincides exactly with the position of HvAPO2/BFL. This locus was found, in previous studies, to control plant height [67]. BFL has been involved in phase duration and tillering by reducing the tiller number by one in barley [17,58]. On 3HL, one locus reduces plant height by 1.32 cm. This region harbors the uzu gene, the first dwarfing gene reporting an important mutation related to the brassinosteroid receptor [68], and also the site sdw1/denso gene controlling dwarfism by a mutation in the HvGA20ox2 gene [69]. On 4H, the candidate gene of the SNP (11_10639) locus is HORVU.MOREX.r3.4HG0394090. This locus was found 660 kb away from the location of HORVU.MOREX.r3.4HG0394170, the ortholog of OsDWARF4 in rice, and known as the HvD4 gene in barley [17]. The locus on 6HS was reported as a novel QTL lacking known candidate genes. On 7H, one locus showed a decrease in plant height, including HvMAX2 and Hv soluble starch synthase 1 (HvSS1).

4.2.3. Spike-Related Traits

Only one locus, located on 2HS, was identified for awn length, corresponding to the position of the previously identified AL2.1 QTL [70]. The GrainGenes database (https://wheat.pw.usda.gov accessed on 1 Jan 2022) lists the candidate genes in this region as HORVU.MOREX.r3.2HG0114290, which codes for a xylose isomerase, and HORVU.MOREX.r3.2HG0114320. Three SNP associations were identified for the GPS trait on 1H, 2H, and 6H. On 1H, the identified locus corresponds to the gene Vrs3, which encodes a histone demethylase that controls lateral spikelet development in barley [71,72]. On 2H, another locus controlling the development and fertility of lateral spikelets is called Vrs1 [73]. The SNP (SCRI_RS_170058) on 6H coincides with the position of the candidate gene HORVU.MOREX.r3.6HG0569710. For spike length, the putative candidate gene for the 3H locus is HORVU.MOREX.r3.3HG0297290, which encodes a protein kinase. The detected SNP on 7H is close to the HORVU.MOREX.r3.7HG0749940 gene, also known as HvBRD2/HvDIM, which is believed to affect plant height and tillering [17].

5. Conclusions

The results of GWAS analysis in barley under a semi-arid climate revealed important genetic markers associated with yield, adaptation, and important agronomic traits such as heading date and plant height. The identification of these markers provides a valuable resource for breeding programs aimed at developing new barley varieties that are better adapted to semi-arid conditions and have improved yield and quality characteristics. The use of these markers in breeding programs can accelerate the development of new varieties, reduce the cost of breeding, and ultimately contribute to sustainable agricultural production in semi-arid regions. In this study, several important loci were identified, especially for heading date and plant height. In barley breeding, the development of reduced plant height varieties with short growing cycles has been a significant objective. This has led to the production of cultivars that are suitable for high-density planting and are more resistant to lodging. These traits are particularly important during the grain-filling stage, when the plant requires adequate photosynthesis and nutrient uptake for optimal yield. However, heat stress can limit photosynthesis and reduce grain yield during this stage, making it challenging for growers to achieve high yields. To mitigate this, we will focus on developing varieties with enhanced heat tolerance and early heading, enabling the crop to flower and set seed before the onset of high temperatures. Through the development of such varieties, barley growers can increase their yields and produce more resilient crops in the face of changing climate conditions.