Landscape and Climatic Variations Shaped Secondary Contacts amid Barn Owls of the Western Palearctic.

The combined actions of climatic variations and landscape barriers shape the history of natural populations. When organisms follow their shifting niches, obstacles in the landscape can lead to the splitting of populations, on which evolution will then act independently. When two such populations are reunited, secondary contact occurs in a broad range of admixture patterns, from narrow hybrid zones to the complete dissolution of lineages. A previous study suggested that barn owls colonized the Western Palearctic after the last glaciation in a ring-like fashion around the Mediterranean Sea, and conjectured an admixture zone in the Balkans. Here, we take advantage of whole-genome sequences of 94 individuals across the Western Palearctic to reveal the complex history of the species in the region using observational and modeling approaches. Even though our results confirm that two distinct lineages colonized the region, one in Europe and one in the Levant, they suggest that it predates the last glaciation and identify a secondary contact zone between the two in Anatolia. We also show that barn owls recolonized Europe after the glaciation from two distinct glacial refugia: a previously identified western one in Iberia and a new eastern one in Italy. Both glacial lineages now communicate via eastern Europe, in a wide and permeable contact zone. This complex history of populations enlightens the taxonomy of Tyto alba in the region, highlights the key role played by mountain ranges and large water bodies as barriers and illustrates the power of population genomics in uncovering intricate demographic patterns.

Introduction

Species distribution patterns fluctuate in response to climatic variations, as populations relocate to follow their shifting niches (Huntley and Webb 1989). When organisms colonize new areas, obstacles in the landscape may lead populations to split with varying degrees of geographic isolation. Evolution, via mutation, drift, local adaptation, and gene flow, will then act independently on each of the isolated populations. If two allopatric populations are later geographically reconnected, after a certain amount of time and divergence, it can create a secondary contact zone. For example, when populations on both sides of an obstacle meet at the end of it in a ring-like fashion, as described in birds (Irwin et al. 2001), amphibians (Devitt et al. 2011; Pereira et al. 2011), or plants (Cacho and Baum 2012). Likewise, climate oscillations can lead to cyclical isolation and secondary contacts between populations as regional suitability varies (Stöck et al. 2012).

This complex interplay between climatic variations and landscape has been extensively studied in the Western Palearctic (Hewitt 1999, 2000), specifically in light of the cycles of glacial and interglacial periods that characterized the Quaternary (Willeit et al. 2019). During the last glaciation, colder temperatures and the expansion of the ice sheets in the north rendered large areas unsuitable for many species, which led them to follow their niches southward. Species found refuge in the Mediterranean peninsulas and in northern Africa where climatic conditions were more amenable, forming isolated populations. At the end of the last glacial maximum (LGM) (i.e., approximatively 20k years ago; Clark et al. 2009), this process was reversed as the climate warmed and the melting of the continental ice caps exposed free land that could be recolonized. An extensive literature addressing the postglacial history of European organisms describes how the complex landscape of the continent, combined with the distribution of species during the glaciation, conditioned their recolonization processes (reviewed in Hewitt [2000]). However, the low-resolution of genetic data used before the genomic era was often insufficient to resolve the intricate and often fine-scale evolutionary processes that occurred during recolonization.

Rapid development of high-throughput sequencing technologies and corresponding methodological tools during the last decades has opened new avenues to study natural populations with high precision. In particular, it has allowed biologists to reconstruct the evolutionary history of species and highlight the diversity of processes acting when populations or subspecies interact in secondary contact. These processes have been found to result in a variety of situations. Although the prolonged isolation of populations may lead to allopatric speciation (many examples in plants, Boucher et al. 2016; Tomasello et al. 2020); amphibians, Devitt et al. 2011; insects, Capblancq et al. 2015; mammals, Hey 2010; and birds, Linck et al. 2020), secondary contact tends to show a broad range of admixture patterns. When admixture occurs, it may vary from narrow hybrid zones between lineages (Poelstra et al. 2014), to the complete dissolution of a lineage (Rougemont et al. 2020), through a gradual level of admixture along a gradient of mixing populations (Duranton et al. 2018).

Microsatellite and mitochondrial data suggested that the barn owl (Tyto alba), a nonmigratory raptor, colonized the Western Palearctic in a ring-like fashion around the Mediterranean Sea after the last glaciation (Burri et al. 2016). Under this scenario, a postglacial expansion from the glacial refugium in the Iberian Peninsula to northern Europe formed the western branch of the ring, with the eastern branch present across the Levant and Anatolia. Although these observations led to conjecture of a potential admixture zone in the Balkans, the available genetic data at the time combined with the overall low genetic differentiation in this species did not allow to fully resolve this question. Moreover, the peculiar genetic makeup of populations in the presumed contact zone brought into question the possibility of a cryptic glacial refugium in the eastern Mediterranean peninsulas.

The difficulty in resolving the postglacial expansion of barn owls is mirrored by their convoluted taxonomy in the Western Palearctic. In this region, Tyto alba is classified into different subspecies based on geography and plumage coloration. First, T. a. erlangeri (Sclater, WL, 1921) reported in Cyprus and Middle East, may match the Levant lineage. Second, T. a. alba (Scopoli, 1769) is white-colored and supposedly present in western Europe and western Canary Islands, and could represent the western arm of the ring colonization. T. a. guttata (Brehm, CL, 1831), the third subspecies is a dark rufous morph allegedly found in Northern and Eastern Europe, in the Balkans and around the Aegean Sea. This taxonomy does not match any known genetic lineage identified so far and overlaps with the area where admixture between the two lineages that colonized Europe supposedly happens, making the presence of a subspecies in this area puzzling in light of the history known so far.

Here, taking advantage of whole-genome sequences of 94 individuals from all around continental Europe and the Mediterranean Sea, we elucidate the demographic history of barn owls in the Western Palearctic. Combining descriptive and modeling approaches based on genomic and ecological data, we identify how the climatic variations and landscape of the region shaped the history of this species. We also investigate how previously isolated populations of barn owls interact at secondary contacts between different lineages, and discuss the convoluted taxonomy with regards to their history.

Results

History of Barn Owls around the Mediterranean Sea

Genetic Diversity and Population Structure in the Western Palearctic

