The role of divergent host use and geography in the evolution of habitat isolation and sexual isolation among sister species of Belonocnema gall wasps Free

Abstract

Ecology and geography can play important roles in the evolution of reproductive isolation across the speciation continuum, but few studies address both at the later stages of speciation. This notable gap in knowledge arises from the fact that traditional ecological speciation studies have predominantly focused on the role of ecology in initiating the speciation process, while many studies exploring the effect of geography (e.g., reinforcement) concentrate on species pairs that lack divergent ecological characteristics. We simultaneously examine the strength of habitat isolation and sexual isolation among three closely related species of Belonocnema gall-forming wasps on two species of live oaks, Quercus virginiana and Q. geminata, that experience divergent selection from their host plants and variable rates of migration due to their geographic context. We find that the strength of both habitat isolation and sexual isolation is lowest among allopatric species pairs with the same host plant association, followed by allopatric species with different host plant associations, and highest between sympatric species with different host-plant associations. This pattern suggests that divergent selection due to different host use interacts with geography in the evolution of habitat isolation and sexual isolation during the later stages of speciation of Belonocnema wasps.

Open in new tabDownload slide

Introduction

A major goal of speciation research is to understand the evolutionary forces that are responsible for generating reproductive barriers between diverging lineages (Coyne and Orr, 2004; Nosil, 2012; Wright, 1978). In this pursuit, many studies have documented an important role for divergent selection associated with differences in ecology in promoting reproductive isolation (RI), a process termed ecological speciation (Rundle and Nosil, 2005; Schluter, 2000). However, ecology may also interact with other factors, including geography, to initiate, promote, or complete the speciation process by influencing the opportunity for gene flow (Doellman et al., 2020; Nosil, 2007; 2012; Schwartz et al., 2010). When lineages are in geographic contact where gene flow is greater than zero (e.g., sympatry, parapatry, or peripatry), selection can favor increased prezygotic isolation (Butlin and Smadja, 2018; Nosil, 2007) if there is a fitness cost to hybridization (e.g., hybrid inviability) through a process known as reinforcement (Servedio and Noor, 2003). In addition, direct selection against immigrants from alternative environments can also promote increased habitat isolation by favoring individuals with high fidelity for their native environment (Nosil et al., 2006). These two mechanisms can occur in lineages experiencing migration and/or gene flow, and under these conditions, lineages are predicted to exhibit higher prezygotic isolation than allopatric lineages with little or no gene flow (Brown and Wilson, 1956; Butlin, 1987; Servedio and Noor, 2003). Alternatively, gene flow in sympatry can lead to the erosion of RI by homogenizing differences between lineages (Nosil et al., 2003; Slatkin, 1987).

RI between lineages with divergent ecology, despite gene flow due to geographical contact, has been observed in numerous systems, providing evidence for the role of ecology in promoting speciation (Craig et al., 1993; Funk, 1998; Lowry et al., 2008; Nosil, 2012; Nosil et al., 2002; Rundle et al., 2000; Via et al., 2000). However, the majority of studies have primarily focused on the early stages of speciation (as reviewed by Nosil, 2012; Rundle and Nosil, 2005). How divergent ecology interacts with geography to affect RI during the later stages of speciation remains surprisingly unclear (Nosil, 2012). This knowledge gap arises because traditional ecological speciation studies have mainly emphasized how ecology initiates the speciation process, whereas reinforcement studies have primarily examined species pairs lacking divergent ecology (reviewed by Servedio and Noor, 2003). One way to assess the role of ecology (i.e., divergent natural selection) and geography (i.e., the likelihood of gene flow) during the later stages of speciation is to compare the strength of RI between lineages characterized by different combinations of ecology and geography (Nosil, 2007).

