1. Introduction

Brachypodium (L.) P. Beauv. is a small, cool-climate genus of about 23 species distributed worldwide and belongs to the monotypic tribe Brachypodieae (Catalán et al., 2016, 2023). Phylogenetically, Brachypodium occupies an intermediate position between the ancestral Pooideae lineages and the recently evolved core pooid clade (Catalán et al., 2016). Twenty of its species are perennial, while the three best-investigated annual taxa, the diploid species B. distachyon and B. stacei and their allotetraploid derivative, B. hybridum, were selected as a functional model system for the study of allopolyploidy (Scholthof et al., 2018; Gordon et al., 2020). Approximately half of Brachypodium species are diploid, while the remaining ones are polyploid, ranging from tetraploid to hexaploid (Díaz-Pérez et al., 2018; Sancho et al., 2022; Decena et al., 2024). Brachypodium displays a remarkable descending dysploidy, with a karyotypic evolution from an ancestral karyotype of x = 10 chromosomes to a more recent one of x = 9 chromosomes, and distinct reductions to different karyotypes of x = 8 and x = 5 chromosomes, some of which are only found in orphan subgenomes of polyploid species (Lusinska et al., 2019; Sancho et al., 2022).

Brachypodium sylvaticum, a perennial Palearctic grass, is one of the most widely distributed species of the genus Brachypodium. B. sylvaticum is native to Europe, Asia, and North Africa (Catalán et al., 2023) and was introduced in western North America, becoming an invasive species (Rosenthal et al., 2008; Roy et al., 2011). This grass has gained attention in plant research as a model species for perennial grasses due to its adaptability and wide distribution across various habitats (Carroll and Somerville, 2009; Dohleman and Long, 2009; Glover et al., 2010). Its significance as a model plant for perenniality-to-annuality transitions has been underscored by extensive genomic (reference genomes Ain1 and Sin1) and biological resources (Steinwand et al., 2013; Gordon et al., 2015; Gordon et al. 2016; Lei et al., 2024; Phytozome v.13, https://phytozome-next.jgi.doe.gov/).

Brachypodium sylvaticum is characterized by its short and slender rhizomes, nodding panicles, densely hairy habit, and long-awned lemmas (Schippmann, 1991). The species is self-compatible and thrives in mesic and humid forest habitats (Khan and Stace, 1999; Steinwand et al., 2013). Alongside typical B. sylvaticum, four microtaxa (B. breviglume, B. kurilense, B. spryginii, B. miserum) have been identified as close relatives, sharing certain morphological characteristics and distribution areas with B. sylvaticum sensu stricto (s.s.) (Tzvelev and Probatova, 2019; Catalán et al., 2023). The entire B. sylvaticum complex is constituted by diploid species, with a consistent chromosome number of 2n = 2x = 18 (x = 9) (Catalán et al., 2016; Díaz-Pérez et al., 2018). Another closely related taxon is the perennial Mediterranean B. glaucovirens, also diploid but with a different chromosome number of 2n = 2x = 16 (x = 8). B. glaucovirens exhibits overall features similar to those of B. sylvaticum (short rhizome, large awns) and B. pinnatum (bright green leaves, broad leaf ribs, erect panicle) (Schippmann, 1991). Previous studies suggested that B. glaucovirens may have a hybrid origin (Scholz, 2007), and plastome-based phylogenies confirmed that B. sylvaticum s.s. served as the maternal progenitor for this taxon (Catalán et al., 2023). Despite its wide distribution and significance as a model organism, the B. sylvaticum complex remains an active area of research (Keng, 1982; Tzvelev, 2015; Catalán et al., 2023). Some of these microtaxa exhibit large disjunct geographic distributions in their native Palearctic area, suggesting genetic isolation and evolutionary divergence, warranting further phylogenetic and genomic investigation.

Phylogenetic studies based on nuclear and plastid loci consistently positioned the Brachypodium sylvaticum complex within the recently evolved core-perennial clade of Brachypodium (Catalán and Olmstead, 2000; Catalán et al., 2016; Díaz-Pérez et al., 2018). Molecular dating suggested that the B. sylvaticum s.s. lineage originated in the Mid-Late Pleistocene (Díaz-Pérez et al., 2018), coinciding with major climatic oscillations that likely influenced its evolutionary trajectory. Comparative chromosome barcoding has further detailed the karyotypic profiles of B. sylvaticum and its close congeners, providing insights into their intermediate karyotype evolutionary history and chromosome dynamics (Lusinska et al., 2019; Sancho et al., 2022).

Brachypodium sylvaticum maintains a ubiquitous and highly stable symbiosis with its Epichloë sylvatica fungal endophyte (Leuchtmann and Schardl, 2022). This association is predominantly characterized by vertical transmission of the endophyte through the host seeds, ensuring the systemic colonization of new generations of host plants without visible symptoms (Meijer and Leuchtmann, 1999). However, E. sylvatica also retains the capacity for sexual reproduction via stromata formation and horizontal transmission, which can inhibit host flowering, though such events are rare and typically confined to specific microhabitats (Brem and Leuchtmann, 2001). This dual transmission strategy underscores the evolutionary flexibility of E. sylvatica, oscillating between predominantly mutualistic and, under certain conditions, antagonistic interactions with its host grass.

The ecological advantages of this mutualistic symbiosis are of great importance for the holobiont. Experimental feeding trials have demonstrated that spontaneous Epichloë sylvatica infection enhances herbivore resistance in Brachypodium sylvaticum (Brem and Leuchtmann, 2001). While the precise chemical basis of deterrence to herbivores remains unclear, Epichloë endophytes in related grasses are known to influence host physiology through metabolic alterations or secondary metabolite (alkaloids) production (Saikkonen et al., 2016; Gundel et al., 2017; Vikuk et al., 2019). Moreover, Epichloë endophytes in other grass species confer additional adaptive advantages, such as increased drought tolerance and pathogen resistance, raising the possibility that B. sylvaticum benefits similarly (Vázquez de Aldana et al., 2015). Given the nearly ubiquitous presence of E. sylvatica in natural populations of B. sylvaticum (100% of tested individuals infected) (Leuchtmann and Schardl, 2022), its role as a key ecological and evolutionary factor in the phylogeography and adaptive strategies of its host species should not be overlooked.

Environmental niche modeling (ENM) has been extensively used to predict the potential distribution ranges of species and to evaluate their ecological niches under different evolutionary and climatic scenarios (Purcell and Stigall, 2021; Žerdoner Čalasan et al., 2021; Zhang et al., 2024; Vásquez-Cruz et al., 2025). These models have provided valuable insights into processes such as ecological speciation, hybridization niches, species isolation, and the impacts of climate change on biodiversity (Arenas-Castro and Sillero, 2021; López-Alvarez et al., 2015). The Pleistocene climatic cycles, characterized by alternating glacial and interglacial periods, profoundly influenced the demographic history and distribution patterns of Palearctic species (Friesen et al., 2016; Hewitt, 2004; Zhang et al., 2024).

Despite the considerable genomic and evolutionary knowledge gained from Brachypodium sylvaticum, no phylogeographic and environmental niche modeling analyses have been performed for this species complex, nor have grass-endophyte holobionts been considered as an evolutionary unit. In this study, we sampled B. sylvaticum accessions across its native Eurasian and North African range, retrieved genetic data from its ubiquitous Epichloë sylvatica endophytes, and integrated phylogenomic, population genetic, coevolutionary, and environmental niche modeling approaches to address key evolutionary questions. We reconstructed the environmental niches of the B. sylvaticum complex from two geoclimatic scenarios, the Last Glacial Maximum (LGM, ~20 thousands years ago (kya)) and the present. These models aimed to identify potential glacial refugia, evaluate oceanic vs. continental occurrence patterns, and infer postglacial range shifts that may have shaped the current distribution of the B. sylvaticum complex species. In this study, we aimed to investigate the evolutionary diversification and ecological dynamics of the B. sylvaticum complex and its ubiquitous fungal symbiont Epichloë sylvatica in the Palearctic. We hypothesize that Quaternary climatic oscillations possibly shaped the divergence of the eastern and western B. sylvaticum lineages, with niche conservatism restricting both lineages to mesic and oceanic habitats, and with possible asymmetric patterns of postglacial expansion. We also hypothesize the plausible coevolution of the B. sylvaticum host lineages with those of its endophyte E. sylvatica, due to the predominant vertical transmission of the fungus through seeds. Our objectives were to: (i) characterize the genomic diversity and phylogenetic relationship among B. sylvaticum complex populations using nuclear genome and plastome data, (ii) assess the coevolutionary dynamics between B. sylvaticum and its endophyte E. sylvatica, (iii) reconstruct historical and contemporary niche distribution patterns to identify glacial refugia and postglacial expansion routes, and (iv) synthesize these findings to elucidate the interplay between niche dynamics and grass–fungal interactions in shaping the plant species’ phylogeography.

2. Materials and methods 2.1. Sampling and genome skimming sequencing

A total of 94 samples belonging to the Brachypodium sylvaticum complex were selected for the genomic study. Sampling included 85 individuals of B. sylvaticum s.s., four B. miserum, two B. breviglume, and one B. kurilense, B. spryginii, and B. glaucovirens. The individuals were strategically collected across the respective native ranges, aiming to sample a comprehensive representation of the species’ genetic, geographical, and ecological diversities (Fig. 1; Appendix A). Six samples of other diploid Brachypodium species were selected as appropriate outgroups (two B. stacei, two B. distachyon, and two B. arbuscula) and three samples of the closely related B. pinnatum were selected as close ingroups, to adequately reconstruct the evolutionary framework of the B. sylvaticum complex samples under study (Catalán et al., 2023). Total genomic DNA was extracted from silica gel-dried leaf tissue from living plants or herbarium specimens using the CTAB method (Doyle and Doyle, 1987), with modifications to improve yield and purity (Decena et al., 2024). High-quality genomic DNA was subjected to genome skimming, a rapid and cost-effective sequencing method ideal for obtaining high-coverage genomic data. This method was performed from PCR-free libraries using Illumina technology at the Centro Nacional de Análisis Genómicos (CNAG, Barcelona) and Macrogen (Madrid). The quality of the sequencing reads was initially assessed using FastQC (Andrews, 2010), which provided detailed summaries of read quality metrics. Low-quality reads and adapter sequences were trimmed using Trimmomatic (Bolger et al., 2014). Reads that failed to pair were discarded to ensure data integrity, resulting in a total of 6.365–48.966 million reads.

2.2. Plastome assembly and annotation

Whole plastome sequences of the studied Brachypodium sylvaticum complex samples were assembled from the filtered genome skimming reads using Novoplasty 4.3.3 (Dierckxsens et al., 2017) and the B. sylvaticum Ain-1 reference plastome (https://phytozome-next.jgi.doe.gov/) as a seed to initiate the extension and assembly of the plastome samples under study. De novo annotation of plastome genes was performed using GeSeq (Tillich et al., 2017) implemented in Chlorobox (https://chlorobox.mpimp-golm.mpg.de/).

2.3. Brachypodium sylvaticum complex and Epichloë sylvatica read-filtering from Illumina genome skimming raw reads

To analyze genomic variation and gene sequences in the Brachypodium sylvaticum complex samples and their associated Epichloë sylvatica fungal endophytes, we employed a customized bioinformatic pipeline (SplitReads.sh). This workflow facilitated the read mapping to both grass and endophyte reference genomes, the extraction of gene sequences, and the filtering of genomic variants (SNPs). Initially, raw reads were mapped to a concatenation of two reference genomes, those of B. sylvaticum Ain1 (https://phytozome-next.jgi.doe.gov/info/BsylvaticumAin_1_v2_1, Lei et al., 2024) and Epichloë sylvatica ASM100826v1 (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_001008265.1/), using Minimap2 (Li, 2018). The output alignment files were sorted, indexed, and filtered into host and fungal alignments using Samtools (Danecek et al., 2021).

2.4. Gene sequence extraction and SNP detection in host Brachypodium sylvaticum complex samples

Mapped reads for host Brachypodium sylvaticum complex taxa were aligned to a concatenated CDSs reference from B. sylvaticum Ain1 to extract gene sequences for our Brachypodium samples. Consensus sequences were generated using Samtools, ensuring a minimum depth threshold of 10 reads per position to account for sequencing errors, obtaining relatively high-confidence host gene coding sequences for subsequent analyses. To identify genetic variants, a SNP-calling workflow was implemented using BCFtools (Danecek et al., 2021). Reads were processed with bcftools mpileup, using the reference genome to generate per-base sequencing depth and variant likelihoods. SNPs and small indels were called using bcftools and individual VCF files were merged using bcftools merge to create a consensus dataset of variants. We then used VCFtools (Danecek et al., 2011) to retain only biallelic SNPs, imposing a MAF (minimum allele frequency) of 0.05 and a maximum missing data threshold of 10%. The combined VCF file was subsequently filtered for linkage disequilibrium (LD) using plink2 (Chang et al., 2015), with an LD threshold of 0.2, and a sliding window of 10,000 bp.

2.5. ITS sequence extraction for Epichloë sylvatica endophytes from the holobiont samples

The pre-filtered Epichloë sylvatica reads were further mapped against a concatenated CDSs reference from E. festucae NCBI ID:35,717 to obtain gene annotations for comparative analysis. However, given the higher recovery of ribosomal genes over single-copy genes, an additional round of mapping was performed using only the complete rDNA 35S cistron (18S + ITS1 + 5.8S + ITS2 + 28S) to refine the detection and characterization of the E. sylvatica 35S sequences. To extract high-confidence fungal sequences, consensus sequences were generated using Samtools. Furthermore, the obtained data for the extensively studied ITS region (ITS1 + 5S + ITS2) were subjected to BLAST analysis to verify that the sequences belonged to Epichloë rather than representing potential contaminations, other endophytes, or chimeric sequences.

2.6. Population genomics and phylogeographic analyses of host Brachypodium sylvaticum complex plants

Phylogenomic trees of host plants were constructed using maximum likelihood (ML) methods to infer the evolutionary relationships among the studied Brachypodium sylvaticum complex taxa and populations. From an initial variant matrix comprising all high quality mapped sites (8,372,248), we pruned for maximum missing data per site ≤ 10%, minor allele frequency ≥ 0.05, and linkage disequilibrium pruning (r² < 0.2), yielding a final panel of 160,124 SNPs used to build an ML tree using IQ-TREE2 (Nguyen et al., 2015) with the best substitution model selected by the program, and 1000 ultrafast bootstrap replicates to assess branch support. Similarly, a ML phylogenomic tree was constructed for the whole plastome sequence dataset using IQ-TREE2. To investigate potential cytonuclear topological discordances between the nuclear and plastome genomes, we tested for potential cophylogenetic trees using the R packages ape and phytools. To quantify the extent of discordance, we employed the generalized Robinson–Foulds (gRF) distance, as implemented in the TreeDist R package (Smith, 2020), which quantifies topological discordance by measuring the proportion of bipartitions (splits) that differ between two phylogenetic trees, providing a normalized value between 0 (identical topologies) and 1 (completely incongruent).

A divergence dating analysis of the nuclear Brachypodium sylvaticum complex tree was performed using TreePL (Smith and O'Meara, 2012), a penalized likelihood method for estimating divergence times. The input topology was the maximum-likelihood phylogeny constructed from the concatenated nuclear SNP dataset (160,124 SNPs). Secondary calibrations were imposed on key nodes, following Catalán et al. (2023) and Sancho et al. (2018) estimations for the Brachypodium crown ancestor split (11.6 millions years ago (Ma)) and the core-perennial ancestor split (2.71 Ma). The optimal smoothing parameter was determined using cross-validation. Additionally, we generated 1000 bootstrap trees from the SNP dataset and dated each with TreePL to assess temporal uncertainty.

To infer the phylogeographic patterns of the Brachypodium sylvaticum complex lineages, we applied the BioGeoBEARS R package (Matzke, 2013), which implements likelihood-based models to reconstruct ancestral areas and infer biogeographic events such as dispersal, extinction, and vicariance. We performed the analysis on the dated phylogeny allowing null ranges and setting the maximum range size to nine. Geographic areas were manually defined based on the current distribution of the sampled populations, and their main evolutionary lineages (see Results), covering an East-to-West gradient: East Palearctic group (Kuril Islands, Japan, Himalayas), and West Palearctic group (East Europe, Middle East, Eastern Mediterranean, British Islands, Pyrenees + Alps, and Morocco + Southern Spain). To evaluate the fit of alternative biogeographic models, we performed model selection among the standard DEC, DIVALIKE, and BAYAREALIKE models, as well as their respective + J variants that incorporate founder-event speciation (jump dispersal). Model comparisons were based on log-likelihood and Akaike Information Criterion (AIC) scores, which were extracted and compared across models.

To investigate population structure and genomic ancestry, we employed ADMIXTURE v.1.3.0 (Alexander et al., 2009) using the 160,124 nuclear SNPs. This dataset was reduced to include only Brachypodium sylvaticum complex samples. This approach estimates individual ancestry by iteratively updating allele frequencies and ancestry fractions. We tested hypothetical genomic groups ranging from K = 1 to K = 10 and selected the best K using Structure Selector (Li and Liu, 2018), The selection of the best K was based on Puechmaille's method (Puechmaille, 2016).

2.7. Phylogenetic analysis of Epichloë sylvatica endophytes and potential plant-endophyte co-evolution

To investigate the potential co-evolution between the Brachypodium sylvaticum host plant and their Epichloë sylvatica fungal endophytes, we also constructed an Epichloë ML phylogenetic tree using the E. sylvatica ITS region (ITS1 + 5.8S + ITS2) alignment and IQ-TREE2. We performed a statistical ParaFit analysis of coevolution (Legendre et al., 2002), which test globally and individually for specific host-endophyte interactions (links). Parafit assumes the null hypothesis of independent evolutionary pathways versus the alternative hypothesis of coevolution; the analysis is based on patristic distances of both host and endophyte trees, and a host-endophyte interaction matrix over 1000 replicates is conducted using the ape package and the Parafit script (ParaFit.R) in R. We also employed the Procrustes Application to Cophylogenetic (PACo; Balbuena et al., 2013) analysis to evaluate host-endophyte phylogenetic congruence using the PACo package and a custom PACo.R script. For each 1000 ultrafast-bootstrap trees of hosts and endophytes we computed patristic distance matrices, then applied PACo with 1000 permutation replicates per tree to extract global goodness-of-fit statistics (m² and p-values) and squared residuals for each host–endophyte link. Residuals were normalized within each replicate, and for each link, we examined the distribution of its 1000 normalized residuals by binning into quartiles. Discordance for a given link was quantified as the combined proportion of replicates falling in the third and fourth quartiles (i.e., > 50th and > 75th percentiles).

2.8. Genomic diversity of Brachypodium sylvaticum complex taxa and genomic groups and impact of spatial isolation and environment isolation

Genomic diversity metrics, including observed heterozygosity (Ho), expected heterozygosity (He), and the inbreeding coefficient (F_IS), were estimated for the Brachypodium sylvaticum complex taxa and genomic groups retrieved from Admixture using the 160,124 SNP dataset with plink2 (Chang et al., 2015). The selfing rate (s) was further calculated as s = 2F_IS/(1 + F_IS) following Ritland (1990). To explore patterns of isolation by distance (IBD) and isolation by environment (IBE) in the Brachypodium sylvaticum complex populations, we first classified individuals into their corresponding populations (Table S1) using the dartR v.2.9.7 library (Mijangos et al., 2022). We then computed pairwise geographic distances between individuals using the Haversine method with the geosphere v.1.5–20 package, and pairwise bioclimatic distances based on a set of 19 bioclimatic variables retrieved from worldclim (current climate) using stats package in R. Pairwise genetic distances were estimated using Manhattan distances between individuals using the vegan R package.

To assess the relative contributions of geographic and environmental factors to the genomic differentiation of Brachypodium sylvaticum complex populations, we applied distance-based redundancy analyses (dbRDA) using the capscale function from the vegan R package. Genetic, geographic, and bioclimatic distance matrices were computed and transformed via Principal Coordinates Analysis (PCoA); axes of each variable were used as predictors. We fitted three models: (i) a full model including both environmental and geographic predictors, (ii) an environment-only model (IBE), and (iii) a geography-only model (IBD). The significance of each model and of individual predictors was assessed using permutation tests (n = 10,000). To disentangle the individual effects of IBE and IBD, we compared nested models using partial dbRDA, evaluating the contribution of environmental distances while controlling geography, and vice versa. This approach allowed us to test whether bioclimatic differences significantly explained genetic structure after accounting for spatial autocorrelation (indicative of IBE), and whether geographic distances explained structure independently of environmental effects (indicative of IBD). Variance partitioning and ANOVA comparisons further quantified the unique and shared contributions of both types of distances to population genomic differentiation.

2.9. Environmental niche modeling analysis of the Brachypodium sylvaticum complex

We performed environmental niche modeling (ENM) with occurrence data from the 94 Brachypodium sylvaticum complex samples used in the genomic study. To further increase the coverage of the species distributions, we also incorporated occurrence records downloaded from GBIF of taxonomically studied vouchers (GBIF.org −25 November 2020- GBIF Occurrence Download https://doi.org/10.15468/dl.jumys7): 168 for B. breviglume, three for B. glaucovirens, 28 for B. kurilense, 921 for B. sylvaticum s.s., and 48 for B. miserum (Fig. 1 and Table S1; Appendix A). To minimize sampling bias, occurrences with pair-distance < 50 km were excluded using the ‘thin’ function from the R package spThin (Aiello-Lammens et al., 2015). The final dataset included 1262 confirmed occurrence sites for B. sylvaticum complex samples classified into their respective taxa. Niche models were constructed for the whole B. sylvaticum complex dataset and separately for the Eastern Palearctic and Western Palearctic groups.

A set of 19 bioclimatic variables was retrieved from the WorldClim database (http://www.worldclim.org) at a spatial resolution of 2.5 arc minutes (~4.5 km² at the equator) for two different geoclimatic scenarios: present and Last Glacial Maximum (LGM, ~21 kya). The ENM analysis was conducted using the Maxent 3.4.1 algorithm (Phillips et al., 2017a) implemented in R through the ENMTools (Warren et al., 2021) and MaxEnt (Phillips et al., 2017b) packages. The most suitable climatic variables were selected by removing highly correlated variables based on the Variance Inflation Factor (VIF). Among correlated variables, we retained those with the greatest biological relevance for B. sylvaticum. The final predictor set comprised Annual Mean Temperature, Temperature Annual Range, Mean Temperature of Wettest Quarter, Annual Precipitation, Precipitation Seasonality, and Precipitation of Driest Quarter. To further mitigate sampling bias, we overlapped 100 independent MaxEnt models with climate data through bootstrap analysis. In each model iteration, 75% of the occurrence records were used for training, while 25% were used for testing. Model performance was assessed using the area under the curve (AUC) metric, which evaluates the discrimination ability of a model against a random one (Lawson et al., 2014). To eliminate areas with low environmental suitability and reduce overprediction, a presence threshold based on a 90% sensitivity criterion was applied, following model evaluation with presence–background data using the dismo package.

To evaluate the potential effect of the Quaternary glacial on the expansions or contractions of the Brachypodium sylvaticum complex niche, we constructed ENMs for the Last Glacial Maximum (LGM, ~22 kya), using paleoclimatic data from the ‘Community Climate System Model’ at 2.5 min resolution (CCSM; Gent et al., 2011) and assuming that the species' climate-driven biological requirements were the same in the past as today.

2.10. Comparing environmental niches and testing niche overlaps of the main Brachypodium sylvaticum groups

We compared the current environmental conditions of the Brachypodium sylvaticum complex, Eastern and Western groups, aiming to detect significant differences among them across 19 bioclimatic variables. In addition, we explored potential environmental variation within the B. sylvaticum complex using the same set of variables. Statistical analyses were performed using non-parametric Kruskal–Wallis tests followed by pairwise multiple comparisons to identify climatic differences among groups and Principal Component Analyses (PCA) to discriminate major environmental gradients, using the ecospat package in R (Broennimann et al., 2012).

We used a niche comparison framework to distinguish between niche overlap (the quantitative measure of the degree to which two groups occupy similar climatic space) and niche conservatism/divergence, which reflects the evolutionary processes underlying whether ecological tolerances have been retained or shifted. Niche overlap between the Eastern and Western groups was quantified using Schoener's D metric and the Hellinger-based I index, both ranging from 0 (no overlap) to 1 (complete overlap). Additionally, niche breadth was calculated using Levins' B2 standardized index. To assess whether observed overlap values indicated conservatism or divergence, we performed niche equivalency and niche similarity tests, based on 1000 random permutations. Following Warren et al. (2021), we interpreted significantly low overlap in the similarity test (p < 0.05 or p < 0.01) as evidence of ecological niche divergence, while non-significant results were associated with niche conservatism, i.e., the retention of ancestral climatic tolerances. We further compared niches from the same Last Glacial Maximum (LGM) group with present-day ones to evaluate temporal changes in both niche breadth and overlap.

3. Results 3.1. Nuclear dataset and phylogenomics of the Brachypodium sylvaticum complex

The total number of raw Illumina short reads across all samples ranged from 10,904,196 to 48,966,188 (Table S2). After filtering and mapping the reads to the Brachypodium sylvaticum Ain-1 nuclear reference genome, the number of mapped reads varied between 10,246,338 and 45,697,368. The total number of genes recovered with at least 90% of their sequence information ranged between 8020 for B. kurilense to 24,287 for B. sylvaticum Bsyl36H (Table S2). The initial whole-genome data matrix comprised 8,372,248 positions; after discarding gaps and applying a linkage disequilibrium (LD) filtering, we obtained a final dataset of 160,124 nuclear SNPs, uniformly distributed across the nuclear chromosomes of the B. sylvaticum complex samples, that were used in downstream analyses. Syntenic SNPs from other diploid Brachypodium species were retrieved from their respective genomic resources (Table S1). The strongly supported B. sylvaticum complex nuclear ML phylogenomic tree based on the nuclear SNP dataset revealed two sister lineages, the Eastern and Western Palearctic clades (Fig. 2). In the Eastern clade, the B. breviglume lineages were ancestral and paraphyletic, and the more recently evolved B. kurilense lineage was sister to the B. miserum clade. Within the Western clade, B. glaucovirens was recovered as the earliest diverging lineage, followed by two well-supported sister subclades. Subclade A included B. spryginii and B. sylvaticum s.s. lineages distributed across western, central, and eastern Europe, as well as Afghanistan, while subclade B comprised B. sylvaticum s.s. lineages from Greece, Denmark, British Islands, Spain and France (Pyrenees and Alps), Morocco, Turkey, and Iran. Most individuals from the same population or region clustered together, and the topology retrieved some geographical structure, although some individuals from the same population nested in separate clades (Fig. 2).

3.2. Biogeographic patterns and ancestral area reconstruction

The BioGeoBEARS ancestral range reconstructions inferred from the TreePL dated phylogeny, supported DIVA + J (Dispersal-Vicariance Analysis with jump dispersals) as optimal model. Divergence dating indicated that the Brachypodium sylvaticum complex lineage split from its stem ancestor in the early Pleistocene (Gelasian, ~1.99 Ma), and the crown node of the complex split into the sister Eastern and Western Palearctic clades in the Calabrian (~1.74 Ma). The estimated ages of the main Eastern clade lineages were slightly older than those of the main Western clade lineages (1.62 vs 1.47 Ma, respectively) (Fig. 2). The inferred geographic origin of the most recent common ancestor (MRCA) of the B. sylvaticum complex indicated a broad Palearctic area, including the Himalayas and the Mediterranean ranges (ADEGH), as the most likely place for it (Fig. 3).

Within the Eastern clade, inferred biogeographic patterns indicated that the crown ancestor and later ancestor of the Eastern group (Brachypodium breviglume) were distributed in the Himalayas (A; 1.62–1.34 Ma). This was followed by dispersal from the Himalayas to Japan (B) in the middle Calabrian (1.34–1.05 Ma) giving rise to the ancestor of the Pacific taxa. Then, in situ speciation gave rise to B. miserum in Japan (1.05–0.86 Ma), while a further dispersal from Japan to the Kuril Islands (C) and subsequent speciation gave rise to B. kurilense in the latter area (Figs. 2 and 3).

Within the Western clade, the most likely inferred ancestral area was a composite range including DEGH areas (Pyrenees-Alps, Middle East, eastern Mediterranean, and Morocco-southern Spain; Fig. 3). From this wide ancestral range, an estimated peripheral isolation event separated the B. glaucovirens lineage in the eastern Mediterranean (G) (~1.47 Ma; Figs. 2 and 3). Subsequently, a vicariance event at ~1.41 Ma likely split the subclade A, distributed in the Pyrenees-Alps region (D), from its sister subclade B, distributed in a large Middle East - Western Mediterranean range (EH). Subclade A comprised lineages predominantly associated with Pyrenean populations, which likely experienced three independent dispersal events: ⅰ) a first long-distance dispersal (LDD) to the Middle East (E) ~1.09 Ma; ⅱ) a second LDD to Morocco-southern Spain (H) ~0.70 Ma; and ⅲ) a third LDD to the eastern Mediterranean (F), which led to the divergence of B. spryginii and related eastern Mediterranean lineages at ~0.56 Ma (Figs. 2 and 3). Subclade B showed a more robust geographic structure, with an inferred early vicariance event separating a Middle Eastern lineage (E) from its sister Morocco-Southern Spain lineage (H) ~1.20 Ma. In the first group, three independent dispersals were estimated to have occurred between nearby geographic areas (Middle East to Eastern Mediterranean ~0.85 Ma, Eastern Mediterranean to Middle East ~0.61 Ma, and Eastern Mediterranean to Western Russia and Central Europe ~0.55 Ma). Within the second group, two independent colonization events of the Pyrenees-Alps from Morocco-southern Spain were estimated to have occurred around 0.95 Ma and 0.91 Ma, respectively, and more recent independent dispersals from the Pyrenees-Alps to the British Isles ~0.60 Ma and back to Morocco-southern Spain (~0.64 Ma, ~0.58 Ma and ~0.44 Ma) (Figs. 2 and 3).

3.3. Plastome dataset and nuclear discordance

The Brachypodium sylvaticum complex whole-plastome alignment included 137,356 positions. The total length of the assembled plastomes ranged from 136,073 to 136,571 bp (Table S2), with an average organelle coverage of 470×. The plastome sequences of the B. sylvaticum complex were highly conserved in terms of synteny and gene number, containing 133 genes (76 protein-coding genes, 20 non-redundant tRNAs, four rRNAs in both inverted repeats, four pseudogenes, and two hypothetical open reading frames) (Fig. S1). The plastomes of B. sylvaticum complex studied exhibited a 1143 bp insertion between the psaI and rbcL genes in the LSC region, corresponding to a rpl23 pseudogene. Additionally, they also showed a 301 bp insertion that aligns with a duplicated rps19 gene located between the psbA and trnH genes in the IRb repeat (Fig. S1).

The plastome phylogenetic tree of the Brachypodium sylvaticum complex (Fig. S2) showed a relatively similar topology to that of the nuclear tree (Fig. 2) for the major lineages, recovering the ancestral split of the Eastern and Western Palearctic lineages. However, the B. pinnatum lineage was resolved here as sister to the Western Palearctic lineage of B. sylvaticum, indicating that the B. sylvaticum complex is not monophyletic in the plastome-based reconstruction (Fig. S2). The Eastern Palearctic clade displayed the same resolution for its lineages as the nuclear tree, whereas successive divergences within lineages of the Western Palearctic clade differed from those of the nuclear tree (Figs. S2 and 4). The most ancestral splits of the plastome tree displayed complete support; however, the most recent splits showed poor support (BS < 50%) due to the similarity in plastome sequences between these closely related lineages (Fig. S2). A cytonuclear discordance analysis performed between the nuclear and plastomic phylogenies using the generalized Robinson-Foulds distance (Smith, 2020) indicated substantial incongruence between both trees (gRF = 0.7286; Fig. 4).

3.4. Detection of fungal endophytes and phylogeny of the Epichloë sylvatica symbionts

Evidence of symbiosis with Epichloë sylvatica was found in 83 of the 94 samples of the Brachypodium sylvaticum complex studied (Table S2). The total number of E. sylvatica reads ranged from 15,224 to 564,201, with the percentage coverage of the 35S rDNA cistron ranging from 63.82% to 89.44% in these samples. All Blast searches indicated that the assembled cistron sequences aligned with the E. sylvatica reference genome (Table S2). Among the 11 excluded individuals, six belonged to the Eastern Palearctic clade, specifically B. breviglume (Bbrev33H, Bbrev34H) and B. miserum (Bmis66-3, Bmis66-4, Bmis67-1, Bmis67-2), and eight samples belonged to the Western Palearctic clade (B. glaucovirens, Bsyl554-3B, Bsyl555-8, Bsyl59-2, Bsyl72-1). However, the failure to detect E. sylvatica in these samples could be due to the quality of their DNA or a lack of reliable detection, so they were discarded from the analysis.

The ITS phylogenetic tree of Brachypodium sylvatica revealed a well-to-moderately supported topology composed of two major clades (Fig. S3). The first clade (Clade Ⅰ) included endophytes found on host plants of B. kurilense, the only Eastern Palearctic species where the symbiont was detected, and lineages from Subclade B of B. sylvaticum from the Western Palearctic (plus three individuals from Subclade A). The second clade (Clade Ⅱ) consisted of endophytes found on B. sylvaticum from a wide geographic distribution, corresponding to lineages from Subclades A and B in the Western Palearctic (Fig. S3). E. sylvatica lineages with very low support values (BS < 20%) were collapsed into polytomies; these corresponded to samples with very similar or identical ITS sequences (Fig. S3).

Potential coevolutionary relationships between Brachypodium sylvaticum host plants and their E. sylvatica endophytes were assessed using the ParaFit and PACo methods. The global ParaFit test, based on the B. sylvaticum nuclear consensus phylogenies and E. sylvatica ITS tree over 1000 replicates, yielded a significant result, indicating overall phylogenetic congruence between hosts and symbionts (Table S3). However, only 11 individual host–endophyte interactions out of 83 were significant (Fig. 5 and Table S3), suggesting that congruence at the global level does not extend evenly to all pairwise associations. Interestingly, the majority of significant interactions corresponded to the Western Subclade A of B. sylvaticum, which included B. spryginii and B. sylvaticum s.s. individuals distributed across western, central, and eastern Europe, as well as Afghanistan (Figs. 2 and 5). In contrast, PACo analysis using 1000 bootstrap phylogenies also indicated statistically significant congruence, although the relatively high residual suggests a moderate level of phylogenetic fit (Table S4). Notably, PACo identified a higher number of consistent individual congruent assemblages, including 50 individual holobionts of B. sylvaticum and E. sylvatica from the Western Palearctic, 17 of 37 in subclade A and 33 of 46 in subclade B (Figs. S4 and S5).

3.5. Population genomics and genomic diversity of Brachypodium sylvaticum complex taxa and groups, and impact of IBD and IBE

The optimal number of population groups (K), determined using the Admixture and Structure Selector method, was K = 4 (Figs. 2 and S6). The first group corresponded to the Eastern Palearctic group -B. breviglume, B. miserum, B. kurilense (Figs. 2 and S6A, blue), with one sample (Bbrev34H, Pakistan) showing some admixture with group 2. The second group included B. spryginii and a subset of B. sylvaticum individuals belonging to the Western Palearctic A subclade (Figs. 2 and S6A, orange); these individuals were predominantly from Pyrenean populations, with one individual showing admixture with group 4 (Bsyl477-11, Spain). The third cluster included B. glaucovirens and B. sylvaticum samples from Morocco, southern Spain, the Middle East, and the eastern Mediterranean, corresponding to the Western Palearctic B subclade (Figs. 2 and S6A, magenta); it contained two samples showing admixture mainly with clusters 1 (Bsyl37H, Iran) and 2 (Bsyl39H, Turkey). The fourth cluster consisted of B. sylvaticum individuals of the Western Palearctic B subclade from the Pyrenees, the Alps, southern Spain, and the British Isles (Figs. 2 and S6A, green), with one individual (Bsyl550-10, Spain) showing some admixture with cluster 2.

All taxa within the Brachypodium sylvaticum complex studied exhibited low overall levels of genomic diversity (Tables 1 and S5), consistent with the predominantly selfing nature of the group (Khan and Stace, 1999; Catalán et al., 2016). Observed heterozygosity values ranged from 0.067 in B. miserum to 0.114 in B. glaucovirens, with correspondingly high inbreeding coefficients and selfing rates (F_IS = 0.113–0.553; s = 0.203–0.712) (Tables 1 and S5). All taxa except B. glaucovirens showed low heterozygosity and high F_IS values, consistent with this autogamous reproductive system. The comparatively higher heterozygosity and lower inbreeding in B. glaucovirens could reflect its putative hybrid origin (Catalán et al., 2023), although the small sample size limits a solid inference. In the ADMIXTURE group (Table 1), similar trends were observed. The Eastern Palearctic group showed the lowest observed heterozygosity (Ho = 0.071) and one of the highest F_IS values (0.534), translating to a selfing rate of 0.696. The groups combining B. sylvaticum with either B. glaucovirens (Cluster 3) or B. spryginii (Cluster 2) exhibited intermediate diversity and inbreeding values (F_IS = 0.428–0.438; s = 0.599–0.609), whereas the group composed exclusively of B. sylvaticum samples (Cluster 4) showed the highest heterozygosity (Ho = 0.090) and the lowest selfing estimate (s = 0.492), although still within the range expected for a selfing species. Overall, several F_IS values deviated significantly from Hardy–Weinberg expectations (Table 1), further supporting the strong effect of selfing in shaping genetic diversity patterns in this complex.

To disentangle the relative contributions of geographic distance (IBD) and environmental factors (IBE) to the genomic structure of Brachypodium sylvaticum complex populations, separate analyses were conducted for the Eastern Palearctic and Western Palearctic groups (Table 2). In the Eastern group (Table 2), neither geographic nor environmental distances explained a significant portion of the genetic variation, whether assessed independently or conditionally. These results suggest that the genomic structure within this clade may reflect other factors, such as demographic history, genetic drift, or recent divergence. In contrast, in the Western group (Table 2), although overall explained variances were low, the effects of IBD and IBE remained significant (p < 0.05), each accounting for approximately 13% of the genetic variation independently. After controlling for geographic distance, environmental factors still explained 10.79% of the variation, indicating that niche-related differentiation may also contribute to genomic structure in this region.

3.6. Bioclimatic analysis and environmental niche shifts of the Brachypodium sylvaticum complex taxa and its Eastern and Western Palearctic groups from the LGM to the present

Taxa and populations of the Brachypodium sylvaticum complex inhabit regions characterized by relatively mild thermal conditions and marked seasonal variability (Table S6). Across the entire distribution of the B. sylvaticum group (Fig. 1), the mean annual temperature averaged 9.43 ℃ (±4.58), with a marked temperature seasonality of 686.46 ℃ (±217.06). Thermal extremes ranged from a maximum temperature of 24.03 ℃ (±5.74) during the warmest month to a minimum of −3.39 ℃ (±6.48) in the coldest month. The annual temperature range was 27.42 ℃ (±7.39), with a mean diurnal variation of 8.88 ℃ (±2.50). In terms of humidity, the B. sylvaticum complex is associated with an annual precipitation of 851.56 mm (±415.91), with substantial contrasts between the wettest (116.98 ± 66.41 mm) and driest (34.74 ± 24.60 mm) months.

When the environmental data were divided into the Eastern Palearctic (including Brachypodium breviglume, B. miserum, B. kurilense) and Western Palearctic (including B. sylvaticum, B. glaucovirens, B. spryginii) groups, significant differences were detected for 16 of the 19 bioclimatic variables (Kruskall-Wallis, p < 0.05; Table S6). While both groups showed comparable mean annual temperatures (Bio1), isothermality (Bio3), and maximum temperatures in the warmest month (Bio5), the Eastern group was characterized by significantly more extreme climatic conditions. This included a greater diurnal temperature range (Bio2), greater temperature seasonality (Bio4), lower minimum temperatures during the coldest month (Bio6), and wider annual temperature ranges (Bio7). Furthermore, the Eastern group experienced significantly higher annual precipitation (Bio12) and greater seasonality of precipitation (Bio15), as well as marked differences in precipitation between the wettest and driest months (Bio13 and Bio14) and quarters (Bio16–19) (Table S6).

Principal Component Analysis (PCA) revealed that the first three axes explained 30.7%, 28.44%, and 17.66% of the total variance, respectively. In the 3D PCA plots (Fig. S7A and S7B), Brachypodium sylvaticum populations largely overlapped with those of the other microtaxa, indicating that overall niche differentiation within the B. sylvaticum complex is relatively subtle, although distinct clusters were observed for B. miserum, B. glaucovirens, and B. kurilense. Ecological niche modeling for the B. sylvaticum complex yielded highly accurate predictions, indicating excellent model performance (AUC = 0.941, SD = 0.001; Fig. 6A). The predicted current species distribution represented a fragmented, ocean-driven climatic niche, with greatest suitability along the northeast Atlantic margins as well as along northwest Pacific coasts, including Japan, the Kuril Islands, and an interior fringe of the Himalayas.

Paleoclimatic projections for the Last Glacial Maximum (LGM) indicate the persistence of glacial refugia in coastal, island, and lowland regions of western Europe (Atlantic margin), eastern Asia (Pacific margin), and in mountain ranges of western and southern Asia (Fig. 6B). These results support a longitudinal shift in distribution, with species likely shifting from glacial refugia in eastern and western coastal areas to more continental interiors as postglacial climates became progressively milder and wetter (Fig. 6). Separate models for the Eastern and Western Palearctic groups of the Brachypodium sylvaticum complex further highlighted their distinct ecological preferences and possible historical distributions (Fig. 6C–F). The Western Palearctic group displayed a primarily European distribution, with high suitability extending eastward to the Caspian region and scattered areas of moderate suitability in central Siberia, the Himalayas, and the Kuril Islands (Fig. 6C). In contrast, the Eastern Palearctic group exhibited a niche centered on China and the Himalayan region, with high probabilities of occurrence extending to Japan and the Pacific coast of East Asia, including the Kuril Islands, and lower suitability inferred for parts of Eastern Europe (Fig. 6E). LGM reconstructions suggest that both groups persisted in distinct glacial refugia: the western group along the Atlantic coast of Europe (Fig. 6D) and the eastern group in the Himalayas and East Asia (Fig. 6F). The present-day distributions of both groups (Fig. 6C–E) reflected partial postglacial expansions from these refugia.

Niche overlap between the Western and Eastern Palearctic groups was moderate according to Schoener's D (0.473) and Hellinger's I (0.702), with a non-significant result in the equivalence test (p = 1), indicating that the two niches were statistically indistinguishable when group identity was randomized. In contrast, the similarity test yielded a significant result (p = 0.027), suggesting that the observed overlap was greater than expected by chance given the environmental context (Table S7A). According to Warren et al. (2021), this pattern supports niche conservatism, implying that despite their deep genomic divergence and geographic separation, the two clades retain broadly similar climatic preferences. Standardized niche breadth values (Levins' B2) remained similarly low across time periods within each clade, with the western group showing a slight increase from 0.109 under present-day conditions to 0.117 during the LGM, and the eastern group increasing from 0.107 to 0.111 (Table S7B). Temporal niche overlap between present-day and LGM projections, as measured by Schoener's D, was moderate in both cases, with values of 0.625 for the western group and 0.487 for the eastern group, indicating partial shifts in climatic suitability over time.

4. Discussion 4.1. Evolutionary diversification of the main Brachypodium sylvaticum complex lineages and microtaxa

Our phylogenomic analyses, integrating genome-wide SNP and plastome data from 94 accessions, resolved the Brachypodium sylvaticum complex into two deeply divergent sister lineages, the Eastern Palearctic clade (Himalayas to the Pacific) and the Western Palearctic clade (Europe, North Africa, and Siberia). These clades diverged approximately 1.74 Ma (95% HPD: 1.91–1.69 Ma; Fig. 2), coinciding with Middle Pleistocene climate shifts that fragmented widely distributed ancestral angiosperm populations into refugia (Hewitt, 2004; Díaz-Pérez et al., 2018). This split partially corroborates previous findings (Catalán et al., 2023), but with unprecedented resolution, revealing reciprocal monophyly and distinct evolutionary trajectories for the Eastern and Western Palearctic lineages (Figs. 2–4 and S2). Importantly, however, the plastome-based phylogeny reconstructed the B. sylvaticum complex as paraphyletic, with B. pinnatum resolved as sister to the Western Palearctic lineage (Figs. 4 and S2), highlighting the proximity of the core perennial lineages B. sylvaticum and B. pinnatum and their high intercrossability (Khan and Stace, 1999), resulting in their shared matrilineal plastomes (Catalán et al., 2016).

Our broader sampling has also refined the geographic and taxonomic boundaries of the studied taxa within each lineage. The East Palearctic clade, comprising Brachypodium breviglume (Himalayas), B. miserum (Japan), and B. kurilense (Kuril Islands), exhibits a robust geographic structure with no evidence of hybridization (Figs. 2, 4 and S6A). The divergence dates (~1.74–1.62 Ma; Fig. 2) predate the Last Glacial Maximum and could correspond to earlier Pleistocene climatic and geological events. The emergence of the Great Kuril Range during the Early Pleistocene (Bulgakov, 1996) and subsequent glacial cycles could have promoted geographic isolation and speciation. Biogeographic reconstruction (DIVA + J model; Fig. 3) supports an ancestral distribution in the Himalayas for B. breviglume, followed by eastward dispersal and allopatric speciation. The paraphyly of B. breviglume and the sister relationship of B. miserum and B. kurilense (Figs. 2 and 3) suggest a stepwise colonization from the Himalayas to the Pacific coast, with subsequent isolation likely driven by Pleistocene sea-level changes that fragmented the temperate biota across East Asia (Qiu et al., 2011).

In contrast, the Western Palearctic clade (Brachypodium sylvaticum, B. glaucovirens, B. spryginii) revealed phylogeographic patterns driven by Pleistocene climatic oscillations. The SNP nuclear tree separated two well-supported subclades, Subclade A, which includes Pyrenean lineages and Caucasian B. spryginii, and Subclade B, a more heterogeneous assemblage including Mediterranean, North African, and British Islands populations. These two subclades showed different patterns of diversification and geographic expansion, with Subclade B exhibiting higher levels of admixture and secondary contact (Figs. 2 and S6A), suggesting contrasting demographic histories within the Western clade. Cytonuclear discordance (RF = 0.7286; Fig. 4) highlighted repeated chloroplast capture events, particularly in populations of B. spryginii (Caucasus) and Mediterranean B. sylvaticum, suggesting historical connections between glacial refugia in Anatolia and the Caucasus (Fig. S2). These cytonuclear discordances could have resulted from incomplete lineage sorting, historical introgression, and repeated plastid capture, processes that are well documented in other species of the genus (Sancho et al., 2018; Campos et al., 2024). The continuity of temperate deciduous forests in these regions during Quaternary climate oscillations likely facilitated gene flow or historical dispersal (Novák et al., 2023). However, phylogeographic patterns in the western A and B subclades reflected the survival of populations in warmer circum-Mediterranean refugia in the Early Pleistocene (> 1 Ma) and recent colonization of cooler high-latitude areas in the Late Pleistocene (< 1 Ma) (Figs. 2 and 3), consistent with similar survival and dispersal patterns in other angiosperms (Hewitt, 2004; Nieto Feliner, 2014). Genomic admixture analysis (K = 4; Figs. 2 and S6A) further revealed genomic admixture in Pyrenean and Balkan populations, reflecting secondary contact zones in these southern European refugia, where divergent lineages were found during the Holocene warming (Hewitt, 2004). Similarly, the topological incongruence observed between nuclear and plastome trees (Fig. 4) reflected the high level of introgression between geographically distant populations of B. sylvaticum, facilitated by the recurrence of LDD in the Western Palearctic regions during the Pleistocene (Figs. 2 and 3).

In this context, Brachypodium glaucovirens, a putative hybrid between B. sylvaticum and B. pinnatum (Scholz, 2007; Catalán et al., 2023), stands out for its high heterozygosity (Ho = 0.114 vs. 0.067–0.093 in other taxa; Table 1), supporting its hybrid origin. This is further corroborated by experimental and natural hybridization studies reporting intermediate individuals between B. glaucovirens and B. sylvaticum with 2n = 17 chromosomes and reduced pollen fertility, indicating spontaneous hybrid formation (Khan and Stace, 1999). Despite displaying partial sterility, these hybrids were able to backcross, suggesting the potential for introgression into sympatric populations. Notably, B. glaucovirens displays lower self-compatibility than B. sylvaticum, implying a tendency towards greater allogamy, which could further promote genetic exchange where their ranges overlap (Khan and Stace, 1999). The predominantly autogamous mating system of B. sylvaticum explains its low nuclear diversity (Ho = 0.067–0.090; Table 1), in contrast to allogamous relatives, such as B. rupestre, which displays higher heterozygosity (Ho = 0.13) and low inbreeding coefficients, even in clonal populations, due to mixed sexual and asexual reproduction strategies (Durán et al., 2025). Despite this, selfing facilitated rapid postglacial colonization events of B. sylvaticum, such as LDDs from western Europe to Iceland and Siberia (Figs. 2, 3 and 6) and invasive success in North America (Rosenthal et al., 2008), mirroring the “colonizer syndrome” of the selfing annuals B. distachyon and B. hybridum (Catalán et al., 2016; Wilson et al., 2019), and closely resembling the rapid spread of the highly selfing annual B. stacei across the Mediterranean through repeated LDDs and bottleneck events (Campos et al., 2024).

4.2. Coevolution and co-speciation of Brachypodium sylvaticum complex and its fungal endophyte Epichloë sylvatica

Our analysis of fungal ITS reads derived from total holobiont DNA successfully recovered Epichloë sylvatica sequences from nearly all Brachypodium sylvaticum complex accessions (> 95% infection rate), yielding nearly complete, high-quality ITS sequences across geographically diverse samples. However, because ITS sequences displayed very low levels of variation among individuals, the resulting phylogeny exhibited limited resolution, which constrains the robustness of deeper phylogenetic inferences and may partly explain the observed incongruences with the host phylogenies. Notably, this study represents one of the first to recover endophyte genetic data directly from host plant genome skimming datasets, highlighting the utility of this approach for symbiont characterization in untargeted genomic workflows.

Strikingly, Epichloë sylvatica was detected in the Eastern Palearctic taxon Brachypodium kurilense. It suggests E. sylvatica may occur more widely across the B. sylvaticum complex (Miwa et al., 2017), though its absence in other Eastern taxa (e.g., B. breviglume, B. miserum) could reflect sampling gaps or ecological exclusion. Notably, E. sylvatica has not been detected in related Western species like B. pinnatum and B. phoenicoides, which host distinct endophytes (e.g., E. typhina; Leuchtmann and Schardl, 2022), underscoring its unique specificity to the B. sylvaticum lineage. The ecological preferences of the B. sylvaticum complex taxa to forested humid habitats across the Palearctic (Fig. 6 and Table S6) probably fostered the infection of their presumably humid-adapted crown grass ancestor by the endophyte's ancestor since the early Pleistocene (Figs. 2–5, S3 and S4), resulting in the co-speciation of both organisms after 1.74 Ma of shared symbiosis. The phylogeny revealed two well-supported E. sylvatica clades, Ⅰ and Ⅱ (Fig. S3), demonstrating genetic structuring within this ubiquitous endophyte. Interestingly, the co-evolutionary analyses detected significant global concordance from both Parafit and PACo approaches (Figs. 5 and S4; Tables S3 and S4) while PACo further revealed a relatively high number of host-endophyte individual interactions (Fig. S4 and Table S4), thus supporting a general common evolutionary framework of the B. sylvaticum population and their respective endophytes. This pattern suggests a long-term association with the co-evolution reinforced by vertical transmission and the host's selfing mating system.

While ITS data could not distinguish mating (sexual) and non-mating (asexual) strains of Epichloë sylvatica (Leuchtmann and Schardl, 2022), we observed geographic structuring within fungal clades. For instance, identical E. sylvatica haplotypes occurred in distant populations (e.g., Pyrenees vs. Caucasus), aligning with prior reports of genetic homogeneity in this endophyte. Phylogenetic congruence analyses revealed complex dynamics. Both Parafit (p = 0.006; Table S3) and PACo (p < 0.01; Table S4 and Fig. S4) detected a significant signal of coevolution at a global level. However, while Parafit identified only localized associations, mainly related to holobionts of the western subclade A (12% significant links; Table S3 and Fig. 5), PACo detected a larger, albeit diffuse, number of individual associations (60.24%; Table S4 and Fig. S4). Although Parafit and PACo are based on different algorithmic principles, both have demonstrated the ability to detect complementary aspects of host-symbiont congruence, reflecting the fidelity imposed by vertical transmission in the highly self-fertilizing host (s = 0.492–0.712; Table 1). The algorithmic basis of PACo likely explains the broader host-endophyte signal recovered by this method, without necessarily implying differences in their sensitivity to horizontal transmission. The almost ubiquitous presence of E. sylvatica in B. sylvaticum holobionts suggests mutualistic benefits, potentially enhancing the ecological success of the host.

4.3. Past climate niche conservatism, glacial refugia, and postglacial spread

Our integrated phylogenomic and environmental niche modeling (ENM) analyses revealed distinct Quaternary responses of Brachypodium sylvaticum Eastern and Western Palearctic clades, which diverged ~1.74 Ma (Fig. 2). While both clades share adaptations to humid, oceanic habitats, their evolutionary trajectories and niche dynamics reflect independent histories shaped by contrasting glacial-interglacial pressures. The Western Palearctic lineage likely persisted in western Europe and Mediterranean LGM refugia -Pyrenees + Alps, Morocco + South Spain, Middle East, SW Europe- (Figs. 3 and 6D), in agreement with general refugial patterns identified for temperate taxa across the Mediterranean and Eurasia (Taberlet et al., 1998; Svenning et al., 2008; Médail and Diadema, 2009). Post-LGM expansion of the Western clade into northern Europe and western Asia—particularly into Scandinavia, the British Islands, and the Caucasus—is supported by phylogenomic and genetic patterns of genomic admixture, cytonuclear discordance, and shallow node divergence times (Figs. 2, 4 and 6), and by the biogeographic models (Fig. 3), which may oversimplify the complexity of recent postglacial migration. This expansion likely reflects a latitudinal migration from southern refugia to the North under increasingly oceanic postglacial climates, a trend observed in other mesic-adapted taxa (Leipold et al., 2017). In contrast, the Eastern Palearctic clade (Himalayas to Asian Pacific coast) exhibited spatial niche persistence since the LGM (Fig. 6E and F), maintaining its distribution within climatically stable mountainous and coastal regions. This pattern reflects long-term range continuity, rather than climatic niche stability, and is consistent with its long-term isolation, low genetic diversity (Table 1), and ancestral origin in East Asia–Himalayas inferred from dating phylogenies (Fig. 2; Catalán et al., 2023) and ancestral area reconstruction (Fig. 3). Oceanic islands (e.g., Japan, Kurils) and the Himalayan corridor likely acted as climatically stable refugia, buffering the Eastern Palearctic lineage against Quaternary climatic upheavals. Similar long-term persistence has been reported in other parts of temperate Asia, including periglacial lowlands, as supported by ancient DNA analyses and plant regeneration studies from ancient permafrost (Yashina et al., 2012; Thomsen and Willerslev, 2015).

Niche overlap between clades was moderate (Schoener's D = 0.473; Hellinger's I = 0.702), with statistically significant similarity (p-value = 0.027) and equivalent niches (p-value = 1.00), indicating that the niches occupied by the Eastern and Western Palearctic lineages are more similar than expected by chance. This pattern is best interpreted as evidence of niche conservatism, suggesting that divergence between clades was likely due to geographic and historical isolation, rather than adaptive ecological speciation. Despite profound genomic divergence and long-term separation, both clades converged on similar mesic climatic niches, likely constrained by ancestral ecological preferences.

Temporal comparisons of niches between the LGM and the present revealed asymmetric niche shifts, the Western Palearctic clade maintained moderate stability (present-day LGM Schoener's D = 0.625), whereas the Eastern Palearctic clade experienced greater climatic niche turnover (present-day LGM Schoener's D = 0.487), likely due to increased climatic fragmentation and heterogeneity in high-altitude refugia in East Asia. This is consistent with phylogeographic evidence from the region showing that Quaternary climatic oscillations promoted population isolation and profound lineage divergence in the Himalayas, eastern China, and Japan due to mountain barriers and habitat fragmentation (Qiu et al., 2011). Both lineages showed a slight expansion of niche breadth (Levins B2: western clade 0.109–0.117; eastern clade 0.107–0.111), suggesting gradual postglacial diversification rather than abrupt ecological changes.

Together, these ENM results provide an independent climate perspective that supports our phylogeographic interpretations, as the spatial patterns of refugia, postglacial expansions, and long-term persistence inferred from the models are consistent with the genomic signals of divergence, admixture, and isolation observed in Brachypodium sylvaticum lineages.

Acknowledgements

We thank the handling editor and three anonymous reviewers for their valuable comments on an earlier version of the manuscript, L.A. Inda Aramendía and J. Ascaso Martorell for their assistance with field sampling and collection of Brachypodium sylvaticum specimens, M.F. Moreno-Aguilar for helpful support with the implementation of custom scripts, the herbaria of B (Berlin), M (München), LD (Lund), VLA (Vladivostok), RIKEN (Japan) for providing access to voucher material, and to these herbaria and to Missouri (MO), gbif.org and Hokkaido University, for providing specimens and images that were used to generate the drawings of the studied taxa shown in Fig. 1. This study was funded by the Spanish Ministry of Science and Innovation PID2022-140074NB-I00, TED2021-131073B–I00, and PDC2022-133712-I00, and the Spanish Aragon Government and European Social Fund Bioflora A01-23R research grants. MD and MC were supported by their respective Spanish Ministry of Science and Innovation Ph.D. fellowships and VS by a Tomsk State University Ph.D. fellowship. DC was supported by a young research contract grant.

Data availability

Supplementary data to this article, all the data matrices assembled and related phylogenomic and biogeographic results in this study can be found online at https://github.com/Bioflora/BsylvaticumPhylogeography. The raw sequencing data have been submitted to the European Nucleotide Archive (ENA) under the project numbers PRJEB64465 and PRJEB89851.

CRediT authorship contribution statement

Ma Ángeles Decena Rodríguez: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Data Curation, Writing – Original Draft, Writing – Review and Editing, Visualization. Miguel Campos Cáceres: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – Original Draft, Writing – Review and Editing, Visualization. Diana Isabel Calderón Pardo: Validation, Writing – Original Draft, Writing – Review and Editing, Visualization. Valeriia Shiposha: Writing – Review and Editing. Marina Olonova: Writing - Review and Editing. Ernesto Pérez Collazos: Conceptualization, Methodology, Investigation, Resources, Writing – Original Draft, Writing - Review and Editing, Visualization, Supervision. Pilar Catalán Rodríguez: Conceptualization, Investigation, Resources, Writing – Original Draft, Writing – Review and Editing, Visualization, Supervision, Project administration, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.12.015.