The Pleistocene species pump past its prime: Evidence from European butterfly sister species

The Pleistocene glacial cycles had a profound impact on the ranges and genetic make‐up of organisms. While it is clear that the contact zones that have been described for many sister taxa are secondary and have formed in the current interglacial, it is unclear when the taxa involved began to diverge. Previous estimates based on small numbers of loci are unreliable given the stochasticity of genetic drift and the contrasting effects of incomplete lineage sorting and gene flow on gene divergence. Here, we use genome‐wide transcriptome data to estimate divergence for 18 sister species pairs of European butterflies showing either sympatric or contact zone distributions. We find that in most cases, species divergence predates the mid‐Pleistocene transition or even the entire Pleistocene period. We also show that although post‐divergence gene flow is restricted to contact zone pairs, they are not systematically younger than sympatric pairs. This suggests that contact zones are not limited to the initial stages of the speciation process, but can involve notably old taxa. Finally, we show that mitochondrial divergence and nuclear divergence are only weakly correlated and mitochondrial divergence is higher for contact zone pairs.


| INTRODUC TI ON
Divergence in allopatry provides a simple null model of speciation (Mayr, 1947). Following geographic isolation and given enough time, reproductive isolation is inevitable as incompatibilities will eventually become fixed as a result of genetic drift and/or selection (Bateson, 1909;Dobzhansky, 1937;Muller, 1942). Taxa that evolved partial reproductive isolation in allopatry may come into secondary contact as a result of range shifts and-depending on their degree of reproductive isolation and niche overlap-either form a contact zone or invade each other's range (Barton, 1985;Pigot, 2013). If allopatric divergence dominates speciation, then local alpha diversity for a given clade cannot accrue until secondary sympatry is achieved (Weir & Price, 2011). Thus, the forces that facilitate or hamper secondary sympatry and the timescale over which this occurs have profound consequences both for speciation and for the spatial distribution of species diversity. While modern ranges only provide a snapshot of the dynamic history of range shifts, understanding the extent to which current range overlap between closely related species can be explained by their speciation history and vice versa has been at the core of speciation research (Coyne & Orr, 2004).
The glacial cycles of the Pleistocene had a profound effect on current diversity of temperate ecosystems (Hewitt, 1996(Hewitt, , 2001Hofreiter, 2009). Populations of temperate taxa in Europe were isolated in icefree refugia around the Mediterranean basin (Iberia, Italy, the Balkans and the larger Mediterranean islands) as glaciers encroached. The observation that the geographic ranges of many young taxa are restricted to individual glacial refugia in southern Europe (Dennis et al., 1991;Hewitt, 1996Hewitt, , 1999Hewitt, , 2011Schmitt, 2007) suggests that this repeated separation into and expansion out of glacial refugia has played a major role in their origin. The availability of allozyme and mitochondrial (mt) data in the 80s and 90s has spurred an abundance of case studies on intra-and interspecific diversity of European taxa including detailed investigations of hybrid zones in taxa ranging from fire-bellied toads (Kruuk et al., 1999), the house mouse (Boursot et al., 1996), grasshoppers (Barton, 1980 andButlin &, to plants (Bacilieri et al., 1996) and marine mussels (Skibinski & Beardmore, 1979). The pervading evidence from these studies is that genetic diversity within and in, many cases, divergence between species is structured by refugia Hewitt, 1996;Schmitt, 2007).

| When was divergence between sister species initiated?
While it is clear that the hybrid zones we observe today are secondary contacts that formed after the last glacial maximum and may have formed many times over throughout the Pleistocene, it is far from clear when divergence between the sister taxa involved was initiated. One possibility is that the Pleistocene glacial cycles initiated species divergence directly by separating populations into allopatric refugia (i.e. a 'species pump' sensu Haffer, 1969). Another possibility is that the initial divergence between sister species predates the Pleistocene, and so, any build-up of reproductive isolation during the Pleistocene (e.g. via the fixation of intrinsic incompatibilities and/or reinforcement) occurred in populations that were already partially diverged. If the Pleistocene species pump hypothesis is correct, we would expect sister species divergence times to be concentrated during or at the beginning of the mid-Pleistocene transition 0.8-1.2 million years ago (MYA), which marks the onset of continent-wide glacial cycling (Bishop et al., 2011). The idea that Pleistocene divergence acted as a species pump was first proposed in the context of American faunas (Avise et al., 1998;Bernatchez & Wilson, 1998;Haffer, 1969), but has dominated phylogeographic studies on European sister taxa (e.g. Habel et al., 2008;Hewitt, 1996Hewitt, , 2000Schmitt, 2007;Schoville et al., 2012). In contrast, other studies including some of the early work on European contact zones (Barton & Hewitt, 1985;Butlin & Hewitt, 1985) conclude that the taxa involved in such secondary contacts may substantially predate the Pleistocene (Abbott et al., 2000;Hewitt, 1996;Klicka & Zink, 1997;Spooner & Ritchie, 2006). Similarly, Pleistocene climate forcing is insufficient in explaining divergence in an Amazonian butterfly suture zone (Dasmahapatra et al., 2010). Thus, it remains unclear to what extent divergence between sister taxa was initiated by 'Pleistocene species pump' dynamics or has an older, deeper origin?
A corollary for the hypothesis of allopatric speciation in different refugia is that range overlap is secondary. Since species can more easily invade each others' ranges once sufficient premating barriers and ecological differentiation have developed, we would expect species pairs with overlapping ranges to be older overall than those without range overlap, all else being equal (Coyne & Orr, 2004). Support for this prediction comes from comparative studies showing that the proportion of range overlap (degree of sympatry (Chesser & Zink, 1994)) is positively (albeit weakly) correlated with genetic divergence (Barraclough & Vogler, 2000;Pigot & Tobias, 2013).
However, a recent study in Chorthippus grasshoppers shows that subspecies that hybridize across contact zones can be older than currently sympatric species (Nolen et al., 2020).

| Mitonuclear discordance
Age estimates for recently diverged taxa have largely relied on singlelocus phylogenies and ignored incomplete lineage sorting. Hewitt (2011) summarizes age estimates for European hybrid-zones taxa including mammals, insects, amphibians and reptiles, which range from hundreds of thousands to several million years ago. However, given that these estimates are based on different markers and calibrations, the extent to which glacial cycles have initiated speciation events remains unknown. Estimates based on mitochondrial (mt) data are particularly unreliable for at least three reasons. First, the mutation rate of mtDNA is highly erratic (Galtier et al., 2009). Second, given the stochasticity of coalescence, the ancestry of a single locus (however well resolved) is a very poor measure of species divergence. In the absence of gene flow, divergence at a single locus may substantially predate the onset of species divergence, while it may be much more recent in the presence of gene flow (Knowles & Carstens, 2007;Wang & Hey, 2010). Mitonuclear discordance in both directions has been found in a large number of animal systems (Toews & Brelsford, 2012) including several closely related species of European butterflies Hinojosa et al., 2019;Wiemers et al., 2010). Finally, mtDNA does not evolve neutrally since transmission of mitochondria is completely linked to maternal inheritance of endosymbionts such as Wolbachia and Spiroplasma and, in organisms with Z/W sex determination, of the W chromosome. Thus, mt diversity and divergence may be driven largely by selective sweeps (including introgression sweeps) rather than neutral gene flow and genetic drift (Galtier et al., 2009;Hurst & Jiggins, 2005;Jiggins, 2003;Martin et al., 2020).