Host-plant-associated insect herbivores offer a powerful experimental system to address these questions, as their ecology and life history are commonly restricted to the discrete space of the host plant where they feed and develop to adulthood, find mates, and lay eggs (Funk et al., 2002). By comparing different allopatric and sympatric populations living on both similar and different host plants, one can determine the individual and combined roles of divergent host plant use (ecology) and geographic overlap (gene flow) in the evolution of RI and speciation (e.g., Egan, Hood, Feder, et al., 2012; Funk, 1998; Nosil et al., 2002; Rosenblum & Harmon, 2011). Despite a handful of studies addressing this topic, the outcomes are still unclear in many natural systems.

We apply this comparative approach to determine the roles of divergent host use and geography by comparing the strength of both habitat isolation and sexual isolation among a complex of three host-plant-specific species of gall-forming wasps in the genus Belonocnema (Hymenoptera: Cynipidae: Cynipini). These three species are each locally adapted to different live oak species in the genus Quercus in the subsection Virentes (Zhang, Hood, Carroo, et al. 2021; Zhang, Hood, Roush, et al., 2021). Two species, B. treatae and B. fossoria, co-occur over a large portion of their geographic ranges in the coastal southeastern United States on their respective host plants, the southern live oak, Q. virginiana (Qv), and the sand live oak, Q. geminata (Qg) (Driscoe et al., 2019; Zhang, Egan, et al., 2021). The third species, B. kinseyi, forms galls on Qv in the western-most portion of its range along the Gulf coast where it is allopatric to both B. treatae and B. fossoria (Figure 1). The two sympatric species B. treatae and B. fossoria are in the later stages of the speciation process based on the following phenotypic and genomic evidence. First, this sympatric species pair exhibits near complete prezygotic RI (range among populations = 0.95–0.99) due to multiple reproductive barriers, including temporal isolation (Hood et al., 2019), sexual isolation (Egan, Hood, Feder, et al., 2012), and habitat isolation (Egan, Hood, & Ott, 2012). Second, the genomic composition of the two species is distinctive, with no evidence of contemporary admixture; however, migration still occurs between the sympatric species at low levels (Driscoe et al., 2019). Thus, both the geographic distribution and patterns of host plant use by the three Belonocnema lineages allow us to simultaneously test the mechanisms underlying the evolution of two prezygotic reproductive barriers: sexual isolation and habitat isolation.

Herein, we assess the strength of each barrier between the following species pairs of Belonocnema to evaluate patterns of RI as a function of differing environments and geographic contexts (Funk et al., 2002; Nosil, 2007): (a) “allopatric, same-host”—B. kinseyi × B. treatae (geographically separate species pair on the same host plant), (b) “allopatric, different-host”—B. kinseyi × B. fossoria (geographically separate species pair on different host plants), and (c) “sympatric, different-host”—B. treatae × B. fossoria (geographically overlapping species pair on different host plants). Using these three types of comparisons, we test whether allopatric lineages on different host plants will exhibit higher RI than allopatric lineages on the same host plant due to the role of divergent selection in promoting speciation and test whether sympatric lineages on different host plants will exhibit higher RI than allopatric lineages on different host plants due to the role of selection favoring reduced migration and reduced hybridization in sympatry.

Methods

Study system

In the spring, sexual generation adult males and females in the genus Belonocnema, which develop separately in single-sex, multichambered root galls, emerge and mate. (Lund et al., 1998). Upon mating, females immediately oviposit into newly flushed leaves. This study focuses on estimating RI among sexual generation individuals, where habitat isolation and sexual isolation can directly affect gene flow among the three Belonocnema species (Egan, Hood, Feder, et al., 2012, Egan, Hood, & Ott, 2012). To isolate the effect of host preference on mating probability, we measured the degree of sexual isolation in the absence of host plant as shown by Egan, Hood, and Ott (2012).

