Semipermeable species boundaries create opportunities for gene flow and adaptive potential

Hybridisation and gene flow can have both deleterious and adaptive consequences for natural populations and species. To better understand the extent of hybridisation in nature and the balance between its beneficial and deleterious outcomes in a changing environment, information on naturally hybridising nonmodel organisms is needed. This requires the characterisation of the structure and extent of natural hybrid zones. Here, we study natural populations of five keystone mound‐building wood ant species in the Formica rufa group across Finland. No genomic studies across the species group exist, and the extent of hybridisation and genomic differentiation in sympatry is unknown. Combining genome‐wide and morphological data, we demonstrate more extensive hybridisation than was previously detected between all five species in Finland. Specifically, we reveal a mosaic hybrid zone between Formica aquilonia, F. rufa and F. polyctena, comprising further generation hybrid populations. Despite this, we find that F. rufa, F. aquilonia, F. lugubris and F. pratensis form distinct gene pools in Finland. We also find that hybrids occupy warmer microhabitats than the nonadmixed populations of cold‐adapted F. aquilonia, and suggest that warm winters and springs, in particular, may benefit hybrids over F. aquilonia, the most abundant F. rufa group species in Finland. In summary, our results indicate that extensive hybridisation may create adaptive potential that could promote wood ant persistence in a changing climate. Additionally, they highlight the potentially significant ecological and evolutionary consequences of extensive mosaic hybrid zones, within which independent hybrid populations face an array of ecological and intrinsic selection pressures.


