Triticum monococcum L. subsp. aegilopoides, accession G52, an ancient and resilient crop with considerable potential for improving stress tolerance in modern wheat, was employed to investigate the role of ADF genes in response to cold stress. Seeds were kindly provided by Professor Lianquan Zhang (Sichuan Agricultural University). To break dormancy, seeds were stratified at 4°C for 24 hours, followed by germination at room temperature for 7 days. Subsequently, ten uniform seedlings were subjected to cold treatment at 4°C for 24 hours, while another set of ten seedlings was maintained under control conditions. Leaf tissues from both treated and control groups were collected, with three biological replicates per group. Samples were immediately flash-frozen in liquid nitrogen and stored at −80°C for subsequent analyses.

The reference genome of Triticum monococcum L. subsp. aegilopoides TA299 was obtained from the DRYAD repository (https://doi.org/10.5061/dryad.v41ns1rxj). The ADF domain (Pfam accession: PF00241) was retrieved from the Pfam database (http://pfam.xfam.org/). A Hidden Markov Model (HMM) search was performed using HMMER v3.0 (Mistry et al., 2013) with a stringent E-value threshold of ≤ 1e−10 to identify candidate ADF protein sequences. The presence of the ADF domain in these sequences was further validated using the NCBI Conserved Domain Database (CDD; https://www.ncbi.nlm.nih.gov/cdd/), SMART (http://smart.embl-heidelberg.de/), and Pfam. Physicochemical properties of the identified TbADFs, including molecular weight, isoelectric point, and instability index, were calculated using the ExPASy ProtParam tool (https://web.expasy.org/protparam/). Subcellular localization of TbADFs were predicted using WoLF PSORT (https://wolfpsort.hgc.jp/) and Cell-PLoc 2.0(http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/).

Positional information for ADF genes was extracted from the GFF3 annotation file of the Triticum monococcum L. subsp. aegilopoides reference genome. Chromosomal locations of the identified ADF genes were visualized using TBtools v2.097 (Chen et al., 2023). A phylogenetic tree was constructed based on full-length ADF protein sequences from eight plant species—Arabidopsis thaliana, Aegilops tauschii, Oryza sativa, Zea mays, Triticum urartu, Triticum turgidum, Triticum aestivum, and Triticum monococcum L. subsp. aegilopoides—using the Maximum Likelihood (ML) method in IQ-TREE v2.0 (Nguyen et al., 2015). The best-fit amino acid substitution model was determined automatically using ModelFinder within IQ-TREE, and tree robustness was assessed with 1,000 bootstrap replicates. To investigate evolutionary relationships and genomic conservation, collinearity analysis between Triticum monococcum L. subsp. aegilopoides and three related wheat species (T. urartu, T. turgidum, and T. aestivum) was performed using MCScanX (Wang et al., 2012) with default parameters.

Gene structure information for the ADF genes was extracted from the GFF3 annotation file of the Triticum monococcum L. subsp. aegilopoides reference genome. Conserved motifs within the TbADFs protein sequences were identified using the MEME Suite (http://meme-suite.org/meme), with the maximum number of motifs set to 10 and other parameters kept at default settings. To investigate potential regulatory mechanisms, the 2,000 bp upstream regions from the transcription start sites of the TbADFs were extracted as putative promoter sequences. These promoter regions were analyzed for cis-acting regulatory elements using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

The tertiary structure of the TbADFs were predicted using AlphaFold2. The predicted 3D structure was subsequently visualized and analyzed using PyMOL v2.5.5 (DeLano, 2002), facilitating detailed examination of the protein’s spatial conformation. To further investigate the surface topology and potential functional sites, the CASTp 3.0 web server (http://sts.bioe.uic.edu/castp/index.html?2pk9) (Tian et al., 2018) was utilized to identify and characterize surface pockets and cavities at multiple scales, offering insights into possible ligand-binding regions and structural features relevant to TbADFs.

Codon usage bias was assessed using CodonW (version 1.4.2; http://codonw.sourceforge.net/), with several key indices calculated to characterize codon usage patterns. These included the effective number of codons (ENC), codon adaptation index (CAI), relative synonymous codon usage (RSCU), overall genomic GC content, GC content at the third codon position (GC3s), and nucleotide frequencies at the third synonymous codon position (T3s, C3s, A3s, G3s). To further investigate the factors shaping codon usage bias, Parity Rule 2 (PR2) analysis was performed to explore the relative contributions of mutational pressure and natural selection. Additionally, ENC-GC3s plots were generated to evaluate the extent of codon bias across genes. All visualizations, including the ENC and PR2 plots, were created using the Matplotlib library (Barrett et al., 2005) to ensure clear graphical representation of the results.

The TbADFs amino acid sequences from Triticum monococcum L. subsp. aegilopoides were submitted to the STRING database (http://www.string-db.org/) (Szklarczyk et al., 2021) to construct a protein-protein interaction (PPI) network, with A. thaliana selected as the reference organism due to its well-characterized protein interaction data. The analysis was performed using default parameters, including a medium confidence score threshold (0.6) to predict functional associations based on evidence such as co-expression, experimental data, and conserved co-occurrence. The resulting PPI network was visualized to identify potential interacting partners of TbADFs and to infer their functional roles in cellular processes.

To evaluate DNA methylation levels in TbADFs, we quantified methylation frequencies in both promoter and gene body regions using publicly available whole-genome bisulfite sequencing data for Triticum monococcum L. subsp. aegilopoides accession TA299, obtained from the Dryad Digital Repository (https://datadryad.org/dataset/doi:10.5061/dryad.v41ns1rxj). Promoter regions were defined as the 2 kilobases upstream of the transcription start site (TSS), while gene body regions included the entire transcribed sequence excluding the promoter.

For each TbADFs and genomic region, methylation frequency was calculated as the ratio of methylated CpG sites to the total number of CpG sites (methylated + unmethylated), using the following formula:

Methylation calls were extracted from the downloaded data files, and region-specific methylation frequencies were computed using custom Python scripts. These scripts parsed CpG methylation count data and calculated the average methylation level across all CpG sites within each defined region for each gene.

Transcriptome data for Triticum monococcum L. subsp. aegilopoides from various tissues and cold stress treatments were retrieved from the NCBI Sequence Read Archive (SRA) under accession numbers PRJEB61155 and PRJEB21284. Adapter sequences and low-quality reads were removed using Trimmomatic (v0.39) (Bolger et al., 2014). The resulting high-quality reads were pseudo-aligned to the reference genome and quantified using Kallisto (v0.46.2) (Bray et al., 2016) to estimate gene expression levels. A heatmap was generated using TBtools (v2.097) (Chen et al., 2023) to visualize expression patterns. To explore gene co-expression relationships, Pearson correlation analysis was performed between four upregulated TbADFs and all other expressed genes. Raw p-values from these pairwise correlation tests were adjusted for multiple comparisons using the Benjamini-Hochberg procedure to control the false discovery rate (FDR) at 5%. Genes with a correlation coefficient greater than 0.9 and a p-value less than 0.001 were first identified as candidate highly co-expressed, these were subsequently subjected to Benjamini-Hochberg correction, and pairs with an FDR-adjusted p-value less than 0.05 were considered highly co-expressed. Co-expression networks and conserved motif distributions were visualized using Cytoscape (v3.10.2) (Shannon et al., 2003). Gene Ontology (GO) enrichment analysis was conducted to infer the potential biological functions of the co-expressed genes, with GO annotations obtained from the KOBAS database (http://bioinfo.org/kobas) (Xie et al., 2011). The enrichment results were visualized using appropriate R packages.

Total RNA was extracted from Triticum monococcum L. subsp. aegilopoides using the RNAprep Pure Plant Kit (Tiangen, China), following the manufacturer’s instructions. First-strand cDNA was synthesized from 1 µg of total RNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA). Gene-specific primers for the TbADF gene family ( Supplementary Table S1 ) were designed using Primer 5 software. Quantitative real-time PCR (qRT-PCR) was conducted on a Veriti 96-Well Thermal Cycler (Applied Biosystems, USA), with actin used as the internal reference gene for normalization. Each 10 µl reaction contained 5 µl of 2× TB Green Premix Ex Taq II (Tli RNaseH Plus; Takara, Japan), 0.4 µl of each primer (10 µM), 1 µl of cDNA, and 3.2 µl of sterile double-distilled water. The amplification protocol consisted of an initial denaturation at 95°C for 3 min, followed by 39 cycles of 95°C for 10 s and 58°C for 30 s. A melt curve analysis was performed to verify amplification specificity, using a ramp from 65°C to 95°C in 0.5°C increments after denaturation at 95°C for 5 s. Relative expression levels of TbADFs were calculated using the 2^−ΔΔCt method.

Clone the CDS sequence of TbADF6 from Triticum monococcum L. subsp. aegilopoides cDNA. Digest the pCAMBIA2300-GFP vector with BamHI, and insert the amplified product into the linearized pCAMBIA2300-GFP vector using the pEASY^®-Basic Seamless Cloning and Assembly Kit. Verify the recombinant construct by DNA sequencing. Transform the construct into Agrobacterium. Resuspend and mix the Agrobacterium cultures, then infiltrate them into tobacco leaves. After 36 hours of incubation, observe the fluorescence using a confocal laser scanning microscope.

To identify ADF genes in Triticum monococcum L. subsp. aegilopoides, BLASTP searches were performed using known ADF protein sequences from A. thaliana, O. sativa, and Zea mays as queries against the Triticum monococcum L. subsp. aegilopoides. Redundant sequences were filtered using HMM profiling, resulting in the identification of nine unique ADF proteins, designated TbADF1 through TbADF9 ( Supplementary Table S2 ). These proteins ranged in length from 138 to 270 amino acids, with molecular weights between 16.08 and 30.19 kDa and isoelectric points (pI) from 4.81 to 9.44. Six of the proteins were acidic (pI< 7). The instability index ranged from 33.17 to 56.08, indicating that six proteins may be unstable (instability index > 40). The aliphatic index varied from 58.17 to 80.52, and all nine proteins had negative GRAVY (grand average of hydropathicity) scores, suggesting a hydrophilic nature. Subcellular localization predictions indicated that all TbADFs are cytoplasmic.

Chromosomal mapping revealed that the nine TbADFs are unevenly distributed across five chromosomes. Specifically, chromosomes 2A, 4A, and 6A each contain one gene, while chromosome 1A has two genes, and chromosome 5A carries four ( Figure 1A ). Notably, TbADF6, TbADF7, and TbADF8 are clustered in adjacent genomic regions, suggesting possible tandem duplication events or localized selective pressures. To assess the evolutionary relationships among TbADFs, a phylogenetic tree was constructed using 97 ADF protein sequences from eight plant species. These ADFs clustered into four subfamilies (subgroups A–D), each containing 13 to 46 members. The TbADFs are represented in all four subgroups: one gene in subgroup A, two genes each in subgroups B and C, and four genes in subgroup D ( Figure 1B ).

To further explore ADF gene conservation, synteny analysis was performed between Triticum monococcum L. subsp. aegilopoides and three wheat relatives of varying ploidy levels: diploid Triticum urartu, tetraploid Triticum turgidum, and hexaploid Triticum aestivum. This analysis identified 3, 12, and 16 orthologous gene pairs, respectively ( Figures 2A–C ), suggesting that ADF genes have been conserved through wheat evolution and polyploidization, likely under purifying selection.

Gene structure analysis revealed a relatively conserved exon-intron architecture among the nine TbADFs. Three genes contained two exons, five had three exons, and TbADF7 exhibited the most complex structure with five exons ( Figure 3A ). This conserved structural architecture suggests evolutionary stability within the ADF gene family. Motif analysis using the MEME Suite identified four conserved motifs shared by all TbADFs, predominantly located within the canonical ADF domain ( Figure 3A ). Multiple sequence alignment and tertiary structure modeling further confirmed the conserved positioning of the ADF domain across all family members ( Figure 3B ). Structural predictions generated by AlphaFold2 and visualized with PyMOL reinforced this conserved spatial configuration, highlighting the structural integrity of the ADF domain among TbADFs ( Figure 3C ).

To explore regulatory potential, the 2,000 bp upstream promoter regions of TbADFs were analyzed using PlantCARE. Numerous cis-acting regulatory elements were detected, including those associated with growth, hormone signaling, and stress response ( Figure 4 ; Supplementary Table S3 ). Notable elements included hormone-responsive motifs such as the abscisic acid response element (ABRE), jasmonic acid-responsive motifs (CGTCA-motif and TGACG-motif), auxin-responsive TGA-element, and gibberellin-responsive P-box. Additionally, stress-associated elements such as the low-temperature response element (LTR), drought-responsive DRE, and hypoxia-responsive motifs (GC-motif, ARE) were also identified. Collectively, these findings suggest that TbADFs are not only structurally conserved but are also likely involved in complex regulatory networks mediating hormonal and environmental stress responses. This highlights their potential roles in enhancing wheat adaptability to adverse environmental conditions.

The detection of protein binding pockets is critical for understanding structural diversity and interaction specificity, particularly in predicting functional sites such as ligand-binding regions. To explore functional diversity among TbADFs, potential catalytic and ligand-binding sites were predicted using CASTp 3.0. Red-highlighted regions indicated the predicted binding pockets ( Figure 5 ), and their surface areas and volumes were calculated ( Supplementary Table S4 ). These binding sites are surrounded by key amino acid residues such as tyrosine (Tyr), phenylalanine (Phe), valine (Val), asparagine (Asn), lysine (Lys), alanine (Ala), aspartic acid (Asp), glutamic acid (Glu), isoleucine (Ile), and leucine (Leu), which are likely involved in mediating protein–ligand interactions and catalytic activity. The conservation of these residues supports their functional importance and potential roles in the stress-responsive behavior of ADFs.

To elucidate the protein interaction network of the ADF gene family in Triticum monococcum L. subsp. aegilopoides, nine TbADF protein sequences were analyzed using the STRING database with A. thaliana as the reference organism. ( Figure 6A ). These interacting proteins are implicated in critical cellular processes, including cytoskeletal dynamics, cell elongation, shape modulation, division, and migration, which collectively influence transpiration regulation, pathogen defense responses, and overall plant growth and development. The high connectivity of TbADFs with these partners, supported by a confidence score threshold of 0.7, suggests a conserved functional module that may underpin stress adaptation and developmental plasticity in wild einkorn wheat. This interaction network provides a foundation for future functional studies.

To further explore transcriptional regulation, promoter regions of the TbADFs genes were scanned for TF binding motifs using the FIMO ( Figure 6B ). This analysis revealed a complex regulatory network characterized by diverse TF interactions. Notably, members of the ERF family exhibited the highest binding site frequency, with a total of 2516 sites identified across multiple TbADFs, particularly in TbADF1 and TbADF7. This suggests that ERF transcription factors may play a dominant role in modulating ADF gene expression, likely in response to abiotic stress signals, given the GC-rich nature of their target motifs. In addition, several other TF families, including LBD, MYB, TCP, and C2H2, also contribute to the regulatory architecture. Although these families are less abundant than ERFs, their presence suggests auxiliary roles in fine-tuning ADF gene expression under specific developmental or stress-related conditions.

The codon usage bias of the ADF gene family in Triticum monococcum L. subsp. aegilopoides was analyzed to explore the evolutionary pressures shaping gene sequence composition. The ENC plot ( Figure 7A ) showed that ADF genes exhibit ENC values ranging from 39.7 to 52.3. Most values fall below the expected standard curve under neutral evolution, indicating a moderate codon usage bias predominantly influenced by GC content at the GC3s, which ranged from 0.644 to 0.829 ( Supplementary Table S5 ).

The PR2-bias plot ( Figure 7B ) demonstrated an approximately balanced distribution between A3 and T3, and between G3 and C3, with most data points clustering near the central point (0.5, 0.5). This pattern suggests that both mutational pressure and natural selection contribute to shaping codon usage in TbADFs.

A heatmap of RSCU values ( Figure 7C ) further highlighted a preference for G/C-ending codons across the ADF family, consistent with the GC-rich genomic context of T. monococcum subsp. aegilopoides. Notably, codons such as GGC (Glycine) and UGC (Cysteine) were overrepresented, implying potential optimization for translation efficiency. This bias may be linked to tRNA abundance or a selective advantage for efficient protein synthesis, particularly under stress conditions, aligning with the known role of ADF proteins in cytoskeletal dynamics and stress response.

This analysis was based on high-quality, genome-wide methylation data for Triticum monococcum L. subsp. aegilopoides accession TA299, which we obtained from Ahmed et al. (2023). The final dataset comprised 221,159,329 CpG sites with a total sequencing depth of 19,524,541,735 and an average coverage of 88.28× per site. Using this high-coverage dataset, we calculated average methylation frequencies for all CpG sites within each defined gene region. The analysis of CpG methylation frequencies in the promoter and gene body regions of nine ADF genes (TbADF1 to TbADF9) revealed significant variation in methylation levels ( Supplementary Table S6 ), with promoter methylation ranging from 0.0907 (TbADF2) to 0.3053 (TbADF9) and gene body methylation from 0.0797 (TbADF2) to 0.2740 (TbADF1). The mean promoter methylation (0.1923) was slightly higher than the mean gene body methylation (0.1823), suggesting a potential regulatory role of promoter methylation in gene expression. Notably, TbADF9 exhibited the highest methylation levels in both regions (0.3053 and 0.2644, respectively), which may indicate a repressive epigenetic state, while TbADF2 displayed the lowest levels (0.0907 and 0.0797).

Transcriptomic analysis across six tissues of Triticum monococcum L. subsp. aegilopoides revealed distinct tissue-specific expression patterns among the nine TbADFs ( Figure 8A ). TbADF2 and TbADF7 exhibited predominant expression in the glumes, whereas TbADF3 and TbADF9 were highly expressed in the grains, with low expression in other tissues. The remaining genes showed broader expression profiles: TbADF1 was mainly expressed in both flag leaves and grains; TbADF4 displayed elevated expression in flag leaves, glumes, and grains; and TbADF5 and TbADF8 were primarily expressed in flag leaves, glumes, and roots. These expression patterns indicate that TbADFs are functionally diversified, contributing to tissue-specific roles during the growth and development of wild einkorn wheat.

Transcriptomic profiling following a 24-hour cold treatment at 4°C revealed distinct differential expression patterns among the TbADFs. Notably, TbADF2, TbADF3, and TbADF9 exhibited no detectable expression under either control or cold-stress conditions, suggesting their non-involvement in the immediate cold response ( Figure 8B ). In contrast, the remaining six ADF genes were cold-responsive: TbADF5 and TbADF8 showed downregulation, whereas TbADF1, TbADF4, TbADF6, and TbADF7 were significantly upregulated. Among them, TbADF6 displayed the most dramatic induction, with TPM increasing from 29.07 to 300.01, highlighting its potential central role in cold stress adaptation. To validate these transcriptomic findings, expression analysis of TbADF1, TbADF4, TbADF6, and TbADF7 was conducted, confirming consistent upregulation of all four genes under cold stress conditions ( Figure 8C ). This concordance between transcriptome data and expression validation strongly supports their involvement in the cold response of wild einkorn wheat.

To further dissect the molecular mechanisms underlying this response, we focused on the four upregulated genes, TbADF1, TbADF4, TbADF6, and TbADF7. Pearson correlation analysis (r > 0.9, p< 0.001) identified strongly co-expressed genes under cold stress: 103 for TbADF1, 1,269 for TbADF4, 1,706 for TbADF6, and 1,657 for TbADF7 ( Figure 9A ). To elucidate the biological relevance of these co-expression networks, we performed GO enrichment analysis ( Figure 9B ).

In the MF category, terms such as purine ribonucleoside triphosphate binding (GO:0035639) and phosphoprotein phosphatase activity (GO:0004721) were significantly enriched, implicating roles in energy metabolism and phosphorylation-dependent signal transduction during cold stress. In the CC category, enrichments for membrane (GO:0016020) and organelle outer membrane (GO:0031968) suggest localization of co-expressed proteins to structural membrane sites, possibly mediating signal transmission and maintaining cellular integrity under stress. BP terms such as maturation of LSU-rRNA (GO:0000470) and intracellular chemical homeostasis (GO:0055082) point to enhanced ribosome biogenesis and metabolic homeostasis, key mechanisms for sustaining protein synthesis and stability during stress. Furthermore, enrichments in response to fatty acid (GO:0070542) and cellular response to endogenous stimulus (GO:0071495) underscore the involvement of lipid signaling and hormone-regulated pathways, both critical for cold adaptation through modulation of membrane fluidity and internal signaling cascades. Collectively, these findings highlight the intricate regulatory and functional roles of cold-inducible ADF genes, particularly TbADF6, in orchestrating a coordinated molecular response to low-temperature stress in wild einkorn wheat.

To investigate the subcellular localization of TbADF6, a GFP-tagged fusion construct was generated and transiently expressed in tobacco leaves via Agrobacterium-mediated transformation. Confocal laser scanning microscopy showed that the GFP-ADF fusion protein was mainly localized in the nucleus and cytoplasm ( Figure 10 ), which is consistent with the subcellular localization prediction ( Supplementary Table S2 ). These results suggest that TbADF6 protein may function in both compartments.