| European butterflies as a model group
Testing whether climate-induced Pleistocene range shifts have triggered speciation or patterned older splits between species requires replication both at the level of genetic loci and at the level of speciation events. Although we can now generate WGS data for any species, there are surprisingly few reliable estimates for the onset of divergence between European sister species and such estimates are lacking even for well-studied contact zone taxa (but see Nolen et al., 2020;Nürnberger et al., 2016).
Lepidoptera are arguably the best-studied arthropod family: European butterflies provide a unique opportunity to investigate divergence and speciation processes comparatively . Near-complete information on geographic ranges and key life history traits (e.g. voltinism and host plant range) is available (Kudrna, 2019;Tolman & Lewington, 2013). Additionally, the taxonomy of all 496 European species is well resolved (Wiemers et al., 2018) and a complete, multilocus phylogeny of all European taxa exists . This, combined with extensive DNA barcode reference libraries Dincӑ et al., 2021), facilitates the identification of species (especially in the case of cryptic taxa) and provides extensive sampling of sister species pairs, many of which abut at narrow contact zones (Dennis et al., 1991;Platania et al., 2020;Vodӑ et al., 2015) (Figure 1). Secondary contact zones have been described in detail for several European taxa, including Spialia orbifer and S. sertorius (Lorkovic, 1973), the Italian Pontia hybrid zone (Porter et al., 1997) and the contacts between Iphiclides podalirus and I. feisthamelii and between Melanargia galathea and M. lachesis along the Pyrenees (Habel et al., 2017;Wohlfahrt, 1996, Gaunet et al., 2019. Here, we use European butterflies as a model system to investigate to what extent the divergence times between sister species in this group are concentrated in the Pleistocene, as predicted by the Pleistocene species pump hypothesis, and test how well recent sister species fit a null model of divergence in allopatry. Although European butterflies have been studied intensively, with few exceptions (see Talla et al., 2017), robust estimates of divergence required for any systematic comparison of speciation are lacking. We generate RNA-seq data for 18 sister species pairs and ask the following specific questions: (i) Has speciation been initiated during the Pleistocene, as envisaged by the species pump hypothesis, or did the glacial cycles pattern pre-existing, older subdivisions?
(ii) Are sister species pairs that form contact zones younger than pairs that overlap in range?
(iii) Is there evidence for gene flow between contact zone species?
(iv) How strongly correlated are mitochondrial and nuclear divergence and do contact zone pairs show increased mitonuclear discordance?
F I G U R E 1 Nine of the 18 sister species pairs of butterfly in which we quantified genome-wide divergence meet at contact zones in southern Europe. In the left group, from left to right across northern Iberia are Satyrium, Pseudophilotes, Melanargia and Iphiclides. In the centre group, from bottom to top across the Alps are Pontia, Euchloe, Pyrgus, and Zerynthia. Finally, on the right across the Balkans is the genus Spialia 2 | ME THODS

| Sampling and molecular work
We identified true sister species pairs in the European butterfly phylogeny (Wiemers et al., 2018;. Species pairs involving island and mountain endemics were excluded, as these cannot achieve secondary sympatry. We also did not consider species pairs that are unlikely to have originated in Europe, for example sister pairs involving North American taxa. Following these criteria, we sampled 18 sister species pairs (Table 1). Our sampling includes 7.3% of European butterfly species (Wiemers et al., 2018) and almost all 'good' butterfly sister species pairs in Europe (Descimon & Mallet, 2009). Specimen identifications were confirmed for 14 species that are difficult to identify based on morphology (morphological characters are subtle and/or internal) but for which COI barcodes are diagnostic of morphoID using LepF/R primers (Hajibabaei et al., 2006) and existing reference databases (Dincӑ et al., 2021). We were unable to obtain fresh material for Erebia euryale and E. ligea, and Fabriciana adippe and F. niobe (two remaining sister pairs meeting our sampling criteria).  (Table S1).
To generate comparable data sets across all samples, Orthofinder v2.3.3 (Emms & Kelly, 2015) was used to cluster proteins into orthogroups. Orthogroups were labelled single-copy orthologues (SCOs) if one protein of each taxon was present. Genus single-copy orthologues (GSCOs) were diagnosed based on the presence of single-copy proteins within the focal pair. Protein sequences from each orthogroup were used to align corresponding DNA sequences using Translatorx v12.0 (Abascal et al., 2010).
Data were generated for 36 species (18 sister pairs) from five families. For 17 pairs, data were generated from 665 SCOs from high-quality transcriptomes (BUSCO scores >90%). For the pair of Zerynthia species (one of which, Zerynthia polyxena, was sampled as a larva), GSCOs (5000 orthologues) were used to avoid restricting the SCOs for other pairs. With the exception of the Zerynthia pair, all analyses are based on SCO to enforce consistent comparisons across pairs. While the SCO data set is much smaller than the pair GSCO data sets and likely enriched for conserved and highly expressed genes, this has very little impact on estimates of divergence and diversity at fourfold degenerate (4D) sites, as these are highly correlated ( Figure S1 and Mackintosh et al., 2019).

| Estimating gene and population divergence
For each species pair, we calculated the average pairwise gene divergence d xy at fourfold degenerate (4D) sites using sequence alignments for one or two diploid samples from each species. This calculation is implemented in the script orthodiver.py (www.github. com/sameb don/ortho diver).
Species split times were estimated using two different approaches.
First, we used a simple rescaling of mean genetic divergence and diversity to obtain a lower bound of the species divergence time. Assuming a simple null model of divergence without gene flow, neutrality and an infinite site mutation model, net mean divergence d a = d xy − (Nei & Li, 1979) is directly proportional to species divergence time Here, is the de novo mutation rate per generation (per base). We assumed = 2.9 * 10 − 9 , an estimate of the spontaneous mutation rate obtained from parent-offspring trios of South American Heliconius melpomene butterflies (Keightley, Pinharanda, et al., 2014).
Since both violations of the mutation model (back-mutations) and the demographic model (gene flow) reduce d a , this time estimate is a lower bound of the true species divergence time. We converted estimates of species divergence time (T) into years ( ) using the mean generation time of each pair (Table 1). Information on generation times was compiled from Collins Butterfly Guide (Tolman & Lewington, 2013) ( Table 1). For species in which generation times vary with latitude, we assumed the minimum generation time of the southern part of the range. This is a reasonable long-term average, given that European glacial refugia are located around the Mediterranean.
Second, we estimated species split times using a multilocus coalescent approach. We considered the distribution of pairwise differences in blocks of a fixed length of 4D sites to fit an isolation-with-migration (IM) model and a nested history of strict divergence to each species pair. In the absence of recombination within blocks, the distribution of pairwise differences has been derived analytically (Lohse et al., 2011;Wilkinson-Herbots, 2012 However, given the high rate of recombination (relative to mutation) in butterflies (Martin et al., 2016(Martin et al., , 2019 and the substantial span of 4D blocks, we expect multilocus inference (assuming no withinblock recombination) to be biased. In particular, recombination is known to lead to (upwardly) biased estimates of divergence time (Wall, 2003).
Given this and other limitations (see Discussion), we will focus on the more simple and robust estimates of species divergence based on d a throughout.
As an additional test for gene flow, we compared the observed distributions of pairwise differences with analytic expectations under a model of strict divergence without gene flow given the estimates of T and ancestral N e obtained from d a and mean (which cannot be affected by recombination).
Thus, in the absence of gene flow, we would expect the empirical distribution of pairwise differences to be narrower than the analytic expectation (Wall, 2003) due to recombination. In contrast, wider-than-expected distributions are indicative of post-divergence gene flow. We re-sampled (without replacement) 10,000 data sets of equal size as the observed data from the expected distribution of each species and tested whether the likelihood of the observed distribution of pairwise differences falls within the distribution of likelihoods expected under a null model of strict divergence.

| Estimating range size and overlap
Geographic ranges were quantified as follows: we obtained occurrence data over Europe for all focal species with a resolution of 60′ latitude and 30′ longitude by critically revising the data from the Distribution Atlas of European Butterflies and Skippers (Kudrna et al., 2011) and by adding data from Roger Vila's collection stored at Institut de Biologia Evolutiva (Barcelona) (Figures S2-S4). To calculate range overlap, we applied the biodecrypt function (Platania et al., 2020) of the recluster R package . This function computes alpha hull with a given concavity ( ) and evaluates the area of overlap among pairs of species. We used = 2 and = 3 for species with discontinuous and continuous distributions in Europe, respectively. We quantified the range overlap of each species pair and calculated the degree of sympatry as: that is the fraction of the smaller range involved in the overlap. In the following, we consider sister pairs with a degree of sympatry ≤0.2 contact zone pairs and those with a degree of sympatry >0.2 sympatric.
Based on this, we classified nine pairs as contact zone taxa. However, since there are only two species pairs with intermediate levels of sympatry (>0.2 and <0.7), our comparisons of contact zone and sympatric pairs are robust to a wide range of thresholds.

| Mitochondrial diversity and divergence
Sequence alignments for the COI barcode locus were obtained from the BOLD database (Ratnasingham & Hebert, 2007) for all 18 sister species pairs. These sequences, together with associated information, are publicly available in the data set DS-EUGENMAP (dx.doi. org/10.5883/DS-EUGENMAP) on BOLD at www.bolds ystems.org and were originally produced by Dincӑ et al. (2021). For each species, we included all available sequence records from Europe (ranging from 21 in E. crameri to 429 in L. sinapis (Table S1)). Mean pairwise diversity ( ) within species and divergence (d xy ) across all sites were computed using DnaSP (Librado & Rozas, 2009).
We obtained the average gene divergence time for each pair from the multilocus-calibrated phylogeny of European butterflies of Wiemers et al., (2020) as half of patristic distances calculated with distTips function of the adephylo R package (Jombart & Dray, 2010).
The correlation between our estimates of species divergences and these node ages was explored with standardized major axis (SMA) regression, using the 'sma' function of the 'smatr' R package. SMA estimates slope and intercept and tests whether slope differs from one.

| Most European butterfly sister species predate the Pleistocene
Mean gene divergence (d xy ) at 4D sites between sister species ranged from 0.015 to 0.085, with a mean of 0.047 (Table 1, Figure 2 (Table S3). However, we find that sampling blocks of a fixed length resulted in a consistent downward bias in d xy of on average 10%. This is likely a result of selecting for long and likely conserved transcripts. In the light of this, and given that blocks based on 4D sites violate other key assumptions of multilocus inference, in particular, no recombination within loci and known phase, we caution against over-interpreting these model-based estimates (see Discussion) and focus on the simpler and more robust estimates of sister species divergence based on d a throughout.

| Sister pairs that form contact zones are not significantly younger than sympatric pairs
There are two reasons to expect species pairs that form contact zones to be younger than sympatric pairs: first, if speciation under a null model of divergence in allopatry is initiated by periods of vicariance, the formation of a contact zone (parapatry) represents an earlier stage in the transition to complete reproductive isolation and substantial range overlap (sympatry) (Coyne & Orr, 2004

| Evidence for recent gene flow in some contact zone pairs
The empirical distribution of pairwise differences deviates significantly from the expectation under a strict divergence model in a

| Pervasive mitonuclear discordance in contact zone species pairs
Our estimates of species divergence are based on average net divergence (d a ) across many hundreds of genes and are robust to how orthologues are filtered ( Figure S1). Given that previous studies on European butterflies have been largely based on mitochondrial (mt) phylogenies, an obvious question is to what extent mt divergence is correlated with mean nuclear divergence. We find that both d a and d xy at COI are positively but only weakly correlated with mean nuclear divergence (Figure 4). The correlation is weaker for d a than for d xy (R 2 = 0.27 and 0.31, respectively), which is unsurprisingly given that mitochondrial diversity (and hence d a ) is both inherently random and may be affected by selective sweeps. Comparing the relation between mt and nuclear d a between contact zone and sympatric pairs, we find a shallower slope for contact zone pairs (0.29 compared to 0.99; Figure 4). This difference, although not significant (p = 0.09), appears to be driven by the reduced mt diversity in contact zone compared with sympatric pairs (mean = 0.0030, SD = 0.0014 and = 0.0047, SD = 0.0031, respectively; t = 1.5763, df = 11.324, p = 0.0712). This suggests that mt diversity may be more strongly affected by selective sweeps in contact zone species than in sympatric pairs. We find no corresponding signal for any difference in nuclear diversity between contact zones and sympatric pairs (t = −0.0139, df = 31.539, p = 0.506) and, in general, no correlation between nuclear and mt diversity ( Figure S7 and Mackintosh et al., (2019)).
Our estimates for the lower bound of sister species divergence differ substantially from the ages of the corresponding nodes in the Wiemers et al., (2020) phylogeny for individual pairs ( Figure S8). This is unsurprising given that the latter are largely informed by mtDNA data ( Figure S9). However, perhaps surprisingly (given the difference in calibration, data and inference approach), our estimates are not consistently older or younger than the node ages of Wiemers et al.,

| Genetic diversity does not correlate with relative range size
Genetic diversity at 4D sites within all 36 species ranged from 0.32% to 4.2% with a mean of 1.5%. Given the H. melpomene mutation rate of = 2.9 × 10 − 9 (Keightley, Pinharanda, et al., 2014), these correspond to effective population sizes ranging from 280,000 to 3,600,000 with a mean of 1,300,000 (assuming = 4N e mu).

Mackintosh et al. (2019) tested whether neutral genetic diversity
across European butterflies correlates with geographic range and found no significant relation across 38 taxa. Our sampling of species pairs allows for a simpler, alternative test of the potential relationship between diversity and range size using sister-clade comparisons, which is less sensitive to potential phylogenetic correlates and uncertainty in current range estimates. If diversity is a function of range size, we expect the species in a pair with the larger range to have higher genetic diversity than the species with the smaller range. We indeed find a difference in the expected direction, 0.0167 (SD = 0.0114) vs 0.0139 (SD = 0.00865), although the effect of relative range size is not significant (t paired = 1.127, df = 17, p = 0.138; Figure S10).

| DISCUSS ION
We quantify and compare genome-wide divergence across 18 sister species pairs of European butterfly. Simple estimates for the onset of species divergence based on net gene divergence (d a ) and a direct mutation rate estimate for butterflies suggest that the majority of pairs have diverged before the onset of the major Pleistocene glacial F I G U R E 2 Species divergence time estimates ( ) plotted against mean genetic divergence (d xy ) for 18 European butterfly sister species pairs. Pairs which abut at contact zones (degree of sympatry ≤ 0.2) are shown in yellow, sympatric pairs with substantial range overlap (> 0.2) in blue. The vertical dashed line represents the beginning of the Pleistocene (2.6 MYA), and the vertical grey bar indicates the mid-Pleistocene transition (0.8-1.2 MYA).
The horizontal dotted lines represent the 95% confidence intervals of estimates given the uncertainty in the mutation rate used for calibration (1.3-5.9 × 10 − 9 (Keightley, Pinharanda, et al., 2014)). The temperature data (5-point running means of global surface temperature) are taken from Hansen et al., (2013) F I G U R E 3 The distribution of pairwise differences, S in blocks of a fixed length of 4D sites in contact (upper box) and sympatric (lower box) pairs. The observed distribution in single-copy orthologues is shown in orange, and the expectation under a history of strict divergence (estimated from and d a ), in grey. Pairs that show wider-than-expected distributions are marked with an asterisk (*), and species that show narrower-than-expected distributions are marked with a plus (+) cycles. Our results support the notion that the modern contact zones are secondary between species that began to diverge much earlier, in the Pliocene or early Pleistocene. Thus, even though the current ranges of many taxon pairs reflect glacial refugia, their initial divergence during the Pliocene or early Pleistocene is unlikely to have been triggered by repeated cycles of range connectivity and vicariance into refugia, as envisaged by the species pump hypothesis and phylogeographic studies based on mt and allozyme data (e.g. Habel et al., 2008;Lai & Pullin, 2004;Schmitt, 2007;Todisco et al., 2010;Zinetti et al., 2013) because substantial glaciation across continents did not develop (Bishop et al., 2011) until the 'mid-Pleistocene transition' from 0.8 to 1.2 MYA and the shift from 41,000-to 100,000- year glacial cycles. Given the antiquity of most sister species, it is perhaps unsurprising that we do not find any relationship between current range overlap and the time since divergence. Specifically, species pairs that form contact zones are not significantly younger than pairs that broadly overlap in range. However, we do find that strong signals of post-divergence gene flow are restricted to contact zone pairs. It is likely that the absence of sympatric pairs with significant gene flow reflects a simple survivorship bias: any incipient species pairs with significant gene flow might have already collapsed.
Likewise, we are more likely to observe old contact zones pairs that have survived repeated glacial cycles.
Our finding that mt divergence between sister species is only weakly correlated with mean nuclear divergence and the possibility that net mt divergence may be greater for contact zone than sympatric species pairs as a result of reduced genetic diversity (note that the differences in d a and mean between contact zone and sympatric pairs are marginally non-significant, p = 0.09 and 0.07, respectively) could suggest that contact zone species may be subject to more frequent selective sweeps linked to mitochondria. Such sweeps may be acting on mt variation directly or indirectly through maternally inherited genomes or chromosomes (e.g. Wolbachia (Jiggins, 2003) and the W chromosome) and have been documented in a number of Lepidopteran systems (Gaunet et al., 2019;Graham & Wilson, 2012;Kodandaramaiah et al., 2013;Martin et al., 2020;Ritter et al., 2013) (see Table 1 for species in this study with confirmed Wolbachia presence). Our results raise the intriguing possibility that such sweeps could play a role in the build-up of reproductive isolation (Giordano et al., 1997;Rokas, 2000;Shoemaker et al., 1999).  Habel et al., 2008;Lai & Pullin, 2004;Schmitt, 2007;Todisco et al., 2010;Zinetti et al., 2013), our divergence estimate for Leptidea reali and L. sinapis, the youngest and only pair for which divergence has been estimated based on genome-wide data before, is lower than previous estimates (Talla et al., 2017).

| Glacial cycling and the Messinian salinity crisis
Taking our estimates of species splits at face value, the divergence of eleven species pairs predates the onset of Pleistocene glacial cycling >2.6 MYA (Gibbard & Head, 2009 (Nürnberger et al., 2016) and reptiles (Kaliontzopoulou et al., 2011) and butterflies in the Melitaea radiation (Leveneu et al., 2009), this remains of course speculation. Roux et al., (2016) (Bunnefeld et al., 2018), and quantify the overlap between inter-and intraspecific divergence times. Nevertheless, our contrasting finding of both gene flow signals in old contact zone pairs (e.g.

| Do European butterfly species fall within the grey zone of speciation?
Pontia) and no evidence for gene flow (and complete sympatry) in the youngest pair (Leptidea) suggests that the 'grey zone of speciation' may be very wide indeed for European butterflies.

| Outlook
Given the challenges of demographic inference from transcriptome data (in particular the high relative recombination rate in butterflies), we have deliberately resisted the temptation to overinterpret models of demographic history. Our goal was instead to establish robust and comparable lower bounds for the age of butterfly sister species in Europe. Being based on mean divergence at 4D sites, these lower F I G U R E 4 Mitochondrial d xy (left) and d a (right) are weakly correlated with mean nuclear divergence (R 2 = 0.3565 and 0.2732, respectively). The coloured lines show the interactions for pairs that form contact zones and sympatric pairs. The two highlighted pairs (Iphiclides and Leptidea) have known Wolbachia-associated sweeps and low mt (and so high mtd a ) bounds for species ages make minimal assumptions and are unaffected by recombination. Likewise, we have decided to focus on a simple and conservative diagnostic for introgression.
Resolving these speciation processes in greater detail will require examination of whole-genome data from larger samples under realistic models of speciation history. Fitting explicit models of speciation to whole-genome sequence data, ideally including both selection and gene flow will undoubtedly refine estimates for the onset of divergence between these species pairs and overcome many of the biases inherent in basing such inferences on transcriptome data. Perhaps even more importantly, it would allow us to quantify the likely endpoints (if present) of speciation processes.
While it is straightforward to determine lower bounds for the onset of divergence under simple null models that assume no gene flow, as we have done here, estimating upper bounds of species divergence in the presence of gene flow is a much harder inference problem.
As pointed out by Barton and Hewitt (1985), the initial time of divergence may be unknowable given that post-divergence gene flow eventually erases all information about this parameter. Although current and historic levels of gene flow between European butterfly sister species remain to be determined, our results already suggest that their speciation histories are older and potentially slower than had been assumed by previous phylogeographic studies based on mt data. It will be fascinating to understand the evolutionary forces that drive both this general pattern and its exceptions, in particular, the selection responsible for the origin of very young but complete (in terms of reproductive isolation) cryptic species such as Leptidea (Talla et al., 2019) and the recently discovered Spialia rosae (Hernández-Roldán et al., 2016).

DATA AVA I L A B I L I T Y S TAT E M E N T
Read data are available from the ENA at PRJEB43082. Sequence alignments for the COI barcode locus were obtained from the dataset DS-EUGENMAP (dx.doi.org/10.5883/DS-EUGENMAP) on BOLD at www.bolds ystems.org and were originally produced by Dincӑ et al., (2021). The script used for calculating diversity and divergence is available at https://github.com/sameb don/ortho diver/ blob/maste r/ortho diver.py.