This study was based on the Ethiopian and Eritrean Barley Collection (EEBC), a germplasm assembled by Mamo et al. (2014) from accessions held at USDA-ARS National Small Grains Collection, the International Centre for Agricultural Research in the Dry Areas (ICARDA) and N.I. Vavilov All-Russian Institute of Plant Genetic Resources (VIR). The EEBC was previously used to identify novel genes for disease resistance and other agronomic traits including grain Zn (GZC) and Fe (GFeC) concentrations. All genotypes were originally phenotyped in 2011 glasshouse and 2012 field trials in Saint Paul, Minnesota (Mamo et al., 2014; Mamo and Steffenson, 2015). The datasets were downloaded from the online archive [T3/Barley (triticeaetoolbox.org)]. Before applying any genetic analysis, missing phenotypic values (5%) were imputed using the MICE (Multivariate Imputation via Chained Equations) package in R software [CRAN - Package mice (r-project.org)] with a maximum of 50 iterations and 5 replications. MICE assumes that the missing data are Missing at Random (MAR), which means that the probability that missing data are related to observed values and can be predicted using linear regression for continuous missing values and logistic regression for categorical missing values. MICE imputes data on a variable-by-variable basis by specifying an imputation model per variable. Data were then checked for a normal distribution and descriptive statistics were used to show the spread. Assumptions of ANOVA were subsequently checked while fitting the data and estimated Best Linear Unbiased Predictor (BLUP) and Best Linear Unbiased Estimate (BLUE) models were assessed when the assumptions were violated. Biological replicates of three seeds per plot were treated as environments to estimate BLUEs and BLUPs. The BLUEs are solutions (or estimates) associated with the fixed effects and BLUPs are the solutions (identified as predictions) associated with the random effects of a model. The BLUEs, BLUPs and correlations (genetic and phenotypic) between genotypes were analysed using META-R software (Alvarado et al., 2015; 2020). Descriptive statistics such as mean, minimum, maximum, and standard deviation were obtained using the ‘stat.desc’ function of ‘pastecs’ package (Grosjean et al., 2014) in the latest version of R software (R Core Team, 2020). Analysis of variance (ANOVA) was conducted using the ‘aov’ function, available in the stats package of R software. Normality and homogeneity of variance were checked by Shapiro–Wilk (Shapiro and Wilk, 1965) and Levene’s test (Levene, 1960), respectively.

Genetic analyses were performed on the DNA extracted from 7-10 days old leaf samples using QIAamp 96 DNA kit (Qiagen GmbH) on an automated nucleic acid purification robot. DNA quality was assessed using a Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, United States) with a requirement for 260/280 and 260/230 ratios > 1.8. DNA was then quantified using Picogreen (Thermo Fisher Scientific, Waltham, MA, United States). 300ng lyophilised DNA per sample was sent to Geneseek (Neogen Europe, Ltd., Auchincruive, United Kingdom) for genotyping using a high throughput Barley 50K SNP array (Bayer et al., 2017) by Illumina HTS processing and HiScan array imaging (Illumina, San Diego, CA, United States). R and Theta scores were extracted from resulting idat files using Genome Studio Genotyping Module v2.0.2 (Illumina, San Diego, CA, United States) and allele scores were created using paRsnps (an R package for clustering, visualising, and comparing Illumina SNP genotyping data). The data were visualised with, and exported from, Flapjack, a genotype visualisation software for further analysis (Milne et al., 2010). The genotypic data were assessed for missing, heterozygous, and monomorphic markers and markers with more than 10% missing data were removed. Polymorphic markers were maintained for cluster analysis, population structure and GWAS. For cluster analysis, the R package ‘SelectionTools’ [Index of/~software/(uni-giessen.de)] was used to read in marker data to construct Principal Components Analysis (PCA) for population structure. Euclidean distance was used to differentiate between genotypes on the PCA axes based on their genetic distance. The optimum number of clusters (K) was determined using the Elbow method which is inferred at the highest value beyond which the total within sum of squares (WSS) is no longer significantly reduced, and the gap statistic (GS) which is increasing at decreasing rate. Population structure analysis was also performed using STRUCTURE software (Pritchard et al., 2000) with a hypothetical number of sub-groups (K) ranging from 1-10 and 5000 burn-in iterations followed by Markov Chain Monte Carlo, MCMC) iterations of 5000. Five independent runs were performed for each value of K. Molecular variances (within and among population variances) were analysed in GenAlex software (Peakall and Smouse, 2012) using 10% of the markers (5,000) randomly selected from all chromosomes. The fixation index (Fst) among the EEBC panel, also known as Wright’s population differentiation statistics (Wright, 1951) was analysed in R software.

