Auxin response factors (ARFs) are a family of plant-specific transcription factors that mediate auxin signaling by binding to auxin-responsive elements in gene promoters, thereby regulating a broad array of developmental and physiological processes. These include organogenesis, embryogenesis, apical dominance, root and shoot development, and reproductive differentiation (Cancé et al. 2022). Typically, ARF proteins harbor three functional domains: a conserved N-terminal B3 DNA-binding domain (DBD), a variable middle region that acts as either a transcriptional activator (glutamine-rich, Q-rich) or repressor (serine/threonine/proline-rich, STP-rich), and a C-terminal PB1 (Phox and Bem1) domain that facilitates homo- and hetero-dimerization with other ARFs or Aux/IAA proteins (Chandler 2016). This domain architecture is remarkable conserved from early land plants to higher angiosperms, underscoring the fundamental roles ARFs play across plant evolution (Li et al. 2023).

In numerous species—such as Arabidopsis thaliana, rice (Oryza sativa L.), maize (Zea mays L.), and wheat (Triticum aestivum L.)—ARF genes have been systematically identified and functionally characterized (Gidhi et al. 2021; Man et al. 2025; Okushima et al. 2005; Wang et al. 2007; Xing et al. 2011). These studies revealed the involvement of specific ARF members in regulating key agronomic traits such as grain size, inflorescence architecture, and stress responses. For example, OsARF3, OsARF4, OsARF6, ZmARF24, ZmARF2, ZmARF17, and TaARF25 have been reported to regulate grain/kernel traits (Gao et al. 2023; Gu et al. 2023; Hu et al. 2018; Jia et al. 2021; Qiao et al. 2021; Wang et al. 2024). Notably, OsARF3 is also involved in balancing the trade-off between abiotic and biotic stress responses. OsARF12, OsARF17, OsARF25, and OsARF10 have been shown to affect stigma development (Guo et al. 2022; Zhao et al. 2022). ZmARF28 can pleiotropically modify plant architecture by single amino acid substitution, its mutant shows reduced plant height, increased leaf number and changed inflorescence identity (Prigge et al. 2025). These functional insights demonstrate the evolutionary versatility and their potential for agronomic improvement of ARF genes across cereals. Despite its agronomic and historical significance, as one of the first domesticated crops (Lister et al. 2018), barley (Hordeum vulgare L.) has not yet seen comprehensive systematic characterization of its ARF gene family. Although ARF genes have been implicated in regulating key yield and fitness traits in related cereals, their characterization in barley remains limited, highlighting a significant knowledge gap. Given the known roles of ARF genes in other species, elucidating the diversity and function of the ARF protein family in barley could reveal previously unknown regulators of important developmental processes and contribute directly to crop improvement efforts.

Recent advances in high-throughput sequencing and the development of pan-genome and pan-transcriptome resources have transformed plant genomics by enabling species-wide resolution of genetic diversity (Marks et al. 2021). Traditional single-reference-based analyses are insufficient to capture the full extent of gene family variation, as they often miss presence/absence variants (PAVs), copy number variations (CNVs), or population-specific alleles. Pan-genomic resources—such as those recently generated for barley—offer a more comprehensive view of structural and functional diversity at both the gene and genome levels (Jayakodi et al. 2024; Jayakodi et al. 2020). In parallel, transcriptome datasets across tissues, genotypes, and even at single-cell resolution now make it possible to explore gene expression divergence, regulation, and potential adaptation at an unprecedented scale.

In this study, we performed a comprehensive, population-scale analysis of the ARF protein family in barley using integrated multi-omics data. We systematically identified and annotated genes and proteins of the ARF family across the 76 barley pan-genomes, revealing approximately 30% more members than in previous studies. By constructing high-resolution phylogenies and characterizing structural variation, expression patterns, cis-regulatory elements, and selective signals, we explored the evolutionary and functional dynamics of ARFs in barley. Our findings revealed significant differences in natural selection pressures and genomic diversity across distinct barley populations. Notably, we identified a superior allele of HvARF3 that appears to have undergone selection in European cultivars, with significant association to increased grain size and weight. These findings demonstrate how integrative pan-genome and transcriptome analyses can illuminate both conserved gene family evolution and novel, population-specific adaptations with agronomic relevance. Our results provide a comprehensive panoramic view of the ARF protein family and establish a foundation for future studies on auxin signaling, inflorescence development, and molecular breeding in barley.