In addition to the overlap of geographic distributions and host associations described above, additional evidence highlights the suitability of the Belonocnema-live oak study system to isolate the role of ecology and geography in the evolution of RI. Despite general geographic sympatry, the two different host plant species occupy different microhabitats, where Qg occurs in drier, higher pH, sandy-like soil compared with Qv that occurs in moist, lower pH soil. Despite this difference, two host plant species can co-occur only meters apart (L. Zhang, personal observation), creating opportunities for Belonocnema wasps to encounter different host plants. In addition, Qv and Qg also differ in the phenology of spring leaf flush, physiological traits, tree height, and leaf morphology (Cavender‐Bares and Pahlich, 2009; Cavender‐Bares et al., 2015). As a result, Qv and Qg exert strong divergent selection on the Belonocnema species, as supported by the observation that both immigrant and hybrid gall-forming wasps suffer significant reductions in fitness dependent on the host plant species (Zhang et al., 2017; Zhang, Hood, Carroo, et al., 2021; Zhang, Hood, Roush, et al., 2021; Supplementary Figure S1).

Sample collection

In the spring of 2016–2019, we collected mature root galls containing penultimate stage pupae of each Belonocnema species from their respective host plants: B. fossoria (three sites), B. treatae (six sites), and B. kinseyi (three sites) (Figure 1; Supplementary Table S1). To secure unmated females, individual root galls were each placed in a 30-ml clear-plastic vial and stored at room temperature. The timing of peak emergence of B. treatae occurs 2–3 weeks earlier than B. fossoria (Hood et al., 2019). Therefore, to synchronize emergence for experimental purposes, we placed one-half of the root galls collected from earlier emerging Qv at each site into a refrigerator at 4 °C for 1 week, which effectively delayed emergence of B. treatae. To minimize age variation among adults used for host preference and mate choice, all root galls were monitored every 2 days. Upon emergence, the 1- to 2-day-old adults were collected for use in the experiments described below. All tests of RI were conducted in a greenhouse under ambient environmental conditions at Rice University, Houston, TX, in 2016, 2018, and 2019, or in a screened-in outdoor enclosure at Archbold Biological Station, Venus, FL, in 2017.

Estimating habitat isolation

To evaluate the potential role of differences in ecology (host plant use) and differences in geography (opportunities for migration and/or gene flow) on habitat isolation, we measured host plant preference for the two sympatric species, B. treatae and B. fosoria, and the allopatric species, B. kinseyi, using a controlled two-choice experimental approach following Egan, Hood, and Ott (2012) (Supplementary Table S1). As female host choice largely determines the environment in which mating, oviposition, and larval development proceed the following generation, and male Belonocnema display lower host fidelity (Egan, Hood, & Ott, 2012), we used females to quantify host preference. In 2016 and 2017, host preference trials were conducted within sealed 500-ml clear-plastic cups that contained a fresh 15- to 20-cm stem cutting of both Qv and Qg from Florida with a similar number and size of leaves to minimize potential effects of differences in plant biomass on host choice. In 2019, trials were conducted within 60 × 15 mm Petri dishes stocked with a single newly flushed leaf of equivalent size of each host plant positioned at opposing sides. For all trials in each year, a single unmated female was aspirated into each container and then observed at 2-min intervals for 30 min for a total of 15 observations. At each interval, we recorded the location of the female as being on Qv, on Qg, or on the surface of the experimental arena. Research in other study systems (Abrahamson and Weis, 1997) and observational data from this system both support the underlying assumption of our current measurement that the proportion of time spent on a specific host plant is correlated with ultimate female oviposition decisions (r = 0.85, P < 0.001, Supplementary Figure S2). Therefore, quantifying the proportion of time an individual spends on its native host plant is an effective way to quantify the degree of host preference. In total, 557 host preference assays were conducted using 126 B. treatae, 130 B. fossoria, and 301 B. kinseyi. Only those trials in which individuals were observed on a host plant at least once were retained for analysis (76.5% of trials, Supplementary Table S1). In total, we made 8,355 individual observations over 278.5 testing hours for the assessment of host preference.

