
==== Front
Sci Adv
Sci Adv
sciadv
advances
Science Advances
2375-2548
American Association for the Advancement of Science

adm8596
10.1126/sciadv.adm8596
Research Article
Earth, Environmental, Ecological, and Space Sciences
SciAdv r-articles
Ecology
Life Sciences
Genetics
The genetic legacy of the Quaternary ice ages for West Palearctic butterflies
The genetic legacy of ice ages for butterflies
https://orcid.org/0000-0001-7129-4526
Dapporto Leonardo Conceptualization Data curation Formal analysis Funding acquisition Investigation Methodology Project administration Resources Software Supervision Validation Visualization Writing - original draft Writing - review & editing 1 *†
https://orcid.org/0000-0002-0707-7495
Menchetti Mattia Conceptualization Data curation Investigation Methodology Project administration Resources Validation Visualization Writing - review & editing 2 †
https://orcid.org/0000-0003-1791-2148
Dincă Vlad Conceptualization Data curation Formal analysis Funding acquisition Investigation Resources Validation Writing - review & editing 3
https://orcid.org/0000-0003-1112-1345
Talavera Gerard Data curation Formal analysis Investigation Methodology Resources Software Validation Writing - review & editing 4
https://orcid.org/0000-0002-2419-2516
Garcia-Berro Aurora Data curation Formal analysis Investigation Methodology Resources Software Writing - review & editing 4
https://orcid.org/0000-0002-3089-7520
D'Ercole Jacopo Conceptualization Methodology Validation Writing - review & editing 5 6
https://orcid.org/0000-0002-3081-6700
Hebert Paul DN Data curation Funding acquisition Investigation Project administration Resources Writing - review & editing 5 6
https://orcid.org/0000-0002-2447-4388
Vila Roger Conceptualization Data curation Funding acquisition Investigation Methodology Project administration Resources Supervision Visualization Writing - original draft Writing - review & editing 2 *
1 ZEN Lab, Department of Biology, University of Florence, Florence, Italy.
2 Institut de Biologia Evolutiva (CSIC-Universitat Pompeu Fabra), Barcelona, Spain.
3 Ecology and Genetics Research Unit, University of Oulu, Oulu, Finland.
4 Institut Botànic de Barcelona (IBB, CSIC-CMCNB), Barcelona, Spain.
5 Centre for Biodiversity Genomics, Guelph, Canada.
6 Department of Integrative Biology, University of Guelph, Guelph, Canada.
* Corresponding author. Email: leonardo.dapporto@unifi.it (L.D.); roger.vila@csic.es (R.V.)
† These authors contributed equally to this work.

20 9 2024
18 9 2024
10 38 eadm859609 11 2023
12 8 2024
Copyright © 2024 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).
2024
The Authors
https://creativecommons.org/licenses/by-nc/4.0/ This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

The interplay between geographic barriers and climatic oscillations over the past 2.6 million years structured genetic variation at the continental scale. The genetic legacy of the Quaternary ice ages (GLQ) hypothesis outlines this phenomenon for Europe, but a comprehensive data-driven assessment is lacking. Using innovative genetic landscape methods, we model the GLQ in the West Palearctic based on 31,653 Cytochrome c oxidase subunit 1 (COI) sequences from 494 butterfly species and three functional traits. Seven distinct bioregions with varying levels of genetic endemicity emerge, revealing a latitudinal gradient in variation that confirms the “southern richness, northern purity” hypothesis. Through shift from case studies to a comparative approach, we objectively identify the main glacial refugia, colonization routes, and barriers to dispersal. Our findings offer a quantitative model of the GLQ across Europe, North Africa, and neighboring Asia, with broader applicability to other taxa and potentially scalable to encompass life on Earth.

A genetic analysis of 494 butterfly species and 31,653 DNA sequences reveals the effects of past climatic cycles on biodiversity.

http://dx.doi.org/10.13039/501100002809 Generalitat de Catalunya 2021-SGR-00420 http://dx.doi.org/10.13039/501100002809 Generalitat de Catalunya Grant 2021-SGR-01334 Academy of Finland 324988, 328895, 352652 MCIN/AEI/ 10.13039/501100011033 PID2019-107078GB-I00 and PID2022-139689NB-I00 MCIN Grant PID2020-117739GA-I00 MCIN / AEI / 10.13039/501100011033 MICIU PID2022-139689NB-I00 Departament de Recerca i Universitats, Generalitat de Catalunya 2021-SGR-00420 Horizon project “Biodiversity Genomics of Europe” N/A Direttiva Biodiversità projects 2019-2020 of the Italian Ministero della Transizione Ecologica N/A NBFC (Italian Ministry of University and Research, PNRR) Project CN00000033930 ”la Caixa” Foundation 100010434, grant LCF/BQ/DR20/11790020 Canada Foundation for Innovation and NSERC N/A ERDF EU N/A Fellowship FPU19/01593 N/A
==== Body
pmcINTRODUCTION

The Quaternary period (2.6 Ma to present) has been characterized by climatic pulses that transformed the composition and distribution of the global biotas (1–3). Seminal papers by Taberlet and co-workers (2) and by Hewitt (1) described the effects of these glacial cycles on the genetic structure of European biota, which have been termed “the genetic legacy of the Quaternary ice ages” (GLQ hereafter). This framework envisages the existence of distinct southern refugia where genetic lineages diverged during long Quaternary glacial cycles. During interglacials, species have likely recolonized northern regions following shared routes. According to the GLQ hypothesis, these phenomena have left deep fingerprints, still visible at the faunistic and genetic levels in current populations. These fingerprints may vary among species characterized by different functional traits, especially those determining their mobility and colonization abilities (4–6). Although on the basis of case studies, each focusing on a few species, the GLQ has influenced phylogeography—the study of the principles governing the geographic distribution of intraspecific genetic variation (IGV) (6–8)—for more than two decades (3, 9–11). However, with the emergence of a dominant paradigm, there is a tendency to favor data that support it while disregarding conflicting results (12). To objectively assess the GLQ, a comparative approach (6) is needed. The ideal dataset should: (i) involve enough species of a diverse taxonomic group to allow a statistically meaningful comparative approach; (ii) include specimens classified according to a well-established taxonomic framework; (iii) use sufficient samples to ensure high spatial resolution; and (iv) cover all Europe and, ideally, neighboring regions of Africa and Asia, which may contribute additional genetic variation.

To put the GLQ paradigm to the test across an entire animal group, we gathered 31,653 mitochondrial DNA (mtDNA) Cytochrome c oxidase subunit 1 (COI) sequences (barcoding marker) for 494 butterfly species from the Western Palearctic region (13). Using this dataset, we aimed to verify the following paradigms of GLQ: (i) Glacial refugia: Identifying recurring distribution patterns among genetic lineages and endemic haplotypes should reveal glacial refugia—both within and beyond the Mediterranean region—and the areas subsequently recolonized from these refugia. (ii) Physical barriers: Mountains and sea straits likely impeded postglacial expansion, mostly in species with low dispersal capacity, resulting in the presence of areas with pronounced phylogeographic breaks. (iii) A latitudinal trend of IGV: Postglacial expansion is expected to produce genetically impoverished northern populations. These patterns could partly be explained by species’ functional traits, and we test for the influence of three functional traits in different aspects of the GLQ.

RESULTS

Regionalization based on butterfly IGV

We divided the study region in areas represented by mainland cells of 250 km by 250 km plus islands (fig. S2). We obtained overall dissimilarity matrices among these areas by calculating genetic distances for all available species pairwise, using three substitution models separately and then by averaging these distances among species using six weighting and averaging methods, thus obtaining 18 dissimilarity matrices. We applied to these dissimilarity matrices a widely used and robust regionalization framework with two different aggregative and divisive clustering methods (14, 15), so that we obtained 36 clustering solutions (Materials and Methods and fig. S1A). The solutions were averaged with a Gower dissimilarity, generating a final dissimilarity matrix to which a 20% of spatial contiguity constraint (16) was added. A final clustering procedure on this matrix produced biogeographic regions (bioregions) characterized by a low IGV within areas of the same bioregion and a high IGV among bioregions.

The solutions identifying an increasing number of clusters of areas (k; fig. S3) produced an expected decrease in the minimum fraction of species showing significant genetic differences among pairs of bioregions in analyses of molecular variance (AMOVAs) (the strength of the weakest geographic division) (table S1 and fig. S4). The majority of species (>50%) showed significant differentiation among pairs of bioregions from k = 2 to k = 8, and two main drops occurred at k = 5 and k = 9. The explained genetic dissimilarity (14) encompassed by pairs of cells from different bioregions showed a rapid rise with values higher than 85% for k > 5. We also measured the total spatial overlap among bioregions using convex hulls, which was lowest for k = 2 and k = 7, indicating a stronger spatial structuring for these solutions (fig. S4). Considering these values, we selected k = 7 as the best compromise; it recovered the following bioregions: (i) Maghreb, (ii) Iberia (including the Balearic islands), (iii) Central-Western Europe (including part of the British Isles, Corsica, Sardinia, and Sicily), (iv) Italian Peninsula and Alps (including most circum-Italian islands), (v) Balkans-Anatolia (including Crete, Cyprus, and smaller Aegean islands), (vi) Fennoscandia (including a part of the British Isles), and (vii) Carpathians (Fig. 1A). The solution for k = 8 (fig. S3) separating the Italian Peninsula from the Alps (17, 18) showed similar support despite the slightly higher spatial overlap. A solution without soft contiguity constraint showed a very similar solution for guide areas (fig. S5).

Fig. 1. Regionalization of the West Palearctic based on butterfly IGV.