Despite an overall low differentiation (overall FST = 0.047), comparable with the overall FST = 0.045 estimated using microsatellites by Burri et al. (2016), the data set revealed a structuration of the genetic diversity among barn owls of the Western Palearctic. The first axis of the genomic principal component analyses (PCA) (explaining 3.32% of the total variance) contrasted individuals from the Levant populations (Israel—IS and Cyprus—CY) to all other individuals (fig. 1), consistent with K = 2 clusters being the best estimate in sNMF (supplementary figs. S4 and S5). For K = 3 (fig. 1), the Canary population (WC) formed an independent genetic cluster, and this was confirmed by the second axis of the PCA (explaining 2.6% of the variance) opposing it to all other individuals (fig. 1 and supplementary fig. S2). This isolation of WC was also observable in table 1, with a lot of private and rare alleles in Tenerife Island, it is higher FIT and the highest population-specific FST of all sampled populations. On the same PCA (fig. 1), individuals from European populations (France—FR, Switzerland—CH, Denmark—DK, Italy—IT, Serbia—SB, Greece—GR, Aegean Islands—AE) formed a third distinct cluster, matching their grouping in a single cluster at K = 3 with sNMF (fig. 1). The Iberian individuals (Portugal—PT) occupied a central position on the PCA (around 0 on both axes) and a mixed composition in sNMF (fig. 1). This central position of the Iberian population was also visible in the pairwise FST, where the highest value in the pairs involving PT is 0.055 (with both CY and WC), whereas all other pairwise comparisons involving populations from two of the distinct groups identified before (Levant, Canary Islands, and Europe) have values equal or higher (fig. 1). The TreeMix analysis (fig. 1) was also consistent with these results, since it also identified these major lineages, isolating the Canary population in a specific lineage, then grouping the Levant populations (IS and CY) in a second lineage and finally all European populations in a third lineage. In Europe, the Iberian population was basal to all other populations.

Fig. 1.

Genetic structure of barn owl populations in Western Palearctic. (a) Population structure for K = 6. Pie charts denote the individual proportion of each of lineages as determined by sNMF and are located at the approximate centroid of the sampled population. (b) Population structure for K = 3. Each bar denotes the individual proportion of each of the three lineages as determined by sNMF. (c) Matrix of pairwise FST between barn owl populations in Western Palearctic. The heatmap provides a visual representation of the FST values given in each cell. (d) PCA based on full set of 94 individuals. Point shape denotes populations and colored circles enclose sample clusters observed in sNMF (K = 3). Values in parentheses indicate the percentage of variance explained by each axis. (e) PCA based on of the 66 European individuals. (f) Population tree and the first migration event in Western Palearctic populations inferred by Treemix, rooted on WC for representation.

Table 1.

Population Genetic Diversity, Inbreeding, and Divergence Estimates for 11 Populations of Barn Owls from the Western Palearctic.

Pop	N	#Pl	#PA	#Rare	F IT	F IS	Pop FST
WC	9	2,217,235 (37,007)	123,016 (2,465)	407,403 (10,568)	0.067 (0.054)	−0.022 (0.057)	0.116 (0.006)
PT	9	2,639,343 (22,857)	102,252 (2,305)	539,157 (10,106)	−0.018 (0.042)	−0.008 (0.042)	0.003 (0.003)
FR	4	2,151,627 (0)	30,290 (581)	287,746 (1,468)	0.039 (0.124)	0.043 (0.131)	0.050 (0.005)
CH	10	2,494,462 (10,723)	47,532 (730)	386,654 (3,117)	0.025 (0.019)	−0.011 (0.019)	0.036 (0.002)
DK	10	2,410,615 (15,558)	40,650 (1,378)	349,800 (5,434)	0.026 (0.020)	−0.02 (0.021)	0.049 (0.002)
IT	9	2,404,069 (7,267)	66,297 (1,638)	401,842 (4,483)	0.035 (0.012)	−0.022 (0.012)	0.052 (0.005)
SB	5	2,336,060 (0)	32,096 (1,515)	326,906 (2,519)	0.025 (0.011)	−0.038 (0.011)	0.056 (0.004)
GR	9	2,454,653 (7,365)	44,996 (1,146)	378,422 (4,152)	0.018 (0.028)	−0.016 (0.027)	0.039 (0.002)
AE	10	2,460,422 (20,650)	56,260 (3,439)	403,259 (10,465)	0.018 (0.060)	−0.001 (0.062)	0.030 (0.002)
CY	10	2,338,318 (48,463)	113,377 (2,229)	480,021 (13,581)	0.022 (0.047)	−0.034 (0.049)	0.059 (0.004)
IS	9	2,509,099 (9,500)	172,624 (4,018)	608,944 (3,597)	−0.019 (0.015)	−0.036 (0.016)	0.021 (0.003)

Note.—SDs of the mean are provided between brackets for each parameter, see Materials and Methods section for details. N, number of individuals in the population; #Pl, number of polymorphic sites per populations; #PA, number of private alleles per population; #Rare, number of rare alleles per population; FIT, mean individual inbreeding coefficient relative to the meta-population; FIS, population level inbreeding coefficient; FST, population-specific FST as in Weir and Goudet (2017). Populations: WC, Canary Islands; PT, Portugal; FR, France; CH, Switzerland; DK, Denmark; IT, Italy; SB, Serbia; GR, Greece; AE, Aegean Islands; CY, Cyprus; IS, Israel.

Population Structure in Continental Europe

Focusing only on European samples, southern populations (PT, IT, GR, and AE) harbored a higher genetic diversity than northern populations (CH, FR, DK, SB) (table 1). The first axis of the genomic PCA based on samples from European-only populations opposed the Italian individuals to all others (fig.1 and supplementary fig. S3). This isolation of Italian samples was also apparent in the pairwise FST within Europe, with all the largest values involving IT. These results were consistent with sNMF on all samples for K higher than 3, where IT individuals formed an independent genetic cluster (supplementary fig. S5, illustrated in fig. 1), as well as TreeMix, where the Italian population was the first to split and had the longest branch within European lineage (fig. 1). Consistently with the distribution of the individuals in the European PCA (fig. 1 and supplementary fig. S3), the ancestry coefficients in the sNMF analyses of K = 4 and K = 5 (supplementary fig. S5) revealed the genetic differentiation of northern populations (FR, CH, and DK) compared with the Italian one, also opposed in the first axis. Individuals from GR and AE individuals shared ancestry with both western and Italian population, in line with their central position along this first axis of the European PCA. For K = 5, a third European component was distinguished in the Aegean individuals. This component was the majoritarian in Greek samples with contributions of both northern and Italian component; and Serbian samples appeared as a mix of the northern and the Aegean component. The eastern lineage (CY and IS), grouped at previous K, were split into two distinct ancestry pools at K = 6 (fig. 1). AE individuals harbored low amounts of CY ancestry, absent in all other European populations, and two CY individuals carried a large contribution of the Aegean component. Consistently, the first migration event detected by TreeMix was from CY, a population from the Levant lineage, to AE, a population from the European lineage (fig. 1 and supplementary fig. S6).