GWAS was initially run to compare with published results on grain colour and the QTL interval was set based on the extent of linkage disequilibrium (LD). Genome-wide and within chromosome LD (r²) decay were analysed in Trait Analysis by aSSociation, Evolution, and Linkage (TASSEL) software (Bradbury et al., 2007) with a sliding window of 500 markers around the significant marker to obtain the extent of inter-chromosomal LD decay. All possible candidate genes were identified based on an LD score of r² ≥ 0.2 between the most significant marker associated with the trait and neighbouring SNPs. The most likely candidate genes from the list were highlighted according to their predicted functions. GWAS was performed using the grain Zn concentration data from Minnesota as described in 2.1 above and 36,431 polymorphic SNPs, and the result compared to the GWAS in the study by Mamo et al. (2014) that used 7,842 polymorphic markers. Results from several GWAS (single locus, SL_GWAS versus multi-locus, ML-GWAS) models were compared to validate the observed results. A Q-Q plot produced by the best fit model was used to determine whether the GWAS model was robust and detected only the most significantly associated SNP markers with the highest observed versus expected [-log₁₀P] values. Too many markers below the Q-Q proportional line means the model is overfitting the data with false negative results, while those above the line indicate false positive results. The first three principal components (PCs) were added as covariates in the model to reasonably capture and correct for the confounding effect of population structure. GWAS was performed with the GAPIT3 package in R software (Wang and Zhang, 2021) using the BLUE data for the phenotypic trait. The results were compared between the multi-locus BLINK and the single-locus MLM models. QTL intervals were determined based on significant LD of squared correlation value r² ≥ 0.20 and possible candidates were prioritised based on genes annotated with relevant functions in other species.