To compare the strength of habitat isolation among species, we conducted two complimentary analyses. First, we calculated individual host preference values by dividing the number of time points where a female was observed on its natal host plant by the total number of time points observed on either host plant during each trial as follows: host preference for Qv = (number of observations on Qv)/(number of observations on Qv + number of observations on Qg). These calculations yield values of host preference range from 0 to 1, with values less than 0.5 indicating increasing preference for the non-natal host plant, values of 0.5 indicating no preference, and values greater than 0.5 indicating preference for the natal host plant. The host preference values were then compared using a generalized linear mixed model (GLMM) with the response variable “host preference” assigned a beta-binomial distribution to control for overdispersion of the binomially distributed data (Kim and Lee, 2017; Length, 2023). The independent variables “gall wasp species” and “arena type” (plastic cup or Petri dish) were included as fixed effects, with “collection site” treated as a random effect. This analysis allowed us to explore species-level differences while controlling for population-level variation.

Second, we converted the point estimates of host preference for each species to pairwise estimates of the difference in host preference between species, a variable we termed “habitat isolation.” Following the approach described in Nosil et al. (2006), the strength of habitat isolation was quantified as the absolute value of the difference in mean host preference towards Qv between each species as follows: difference in preference between B. treatae and B. fossoria = absolute value of (host preference for Qv in B. treatae—host preference for Qv in B. fossoria). This species-pair metric ranges from 0 to 1, with 0 indicating similar habitat preferences and no habitat isolation and 1 indicating divergent habitat preferences and complete habitat isolation between species. To guide statistical inference, we then calculated 95% confidence intervals for the strength of habitat isolation from 10,000 bootstrap values for each of the three species pairs, B. kinseyi × B. treatae, B. kinseyi × B. fossoria, and B. treatae × B. fossoria. The bootstrap values were obtained by randomly resampling individual host preference values from the original data and recalculating habitat isolation (St. John and Fuller, 2021).

Estimating sexual isolation

To evaluate the potential role of differences in ecology (host plant use) and differences in geography (opportunities for migration and/or gene flow) on sexual isolation, we observed courtship and mating interactions between males and females for the two sympatric species, B. treatae and B. fossoria, as well as for the allopatric species, B. kinseyi, using controlled no-choice preference trials following the methods of Egan, Hood, Feder, et al. (2012) (Supplementary Table S3). Similar to the host preference tests, trials were conducted in 2016, 2017, and 2018 within 60 × 15 mm Petri dishes lined at the base with damp filter paper. For each trial, one male and one female were aspirated into the Petri dish and observed at 2-min intervals for 30 min for a total of 15 observations. At each interval, we recorded three courtship and mating-related interactions: (a) male wing buzzing (or fanning), an important male courtship display in insects (Villagra et al., 2011), (b) male mounting the female, a behavior that precedes copulation, and (c) copulation, defined as direct contact of male and female abdomens. In total, 1,123 mating assays were conducted using 225 female and 153 male B. treatae from Qv, 563 female and 455 male B. fossoria from Qg, and 335 female and 515 male B. kinseyi from Qv (see Supplementary Table S3 for pairings and replication). In total, we recorded 16,845 observations over 561.5 testing hours for the assessment of mate preference.

First, we calculated and then tested whether the probability of mating differed between conspecific and heterospecific pairs for each set of Belonocnema species with larger differences in the probability of mating equating to greater sexual isolation between species pairs. During our observations, the probability of mating for each individual pair was quantified as 0 if copulation did not occur and 1 if copulation occurred. Analysis of male-specific (wing buzzing) and female-specific (allowing male to mount) mate preferences are examined in Supplementary Material. The mating probability of the same species vs. different species treatments for each species pair was then compared using GLMM where the response variable “mating probability” was considered to be binomially distributed, with the independent variables “cross type” (conspecific or heterospecific) and “collection year” (2016, 2017, 2018, or 2019) included as fixed effects and “collection site” of males and females treated as a random effect. The three different heterospecific species pairs, B. kinseyi × B. treatae, B. kinseyi × B. fossoria, and B. treatae × B. fossoria, were compared with the corresponding conspecific mating trials for each species pairing (B. kinseyi × B. kinseyi, B. treatae × B. treatae, or B. fossoria × B. fossoria). Additionally, we also analyzed the degree of mate preference for males and females separately. These analyses and results can be found in Supplementary Materials.