Fine Structure and Haplotype Sharing

The clustering of individuals by FineStructure, based on shared haplotypes between individuals was consistent with previous results (supplementary fig. S7). Individuals from the different populations sampled clustered together, except for CH and FR individuals, mixed in the same population. Consistently with this grouping, haplotypes from any given population were more likely to be found in individuals from the same population, followed by its most related populations (fig. 2). In the Levant lineage, IS haplotypes mostly painted IS individuals but also CY individuals and vice versa. Iberian haplotypes mostly painted PT individuals, but also contributed greatly to the painting of all European individuals, decreasing with distance. Western European haplotypes (from FR and CH) mostly painted western European individuals, followed by northern individuals (from DK) and finally eastern individuals (from SB, GR, and AE). The reverse pattern was observed for haplotypes from eastern Europe (GR, AE), with a gradient of contribution decreasing from east to west. Haplotypes from DK and SB mostly painted individuals from their own population, but also their respective geographical neighbors in both eastern and western European populations. Italian haplotypes were the most distinct haplotypes among European populations, mostly painting Italian individuals, followed by Greek individuals, and being painted by other populations at a lower rate than expected given its geographic position. Finally, AE haplotypes also painted more often CY individuals than IS individuals. This painting of levant individuals by AE haplotypes was higher than the contribution from any other European individual, and both CY and IS haplotypes painted more AE individuals than any other European individual.

Fig. 2.

Individual haplotype sharing between barn owl populations. Part of the total length of ChromoPainter chunks inherited from other genomes. Each graph summarizes the information of all the genomes from a given population, indicated on the top-left corner. Background colors match the lineages identified in figure 1.

Modeling of the History of European Barn Owls

Species Distribution Modeling

Habitat suitability projections showed that, from a climatic point of view, there were suitable regions for barn owls all around the Mediterranean Sea during the glaciation (20,000 years BP; fig. 3). Large areas were suitable in northern Africa and the Iberian Peninsula, but also in the two eastern Mediterranean peninsulas (current Italy and Greece). At this point, the sea levels were lower than today’s and the two eastern peninsulas were more connected, allowing for a continuous region of suitable barn owl habitat. At the mid-Holocene (6,000 years BP), major changes in sea level revealed a coastline very similar to nowadays. Our projections revealed a reduction of habitat suitability in northern Africa at this time, whereas the suitability of western and northern Europe increased (fig. 3). Finally, today, nearly all continental Europe is suitable for the barn owls, with the notable exception of mountain areas (fig. 3).

Fig. 3.

Modeling of the history of the barn owl in Europe. (a) Schematic representation of the five demographic scenarios tested for the colonization of the Europe by barn owls. Three models included one refugium in the Iberia during the LGM, whereas the last two included two refugia, one in Iberia and the second in Italy. Gray bars with snowflakes represent the last glaciation. (b) Best supported demographic model for the history of European barn owl populations as determined by fastsimcoal2. Time is indicated in thousands of years, determined using a 3-year generation time, CIs at 95% are given between brackets. Population sizes (haploid) are shown inside each population bar; arrows indicate forward-in-time number of migrants (2Nm) and direction. (c) SDM of barn owls based on climatic variables, projected into the past (last glacial maximum, 20 ka; mid-Holocene, 6 ka) and today’s condition. Locations in dark gray were highly suitable in 90% of the models. Below that threshold cells were considered as unsuitable (lightest gray shade on the graph). The present coastline is outlined in blue in all graphs.

Demographic Inference

Considering the position of Italy on the PCA (fig. 1), its high pairwise FST (fig. 1) and lower haplotype sharing with the rest of European populations (fig. 2) and its suitability for the barn owl during LGM (fig. 3), we tested whether the Italian peninsula could have been a cryptic glacial refugium during the last glaciation. We thus modeled five different demographic scenarios of recolonization of Europe using fastsimcoal2 (Excoffier and Foll 2011; Excoffier et al. 2013). Three models included only one refugium in the Iberian Peninsula and various possibilities of colonization scenarios, whereas two models included two refugia during the LGM, a western refugium in the Iberian Peninsula and an eastern refugium, in the Italian Peninsula (fig. 3, supplementary fig. S8, see Materials and Methods for details).

Akaike’s information criterion (AIC) and raw likelihood comparisons showed that the two refugia model 2R-1 explains best the site frequency spectra (SFS) of our data set (supplementary table S4 and fig. 3). In this model, the ancestral Italian lineage (IT) split from the Iberian lineage (PT) before the last glaciation, estimated at approximately 69,000 years BP (95% CI: 24,000–90,000 years BP; calculated with 3-year generation time). After its initial expansion, the ancestral population is estimated to have been larger in the Italian Peninsula than in Iberia (respectively, 189K (11K–320K) and 11K (10K–73K) haploid individuals), although the credibility intervals overlap. During the LGM (fixed between 24K and 18K years BP), both populations experienced a bottleneck, with a population size reduced to 6.8K (1.2K–141K) individuals in the Iberian lineage and 62 (36–116K) in the Italian. After the glaciation, the size of both populations increased to their current size, estimated at 44K (20K–380K) in the Iberian Peninsula and 1.3K (1K–326K) in the Italian Peninsula, both smaller but consistent with their estimated census size (55K–98K for the Iberian, Lourenço et al. 2015; Martí and del Moral 2003 and 6K–13K for the Italian, Pierandrea and Giancarlo 2005). The Greek population split from the Italian branch around 5,700 years BP, whereas the Swiss (CH) population split slightly later from Iberia (5,000 years BP) and maintain a high level of gene flow (estimated to 90 (46–1.4K) from CH to PT and 8 (3–62) in the reverse direction). Current effective population sizes of the CH and GR populations are estimated to 3.4K (1K–205K) and 1.4K (1K–208K), respectively (fig. 3), in line with census results (200–1,000, Knaus et al. 2018 and 3,000–6,000, BirdLife International 2004, respectively). Migration between these populations is estimated to be highest from IT and GR to CH (respectively, 27 (0.2–157) migrants from IT and 42 (0.1–156) from GR) and lowest in the opposite direction (respectively, 1.3 (0.1–96) migrants from CH to GR and 0.02 (0.2–38) from CH to IT) (supplementary table S5). Point estimates with 95% confidence intervals (CI) for all parameters of the best model are given in (supplementary table S5), as well as single point estimates for all models (supplementary table S3).