(A) Map of the bioregions identified for the k = 7 solution. Large mainland squares and islands marked with circles are those with most data and represent guide areas used to obtain the bioregions; the small mainland squares and islands marked with triangles represent nonguide areas attributed a posteriori to their bioregion. (B) The Ward tree obtained from mean PhiST values of species shared among pairs of bioregions identified by the regionalization procedure. The inset radar graphs represent the mean similarity (1 − PhiST) between each bioregion (the color of the central dot) and the other bioregions (circular sectors) used to construct the tree. (C) Level of intraspecific endemicity measured as number of EMs per square-rooted number of specimens (dots represent individual species, left axis), and estimated marginal means and their SE (whiskers) for the same measure obtained after a GLMM comparing bioregions (right axis). Three groups of bioregions (N. Africa, S. Peninsulas, and C-N Europe, indicated in boxes) were significantly different to each other and not differing within groups.

The mean PhiST values from AMOVAs among species and the resulting a posteriori Ward tree revealed that the Maghreb harbors the most differentiated butterfly fauna, followed by Iberia (Fig. 1B). The following division separated the group represented by Italian Peninsula-Alps and Central-Western Europe from the group including Balkans-Anatolia, Fennoscandia, and Carpathians.

Patterns of regionalization were further explored for the number of endemic mutations (EMs) per species. EMs for each species were obtained on haplotype networks as the number of mutations characterizing each exclusive haplotype from the closest haplotype (fig. S6). EMs have been identified, standardized by square-rooting the number of sequenced individuals, for each species around each cell of 25 km by 25 km and for each bioregion.

A generalized linear mixed model (GLMM) was applied to highlight differences in EMs for each species in each bioregion (response variable) in comparison with bioregion membership and three species traits (phenology, wingspan, and number of used host plant genera); species membership nested with genus was included as a random variable. Significant differences in EMs among bioregions were retrieved (GLMM, n = 1517 species × bioregions, χ2 = 167.841, df = 6, P < 0.001). Among the three functional traits, phenology had a significant effect, with species showing longer flight periods expressing fewer EMs (χ2 = 11.588, df = 1, P < 0.001), while wingspan and number of host plant genera did not show a significant effect (wingspan, χ2 = 0.308, df = 1, P = 0.579; host plant genera, χ2 = 2.677, df = 1, P = 0.102). This result agrees with previous studies where a shorter flight period is often linked to higher genetic differentiation in butterflies, because adults mediate dispersal and gene flow (5, 19). Pairwise comparisons (Fig. 1C and table S2) showed that the Maghreb has significantly higher haplotype endemicity than all other bioregions. In addition, Iberia, Italian Peninsula-Alps, and Balkans-Anatolia were all significantly richer in EMs than Central-Western Europe, Carpathians, and Fennoscandia. In the GLMM, most variance was explained by differences among species (full model conditional R2 = 0.446), while differences among bioregions and species traits had a smaller effect (marginal R2 = 0.069). This confirms the existence of a strong idiosyncrasy among taxa in phylogeographic patterns (20, 21).

In the analysis of EMs for individual cells, we considered endemic haplotypes those occurring exclusively in an area of 325-km radius around a cell. EM values for each species were further standardized by forcing the mean to zero, so that cells with lower-than-mean values scored a negative contribution for a given species. Summing up the standardized EMs for each cell allowed to identify cells with positive values (several species showing concordant positive values). We lastly visualized the cells in the highest 10% quantile of this value (main hotspots) (Fig. 2A). The cells within the top 25% quantile (including also minor hotspots) are visualized in fig. S7. There was a positive correlation (Pearson R = 0.295, P < 0.001; Fig. 2B) between the number of species sequenced in a 125-km area around a given cell (regional richness and sequencing effort) and their cumulative standardized EM, but the cells falling in the top 10% quantile occurred all along the species richness variance and were mainly located around the Mediterranean region.

Fig. 2. The genetic endemicity hotspots.

(A) Map showing the regions within the highest 10% quantile for cumulative EMs. (B) Relationship between the EMs of 25 km by 25 km cells and the number of species sequenced in the 125-km-radius buffers around them. Colors differentiate three types of genetic endemicity hotspots, based on the number of species sequenced.

Main hotspots occurred in mountain-hill areas and islands, and we identified large species-rich hotspots with >140 analyzed species in Iberia and central Europe (Alps, Pyrenees, and Cantabrian chains; red in Fig. 2), intermediate species-rich hotspots with species ranging between 70 and 140 (Baetic chain in Iberia, Southern Italian Peninsula, Sicily, Peloponnesus in the Balkans, and Atlas chain in the Maghreb; blue in Fig. 2) and species-poor hotspots with <70 species (Corsica, Crete, Cyprus, and some areas of the Middle East, the latter highly underestimated for actual species richness because of undersampling; green in Fig. 2). The 25% top quantile expanded the main hotspots with minor hotspots with a special reference to Sardinia, Central Apennines, Northern Balkans, and Carpathians and to Northern and Central Europe with special reference to marginal areas in British Isles and Scandinavia (fig. S7) (22).

Identification of barriers to dispersal

Genetic landscapes were performed (fig. S1B) to test for the effect of key geographic barriers in producing phylogenetic breaks. We created a Delaunay triangulation for sequences of each species (with at least six specimens, three areas sampled, and one mutation differentiating the sequences) and computed minimum cost paths connecting them. For each path, we recorded its length, the maximum altitude and the length of sea straits encountered; for each species, we also included the three aforementioned traits. Then, we carried out a generalized additive mixed model (GAMM) to understand how p-distances among pairs of individuals are influenced by these predictors, including possible curvilinear trends, and using species as a random effect. As expected under a model of isolation by distance, a GAMM (n = 63,445 midpoints, 328 species) revealed a highly significant relationship between the genetic p-distances among pairs of specimens and the log-transformed length of the overground paths separating them (Fig. 3B and table S3). The length of sea straits encountered along the paths also showed a significant curvilinear trend (Fig. 3D). Last, the paths crossing higher maximum altitudes showed stronger phylogeographic breaks, with a significant curvilinear trend (Fig. 3E). Latitude also showed a significant curvilinear effect, with the highest values corresponding to southern areas at 30°N to 45°N (Mediterranean) and then declining to Central and Northern Europe, corresponding to the occurrence of cold-period ice sheets (Fig. 3C). As found for the number of EMs, only phenology had a significant effect among the three functional traits, with species showing longer flight periods expressing lower genetic p-distances across pairs of sites (table S2).

Fig. 3. General pattern of phylogeographic breaks for West Palearctic butterflies.

(A) Map illustrating the mean effective migration surfaces across species interpolated by inverse distance weighting. Blue and red represent surfaces of high (corridors) and low (barriers) effective migration, respectively. (B to E) Effects of different variables on the genetic p-distances among individuals in a GAMM analysis: logarithm of path length (B), latitude of the path midpoint (C), logarithm of length of sea straits encountered along the path (D), and maximum altitude encountered along the path (E).

The spatial distribution of phylogeographic breaks and corridors was mapped using Fast and Flexible Estimation of Effective Migration Surfaces (FEEMS) (23). Effective migration surfaces were obtained for 320 species with minimum requirements as above for genetic landscapes (FEEMS analysis failed in eight species) over a triangle grid with a spacing of ≈100 km between centroids. We then exported the effective migration surface layers and averaged the results among species (24). We interpolated the averaged values by inverse distance weighting on QGIS 2.18. The resulting map (Fig. 3A) shows that main phylogenetic breaks occur along the Mediterranean Sea. Mountain areas appear as weaker dispersal barriers compared to sea, limited to the main chains (Pyrenees and Alps) delimiting two main refugia hypothesized by the GLQ: Iberian and Italian peninsulas. In the third refugium, the Balkan Peninsula, a relatively weak barrier to dispersal is found in Northern Greece. River valleys are not recovered as barriers to dispersal, with the remarkable exception of the Ebro valley in Iberia, which emerges as the main barrier within this peninsula. In central Europe, the only phylogeographic break is found in the Alps, extending northward into Germany (1, 2). Northern Europe generally presents much higher effective migration values, apart from a barrier highlighted in northern Fennoscandia and between Britain and Hebrides (22). Other areas of low effective migration values along the margins of the study area may be due to a FEEMS methodological artifact and to a reduced number of specimens and species sampled for those areas.

Latitudinal trends in IGV

Latitudinal trends in different aspects of biodiversity have been found or proposed for a variety of organisms (11, 20, 25–28). The “southern richness-northern purity” hypothesis states that, in temperate areas of the Northern Hemisphere that were partly covered by ice sheets during the glacial maxima, the IGV is inversely correlated with latitude. This phenomenon is postulated to be produced by gene surfing during postglacial recolonization: repeated and spatially oriented founder effects that resulted in an increase of particular variants and, thus, a gradual impoverishment of genetic diversity (1, 29). We applied a widely used band-wise method (26, 30) (fig. S1C) by dividing the study area in 42 bands, each 2° wide, with 1° overlap the nearby ones. For each species in each latitudinal band, we computed the mean p-distances across all pairs of sequences (25). To obtain insights about the mechanisms that generated the observed latitudinal gradient of IGV, we ascertained whether each pair of sequences within each band belonged to the same bioregion (intraregion IGV) or to different bioregions (inter-region IGV). In particular, we tested three main hypotheses regarding whether postglacial recolonization resulted in secondary sympatry and whether glacial refugia contributed equally or unevenly to the recolonization process (illustrated in Fig. 4, A to C).

Fig. 4. Three main hypotheses for a latitudinal trend in IGV.

