Evidence for spatial clines and mixed geographic modes of speciation for North American cherry‐infesting Rhagoletis (Diptera: Tephritidae) flies

Abstract An important criterion for understanding speciation is the geographic context of population divergence. Three major modes of allopatric, parapatric, and sympatric speciation define the extent of spatial overlap and gene flow between diverging populations. However, mixed modes of speciation are also possible, whereby populations experience periods of allopatry, parapatry, and/or sympatry at different times as they diverge. Here, we report clinal patterns of variation for 21 nuclear‐encoded microsatellites and a wing spot phenotype for cherry‐infesting Rhagoletis (Diptera: Tephritidae) across North America consistent with these flies having initially diverged in parapatry followed by a period of allopatric differentiation in the early Holocene. However, mitochondrial DNA (mtDNA) displays a different pattern; cherry flies at the ends of the clines in the eastern USA and Pacific Northwest share identical haplotypes, while centrally located populations in the southwestern USA and Mexico possess a different haplotype. We hypothesize that the mitochondrial difference could be due to lineage sorting but more likely reflects a selective sweep of a favorable mtDNA variant or the spread of an endosymbiont. The estimated divergence time for mtDNA suggests possible past allopatry, secondary contact, and subsequent isolation between USA and Mexican fly populations initiated before the Wisconsin glaciation. Thus, the current genetics of cherry flies may involve different mixed modes of divergence occurring in different portions of the fly's range. We discuss the need for additional DNA sequencing and quantification of prezygotic and postzygotic reproductive isolation to verify the multiple mixed‐mode hypothesis for cherry flies and draw parallels from other systems to assess the generality that speciation may commonly involve complex biogeographies of varying combinations of allopatric, parapatric, and sympatric divergence.

tions experience periods of allopatry, parapatry, and/or sympatry at different times as they diverge. Here, we report clinal patterns of variation for 21 nuclear-encoded microsatellites and a wing spot phenotype for cherry-infesting Rhagoletis (Diptera: Tephritidae) across North America consistent with these flies having initially diverged in parapatry followed by a period of allopatric differentiation in the early Holocene.
However, mitochondrial DNA (mtDNA) displays a different pattern; cherry flies at the ends of the clines in the eastern USA and Pacific Northwest share identical haplotypes, while centrally located populations in the southwestern USA and Mexico possess a different haplotype. We hypothesize that the mitochondrial difference could be due to lineage sorting but more likely reflects a selective sweep of a favorable mtDNA variant or the spread of an endosymbiont. The estimated divergence time for mtDNA suggests possible past allopatry, secondary contact, and subsequent isolation between USA and Mexican fly populations initiated before the Wisconsin glaciation.
Thus, the current genetics of cherry flies may involve different mixed modes of divergence occurring in different portions of the fly's range. We discuss the need for additional DNA sequencing and quantification of prezygotic and postzygotic reproductive isolation to verify the multiple mixed-mode hypothesis for cherry flies and draw parallels from other systems to assess the generality that speciation may commonly involve complex biogeographies of varying combinations of allopatric, parapatric, and sympatric divergence.