Barriers and Corridors

Migration Surface Estimate in the Western Palearctic

Estimated Effective Migration Surface identified large water bodies, especially in the eastern Mediterranean and around Cyprus, as regions resisting to gene flow (supplementary fig. S12). On the mainland, barriers to gene flow matched the main formations of the Alpide belt in the region, an orogenic formation spanning from western Europe to eastern Asia. From west to east, a light barrier overlapped with the Pyrenees, a strong barrier spanned the Alps to the Balkans, and a third obstacle matched the Taurus mountains in Anatolia. A region with high gene flow was identified in continental Europe above the Alps, spanning from western Europe to the Balkan Peninsula.

Isolation by Distance in Europe

In continental Europe, the shortest path overland did not correlate significantly with genetic distance (supplementary fig. S13) (mantel test, P = 0.193, R = 0.20/linear model, P = 0.26, R2 = 0.012). On the contrary, when the geographic distance between populations included the barrier formed by the Alps (i.e., the Italian population was connected to other populations via the Greek Peninsula, itself connected to western Europe via northern Europe), both tests were significant (mantel test, P = 0.002, R = 0.68/linear model, P = 1.3×10−5, R2 = 0.507).

Discussion

The history of natural populations is shaped by the combination of landscape barriers and climatic variations that isolate and mix lineages through their combined actions. Consistently with previous work using microsatellites (Burri et al. 2016), we show that barn owls colonized the Western Palearctic in a ring-like fashion around the Mediterranean Sea, with one arm around the Levant and a second throughout Europe. However, using whole-genome sequences, we found this colonization actually predates the last glaciation and pinpoint a secondary contact zone between the two lineages in Anatolia rather than in the Balkans. In addition, we provide evidence that barn owls recolonized Europe after the LGM from two distinct glacial refugia—a western one in Iberia and newly identified eastern one in Italy. As temperatures started rising, western and northern Europe were colonized by owls from the Iberian Peninsula whereas, in the meantime, the eastern refugium population of Italy had spread to the Balkans (fig. 4). The western and eastern glacial populations finally met in eastern Europe. This complex history of populations questions the taxonomy of the multiple T.alba subspecies, highlights the key roles of mountain ranges and large water bodies as barriers to gene flow for a widespread bird, and illustrates the power of population genomics in unraveling intricate patterns.

Fig. 4.

Schematic representation of the history of barn owls in the Western Palearctic and the main barriers in the region. Orange arrows depict the colonization of the region by the three main lineages (Levant, Canary Islands, European). Yellow arrows represent the modeled postglacial recolonization scheme of Europe, with two distinct refugia (yellow dots). Blue lines represent the main barriers identified in this work, namely the Alps in Europe (dashed line) and the Taurus and Zagros mountains in Anatolia (solid line).

Colonization of the Western Palearctic and Gene Flow in Anatolia

Our results show that two distinct barn owl genetic lineages surround the Mediterranean Basin: one in the Levant, and a second in Europe (fig. 1), likely connected via northern Africa. Supported by the higher and specific diversity of the basal population of each arm (namely, Israel—IS and the Iberian peninsula—PT; table 1), these observations are consistent with the ring colonization scenario hypothesized by Burri et al. (2016). However, we show that a barn owl population survived the last glaciation in Italy (see next section for details) and that its genetic makeup resembles the European lineage (fig. 1). Therefore, the ring colonization of Europe around the Mediterranean appears to have been preglacial, whereas the postglacial history is more convoluted (see next section).

In previous studies, the ancestry of Greek and Aegean populations was unclear, with a hypothesized mixed origin between the European and Levant lineages (Burri et al. 2016). This uncertainty was mostly likely due to the low resolution of genetic markers (Mitochondrial DNA and microsatellites), as the genomic data reported here clearly show that Greek and Aegean owls are genetically much closer to European than to Levant ones (fig. 1). This observation indicates that the European lineage reached further east than previously assumed, allowing us to pinpoint the secondary contact zone between the European and Levant lineages to Anatolia, instead of the Balkans as it had been proposed. In Anatolia, the Taurus and Zagros mountain-ranges form an imposing barrier that appears to have stopped the expansion of the Levant lineage both during the ring colonization and nowadays.

Despite the barrier, and although we do not see a complete admixture of the two lineages, there is evidence of some gene flow. Indeed, the first migration in Treemix (fig. 1) pointed to a secondary contact between Cypriot and Aegean populations, consistent with the signals of admixture between them (figs. 1, 1, and 2). The admixture pattern however is restricted geographically to this narrow region and does not permeate further into either of the lineages, as surrounding populations (Israel in the Levant, and Greece and Serbia in Europe) do not show signals of admixture (figs. 1 and 2). Thus, the migration between populations on both sides seems limited and possibly only occurs along a narrow corridor along the Turkish coast where only a few barn owls have been recorded (Göçer and Johnson 2018). Further analyses with samples from Anatolia should allow to characterize in high resolution how and when admixture occurred in this region.

Glacial Refugia and Recolonization of Europe

Previous studies showed that barn owls survived the last glaciation by taking refuge in the Iberian Peninsula and maybe even in emerged land in the now immerged Bay of Biscay (Antoniazza et al. 2014; Machado et al. 2021). The observed distribution of diversity in Europe, and especially the specific makeup of the Italian populations, is best explained by a demographic model with two glacial refugia—one in Iberia and a second in Italy, derived from the Iberian population before the glaciation (fig. 3; supplementary table S6). Environmental projections not only support this model, as Italy was highly suitable for the species at the time, but also show that, due to the low levels of the Adriatic Sea, the suitable surface extended to the west coast of the Balkans (LGM—fig. 3). Crucially, three of the barn owl’s key prey also had glacial refugia in the Italian and Balkan peninsulas, namely the common vole (Microtus sp.) (Stojak et al. 2019), the wood mouse (Apodemus sp.) (Herman et al. 2016), and shrews (Crocidura sp.) (Dubey et al. 2007).