(A to C) Representation of a system characterized by glacial refugia (triangles at the bottom) and postglacial expansion toward the smaller polar area (top). Genetic variability is indicated by different colors of dots (individuals) with a hypothetical haplotype network. The inset graphs (top right) summarize the expected trends of inter-region (red lines) and intraregion (blue lines) IGV versus latitude. (A) Distinct and genetically diverse glacial refugia similarly contribute to colonization of northern areas by gene surfing without secondary sympatry in recolonized areas; this should produce northern decrease in intraregion IGV and latitudinal maintenance of inter-regions IGV. (B) As in (A), but lineages come into secondary sympatry in recolonized smaller northern areas; this should lead to a northern increase in intraregion IGV and decrease in inter-regions IGV. (C) As in (B), but refugia contribute differently to postglacial recolonization; this should determine decreasing intraregion and inter-regions IGV with latitude. (D) The significant effects of latitude, type of IGV, and their interaction on IGV obtained after a GLMM. The continuous sections of regression lines represent the latitude interval for which data were available, and discontinuous extremes represent extrapolations.

A GLMM for 310 species with minimum requirements (n = 6206 species × band, species nested in genus and band as random factors) revealed that IGV increased as the sampling area expanded (distances among specimens) but that it was not affected by the number of sequences available in each band (table S4). We found an overall decrease of IGV in northern areas and a significantly higher inter-region IGV than intraregion IGV; moreover, a significant interaction between latitude and IGV type (inter- or intraregion) (table S4) revealed that the latitudinal trend was mostly due to an increased inter-region IGV in southern areas (Fig. 4D). Independent GLMMs for intra (n = 2361 species × band) and inter-region (n = 3845 species × band) IGV showed that latitude had a significant negative effect in both cases (table S5). This evidence strongly supports the third hypothesis for the generation of a latitudinal trend of IGV (Fig. 4C): an imbalanced contribution of different refugia to the northward postglacial colonization, producing gene surfing and founder effects and with potential secondary sympatry after contact.

As found in previous analyses, species characterized by short flight periods (fig. S8A) had elevated IGV. The number of host plants used also showed a strongly significant interaction with latitude, with specialized species showing a steeper decrease in IGV with latitude than generalists (fig. S8B). Most variance in the GLMM was explained by differences among species (full model conditional R2 = 0.518), while geographic determinants (distance and latitude) and species traits had a smaller effect (marginal R2 = 0.176).

An objective model of GLQ

We quantitatively inferred colonization routes and barriers to dispersal and produced an objective model of GLQ for West Palearctic butterflies (Fig. 5A). Colonization routes were established on the basis of PhiST values among pairs of bioregions lower than 0.3 (the most similar bioregions in Fig. 1B and table S1), which were all located on the European mainland. The direction of these routes was based on significant differences in the number of EMs between the bioregions and on the occurrence of endemicity hotspots in source areas (Fig. 1B): higher endemicity for the source, lower endemicity for the region colonized. All the connections that started in the bioregions representing the three south European peninsulas were oriented northward. Thus, these were interpreted as reflecting postglacial northward expansion routes. The rest of connections could not be oriented and are considered bidirectional, on the basis of the lack of significant EM trends. These include the only longitudinal connection, between Central-Western Europe and the Carpathians, and two latitudinal connections between these areas and Fennoscandia. Worth noting, none of these three regions showed main endemicity hotspots.

Fig. 5. The model of the GLQ for West Palearctic butterflies.

(A) Barriers to dispersal (red) are the lowest (10%) quartile cells for average effective migration surfaces. Dispersal routes (black arrows) connect the bioregions showing PhiST values lower than 0.3, oriented from significantly higher EMs to lower. Thick solid arrows connect bioregions showing mean PhiST values <0.20, thin solid arrows represent values between 0.20 and 0.25, and dashed arrows represent values between 0.25 and 0.30. (B and C) Main routes of postglacial recolonization and suture zones proposed by Taberlet and co-workers (2) (B) and Hewitt (1) (C).

Main barriers to dispersal were revealed by mapping the cells in the lowest 10% quantile for effective migration surfaces (Fig. 3A). These phylogenetic breaks show a remarkable difference to suture zones inferred by the GLQ paradigm: They are located both over sea and mainland, while the GLQ paradigm highlighted only mainland regions where lineages converge in narrow areas of secondary sympatry (Fig. 5, B and C). The use of migration surfaces and the inclusion of North Africa made it possible to document strong breaks along Mediterranean Sea straits (Fig. 5A), which confirms that the sea represents a main barrier for butterfly dispersal (31–33).

DISCUSSION

The GLQ paradigm has fascinated biologists for its clear predictions about the distribution of genetic variation. While predominant patterns can be extracted, we acknowledge that phylogeography should be still considered a “tale of tales” (21) because of the strong species idiosyncrasy. The genetic fingerprints of past events may be differentially preserved: The most mobile and generalist species—those showing a long adult flight period and using many host plants—had a lower incidence of endemic lineages, as well as a lower tendency to establish spatial differentiation and to show latitudinal trends of genetic diversity (5, 34).

Despite the complexity and stochasticity of processes involved in the genesis of IGV, macrogenetics is an effective approach for inferring general phylogeographic theories, taking advantage of the increasing availability of DNA sequences (5, 11, 18, 20, 30). Scaling up classical analyses and using state-of-the art methods on a large genetic dataset comprising virtually all West Palearctic butterfly species, we could identify main and common processes determining the distribution of IGV in this region. Current efforts dedicated to mtDNA sequencing, mostly in the form of DNA-barcoding libraries, will enable direct comparisons among different taxa and regions. Because butterflies are highly mobile organisms (35), it is possible that the IGV patterns generated by GLQ could be even stronger in taxa with lower dispersal capacity. The use of multiple nuclear markers will test the generality of the patterns found, including patterns produced by male dispersal, and expand the focus to cover adaptation and hybridization among lineages (6, 36).

Regionalization based on butterfly IGV

The quest for a shared distribution of genetic lineages that proves the existence of a common set of historical vicariant events and postglacial recolonization routes is the grail of phylogeography (1, 2, 5, 8, 37). In general, the importance of the Mediterranean as both a glacial refugium and a center of diversification is illustrated by the fact that bioregions in this area showed the highest distinctiveness and/or frequency of EMs. This pattern also indicates that, during glacial periods, the number of species that survived in extra-Mediterranean refugia (38), likely located in the central European bioregions, was comparatively low.

The Mediterranean bioregions identified in this study have a clear correspondence to the southern European peninsulas identified as refugia by the GLQ (1, 2). To these, our analysis adds the Maghreb, a region generally suitable for the survival of butterfly populations during glacial maxima (39). Central-Western Europe and Carpathians (and Alps identified as a distinct bioregion for k = 8) comprise potential extra-Mediterranean refugia (38, 40) but also represent recent ecological centers of endemism recolonized by several species during postglacial population expansions (17, 41). Last, Fennoscandia, most of Britain, and North-East Europe, were largely covered by ice sheets during glacial periods and, hence, devoid of butterfly fauna. Current populations in this entire region, except perhaps a handful of arctic species, represent postglacial colonizers. Sicily, Sardinia, and Corsica were grouped with Central-Western Europe, which is not unexpected given the recurrent evidence of similarity among these areas in butterfly genetics (5, 42) and the strong contrast between these islands and the Italian mainland (32).

We identified the Maghreb as the most differentiated bioregion. There is no influence on this result of the fauna typical of this area and belonging to the Afrotropical region, because in the analysis of IGV and for EMs, only genetic differences among species shared among bioregions were considered. The long-term separation between Europe and Africa has likely generated the observed pattern, despite the short distance (14 km) between the two continents at the Strait of Gibraltar. The Maghreb is possibly the area currently facing the highest conservation threat: Extinctions of butterfly populations are rampant because of anthropogenic pressure (e.g., habitat fragmentation, overgrazing, and climate change) (43) undoubtedly making it a genetic erosion hotspot. These phenomena are also important throughout the South European peninsulas—here confirmed as three distinct bioregions (2, 3, 44)—mostly in isolated mountain areas, here identified as main endemicity hotspots, where small relict populations are threatened by current climatic changes (45, 46). Among the European bioregions, Iberia emerged as the most divergent unit. This is confirmed by the observation that three distinct 90%-quantile genetic hotspots were retrieved for this bioregion. The peripheral position of Iberia and the presence of a series of isolated mountain chains are known to have produced a complex system of refugia within a refugium (47), where Iberian lineages have diversified into different endemic subgroups (48, 49). In addition, the Ebro valley was retrieved as a main barrier to dispersal, which further helped isolating and subdividing this refugium. Moreover, the westernmost marginal location of Iberia in the Palearctic has also probably preserved ancestral and endemic populations from selective sweeps (22) by eastern lineages.

Latitudinal trends in IGV

Our results support the gradual accumulation of genetic variants in areas where species persisted over long intervals. The higher IGV shown by inter-region comparisons and the steeper relationship with latitude confirm the primacy of divergence among allopatric populations in different refugia as the main driver for genetic differentiation in the West Palearctic. By comparison, the accumulation of mutations in the same refugium was a weaker mechanism. While the decrease of intraregion IGV with latitude agrees with well-established models of gene surfing and founder effect during recolonization of formerly glaciated areas (29, 30), the strong trend in inter-region IGV has been less investigated. A first explanation derives from our regionalization results demonstrating that central European bioregions are genetically less diversified to each other than southern ones. Moreover, despite the common observation that suture zones among lineages tend to be narrow (29), dispersal to northern regions forces the lineages to converge in relatively smaller circumpolar areas, facilitating secondary sympatry. Secondary sympatry in northern areas can increase intraregion IGV, thus counterbalancing genetic impoverishment caused by gene surfing. This is evidenced by a flatter intraregion IGV trend with latitude and by the fact that intraregion IGV was extrapolated to exceed inter-region IGV at latitudes above 55°N (Fig. 4D). The occurrence of subarctic melting pots of genetic diversity is a well-known phenomenon (9, 50) also supported in our study by the presence of a phylogeographic break in northern Fennoscandia (Fig. 3A).