| INTRODUC TI ON
Biogeography is a key consideration for understanding speciation, as it affects the relationship between geographic isolation and gene flow under which reproductive isolation (RI) evolves and new biodiversity forms (Coyne & Orr, 2004). The geographic context of speciation has traditionally been divided into allopatric, parapatric, and sympatric modes, depending upon the spatial distribution of diverging populations (Bush, 1975;Coyne & Orr, 2004). Under classic allopatric or vicariant speciation ( Figure S1a), populations are physically separated (become bifurcated) into two large subpopulations by an uninhabitable area (Mayr, 1963). These isolates subsequently diverge in allopatry to form new species. Parapatric or clinal speciation ( Figure S1b) occurs when an ancestral population diverges across a spatial or ecological gradient (Endler, 1977). In this case, populations partly overlap with one another or can be arrayed as a series of stepping-stones partially connected by gene flow, as they evolve into different species. Finally, sympatric speciation (Figure S1c) occurs when diverging populations completely overlap, such that individuals of the evolving taxa remain in the cruising range of one another (Bush, 1975;Mayr, 1963).
The different modes of speciation are associated with different levels of gene flow among populations during the divergence process. Under allopatric speciation ( Figure S1a), gene flow between populations ceases due to a physical barrier of inhospitable habitat preventing migration (Mayr, 1963). The isolated populations are consequently free to independently accumulate genetic differences by natural selection and/or genetic drift to evolve RI, eventually leading to speciation. In parapatric speciation (Figure S1b), the partial overlap of populations results in a reduction but not elimination of migration between diverging demes. In this case, primary clines of changing gene and phenotype frequencies can evolve across the landscape generating a pattern of isolation by distance (IBD), which can potentially result in populations on different sides of the clines evolving into reproductively isolated species (Endler, 1977). Sympatric speciation ( Figure S1c) assumes that gene flow is possible throughout the divergence process. In this case, ecologically based divergent selection associated with biotic or abiotic environmental conditions or sexual selection is greater (stronger) than the level of gene flow, permitting populations to differentiate in the absence of geographic isolation (Bush, 1975;Mayr, 1963;Nosil, 2012). Distinguishing among these biogeographic possibilities is important because it can shed light on the factors and chronology whereby different types of RI evolve during speciation.
The classification of speciation into distinct allopatric, parapatric, and sympatric modes is to some degree one of convenience, distinguishing different general relationships between geographic isolation and gene flow during speciation. While it is possible for speciation to occur from start to finish solely in allopatry, parapatry, or sympatry, in reality, this may not always be the case. In other words, mixed or "hybrid" modes of speciation may be common in which the geographic context and levels of gene flow temporally vary during the divergence process (Coyne & Orr, 2004; 0.01 One example of a mixed mode of speciation involves secondary contact ( Figure S2a). In this case, a population initially becomes subdivided into isolates that diverge in allopatry. Through time, the physical barrier separating taxa dissipates and populations reestablish contact. If diverged populations have not evolved complete RI, then secondary genetic clines (a hybrid zone) can form, reflecting rates of migration, interbreeding, and selection against hybrids (Barton & Hewitt, 1981, 1985, 1989. Moreover, when a degree of postzygotic isolation has evolved in allopatry, a process of "reinforcement" can occur following secondary contact in which increased prezygotic RI is selected for to reduce the production of less fit hybrids (Servedio & Noor, 2003). Thus, reinforcement following secondary contact represents a mixed mode of speciation in that RI does not evolve solely in allopatry in the absence of gene flow. The same applies to other mechanisms in addition to reinforcement that can evolve following secondary contact to increase RI between populations in the face of gene flow (e.g., divergent ecological adaptation to different habitats). Such a mixed mode has been termed "allo-parapatric" or "allo-sympatric" by Coyne and Orr (2004). Populations in secondary contact may reflect a single cycle of allopatry, divergence, and subsequent overlap ( Figure S2a). However, populations can potentially undergo repeated cycles of allopatry and secondary contact during their history, with divergence accruing at a variety of different times when populations were isolated and overlapped (Hewitt, 2000).
Thus, there can also be various types of multiple mixed modes of speciation, such as "allo-para-allopatric," as the geographic distributions of diverging taxa change repeatedly over time ( Figure S2b).
Although there are many examples of secondary contact and hybrid zones in nature (see Barton & Hewitt, 1985;Coyne & Orr, 2004 for synopsis), some of which may involve the evolution of reinforcement (Servedio & Noor, 2003), there are other possible chronologies in which a mixed mode of speciation may unfold. One alternative mixed mode involves divergence beginning in parapatry, with populations initially forming a "primary" cline (Barton & Hewitt, 1981, 1989; Figure S2c). Over time, the populations on the ends of the clines may become geographically isolated, with intermediate populations potentially going extinct. In allopatry, the two populations can subsequently continue diverging and evolve additional RI. The above scenario is mixed because divergence did not accumulate only after the complete geographic isolation of populations but accrued to some degree in the face of gene flow in parapatry (or even sympatry) prior to allopatry. Such a mixed mode has been termed "para-allopatric" by Coyne and Orr (2004).
Envisioning different types of mixed or multiple mixed modes of speciation are not difficult and population divergence may often occur by such means in nature (Coyne & Orr, 2004). The problem is demonstrating their occurrence, as finding evidence supporting a mixed mode can be challenging, even for the seemingly straight-forward case of reinforcement following secondary contact (Barton & Hewitt, 1985;Butlin, 1989). In this regard, testing requires a study system having a known biogeography with features allowing inferences to be made that populations diverged under differing spatial contexts (e.g., parapatric and allopatric), experiencing varying levels of gene flow through time in the process. The problem is twofold.
First, one must characterize the biogeography of a pair of taxa to demonstrate that it has changed through time. Second, one must show that during times when taxa were both in contact and experiencing gene flow and geographically separated and isolated, they evolved RI.
Here, we focus on the first question of changing biogeography by examining two sister taxa of cherry-infesting fruit flies in the genus Rhagoletis (Diptera: Tephritidae), Rhagoletis cingulata, and Rhagoletis indifferens that possess attributes conducive for undergoing and testing for a mixed mode of divergence (Bush, 1966 Figure 1a). Thus, if cherry flies have or are speciating via a single mode of divergence, the appropriate null hypothesis would be an allopatric model reflecting their current geographic isolation and not parapatric or sympatric alternatives.
However, the geographic distribution of cherry flies is more complicated than a simple PNW versus ENA bifurcation or PNW, ENA, and SoM trifurcation. Flies described as R. cingulata and their black cherry host plants also exist in a semicircular arc of isolated populations found at higher elevations through the southwestern USA (SW) states of Texas, New Mexico, and Arizona that could have previously connected or been the source of the populations in the PNW and ENA (Figure 1a; Bush, 1966). In these isolated "sky islands" of the SW (Dodge, 1943;Heald, 1967) and northern Mexico (NoM), flies and their black cherry hosts persist in restricted pockets in certain mountain canyons experiencing higher rainfall and lower temperatures than surrounding low-lying deserts (Bush, 1966). Black cherries and R. cingulata also extend southerly from the SW through isolated mountains in the Chihuahuan Desert in NoM down into the Sierra Madre Oriental Mountains (SMOr) and SoM (Figure 1a; Bush, 1966;Rull, Aluja, & Feder, 2011).
The current distribution of R. cingulata and R. indifferens may still be explained by a null vicariant model of allopatric divergence (Figure 2a).
In this case, rather than being isolated by the northern plain states, R. cingulata and R. indifferens may instead have initially been separated between the SW and PNW into two large allopatric populations, with R. cingulata subsequently becoming subdivided across the SW, NoM, SoM, and ENA (Figure 2a). Such a modified vicariance hypothesis predicts that sharp genetic and phenotypic breaks should be observed associated with the geographic discontinuity between R. cingulata and R.
indifferens in the SW versus PNW, respectively (Figure 2a). A degree of divergence may also potentially be observed among R. cingulata populations across the SW, NoM, SoM, and ENA, depending upon how long ago they became isolated. However, the arc of cherry-infesting fly populations across the SW and NoM, as well as the existence of populations in SoM, also suggests the possibility of various alternative mixed modes of divergence as opposed to the null vicariance model.
For example, populations in the SW and NoM could represent recent secondary contact between previously isolated ENA, PNW, and SoM populations (i.e., allo-parapatry; Figure 2b). Alternatively, the current range of cherry flies could reflect a sequence of expansion from and fragmentation of what was once a continuous ancestral population across the SW and NoM (i.e., para-allopatry; Figure 2c). Finally, more complex mixed modes are also possible in which the biogeography involves past allopatry of ENA, PNW, and SoM populations, followed by secondary contact across the SW and NoM, and then recent isolation In the current study, we characterize cherry fly populations across North America for variation in nuclear-encoded microsatellite loci, maternally inherited mitochondrial DNA (mtDNA), and the wing spot morphological trait. Our primary objective is to distinguish between vicariance and a mixed mode of divergence by determining whether clinal variation and a pattern of IDB exist for current sky island isolates of cherry flies across the SW and NoM connecting R. cingulata in ENA with R. indifferens in the PNW, as well as potentially with R. cingulata in SoM. A secondary objective, given the existence of clines, is to evaluate whether the genetic data provide any insight into which of the various mixed-mode scenarios may be occurring. In this regard, much is known about the history of climate change and its effects on the biogeography of the flora and fauna of the SW and NoM (Dansgaard et al., 1993;Martin & Harrell, 1957;McDonald, 1993). Specifically, this region is known to have undergone a series of dramatic climatic oscillations during the Pleistocene driven by major ice ages occurring at ~100,000 year-intervals (Hewitt, 1996). These ice ages were punctuated by interglacial periods, the latest one beginning ~11,700 years ago (ya), demarcating the start of the Holocene (Hewitt, 1996). During the cooler and wetter ice ages, forests and woodlands, and their associated flora and fauna, were more contiguous across the SW and NoM (Betancourt, Van Devender, & Martin, 1990), fostering migration and gene flow. In contrast, during the intervening interglacial periods, warmer and drier conditions drove the forests and woodlands to higher elevation sky islands separated by low-lying deserts, restricting gene flow.
The information on climate change can be used in conjunction with the current distribution of cherry flies to generate predictions concerning expected divergence times of populations under the differing mixed-mode hypotheses. As cherries, and hence cherry flies, would be associated with the changing distributions of forest and woodland habitats, ice ages during the Pleistocene would represent times when fly populations were more contiguous through the SW and NoM. Thus, during glacial periods, gene flow would be possible, fostering the potential formation of clines. In contrast, interglacial periods, when cherries and cherry flies were restricted to isolated pockets in sky islands, as is currently the case, would correspond to times of limited or no gene flow. Consequently, during interglacial periods, populations could diverge in allopatry. Given that the most recent warming and drying started 8,000 to 9,000 ya in the early Holocene (Allen, Betancourt, & Swetnam, 1998), if the allo-parapatry hypothesis is true, then this would imply recent allopatric separation of ENA, PNW, and F I G U R E 2 Alternative geographic divergence hypotheses for cherry-infesting Rhagoletis in North America. (a) Vicariance model; a large, ancestral cherry fly population harboring little to no spatial variation bifurcated into R. cingulata and R. indifferens along a spatial divide between the southwestern USA (SW) and Pacific Northwest (PNW). Subsequently, R. cingulata became recently subdivided into Southwestern, northern Mexican (NoM), southern Mexican (SoM), and eastern USA (ENA) subpopulations. Note that the classic allopatric/ vicariant model predicts sharp genetic and phenotypic breaks should currently be seen between cherry flies in the PNW and elsewhere, as indicated by the orange compared to green color, respectively, of populations in these regions; (b) Allo-parapatric model; a large, ancestral population harboring little to no spatial variation separated into ENA, PNW, and SoM isolates. These isolates subsequently diverge in allopatry. Most recently in the Holocene, flies were reintroduced into the SW and NoM, forming secondary clines; (c) Para-allopatric model; a primary cline of isolation by distance forms through the ancestral cherry fly population across the SW and NoM sometime during the Wisconsin glaciation. As ice sheets retreated following the last glacial maximum, cherry flies migrated into ENA and the PNW. Subsequent warming and drying in the early Holocene resulted in these populations and those in the SW and NoM becoming geographically isolated and differentiating; (d) Allo-paraallopatric model; an ancestral cherry fly population becomes geographically isolated and diverges into ENA, PNW, and SoM subpopulations during the interglacial period preceding the Wisconsin glaciation. Populations then reestablish contact during the Wisconsin glaciation, creating secondary clines through the SW and NoM. In addition, introgression from SoM introduced the diverged mtDNA haplotype II into NoM and the SW. Subsequent warming and drying in the early Holocene resulted in the cherry fly populations becoming geographically separated into PNW, SW, NoM, SoM, and ENA and differentiating, with the mtDNA haplotype II not reaching the PNW and ENA prior to regional subdivision. A multiple mixed model involves a combination of hypotheses c and d, with primary clines forming across the SW and NoM being augmented by allopatry, secondary contact, and gene flow from SoM, especially for mtDNA haplotype II

| Study sites and specimens
All flies were directly collected as larvae in infested cherry fruit from the field and reared to adulthood in the laboratory using standard Rhagoletis husbandry methods (Feder, Chilcote, & Bush, 1989) before being frozen at −80°C for later genetic analysis. Only one egg is typically deposited per cherry, with a maximum of one larva emerging per each infested fruit (Bush, 1966). Thus, there is not an issue with having to choose to genotype one individual per cherry to avoid possible problems in biasing estimates of gene or haplotype frequencies due to scoring of related individuals.
Flies from 17 populations sampled between 2004 and 2011 were genotyped for microsatellites and/or mtDNA ( Figure 1a; Table S1).
Previous surveys of R. cingulata in Mexico indicated that while the host black cherry is contiguous in the SMOr, the fly is now found only in two northern and southern areas of this mountain range (Rull et al., 2011), both sampled here ( Figure 1a). Similarly, while P. serotina is present in the Sierra Madre Occidental Mountains of Mexico, flies have not been found in the area (Rull et al., 2011).
Bitter and black cherries co-occur in areas of the states of Arizona and New Mexico ( Figure 1a). However, only black cherries were sampled from the SW, as bitter cherry is rare. The presence of both black and bitter cherries in the SW suggests that these two host plants and their associated fly populations cooccurred in the region in the past. Rhagoletis cingulata and R. indifferens are known to use other cherry species in common as secondary hosts (Bush, 1966). Thus, host affiliation is not likely to be a strong barrier to gene flow. In the absence of nonhost related pre-or postzygotic RI, if these two taxa overlapped in the past, or were to cooccur in the future, then they could potentially hybridize, thereby generating clines. In this regard, mating trials between R. cingulata from ENA and R. indifferens from the PNW showed moderate prezygotic isolation (isolation index = 0.27 = the percentage of intraspecific subtract interspecific matings in trials divided by the total percentage of matings; Hood, Egan, & Feder, 2012). In addition, Doellman et al. (2019) reported no evidence for strong postzygotic isolation in crosses between PNW and ENA cherry flies, further supporting the likelihood that R. cingulata from ENA and R. indifferens from the PNW would form a hybrid zone and even possibly fuse if they were to come back into contact.

| Multilocus nuclear data
A total of 364 adult flies from populations 1-15 ( Figure 1a; Table S1) were genotyped for 21 microsatellite loci (see Supplemental Methods for details). References for primers and conditions used for PCR amplification can be found in Maxwell, Thistlewood, and Keyghobadi (2013) and Michel et al. (2010). Due to their relatively high rates of mutation and levels of polymorphism (Estoup & Angers, 1998) We first constructed a genetic distance network useful to graphically summarize the genetic relationships among cherry fly populations across North America and visualize possible patterns of IBD using PowerMarker v3.25 (Liu & Mus, 2005). This was followed by Mantel tests of IBD using vegan (R package version 2.2-1; Oksanen et al., 2015). We also quantified levels of genetic divergence within and between regions by calculating Jost's D values (D Jost ) with the R package diverSity, as suggested for highly polymorphic microsatellite loci (Jost, 2008;Keenan, McGinnity, Cross, Crozier, & Prodöhl, 2013). Tests for evidence of genetic subdivision among fly populations were next performed using Structure v2.3.4 (Evanno, Regnaut, & Goudet, 2005;Pritchard, Stephens, & Donnelly, 2000).
Lastly, we estimated divergence times between populations based on the Metropolis-coupled Markov Chain Monte Carlo (MCMC) sampling algorithm in iMa2P (Sethuraman & Hey, 2016). We focus on reporting results assuming an average microsatellite mutation rate of 6.3 × 10 -6 per meiosis, which generated divergence time estimates for cherry flies corresponding well to known climatic changes in Holocene (see Sections 3 and 4).

| Mitochondrial DNA
Mitochondrial DNA sequence data were generated for a total of 97 flies for a 1,482 bp fragment that included parts of the cytochrome oxidase I (COI), leucine transfer RNA (tRNA-Leu), and cytochrome oxidase II (COII) genes (Table S1; GenBank accession numbers: KT221476-KT22179), using the overlapping primer pairs UEA5 and UEA10 (Lunt, Zhang, Szymura, & Hewitt, 1996) and C2-J-3138 and TK-N-3782 (Simon et al., 1994). Due to its lack of recombination and maternal inheritance, we anticipated that analysis of mtDNA haplotypes could provide useful information to detect and date past periods of geographic isolation for cherry flies. Methods for Sanger sequencing of mtDNA are described in Feder et al. (2003). Sequences were edited manually and assembled using codoncode aligner v3.7 (Codon Code Corp). A haplotype network was constructed for the mtDNA data using PegaS (R package version 0.11; Paradis, 2010).

| Clinal patterns in multilocus nuclear data
The genetic distance network for microsatellites shows a pattern of IBD in extant cherry fly populations across North America (Figure 1b).
Populations of R. cingulata and R. indifferens were arrayed in the network in relation to their geographic location from ENA to Mexico (north to south), to the SW, to the PNW. As a result, a significant correlation was observed between geographic and overall genetic distance (IBD) for all 21 microsatellite loci across the 15 populations genotyped in North America (r = .81; p < .0001, as determined by nonparametric Monte Carlo simulations; Table 1). Significant genetic IBD was not only observed for the 21 microsatellites analyzed together but also for 19 of the 21 loci (90.5%) considered individually (Table 2). Although microsatellite genetic distance between populations was strongly correlated with geographic distance overall, the pattern was largely driven by sites in the SW and Mexico (Table 1; Figure 1b). Thus, among the six SW and Mexican populations of R. cingulata, significant IBD was observed (r = .86; p < .0001; Table 1).
However, significant IBD was not detected among the six R. indifferens populations in the PNW (r = .48, p = .071) or among the three R. cingulata populations in ENA (r = −.22, p = .50).
Estimates of D Jost within and among regions displayed patterns of differentiation matching those seen in the IBD analysis (see Table S2 for  (Table S3). The situation differed slightly for the SoM sites, as two microsatellites (WCFF24 and WCFF105) possessed unique alleles in these two populations at combined high frequencies not found elsewhere (Table S3). In addition, the two SoM populations had more alleles unique to the region not found at any other site (n = 45) than both the three ENA populations (n = 25) and six PNW populations (n = 21) surveyed (Table S3).

| Population structure from multilocus nuclear data
Structure analysis implied that the clines were due to fly populations differing in the frequencies of alleles they possess rather than representing differences in the proportions of parental R. indifferens and R. cingulata genotypes across sites (Figure 3c). The largest change in TA B L E 1 Correlation coefficients (r) and significance levels (p-value) for Mantel tests of Isolation by distance (IBD) for the combined 21 microsatellite dataset scored across cherry-infesting fly populations spanning indicated geographic regions (North America = sites 1-15; Pacific Northwest = sites 1-6; southwest USA & Mexico = sites 7-12; eastern USA = sites 13-15) log likelihood in the Structure analysis was observed around K = 2 (Table 3)

| Divergence estimates from multilocus nuclear data
The estimated divergence time with the highest posterior probability between R. indifferens in the PNW (site 4) and R. cingulata in Arizona (site 7), assuming a mutation rate of 6.3 × 10 -6 , was 5,556 in the iMa2P analysis. The results above were based on a rate of 6.3 × 10 -6 that corresponds to that of Drosophila (Schug, Mackay, & Aquadro, 1997). However, using higher mutation rates in the analysis resulted in proportionately more recent estimates of population divergence ( Figure S3).

| Mitochondrial DNA divergence
Mitochondrial DNA haplotypes displayed a different pattern of geographic variation compared to the microsatellites (Figure 4) (Brower, 1994).

| Clinal patterns in morphological variation
The

| Null vicariant versus mixed-mode divergence
With respect to our primary objective, the results for the microsatellites reveal a clinal pattern of allele frequency differences in cherry fly populations across the PNW, SW, ENA, and NoM (Figures 1b and 3b,c), as predicted for a mixed mode of divergence.
Populations from SoM also fit the general IBD pattern for microsatellites ( Figure 1b). However, there appears to be a geographic  Table S3) and vice versa (see orange boxes in Table S3). The pat-

| Assessing different mixed-mode models
Given the finding of clinal variation, our secondary objective was to assess whether patterns of genetic differentiation and estimated divergence times for populations, when coupled with the known climatic history of the SW and NoM, might provide insight into the likelihoods of different mixed modes of speciation for cherry flies.
As discussed in the introduction, there are three major possibilities

| Hypothesis one: allo-parapatric model
Hypothesis one contends that the microsatellite clines observed across the SW and NoM were formed by recent secondary contact during the Holocene of isolated fly populations from the PNW, ENA, and SoM that diverged previously in allopatry (Figure 2b). We consider this scenario unlikely for several reasons. First, given the current physical separation of cherry flies on sky islands across the SW and NoM, there is little likelihood of migration and gene flow among these populations at the present time. Second, the current isolation of cherry flies corresponds climatically and dates genetically to the Holocene, when warmer and drier conditions starting 8,000 to 9,000 ya forced forest and woodland associated plants and animals to higher elevations (Allen et al., 1998). We assumed a microsatellite mutation rate of 6.3 × 10 -6 for our estimates of divergence time for cherry flies, which corresponds to that of Drosophila (Schug et al., 1997) and represents the lower end of rates for insects (Estoup & Angers, 1998). Thus, divergence time is not likely to be much greater than those we report, as this would require unrealistically low mutations rate. However, it is possible that actual mutation rates in cherry flies are higher than 6.3 × 10 -6 , which would proportionately reduce estimated divergence times. For example, a mutation rate of 1 × 10 -3 would change the estimated divergence time between R.
indifferens in the PNW and R. cingulata in Arizona from 5,556 to 50 ya ( Figure S3). Such a timeline would seem to make the hypothesis of recent secondary contact more plausible. However, such recent divergence times would require that cherry flies were independently introduced into the SW and NoM from allopatric PNW, ENA, and SoM populations within historical times, and there is no record of this (Bush, 1966). Also, following these introductions, cherry flies from the different regions would have subsequently had to spread rapidly in an island-hopping fashion among mountain isolates to generate the observed microsatellite clines, which seems highly unlikely. Third, and finally, a scenario of recent introductions would have difficulty in accounting for the lack of mtDNA variation throughout the SW and Mexico. Given recent gene flow from the PNW, ENA, and SoM posited by hypothesis one, we would expect to observe polymorphism for the two major mtDNA haplotypes I and II across the SW and NoM, which is not the case, as only haplotype II is present.

| Hypothesis two: para-allopatric model
Setting aside consideration of the diverged mtDNA haplotype II in the SW and Mexico, hypothesis two (para-allopatric model) appears to best explain the pattern of microsatellite and wing spot variation observed across North America, particularly with respect to the PNW and ENA (Figure 2c) share the same mtDNA haplotype I in common and there is limited mtDNA variation within either region (Figure 4), which again implies a recent origin and lack of previous allopatry between flies in these two regions. It may seem odd that mtDNA haplotype I is identical between the ENA and the PNW given their estimated di-  99949 1,482 ). There is a major difficulty with the para-allopatric hypothesis, however, which concerns the presence of the diverged mtDNA haplotype II in the SW and Mexico.
As alluded to above and discussed below, we hypothesize that the mtDNA discrepancy may reflect a cycle of allopatry and secondary contact involving flies from SoM dating to the interglacial period 100,000 to 150,000 ya preceding the Wisconsin glaciation.

| Hypothesis three: allo-para-allopatric model
We argue that hypothesis three (Figure 2d) involving an allo-paraallopatric sequence of events is not the most parsimonious explanation for the pattern of differentiation in the PNW and ENA but that a variation of the scenario may help account for the diverged mtDNA haplotype II observed through the SW and Mexico. With respect to the PNW and ENA, there is no genetic evidence for these populations having been allopatrically isolated during the interglacial period 100,000-150,000 ya prior to the Wisconsin glaciation. The estimated divergence time between PNW and ENA based on the microsatellites is 15,079 ya (95% credible interval 7,143 to 31,270 years), which would place any past allopatry, if it occurred, in the Wisconsin glaciation.
One possible scenario that could account for allopatry during the Wisconsin glaciation is that populations in the PNW and ENA resided in one or more northern refugia during this time. Thus, instead of expanding northward following retreating ice sheets, current PNW and ENA populations migrated southward from refugia as conditions warmed to contact one another and form clines across the SW and NoM by the early Holocene. While seemingly plausible, the refugium/refugia scenario of secondary contact has several problems.
For example, if two refugia were involved (possible candidates being the west coast of Alaska/British Columbia and off the east coast of New England), then it would be difficult to explain how current populations in the PNW and ENA that emanated from these refugia share the same mtDNA haplotype I in common. A potential solution to this issue may be that PNW and ENA were actually connected at some time by a conduit through the North, permitting gene flow and, thus, the homogenization of mtDNA haplotypes between the two regions. In this regard, the pin cherry, Prunus pensylvanica, currently has a distribution through Canada connecting bitter cherry, the main host of R. indifferens in the PNW, with black cherry, the major host of R. cingulata in ENA ( Figure S4). Maxwell et al. (2013) have reported rearing R. indifferens from pin cherry fruit collected in British Columbia, Canada. Pin cherry has also been cited as an alternate, although rare, host for R. cingulata in ENA (Bush, 1966).
If true, then cherry flies may actually have formed a ring around North America at some time in the past (Irwin, Irwin, & Price, 2001).
However, pin cherry is also the primary native host for R. fausta, a fly that outcompetes R. cingulata on P. pensylvanica (Bush, 1966).
Moreover, our attempts at rearing R. cingulata from pin cherries in the states of Wisconsin and Minnesota, USA have to date been unsuccessful. Thus, there are reasons to doubt that pin cherry is or was an effective bridge between PNW and ENA populations of cherry flies now or in the past.
A single refugium scenario also has difficulties. For example, given that it does not appear likely for pin cherry to have served as a conduit for east/west movement of flies, some alternative and not readily apparent means would need to be discerned for how migrants spread from a hypothesized single refugium to subsequently come into secondary contact in the SW and NoM. In addition, while the existence of a single refugium could explain the lack of mtDNA haplotype I divergence between the PNW and ENA, it would also predict that a genetic signature should be observed for microsatellites reflecting the period of isolation for these flies in the refugium. However, this is not the case for the microsatellites, as there are no alleles of consequential frequency shared uniquely among all fly populations in the PNW and ENA and absent from the SW and Mexico (Table S3).
While the current data do not provide evidence for glacial refugia/refugium in the North and are most consistent with a primary origin for the microsatellite clines connecting PNW and ENA populations, this may not be the case for the diverged mtDNA haplotype II found in the SW, NoM, and SoM. In this regard, the estimated divergence time between mtDNA haplotypes I and II is 100,000-150,000 ya, depending upon the rate considered for the molecular clock.
This timeline corresponds with the interglacial period preceding the Wisconsin glaciation when, similar to the current time, cherry fly populations in SoM may have been isolated from those elsewhere.
Thus, it is conceivable that the diverged mtDNA haplotype II may reflect a period of past allopatry for cherry flies. In accord with this scenario, flies from SoM do possess a greater number of microsatellite alleles unique to the region (n = 45) than are found in the PNW (n = 21) and ENA, despite fewer numbers of sites and flies being genotyped from SoM (see blue boxes in Table S3). Indeed, for the five loci WCFF24, WCFF93, WCFF105, WCFF11, and P27, the alleles unique to SoM populations are either individually or collectively present at appreciable frequencies >0.10. Thus, there is better evidence for a nuclear genetic signature of past allopatry for the SoM compared with PNW and ENA populations.
The allo-parapatric scenario for SoM flies, if true, still must explain certain anomalies in the data. Most importantly, if secondary contact and gene flow occurred between SoM, SW, and NoM flies, then why is mtDNA haplotype II fixed, while microsatellites display IBD among these populations (Tables 1 and 2)? Several factors alone or in combination may potentially account for differences in the apparent introgression dynamics between mtDNA and microsatellites, but still require testing to verify. These include differential lineage sorting related to the smaller effective population size for mtDNA (Charlesworth, Bartolome, & Noel, 2005), sex-related differences in migration rates, the selective sweep of a favorable mtDNA variant (Galtier, Nabholz, Giémin, & Hurst, 2009;Silva, Lima, Martel, & Castilho, 2014), or the presence of a cytoplasmically inherited endosymbiont (e.g., Wolbachia; Brucker & Bordenstein, 2012;Jaenike, Dyer, Cornish, & Minhas, 2006;Werren, 1998). We return to the possibility of an endosymbiont below when we discuss what is currently known concerning patterns of RI among cherry flies across North America.
In summary, it appears that a combination of hypotheses two and three may best explain the history of cherry flies in North America based on the current genetic and phenotypic data. Thus, while the presence of clines was confirmed in the present study, the disjunct pattern for mtDNA complicates the story. Additional DNA sequencing studies could help to further resolve and complement the microsatellite and mtDNA results. In particular, sequences for nuclear-encoded genomic regions experiencing low rates of recombination may confirm or discount past allopatry.
In other Rhagoletis flies showing evidence of mixed geographic modes of speciation (Xie et al., 2007), sequencing of inversion polymorphisms has proven useful for inferring secondary contact when rearrangements mapping to different chromosomes were shown to have similar divergence times (Feder et al., 2003(Feder et al., , 2005.

| Evolution of RI
As outlined in the Introduction, the problem of verifying a mixed mode of speciation is two-fold. First, one must characterize the biogeography of a pair of taxa to demonstrate that it has changed through time. Second, one must show that during times when taxa were both in contact and experiencing gene flow, and geographically separated and isolated, they evolved RI. We have accomplished the first goal in the current study. The genetic and phenotypic patterns of clinal variation found for cherry flies imply that they have a biogeography consistent with a mixed mode of speciation and our results shed light on which of the various types of mixed-mode models appear most likely to be responsible for generating the observed patterns. However, our findings do not directly address the second issue of the evolution of RI. Thus, we may be able to infer that genetic differences likely accrued between cherry flies during periods of parapatry and allopatry but this does not demonstrate that RI also evolved at these times. shown for other insects (Brucker & Bordenstein, 2012;Jaenike et al., 2006;Werren, 1998). While speculative, Schuler et al. (2013) reported the presence of Wolbachia in native populations of R. cingulata in the USA and in introduced populations of the species in Europe. In addition, different strains of Wolbachia have been associated with specific mtDNA haplotypes and the occurrence of cytoplasmic incompatibility between different races of the European cherry fruit fly, R. cerasi (Arthofer et al., 2009;Riegler & Stauffer, 2002). Thus, there is precedent for Wolbachia sweeping through Rhagoletis populations (Bakovic, Schebeck, Telschow, Stauffer, & Schuler, 2018;Schuler et al., 2013Schuler et al., , 2016. With regard to the results of Tadeo et al. (2015), if Wolbachia is involved, then we would predict that at least two strains of the endosymbiont differentiate Mexican from ENA cherry flies in order to account for the apparent bidirectionality of postzygotic RI.

| CON CLUS ION
While potentially common, our study highlights how challenging it can be to assess the mixed-mode hypothesis. Thus, aside from reinforcement evolving following secondary contact, evidence for mixed-mode speciation is generally lacking. Our findings for R. indifferens and R. cingulata are consistent with a mixed mode of divergence for these flies. However, further DNA sequencing studies and analyses of pre-and postzygotic isolation among cherry fly populations are needed to strengthen the biogeographic implications of the current microsatellite and mtDNA data and to show that in addition to population divergence, RI also accrued during periods of parapatry and allopatry.
Moreover, our findings suggest that different mixed modes of divergence may apply to different parts of the cherry fly range, further complicating the problem. Studies of other systems are also clearly needed to extend the results for cherry flies more generally. However, many species may not possess attributes ideal for detecting a mixed mode of divergence despite its occurrence. For example, if climate change continues and causes the extinction of cherry fly populations through the SW and NoM in the near future, then things would change. Instead of detecting clinal variation and inferring a mixed mode of divergence, a future study would likely conclude that PNW, ENA, and SoM populations represent allopatric isolates undergoing vicariant speciation. Thus, the dissimilarities between cherry flies and other systems may reflect the particular snapshot in time that they happen to capture rather than any inherent differences. In this regard, suture zones in which populations have been hypothesized to have undergone repeated cycles of allopatry, secondary contact, and gene flow associated with glaciation potentially represent rich systems to extend the results for cherry flies more distantly in time. Studies of more recently diverged species may help reveal if signatures of parapatric divergence can be found prior to allopatry. Regardless, our results also underscore the importance of analyzing different types of genetic markers and genomes with differing demographic and evolutionary dynamics to accurately and fully resolve the biogeography of a system and its consequences for speciation.
Moreover, the possibility of cytoplasmic sweeps, perhaps due to endosymbionts, urges caution about inferring biogeographic history based on the assumption that cytoplasmically inherited genomes like mtDNA represent passive, neutral markers of population structure. Thus, while categorizing speciation into different modes and assuming neutral gene flow dynamics be can be helpful for understanding population divergence, life may often produce additional complications.

CO N FLI C T O F I NTE R E S T
All authors declare no conflict of interest.