The inferred size of the preglacial population of barn owls inhabiting current Italy was larger than any other (189K; 11–320K), prior to a very strong bottleneck during the glaciation (population reduced to 62 individuals; 36–116K). These extreme values should be considered with caution as they are likely a reflection of the limitations of our modeling approach which itself may have broader implications. Firstly, parameter estimation can be affected by modeling simplifications, such as constant population sizes and instant bottlenecks. Secondly, although we attempted to remove the influence of direct selection by removing coding regions from the data set, we might have been unable to account for linked selection which can affect noncoding regions of the genome and impact the results of the inference (Pouyet et al. 2018). Finally, the software (fastsimcoal2) cannot account for potential gene flow from unmodeled or unsampled populations. For example, we purposefully excluded the Levant lineage from the models, in an effort to reduce the computational burden of what would likely be an unreliable estimation considering the longtime of divergence of this population. However, there is evidence of gene flow from the levant lineage into the European populations, and its absence from our modeling could have impacted parameter estimates. Further, unsampled regions may have contributed to the genetic composition of the modeled populations. For example, gene flow from northern Africa to Italy via Sicily might have taken place, facilitated by the increased connectivity between Italy and North Africa during the glaciation (fig. 3) (Dapporto and Bruschini 2012; Husemann et al. 2014; Lodolo et al. 2020).

With the warming following the LGM, Europe became gradually more suitable and, by the mid-Holocene (6,000 years BP), most of western and northern Europe were appropriate for barn owls (fig. 3) as well as the common vole (Microtus arvalis) (Heckel et al. 2005). The genetic similarity between Iberian (PT) and northwestern populations (CH, FR, DK; fig. 1) indicates that barn owls colonized these newly available regions from the Iberian refugium as previously thought (Antoniazza et al. 2014). The contribution from the Italian refugium to northern populations appears to have been hindered by the Alps (see next sections), as suggested by the higher genetic distance between them (figs. 1, 1, 1, and 2). Instead, at this time the Adriatic Sea had neared today’s levels, isolating genetically and geographically the Italian refugium from its component in the Balkan Peninsula (fig. 3IT-GR split ∼6K years BP fig. 3). Only more recently did the rise of temperatures allow for areas in the east of Europe to become suitable, finally connecting the southeastern populations near the Aegean Sea (GR and AE) with populations in the northeastern part of Europe (fig. 3). In particular, the high heterozygosity and admixed ancestry of Serbian individuals (table 1and fig. 1) suggest that the suture between the Iberian and the Italo-Greek glacial lineages took place in eastern Europe. This newly identified postglacial recolonization scheme of continental Europe by the barn owl matches the general pattern described for the brown bear (Ursus arctos) and shrew (Sorex sp.) (Hewitt 1999).

The isolation by distance pattern observed between European populations (supplementary fig. S13) highlights a diffusion of alleles in the European populations rather than a narrow hybrid zone (Bertl et al. 2018). Further, the inferred migration rates support high current gene flow in the region (fig. 3 and CH and GR in supplementary table S5), and we found signals of each ancestry in populations far from the suture zone (dark blue in GR and light blue in CH and DK, fig. 1). Finally, the measure of haplotype sharing decreases consistently with distance between populations around the northern side of the Alps between populations from Iberia to Greece, excluding Italy (fig. 2). Surrounded by the sea and the Alps, Italy is the exception and appears to have avoided incoming gene flow (supplementary fig. S12, isolation in fig. 2), thus being a better-preserved relic of the refugium population (own cluster in sNMF K > 4, supplementary fig. S5). In contrast, the Balkan component admixes smoothly with the other European populations. Such seamless mixing of the two glacial lineages from southern refugia where barn owls are mostly white (Antoniazza et al. 2010; Burri et al. 2016), brings further into question the subspecies T. a. guttata.

The Case of Tyto alba guttata

Traditionally, in Europe, the eastern barn owl (T. a. guttata) is defined by its dark rufous ventral plumage in contrast to the white western barn owl (T. a. alba) (Mátics et al. 2001). With a wide distribution, the former is recorded from The Netherlands to Greece, including most of northern and eastern Europe (Clements et al. 2019). However, this repartition does not match the history of any specific glacial lineage identified above, nor any genetically differentiated population. The dark populations of northern Europe (DK) are genetically as similar to lighter western populations (FR, CH) than to the dark ones in the east (SB). This color variance within European populations has been shown to be maintained through local adaptation (Antoniazza et al. 2010), and while the genomic basis and history of this trait remain worthy of future investigations, we suggest that all European barn owls form a single subspecies (T. a. alba), reflecting the entire European population, regardless of their color.

Barriers and Corridors Shape the Connectivity of the Western Palearctic Meta-Population

The partition of genetic diversity among barn owls in the Western Palearctic allowed us to identify barriers and corridors to gene flow. Populations isolated by large water bodies have accumulated substantial genetic differences as, for example, the higher FST in the Canary and Cyprus islands (fig. 1 and table 1), and reflect the importance of water as a barrier to dispersion in this species (Machado et al. 2021). On the mainland, and as described for the American barn owl (Machado et al. 2018), major mountain ranges act as significant obstacles to migration for European barn owls and can generate genetic structure. First, the high mountain ranges of Taurus and Zagros coincide with the contact zone between the Levant and the European lineages both nowadays and potentially at the time of the preglacial ring colonization of Europe (see above). Second, the Alps and the Balkan Mountains slowed the northward expansion of the glacial populations of Italy and Greece after the LGM and still constrain migration between populations on both sides of their ranges. If these results remain to be confirmed with observational data (i.e., ringing data not available for all countries), they emphasize that, despite its worldwide repartition and its presence on many islands, the connectivity of barn owl populations is heavily driven by biogeographical barriers.

Conclusion