An objective model of GLQ

The postglacial colonization routes for butterflies are unlike any previously proposed paradigm in the GLQ hypothesis: foremost the out-of-Italy route, then out-of-the-Balkans and, lastly, out-of-Iberia. Hence, the Italian peninsula and the Alps appeared as the main sources of northward expanding populations for Western and central Europe, and the Carpathians bioregion was equally sourced by Italy and the Balkans. This pattern aligns with that described for the butterfly Melanargia galathea and related taxa, for which the term “fourth paradigm” was coined (51) as an additional exemplary pattern to the three described by Hewitt (1). We can here generalize this paradigm as the main phylogeographic pattern of West Palearctic butterflies.

The absence of physical barriers has been proposed as one possible reason for the large contribution of lineages made by the Balkan refugium to central European biotas (44). This conclusion is only partially supported by our results, because the Balkans bioregion is quite different from the Carpathians and Fennoscandia and because a phylogeographic break separating the Peloponnesus and Anatolia genetic hotspots from rest of the Balkans was obtained. Likely, several lineages that diverged in the southern Balkans during the glacial periods did not disperse northwards across the Danube river basin and are now altitudinally restricted to mountain areas that function as an interglacial refugium. In turn, the Carpathians, a minor endemicity hotspot revealed by the 75% quantile (fig. S7), exchanged genetic lineages with Northern Europe through the wide migration surface represented by Northern-Eastern Europe (Fig. 3A). The substantial gene flow with multiple bioregions and the interglacial immigration of lineages from southern peninsulas explain why Central-Western Europe, Carpathians, and Fennoscandia are less distinct, as evidenced by their lower PhiST values and few EMs. Nevertheless, the differential input from southern refugia and the putative existence of extra-Mediterranean refugia bordering the Alps and Carpathians (38, 52), at least for a reduced number of species, could have contributed to the distinctiveness of these bioregions (17, 53, 54).

Although the Western Palearctic is bounded to the south by the Sahara desert and by the Atlantic Ocean to the North-West, it lacks clear boundaries to the East. The repeated colonization of Northern Europe by populations from Siberia during the interglacial, a recurrent pattern in butterfly phylogeography (49, 55, 56) not tested here, was likely the main factor leading this formerly glaciated area to gain status as a distinct bioregion compared to other European areas.

Methodological considerations and limitations

Although the mitochondrial COI gene is the primary marker used in phylogeography, mitonuclear discordance is widely recognized (6, 57–61). First, mtDNA is maternally inherited and reflects exclusively the biogeographic and population history of females (6, 58). In addition, a relatively high incidence of mtDNA introgression across species (57) may lead to conflicts with species trees. Last, selective sweeps are relatively common in mtDNA, sometimes as a consequence of infections by endosymbionts like Wolbachia (58). Because of these phenomena, there is contrasting evidence about the correlation between the patterns exposed by nuclear DNA and mtDNA (6, 59). Nevertheless, macrogenetic studies based on mtDNA revealed large-scale patterns in IGV (11, 26, 60, 61), and, up to now, there is no alternative marker to COI with sufficient data available to carry out comparative analyses at broad taxonomic and spatial scope (11, 25, 59, 60). We hope that the availability of sequencing data obtained from next-generation sequencing on a large number of species will allow, in the near future, for the provision of more refined data to replicate similar studies based on multiple nuclear markers (6). Last, taxonomic decisions can influence the evaluation of IGV and, while butterflies are relatively well studied, some cryptic species are still being described requiring continuous updates in the checklists (13, 62). Nevertheless, unfound cryptic species are expected to comprise a low proportion of this dataset and should have a small impact on the regionalization results.

General implications

In line with the GLQ paradigm, we found that IGV was mostly generated by the lineages that diverged in allopatry in glacial refugia, also representing hotspots of genetic endemicity. Current and past geography together with climate set the barriers to gene flow. Our analyses identified sea straits in the Mediterranean Sea as the main obstacles to dispersal during the Quaternary. Although phylogeographic breaks were stronger along paths crossing high-altitude areas, we found that mountain chains represent not only phylogeographic breaks for butterflies but also centers of genetic endemicity. It has to be noted that European mountains represent the richest regions for butterfly alpha diversity, because many cold-tolerant species (and populations) widely distributed during glacial periods have recently converged on these mountains (41, 63). In addition to being barriers during the glacial periods and currently for a subset of warm-adapted taxa, they function as interglacial refugia for many species, therefore representing reservoirs of genetic diversity (hotspots) from where propagules may disperse (17, 18).

Last, the study of latitudinal trends revealed that three neutral processes occurred simultaneously during the northward expansions of populations in interglacials: (i) gene surfing decreasing intraregion IGV from endemicity hotspots in the south to northern areas; (ii) expansions from only a subset of refugia decreasing inter-region IGV in the north; and (iii) secondary sympatry in the north increasing intraregion IGV and decreasing genetic differentiation among region with latitude. These three processes jointly produced the observed pattern of latitudinal trend in IGV but also strongly contributed to the intraspecific regionalization and to the location of phylogeographic breaks. Thus, the spatial pattern of IGV in Europe can be summarized as longitudinally differentiated and latitudinally impoverished.

Knowledge of processes that transpired during past climate changes provides insights into the potential impacts of the ongoing climatic crisis on biodiversity. In particular, our results suggest that, while the expansion of taxa northwards in response to global warming (64, 65) will probably favor the survival of species, populations occupying these regions will likely be genetically impoverished. Because most genetic diversity involves the differentiation of the Maghreb from Eurasia, as well as differentiation among southern European refugia, specific conservation efforts should focus on preserving heterogeneous landscapes in southern regions, with a special reference to mountain chains, to minimize the loss of intraspecific genetic lineages. In central and Northern Europe, maintaining multiple biological corridors can foster northward dispersal while minimizing genetic filtering.

MATERIALS AND METHODS

Sampling and COI data filtering

We obtained 31,653 DNA barcode sequences of >400 bp for 494 species of butterflies from the Western Palearctic region from the recently published iodatabase dataset (13). All sequences derived from specimens with information on the site of their collection in decimal degrees coordinates [World-Geodetic-System 1984 (WGS84)] and with identifications based on the most recent taxonomic checklist for Western Palearctic butterflies (13). Sequences were selected within a rectangle delimited by the following Mercator coordinates: −13.25° to 40° longitude, 27° to 73° latitude. An important limitation in the use of massive mitochondrial COI data in IGV assessments stems from the fact that, in some published studies, only single or newly found haplotypes are deposited in data repositories instead of the full datasets (59) thus providing a false signal of genetic diversity. The iodatabase dataset is not affected by this limitation (13).

Biogeographic regionalization

Regionalization was generated with a modified version of the algorithm proposed by Dapporto and co-workers (5) (fig. S1A), which is an application to IGV of the typical framework adopted for regionalization of faunistic data (14). To obtain equal-area mainland cells, WGS84 coordinates were transformed into Lambert azimuthal equal-area projection. COI sequences were assigned to two kinds of areas: islands and mainland grid cells of 250 km by 250 km and largely corresponding to the cell size used in macrogenetic studies (25). Because of their large size, Britain and Ireland were distinguished from the European continent and also divided into grid cells. Some cells that covered a very small area of land and islands very close to each other (less than 3 km) were aggregated with the nearest area at the same latitude (see R script in Dryad).

Three distance matrices of pairwise interarea mean genetic distances using p-distances and Felsenstein’s (1981) and Tamura and Nei (1993) models were subsequently assembled for each species using the COI sequences in each area using the dist.dna function of the “ape” R package (66) (three series of models). The average genetic differentiation among any pair of areas (Di,j) has been calculated using the following formulaDi,j=∑s=1ri,jdi,j,s×wsn∑s=1ri,jwsnn(1)

where D is a mean genetic distance among specimens from areas i and j calculated for each species (s) using the three different substitution models. The n exponent was set to 1 and 2 to provide a linear or quadratic importance to species with higher differentiation (two series of models). The contribution of species to the mean genetic distance is weighted (ws) with three different approaches: (i) no weight with ws equal to 1 for all species (series 1); (ii) the distribution with ws equal to the number of areas they occupy to provide a higher contribution to species showing wider distribution at the continental level (series 2); and (iii) the spatial structure in genetic differentiation with ws equal to the species GST, to increase the effect of species showing high interarea diversity as opposed to sympatric genetic diversity (series 3). GST was defined as (HT − HS)/HT (7), where HT is the average of p-distances for all the barcode records belonging to each pair of areas, and HS is the average of p-distances within areas. This index ranges between 0 and 1, with 0 implying a species showing panmixia and 1 indicating a complete spatial structuring among areas. GST depends on grain (cell size) and on the arbitrary subdivision of the mainland cells. For this reason, we computed, for each species, the GST among cells of 100 km by 100 km, 200 km by 200 km, 300 km by 300 km, 400 km by 400 km, and we replicated the four analyses by placing the center of each cell on the vertex among cells in the previous analysis (13). The average GST for each species among the 10 values obtained has been used as the generalized tendency for that species to show spatial structuring.

We then selected for each combination of models (18; 3 distances × 2 exponent × 3 weights) the pairwise average genetic differentiation (Di,j) among areas that contained at least 14 species each and sharing at least half this minimum requirement (seven species). These 18 among-areas distance matrices based on all the examined species had a large number of missing pairwise values due to incomplete sampling in some areas. Following a recent paper (5), we ordered the areas on the basis of the decreasing number of empty cells in the dissimilarity matrices, and we removed them in this order unless a set of areas showing a complete dissimilarity matrix was obtained (by use of the guide areas).