Genomic sequences of the HvZIFL1 candidate gene identified in the genome-wide association (GWA) analyses were downloaded from the genome assembly of Morex v1 (Mascher et al., 2017)(HORVU4Hr1G081570.1). Both version 2 (Monat et al., 2019) (HORVU.MOREX.r2.4HG0341120.1) and version 3 (Mascher et al., 2021) (HORVU.MOREX.r3.4H409580) were released during this study and used for subsequent genetic characterisation of the gene. The latest reference transcriptome database [Morex Gene Atlas - Hv_Mx_chr4HG27982 (hutton.ac.uk)] was used to identify potential transcripts and compared to EnsemblPlants which maintains 5 transcripts of HvZIFL1 (https://plants.ensembl.org/Hordeum_vulgare/Gene/Summary?g=HORVU.MOREX.r3.4HG0409580;r=4H:586550476-586579926;t=HORVU.MOREX.r3.4HG0409580.1;db=core). A Kompetitive Allele Specific PCR (KASP) assay was designed ( Supplementary Table 6 ) to the most significantly associated marker for technology independent validation across the EEBC germplasm as well as for other barley germplasm. Haplotype analysis was performed for 194 georeferenced EEBC landrace collections using 19 SNPs from the 50K SNP platform and a locus haplotype was also plotted using 80 polymorphic SNPs across the gene obtained from exome capture data of the European Wheat and Barley Legacy for Breeding Improvement (WHEALBI) project (Russell et al., 2016; Bustos-Korts et al., 2019).

Amino acid sequences of ZIFL1 paralogues were downloaded from EnsemblPlants and BLAST searched in the National Centre of Biotechnology institute (NCBI) database. The top 100 sequences were downloaded and used for multiple alignment in Jalview 2.11.2.1 software (Waterhouse et al., 2009) and sequences were uploaded in the Multiple Expectation maximisation for Motif Elicitation (MEME) online database [ https://memesuite.org/meme/tools/meme] to explore protein motifs and hence possible similar functions (Bailey and Elkan, 1994). Known conserved motifs and new signatures were discovered using detailed analysis in the latest version of Jalview (Waterhouse et al., 2009) and by multiple alignment using Muscle (Edgar, 2004).

Population structure is one of the confounding factors limiting efficiency of detecting quantitative trait loci in genome-wide association studies. It is hence crucial to characterise and determine the diversity and genetic structure existing within the germplasm panel. To this end, we used the same extensive EEBC panel assembled by University of Minnesota (Mamo et al., 2014) to determine the level of population structure. We genotyped the accessions using the barley 50K SNP assay (Bayer et al., 2017) and after filtering to remove failed, missing (>10% removed), monomorphic and low frequency (< 5%) markers, 36431 polymorphic gene-based and physically mapped SNPs remained. We observed four major clusters of the EEBC germplasm with 72% present in subpopulation 1 (SP1) with no partitioning by row type or by geography ( Figures 1A, C ). Molecular variance explained 69% and 31% ( Figure 1B ) of genetic diversity within and among subpopulations of the EEBC, and the first two principal components (PC1 and PC2) explained 82% and 7% of the genetic variation between the sub-populations, respectively. To examine diversity of the EEBC in a global context 50K SNP data was combined with the 1000 accession core barley collection in the German Federal Gene bank IPK Gatersleben which had been genotyped using the same platform (Darrier et al., 2019). Principal Component Analysis (PCA) uncovered the level of diversity and population structure with the first (PC1) and second (PC2) axes explaining 41% and 25% of the variation, respectively identifying four separate genotypic clusters ( Supplementary Figure 1A ). The Ethiopian (EEBC) germplasm was shown to form tight clusters with 238 lines (80%) of the 296 accessions forming a single main cluster (highlighted in green in Supplementary Figure 1A ). Although 46 Ethiopian lines were among the core panel, only half grouped within the EEBC cluster. Molecular variance depicted 57% of the variation within and 43% between sub-populations ( Supplementary Figure 1B ).

The phenotypic data (Mamo et al., 2014) was analysed and descriptive statistics used to summarise grain Zn concentration of the EEBC panel in both Minnesota glasshouse and field trials ( Supplementary Table 1 ). Grain Zn concentration ranged between 6-72 mg/kg in the glasshouse and 20-87 mg/kg in the field. Analysis of variance (ANOVA) depicted significant difference between genotypes (G), locations/environment (E) and genotype-by-location/environment interaction [GEI] ( Supplementary Table 2 ). Genotype, location/environment (E) and genotype-by-environment interaction (GEI) contributed 32%, 36%, and 21% to the total variance (% total sum of squares, %TSS) ( Supplementary Table 2 ). In contrast, within-population analysis revealed no statistical difference among the sub-populations ( Figure 1D ) in the 2011 glasshouse trial with a mean GZC ranging from 28 mg/kg for SP4 to 32 mg/kg for SP2. The 2012 field trial, however, depicted significantly higher GZC of 51 mg/kg for the SP1 compared to the lowest grain Zn concentration of 41 mg/kg for SP4. The two sub-populations SP2 and SP3 had statistically similar grain Zn concentration of 44 mg/kg and 45 mg/kg, respectively. The SP4 is comprised of 20 two-rowed elite cultivars or breeding lines including the UK standard check variety known as Concerto.

We then performed GWAS, accounting for the observed population structure, and compared several single locus (SL-GWAS) and multi-locus (ML-GWAS) models with the objective of identifying consistent associations regardless of the statistical model. GWAS for grain zinc concentration (BLUE data, Mamo et al., 2014) from glasshouse and field trial was analysed using three different statistical models ( Figures 2A, B ; Supplementary Table 4 ; Supplementary Figures 2A, B ). For the glasshouse data GWAS consistently revealed a significant association at the bottom of chromosome 4HL. From the field data, only BLINK detected a significant association at this locus ( Supplementary Figure 2B ). Given grain Zn concentration is likely to be under complex genetic control, on weight of evidence we decided to investigate the terminal 4H region in more detail. High SNP marker density allowed the association boundaries to be defined accurately based on the drop in LD ( Supplementary Figure 3 ) between the most highly associated and physically adjacent non-significant markers. The most highly associated SNP (JHI-Hv50-2016-265280) was located within a gene annotated as a member of the Major Facilitator Superfamily (MFS) of transporters that includes Zinc Induced Facilitator (ZIF) and ZIF-Like (ZIFL) proteins. As no other annotated genes within the region appeared strong candidates, we prioritised HvZIFL1 as a strong causal gene candidate for more detailed analyses.

We extracted the complete HvZIFL1 (HORVU4Hr1G081570.1) gene sequence from Morex reference genome assembly version 1 (Mascher et al., 2017) to investigate molecular polymorphism, and again when version 2 became available (Monat et al., 2019). The protein encoded by HvZIFL1 consists of 493 amino acid residues and its secondary structure was predicted to contain a typical MFS protein structure with 12 membrane spanning domains known as transmembrane helices (TMH) that facilitate substrate transport across the biological membrane ( Supplementary Figure 4 ). We predicted that a mutation changing Glutamine to Histidine (Q330H) and an adjacent Valine to Phenylalanine (V331F) respectively) located at the gate of TMH8 ( Figure 3B ; Supplementary Figure 4 ) may be functionally disruptive (positive mutation because it increase grain GZC). Further examination of the 50K SNP data revealed that 19 polymorphic SNPs were located within HvZIFL1 and haplotype analysis revealed eight haplotypes ( Table 1 ) across the 297 (296 EEBC plus Concerto) accessions with frequencies ranging from 0.007 (2 accessions) to 0.581 (172 accessions). In the combined EEBC and IPK dataset (Darrier et al., 2019), we observed 10 unique haplotypes with frequencies ranging from 0.003 to 0.536. The most significantly associated GWAS SNP was a T to G transversion that was unique to Haplotype 1 ( Tables 1 , 2 ) which had significantly higher grain zinc concentration (47.6 mg/kg) when compared to all other haplotypes (42.8 mg/kg for haplotype 5 to 34.1 mg/kg for haplotype 4).

We designed 12 pairs of primers ( Supplementary Table 5 ) to amplify the complete HvZIFL1 based on Morex version 1 and PCR-sequenced the gene in twelve genotypes that represented high and low grain Zn concentration. A total of 95 SNP polymorphisms were identified, almost all of which are silent/synonymous changes occurring within non-coding regions (black inverted triangle marks, Figure 3A ). The two exceptions were the diagnostic SNP marker JHI-Hv50k-2016-265280 identified by GWAS, and an adjacent SNP discovered by Sanger sequencing using primer pair 9 ( Supplementary Table 5 ; indicated with branching lines in Figure 3 ). Sequencing and alignment also revealed a variable INDEL of either 162 bp or 153 bp in length within the non-coding region between exons six and seven ( Figure 3B ). This INDEL was confirmed using agarose gel electrophoresis and was present in 10 of the 12 Ethiopian barley sequences, with EEBC_005 having a 162 bp INDEL and EEBC_114 the same 153bp INDEL as the Morex reference annotation ( Figure 3B ). Assaying the INDEL across the EEBC identified 51 lines with the 153bp insertion identical to the reference cultivar Morex ( Supplementary Table 8 ). The highly significant markers from GWAS were found in Exon 11 changing a CAA (Glutamine, Q) to CAC (Histidine, H) with the adjacent mutation changing GTT (Valine, V) to TTT (Phenylalanine, F). Somewhat surprisingly, the most significantly associated GWAS SNP was not in complete LD with the 153 bp deletion. Further investigation of the HvZIFL1 gene sequence revealed duplication of ~7 kilobase pairs with 15.5 kbp interspaced sequence ( Figure 4 ).

Because the 50K SNP assay has been widely used by the barley research community as a platform for genetic analysis and marker assisted selection (MAS), we were able to examine the distribution of the highly associated marker JHI-Hv50k-2016-265280 across a wide range of germplasm including 575 cultivars, 197 landraces and 258 wild barley accessions ( Table 2 ). The frequency of the minor allele associated with increased grain zinc is 0.202 in the EEBC compared to 0.047 in cultivars and 0.004 in wild barley accessions. These frequencies observed in both the landrace and wild germplasm are below 0.05 which is usually chosen as the threshold for robust marker:trait detection in association studies and would therefore not have been significant in any of these germplasm sets. We also assembled haplotype sequences from the genetically diverse WHEALBI collection (N= 432) (Bustos-Korts et al., 2019) using 80 polymorphic SNPs discovered in exome capture data. We identified 38 unique haplotypes where the favourable JHI-Hv50k-2016-265280 ‘C’ allele was found only in Haplotype 18 (Hap_18), comprised of 40 lines of which 30% were from Ethiopia and Eritrea ( Supplementary Tables 7 , 9 , 10 ). Analysis of grain Zn concentration showed a significant effect associated with both the 153bp INDEL and the diagnostic marker, individually and in combination ( Figures 5A–E ). The boxplots and t-tests showed that genotypes with the 153 bp deletion in ZIFL1 had a mean grain Zn concentration of 31 mg/kg compared to 28 mg/kg for those without the deletion (t=2.71**, Figure 5A ) using the Minnesota (glasshouse) data. The difference between these groups in the field was 8 mg/kg (t= 5.23**, Figure 5B ). The JHI-Hv50k-2016-265280 SNP effect for the glasshouse and the field trials showed an 8 mg/kg (t= 5.94***) and 12 mg/kg (t= 7.55***) difference between alleles ( Figure 5C ). Variability in grain Zn concentration was also apparent in the data across the two locations/years ( Figures 5D, E ). In the 2011 glasshouse trial, deletion had an average effect of 1 mg/kg (t-statistics = 1.12 ns) while the SNP marker had a mean effect of 8 mg/kg (t=5.67***, Figure 5D ). In the 2012 field trial, deletion had an average effect of 5 mg/kg (t= 3.50***) while the SNP effect was 11 mg/kg (t= 6.53***). Interestingly, the combined effect of the 153bp INDEL and the favourable allele increased grain Zn concentration level by 9 and 16 mg/kg in the glasshouse and field conditions respectively in Minnesota.

A BLAST search of the amino acid sequence of the candidate HvZIFL1 gene identified 33 paralogues distributed across the barley genome ( Supplementary Table 3 ). For most of these genes, SNPs were present on the 50K chip and almost all were polymorphic with minimum allele frequency (MAF) values ranging from 0.021 to 0.476. However, none of the SNPs were associated with grain zinc concentration except the HvZIFL1 on 4H with an MAF of 0.202. We then compared the HvZIFL1 protein sequence with its cereal orthologs and identified 30 short conserved sequence motifs that may be related to both conserved and unique functions of the protein. The motifs could be categorised into three major and three minor groups ( Supplementary Figure 5 ). The first major group (G-I) contained HvZIFL1 and 10 other ZIFL genes and exhibited 11 conserved motifs. The second major group (G-II) consisted of 9 ZIFL proteins with 10 unique motifs while the third group (G-III) contained 9 ZIFL genes with 12 unique sequence motifs including motif 6 and 15 which were conserved in all groups. The three minor groups consisted of 1-2 ZIFL genes with four conserved motifs.

We used the cv. Morex EoRNA gene expression atlas to explore the temporal and spatial transcript abundance and post transcriptional variation of HvZIFL1 (https://ics.hutton.ac.uk/morexgeneatlas/index.html). HvZIFL1 (transcript Hv_Mx_chr4HG27982 in MorexGeneAtlas) (Guo et al., 2025) was detected as multiple alternative transcripts ( Supplementary Figure 6 ) with all, including the most abundant, expressed exclusively in root tissues with no indication of expression in any above ground tissues, including the developing grain ( Supplementary Figure 7 ). These observations are consistent with HvZIFL1 contributing to variation in total Zn uptake from soil. Whether subsequent transport to the developing grain is passive (e.g. via the transpiration stream) or active (involving other Zn transporters) could not be determined but does not apparently involve HvZIFL1. Furthermore, there was little variation in the level of HvZIFL1 expression in roots when Zn (or other heavy metals) were added to the growth medium suggesting its expression is not responsive to Zn availability (Kintlová et al., 2017) but rather is genetically pre-determined in different genotypes.

The diagnostic SNP marker (JHI-Hv50k-2016-265280) identified by GWAS was selected for validation as a tool for marker-assisted selection (MAS) using a KASP assay. Alternative alleles in all EEBC accessions were then classified using the KASP assay. 237 lines had the JHI-Hv50k-2016-265280 ‘A’ allele associated with lower grain zinc concentration and 59 lines the JHI-Hv50k-2016-265280 ‘C’ allele associated with higher grain zinc concentration ( Figures 5A, B ). Comparing the data from the KASP assay and the 153bp deletion PCR-based size marker, the JHI-Hv50k-2016-265280 ‘C’ allele was only associated with the deletion whereas the ‘A’ allele was associated with and without the deletion in a frequency of 0.212 (‘A’ allele, no deletion) and 0.787 (‘A’ allele, with deletion). The impact of using both markers increased the differences of grain zinc concentration to 11 mg/kg and 16 mg/kg for glasshouse and field data, respectively ( Figures 5C, D ).