The combination of whole-genome sequencing and sophisticated modeling methods has revealed and completed the evolutionary history of many species, starting with humans (Garrigan and Hammer 2006; Auton et al. 2009) and domestic species (vonHoldt et al. 2010; Rendón-Anaya et al. 2017) often revealing intricate complexities in their past histories. Applied to nonmodel species, these techniques have also shown that theoretical models such as the ring speciation theory by Mayr (1942) can only partially explain the observed data. For example, populations of the western fence lizard (Sceloporus occidentalis), that spread in a ring-like manner around the Sierra Nevada, actually show strong evidence for gene flow at the northern terminus of the ring where they were thought to have become isolated gene pools (Bouzid et al. 2021). Similarly, genomic analyses of Greenish warblers around the Himalayan plateau showed that admixture also occurs at the end of the ring and historical breaks in gene flow have existed at multiple location around the ring (Alcaide et al. 2014). Here, we show that the barn owl follows a similar trend with a secondary contact zone at the end of the ring, discontinuities in the ring, as well as the discovery of a cryptic glacial refugium within the ring itself, thus completing the history of the species around the Mediterranean Sea. However, several questions remain unanswered, waiting for relevant samples to be collected and analyzed: what role did northern African populations play in connecting the Levant and European lineages? Did they contribute to the diversity observed in Italy? How narrow is the contact zone between the Levant and European lineages in Anatolia? Lastly, the origin of barn owls from the Western Palearctic as a whole also deserves further investigation, as they are believed to have colonized the Western Palearctic from the East, given the species’ supposed origin in Africa or southeastern Asia approximately 4 My BP (Aliabadian et al. 2016). But this is at odds with the higher genetic diversity of the Iberian population compared with the Levant one. Such inconsistency points to the need for samples from around the world, to understand how this charismatic group of nocturnal predators conquered the entire planet.

Moreover, research on postglacial recolonization and the subsequent phylogeographic patterns peaked at the turn of the century, with many studies providing an overview of the history of a wide variety of organisms (reviewed in Hewitt [1999]). The rise in availability of genomic data for these species will rewrite the history of many of them and will probably challenge paradigms about recolonization processes. Furthermore, such data will allow to detail the genomic consequences of the history; both from a neutral and selective perspective. Applied to several species, these approaches will redefine with greater clarity the broad phylogeographical patterns in the Western Palearctic and elsewhere, to rethink taxonomic classifications and to better understand how organisms might adapt to a changing environment in a complex, fragmented and rapidly changing landscape.

Materials and Methods

Samples and Data Preparation

Sampling, Molecular, and Sequencing Methods

The whole genomes of 96 individual barn owls (T.alba) were used in this study (supplementary table S1): 94 individuals were sampled in 11 Western Palearctic localities: Canary Islands (Tenerife Island—WC), Portugal (PT), France (FR), Switzerland (CH), Denmark (DK), Serbia (SB), Greece (GR), Italy (IT), Aegean islands (AE), Cyprus (CY), and Israel (IS). In addition, one Eastern (Tyto javanica from Singapore) and one American barn owl (Tyto furcata from California) were used as outgroups. Illumina whole-genome sequences of individuals from PT, FR, CH, DK, and the outgroups were obtained from the GenBank repository (BioProject PRJNA700797). For the remaining 61 individuals, we followed a similar library preparation and sequencing protocol as outlined in Machado et al. (2021). Briefly, genomic DNA was extracted using the DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany), and individually tagged. About 100 bp TruSeq DNA PCR-free libraries (Illumina) were prepared according to manufacturer’s instructions. Whole-genome resequencing was performed on multiplexed libraries with Illumina HiSeq 2500 PE high-throughput sequencing at the Lausanne Genomic Technologies Facility (GTF, University of Lausanne, Switzerland), producing a mean number of 121,529,696 paired-end reads per individuals (2×150 bp long, SD: 21,012,431 reads).

Data Processing, SNP Calling, and Technical Filtering

The bioinformatics pipeline used to obtain analysis-ready SNPs was adapted from the Genome Analysis Toolkit (GATK) Best Practices (Van der Auwera et al. 2013) to a nonmodel organism following the developers’ instructions, as in Machado et al. (2021). Raw reads were trimmed with Trimommatic v.0.36 (Bolger et al. 2014) and aligned to the reference barn owl genome (Machado et al. 2021) with BWA-MEM v.0.7.15 (Li and Durbin 2009). Base quality score recalibration (BQSR) was performed using high-confidence calls obtained from two independent callers—GATK’s HaplotypeCaller and GenotypeGVCF v.4.1.3 and ANGSD v.0.921 (Korneliussen et al. 2014)—as a set of “true variants” in GATK v.4.1.3 following the procedure described in Machado et al. (2021).

Genotype calls were filtered for analyses using a hard-filtering approach as proposed for nonmodel organisms, using GATK and VCFtools (Danecek et al. 2011). Calls were removed if they presented: low individual quality per depth (QD < 5), extreme coverage (1,100 > DP > 2,500), mapping quality (MQ < 40 and MQ > 70), extreme hetero or homozygosity (ExcessHet > 20 and InbreedingCoeff > 0.9), and high read strand bias (FS > 60 and SOR > 3). Then, we removed calls for which up to 5% of genotypes had low quality (GQ < 20) and extreme coverage (GenDP < 10 and GenDP > 40). We kept only biallelic sites, excluded SNPs on the sexual chromosome (super scaffolds 13 and 42; Machado et al. 2021) and an exact Hardy–Weinberg test was used to remove sites that significantly departed (P = 0.05) from the expected equilibrium using the package HardyWeinberg (Graffelman and Camarena 2008; Graffelman 2015) in R (R Core Team 2020), yielding a data set of 6,448,521 SNP (mean individual coverage: 19.99×; SD: 4.38). Lastly, we discarded singletons (minimum allelic count [mac] <2), yielding to a total of 5,151,169 SNP for the population genomic analyses.

SNP Phasing and Quality Control

The set of 6,448,521 variants was phased in two steps. First, individual variants were phased using a read-based approach in which reads covering multiple heterozygous sites were used to resolve local haplotypes. To do so, WhatsHap v1.0 (Martin et al. 2016) was run independently for each individual with default parameters. Secondly, variants were statistically phased with Shape-It v4.1.2 (Delaneau et al. 2019). This algorithm integrates local individual phase and applies an approach based on coalescence and recombination to statistically phase haplotypes and impute missing data. Shape-It was run following the manual instructions for a better accuracy: the number of conditioning neighbors in the PBWT was set to 8, and the MCMC chain was run with ten burn-in generations, five pruning iterations, each separated by one burn-in iteration, and ten main iterations.