The complete distance matrix was subjected to regionalization analyses using the recluster.region function of the recluster R package with two different agglomerative and divisive clustering algorithms: (i) ward.D2 as implemented in the hclust function of the stat R package and (ii) the Divisive ANAlysis Clustering (DIANA) as implemented in the cluster R package (two series of models) obtaining clustering solutions ranging from k = 2 to k = 15 bioregions. The total number of models obtained was 36 (18 distance matrices × 2 clustering algorithms) for each k value. A consensus, for each k value, was obtained by calculating a Gower dissimilarity among the 36 vectors defining bioregion membership of areas in different models. This distance matrix indicated the number of times each pair of areas did not belong to the same cluster. To this matrix, we applied a soft contiguity constraint as done in the hclustgeo function of the “Clustgeo” R package (16). This algorithm combined the distance matrix obtained by Gower dissimilarity with the matrix of geographic distances among area centers to attribute undecided areas to spatially close clusters. To allow genetic differentiation to be dominant over spatial contiguity, we set α to 0.2 producing a relative importance of 80% for the genetic dissimilarity and 20% for geographic distances. The resulting dissimilarity matrices, one for each k value, were subjected to a final recluster.region clustering (ward.D2 method) asking for the formation of number of groups corresponding to k values.

After creating a clustering solution for the guide areas showing a complete average genetic dissimilarity matrix, the areas excluded because of the presence of NA values in the dissimilarity matrix but including at least 14 species (nonguide areas) have been re-entered in the analysis iteratively as follows. An average dissimilarity matrix among the 18 original dissimilarity matrices (containing all the areas with more than 14 species) is obtained and weighted by geographic distance as done before. Then, nonguide areas are inversely ordered by number of missing values in the resulting dissimilarity matrix, and we calculated for each area the mean of the dissimilarities of that area with respect to those previously attributed included in bioregions (both guide and nonguide areas). We then attributed the area under examination to the bioregion showing the lowest mean distance. We replicated the procedure unless all the nonguide areas are attributed.

To evaluate distinctiveness among bioregions, we carried out pairwise AMOVAs between pairs of bioregions for each shared species with at least five specimens in each of the two bioregions (using the function pairPhiST of the “haplotypes” R package).

The strength of the regionalization solutions from k = 2 to k = 15 were compared using three parameters: (i) the spatial overlaps among areas, (ii) the explained dissimilarity [the fraction of inter-region overall genetic distances divided by the total distances; a strong regionalization solution should explain around 90% of the dissimilarity variance (14)], and (iii) the minimum fraction of species showing a significant geographic structure in pairwise AMOVAs among bioregions representing the strength of the least supported pair of bioregions. The overlap among bioregions was obtained by using the biodecrypt function of the recluster R package. By using this algorithm, we computed convex hulls joining centers of areas belonging to the same cluster, and we obtained the overlaps among pairs of areas (in square kilometers). For each k value, we obtained the total overlap and we considered as more spatially structured the solution(s) showing the lowest overlap. Explained dissimilarity was obtained using the recluster.explain R function in the “recluster” package (67).

Intraspecific genetic endemicity

We lastly quantified intraspecific endemicity of bioregions as the number and extent of genetic divergence of haplotypes occurring in single area units. As area units, we used two different resolutions: the bioregions identified using regionalization and circular areas of 125-km radius around cells from a grid of 25 km by 25 km. Endemic haplotypes for bioregions were those found in a single bioregion and nowhere else in a dataset including all sequences for the West Palearctic, Asia (68), and North America (69). Similarly, endemic haplotypes among the specimens sequenced in the 125-km radius units were those occurring exclusively in an area around the cell of 325-km radius. We used species occurring in at least two bioregions or occurring both within and outside the 325-km radius to exclude mixing intraspecific endemicity and divergence with interspecific endemicity (67, 68).

Once the haplotypes exclusive to each cell and bioregion were identified, the degree of divergence among the endemic haplotypes in a given bioregion for each species was scored as minimum number of mutations from the nearest haplotype(s), EMs. EMs were calculated on haplotype networks using the minimum p-distances from the nearest haplotype (fig. S5). The sum of the EMs was obtained for each species in each bioregion or cell. This procedure was implemented in the function “haplotype.endemicity,” which requires as input a haplotype network obtained with the pegas R package, the original distance matrix, and the list of endemic haplotypes. The number of detected haplotypes does not increase linearly with the number of sampled specimens, being better represented by accumulation curves (30, 69). To obtain a standardized frequency of EMs per species in each area, we divided the sum of the EMs obtained for each species by the square-rooted number of specimens analyzed for that species in that area, thus obtaining the index EMs per square-rooted specimens. Species traits can affect the general tendency of species to show high endemicity, notably because a low dispersal capability can reduce gene flow, thus allowing for genetic differentiation into endemic lineages. We selected three functional traits largely used to describe dispersal and colonization capabilities of butterfly species (70): (i) wing index (a composite measure for size); (ii) phenology (obtained as a first principal components analysis component among length of flight period, number of generations, and first and last month of flight); and (iii) the number of host plant genera used by the larvae. The traits have been updated considering recent shifts in taxonomy, and an updated table is presented in the GLQ package. To highlight areas with excess of EMs at the continental scale, we proceeded using the data obtained for the 25 km by 25 km grid cells as follows. For each species, we scaled and centered the EMs per square root specimens among grid cells to obtain a mean value equal to zero and SD equal to 1. This allowed us to (i) identify cells with EMs lower than the mean value as negative and those higher than the mean as positive and (ii) standardize species among them for their mutation rate and time since the last sweep, thus providing similar absolute values.

Then, values for the different species in each cell have been summed. As a result, if a cell has an excess of species below the mean (likely in area of recent colonization), then they will score negative values as opposed to positive values in refugial areas from where colonization originated.

Highest (or lowest) values in around a cell occur if two conditions occur: (i) Species are in most cases positive (or negative) and (ii) the area around the cell is rich in species. In this optic, the analysis produces hotspots of genetic endemicity.

Phylogeographic breaks

To test for the effect of isolation by distance, of main geographic barriers (mountain chains and sea channels), and of species functional traits (size, degree of generalism, and phenology) on mtDNA IGV, we applied genetic landscapes using a modified published algorithm (71). The procedure is shown in a schematic flowchart in fig. S1B, while the recluster.landscape function used to carry out the analysis is provided in the GLQ R package. A p-distance matrix was computed for each species with at least five COI sequences, and the geolocation of each specimen was used to build a Delaunay triangulation [deldir function of the “deldir” R package (72)]. A least-cost path for the records connected in each triangulation was then calculated under the assumption that dispersing over ground is three times less costly than crossing a sea barrier. This was done by obtaining ground and sea cells in the study area from a climatic layer for mean annual temperature (www.worldclim.org) with 2.5-min spatial resolution. The ground and sea cell layer was transformed into a permeability layer by setting ground cells to a value of 3 and sea cells to 1. A transition layer was then obtained [transition and geoCorrection functions of the “gdistance” R package (73)], and the least-cost path between each pair of records was computed on the transition layer (shortestPatch function of the “gdistance” R package).

A depth-elevation layer (2.5-min spatial resolution) was downloaded from the General Bathymetric Chart of the Oceans GEBCO (www.gebco.net) and clipped to the study area. For each least cost path connecting two specimens in the Delaunay triangulation, we obtained the following: (i) midpoint over the path in decimal latitude and longitude, (ii) length of the path, (iii) length of the sea part of the path, (iv) maximum elevation, and (v) p-distance between COI sequences of connected specimens.

To map dispersal surfaces and phylogeographic breaks, we applied FEEMS (23). We used the same species selected for the genetic landscapes and organized the sequenced specimens over a grid a grid spacing ≈100 km between centroids. We selected the smoothing regularization parameter (λ) for each species separately by selecting the lowest values in a range between 1 × 10−6 and 100, applying a leave-one-out cross-validation procedure. Thereafter, we integrated all results in a single layer by averaging among species the effective migration values (standardized forcing mean equal to zero for each species) in each segment connecting grid centroids (24). QGIS 2.18 was used to interpolate the average layer over the study area by inverse distance weighting with default run settings.

Latitudinal trend of IGV

The latitudinal trend of IGV was assessed using a band-wise method following Miraldo et al. (26) and Dincă et al. (30) (fig. S1C). By using the function “recluster.latitudinal.bands” included in the GLQ package, we divided the study area (28°N to 70°N) into 42 bands, each 2° wide. As the band centers were offset 1° to each other, each band overlapped by 1° with the nearby ones, an approach that smoothed the resulting trend (30). For each latitudinal band, only species represented by at least five specimens were included in the analysis. Genetic diversity for each species was calculated following the approach used in most macrogenetic studies (25): by computing the average number of nucleotide differences per site across all pairwise sequence comparisons for that species, which equates to the mean of a p-distance matrix calculated by using pairwise deletion.

Increased IGV at lower latitudes can result from long species persistence and warmer climate that could produce high genetic diversity within southern regions followed by (re)colonization of northern areas during the interglacial by some lineages only. Another not mutually exclusive mechanism can involve allopatric genetic differentiation among bioregions produced by long persistence in isolated refugia followed by re(colonization) of northern areas by a subset of populations from different refugia (29, 44, 74). To separate and directly compare the role of these two processes, we ascertained for each band whether each pair of sequences belonged to the same bioregion previously identified (intraregion IGV) or to two different bioregions (inter-region IGV). To control for the effects of distance among sequenced specimens that should be greater at low latitudes due to the spherical shape of earth (75), we also evaluated the distances among sequences by using the formulaDist=∑i=1ngeoDistseqi,Bar2n(2)