Due to the lack of independence, the interaction between species pair type and individual pair cannot be directly tested using GLMM. Thus, we also quantified the strength of sexual isolation (SI) between species pairs using the method described by Sobel and Chen (2014): SI = 1 − 2 * [heterospecific mate frequency/(heterospecific mate frequency + conspecific mate frequency)]. In this metric, a value of zero indicates no sexual isolation, whereas a value of 1 indicates complete sexual isolation. Similar to the analysis of habitat isolation, we then generated 95% confidence intervals from 10,000 bootstrap values of sexual isolation between each of the three species pairs by randomly resampling mated and unmated pairs from the original data and recalculating sexual isolation.

All GLMM analyses were followed by multiple pairwise comparisons among Belonocnema lineages using a Tukey’s post hoc test with the function lsmeans in package “emmeans” (Length, 2023). All analyses including bootstrapped simulations were conducted in R version 4.0.2 (R Core Team, 2021).

Results

Strength of habitat isolation

Sympatric Belonocnema on Qv exhibited elevated levels of habitat isolation compared with the other sympatric or allopatric species, which is consistent with predictions that selection would favor higher host fidelity under sympatric conditions where immigration resulted in lower fitness. Specifically, populations of B. treatae developing on Qv that are sympatric with B. fossoria displayed significantly higher host plant fidelity (mean ± SE, 0.778 ± 0.042) than allopatric populations of B. kinseyi developing on Qv (0.647 ± 0.023, t = 2.49, P = 0.035, Figure 2A) and sympatric populations of B. fossoria, specialized on Qg (0.626 ± 0.037, t = 2.56, P = 0.029, Figure 2A). Correspondingly, we found the strength of habitat isolation to be the lowest between allopatric populations of B. kinseyi and B. treatae that share the same host plant species (mean ± 95% CI: 0.157, 0.064–0.248), followed by an intermediate value for allopatric populations of B. kinseyi and B. fossoria (0.183, 0.099–0.265) that use different host plant species, with the highest habitat isolation levels observed between sympatric populations of B. treatae and B. fossoria (0.340, 0.235–0.442) that use different host plant species (Figure 2B).

Interestingly, the trials performed in the smaller Petri dishes exhibited significantly higher host preference than trails conducted in the larger cups across our study (χ² = 7.376, P < 0.01), suggesting proximity to the two host plant options increased the expression of host preference, which is an important element that will inform future studies in this system. However, the testing environment did not affect the overall difference in host preference among different Belonocnema species, as indicated by the interaction term in the GLMM (χ² = 2.065, P = 0.356, Supplementary Table S4), suggesting that this difference did not affect our qualitative interpretations.

Strength of sexual isolation

Sympatric Belonocnema on Qg exhibited elevated levels of sexual isolation compared with the other sympatric or allopatric species, which is consistent with the prediction that selection favors stronger conspecific mate preference under sympatric conditions where hybridization results in lower fitness. The three Belonocnema species exhibited significantly higher probability of mating with conspecifics compared with heterospecifics that used different host plants (B. kinseyi × B. fossoria: Z = 2.511, P = 0.012; B. treatae × B. fossoria: Z = 2.279, P = 0.023; Figure 3A). In contrast, the probability of mating did not differ among allopatric conspecific and heterospecific B. kinseyi and B. treatae, which share the same host plant species (Z = 0.877, P = 0.380; Figure 3A).

Correspondingly, the strength of sexual isolation was lowest between allopatric B. kinseyi and B. treatae using the same host plant species (mean SI, 95% CI: 0.081, −0.137 to 0.302). Sexual isolation was ~1.7 times greater between allopatric B. kinseyi and B. fossoria using different host plant species (mean SI, 95% CI: 0.138, 0.013 to 0.259). The strength of sexual isolation was highest between sympatric populations of B. treatae and B. fossoria using different host plants (mean SI, 95% CI: 0.337, 0.068 to 0.598). Additionally, the degree of sexual isolation was significantly different among species pairs B. kinseyi × B. fossoria and B. treatae × B. fossoria where the 95% bootstrap CIs did not overlap with zero (Figure 3B).