Raw transcriptome data from five tissues—caryopsis, coleoptile, inflorescence, root, and shoot—across 20 barley cultivars were retrieved from the recently published barley pan-transcriptome dataset (accession: PRJEB64639). Transcript abundance was quantified using Kallisto v0.46.1 (Bray et al. 2016), yielding TPM values and raw read counts. To minimize genotype-specific effects in the functional analysis of ARF genes, genotype-corrected residual expression values ( ) were calculated (see below) and used for downstream analyses. Weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA R package with default parameters and a correlation threshold of 0.7. Gene Ontology (GO) enrichment was conducted using clusterprofiler v3.2.1 (Yu et al. 2012), with GO annotations sourced from the IPK barley database (https://galaxy-web.ipk-gatersleben.de/). Single-cell expression patterns of ARF genes were extracted from a recently published barley inflorescence single-cell transcriptomic resource BARVISTA (https://www.plabipd.de/projects/hannah_demo/BARVISTA/tissue_no_download.html) (Demesa-Arevalo et al. 2025).

To quantify the contribution of genotype and tissue to gene expression variation, a linear mixed model (LMM) was fitted for each ARF gene using the following equation:

Where represents the expression value (TPM), is the random genotype effect, is the tissue effect, and is the residual error. The proportion of variance explained by genotype and tissue was calculated as:

Genotype- and tissue-corrected residuals were then extracted using:

All calculations were performed in R using lme4 package (Bates et al. 2015).

TE identification and cis-eQTL analysis

TEs were identified in 76 barley genomes using panHiTE (Hu et al. 2025), and further confirmed using TE annotations from the IPK panBarlex database (https://panbarlex.ipk-gatersleben.de/#). Only high confidence TEs present in both of datasets were retained. To determine TE insertions near ARF genes, 100 kb flanking regions were analyzed using bedtools. For cis-eQTL analysis, 1 Mb upstream and downstream of each ARF gene was defined as the cis-regulatory region. values were used as gene expression input in MatrixEQTL v2.3 (Shabalin 2012) under default settings.

To explore potential domestication-related selection signals, raw genomic data were obtained from a recent barley population study (Guo et al. 2024). A 1 Mb window surrounding each ARF gene was analyzed. Genomic regions were extracted using bcftools v1.15.1, and downstream processing and visualization were carried out using Perl under the guidance of the original study’s first author (Guo et al. 2024).

Whole-genome assemblies of the 76 barley accessions were re-mapped to the Morex V2 reference genome (Figure, S4). Genomes were split into 1 kb fragments using seqkit v0.9.1 (Shen et al. 2016) to simulate FASTQ reads. Variant calling was performed using BWA v0.7.17 (Li 2013), SAMtools v1.16.1(Li et al. 2009), and BCFtools v1.15.1. SNPs and InDels were annotated using SnpEff (Cingolani et al. 2012). Amino acid substitutions in ARF proteins were predicted with PPVED v1.0 (Gou et al. 2022). Correlation analysis between substitution scores and domestication degrees were calculated using Pearson’s correlation coefficient in R.

To genotype the functional HvARF3 R475G variant, a dCAPS marker was developed targeting a SNP located at chr3H:494973445 (based on the MorexV2 reference genome). Primers (Supplementary data 5) were designed using the online tools dCAPS Finder 2.0 and Primer3, incorporating a mismatch to generate a specific restriction site for the enzyme DdeI (New England Biolabs Inc). PCR amplification was carried out using Taq DNA Polymerase (Qiagen, Cat. No. 201203) following the manufacturer’s recommended protocol. Amplified products were subsequently digested with DdeI and separated on a 2.5% agarose gel.

Identification and validation of ARF genes and their PAVs/CNVs across 76 barley pangenome assemblies

To comprehensively catalog ARF-coding genes in barley, we employed a domain-based gene discovery approach by integrating three major protein family databases—Pfam (PF02362, PF06507, PF02309), ProSiteProfiles (PS50863, PS51745), and PANTHER (PTHR31384) —using the recently released barley pan-genome v2 database (BPGv2) (Jayakodi et al. 2024). All candidate sequences were further validated against the PanBARLEX online resource (https://panbarlex.ipk-gatersleben.de/#) to ensure annotation accuracy. Given that the conventional reference genome (Morex) lacks two ARFs due to genotype-specific presence/absence variations (PAVs), we selected the European modern spring barley cultivar Barke for gene nomenclature. Barke harbors all identified ARF-coding genes and has been extensively characterized at the transcriptomic level (Coulter et al. 2022; Guo et al. 2025). To maintain cross-species consistency, gene names were assigned based on the wheat gene nomenclature system (Boden et al. 2023) and aligned with orthologous ARF genes in rice (Wang et al. 2007).

Compared to previous single-genome based analyses (Tombuloglu 2019), our pan-genomic approach uncovered approximately 30% more ARF genes, underscoring the power of multi-genotype analysis to capture greater gene diversity. We retained 26 full-length ARFs containing at least one conserved ARF-related domain in ≥ 10% of accessions for downstream analysis. Based on their frequency across the 76 barley genomes, these genes were classified into core, soft-core, and shell categories following established dispensability thresholds (Marroni et al. 2014) (Fig. 1A, E). Notably, approximately ~ 80.7% of ARFs belonged to the core or soft-core categories, being present in 76 or 75 genomes, respectively. This high evolutionary conservation highlights the essential roles of ARFs for barley development and fitness.

An exception was the HvARF13 clade, which exhibited extensive structural and presence variation. While HvARF13-3 was highly conserved (found in 74 genotypes), HvARF13-1 and HvARF13-2 were present in only 47 and 59 accessions, respectively. Although HvARF13-1 and HvARF13-2 are shorter than HvARF13-3, all three retain an intact B3 DNA-binding domain in at least some genotypes and are expressed at transcriptome level, suggesting potential functionality (Fig. 1B). This structural variation reflects dynamic gene duplication events and copy number variations (CNVs) within the barley ARF family. To validate these PAVs, we extracted 5–10 kb flanking genomic sequences and performed local alignment. A 479 bp deletion was identified at the 3′ end of HvARF13-2 in several genotypes, including Morex (Figure S1B), likely explaining its absence in previous annotations. A larger deletion was detected for HvARF13-1 in some accessions, leading to complete gene loss (Figure S1A). Additionally, a truncated gene (HORVU.BARKE.PROJ.7HG00820280) lacking conserved domains was likely misannotated and was excluded after sequence and domain verification.

To explore functional implications of ARF protein diversity, we examined the composition of the middle region in each Barke ARF protein, combining domain annotation with phylogenetic analysis using rice homologs (Fig. 1C). We specifically quantified glutamine (Q) and serine/threonine/proline (STP) residue enrichment, which are indicative of functional differentiation (Chandler 2016; Shen et al. 2010). ARFs with Q-rich middle regions were classified as Class I (typically transcriptional activators; e.g. HvARF19), whereas those enriched in STP residues were categorized as Class II or Class III (transcriptional repressors; e.g. HvARF1) (Fig. 1D). In total, 26 well-characterized ARFs with high-confidence annotations were selected for our subsequent analyses.

Phylogenetic and Ka/Ks analyses reveal evolutionary dynamics and selection pressures on ARF proteins in barley

To elucidate the evolutionary relationships within the ARF gene family in barley, we performed a comprehensive phylogenetic analysis based on full-length ARF protein sequences derived from 76 barley genotypes. Despite extensive sequence variation among accessions, all ARF proteins could be reliably assigned to their respective clades, consistent with ARF classifications established in other plant species (Wang et al. 2007) (Fig. 2A). Comparative genome analyses across Arabidopsis, rice, and barley revealed substantial changes in ARF copy number, with both gene losses and expansions observed. Notably, several Arabidopsis ARFs, including AtARF9, AtARF13, AtARF14, AtARF23, AtARF12, AtARF22, AtARF15, AtARF10, AtARF20, and AtARF21, were completely absent in both rice and barley (Fig. 2A), suggesting lineage-specific gene loss events following the divergence of grasses and dicots (Aravind et al. 2000; Finet et al. 2013). This pattern likely reflects evolutionary streamlining of the ARF repertoires in grasses and possibly the adoption of distinct auxin signaling strategies (McSteen 2010).

Interestingly, the AtARF5 clade retained a single-copy ortholog across all three species tested here, indicative of strong purifying selection and a critical, possibly irreplaceable, role in plant development. In contrast, barley showed expanded ARF families compared to rice, particularly in HvARF4, HvARF8, HvARF10, and HvARF13. Among these, HvARF4-1 and HvARF4-2 represent a rare example of full-length gene duplication, with both paralogs retaining all key functional domains, including the B3 DNA-binding domain and PB1-like domain (Fig. 2A, 2B). Other duplicated genes, such as members of the HvARF13 clade, often exhibited partial structures lacking one or more conserved domains, suggesting potential subfunctionalization or neofunctionalization. Notable, most of the observed barley-specific ARF duplications clustered within the AtARF17/16/10-like clade, which includes AtARF10, AtARF16, AtARF17, and OsARF18—a group known to function as negative regulators of plant development and seed germination (Huang et al. 2016; Liu et al. 2007; Liu et al. 2013; Mallory et al. 2005). The increasing copy number of this clade from Arabidopsis to rice and barley suggests an adaptive expansion, potentially driven by environmental selection pressures or functional diversification. To investigate the variation among ARF protein in barley, all ARF protein sequences were analyzed and are presented in a schematic diagram (Fig. 2B) (Cancé et al. 2022). The DBD domain is relatively conserved in class I ARFs, whereas almost all class III ARFs contain a short insertion within this region, suggesting greater structural stability and conservation in class I. Notably, HvARF25, although belonging to class I, has a relatively high content of glutamine (Q) residues but does not form a poly-Q structure. Most class III ARFs possess the B3 repression domain (BRD, motif pattern: R/KLFG) (Choi et al. 2018), a characteristic indicator of repressor function. However, HvARF4-1 and HvARF4-2 are exceptions: instead of the canonical R/KLFG motif, each contain a single amino acid substitution that produces a modified KIFG motif, potentially indicating a functional divergence.

To assess selection constraints acting on ARF proteins, we calculated the ratio of non-synonymous (Ka) to synonymous (Ks) substitution rates (Ka/Ks) across different barley populations. Ka/Ks values below, equal to, or above 1 are generally interpreted as indicative of purifying, neutral, or positive selections, respectively (Kimura 1985). Overall, most ARF proteins exhibited Ka/Ks ratios well below 1, indicating strong purifying selection acting to conserve gene function across genotypes (Fig. 2C). This evolutionary constraint underscores the essential regulatory roles of ARF proteins during barley development and adaptation. However, HvARF8-2 was a notable exception, with Ka/Ks ratios frequently exceeding 1. The lack of a B3 domain in this gene may have relaxed functional constraints, permitting rapid evolution under natural or selective breeding. This gene may represent an example of lineage-specific divergence or a recently neofunctionalized paralog. Further comparison of Ka/Ks ratios among wild barley, landraces, and modern cultivars revealed a statistically significant decline from wild accessions to cultivars, indicative of intensified purifying selection during domestication and breeding. A similar trend was observed across geographic regions: accessions from East Asia exhibited lower Ka/Ks values than those from European and West Asia (Fig. 2D), suggesting differential selection pressures linked to regional adaptation.

Interestingly, HvARF12 deviated from this general pattern, showing elevated Ka/Ks ratios in cultivars compared to wild or landrace accessions. Its rice ortholog, OsARF12, is functionally redundant with OsARF17 and OsARF25, and has been implicated in auxin-mediated regulation of tiller angle, as well as in responses to viral infection and disease resistance (Li et al. 2020) (Zhao et al. 2022). These findings raise the possibility that HvARF12 has undergone adaptive evolution, possibly driven by trade-offs between agronomic performance and stress tolerance.

Taken together, these results reveal diverse evolutionary dynamics within the ARF family in barley, shaped by natural selection, domestication, and breeding. The combination of phylogenetic, structural, and Ka/Ks analyses provides key insights into functional divergence and adaptation of ARF genes in cereal crops.

TE identification, classification, and their regulatory influence on ARF gene expression

Transposable elements (TEs) are mobile genetic elements capable of within the genome, frequently inducing mutations or altering gene expression. Representing a substantial fraction of plant genomes, TEs have been recognized as major contributors to genome evolution, structural variation, and transcriptional regulation (Bennetzen and Wang 2014; Hirsch and Springer 2017; Lisch 2013). To explore the potential cis-regulatory effects of TEs on ARF gene expression, we surveyed TE insertions within 100 kb upstream and downstream of all ARF genes across 76 barley accessions, using a high-confidence TE annotation dataset.

In total, 516 TEs were identified near ARF loci, encompassing diverse TE families, with retrotransposons constituting the predominant class (Fig. 3A). Interestingly, TE insertions were not randomly distributed: significantly enrichment was observed within the 50 kb upstream of many HvARF genes. A secondary enrichment peak occurred around 20 kb upstream of several loci (Fig. 3B). Given their proximity to core promoter regions and their established roles in modulating chromatin accessibility and transcription factor binding (Hirsch and Springer 2017; Lönnig and Saedler 2002), these upstream TEs are strong candidates for modulating ARF gene expression via cis-regulatory mechanisms. The observed spatial enrichment pattern suggests non-random TE insertion or retention, possibly reflecting selection for regulatory function. To investigate the potential regulatory impact of nearby TEs, we conducted expression quantitative trait loci (eQTL) mapping (Druka et al. 2010) in a panel of 20 representative genotypes by integrating high-quality pan-genome and pan-transcriptome datasets (Guo et al. 2025). To disentangle the complex sources of expression variation, we partitioned the variance in gene expression into components attributable to genotype, tissue, and their interaction. Among the 26 ARF genes, 22 displayed sufficient expression variability and were retained for further analysis.

Variance component analysis revealed that tissue and tissue × genotype interactions were the primary drivers of expression variation, consistent with the known tissue-specific expression profiles of ARF genes (Fig. 3C). To mitigate tissue effects in subsequent association analyses, we used residual expression values—obtained by regressing out tissue effects—instead of raw expression levels. Cis-eQTLs were identified using SNPs located within a 1 Mb window flanking each ARF gene. As anticipated, genes with high genotype or genotype-by-tissue explained variance (PVE) —such as HvARF20, HvARF13-3, HvARF18, and HvARF25—displayed strong eQTL signals (–log₁₀P > 3) (Fig. 3D, E). Notably, HvARF25 exhibited a distinct eQTL peak located 23–32 kb upstream of its transcriptional start site, coinciding with the position of a retrotransposon. This suggests that the TE insertion may act as a cis-regulatory element modulating HvARF25 expression in a genotype-dependent manner (Fig. 3B, D).

A particularly interesting outlier was HvARF11, which showed the most significant eQTL signal despite having low genotype and genotype-by-tissue PVE, and lacking annotated TEs within its vicinity (Fig. 3D, E). This implies the presence of a potent, localized regulatory variant—potentially an unannotated TE insertion, structural variation, or epigenetic modification—contributing substantially to HvARF11 expression, independent of broader genomic background effects. Overall, TE located near ARF genes are likely to play a regulatory role contributing to genotype-specific expression differences.

Pan-transcriptomic analysis reveals tissue-specific expression and functional networks of ARF genes

To elucidate the functional roles and tissue-specific regulatory dynamics of ARF genes in barley, we integrated population-scale pan-genomic (Jayakodi et al. 2024) and pan-transcriptomic (Guo et al. 2025) datasets. A total of 22 ARF genes with consistently detectable expression across five major tissues in 20 representative genotypes were retained for downstream analyses, minimizing noise from low-abundance transcripts.

ARF genes displayed pronounced tissue-specific expression patterns, with the inflorescence exhibiting the highest genotype-dependent expression variability, followed by roots (Fig. 4A). Notably, 12 of the 22 ARF genes were highly and specifically expressed in developing inflorescences at Waddington stages W6–W7 (Waddington et al. 1983), a critical window for spikelet and floret meristem differentiation. Seven genes showed peak expression in seedling roots, while HvARF2 and HvARF6 were predominantly expressed in shoots tissues (Fig. 4B). HvARF22 exhibited unique expression in caryopsis, and no ARF genes demonstrated tissue specificity in embryonic structures such as the embryo, mesocotyl, or seminal root. To further resolve the spatial expression dynamics within inflorescence tissue, we utilized a high-resolution single-cell RNA-seq atlas of barley inflorescence (Morex) (Demesa-Arevalo et al. 2025). This analysis uncovered cell-type-specific expression patterns for individual ARF genes (Fig. 4C). HvARF4-1 and HvARF4-2 were broadly expressed across diverse cell types, indicating general developmental roles. In contrast, HvARF14, HvARF15, and HvARF11 were enriched in vascular-associated tissues, with HvARF11 being specifically expressed during the transition from triple-spikelet meristem to floret meristem—an essential developmental checkpoint. Additionally, HvARF6, HvARF12, and HvARF25 were predominantly expressed in the adaxial epidermis, suggesting possible roles in organ polarity or boundary formation. These temporally and spatially resolved expression profiles support the hypothesis that ARF genes operate in multiple, partially distinct regulatory programs during inflorescence development.

To further explore the potential regulatory networks associated with those ARF gene, we conducted a weighted gene co-expression network analysis (WGCNA). Residual expression values (adjusted for genotypic effects) were used as input to focus on tissue-specific transcriptional relationships. For each ARF gene, co-expressed gene modules were identified, and Gene Ontology (GO) enrichment analysis was performed (Fig. 4D). As expected for transcription factors, most HvARF-associated modules were enriched in functions related to RNA/DNA binding and transcriptional regulation. Notably, HvARF3, HvARF14, and HvARF15 were co-expressed with genes involved in fungal defense responses. This observation is consistent with the functional identification of their rice orthologs (OsARF3, OsARF3a, OsARF3b), which have been reported to modulate disease resistance pathways.

Overall, our discovery indicates that ARF genes exhibit tissue and cell-type-specific expression patterns and are involved in distinct regulatory modules that govern a range of biological processes—from developmental patterning in inflorescence to biotic stress responses. These results insight a foundational framework for future functional studies and genetic manipulation strategies for optimizing inflorescence architecture or resistance in barley.

In this study, we present a novel method of gene family analysis by utilizing abundant, publicly available genomic and transcriptomic resources to comprehensively investigate the ARF gene family in barley. Through the integration of multi-omics datasets, we performed a series of systematic analyses including PAV detection, evolutionary dynamics, TE-mediated eQTL mapping, transcriptome analysis, and domestication signal identification. A summary of ARF genes in barley has been created by using different data source and methods, which provide insights and clues for future cloning and analysis. Those analyses, which were previously impossible with a single reference genome, demonstrate the advantage of pan-genome, pan-transcriptome, and resequencing data in mining gene family complexity across diverse genotypes in barley.

Recent advances in barley multi-omics resources, including high-resolution pan-genomes and population-wide transcriptomes, offer unprecedented opportunities for gene family research. However, the challenge lies in effective data integration and interpretation. Our study demonstrates that combining these layers enables the discovery of functional alleles and regulatory mechanisms otherwise inaccessible. For instance, the identification of HvARF3^R475G as a putative domestication-related variant associated with increased grain size and TKW underscores the value of population-scale pan-genomic approaches. Conventional gene family studies relying on single genomes often fail to capture intraspecies diversity— especially in crops like barley with global distribution and rich germplasm pools. Pan-genomics strategies, in contrast, enable the identification of both structural and nucleotide-level polymorphisms, accelerating trait mapping and gene discovery. Previously, genes such as Ppd-H1 and VRN-1 were identified through large-scale germplasm screenings (Santra et al. 2009; Turner et al. 2005). Today, pan-genomic datasets allow for earlier and more targeted discovery. In this study, we revealed geographic stratification of HvARF3 haplotypes—paralleling patterns observed in well-characterized flowering regulators like wheat Ppd-D1—highlighting possible regional adaptation and breeding selection.

The ARF gene family comprises a moderate number of transcription factors with essential roles in auxin-mediated growth and development. While several ARFs have been functionally characterized in Arabidopsis and rice, their roles in barley remain underexplored. Phylogenetic analyses revealed lineage-specific expansions and contractions. Notable, the AtARF17/16/10 clade exhibited multiple duplications in barley, particularly within the HvARF13 subgroup. HvARF13-3 remains highly conserved, while HvARF13-1 and HvARF13-2 show CNVs and sequence divergence—suggesting ongoing sub-functionalization or neofunctionalization (Fig. 1B). These variations are more frequent in wild and landrace accessions, implying adaptive relevance in non-domesticated environments. Comparative genomics also revealed notable differences in ARF gene retention and loss. For example, rice contains two paralogs, OsARF23 and OsARF24 (formerly OsARF1), involved in grain formation via rice morphology determinant (RMD) regulation (Li et al. 2014; Waller et al. 2002). These genes have no direct orthologs in barley (Fig. 2A), indicating species-specific divergence. Conversely, barley harbors a duplicated OsARF4 ortholog pair—HvARF4-1 and HvARF4-2—both retaining full-length domain structures. Given OsARF4’s known role in negatively regulating grain size in rice (Hu et al. 2018), its barley counterparts represent strong candidates for functional studies.

Despite previous findings that several rice ARF genes (e.g., OsARF3, OsARF4, OsARF6, OsARF14, OsARF15) are associated with grain-related traits (Gu et al. 2023; Hu et al. 2018; Qiao et al. 2021), most ARF genes in barley do not show direct expression in the caryopsis (Fig. 4B). Instead, their strong expression in inflorescence tissues suggests that these ARF genes may regulate grain development indirectly by modulating early reproductive organogenesis and floral architecture, which likely affects final grain size or weight. Another example is HvARF3, the selected gene during breeding, it shows high gene expression levels in the developing inflorescence while no expression in grain. In the single-cell transcriptome dataset of barley inflorescences, HvARF3 is specifically expressed in dividing cells, suggesting a potential role in regulating cell division and meristem maintenance during early stages of inflorescence development. But according to GWAS haplotype analysis and performance of its ortholog OsARF3, it most likely affects grain size and grain traits, too. Interestingly, the favorable haplotypes of HvARF3 are rare in Asian barley accessions. This may be attributable to a regional genetic bottleneck or local adaptation, which limited the introgression of the advantageous haplotype. Alternatively, certain HvARF3 alleles may exhibit functional trade-offs similar to its rice ortholog OsARF3, which has been shown to mediate a balance between heat tolerance and pathogen resistance (Gu et al. 2023). Furthermore, the Ka/Ks rate showed that HvARF3 was higher in Asian accessions, indicating that the purifying selection was not as strong as that in Europe. Thus, environmental difference between European and Asian growth conditions could have shaped regional selection pressures acting on different HvARF3 variants and reflect different breeding strategy.

In conclusion, our integrative analysis highlights the critical role of ARF genes in barley development and demonstrates the utility of multi-omics strategies in gene family research. By resolving patterns of structural variation, expression regulation, and selection, we identify promising candidates like HvARF3 for targeted breeding. These findings contribute not only to our understanding of auxin signaling but also to practical efforts aimed at improving grain yield and adaptability in barley. The dCAPS marker developed for the HvARF3^R475 variant will facilitate marker-assisted selection in future breeding programs.