where geoDist(seqi, Bar) is the geodetic distance between the ith sequenced specimen and the geographic barycenter (average latitude and longitude) among the specimens in the same band, and where n represents the number of specimens in the same band. For intraregion IGV, barycenters and distances were computed among specimens belonging to the same bioregion, and the results from different bioregions of the same band were then averaged. When a bioregion only held a single sequence, the analysis generated a missing value (NA) so it was excluded from the intraregion IGV distances.

Tracing the colonization routes

To trace objective colonization routes and barriers to dispersal, we combined multiple pieces of information: (i) the similarity among bioregions as measured by PhiST; (ii) significant differences between bioregions in the number of EMs and presence of genetic endemicity hotspots; and (iii) areas of low and high effective migration surfaces. In particular, connections among bioregions have been identified as arrows connecting the regions with PhiST values <0.3. We directed the arrows toward areas with a significantly lower number of EMs. Barriers to dispersal were drawn on the basis of the cells with the lowest 10% of effective migration from the average layer obtained from the effective migration surfaces analyses.

Statistical analyses

We compared the frequency of EMs per square-root specimens among bioregions and against species functional traits with a mixed generalized linear model [glmmTMB function of the “glmmTMB” R package (76)], we used a tweedie family allowing the search of the most appropriate distribution of the response variable. Species membership was included as a random factor nested with genus membership to account for phylogenetic autocorrelation. Post hoc tests among bioregions were applied using the function emmeans from the “emmeans” R package (77) with the Tukey test option. The proportion of variance explained by bioregions and traits (marginal R2) and that explained by the full model including the species random factor (conditional R2) were obtained by using the r.squaredGLMM function of the “MuMIn” package (78).

For the analysis of phylogeographic breaks, we only used Delaunay triangulation segments showing the least cost path shorter than 1000 km to avoid cases where the midpoints of very long paths were located too far from the area where specimens were collected. The following predictors were tested for their effects on the p-distances between specimens in a GAMM [gamm function of the “mgcv” R package (79)] using Delaunay triangulation segments as cases: (i) the maximum elevation along the path (smooth variable), (ii) the log transformed length of the sea strait (smooth variable), (iii) the log transformed length of the path (smooth predictor), (iv) latitude, and (v) longitude. Prior macrogenetic studies have shown that most of the variance in IGV is not explained by latitude because of strong differences among species in their functional traits (mostly related to life span and mobility) (20). For this reason, we included the three previously described butterfly traits (80): Because many segments (cases) are produced for each species in genetic landscapes, species membership was included as a random factor. Spatial autocorrelation was also considered in the model by including a covariance structure with distances based on midpoints (see R script in Dryad). Because the p-distances were not normally distributed and right skewed, we used a tweedie family density method.

The latitudinal trend of genetic diversity was modeled using GLMM (81) (glmmTMB function). As a response variable, we used the IGV calculated for each band (inter- and intraregion separately), and as predictors, we used the central latitude value of each band, the type of IGV (intraregion or inter-region), the geographic distances, and the three species traits previously described. We also included the number of specimens for each species sequenced from each band to examine possible effects of sample size on IGV values. Four interaction terms between central latitude and the type of IGV and the three species traits were also included. Last, we included two random factors: species membership (because latitudinal effects are expected to vary widely among species) nested within genus membership (to account for phylogenetic autocorrelation) and latitudinal bands (because two measures, intra- and inter-region IGV, are available for each band). We also included a covariance term in the model based on geographic distances among the specimens (calculated using barycenter coordinates) to account for spatial autocorrelation as indicated in the glmmTMB package (see R script in Dryad). Because the distribution of IGV values was highly right skewed, we used a tweedie family. Last, we checked for multicollinearity between the predictors (latitude and species traits) using a Pearson correlation. We also conducted two GLMMs for inter- and intraregion IGV separately to better understand these different component of latitudinal diversity trends. The proportions of variance explained by both fixed and random factors (conditional R2) and that explained by fixed factors only (marginal R2) were obtained by using the r.squaredGLMM function of the “MuMIn” package.

Acknowledgments

We thank many lepidopterists, too numerous to be mentioned individually, who provided samples, as well as many researchers who supported COI sequencing. We also thank the anonymous reviewers for constructive comments that helped improve the manuscript.

Funding: This work was supported by the following: National Biodiversity Future Center (NBFC, Italian Ministry of University and Research, PNRR, Missione 4 Componente 2, “Dalla ricerca all’impresa,” Investimento 1.4, Project CN00000033) (L.D.); Horizon project “Biodiversity Genomics of Europe” (L.D.); Direttiva Biodiversità projects 2019-2020 of the Italian Ministero della Transizione Ecologica (L.D.); grant PID2020-117739GA-I00 MCIN / AEI / 10.13039/501100011033 (G.T.); grant 2021-SGR-01334, Generalitat de Catalunya (G.T.); Fellowship FPU19/01593 (A.G.-B.); grant PID2022-139689NB-I00 by MICIU/ AEI/ 10.13039/501100011033 and by ERDF, EU (R.V.); grant 2021-SGR-00420 by Departament de Recerca i Universitats, Generalitat de Catalunya (R.V.); the “la Caixa” Foundation (ID 100010434, grant LCF/BQ/DR20/11790020) (M.M.); Academy of Finland (Academy Research Fellow, decision nos. 324988, 328895, and 352652) (V.D.); and The Canada Foundation for Innovation and NSERC (P.D.N.H.).

Author contributions: Conceptualization: L.D. and R.V. Methodology: L.D., R.V., M.M., J.D., A.G.-B., and G.T. Data curation: L.D., R.V., V.D., M.M., P.D.N.H., J.D., and G.T. R and Python coding: L.D. and A.G.-B. Investigation: L.D., R.V., M.M., J.D., A.G.-B., and G.T. Visualization: L.D. and M.M. Supervision: L.D., M.M., V.D., and R.V. Writing—original draft: L.D. and R.V. Writing—review and editing: L.D., R.V., V.D., M.M., P.D.N.H., J.D., G.T., and A.G.-B.