Discussion

In this study, we tested for the contributions of divergent host use and geography to the evolution of habitat isolation and sexual isolation (Funk et al., 2002; Nosil, 2005, 2007). These two prezygotic barriers constitute important components of RI among the three closely related species of Belonocnema gall wasps, which constitute an emerging model for speciation research (Driscoe et al., 2019; Egan, Hood, Feder, et al., 2012; Egan, Hood, & Ott, 2012; Hood et al., 2019; Zhang, Hood, Carroo, et al., 2021; Zhang, Hood, Roush, et al., 2021). For both habitat and sexual isolation, we found that the strength of RI between species was lowest among allopatric species that shared the same host plant association, higher between allopatric species with different host plant associations, and highest between sympatric species with different host plant associations as illustrated in Figures 2 and 3. This pattern of RI in relation to geographic context and host plant association is consistent with previous work documenting an important role of both divergent ecology and geography in the process of speciation (e.g., Coyne and Orr, 2004; Nosil, 2007).

Our finding of significant habitat isolation and sexual isolation among allopatric species that share the same host plant association suggests that selection is not directly associated with host plant use and neutral processes associated with genetic drift can also modestly contribute to the evolution of both reproductive barriers. Significant estimates of prezygotic RI between allopatric species with the same host plant association indicate a potentially more complicated evolutionary outcome in a system where divergent selection from different host plant environments has been the main focus (e.g., Driscoe et al., 2019; Egan, Hood, Feder, et al., 2012, Egan, Hood, & Ott, 2012, Hood et al., 2019; Zhang et al., 2017, 2019; Zhang, Hood, Carroo, et al., 2021; Zhang, Hood, Roush, et al., 2021). For instance, RI could be an indirect product of the fixation of different advantageous mutations in separate populations experiencing similar selection environment, a process known as mutation-order speciation (Nosil, 2012; Schluter, 2009). In the case of allopatric B. treatae and B. kinyesi that share the same host plant, each lineage may have developed different phenotypic pathways underlying adaptation to the same host plant environment, and these differing solutions could generate assortative mating within lineages leading to partial habitat and/or sexual isolation (Nosil and Flaxman, 2011; Schluter, 2009). Given the complex genetic pathways involved in gall formation (e.g., Schultz et al., 2019), particularly in Hymenoptera (Martinson et al., 2022), it might not be surprising that lineages differ in their adaptations to manipulate their host.

The estimates of habitat isolation and sexual isolation observed between allopatric species feeding on different host plants compared with allopatric species feeding on the same host plant species strongly suggest that divergent selection can promote RI regardless of the geographic context. Allopatric populations of B. kinseyi sampled from Qv exhibit significant habitat isolation with respect to the alternative host plant Qg (mean host preference and 95% CI: 0.629, 0.581 to 0.676). Such habitat isolation in allopatric species might evolve regardless of the geographic context when there is a fitness cost for searching for hosts (e.g., predation risk or highly variable environments), and selection increases the host search efficiency and/or reduces risk via increased host fidelity (e.g., Egan and Funk, 2006). Similarly, the significantly higher estimate of sexual isolation among allopatric species pairs with different host associations (B. kinseyi × B. fossoria) compared to allopatric species with the same host associations (B. kinseyi × B. treatae) might arise as a by-product of phenotypic divergence (e.g., body size) due to local adaptation to different host plants (Egan, Hood, Feder, et al., 2012, Roush et al., in review).

We also observed that both habitat isolation and sexual isolation were highest among sympatric species with different host plant associations compared with allopatric species with different host plant associations (Figures 2B and 3B). There are three notable mechanisms in the speciation literature suggested to promote this pattern between sympatric B. treatae and B. fossoria: (a) the Templeton effect, (b) reinforcement, and (c) reproductive interference (Hollander et al., 2018; Noor, 1999; Templeton, 1981). The Templeton effect, also known as differential fusion, posits that isolated populations experiencing strong prezygotic RI persist in sympatry, whereas weakly isolated populations fuse (Templeton, 1981). Alternatively, reinforcement is the process where natural selection increases RI between two lineages as a result of selection acting against the production of hybrids with lower fitness (Noor, 1999) and has been demonstrated as an important factor promoting speciation in many systems (Servedio and Noor, 2003), particularly in insects (Yukilevich, 2012). Lastly, reproductive interference posits that the cost of migration or the presence of congeners during mating could promote the elevation of prezygotic isolation between the two sympatric species even without current gene flow (Butlin, 1987; Hollander et al., 2018). Future work in this system will need to address the role that these factors play in the speciation process in Belonocnema.