To assess the quality of the phasing, we examined phase accuracy by using the switch-error-rate metric (Choi et al. 2018). When comparing two phasing for an individual’s variants, a switch error occurs when a heterozygous site has its phase switched relative to that of the previous heterozygous site. Thus, for each individual, we compared the true local phasing inferred form the read-based approach (WhatsHap) and the statistical phasing of this individual’s variants statistically phased by Shape-It, with read-based phase information ignored only for the individual considered (same version and parameters than in the paragraph above). The final estimation of the switch error rate was done using the switchError code to compare both phasing sets (available at https://github.com/SPG-group/switchError) (supplementary fig. S1).

Population Structure and Genetic Diversity

In order to investigate population structure among our samples, sNMF (Frichot et al. 2014) was run for a number of clusters K ranging from 1 to 10, with 25 replicates for each K to infer individual clustering and admixture proportions. For this analysis, SNPs were pruned for linkage disequilibrium in PLINK v1.946 (Purcell et al. 2007) (parameters –indep-pairwise 50 10 0.1) as recommended by the authors, yielding 594,355 SNP. Treemix (Pickrell and Pritchard 2012) was used to calculate a drift-based tree of our populations, using this LD-pruned data set. To detect admixture events between populations, ten Treemix replicates were run for zero to ten migration events, with the tree rooted on the WC population, representative of a nonadmixing population (see Results and supplementary fig. S6).

Population expected and observed heterozygosity, population-specific private alleles, population-specific rare alleles (mac ≤ 5), population-specific total number of polymorphic sites (table 1), population-level genetic diversity (pi and Watterson estimator—supplementary table S2) were estimated on the 5,151,169 variants data set using the package hierfstat v.0.5-10 (Goudet 2005). To account for differences in sample sizes, which ranges from 4 to 10, population-specific statistics were calculated by randomly sampling five individuals from the larger populations (all except FR and SB) ten times in a bootstrap-fashion and estimating the mean and SD. Individual-based relatedness (β) (Weir and Goudet 2017) and inbreeding coefficient for SNP data were calculated with the R package SNPRelate (Zheng et al. 2012). Overall FST, population pairwise FST and population-specific FST (Weir and Goudet 2017) were computed with the package hierfstat. CIs for population-specific FST were computed by dividing the SNPs into 100 blocs, and bootstrapping 100 times 100 blocks with replacement. Finally, PCA were also performed with the R package SNPRelate, first with all individuals and second only with the 66 European ones (excluding WC, CY, and IS).

Haplotype Sharing

To measure shared ancestry in the recent past between individuals, we ran fineSTRUCTURE (Lawson et al. 2012) on the phased data set including all individuals. For this analysis, we initially modeled haplotype sharing between individuals using ChromoPainter to generate a co-ancestry matrix, which records the expected number of haplotypes chunks each individual donates to another. For this ChromoPainter step, we converted phased haps files to ChromoPainter phase files using the impute2chromopainter.pl script provided at http://www.paintmychromosomes.com and generated a uniform recombination map with the makeuniformrecfile.pl script. Using the version of ChromoPainter built into fineSTRUCTURE v.2.0.8, we performed ten EM iterations to estimate the Ne and Mu parameters (switch rate and mutation rate). The model was then run using the estimated values for these parameters (respectively, 570.761 for Ne and 0.0074240 for Mu), and we used default settings to paint all individuals by all others (-a 0 0). We ran fineSTRUCTURE’s MCMC model on the co-ancestry matrix for 500,000 burn-in and 500,000 sampling iterations, sampling every 10,000 iterations to determine the grouping of samples with the best posterior probability.

Modeling of History of European Barn Owl

Maximum-Likelihood Demographic Inference

Data Preparation

To describe the history of barn owls in Europe, we modeled five different demographic scenarios using fastsimcoal2 (Excoffier and Foll 2011; Excoffier et al. 2013). Given the position of Italy on the PCA (fig. 1), its high FST and lower haplotype sharing with the rest of European populations (fig. 2), we tested in particular whether it could have been a cryptic glacial refugium during the last glaciation. To focus on European populations and due to computational constraints, we simplified the data set to model the history of four populations of eight individuals (supplementary table S1): PT as representatives of the known refugium in the Iberian Peninsula (Antoniazza et al. 2014); IT and GR representing the peninsulas of Italy and Balkans, respectively; and CH, a product of the recolonization of northern Europe from the Iberian refugium(Antoniazza et al. 2014). The populations from the Levant were excluded from this European-centered modeling because the inferred convergence between the two lineages would likely be unreliable for an event so far back in the past without calibrating events. Therefore, the computational burden associated with including the Levant lineage was deemed unnecessary.

The set of 6,448,521 autosomal SNPs was filtered to retain only neutrally evolving regions. We first excluded SNPs found in genic regions and CpG mutations (Pouyet et al. 2018). To achieve homogeneity among SNPs, we removed all sites with missing data and excluded positions with a coverage outside two-thirds of the SEM. We employed a parsimony approach based on the Tytonidae phylogenetic tree (Uva et al. 2018) to determine the ancestral state of the SNPs using the genomes of the two outgroups. Sites for which it was impossible to attribute a state based on the available outgroups were discarded. The remaining 770,718 SNPs were used to produce population pairwise SFS.

Demographic Scenarios and Parameters

Five different scenarios were tested to model the history of barn owls in continental Europe, with a special focus on the period since the last glaciation (fig. 3 and supplementary fig. S8). Three models included only one refugium in the Iberian Peninsula and various possibilities of colonization scenarios, thus excluding the persistence of barn owls in a second refugium during the glaciation. The two last models included two refugia during the LGM, a western refugium in the Iberian Peninsula and an eastern refugium, in the Italian Peninsula. Given that descriptive analyses clearly point to a certain degree of gene flow between all mainland populations, we did not model a null scenario with no migration. Instead we allowed the search of migration parameters to start at zero.

The models 1R-1, 1R-2, and 1R-3 included only one refugium in the Iberian Peninsula. 1R-1 model assumed only one colonization route around the north side of the Alps (forming the CH population), and from there move southeast to reach first the Balkans (GR) and then Italy (IT). The two other single-refugium models (1R-2 and 1R-3) assumed two distinct colonization routes, one north of the Alps to CH and the other south of the Alps along the Mediterranean coast to IT. 1R-2 assumes current Greece would have been colonized by owls from the Italian Peninsula, following the route along the Mediterranean coast. In 1R-3, Greece would have been colonized via northern Europe, whereas the Mediterranean expansion would have stopped in Italy.

The last two models included a second, eastern, refugium (2R-1 and 2R-2). In these models, the western lineage expansion from the Iberian population would have colonized Europe before the last glaciation, thus occupying all the Mediterranean peninsulas. During the last glaciation, two distinct populations would have survived, respectively, in the Iberian (western refugium) and Italian (eastern refugium) Peninsulas. Both models assume that northern Europe was recolonized from the Iberian lineage after the glaciation, but they differ in the scenario of recolonization of southeastern Europe. In 2R-1, Greece was recolonized from the Italian refugium, whereas in 2R-2 the expansion from the Iberian Peninsula would have recolonized all eastern Europe, including current Greece. In this last scenario, the Italian population would be the only relic from the eastern refugium. For all scenarios, migrations between populations were allowed (see fig. 3 and supplementary table S3).

Wide search ranges for initial simulation parameters were allowed for population sizes, divergence times, and migration rates (supplementary table S3). Each population split was preceded by an instantaneous bottleneck, in which the founding population size was drawn from a log-uniform distribution between 0.01 and 0.5 proportion of current population sizes.

Demographic Inference

Demographic simulations and parameter inference were performed under a composite-likelihood approach based on the joint SFS as implemented in fastsimcoal2 (Excoffier and Foll 2011; Excoffier et al. 2013). For each of the five scenarios, 100 independent estimations with different initial values were run. For each run, there were 500,000 coalescent simulations (option -n), with 50 expectation–maximization (EM) cycles (-M and -L). As we do not have an accurate mutation rate for barn owls, we fixed the end of the glaciation to 6,000 generations BP (∼18,000 years BP with a 3-year generation time) and scaled all other parameters relative to it using the -0 command option (using only polymorphic sites). The best-fitting scenario out of the five tested was determined based on AIC (Akaike 1998) and confirmed through the examination of the likelihood ranges of each scenario as suggested in Kocher et al. (1989). Nonparametric bootstrapping was performed to estimate 95% CI of the inferred parameters under the best-fitting scenario. To account for LD, a block-bootstrap approach was employed as suggested by the authors (Excoffier and Foll 2011; Excoffier et al. 2013): the SNPs were divided into 100 same-size blocks, and then 100 bootstrap SFS were generated by sampling these blocks with replacement. Due to computational constraints, for bootstrapping, we ran 50 independent parameter inferences per bootstrapped SFS with only ten EM cycles each, instead of 50 cycles used for comparing scenarios above. This procedure has been defined as conservative (Malaspinas et al. 2016), and is expected to produce quite large CIs. We accepted this trade-off as our main goal was to determine the best demographic topology, accepting uncertainty on specific parameter values. The highest maximum-likelihood run for each bootstrapped SFS was used to estimate 95% CI of all parameters.

Niche Modeling

In order to identify the regions of high habitat suitability for barn owls at the LGM (20,000 years BP) and to support the demographic scenarios tested in the previous section, we modeled the past spatial distribution of the species in the Western Palearctic. We built species distribution model (SDM) using Maximum Entropy Modeling (MaxEnt), a presence-only based tool (Elith et al. 2011). Current climatic variables for the Western Palearctic (supplementary fig. S10) were extracted from the WorldClim database at 5 arc min resolution using the R package rbioclim (Exposito-Alonso 2017), and filtered to remove variables with a correlation of 0.8 or higher (supplementary fig. S9). The variables retained were: Mean Diurnal Range (Bio2), Min Temperature of Coldest Month (Bio6), Temperature Annual Range (Bio7), Mean Temperature of Wettest Quarter (Bio8), Precipitation Seasonality (Bio15), Precipitation of Driest Quarter (Bio17), and Precipitation of Coldest Quarter (Bio19). We built models with linear, quadratic, and hinge features, and with a range (1–5) of regularization multipliers to determine which combination optimized the model without over complexifying it. The best combination based on the corrected AIC (as recommended by Warren & Seifert 2011) was achieved with a quadratic model with 1 as regularization multiplier (supplementary table S6). We ran 100 independent maxent models, omitting 25% of the data during training to test the model. To avoid geographic bias due to different sampling effort in the distribution area of the species, we randomly extracted 1,000 presence points within the IUCN distribution map (BirdLife International 2019) for each model run (Fourcade 2016).

Predictive performances of the models were evaluated on the basis of the area under the curve (AUC) of the receiver operator plot of the test data. For all models with an AUC higher than 0.8 (considered a good model; Swets 1988; Li et al. 2020), we transformed the output of Maxent into binary maps of suitability. We assigned a cell as suitable when its mean suitability value was higher than the mean value of the 10% test presence threshold. This conservative threshold allows us to omit all regions with habitat suitability lower than the suitability values of the lowest 10% of occurrence records. Finally, we averaged the values of the models for each cell, and only cells suitable in 90% of the models were represented as such in the map.

We projected the models to the climatic conditions of the mid-Holocene (6,000 years BP) and the LGM (20,000 years BP), which we extracted from WorldClim at the same resolution as current data. When projecting to past climates, the Multivariate Environmental Similarity Surface (MESS) approach (Elith et al. 2011) was used to assess whether models were projected into climatic conditions different from those found in the calibration data. As our goal was to highlight only areas of high suitability for barn owls, cells with climatic conditions outside the distribution used to build the model were assigned as unsuitable (0 attributed to cell with negative MESS). For each timepoint, the results of the models were merged and transformed into a binary map as described for current data (fig. 3).

Migration Surface Estimation in the Western Palearctic

The Estimated Effective Migration Surface (EEMS) v.0.0.9 software (Petkova et al. 2016) was used to visualize geographic regions with higher or lower than average levels of gene flow between barn owl populations of the Western Palearctic. Using the SNP data set pruned for LD produced above, we calculated the matrix of genetic dissimilarities with the tool bed2diff. The free Google Maps api v.3 application available at http://www.birdtheme.org/useful/v3tool.html (last accessd December 08, 2021) was used to draw the polygon outlining the study area in the Western Palearctic. EEMS was run with 1,000 demes in five independent chains of 5 million MCMC iterations with a burn-in of 1 million iterations. Results were visually checked for MCMC chain convergence (supplementary fig. S11) and through the linear relation between the observed and fitted values for within‐ and between‐demes estimates using the associated R package rEEMSplots v.0.0.1 (Petkova et al. 2016). With the same package, we produced a map of effective migration surface by merging the five MCMC chains.

Isolation by Distance in Continental Europe

To investigate how population structure correlated with spatial distances between European populations and to detail the role of the Alps as a barrier to gene flow, we performed Mantel tests as implemented in the ade4 package v.1.7-15 (Dray and Dufour 2007, p. 4) for R. We compared the genetic distances (pairwise FST between populations, see Population Structure and Genetic Diversity for details) with different measures of geographical distance between populations: the shortest distance over land via direct flight and the distance constrained by the presence of the Alps, forcing the connection of the Italian population to the other populations via the Greek Peninsula (supplementary fig. S13). We also tested the linear regression between both variables in R.

Click here for additional data file.