| INTRODUC TI ON
Hybridisation and gene flow between closely related species are common in nature and predicted to increase because of human activities and climatic warming (Chunco, 2014;Mallet, 2005;Rieseberg, 2009;Scheffers et al., 2016), as changing environmental conditions induce range shifts and connections between previously allopatric populations. This may, over longer evolutionary time scales, either collapse species barriers (Owens & Samuk, 2019) or reinforce them (Lemmon & Juenger, 2017;Pfennig, 2016). For some time, the species barriers can also remain semipermeable, a situation where the taxa remain distinct in their natural habitat despite gene flow between them.
Hybridisation and gene flow have both positive and negative effects for fitness. Complete reproductive isolation is not necessarily optimal (Barton, 2020): increased genetic diversity through hybridisation may fuel adaptation Torda & Quigley, 2021), and introgression has repeatedly underlied adaptation into new niches (De-Kayne et al., 2022;Kagawa & Takimoto, 2018;Meier et al., 2017). However, populations can also suffer from hybridisation due to deleterious fitness effects (Ålund et al., 2013) and mortality (Ellison et al., 2008), arising from the inviability of offspring or hybrid breakdown in later generations. Revealing the extent and consequences of hybridisation in natural populations is needed to clarify the fate and adaptive potential of hybrid populations in the face of changing environment, like warming climate.
The geographical extent and type of hybrid zone greatly influence the evolutionary and ecological impacts of hybridisation: In contrast to more commonly studied, geographically restricted clinal zones (Barton & Hewitt, 1985), some taxa exhibit mosaic hybrid zones.
These may be as large as the species' sympatric area. In mosaic zones (Arnold, 1997), hybridisation events are typically independent, and recently admixed populations form a spatial mosaic with nonadmixed and backcrossed populations of at least two species (Arnold, 1997;Harrison, 1986;Rand & Harrison, 1989). Under the mosaic zone model, the species are assumed to be adapted to different environments, and their hybrids may either have lower fitness (due to intrinsic selection) or be favoured (extrinsic, environment-dependent selection) (Arnold, 1997). Mosaic hybrid zones have been reported, for example for various plant species (Abbott & Brennan, 2014;Rieseberg et al., 1999) and molluscs (Fraïsse et al., 2014). Since mosaic zones can be large, and selection can act independently in different hybrid populations, they may result in more heterogeneous patterns and geographically widespread effects than those in clinal hybrid zones. Furthermore, independently evolving hybrid populations allow predictability in the outcomes of hybridisation to be studied (Nouhaud, Martin, et al., 2022).
Hybridisation is commonly investigated between two species.
However, genomic data have facilitated the detection of multispecies hybridisation in many taxa, including white oaks (Reutimann et al., 2020), birds (Grant & Grant, 2020;Natola et al., 2022;Ottenburghs, 2019), fish (Banerjee et al., 2022), and butterflies (Heliconius Genome Consortium, 2012). Multispecies hybridisation is likely to be more frequent in sympatric groups with a low degree of reproductive isolation, and one species may act as a 'conduit' of gene flow between two otherwise isolated species (Grant & Grant, 2020).
Wood ants of the Formica rufa species group (Hymenoptera, Formicidae) represent ideal models of hybridisation and understanding its role in adaptation: they have diverged recently, within the past 500,000 years (Goropashnaya et al., 2004(Goropashnaya et al., , 2012Portinha et al., 2022), have a wide sympatric range in Eurasia including Finland and have different climatic adaptations (Table 1 and Martin-Roy et al., 2021). Furthermore, there are premises for the formation of mosaic hybrid zones: morphological and genetic marker findings indicate extensive hybridisation across the F. rufa group within their sympatric Eurasian ranges (Seifert, 2018(Seifert, , 2021Seifert et al., 2010).
Albeit reproductive isolation is not complete within the group, both pre-and postzygotic isolation, including intrinsic and environmentally dependent selection, have been detected Kulmuni & Pamilo, 2014;Seifert, 2018). Wood ant hybridisation has resulted in male-biased mortality during development, but also female heterosis and potential for thermal adaptation (Kulmuni & Pamilo, 2014;Martin-Roy et al., 2021). Instability in the barriers of gene flow has been reported previously (Kulmuni et al., 2020) in this system, suggesting that the balance between the costs and benefits of hybridisation may vary across space and time. Understanding hybridisation in the F. rufa group is important since these ants are keystone species in boreal and mountain forests across Eurasia (Stockan & Robinson, 2016;Trigos-Peral et al., 2021). They have a significant role in the forest ecosystem: they build stable nest mounds (Seifert, 2018), and populations of polygynous (i.e. multiple queens per nest) species in particular can reach densities of hundreds of nests per km 2 , and impact ecosystem characteristics, ranging from nutrient cycling below ground to aphid farms in the forest canopy. In TA B L E 1 Characteristics of the study species.  (Yarrow, 1955), which is the most common of these species in Finland (Punttila & Kilpeläinen, 2009).
Multiple hybrid populations of F. aquilonia and F. polyctena have been genetically characterised within a small region in southern Finland (Beresford et al., 2017). In addition, morphological and allozyme data support several other observations of F. aquilonia and F. polyctena hybridisation, mainly in southern Finland (Pamilo & Kulmuni, 2022;Seifert, 2021;Sorvari, 2022), and indicate F. aquilonia and F. lugubris hybridisation in northern Finland (Pamilo & Kulmuni, 2022). Since many species pairs within the group hybridise, this raises the possibility of previously undetected multispecies hybridisation. However, outside F. aquilonia and F. polyctena (Nouhaud, Martin, et al., 2022;Portinha et al., 2022), no genomic studies on the F. rufa group exist, and consequently also no information on the congruence between morphological and genomic data is available.
Under a changing climate, semipermeable species boundaries and hybridisation between differentially adapted species may provide an opportunity to reshuffle adaptations and help species to cope with environmental shifts. Hybridisation between the F. rufa group species may have potential for this, since the five species targetted here differ in their life history characteristics and habitat preferences (Seifert, 2018) and, based on their Eurasian distributions, are expected to be adapted to different climatic conditions ( and are assumed to be warm-adapted (Table 1) (Seifert, 2018).
Since successful wood ant hibernation, reproduction and brood development are temperature-dependent, climatic conditions that affect within-nest temperatures are critical for their fitness and survival throughout the year (Frouz, 2000;Frouz & Finer, 2007;Kadochová & Frouz, 2014). However, the ants can enhance their fitness by active nest thermoregulation from early spring until autumn (Horstmann & Schmid, 1986;Kadochová & Frouz, 2014). Active regulation is paused during hibernation, but recovers in the spring, when the species start to reproduce at different but overlapping times (Seifert, 2018), with warmth-dependent brood production and development continuing until summer. In the winter, high hibernation temperatures may increase ant metabolism, leading in contrast to reduced lipid reserves, and thus reduced fitness and increased mortality, as has been shown in the laboratory for cold-adapted F. aquilonia (Sorvari et al., 2011). Specifically, workers need these body fat reserves in early spring for nest heating (Seifert, 2018), and queens for successful sexual reproduction (Sorvari et al., 2011).

| Morphological data
We carried out state-of-the-art morphological species identification with NUMOBAT (Seifert, 2021) using three to six individuals from each nest sampled for genomic data. Since both the genomic analysis and morphological identification required the use of whole individuals, we used different individuals from each nest for these separate purposes, which is justified by the fact that wood ant nests are genetically homogeneous (Kulmuni et al., 2020;Nouhaud, Martin, et al., 2022;Pamilo, 1982). We measured altogether 17 morphological characteristics (See Table S2 for a description of the characteristics, and Table S3 for primary morphological data), producing species assignment probabilities at the nest level. The morphological characteristics measured and their recording is described in detail in Seifert (2021).
All 77 locations with NUMOBAT data (69 nests collected for this study and eight nests analysed in Portinha et al., 2022, see Table 2 and  (Seifert, 2021) and do not show any phenotypical indications of a hybrid identity.
There was only a small amount of reference data for F. lugubris and F. pratensis for six of the morphological characteristics. Hence, we ran first all samples in a 5-class LDA using 11 morphological characteristics for which sufficient reference data was available (see also

| Genomic data
For the genomic characterisation of our study species and individual samples, as well as documenting the extent of hybridisation within the F. rufa group in sympatry, we studied nuclear admixture and mito-nuclear mismatch. To achieve this, we performed whole-genome sequencing for our 69 new samples and analysed these along with the 20 previously sequenced genomes   (Table S4).

Nuclear genomic data
Single nucleotide variants (SNPs) for the nuclear genome were called with the FREEBAYES software (v1.3.1, -k option was used for disabling population priors; Garrison & Marth, 2012) and the VCF was normalised using VT (v2015.11.10; Tan et al., 2015). Sites located at less than two base pairs from indels were excluded, as were sites supported by only Forward or Reverse reads, using BCFTOOLS (v1.10; Danecek Only biallelic SNPs with quality equal or higher than 30 were kept. Sites were then filtered based on individual sequencing depth: individual genotypes that had more than twice the mean depth of the individual in question were set as missing. Genotyping errors that would occur, for example due to misaligned reads were removed by excluding sites that displayed heterozygote excess after pooling all samples (p < .01, VCFTOOLS v0.1.16; Danecek et al., 2011).
Subsequently, individual genotypes with quality <30 and depth of coverage <8 were coded as missing data. Sites with more than 10% missing data over all samples were discarded. Finally, singletons were removed (minor allele count <2) with VCFTOOLS (v0.1.16).
After these steps, 1,829,565 SNPs from 89 individuals with sequencing depths of at least 8× remained. All the genetic analyses (detailed below) were performed by subsampling 9770 genome-wide SNPs. These were obtained by thinning the final dataset with 20 kb distances to minimise linkage disequilibrium, as well as by excluding chromosome number three (social chromosome with reduced recombination due to supergene inversions; Brelsford et al., 2020) and shorter contigs with no known genomic location.

Mitochondrial genomic data
Single nucleotide variant calling for the mitochondrial genome was done with FREEBAYES and a frequency-based approach. The SNP filtering followed the same pipeline that was used for the nuclear genome except that the steps after the retention of biallelic SNPs with a quality greater than or equal to 30 were not applied to the mitochondrial data. We utilised two lines of evidence to examine genome-wide differentiation and to identify admixed individuals and the extent of introgression in our study species: we investigated individual nuclear genomic ancestry proportions with ADMIXTURE, and studied mitonuclear mismatches. These results were supported by information from the phylogenetic network and principal component analyses (PCA). We used the genomic analyses in combination with the morphological identification from NUMOBAT, which enabled the identification of the samples at the species level.

| Admixture analysis
We performed ADMIXTURE (Alexander et al., 2009) with the 9770 nuclear SNP data to assess the extent of hybridisation in the F. rufa group in sympatry, to reveal whether the samples group genomically by species, and to detect nuclear genomic signs of hybridisation. All five species were included, but we aimed at balanced sample sizes of nonadmixed individuals (from six to eight individuals per species) and thus randomly selected six F. aquilonia individuals (see Table S1), since unequal sample sizes could yield biased clustering results (Toyama et al., 2020). For ADMIXTURE, the number of clusters K was chosen to range from K = 2 to 5, with prior knowledge of sampling from five morphologically identified species. Random selection of six F. aquilonia individuals was repeated multiple times to confirm that certain individuals did not drive the results. The most likely value of K was chosen based on cross-validation error (CV, prediction of withheld data; Alexander & Lange, 2011).

| Mitotype network
To detect hybridisation via mitonuclear mismatches, we characterised the mitotypes of our samples using the mitochondrial SNP data.

| Observed heterozygosity
To characterise the hybrids, we also studied the observed heterozygosity of our Finnish samples (N = 77, see Table S1). Specifically, we tested whether admixed individuals between F. aquilonia and the F. polyctena/F. rufa clade, whose adaptive potential we were interested in studying due to the prevalence of their hybridisation, We thus recovered 305 Ancestry-Informative Markers (AIMs) that were used to compute HI and ancestry heterozygosity. For these analyses, we excluded two admixed samples with more than 50% missing data for the AIMs (see Table S1).  Table S1).

|
We generated climatic time series for each of these nests extending back 1, 5, and 20 years from the nest-specific field sampling date.
The time series data were extracted from digital rasters of daily mean air temperature and daily precipitation sum (sourced from the To test whether F. aquilonia and hybrid nest sites diverged climatically within their sympatric range, we fitted Bayesian logistic regression models using the R package brms (Bürkner, 2017) with nest status (hybrid vs. F. aquilonia) as the response variable. We fitted 10 separate models, each with either one of the climatic covariates, or site latitude or longitude, as the explanatory variable.
We did not include multiple explanatory variables in a single model due to high multicollinearity between the climatic variables and the coordinates, and the relatively low sample size. We normalised the explanatory variables to mean zero and unit variance, and used an informative normal(0,2) prior distribution for the model parameters.
We fitted the models using four chains each with 4000 iterations and discarded the first 2000 iterations as warm up resulting in a total of 8000 posterior samples. The convergence and mixing of the fitting algorithm was checked with the Rhat and Effective Sample Size values. The models were ranked according to their expected log pointwise predictive density from leave-one-out cross-validation using the loo R package (Vehtari et al., 2017). We also checked whether the exclusion of three nests that are spatial and climatic outliers, due to their location in the Åland islands in the Baltic Sea, altered the model results.

| Within-nest temperature
Lastly, we studied whether ambient temperature affects withinnest temperatures as experienced by the ants, and hence whether our regional climatic modelling results are likely to have biological relevance. For this purpose, we tested whether within-nest temperatures correlated with ambient temperatures and whether they differed between hybrid and F. aquilonia nests. Three of the five ant populations used in these analyses (for which internal nest temperature data were available) were same, and two different from those used in the above morphological, genomic, and ambient climatic analyses (see Table S1).
In the years 2020 and 2021, we deployed Hobo data loggers to record within-nest temperatures hourly in 42 ant nests located within three admixed and two F. aquilonia populations in southern Finland ( Figure S1). In total, four populations were sampled in 2020 and three populations in 2021; two of these were sampled in both years. In most cases, data logging ran continuously from 1st January until 31st July. However, data logging in one of the 2021 populations ('Svanvik') ran during May-July only. Daily mean within-nest temperatures per nest were calculated from the hourly values. Daily mean ambient temperatures for the corresponding dates were then extracted for each ant population based on the mean coordinates of nests within that population from the ClimGrid gridded daily climatology rasters at a 1-km spatial resolution (Aalto et al., 2016).
We modelled the mean daily within-nest temperatures of these nests in selected winter (January, February), spring (March-May) and summer (June, July) months as a function of nest status (hybrid vs. F. aquilonia), site ambient temperature and sampling year (2020 vs. 2021) using the R package brms (Bürkner, 2017). We compared a simple additive model structure (

| Wood ant species maintain their genomic distinctiveness despite hybridisation
Despite the incomplete reproductive isolation and extensive hybridisation in sympatric F. rufa group wood ants we found support for four nonadmixed F. rufa group species (F. aquilonia, F. lugubris, F. rufa, and F. pratensis) in our Finland-wide sample, indicating genomically distinct species rather than a panmictic gene pool. This result is based on ADMIXTURE, PCA, mitochondrial haplotypes, and the phylogenetic network, detailed below (Figures 1 and 2, for PCA see Figure S2). TA B L E 3 NUMOBAT LDA probabilities of assigning nonadmixed samples to the species in which they belong to according to genomic data.
Formica pratensis was genomically most distinct from all the others and formed a well-separated genepool in both nuclear and mitochondrial analyses (Figures 1a,c and 2, Figure S2A). None of the seven F. pratensis populations showed signs of admixture.

| Hybrid populations are not transient and the hybrids have higher heterozygosity than the nonadmixed individuals
We were interested in the amount of genetic variation as a proxy for adaptive potential especially in the case of individuals admixed between F. aquilonia and the F. polyctena/F. rufa clade ('Admixed 1'), since they are abundant in the landscape. These populations had significantly higher observed overall heterozygosity than either sympatric nonadmixed F. aquilonia (Wilcoxon Rank Sum test, z = 4.92, p < .001) or F. rufa (z = 3.38, p < .001) individuals ( Figure 3).

Furthermore, Ancestry Informative Markers (AIMs) indicate that
these were late-generation hybrids, since heterozygosity based on AIMs was well below 0.5, suggesting hybridisation leads to persistent hybrid populations. Furthermore, hybrid individuals from different populations vary in ancestry, supporting the idea that the hybrid populations evolve independently of one another ( Figure 4).

| Hybrids occur at sites with warmer temperatures and less spring precipitation than F. aquilonia populations
Our results indicate that the hybrid populations between F. aqui-  The absolute values of these effect estimates are highly correlated with the difference in the expected log pointwise predictive density (elpd-diff) calculated with loo, which was used to rank the models according to their predictive performance ( Figure S4).

| Ambient temperature impacts within-nest temperature and hybrid nests are warmer than F. aquilonia nests
Monitoring of three hybrid and two F. aquilonia populations shows that ambient temperature has a linear effect on within-nest temperatures in most months ( Table 4, Table S6). Furthermore, even when controlling for the ambient temperature and inter-year effects, hybrid nests were internally warmer, on average, than F. aquilonia nests.
This result is consistent in models with and without an interaction between nest status and ambient temperature ( Table 4, Table S6) from January to March and in June. In addition, during January and April-May, hybrid nests warmed up more efficiently with respect to increases in ambient temperatures (see interaction model, Table 4).
In June, the hybrid nests tended to warm up more slowly than the cooler F. aquilonia nests ( Table 4).

| Genomic distinctness in sympatry, and mosaic hybrid zone with broad implications between differentially adapted F. rufa group species
Despite our results on extensive hybridisation and the finding that all the studied species hybridise with one or more other species in sympatry, we detected both whole-genome and mitochondrial divergence among most of them in Finland. Specifically, four out of the five studied species-F. aquilonia, F. lugubris, F. rufa and F. pratensis-were found in nonadmixed populations in Finland. Our current dataset, in which all the F. polyctena populations are admixed in sympatry, does not allow testing whether F. polyctena and F. rufa represent different species or ecotypes of the same species. The mosaic hybrid zone that we detected (ca 60°N to 63°N) greatly expands the known area and extent of wood ant hybridisation in Finland (Beresford et al., 2017;Pamilo & Kulmuni, 2022;Seifert, 2021). We report 10 new populations that are admixed between F. aquilonia and the F. polyctena/F. rufa clade, bringing the total within this dataset, comprising both published and unpublished genomes, to 14. The hybrid zone is likely not restricted to F I G U R E 5 Bayesian logistic regression using climatic data from a total of 33 sympatric F. aquilonia populations and populations admixed between F. aquilonia and the F. polyctena/F. rufa clade ('Admixed 1') shows that they occur in different microclimates. Hybrid nest locations are warmer both in the winter and in the spring, and receive less precipitation in the spring than F. aquilonia nest locations. The climate data spans 20 years backward from the nest sampling year. Covariates that significantly correlate with hybrid occurrence have nonoverlapping credibility intervals (95% credibility intervals are given) with the midline (zero). The winter hibernation season ('Winter') = December-February; spring reproductive season ('Spring') = April-May. Upper quartile temp. = upper quartile of winter temperatures, N days >0°C = number of >0°C winter days, Mean temp. = annual mean winter or spring temperature, Mean prec. = annual mean winter or spring precipitation, Lower quartile temp. = lower quartile of spring temperatures, N days <0°C = number of spring frost days.
Winter Spring our study area, but may continue across the species' sympatric Eurasian area, as morphological evidence suggests (Seifert, 2021).  (Arnold, 1997;Barton & Hewitt, 1985), but also by local conditions within each patch within the broad hybrid zone as expected according to the mosaic hybrid zone model (Arnold, 1997).
Second, in contrast to clinal zones, mosaic zones typically comprise a mixture of environmentally dependent and independent selection pressures, the strength of which can vary among patches (Arnold, 1997). Clinal hybrid zones, in contrast, are usually assumed to be arrayed linearly along a habitat gradient (Arnold, 1997). The association that we observed between hybrids and warmer microclimatic conditions aligns with earlier evidence of thermally dependent selection (Martin-Roy et al., 2021) that may favour hybrids.
This could counteract hybrid breakdown due to intrinsic incompatibilities Kulmuni & Pamilo, 2014) and reduced hatching rates (Beresford, 2021) which have been detected previously in wood ant hybrids Kulmuni & Pamilo, 2014). Third, variation in F. aquilonia ancestry and ancestry heterozygosity, combined with the limited dispersal capabilities of supercolonial wood ants, suggest that hybrid populations may have emerged from independent admixture events. This finding is consistent with previous research on three closely located populations (Nouhaud, Martin, et al., 2022). The mosaic of ancestry proportions across the hybrid zone contrasts with the spatial gradients typical of clinal hybrid zones (Barton & Hewitt, 1985). Together with relatively low ancestry heterozygosity, indicative of later generation or backcrossed hybrids, this attests to long-lived hybrid populations that are not under overwhelmingly negative selection. In summary, the evolutionary outcomes of mosaic hybrid zones may differ from those of clinal zones. This is because hybrid populations may evolve independently and be shaped in various ways by both ecological and intrinsic selection pressures and interactions with the nonadmixed populations leading to a diversity of hybrid lineages and phenotypes (Harrison, 1986;Rand & Harrison, 1989).
Which selection pressures drive the genomic and species composition of the mosaic zone is an interesting question for future studies. They may be, for example, affected by the degree of forest fragmentation, since in Central Europe F. polyctena and F. rufa hybrids are more prevalent in small fragments than large forests . Furthermore, the interplay of ecological factors such as temperature, habitat preferences, differences in mating times (Martin-Roy et al., 2021;Seifert, 2018;Stockan & Robinson, 2016), and intrinsic selection on genetic incompatibilities Kulmuni & Pamilo, 2014) makes predicting the evolutionary dynamics of such a mosaic challenging. Our unique setting of more than 10 independent hybrid populations will provide an unrivalled opportunity to study the dynamics and outcomes of hybridisation in natural populations.
We predict that the northern border of the hybrid zone, that is the northern range limit of F. rufa (and F. polyctena), is likely to shift TA B L E 4 The relationship of within-nest and ambient temperature in F. aquilonia populations and populations that are admixed between F. aquilonia and the F. polyctena/F. rufa clade. and extend northwards in the future following a warming climate, as has been observed in other taxa (Antão et al., 2022). Global warming is markedly faster in Fennoscandia than the global average (Masson-Delmotte et al., 2021), yet biodiversity responses may be complex and not necessarily follow directly climatic changes (Waldock et al., 2018).

| Hybridisation patterns are consistent with the mating strategies of the different species
Our findings are among the few that pertain to a mosaic hybrid zone involving multiple species. We add a new system to the groups of taxa exhibiting natural mosaic hybrid zones (Abbott & Brennan, 2014;Fraïsse et al., 2014;Rieseberg et al., 1999) and multispecies hybridisation (Banerjee et al., 2022;Grant & Grant, 2020;Heliconius Genome Consortium, 2012;Natola et al., 2022;Ottenburghs, 2019;Reutimann et al., 2020).  . This pattern could be due to the life history characteristics of the species and/or selection. The simplest explanation is that F. polyctena or F. rufa queens are forced to mate with F. aquilonia males due to a lack of conspecifics, because of differences in their abundance in Finland (Seifert, 2018), similarly to swordtail fish three-way hybridisation in Mexico (Banerjee et al., 2022). A complementary explanation is that since highly supercolonial F. aquilonia mates on top of or close to its nest (Seifert, 2018), workers protect the queens from nonconspecifics, as has been observed in other ant species (Eyer et al., 2023;Staab & Kleineidam, 2014).
However, F. aquilonia males could disperse further as reported for polygynous F. lugubris (Gyllenstrand & Seppä, 2003) (Seifert, 2019a), Lasius emarginatus × L. platythorax (Seifert, 2019b) and Myrmica scabrinodis × M. vandeli (Yazdi et al., 2012) or, in the case of gyne samples, Lasius umbratus × L. merdionalis (Seifert, 2006).   (Sorvari et al., 2011), or hamper the ants in heating their nests in early spring (Seifert, 2018).  (Frouz & Finer, 2007). This may make ambient spring temperatures less crucial for ant survival than winter temperatures. One nonadaptive explanation for warmer hybrid nests would be differences in colony sizes, which in turn affect colony thermoregulation (Kadochová et al., 2019). We were unable to control for colony size here, but based on field observations this does not markedly differ between our hybrid and F. aquilonia populations. Overall, these results suggest that ambient temperature fluctuations might affect the wood ants' hibernation and reproductive success. It will be of interest in future studies to determine whether and to what extent thermal differences detected between F. aquilonia and hybrid nests are caused by either differences in the worker population size, ambient temperature, relative sun exposure, nest structure, and/or ant activity patterns.

| Hybrids occur in warmer and drier
The longer term adaptive potential of the hybrids in a changing environment is further supported by our finding of more heterozygosity in hybrids compared to the nonadmixed individuals.
Genetic diversity arising from hybridisation may fuel adaptation (Torda & Quigley, 2021) even in the presence of hybrid incompatibilities  and could provide means to adapt to new environmental challenges (Grant & Grant, 2019;Seehausen, 2013).

| CON CLUS IONS
We used combined genomic and morphological approaches to analyse species divergence and hybridisation in natural populations of extensively hybridising (Seifert, 2021), evolutionarily young (Goropashnaya et al., 2004(Goropashnaya et al., , 2012 (Antão et al., 2022). This is because we show that the hybrid populations are abundant in the landscape, are presumably longlived, have increased genetic diversity, and tend to live in warmer habitats than the sympatric nonadmixed F. aquilonia populations.
Since the majority of our samples were collected over 15 years ago, repeated sampling of the same and novel locations could provide insights into the evolution of the mosaic hybrid zone and pinpoint the role of climate as a potential driver of its evolution. Future research is needed to examine the relative roles and patterns of barriers to gene flow between the wood ant species, as well as patterns of introgression, and the consequences of hybridisation for the wood ants and the surrounding ecosystem.  which improved the quality of our manuscript.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interest.

O PEN R E S E A RCH BA D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/ dryad.4f4qr fjhf.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genetic data (filtered nuclear and mitochondrial SNP datasets used in this study), temperature and precipitation data (broadscale climatic data on temperature and precipitation, and smaller-scale within-nest and ambient temperature data), and related scripts have been archived in Dryad at https://doi.org/10.5061/dryad.4f4qr fjhf.