- 1Faculty of Agriculture, Forestry and Food Engineering, Yibin University, Yibin, Sichuan, China
- 2Crop Research Institute, Shandong Academy of Agricultural Sciences, Jinan, Shandong, China
- 3State Key Laboratory of Crop Gene Exploration and Utilization in Southwest China, Sichuan Agricultural University, Chengdu, Sichuan, China
- 4Key Laboratory of Aquatic Genomics, Ministry of Agriculture and Rural Affairs, Beijing Key Laboratory of Fishery Biotechnology, Chinese Academy of Fishery Sciences, Beijing, China
The NAC transcription factor family is pivotal in regulating plant development and stress responses, yet its diversity and evolutionary dynamics in barley (Hordeum vulgare) remain underexplored. In this study, we performed a comprehensive pan-genome analysis to identify and characterize the HvNACs across 20 barley accessions. A ranging from 127 to 149 HvNACs were identified in each genome, in which the Morex genome harbored the highest count. These HvNACs were classified into 201 orthogroups, further stratified into core (102), soft-core (18), shell (25), and lineage-specific (56) categories. Phylogenetic analysis delineated them into 12 subfamilies, of which the core genes have undergone strong purifying selection, by contrast, the shell and lineage-specific genes were under relaxed selection constraint, suggesting functional diversification in barley. Genomic variation, such as PAVs and CNVs, largely driven by TEs, highlighted the dynamic nature of NAC loci. Furthermore, transcriptome profiling of the HvNACs demonstrated diverse tissue expression patterns and different response characteristics under salt stress. These findings elucidate the evolutionary and functional dynamics of HvNACs, offering valuable insights for genetic improvement of breeding programs in barley as well as in other crops.
Introduction
NAC transcription factors, emerged during the transition from aquatic to terrestrial plant life and expanded thereafter, establishing themselves as one of the largest transcription factor families in plants (Mathew and Agarwal, 2018; Chen and Xia, 2025). They modulate an array of biological processes, such as secondary growth, cell division, senescence, hormone signaling, and nutrient uptake (Xiong et al., 2025). Moreover, NAC transcription factors coordinate multiple developmental and environmental signals, creating complex regulatory networks that precisely adjust target gene expression (Li et al., 2011; Nuruzzaman et al., 2013; Parankusam et al., 2017). These proteins possess a conserved N-terminal NAC domain (~150 amino acids) and a highly variable C-terminal transcriptional regulation region, which exhibit functional specificity through their structural diversity The NAC domain consists of five subdomains (designated A to E), each performing distinct molecular functions (Mohanta et al., 2020). Of these, subdomains C and D are highly conserved and essential for DNA binding, while subdomain A supports homo- and heterodimerization, enabling interactions with other NAC proteins. By contrast, the more variable subdomains B and E promote the functional diversification of NAC proteins across plant species (Kim et al., 2016). The C-terminal region shows considerable sequence variation and functions as a transcriptional activation or repression domain, often interacting with other regulatory proteins to modulate gene expression (Hu et al., 2010; Mohanta et al., 2020). This dual-domain structure supports the diverse roles of NAC transcription factors in regulating essential physiological processes and aiding stress adaptation.
NAC transcription factors serve as pivotal regulators of plant responses to both biotic and abiotic stresses (Sun et al., 2012; Shao et al., 2015). These plant-specific transcription factors feature conserved DNA-binding domains and function as either activators or repressors in a variety of cellular processes (Singh et al., 2021). They engage in intricate signaling networks and transcriptional reprogramming crucial for stress adaptation (Nuruzzaman et al., 2013). Notably, NAC transcription factors often exhibit auto-regulation and cross-regulation, enabling them to coordinate responses to multiple stresses simultaneously (Yuan et al., 2019). For example, ANAC019, ANAC055, and ANAC072/RD26 have been demonstrated to enhance drought tolerance in Arabidopsis thaliana (Tran et al., 2004). The ANAC072/RD26 functions as transcription activator in ABA-dependent stress-signaling pathway (Fujita et al., 2004). In rice (Oryza sativa), the stress-responsive SNAC3 boosts heat and drought tolerance by regulating reactive oxygen species homeostasis (Fang et al., 2015). Moreover, ANAC096 and bZIP transcription factors work together to mediate dehydration and osmotic stress responses in Arabidopsis (Xu et al., 2013). Overexpressing of grapevine VvNAC08 and wheat (Triticum aestivum) TaNAC29 has also been shown to enhance drought and salt tolerance in transgenic Arabidopsis (Xu et al., 2013). In addition to abiotic stress tolerance, NAC transcription factors play vital roles in plant immunity. In wheat, the NAC protein GRAB1–2 binds to the wheat dwarf geminivirus RepA protein, regulating viral DNA replication (Xie et al., 1999). Considering the increasing environmental challenges in modern agriculture, the NAC transcription factors offer promising avenues for improving crop resilience. Overexpressing stress-responsive NAC genes has been successfully used in transgenic approaches to improved tolerance to various stresses (Nuruzzaman et al., 2013). However, precise promoter selection and rigorous field validation are essential to mitigate potential trade-offs between stress tolerance and plant growth (Shao et al., 2015). Future research should aim to unravel the intricate regulatory mechanisms of NAC transcription factors and refine their application in breeding stress-resilient crop varieties.
Genome-wide studies have identified a large number of NAC genes across diverse plant species, such as there are 105 NAC genes in Arabidopsis (Ooka et al., 2003), 154 in tobacco (Nicotiana tabacum) (Li et al., 2018), 151 in rice(Oryza sativa) (Nuruzzaman et al., 2010), 559 in wheat (Guerin et al., 2019), and 147 in foxtail (Setaria italica) (Puranik et al., 2013). Furthermore, a comparative genomic analysis of 24 land plant species revealed a total of 3,187 NAC transcription factors, which were clustered into six major groups. These findings highlight the evolutionary complexity of the NAC gene family, marked by variations in gene number and phylogenetic relationships across species. Traditional gene family identification often relies on identifying conserved domains within a single reference genome. However, this approach fails to capture the full extent of genetic diversity within species, particularly gene presence/absence variation (gPAV), which is common among individuals (Tettelin et al., 2005). To gain a more comprehensive understanding of species-wide diversity in genome, the concept of the pan-genome has emerged and gained significant attention. Unlike traditional approaches that rely on a single reference genome, pangenome overlook gene presence/absence variation among individuals, which integrates genomic data from multiple accessions, reconstructing the species’ complete gene repertoire (Jayakodi et al., 2021). There are two categories for genes included in a pangenome, the core genes shared across all individuals and the dispensable/variable genes present only in subsets of accessions (Bayer et al., 2020). Pan-genome analysis provides a robust framework for exploring intra-species genetic variation, adaptive evolution, domestication, and functional diversity (Tong et al., 2025).
Barley (Hordeum vulgare), one of the world’s oldest cultivated cereal crops, ranks the fourth-largest crop after maize, rice, and wheat in terms of yields. It possesses broad adaptability to diverse environments and plays essential roles in feeding, malting, and brewing. Cultivated barley (H. vulgare subsp. vulgare) was domesticated approximately 10,000 years ago from its immediate ancestor wild barley (H. vulgare subsp. spontaneum) (Zohary et al., 2012), which exhibits substantial genetic diversity potentially improving barley cultivars. Barley has been well developed in genetics, genomics, and breeding over the past two decades, particularly in pan-genomics in recent five years. The pan-genome assembly consisting of 19 cultivated accessions and one wild barley accession has been published (Jayakodi et al., 2020), and subsequently released a barley pan-transcriptomic datasets of 20 barley genotypes (Guo et al., 2025). The published barley pan-genome have facilitated functional genomics analyses, several studies have performed a pan-genomic approach for gene-family characterization and exploring structural variations during domestication process in barley (Chen et al., 2022; Jayakodi et al., 2024; Tong et al., 2025). In this study, we systematically identified the NAC gene family members in barley at the pan-genome level. We further analyzed gene structures, physicochemical properties, chromosomal localization, phylogenetic relationships, and tissue-specific expression patterns. These approaches addressed the limitations of single-reference genome studies and provided a more comprehensive understanding of the HvNACs’s evolutionary dynamics and functional diversity in barley. Moreover, we explored the potential roles of HvNACs in salt stress responses, offering new insights into stress tolerance mechanisms and a theoretical basis for genetic breeding improvement in barley.
Materials and methods
Identification of HvNACs
The genomic sequences of 20 barley pangenome accessions were retrieved from the IPK database (https://barley-pangenome.ipk-gatersleben.de) as described by Jayakodi et al. (2020). The genome and protein sequence data of rice (Oryza sativa. IRGSP-1.0) were downloaded from the Ensembl database (http://asia.ensembl.org/index.html). The NAM domain (PF02365) hidden Markov model (HMM) file was downloaded from the Pfam website (https://pfam.xfam.org/). Using HMMER (v3.0) software with an E-value threshold of 1E-5 and a minimum domain score of 20, as described by Chen et al. (2022), HvNAC protein sequences containing the NAM domain were identified. To further confirm NAC domain presence, BlastP (v2.9.0) was used with an E-value threshold of 1E-10 and a minimum query coverage of 70% against the Pfam NAM domain sequence. These potential protein sequences were further validated for the presence of the NAC domain using the NCBI Conserved Domain Database (CDD; https://www.ncbi.nlm.nih.gov/cdd/), the SMART database (http://smart.embl-heidelberg.de/), and the Pfam database.
Orthologous HvNACs clusters across the 20 barley pangenome accessions were inferred using the CD-HIT (v4.6) (Li and Godzik, 2006) with sequence identity threshold of 95% and alignment coverage of 90% for the longer sequence, requiring a minimum alignment length of 100 amino acids. Orthologous gene groups (OGGs) were defined as clusters of NAC genes sharing at least 95% sequence identity and were further classified based on their presence across accessions: core OGGs (present in all 20 accessions), soft-core OGGs (present in 18–19 accessions), shell OGGs (present in 3–17 accessions), and cloud OGGs (present in only one accession). To calculate copy number variation (CNV), we counted the number of gene members in each OGG for every barley accession. The resulting OGG-by-accession matrix was used to compare gene copy numbers across the 20 genomes. The identified pangenome HvNACs were sequentially numbered and sorted according to their conservation category and the order of their corresponding OGGs. The longest protein in each cluster was designated as the representative sequence. The physicochemical properties of the NAC proteins were computed using the ExPASy online tool (https://web.expasy.org/protparam/). To visualize the presence/absence patterns of NAC genes across the 20 barley genomes, we used the R package ComplexHeatmap (Gu et al., 2016) to generate a heatmap. Additionally, ggplot2 (Villanueva and Chen, 2019) was used to create plots summarizing the frequency distribution of NAC gene occurrence across accessions.
Phylogenetic tree construction
The phylogenetic tree of NAC proteins was constructed using the Maximum Likelihood (ML) method implemented in IQ-TREE v2.0 (Nguyen et al., 2015). Full-length amino acid sequences of NAC proteins from barley and rice were first aligned using MUSLE. The resulting multiple sequence alignment was then trimmed using trimAl with the automated option to remove poorly aligned regions and reduce noise. The best-fit substitution model (Dayhoff) was selected automatically by IQ-TREE. To assess the robustness of phylogenetic relationships, bootstrap analysis was performed with 1,000 replicates. Based on the clustering results of the phylogenetic tree and previously reported classifications in rice and miaze (Ooka et al., 2003; Peng et al., 2015), the NAC proteins were assigned to distinct subfamilies. The final phylogenetic tree (*.nwk file) was uploaded to the iTOL platform (http://itol.embl.de/) for enhanced visualization and further annotation.
Analysis of gene structures and conserved motifs
Intron-exon structures for representative HvNACs were obtained from the annotation files corresponding to 20 barley genomes. Additionally, conserved motifs within the HvNACs of barley were identified utilizing the MEME Suite (http://meme-suite.org/meme), specifying a maximum limit of 10 motifs. Subsequently, the resulting gene structures and motif distributions were visualized employing TBtools (v2.097) (Chen et al., 2020).
Analysis of duplication events and collinearity
Synteny analysis for HvNACs was conducted using MCScanX (Wang et al., 2013), utilizing amino acid sequences and corresponding chromosomal positions from the Morex V3 genome assembly (https://wheat.pw.usda.gov/GG3/) (Mascher et al., 2021) These syntenic relationships were subsequently visualized using a circos plot generated with TBtools (v2.097) (Chen et al., 2020). To further explore evolutionary selection pressures, synonymous (Ks) and non-synonymous (Ka) substitution rates were calculated for orthologous gene pairs across 20 barley accessions. This estimation utilized the approximate method available in KaKs Calculator (v1.2) (Zhang et al., 2006). A Ka/Ks ratio greater than 1 indicates positive selection, a ratio equal to 1 suggests neutral selection, and a ratio less than 1 signifies purifying selection.
Analyses of presence/absence variation and transposable element
Datasets for PAV and TE covering the Morex genome and 19 additional barley genomes were obtained from the source indicated by https://doi.org/10.5447/ipk/2020/24. Subsequently, custom Python scripts were employed to identify TE regions overlapping HvNACs. These overlapping TE elements were then categorized based on their position relative to HvNACs as upstream, downstream, or genic. Furthermore, the transposons themselves were classified according to the three-letter code system outlined by Wicker et al. (2007).
Codon usage evaluation
Codon usage bias was analyzed using codonW (v1.4.2) (http://codonw.sourceforge.net/). The following parameters were calculated: effective number of codons (ENC), codon adaptation index (CAI), relative synonymous codon usage (RSCU), overall genomic G+C content (GC), G+C content at the third synonymous codon position (GC3s), and nucleotide composition at the third synonymous codon position (T3s, C3s, A3s, G3s). Optimal codons were identified based on genome-wide ENC values. Genes in the lowest and highest 10% ENC categories were classified as high-expression and low-expression gene sets, respectively. Optimal codons were defined as those with an RSCU > 1 in high-expression genes, RSCU < 1 in low-expression genes, and a ΔRSCU (RSCU difference) ≥ 0.3.
To assess the relative contributions of natural selection and mutational pressure on codon usage bias, Parity Rule 2 (PR2) analysis was performed. PR2 plots were generated by calculating the A3/(A3 + T3) ratio and G3/(G3 + C3) ratio at the third codon position for each HvNAC across the 20 barley accessions. A reference point at (0.5, 0.5) was included to indicate a state of equal mutational pressure and no bias between purines and pyrimidines. ENC-plot analysis was conducted to evaluate codon usage patterns and the effect of GC content at the third synonymous codon position (GC3s). ENC values were plotted against GC3s for each NAC gene, and the expected ENC values under a neutral mutation model were computed using the formula:
Genes with ENC values falling below this curve were inferred to be under selection rather than mutational pressure. All plots, including PR2 and ENC analyses, were visualized using Matplotlib (v3.5.2) (Barrett et al., 2005), with consistent color coding across barley accessions to facilitate comparative analysis.
Transcriptome analysis
Raw transcriptome sequencing data from five distinct tissues (caryopsis, coleoptile, inflorescence, root, and shoot) across 20 barley genomes were utilized in this study, sourced from a recent comprehensive barley transcriptome investigation (Guo et al., 2025). Initially, raw sequencing reads underwent processing with Trimmomatic (v0.39) (Bolger et al., 2014) to eliminate adapter sequences, low-quality reads, and potential contaminants, thereby ensuring high data integrity. Subsequently, the filtered reads were aligned to the Morex V3 reference genome assembly (https://wheat.pw.usda.gov/GG3/) (Mascher et al., 2021) using HISAT2 (v2.2.1) (Kim et al., 2019), noted for its efficient and accurate mapping capabilities. Transcript assembly and quantification were then performed using StringTie (v2.1.3) (Pertea et al., 2015) based on the alignment results. Gene expression variation across the different tissues was assessed by calculating the average FPKM value for each HvNAC across all 20 barley genomes. Finally, these expression patterns were visualized as a heatmap generated with the ComplexHeatmap package in R (v4.4.0) (Gu et al., 2016), facilitating the analysis of tissue-specific expression profiles and potential functional divergence among the HvNACs.
Analysis of gene expression clustering
To investigate the salt stress response of HvNACs in barley, transcriptome data of leaf tissues subjected to varying durations of salt treatment (0 h, 1 h, 24 h, and 10 d) were retrieved from the SRA database (accession ID: PRJNA962512). In this experiment, barley seedlings were cultivated under controlled conditions until the three-leaf stage, at which point they were treated with 300 mM NaCl to induce salt stress. Leaf samples were collected at the designated time points following treatment for RNA sequencing and subsequent expression analysis. Expression patterns across these time points were analyzed using the Mfuzz package (Kumar and Futschik, 2007), as previously described by Chen et al. (2023), resulting in the classification of genes into 10 distinct clusters based on their temporal expression profiles. The fuzzification parameter (m) was optimized at 1.25 to balance cluster membership stringency and flexibility. To explore the co-expression relationships, Pearson correlation analysis was performed to calculate the correlation between HvNACs in each cluster and other genes, with a threshold of correlation coefficient greater than 0.9 and p-value less than 0.01 applied to identify highly correlated genes. Additionally, NAC transcription factor binding motifs were downloaded from the PlantTFDB database (https://planttfdb.gao-lab.org/prediction.php). The FIMO (v5.3.3) (Grant et al., 2011) was employed to search for these motifs within the promoter sequences of genes highly correlated with HvNACs, enabling the identification of potential regulatory targets. Visualization of the co-expression networks and motif analysis results was performed using Cytoscape (v3.10.2) (Su et al., 2014). Gene Ontology (GO) enrichment analysis was subsequently performed to elucidate the potential functional roles associated with these modules. Gene Ontology terms for the enrichment analysis were retrieved from the KOBAS database (http://bioinfo.org/kobas). Visualization of the results was performed using relevant R packages to facilitate a comprehensive understanding of the molecular mechanisms governing the salt stress response of HvNACs.
Results
Identification of HvNACs in 20 barley pangenome
Using HMMER and BlastP, we identified HvNACs across 20 barley genomes. To eliminate redundancy and facilitate downstream analyses, we used CD-HIT to cluster similar sequences. For each cluster, one representative gene per accession was retained. The number of HvNACs in each genome ranged from 127 (HOR_21599) to 149 (Morex), with most genomes (85%) containing 130–135 HvNACs. To further characterize these HvNACs, we conducted a comprehensive analysis of their physicochemical properties (Supplementary Table S1). The sequence lengths varied significantly, ranged from 100 amino acids to 863 amino acids, while their molecular weight spanned from 11,854.43 Da to 96,275.98 Da. The pI ranged from 4.29 to 10.57, indicating a broad spectrum of charge distributions. The instability index ranged from 27.00 to 82.12, reflecting variations in protein stability. Meanwhile, the aliphatic index spanned from 42.46 to 73.75, suggesting differences in thermostability. The GRAVY values, all negative, ranged from -0.243 to -0.951, confirming the hydrophilic nature of these HvNACs. The considerable diversity in physicochemical properties highlights the potential functional versatility of barley HvNACs.
The Morex-v3 assembly contained a significantly higher number of NAC genes than the other 19 barley genomes, prompting an investigation into whether this expansion may have resulted from improvements in genome assembly quality. To explore this, we performed a comparative analysis between Morex-v3 and its predecessor, Morex-v2, which contains 137 NAC genes. Clustering of NAC genes from both versions revealed divergences. A total of 167 clusters were identified, including 63 clusters with single genes, 103 clusters forming 1:1 pairs, and one cluster containing three genes (Supplementary Table S2). These results suggest differences between assembly versions, potentially including resolved fragmentation, merged loci, or novel gene annotations in v3. The analysis indicates that advancements in genome assembly continuity may contribute to more accurate identification of genes, which might help explain the elevated NAC gene count in Morex-v3 relative to both its earlier version.
Core and dispensable HvNACs
To investigate the conservation of HvNACs in the barley pan-genome, we assigned the HvNACs identified from 20 barley accessions into different orthologous gene groups (OGGs). Using a threshold of ≥95% amino acid sequence identity, a total of 2,842 HvNACs from the 20 barley genomes were clustered into 201 orthologous gene groups (OGGs).The 201 OGGs were stratified into four conservation categories (Supplementary Table S3; Figures 1A, B): (1) Core clusters (n=102), exhibiting strict evolutionary conservation with all 20 genes present in every accession, likely maintaining essential cellular functions; (2) Soft-core clusters (n=18), displaying near-ubiquitous distribution (conserved in ≥90% lines) with limited lineage-specific losses; (3) Shell clusters (n=25), showing moderate conservation (present in >10% but <90% of the accessions) and missing in 3–17 genomes, indicative of evolutionary plasticity; and (4) Lineage-specific/cloud clusters (n=56), characterized by extreme diversification (≤10%) and absence in 18–19 accessions, potentially encoding adaptive traits. This conservation-to-divergence continuum not only delineates fundamental versus specialized NAC functions but also provides an evolutionary framework prioritizing candidate genes during barley improvement.

Figure 1. Number and distribution of HvNACs in the barley pangenome. (A) Statistics of the number of different types of NAC genes in 20 barley genomes. (B) Heatmap showing the presence and absence of NAC genes across 20 barley varieties. Each row represents a NAC gene, and each column represents a barley accession. Orange indicates the presence of a gene, while blue indicates its absence.
Phylogenetic analysis of HvNACs
Phylogenetic analysis of HvNACs in the barley pan-genome and rice (based on the Oryza sativa reference genome, IRGSP-1.0) identified 12 major subfamilies (group 1 to group 12), with groups 5 and 6 comprising only rice NAC genes, suggesting potential lineage-specific functional divergence (Supplementary Table S3; Figure 2). In barley, core and soft-core genes, collectively accounting for 81.2% of all NAC entries, were predominantly clustered in group 1 (23.7%), group 2 (15.3%), group 3 (19.5%), group 4 (14.4%), and group 8 (8.3%), suggesting conserved evolutionary roles in essential regulatory processes. By contrast, lineage-specific/cloud and shell genes, which comprised 18.8% of total entries, were enriched in groups 9 to 12. Notably, group 10 (7.1%) and group 11 (6.5%) contained the highest proportions of dynamic, species-specific genes, likely contributing to adaptive diversification. Interestingly, group 7 harbored only a single HvNAC gene. The absence of HvNACs in groups 5 and 6, along with their exclusive presence in rice, underscores the divergent evolutionary trajectories of these two Poaceae species. This pattern illustrates that conserved subfamilies retain essential functions, reflecting a balance between genetic stability and evolutionary innovation within the NAC gene family.

Figure 2. Phylogenetic tree of barley NAC genes. The evolutionary tree was constructed using IQ-TREE with the maximum likelihood method (1,000 bootstrap replicates) based on the Dayhoff model.
Analysis of gene structure and conserved motif
To investigate the evolutionary relationships and structural diversity of HvNACs within the barley pangenome, we performed a comprehensive analysis of their gene structures and conserved protein motifs (Supplementary Table S4). Our analysis reveals that HvNACs within the same subfamily typically display consistent exon-intron architectures and motif compositions, highlighting their close evolutionary relationships. Most HvNACs contain 1 to 6 exons. Motif analysis identified a conserved core set of motifs (Motif_0 to Motif_6) shared by most subfamilies. However, subfamily-specific deviations were evident; for example, group 10 genes (e.g., HvNAC.CR085) frequently lacked Motif_1 or included Motif_9, while group 11 genes (e.g., HvNAC.CR096) often contained Motif_7 or Motif_8, suggesting potential functional specialization. Dispensable HvNACs largely mirrored the motif patterns of core genes within their subfamilies; however, the presence of unique motifs in lineage-specific genes (e.g., HvNAC.CL177 in group 1 with Motif_9) suggests adaptive divergence.
Analysis of CNVs, PAVs, and transposon elements
Analysis of CNV in HvNACs across the barley pangenome revealed substantial variability, with gene copy numbers ranging from 1 to 58 (Supplementary Table S5; Figure 3A). Among the 201 HvNACs identified, 89 were single-copy core genes consistently present across all genomes, whereas some lineage-specific genes exhibited copy number expansion in specific cultivars. For instance, genes such as HvNAC.CR001 to HvNAC.CR017 maintained a uniform single copy across all lines. In contrast, HvNAC.CR021 exhibited up to four copies in cultivars including HOR_3081, HOR_7552, and Akashinriki, while HvNAC.CR057 reached seven copies in Morex. Further analysis revealed additional complexity for CNVs of barley NAC. HvNAC.CL162 and HvNAC.CL164 exhibited five copies exclusively in Morex, suggesting potential cultivar-specific amplification events. This variability reflects the dynamic evolutionary of HvNACs and may indicate functional diversification or adaptive significance associated with specific genomic backgrounds.

Figure 3. Copy number variation, presence-absence variation, and transposable element distribution of HvNACs. (A) Heatmap of NAC gene copy numbers across different barley genomes. (B) PAV overlapping with the NAC gene family, with two panels showing presence and absence compared to Morex. (C) Distribution of different types of TEs identified within 2 kb upstream and downstream of NAC genes; the left panel shows TE positions, and the right panel indicates TE types.
PAVs analysis of NAC gene regions across the barley pan-genome, based on a comparative assessment of 19 barley genomes against the Morex reference, identified a total of 567 distinct PAV events (Supplementary Table S6). These included both insertions and deletions relative to the Morex genome. For example, the locus Morex-HvNAC1 harbors a 58 bp insertion in the Golden Promise genome that is absent in Morex, whereas Morex-HvNAC4 contains an 886 bp region present in Morex but missing in Golden Promise (Figure 3B). Further analysis revealed extensive PAV complexity at specific NAC loci among different barley. The Morex-HvNAC4 locus exhibited substantial insertional variation, including a 6,822 bp segment in Barke and a 9,304 bp segment in Golden Promise, alongside deletions such as a 1,559 bp segment absent in Chiba. Similarly, the Morex-HvNAC7 locus demonstrated remarkable PAV diversity, ranging from a 58 bp deletion in Barke to a large 553,451 bp deletion in Golden Promise, suggesting major genomic rearrangements at this locus.
TEs were identified as key contributors to NAC gene duplication in barley. A comprehensive analysis revealed 1,377 TE insertions within or flanking HvNACs regions across the barley pan-genome (Supplementary Table S7). The positional distribution of these elements showed a prevalence in potential regulatory regions: 526 (38.2%) were found upstream, 502 (36.5%) downstream of the coding sequences, and 349 (25.3%) within genic boundaries (Figure 3C). Notably, DNA transposons constituted the vast majority (85.25%) of these identified TEs, although various retrotransposon types, including LTR elements as Ty1-copia retrotransposons (RLC) and Ty3 retrotransposons (RLG), were also present (Figure 3D). In addition, NAC-associated TEs were predominantly distributed toward the distal ends of chromosomes, with chromosome 7H showing the highest overall TE abundance associated with HvNAC regions (Supplementary Figure S1).
Analyses of natural selection
In the barley pangenome, we performed an evolutionary analysis of 201 OGGs within the HvNACs. The results revealed that core HvNACs exhibited significantly lower Ka, Ks, and Ka/Ks values compared to dispensable genes (p < 0.05) (Figure 4), suggesting that core genes are subject to stronger purifying selection. This intensified selective constraint likely reflects their indispensable roles in essential biological processes, thereby necessitating higher sequence conservation. In contrast, dispensable HvNACs appear to be under relaxed selection pressure, allowing greater evolutionary plasticity and enabling the acquisition of more diverse or specialized functions.

Figure 4. Comparison of Ka, Ks, and Ka/Ks values between core and dispensable NAC genes. P-values were calculated using the Mann-Whitney U test (non-parametric test).
Collinearity analysis
In the Morex barley reference genome, the HvNACs displays an uneven chromosomal distribution, with notable enrichment at the distal ends of chromosomes 2H, 3H, 5H, and 7H (Figure 5). These telomeric and subtelomeric regions—characterized by high recombination frequencies and elevated structural variation—are recognized as hotspots for rapid gene family expansion and functional innovation. A comprehensive analysis revealed that dispersed duplication is the primary mechanism driving NAC gene expansion, accounting for approximately 56.37% of the family members (Figure 5). This suggests that transposon-mediated replication and large-scale segmental rearrangements have played a central role in the relocation, dispersion, and functional diversification of HvNACs. Additionally, tandem duplication contributes to 18.12% of NAC gene expansion, emphasizing the importance of localized gene cluster amplification, likely driven by unequal crossing-over events, increasing gene copy number and supporting functional diversity. Collectively, the observed duplication patterns and chromosomal distribution highlight a complex evolutionary history for the barley NAC gene family, shaped by diverse duplication events and selective pressures.

Figure 5. Intraspecific collinearity analysis of NAC genes in the Morex genome. Gene types are color-coded as follows: Dispersed (red), Proximal (blue), Tandem (pink), and WGD or Segmental (orange). Collinear gene pairs are connected by red lines.
Codon usage evaluation
We observed distinct codon usage patterns of NAC gene lineages across the barley pan-genome. Analyzing 2,666 HvNACs revealed that the effective number of codons (ENC) ranged from 26.99 to 61.00, with an average of 42.79. As most ENC values exceeded 40, this suggests a generally low codon usage bias in HvNACs (Supplementary Table S8). Furthermore, the codon adaptation index (CAI) values ranged from 0.168 to 0.356, with an average of 0.253, significantly lower than the optimal value of 1. This indicates both weak codon preference and relatively low expression potential.
The relative synonymous codon usage (RSCU) values, as presented in Supplementary Table S9, reveal that 29 codons had an RSCU greater than 1. RSCU is a measure of how frequently a particular codon is used to encode an amino acid, relative to the expected usage if all synonymous codons for that amino acid were used equally (Sharp et al., 1986). An RSCU value greater than 1 indicates a preference for that codon in the gene set analyzed. Among these, the preferred codons were predominantly those ending in G (13 codons) and C (16 codons). In contrast, 28 codons were underrepresented, with an average RSCU value of less than 0.6, while 13 codons were overrepresented, with an average RSCU value greater than 1.6. By comparing the RSCU values between the two codon preference libraries in HvNACs, we identified four optimal codons: CCC (Proline), GCT (Alanine), GAT (Aspartic acid), and CTC (Leucine) (Supplementary Table S10). These codons showed a ΔRSCU greater than 0.3, with RSCU values exceeding 1 in high-preference genes and less than 1 in low-preference genes.

Figure 6. Codon usage bias analysis of NAC genes in the barley genome using ENC and PR2 plots. (A) ENC-plot analysis of NAC genes in the barley genome. The red line indicates the expected position of genes when codon usage is determined solely by GC3s. ENC represents the Effective Number of Codons; GC3s represents the G+C content at the third position of synonymous codons. (B) PR2-plot analysis of NAC genes in the barley genome. A3/(A3 + T3) represents the ratio of A to A+T at the third codon position. G3/(G3 + C3) represents the ratio of G to G+C at the third codon position.
To assess codon usage bias (CUB) in the analyzed genomes, a Parity Rule 2 (PR2) plot was generated, where A3/(A3+T3) was plotted on the vertical axis and G3/(G3+C3) on the horizontal axis for each gene (Figure 6B). The results revealed a slight deviation of data points from the expected neutral position (0.5, 0.5), with average values of 0.4951 for A3/(A3+T3) and 0.4610 for G3/(G3+C3), suggesting a moderate codon usage bias. Most G3/(G3+C3) values were below 0.5, indicating a preference for cytosine over guanine at the third codon position. This bias is likely shaped by both mutational pressure and natural selection. This conclusion is further supported by the significantly higher G3+C3 content (0.768) compared to A3+T3 (0.232), reflecting a marked GC preference in the analyzed genomes. In contrast, the average A3/(A3+T3) value of 0.4951, with individual genes reaching as high as 0.7, suggests a relatively balanced usage of adenine and thymine, with a slight preference for adenine in specific cases.
Expression profiling of HvNACs across different tissues
To explore the tissue-specific expression patterns of NAC transcription factors in barley, we analyzed the pan-transcriptomic transcript abundance of 201 HvNACs across five representative tissues consisting of caryopsis, embryonic, inflorescence, root, and seedling epicotyl. These genes were categorized into 12 groups based on their transcriptional signatures among all tissues. The results revealed a pronounced variation in expression profiles among different tissues (Figure 7). Notably, Group 7 (containing only one gene) displayed the highest average expression across all tissues, particularly in inflorescence (mean FPKM=91.34) and seedling epicotyl (mean FPKM=77.37). Group 3 and Group 8 also exhibited high expression levels across multiple tissues, with Group 3 showing high expression in roots (mean FPKM = 46.95). Byn contrast, Group 10 and Group 11 were characterized by low to moderate expression across all tissues. Group 10 genes showed relatively higher expression in roots (4.07) but remained low in other organs. Group 11, although widely distributed in the phylogenetic tree, showed consistently low expression levels (mean FPKM<1.0 in all tissues). Group 2 and Group 4 demonstrated moderate and constitutive expression across all tissues, indicating a more conserved and possibly fundamental role in general plant development. Group 1 genes, despite having a wide range of expression levels (standard deviation of 16.84–32.72), showed high expression in roots (mean FPKM =15.16) and caryopsis (mean FPKM =11.53). Overall, the transcriptomic landscape of HvNACs reveals both conserved and divergent regulatory roles in barley growth and development.

Figure 7. Expression analysis of NAC genes in different barley tissues, based on 20 barley accessions covering five tissue types (caryopsis, coleoptile, inflorescence, root, and shoot). Gene expression values represent the average transcript abundance across the 20 genomes; for non-core genes, the average was calculated only from the genome in which the gene was present. The heatmap was generated using min-max normalization to scale expression values to the range [0, 1] for each gene across tissues.
Mfuzz clustering analysis of HvNACs under salt stress
To characterize the dynamic transcriptional response of HvNACs to salt stress, we analyzed transcriptome data at four time points (0h, 1h, 24h, 10d) after salt treatment. Low-expression genes (average FPKM < 10) were filtered out to enhance analytical robustness. The remaining genes were then subjected to time-series clustering using the Mfuzz algorithm, which grouped them into 10 distinct clusters based on their expression profiles (Figure 8A). A total of 53 HvNACs were unevenly distributed across these 10 clusters. Cluster 2 contained the largest number (26 HvNACs), while Clusters 5 and 7 each included only one (Figure 8B). Further analysis of expression trends highlighted diverse regulatory roles of HvNACs. Notably, two contrasting profiles emerged: Cluster 6 showed gradual upregulation from 0h to 10d, whereas Cluster 8 exhibited progressive downregulation over the same period. These findings suggest that Cluster 2 (due to its high HvNAC enrichment), Cluster 6 (gradual upregulation), and Cluster 8 (progressive downregulation) represent key regulatory modules warranting further investigation into the molecular mechanisms responsible for salt stress response in barley.

Figure 8. Temporal expression patterns and clustering of HvNAC genes in barley leaves under salt stress. (A) Mfuzz clustering analysis in barley leaves under salt treatment at 0h, 1h, 24h, and 10d, with genes showing expression levels greater than 10 into 10 clusters. (B) The number of HvNACs in each cluster.
Using Pearson correlation analysis, we identified genes highly correlated with HvNACs (correlation coefficient > 0.9, p-value < 0.01) within specific clusters under salt stress in barley (Supplementary Table S11). The promoter regions of these genes were analyzed for NAC binding sites using FIMO, employing 21 HvNACs binding motifs retrieved from the PlantTFDB database (Supplementary Table S12). The results revealed 1,075 genes in Cluster 2, 244 genes in Cluster 6, and 362 genes in Cluster 8 with predicted NAC binding sites and high correlation with HvNACs (Figure 9A), indicating potential transcriptional regulatory networks associated with salt stress response.

Figure 9. Gene co-expression network and functional enrichment analysis of salt-responsive HvNACs clusters. (A) Network diagram of genes in clusters 2, 6, and 8 that are highly correlated with HvNACs and contain NAC binding sites in their promoter regions. (B) GO enrichment analysis of genes in clusters 2, 6, and 8.
GO enrichment analysis of the three key clusters uncovered significant biological processes related to salt stress (Figure 9B). In cluster 2, enriched GO terms such as signaling receptor activity (GO:0038023), transferase activity (GO:0016740), carbohydrate homeostasis (GO:0033500), and beta-glucan biosynthesis (GO:0051274) indicate involvement in salt stress signal transduction, carbohydrate metabolism, and cell wall integrity maintenance, with additional roles in abscisic acid signaling and responses to osmotic stress, salinity, and dehydration. Cluster 6 showed enrichment in ion binding (GO:0043167), cation binding (GO:0043169), and response to salt (GO:0009651), highlighting functions in ion transport, homeostasis, and stress adaptation. Conversely, a downregulated cluster 8 exhibited significant enrichment in catalytic activity (GO:0140101), acting on a tRNA (GO:0140101), ligase activity (GO:0016874), and chloroplast (GO:0009507), suggesting reduced RNA processing and photosynthesis-related activities to conserve energy and minimize oxidative damage. Enriched biological processes like response to salt stress (GO:0009651) and carboxylic acid transport (GO:0046942) in this cluster 8 may imply a strategic suppression of non-essential pathways to prioritize stress response resilience, enhancing adaptation to salinity in barley.
Discussion
The NAC transcription factor family member is a key regulator of plant development and abiotic stress responses across diverse species (Shao et al., 2015). Understanding the evolutionary dynamics and functional diversification of NAC genes in cereal crops like barley is critical for developing stress-resilient varieties through targeted breeding (Fuertes-Aguilar and Matilla, 2024). Our systematic investigation of HvNACs at pan-genome level gives novel insights into their evolutionary trajectories and functional specialization patterns in barley.
Pangenome analysis reveals significant diversity and impact of genome assembly quality
Through systematic analysis of NAC transcription factors in 20 barley pangenome datasets, we obtained critical insights into their genomic architecture, evolutionary patterns, and functional characteristics. A total of 127–149 HvNACs were identified in each barley genome, with the majority of accessions containing 130–135 functional members. These variations reflect the natural genetic heterogeneity characteristic of barley pangenome systems, highlighting the rich genetic diversity available for breeding programs aimed at improving specific traits. Compared with other species, the number of NAC genes in barley is higher than in A. thaliana (105 genes) (Ooka et al., 2003), comparable to rice (151 genes) (Nuruzzaman et al., 2010); and foxtail millet (147 genes) (Puranik et al., 2013), but lower than in wheat (559 genes) (Guerin et al., 2019). Significant variations in HvNACs characteristics (sequence length, molecular weight, isoelectric point, instability index) may indicate extensive functional diversity and adaptive regulatory capacity under environmental stresses (Puranik et al., 2012). The Morex V3 assembly showed maximum NAC gene representation (149), highlighting the critical influence of assembly quality on accurate gene prediction. This conclusion was further corroborated by comparative analysis between Morex-v3 and its predecessor Morex-v2. The increased NAC gene complement in Morex-v3 probably stems from improved assembly continuity, resolution of fragmented genomic regions, and enhanced novel gene annotation. These advancements enabled detection of NAC genes that were either fragmented or undetected in previous assemblies.
Evolutionary stratification into core and dispensable gene sets
The 201 NAC OGGs were systematically categorized into four evolutionary classes: core genes (present in all accessions), soft-core (≥90% frequency), shell (5 present in >10% but <90% of the accessions), and lineage-specific/cloud genes (≤10% prevalence), establishing a comprehensive evolutionary classification system. Collectively representing 81.2% of HvNACs, core and soft-core genes were primarily localized to phylogenetic group 1–4 and clade 8, suggesting their fundamental role in maintaining conserved regulatory networks in barley. The strong conservation of these clades, evidenced by lower Ka/Ks values suggestive of purifying selection, confirms their essential involvement in core cellular processes including developmental regulation and stress response pathways (Panchy et al., 2016). In contrast, shell (25 OGGs) and lineage-specific/cloud genes (56 OGGs) constitute an adaptive reservoir characterized by incomplete pan-genome representation and relaxed selection constraints. This gene cohort likely facilitates environmental adaptation and niche-specific stress tolerance, representing promising targets for future functional studies and breeding efforts to enhance stress resilience. These dispensable genes may represent a valuable resource for identifying novel alleles associated with adaptation to specific local environments or emerging stresses, offering promising avenues for targeted breeding strategies to enhance resilience in diverse agricultural settings (Dawson et al., 2015). Comparative phylogenetics revealed lineage-specific divergence, with group 5 and 6 uniquely conserved in rice NACs. These evolutionary characteristics may imply Poaceae-specific functional specialization, possibly driven by differential domestication histories or ecological adaptation requirements (Doebley et al., 2006).
Genomic rearrangements and transposable elements drive HvNACs dynamic expansion
CNVs and PAVs detected in HvNACs demonstrate their evolutionary plasticity. This variability spans conserved single-copy core genes to genotype-specific expansions, exemplified by HvNAC.CR057 and HvNAC.CL162/164 amplifications in Morex, may implying dosage-dependent regulation of NAC-mediated pathways. PAVs involving segmental insertions/deletions indicate ongoing genomic restructuring, a phenomenon also reported in wheat NAC genes (Singh et al., 2021).This CNV and PAV diversity may be partly driven by differential selection pressures associated with barley’s adaptation to diverse agro-ecological environments during domestication and improvement (Pourkheirandish and Komatsuda, 2007; Russell et al., 2016). As barley has been cultivated across a wide geographic range—from arid regions in the Fertile Crescent to temperate zones in Europe and East Asia—local adaptation likely favored the expansion or retention of lineage-specific NAC genes with stress-related functions, particularly those associated with drought, salinity, and cold tolerance (Dawson et al., 2015). Additionally, repeated selection bottlenecks and introgressions during modern breeding may have contributed to gene copy variation in elite cultivars.
A strong association emerged between NAC gene diversification and TE mobilization. DNA transposons (85.25% of TE content) predominantly occupied NAC gene flanking regions, with significant enrichment in 5’/3’ regulatory sequences. Spatial correlation between TEs and CNVs suggests TE-mediated local duplication mechanisms, consistent with Poaceae genome evolution patterns (Bennetzen and Wang, 2014; Wicker et al., 2018). In this study, retrotransposons exhibited lower abundance, while retrotransposons are responsible for increasing genomic flexibility as revealed previously (Chen et al., 2023). Collectively, these data suggest that TEs act as primary architects of NAC gene diversity, which likely contributes to barley adaptation coping with environmental stresses. The observation that dispersed duplication is the dominant mechanism (56.37%) for NAC expansion, coupled with enrichment at recombination-prone distal chromosome ends, further supporting the role of TEs and genome instability in shaping the molecular evolution of HvNACs. Understanding the interplay between TEs and NAC gene evolution provides critical insights into the mechanisms driving adaptive evolution in barley, potentially enabling future strategies for manipulating TE activity to generate beneficial genetic variations for crop improvement.
Functional insights from structure, codon usage, and expression patterns
The consistent exon-intron structures and conserved motifs within NAC subfamilies suggest a common evolutionary origin. However, subtle differences among clades, such as the absence of Motif_1 in Group 10, the presence of Motif_7 and Motif_8 in Group 11, and the emergence of unique domains in lineage-specific genes, may reflect functional specialization and adaptive evolution of HvNACs.
Codon usage analysis revealed a generally weak codon bias among HvNACs, indicating that the translation efficiency of the most genes is not likely to undergo strong selective pressure. The identified preferred codons provide valuable clues for optimizing gene expression in future genetic engineering applications. At the same time, the neutral deviations observed in the ENC and PR2 plots, may implying both mutational pressures, especially GC bias, and possible translational selection contribute to codon preference. These factors may modulate expression efficiency or influence protein folding of certain genes under specific conditions. Furthermore, a set of preferred and optimal codons, predominantly ending in G or C, was identified. As this analysis was based on 20 barley genomes, the limited sample size may constrain the generalizability of the codon usage patterns observed. A larger genomic dataset could enhance the robustness of these findings. This analysis provides preliminary insights into the evolutionary forces influencing codon preferences in HvNACs, with potential applications for gene expression optimization.
Tissue expression profiling revealed that HvNACs show pronounced expression differences across various tissues, and these expression patterns are closely associated with their phylogenetic classifications. Genes belonging to Group 3 and Group 8 showed high expression levels in roots and inflorescences, suggesting their roles in tissue-specific developmental processes. This pattern is consistent with the reported functions of NACs in tissue differentiation and secondary cell wall formation (Olsen et al., 2005; Zhong et al., 2010). Genes in Group 2 and Group 4 exhibited relatively stable expression patterns, suggesting potential housekeeping roles. In contrast, the low expression levels of Group 10 and Group 11 suggest possible involvement in developmental stages or stress conditions not examined in this study, or in regulatory functions that require minimal transcriptional activity. Moreover, the considerable variation in expression within Group 1 highlights the functional diversity in HvNACs.
HvNACs involve in salt stress response
Using Mfuzz clustering analysis, we identified the dynamic transcriptional responses of HvNACs to salt stress, supporting their essential roles in barley adaptation coping with abiotic stress. Interestingly, among all clusters, only cluster 2 and cluster 8 each contained a single gene classified as soft-core, while the remaining HvNACs in these stress-responsive clusters were core genes. This pattern indicates that the transcriptional response to salt stress is primarily mediated by evolutionarily conserved NAC genes, and that core NACs may play dominant roles in the regulation of key stress adaptation pathways. The presence of a few SC genes may reflect lineage-specific modulation or fine-tuning of stress responses in certain accessions.
Cluster 2 (HvNACs high-enriched) and Cluster 6 (upregulated expression after salt stress) are associated with signal transduction, ion homeostasis, and cell wall integrity, all of which are crucial for salt tolerance (Zhu, 2002). Within Cluster 6, we identified two HvNACs, Morex-HvNac55 and Morex-HvNac65, which show high co-expression with multiple salt stress-responsive genes, including regulate signal transduction, ion homeostasis, and cell wall integrity, crucial for salt tolerance. The downregulation of Cluster 8 indicates a possible inhibitory role, potentially preventing excessive activation through modulation of stress responses (Hu et al., 2008). GO enrichment analysis of these clusters revealed that HvNACs are involved in several key stress response pathways, including signal perception and transduction (e.g., receptor activity, ABA signaling), metabolic regulation (carbohydrate homeostasis), ion homeostasis (ion/cation binding), and cellular protective functions (response to stimuli, dehydration, osmotic stress). This finding is consistent with previous studies emphasizing the role of NAC transcription factors in regulating abiotic stress tolerance in various plant species (Nakashima et al., 2009; Tran et al., 2010; Shao et al., 2015). These key HvNACs, especially Morex-HvNac55 and Morex-HvNac65, are strong candidates for molecular breeding and transgenic strategies aimed at enhancing salt tolerance in barley.
Conclusion
This study provides the comprehensive pangenome-scale characterization of the NAC transcription factor family in barley. By integrating analyses of gene content, phylogeny, structural variation, selection pressure, codon usage, and expression, we have unveiled the complex evolutionary history and functional diversification of HvNACs. Our findings indicate that a ranging from 127 to 149 HvNACs are identified across 20 barley accessions, reveal distinct evolutionary trajectories for core versus dispensable genes in HvNACs, underscore the pivotal role of TEs and genomic rearrangements in driving the expansion of HvNACs, and provide strong evidence for the involvement of specific NAC subsets in tissue development and salt stress responses. These findings establish a valuable foundation for future functional genomics studies and may inform efforts to explore the roles of specific HvNACs in stress response and developmental pathways. In addition, this work could support barley improvement programs by identifying candidate genes for enhancing stress tolerance and agronomic performance. Future research directions may include functional validation of key HvNACs using gene editing and transcriptomic approaches, as well as the incorporation of HvNACs diversity into marker-assisted selection to facilitate the development of barley varieties adapted to a range of agro-ecological conditions.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Author contributions
XL: Software, Writing – original draft, Data curation, Conceptualization, Visualization, Writing – review & editing, Formal Analysis. MZ: Data curation, Visualization, Writing – review & editing. JS: Writing – review & editing, Conceptualization, Software. LW: Investigation, Writing – review & editing. MS: Investigation, Writing – review & editing. YZ: Conceptualization, Supervision, Writing – review & editing. QW: Software, Visualization, Writing – review & editing. GC: Supervision, Writing – review & editing, Conceptualization, Project administration, Validation, Methodology.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This study was financially supported by the National Natural Science Foundation of China (32401831 to GC); Youth Taishan Scholar Project of Shandong Province (tsqn202408302 to YZ); Open Foundation of Key Laboratory of Wuliangye-flavor Liquor Solid-state Fermentation, China National Light Industry (2023JJ008); Sichuan Science and Technology Program (2025ZNSFSC1010; MZGC20240023); the Agricultural Scientific and Technological Innovation Project of Shandong Academy of Agricultural Sciences (CXGC2025C01 to GC); The Shandong Postdoctoral Science Foundation (SDCX-ZG-202400117 to CG); and Yibin University Research Projects (2025XJKY027).
Acknowledgments
We thank Jingye Fu from Sichuan Agricultural University for editing the manuscript and providing valuable suggestions.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1635416/full#supplementary-material
References
Barrett, P., Hunter, J., Miller, J. T., Hsu, J. C., and Greenfield, P. (2005). Matplotlib – A Portable Python Plotting Package. In Astronomical Data Analysis Software and Systems XIV, 347, 91.
Bayer, P. E., Golicz, A. A., Scheben, A., Batley, J., and Edwards, D. (2020). Plant pan-genomes are the new reference. Nat. Plants 6, 914–920. doi: 10.1038/s41477-020-0733-0
Bennetzen, J. L. and Wang, H. (2014). The contributions of transposable elements to the structure, function, and evolution of plant genomes. Annu. Rev. Plant Biol. 65, 505–530. doi: 10.1146/annurev-arplant-050213-035811
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Chen, C., Chen, H., Zhang, Y., Thomas, H. R., Frank, M. H., He, Y., et al. (2020). TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol. Plant 13, 1194–1202. doi: 10.1016/j.molp.2020.06.009
Chen, G., Mishina, K., Wang, Q., Zhu, H., Tagiri, A., Kikuchi, S., et al. (2023). Organ-enriched gene expression during floral morphogenesis in wild barley. Plant J. 116, 887–902. doi: 10.1111/tpj.16416
Chen, G., Mishina, K., Zhu, H., Kikuchi, S., Sassa, H., Oono, Y., et al. (2022). Genome-wide analysis of Snf2 gene family reveals potential role in regulation of spike development in barley. Int. J. Mol. Sci. 24, 457. doi: 10.3390/ijms24010457
Chen, Y. and Xia, P. (2025). NAC transcription factors as biological macromolecules responded to abiotic stress: A comprehensive review. Int. J. Biol. Macromol. 242, 142400. doi: 10.1016/j.ijbiomac.2025.142400
Dawson, I. K., Russell, J., Powell, W., Steffenson, B., Thomas, W. T., and Waugh, R. (2015). Barley: a translational model for adaptation to climate change. New Phytol. 206, 913–931. doi: 10.1111/nph.13266
Doebley, J. F., Gaut, B. S., and Smith, B. D. (2006). The molecular genetics of crop domestication. Cell 127, 1309–1321. doi: 10.1016/j.cell.2006.12.006
Fang, Y., Liao, K., Du, H., Xu, Y., Song, H., Li, X., et al. (2015). A stress-responsive NAC transcription factor SNAC3 confers heat and drought tolerance through modulation of reactive oxygen species in rice. J. Exp. Bot. 66, 6803–6817. doi: 10.1093/jxb/erv386
Fuertes-Aguilar, J. and Matilla, A. J. (2024). Transcriptional control of seed life: new insights into the role of the NAC family. Int. J. Mol. Sci. 25, 5369. doi: 10.3390/ijms25105369
Fujita, M., Fujita, Y., Maruyama, K., Seki, M., Hiratsu, K., Ohme-Takagi, M., et al. (2004). A dehydration-induced NAC protein, RD26, is involved in a novel ABA-dependent stress-signaling pathway. Plant J. 39, 863–876. doi: 10.1111/j.1365-313X.2004.02171.x
Grant, C. E., Bailey, T. L., and Noble, W. S. (2011). FIMO: scanning for occurrences of a given motif. Bioinformatics 27, 1017–1018. doi: 10.1093/bioinformatics/btr064
Gu, Z., Eils, R., and Schlesner, M. (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847–2849. doi: 10.1093/bioinformatics/btw313
Guerin, C., Roche, J., Allard, V., Ravel, C., Mouzeyar, S., and Bouzidi, M. F. (2019). Genome-wide analysis, expansion and expression of the NAC family under drought and heat stresses in bread wheat (T. aestivum L.). PloS One 14, e0213390. doi: 10.1371/journal.pone.0213390
Guo, W., Schreiber, M., Marosi, V. B., Bagnaresi, P., Jørgensen, M. E., Braune, K. B., et al. (2025). A barley pan-transcriptome reveals layers of genotype-dependent transcriptional complexity. Nat. Genet. 57, 1–10. doi: 10.1038/s41588-024-02069-y
Hershberg, R. and Petrov, D. A. (2008). Selection on codon bias. Annu. Rev. Genet. 42, 287–299. doi: 10.1146/annurev.genet.42.110807.091442
Hu, R., Qi, G., Kong, Y., Kong, D., Gao, Q., and Zhou, G. (2010). Comprehensive analysis of NAC domain transcription factor gene family in Populus trichocarpa. BMC Plant Biol. 10, 145. doi: 10.1186/1471-2229-10-145
Hu, H., You, J., Fang, Y., Zhu, X., Qi, Z., and Xiong, L. (2008). Characterization of transcription factor gene SNAC2 conferring cold and salt tolerance in rice. Plant Mol. Biol. 67, 169–181. doi: 10.1007/s11103-008-9309-5
Jayakodi, M., Lu, Q., Pidon, H., Rabanus-Wallace, M. T., Bayer, M., Lux, T., et al. (2024). Structural variation in the pangenome of wild and domesticated barley. Nature 636, 1–9. doi: 10.1038/s41586-024-08187-1
Jayakodi, M., Padmarasu, S., Haberer, G., Bonthala, V. S., Gundlach, H., Monat, C., et al. (2020). The barley pan-genome reveals the hidden legacy of mutation breeding. Nature 588, 284–289. doi: 10.1038/s41586-020-2947-8
Jayakodi, M., Schreiber, M., Stein, N., and Mascher, M. (2021). Building pan-genome infrastructures for crop plants and their use in association genetics. DNA Res. 28, dsaa030. doi: 10.1093/dnares/dsaa030
Kim, H. J., Nam, H. G., and Lim, P. O. (2016). Regulatory network of NAC transcription factors in leaf senescence. Curr. Opin. Plant Biol. 33, 48–56. doi: 10.1016/j.pbi.2016.06.002
Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915. doi: 10.1038/s41587-019-0201-4
Kumar, L. and Futschik, M. E. (2007). Mfuzz: a software package for soft clustering of microarray data. Bioinformation 2, 5–7. doi: 10.6026/97320630002005
Li, W. and Godzik, A. (2006). Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659. doi: 10.1093/bioinformatics/btl158
Li, W., Han, L., Qian, Y., and Sun, Z. (2011). Characteristics and functions of NAC transcription factors in plants. Chin. J. Appl. Environ. Biol. 17, 596–606.
Li, W., Li, X., Chao, J., Zhang, Z., Wang, W., and Guo, Y. (2018). NAC family transcription factors in tobacco and their potential role in regulating leaf senescence. Front. Plant Sci. 9, 1900. doi: 10.3389/fpls.2018.01900
Mascher, M., Wicker, T., Jenkins, J., Plott, C., Lux, T., Koh, C. S., et al. (2021). Long-read sequence assembly: a technical evaluation in barley. Plant Cell 33, 1888–1906. doi: 10.1093/plcell/koab077
Mathew, I. E. and Agarwal, P. (2018). May the fittest protein evolve: favoring the plant-specific origin and expansion of NAC transcription factors. BioEssays 40, 1800018. doi: 10.1002/bies.201800018
Mohanta, T. K., Yadav, D., Khan, A., Hashem, A., Tabassum, B., Khan, A. L., et al. (2020). Genomics, molecular and evolutionary perspective of NAC transcription factors. PloS One 15, e0231425. doi: 10.1371/journal.pone.0231425
Nakashima, K., Ito, Y., and Yamaguchi-Shinozaki, K. (2009). Transcriptional regulatory networks in response to abiotic stresses in Arabidopsis and grasses. Plant Physiol. 149, 88–95. doi: 10.1104/pp.108.129791
Nguyen, L.-T., Schmidt, H. A., Von Haeseler, A., and Minh, B. Q. (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274. doi: 10.1093/molbev/msu300
Nuruzzaman, M., Manimekalai, R., Sharoni, A. M., Satoh, K., Kondoh, H., Ooka, H., et al. (2010). Genome-wide analysis of NAC transcription factor family in rice. Gene 465, 30–44. doi: 10.1016/j.gene.2010.06.008
Nuruzzaman, M., Sharoni, A. M., and Kikuchi, S. (2013). Roles of NAC transcription factors in the regulation of biotic and abiotic stress responses in plants. Front. Microbiol. 4, 248. doi: 10.3389/fmicb.2013.00248
Olsen, A. N., Ernst, H. A., Leggio, L. L., and Skriver, K. (2005). NAC transcription factors: structurally distinct, functionally diverse. Trends Plant Sci. 10, 79–87. doi: 10.1016/j.tplants.2004.12.010
Ooka, H., Satoh, K., Doi, K., Nagata, T., Otomo, Y., Murakami, K., et al. (2003). Comprehensive analysis of NAC family genes in Oryza sativa and Arabidopsis thaliana. DNA Res. 10, 239–247. doi: 10.1093/dnares/10.6.239
Panchy, N., Lehti-Shiu, M., and Shiu, S.-H. (2016). Evolution of gene duplication in plants. Plant Physiol. 171, 2294–2316. doi: 10.1104/pp.16.00523
Parankusam, S., Bhatnagar-Mathur, P., and Sharma, K. K. (2017). Heat responsive proteome changes reveal molecular mechanisms underlying heat tolerance in chickpea. Environ. Exp. Bot. 141, 132–144. doi: 10.1016/j.envexpbot.2017.07.007
Peng, X., Zhao, Y., Li, X., Wu, M., Chai, W., Sheng, L., et al. (2015). Genomewide identification, classification and analysis of NAC type gene family in maize. J. Genet. 94, 377–390. doi: 10.1007/s12041-015-0526-9
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T.-C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Pourkheirandish, M. and Komatsuda, T. (2007). The importance of barley genetics and domestication in a global perspective. Ann. Bot. 100, 999–1008. doi: 10.1093/aob/mcm139
Puranik, S., Sahu, P. P., Mandal, S. N., Parida, S. K., and Prasad, M. (2013). Comprehensive genome-wide survey, genomic constitution and expression profiling of the NAC transcription factor family in foxtail millet (Setaria italica L.). PloS One 8, e64594. doi: 10.1371/journal.pone.0064594
Puranik, S., Sahu, P. P., Srivastava, P. S., and Prasad, M. (2012). NAC proteins: regulation and role in stress tolerance. Trends Plant Sci. 17, 369–381. doi: 10.1016/j.tplants.2012.02.004
Russell, J., Mascher, M., Dawson, I. K., Kyriakidis, S., Calixto, C., Freund, F., et al. (2016). Exome sequencing of geographically diverse barley landraces and wild relatives gives insights into environmental adaptation. Nat. Genet. 48, 1024–1030. doi: 10.1038/ng.3612
Shao, H., Wang, H., and Tang, X. (2015). NAC transcription factors in plant multiple abiotic stress responses: progress and prospects. Front. Plant Sci. 6, 902. doi: 10.3389/fpls.2015.00902
Sharp, P. M., Tuohy, T. M., and Mosurski, K. R. (1986). Codon usage in yeast: cluster analysis clearly differentiates highly and lowly expressed genes. Nucleic Acids Res. 14, 5125–5143. doi: 10.1093/nar/14.13.5125
Singh, S., Koyama, H., Bhati, K. K., and Alok, A. (2021). The biotechnological importance of the plant-specific NAC transcription factor family in crop improvement. J. Plant Res. 134, 475–495. doi: 10.1007/s10265-021-01270-y
Su, G., Morris, J. H., Demchak, B., and Bader, G. D. (2014). Biological network exploration with Cytoscape 3. Curr. Protoc. Bioinf. 47, 8–13. doi: 10.1002/0471250953.bi0813s47
Sun, L.-J., Li, D.-Y., Zhang, H.-J., and Song, F.-M. (2012). Functions of NAC transcription factors in biotic and abiotic stress responses in plants. Hereditas 34, 993–1002. doi: 10.3724/sp.j.1005.2012.00993
Tettelin, H., Masignani, V., Cieslewicz, M. J., Donati, C., Medini, D., Ward, N. L., et al. (2005). Genome analysis of multiple pathogenic isolates of streptococcus agalactiae: implications for the microbial “pan-genome. Proc. Natl. Acad. Sci. 102, 13950–13955. doi: 10.1073/pnas.0506758102
Tong, C., Jia, Y., Hu, H., Zeng, Z., Chapman, B., Li, C., et al. (2025). Pangenome and pantranscriptome as the new reference for gene-family characterization: A case study of basic helix-loop-helix (bHLH) genes in barley. Plant Communications, 6, 101190. doi: 10.1016/j.xplc.2024.101190
Tran, L.-S. P., Nakashima, K., Sakuma, Y., Simpson, S. D., Fujita, Y., Maruyama, K., et al. (2004). Isolation and functional analysis of Arabidopsis stress-inducible NAC transcription factors that bind to a drought-responsive cis-element in the early responsive to dehydration stress 1 promoter. Plant Cell 16, 2481–2498. doi: 10.1105/tpc.104.022699
Tran, L.-S. P., Nishiyama, R., Yamaguchi-Shinozaki, K., and Shinozaki, K. (2010). Potential utilization of NAC transcription factors to enhance abiotic stress tolerance in plants by biotechnological approach. GM Crops 1, 32–39. doi: 10.4161/gmcr.1.1.10569
Villanueva, R. A. M. and Chen, Z. J. (2019). ggplot2: elegant graphics for data analysis, 2nd edition. Meas. Interdiscip. Res. 17, 160–167. doi: 10.1080/15366367.2019.1565254
Wang, Y., Li, J., and Paterson, A. H. (2013). MCScanX-transposed: detecting transposed gene duplications based on multiple colinearity scans. Bioinformatics 29, 1458–1460. doi: 10.1093/bioinformatics/btt150
Wicker, T., Gundlach, H., Spannagl, M., Uauy, C., Borrill, P., Ramírez-González, R. H., et al Impact of transposable elements on genome structure and evolution in wheat. (2018) Genome Biol. 19, 103. doi: 10.1186/s13059-018-1479-0
Wicker, T., Sabot, F., Hua-Van, A., Bennetzen, J. L., Capy, P., Chalhoub, B., et al. (2007). A unified classification system for eukaryotic transposable elements. Nat. Rev. Genet. 8, 973–982. doi: 10.1038/nrg2165
Xie, Q., Sanz-Burgos, A. P., Guo, H., García, J. A., and Gutiérrez, C. (1999). GRAB proteins, novel members of the NAC domain family, isolated by their interaction with a geminivirus protein. Plant Mol. Biol. 39, 647–656. doi: 10.1023/A:1006138221874
Xiong, H., He, H., Chang, Y., Miao, B., Liu, Z., Wang, Q., et al. (2025). Multiple roles of NAC transcription factors in plant development and stress responses. J. Integr. Plant Biol. 67, 510–538. doi: 10.1111/jipb.13854
Xu, Z.-Y., Kim, S. Y., Hyeon, D. Y., Kim, D. H., Dong, T., Park, Y., et al. (2013). The Arabidopsis NAC transcription factor ANAC096 cooperates with bZIP-type transcription factors in dehydration and osmotic stress responses. Plant Cell 25, 4708–4724. doi: 10.1105/tpc.113.119099
Yuan, X., Wang, H., Cai, J., Li, D., and Song, F. (2019). NAC transcription factors in plant immunity. Phytopathol. Res. 1, 1–13. doi: 10.1186/s42483-018-0008-0
Zhang, Z., Li, J., Zhao, X.-Q., Wang, J., Wong, G. K.-S., and Yu, J. (2006). KaKs_Calculator: calculating Ka and Ks through model selection and model averaging. Genomics. Proteomics Bioinf. 4, 259–263. doi: 10.1016/S1672-0229(07)60007-2
Zhong, R., Lee, C., and Ye, Z.-H. (2010). Evolutionary conservation of the transcriptional network regulating secondary cell wall biosynthesis. Trends Plant Sci. 15, 625–632. doi: 10.1016/j.tplants.2010.08.007
Zhu, J.-K. (2002). Salt and drought stress signal transduction in plants. Annu. Rev. Plant Biol. 53, 247–273. doi: 10.1146/annurev.arplant.53.091401.143329
Keywords: Hordeum vulgare, NAC, pangenome, evolutionary dynamics, pantranscriptome, salt stress response
Citation: Liu X, Zhang M, Su J, Wu L, Shen M, Zhuang Y, Wang Q and Chen G (2025) The evolution, variation, and expression patterns under development and stress responses of the NAC gene family in the barley pan-genome. Front. Plant Sci. 16:1635416. doi: 10.3389/fpls.2025.1635416
Received: 26 May 2025; Accepted: 14 July 2025;
Published: 07 August 2025.
Edited by:
Prateek Gupta, SRM University, IndiaReviewed by:
Rezwan Tariq, Texas A&M University, United StatesSantisree Parankusam, International Crops Research Institute for the Semi-Arid Tropics (ICRISAT), India
Yong Jia, Murdoch University, Australia
Copyright © 2025 Liu, Zhang, Su, Wu, Shen, Zhuang, Wang and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Gang Chen, Y2cxOTkxMDkwNUBnbWFpbC5jb20=