Interestingly, we observed asymmetry in the direction of habitat isolation and sexual isolation between B. treatae and B. fossoria (Figure 2A, Supplementary Figure S3) and that the direction of asymmetry is consistent with the direction of asymmetry in migration rate among these species (Driscoe et al., 2019). More specifically, extensive population genomics of Belonocnema across its known range has demonstrated lower migration rates of B. fossoria from Qg to the sympatric alterative host plant Qv, in comparison to the migration of B. treatae from Qv to the sympatric alternative host plant Qg (Driscoe et al., 2019; Zhang, Egan, et al., 2021). Correspondingly, the lineage that emigrates more (B. treatae) is predicted to evolve stronger habitat isolation due to the lower fitness of immigrants in the alternative habitat (Figure 2A), while the lineage that receives more immigrants (B. fossoria) is predicted to exhibit stronger sexual isolation due to the lower fitness of hybrids (Supplementary Figure S3). Despite asymmetric migration and asymmetric RI being commonly found across different study systems (Bolnick et al., 2008; Lowry et al., 2008a; Oswald et al., 2017) and the understanding that migration rate is important with respect to reinforcement (Servedio and Noor, 2003), few studies have formally described the relationship between asymmetry in migration and asymmetry in the evolution of RI (but see Suni and Hopkins, 2018; Yukilevich, 2012).

The results presented here serve as a valuable case study to evaluate the relative contributions of ecology and geography to two critical prezygotic reproductive barriers during the late stage of speciation. Although the method of contrasting lineages with different ecologies and geographies is primarily applied among populations at early stages of speciation, where the effects of evolutionary history and neutral genetic divergence are negligible (reviewed by Funk et al., 2002), we argue that this comparative framework remains informative even for diverged speciation pairs. In fact, this comparative framework should be considered an important initial step to explore evolutionary factors contributing to the completion of the speciation process, which is a critical knowledge gap in speciation research (reviewed by Kulmuni et al., 2020). Furthermore, our study is one of the few cases that identify the possible role of sympatry in promoting the evolution of habitat isolation, a common prezygotic barrier that has often been overlooked in studies investigating the occurrence of reinforcement (but see Nosil and Yukilevich, 2008; Nosil et al., 2006; Yukilevich and True, 2006). In future studies, it is essential to apply demographic history analysis to the three species pairs to further isolate the role of evolutionary history in shaping the evolution of RI.

Data availability

All datasets and R code can be found in GitHub repository using this https://doi.org/10.5281/zenodo.10278414.

Author contributions

L.Z. and S.P.E. designed the study and all authors collected data. L.Z. analyzed the data and wrote the first draft of the manuscript. All authors edited subsequent drafts and approved the final version.

Funding

Support was provided to S.P.E. by Rice University, a Rosemary Grant Award from the Society for the Study of Evolution (L.Z.), student research awards from the Society for Integrative and Comparative Biology, and the American Association of Naturalist (L.Z.).

Acknowledgments

We thank the staff and volunteers at Archbold Biological Station, especially Mark Deyrup and Hilary Swain, for field assistance in Florida; Robert Busbee, Amanda Driscoe, and Chelsea Smith for plant husbandry; and Gaston Del Pino, Sean Liu, and Elaine Hu for data collection. The authors acknowledge Lake Lizzie Conservation Area in Osceola County and South Florida Water Management District at Kissimmee River for access to experimental sites.

Conflicts of interest

The authors declare no competing interest.

References