Competing interests: The authors declare they have no competing interests.

Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The data that support the findings of this study (COI sequences, specimens metadata, and layers necessary to carry out the analyses) together with R functions and R scripts to replicate the analysis are available in the Dryad repository (DOI: 10.5061/dryad.73n5tb33n) and in the GLQ GitHub repository (https://github.com/leondap/GLQ). All the sequences belong to BOLD and GenBank repositories, and their accession codes are reported in the metadata.

Supplementary Materials

The PDF file includes:

Figs. S1 to S8

Tables S1 to S5

Legend for data S1

Other Supplementary Material for this manuscript includes the following:

Data S1
==== Refs
REFERENCES AND NOTES

1 G. Hewitt, The genetic legacy of the Quaternary ice ages. Nature 405 , 907–913 (2000).10879524
2 P. Taberlet, L. Fumagalli, A. Wust-Saucy, J. Cosson, Comparative phylogeography and postglacial colonization routes in Europe. Mol. Ecol. 7 , 453–464 (1998).9628000
3 T. Schmitt, Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front. Zool. 4 , 11 (2007).17439649
4 A. Papadopoulou, L. L. Knowles, Toward a paradigm shift in comparative phylogeography driven by trait-based hypotheses. Proc. Natl. Acad. Sci. U.S.A. 113 , 8018–8024 (2016).27432974
5 L. Dapporto, A. Cini, R. Vodă, V. Dincă, M. Wiemers, M. Menchetti, G. Magini, G. Talavera, T. Shreeve, S. Bonelli, L. P. Casacci, E. Balletto, S. Scalercio, R. Vila, Integrating three comprehensive data sets shows that mitochondrial DNA variation is linked to species traits and paleogeographic events in European butterflies. Mol. Ecol. Resour. 19 , 1623–1636 (2019).31325412
6 S. V. Edwards, V. V. Robin, N. Ferrand, C. Moritz, The evolution of comparative phylogeography: Putting the geography (and more) into comparative population genomics. Genome Biol. Evol. 14 , evab176 (2022).34347070
7 J. C. Avise, J. Arnold, R. M. Ball, E. Bermingham, T. Lamb, J. E. Neigel, C. A. Reeb, N. C. Saunders, Intraspecific phylogeography: The mitochondrial DNA bridge between population genetics and systematics. Annu. Rev. Ecol. Syst. 18 , 489–522 (1987).
8 J. C. Avise, Phylogeography: Retrospect and prospect. J. Biogeogr. 36 , 3–15 (2009).
9 R. J. Petit, I. Aguinagalde, J.-L. de Beaulieu, C. Bittkau, S. Brewer, R. Cheddadi, R. Ennos, S. Fineschi, D. Grivet, M. Lascoux, A. Mohanty, G. Müller-Starck, B. Demesure-Musch, A. Palmé, J. P. Martín, S. Rendell, G. G. Vendramin, Glacial refugia: Hotspots but not melting pots of genetic diversity. Science 300 , 1563–1565 (2003).12791991
10 J. R. Stewart, A. M. Lister, I. Barnes, L. Dalén, Refugia revisited: Individualistic responses of species in space and time. Proc. Biol. Sci. 277 , 661–671 (2010).19864280
11 C. M. French, L. D. Bertola, A. C. Carnaval, E. P. Economo, J. M. Kass, D. J. Lohman, K. A. Marske, R. Meier, I. Overcast, A. J. Rominger, P. P. A. Staniczenko, M. J. Hickerson, Global determinants of insect mitochondrial genetic diversity. Nat. Commun. 14 , 5276 (2023).37644003
12 T. Kuhn, The Structure of Scientific Revolutions (Princeton Univ. Press, 2021).
13 L. Dapporto, M. Menchetti, R. Vodă, C. Corbella, S. Cuvelier, I. Djemadi, M. Gascoigne-Pees, J. C. Hinojosa, N. T. Lam, M. Serracanta, G. Talavera, V. Dincă, R. Vila, The atlas of mitochondrial genetic diversity for Western Palaearctic butterflies. Glob. Ecol. Biogeogr. 31 , 2184–2190 (2022).
14 B. G. Holt, J.-P. Lessard, M. K. Borregaard, S. A. Fritz, M. B. Araújo, D. Dimitrov, P.-H. Fabre, C. H. Graham, G. R. Graves, K. A. Jønsson, D. Nogués-Bravo, Z. Wang, R. J. Whittaker, J. Fjeldså, C. Rahbek, An update of Wallace’s zoogeographic regions of the world. Science 339 , 74–78 (2013).23258408
15 L. Dapporto, G. Ciolli, R. L. H. Dennis, R. Fox, T. G. Shreeve, A new procedure for extrapolating turnover regionalization at mid-small spatial scales, tested on British butterflies. Methods Ecol. Evol. 6 , 1287–1297 (2015).
16 M. Chavent, V. Kuentz-Simonet, A. Labenne, J. Saracco, ClustGeo: An R package for hierarchical clustering with spatial constraints. Comput. Stat. 33 , 1799–1822 (2018).
17 M. Menchetti, G. Talavera, A. Cini, V. Salvati, V. Dincă, L. Platania, S. Bonelli, E. Balletto, R. Vila, L. Dapporto, Two ways to be endemic. Alps and Apennines are different functional refugia during climatic cycles. Mol. Ecol. 30 , 1297–1310 (2021).33421216
18 T. Schmitt, U. Fritz, M. Delfino, W. Ulrich, J. C. Habel, Biogeography of Italy revisited: Genetic lineages confirm major phylogeographic patterns and a pre-Pleistocene origin of its biota. Front. Zool. 18 , 34 (2021).34187502
19 L. Dapporto, A. Cini, M. Menchetti, R. Vodă, S. Bonelli, L. P. Casacci, V. Dincă, S. Scalercio, J. C. Hinojosa, H. Biermann, L. Forbicioni, U. Mazzantini, L. Venturi, F. Zanichelli, E. Balletto, T. G. Shreeve, R. L. H. Dennis, R. Vila, Rise and fall of island butterfly diversity: Understanding genetic differentiation and extinction in a highly diverse archipelago. Divers. Distrib. 23 , 1169–1181 (2017).
20 D. M. Leigh, C. B. van Rees, K. L. Millette, M. F. Breed, C. Schmidt, L. D. Bertola, B. K. Hand, M. E. Hunter, E. L. Jensen, F. Kershaw, L. Liggins, G. Luikart, S. Manel, J. Mergeay, J. M. Miller, G. Segelbacher, S. Hoban, I. Paz-Vinas, Opportunities and challenges of macrogenetic studies. Nat. Rev. Genet. 22 , 791–807 (2021).34408318
21 G. N. Feliner, Southern European glacial refugia: A tale of tales. Taxon 60 , 365–372 (2011).
22 S. Arif, M. Gerth, W. G. Hone-Millard, M. D. S. Nunes, L. Dapporto, T. G. Shreeve, Evidence for multiple colonisations and Wolbachia infections shaping the genetic structure of the widespread butterfly Polyommatus icarus in the British Isles. Mol. Ecol. 30 , 5196–5213 (2021).34402109
23 J. Marcus, W. Ha, R. F. Barber, J. Novembre, Fast and flexible estimation of effective migration surfaces. eLife 10 , e61927 (2021).34328078
24 B. Burriel-Carranza, H. Tejero-Cicuéndez, A. Carné, G. Riaño, A. Talavera, S. Al Saadi, J. Els, J. Šmíd, K. Tamar, P. Tarroso, S. Carranza, The origin of a mountain biota: Hyper-aridity shaped reptile diversity in an Arabian biodiversity hotspot. bioRxiv 2023.04.07.536010 [Preprint] (2023). 10.1101/2023.04.07.536010.
25 S. Theodoridis, D. A. Fordham, S. C. Brown, S. Li, C. Rahbek, D. Nogues-Bravo, Evolutionary history and past climate change shape the distribution of genetic diversity in terrestrial mammals. Nat. Commun. 11 , 2557 (2020).32444801
26 A. Miraldo, S. Li, M. K. Borregaard, A. Flórez-Rodríguez, S. Gopalakrishnan, M. Rizvanovic, Z. Wang, C. Rahbek, K. A. Marske, D. Nogués-Bravo, An Anthropocene map of genetic diversity. Science 353 , 1532–1535 (2016).27708102
27 J. C. Habel, W. Ulrich, T. Assmann, Allele elimination recalculated: Nested subset analyses for molecular biogeographical data. J. Biogeogr. 40 , 769–777 (2013).
28 J. C. Habel, C. Drees, T. Schmitt, T. Assmann, “Review refugial areas and postglacial colonizations in the western Palearctic” in Relict Species: Phylogeography and Conservation Biology, J. C. Habel, T. Assmann, Eds. (Springer, 2010), pp. 189–197.
29 J. M. Waters, C. I. Fraser, G. M. Hewitt, Founder takes all: Density-dependent processes structure biodiversity. Trends Ecol. Evol. 28 , 78–85 (2013).23000431
30 V. Dincă, L. Dapporto, P. Somervuo, R. Vodă, S. Cuvelier, M. Gascoigne-Pees, P. Huemer, M. Mutanen, P. D. N. Hebert, R. Vila, High resolution DNA barcode library for European butterflies reveals continental patterns of mitochondrial genetic diversity. Commun. Biol. 4 , 315 (2021).33750912
31 L. Dapporto, C. Bruschini, D. Baracchi, A. Cini, S. F. Gayubo, J. A. GonzÁlez, R. L. H. Dennis, Phylogeography and counter-intuitive inferences in island biogeography: Evidence from morphometric markers in the mobile butterfly Maniola jurtina (Linnaeus) (Lepidoptera, Nymphalidae). Biol. J. Linn. Soc. 98 , 677–692 (2009).
32 S. Scalercio, A. Cini, M. Menchetti, R. Vodă, S. Bonelli, A. Bordoni, L. P. Casacci, V. Dincă, E. Balletto, R. Vila, L. Dapporto, How long is 3 km for a butterfly? Ecological constraints and functional traits explain high mitochondrial genetic diversity between Sicily and the Italian Peninsula. J. Anim. Ecol. 89 , 2013–2026 (2020).32207150
33 L. Dapporto, S. Fattorini, R. Vodǎ, V. Dincǎ, R. Vila, Biogeography of western Mediterranean butterflies: Combining turnover and nestedness components of faunal dissimilarity. J. Biogeogr. 41 , 1639–1650 (2014).
34 C. W. Burney, R. T. Brumfield, Ecology predicts levels of genetic differentiation in Neotropical birds. Am. Nat. 174 , 358–368 (2009).19627230
35 V. Devictor, C. Van Swaay, T. Brereton, L. Brotons, D. Chamberlain, J. Heliölä, S. Herrando, R. Julliard, M. Kuussaari, Å. Lindström, J. Reif, D. B. Roy, O. Schweiger, J. Settele, C. Stefanescu, A. Van Strien, C. Van Turnhout, Z. Vermouzek, M. WallisDeVries, I. Wynhoff, F. Jiguet, Differences in the climatic debts of birds and butterflies at a continental scale. Nat. Clim. Chang. 2 , 121–124 (2012).
36 C. Schmidt, C. J. Garroway, The conservation utility of mitochondrial genetic diversity in macrogenetic research. Conserv. Genet. 22 , 323–327 (2021).
37 D. E. Soltis, A. B. Morris, J. S. McLachlan, P. S. Manos, P. S. Soltis, Comparative phylogeography of unglaciated eastern North America. Mol. Ecol. 15 , 4261–4293 (2006).17107465
38 T. Schmitt, Z. Varga, Extra-Mediterranean refugia: The rule and not the exception? Front. Zool. 9 , 22 (2012).22953783
39 M. Husemann, T. Schmitt, F. E. Zachos, W. Ulrich, J. C. Habel, Palaearctic biogeography revisited: Evidence for the existence of a North African refugium for Western Palaearctic biota. J. Biogeogr. 41 , 81–94 (2014).
40 H. J. B. Birks, K. J. Willis, Alpines, trees, and refugia in Europe. Plant Ecol. Divers. 1 , 147–160 (2008).
41 R. Ohlemüller, B. J. Anderson, M. B. Araújo, S. H. M. Butchart, O. Kudrna, R. S. Ridgely, C. D. Thomas, The coincidence of climatic and species rarity: High risk to small-range species from climate change. Biol. Lett. 4 , 568–572 (2008).18664421
42 L. Dapporto, Satyrinae butterflies from Sardinia and Corsica show a kaleidoscopic intraspecific biogeography (Lepidoptera, Nymphlidae). Biol. J. Linn. Soc. 100 , 195–212 (2010).
43 M. R. Tarrier, J. Delacre, Les Papillons de Jour Du Maroc: Guide d’identification et de Bio-Indication (Biotope, 2008).
44 G. M. Hewitt, Post-glacial re-colonization of European biota. Biol. J. Linn. Soc. 68 , 87–112 (1999).
45 H. Romo, E. García-Barros, R. J. Wilson, R. G. Mateo, M. L. Munguira, Modelling the scope to conserve an endemic-rich mountain butterfly taxon in a changing climate. Insect Conserv. Divers. 16 , 451–467 (2023).
46 G. Sistri, M. Menchetti, L. Santini, L. Pasquali, S. Sapienti, A. Cini, L. Platania, E. Balletto, F. Barbero, S. Bonelli, L. P. Casacci, V. E. Dinca, R. Vila, C. Mantoni, S. Fattorini, L. Dapporto, The isolated Erebia pandrose Apennine population is genetically unique and endangered by climate change. Insect Conserv. Divers. 15 , 136–148 (2022).
47 A. Gómez, D. H. Lunt, “Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula” in Phylogeography of Southern European Refugia, S. Weiss, N. Ferrand, Eds. (Springer, 2007), pp. 155–188.
48 J. L. Hernández-Roldán, L. Dapporto, V. Dincă, J. C. Vicente, E. A. Hornett, J. Šíchová, V. A. Lukhtanov, G. Talavera, R. Vila, Integrative analyses unveil speciation linked to host plant shift in Spialia butterflies. Mol. Ecol. 25 , 4267–4284 (2016).27393640
49 V. Dincǎ, L. Dapporto, R. Vila, A combined genetic-morphometric analysis unravels the complex biogeographical history of Polyommatus icarus and Polyommatus celina common blue butterflies. Mol. Ecol. 20 , 3921–3935 (2011).21883579
50 M. W. Jacobsen, N. W. Jensen, R. Nygaard, K. Præbel, B. Jónsson, N. H. Nielsen, J. M. Pujolar, D. J. Fraser, L. Bernatchez, M. M. Hansen, A melting pot in the Arctic: Analysis of mitogenome variation in Arctic char (Salvelinus alpinus) reveals a 1000-km contact zone between highly divergent lineages. Ecol. Freshw. Fish 31 , 330–346 (2022).
51 J. C. Habel, T. Schmitt, P. Müller, The fourth paradigm pattern of post-glacial range expansion of European terrestrial species: The phylogeography of the Marbled White butterfly (Satyrinae, Lepidoptera). J. Biogeogr. 32 , 1489–1497 (2005).
52 M. Mutanen, A. Hausmann, P. D. N. Hebert, J.-F. Landry, J. R. de Waard, P. Huemer, Allopatry as a Gordian knot for taxonomists: Patterns of DNA barcode divergence in Arctic-Alpine Lepidoptera. PLOS ONE 7 , e47214 (2012).23071761
53 J. Davis, A. Pavlova, R. Thompson, P. Sunnucks, Evolutionary refugia and ecological refuges: Key concepts for conserving Australian arid zone freshwater biodiversity under climate change. Glob. Chang. Biol. 19 , 1970–1984 (2013).23526791
54 J. C. Habel, R. Vila, R. Vodă, M. Husemann, T. Schmitt, L. Dapporto, Differentiation in the marbled white butterfly species complex driven by multiple evolutionary forces. J. Biogeogr. 44 , 433–445 (2017).
55 V. Dincǎ, C. Wiklund, V. A. Lukhtanov, U. Kodandaramaiah, K. Norén, L. Dapporto, N. Wahlberg, R. Vila, M. Friberg, Reproductive isolation and patterns of genetic differentiation in a cryptic butterfly species complex. J. Evol. Biol. 26 , 2095–2106 (2013).23909947
56 M. S. Tahami, V. Dincă, K. M. Lee, R. Vila, M. Joshi, M. Heikkilä, L. Dapporto, S. Schmid, P. Huemer, M. Mutanen, Genomics reveal admixture and unexpected patterns of diversity in a parapatric pair of butterflies. Genes 12 , 2009 (2021).34946956
57 D. P. L. Toews, A. Brelsford, The biogeography of mitochondrial and nuclear discordance in animals. Mol. Ecol. 21 , 3907–3930 (2012).22738314
58 A. Duplouy, E. A. Hornett, Uncovering the hidden players in Lepidoptera biology: The heritable microbial endosymbionts. PeerJ 6 , e4629 (2018).29761037
59 I. Paz-Vinas, E. L. Jensen, L. D. Bertola, M. F. Breed, B. K. Hand, M. E. Hunter, F. Kershaw, D. M. Leigh, G. Luikart, J. Mergeay, J. M. Miller, C. B. Van Rees, G. Segelbacher, S. Hoban, Macrogenetic studies must not ignore limitations of genetic markers and scale. Ecol. Lett. 24 , 1282–1284 (2021).33749962
60 K. L. Millette, V. Fugère, C. Debyser, A. Greiner, F. J. J. Chain, A. Gonzalez, No consistent effects of humans on animal genetic diversity worldwide. Ecol. Lett. 23 , 55–67 (2020).31637822
61 E. P. Lessa, J. A. Cook, J. L. Patton, Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary. Proc. Natl. Acad. Sci. U.S.A. 100 , 10331–10334 (2003).12913123
62 M. Wiemers, E. Balletto, V. Dincă, Z. F. Fric, G. Lamas, V. Lukhtanov, M. L. Munguira, C. A. M. VanSwaay, R. Vila, A. Vliegenthart, N. Wahlberg, R. Verovnik, An updated checklist of the European butterflies (Lepidoptera, papilionoidea). Zookeys 811 , 9–45 (2018).
63 B. A. Hawkins, Multiregional comparison of the ecological and phylogenetic structure of butterfly species richness gradients. J. Biogeogr. 37 , 647–656 (2010).
64 C. Parmesan, N. Ryrholm, C. Stefanescu, J. K. Hill, C. D. Thomas, H. Descimon, B. Huntley, L. Kaila, J. Kullberg, T. Tammaru, W. J. Tennent, J. A. Thomas, M. Warren, Poleward shifts in geographical ranges of butterfly species associated with regional warming. Nature 399 , 579–583 (1999).
65 M. T. Burrows, D. S. Schoeman, L. B. Buckley, P. Moore, E. S. Poloczanska, K. M. Brander, C. Brown, J. F. Bruno, C. M. Duarte, B. S. Halpern, J. Holding, C. V. Kappel, W. Kiessling, M. I. O’Connor, J. M. Pandolfi, C. Parmesan, F. B. Schwing, W. J. Sydeman, A. J. Richardson, The pace of shifting climate in marine and terrestrial ecosystems. Science 334 , 652–655 (2011).22053045
66 E. Paradis, J. Claude, K. Strimmer, APE: Analyses of phylogenetics and evolution in R language. Bioinformatics 20 , 289–290 (2004).14734327
67 L. Dapporto, M. Ramazzotti, S. Fattorini, G. Talavera, R. Vila, R. L. H. Dennis, Recluster: An unbiased clustering procedure for beta-diversity turnover. Ecography 36 , 1070–1075 (2013).
68 V. A. Lukhtanov, A. Sourakov, E. V. Zakharov, P. D. N. Hebert, DNA barcoding Central Asian butterflies: Increasing geographical dimension does not significantly reduce the success of species identification. Mol. Ecol. Resour. 9 , 1302–1310 (2009).21564901
69 J. D’Ercole, V. Dincă, P. A. Opler, N. Kondla, C. Schmidt, J. D. Phillips, R. Robbins, J. M. Burns, S. E. Miller, N. Grishin, E. V. Zakharov, J. R. DeWaard, S. Ratnasingham, P. D. N. Hebert, A DNA barcode library for the butterflies of North America. PeerJ 9 , e11157 (2021).33976967
70 R. L. H. Dennis, P. B. Hardy, L. Dapporto, Nestedness in island faunas: Novel insights into island biogeography through butterfly community profiles of colonization ability and migration capacity. J. Biogeogr. 39 , 1412–1426 (2012).
71 A. G. Vandergast, W. M. Perry, R. V. Lugo, S. A. Hathaway, Genetic landscapes GIS Toolbox: Tools to map patterns of genetic divergence and diversity. Mol. Ecol. Resour. 11 , 158–161 (2011).21429115
72 R. Turner, M. R. Turner, Package ‘deldir’ (2021).
73 J. van Etten, R package gdistance: Distances and routes on geographical grids. J. Stat. Softw. 76 , 1–21 (2017).36568334
74 M. Pontarp, L. Bunnefeld, J. S. Cabral, R. S. Etienne, S. A. Fritz, R. Gillespie, C. H. Graham, O. Hagen, F. Hartig, S. Huang, R. Jansson, O. Maliet, T. Münkemüller, L. Pellissier, T. F. Rangel, D. Storch, T. Wiegand, A. H. Hurlbert, The latitudinal diversity gradient: Novel understanding through mechanistic eco-evolutionary models. Trends Ecol. Evol. 34 , 211–223 (2019).30591209
75 P. Gratton, S. Marta, G. Bocksberger, M. Winter, P. Keil, E. Trucchi, H. Kühl, Which latitudinal gradients for genetic diversity? Trends Ecol. Evol. 32 , 724–726 (2017).28807398
76 M. E. Brooks, K. Kristensen, K. J. Van Benthem, A. Magnusson, C. W. Berg, A. Nielsen, H. J. Skaug, M. Machler, B. M. Bolker, glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9 , 378–400 (2017).
77 R. Lenth, H. Singmann, J. Love, P. Buerkner, M. Herve, Package ‘emmeans’ (2019).
78 K. Barton, MuMIn: Multi-model inference. R package version 1. 0. 0. (2009); https://r-forge.r-project.org/projects/mumin/.
79 S. Wood, M. S. Wood, Package ‘mgcv.’ R package version 1, 729 (2015).
80 J. Middleton-Welling, L. Dapporto, E. García-Barros, M. Wiemers, P. Nowicki, E. Plazio, S. Bonelli, M. Zaccagno, M. Šašić, J. Liparova, O. Schweiger, A. Harpke, M. Musche, J. Settele, R. Schmucki, T. Shreeve, A new comprehensive trait database of European and Maghreb butterflies, Papilionoidea. Sci. Data 7 , 351 (2020).33060594
81 H. de Kort, J. G. Prunier, S. Ducatez, O. Honnay, M. Baguette, V. M. Stevens, S. Blanchet, Life history, climate and biogeography interactively affect worldwide genetic diversity of plant and animal populations. Nat. Commun. 12 , 516 (2021).33483517
