Context-Dependent Reproductive Isolation: Host Plant Variability Drives Fitness of Hybrid Herbivores

The role of divergent selection between alternative environments in promoting reproductive isolation (RI) between lineages is well recognized. However, most studies view each divergent environment as homogenous, thereby overlooking the potential role within-environment variation plays in RI between differentiating lineages. Here, we test the importance of microenvironmental variation in RI by using individual trees of two host plants, each harboring locally adapted populations of the cynipid wasp Belonocnema treatae. We compared the fitness surrogate (survival) of offspring from hybrid crosses with resident crosses across individual trees on each of two primary host plants, Quercus virginiana and Q. geminata. We found evidence of weak hybrid inviability between host-associated lineages of B. treatae despite strong genomic differentiation. However, averaging across environments masked great variation in hybrid fitness on individual trees, where hybrids performed worse than, equal to, or better than residents. Thus, considering the environmental context of hybridization is critical to improving the predictability of divergence under variable selection.


Introduction
Identifying reproductive barriers that reduce gene flow between differentiating lineages is critical to gaining insight into mechanisms that initiate and/or maintain divergence (Ramsey et al. 2003;Coyne and Orr 2004;Rundle and Nosil 2005;Lowry et al. 2008). Rarely, however, do studies explore the basis of variation in estimates of reproductive isolation (RI) among populations. Given that variation in the strength of RI can provide stepping-stones for individuals migrating between lineages and that the effects of subsequent hybridization and introgression can be consequential for divergence patterns (Slatkin 1987;Arnold 2016), estimating variability of RI and understanding its underlying causes are crucial to predicting the evolutionary dynamics of interacting lineages (Zuellig and Sweigart 2018).
To date, studies of variation of RI have focused largely on genetic polymorphism of Bateson-Dobzhansky-Muller incompatibilities of hybridizing lineages as the primary mechanism driving variation of postzygotic RI (Reed and Markow 2004;Cutter 2012;Larson et al. 2018). In contrast, the degree to which selection varies among microenvironments within alternate environments and its concomitant effects on the strength of RI are relatively unexplored. Heterogeneity within alternate environments may generate variation in RI among diverging lineages by changing species interactions or the strength of selection (Wade and Kalisz 1990;Sweigart et al. 2007;Scopece et al. 2010;Suni and Hopkins 2018). Therefore, the connection between microenvironment and strength of RI is particularly important for understanding adaptation to divergent ecological conditions (Schluter 2000;Rundle and Nosil 2005;Nosil 2012). Alternate environments occupied by diverging lineages are typically viewed as homogenous entities (Richards and Ortiz-Barrientos 2016;Lackey and Boughman 2017;Karrenberg et al. 2019). However, as illustrated by Abrahamson and Weis (1997), Egan and Ott (2007), Egan et al. (2011), andMaciejewski et al. (2020), it is crucial to examine variation in selection and local adaptation within ecologically divergent environments. To date, few studies have explored how environmental variation affects RI between lineages (but see Seehausen 1997;Craig et al. 2007;Grant and Grant 2008). Here, we directly assess the link between heterogeneity within alternative environments (individual trees) and variation in the strength of RI in an insect herbivore.
Host-specific phytophagous insects are suitable for exploring the role environmental variation plays in promoting variation in RI because (1) alternative host plant species can represent divergent selective environments (Drès and Mallet 2002;Funk et al. 2002) and (2) interindividual variation in plant susceptibility and quality within each host species represents a form of habitat heterogeneity (Strauss and Karban 1998;Craig et al. 2007;Egan and Ott 2007). Here, we test how the strength of the postzygotic reproductive barrier "hybrid inviability" of the gall-forming specialist Belonocnema treatae varies among individual trees (i.e., different developmental environments) within each of their primary host plants, Quercus virginiana and Q. geminata. Host-associated lineages of this herbivore exhibit ecological divergence in gall and adult morphology (Egan et al. 2012a(Egan et al. , 2013 as well as multiple reproductive barriers (i.e., temporal, habitat, and sexual isolation; immigrant inviability; and reduced immigrant fecundity; Egan et al. 2012aEgan et al. , 2012bZhang et al. 2017Zhang et al. , 2021bHood et al. 2019). While no evidence of genomic admixture is present, low levels of migration between host-associated lineages is observed (Driscoe et al. 2019). This live oak-gall wasp system therefore allows us to test (1) whether these genomically divergent lineages of the herbivore exhibit hybrid inviability and (2) whether variation among individual trees contributes to the magnitude and/or direction of hybrid fitness.

Study System
The gall-forming wasp Belonocnema treatae Mayr (Hymenoptera: Cynipidae) develops on two sister species of live oaks (Fagaceae: Quercus: subsection Virentes), Quercus virginiana Mill (Qv) and Q. geminata Small (Qg), whose geographic ranges broadly overlap across the southeastern United States, creating sympatric and allopatric populations (Cavender-Bares and Pahlich 2009). The geographic distribution of B. treatae parallels the distribution of its hosts (Driscoe et al. 2019). In sympatry, host-associated lineages exhibit multiple prezygotic reproductive barriers that collectively reduce estimated gene flow by 95%-99% (Egan et al. 2012a(Egan et al. , 2012bEgan et al. 2013;Zhang et al. 2017Zhang et al. , 2021bHood et al. 2019). Correspondingly, sympatric host-associated lineages are genomically distinctive (average G ST p 0:375 0:006 confidence interval [CI]), with evidence of migrants but no admixture. Allopatric populations of B. treatae from the western range of Qv exhibit even stronger divergence from Qg-associated lineage (G ST p 0:8150:006 CI; Driscoe et al. 2019) but lower prezygotic isolation (Egan et al. 2013; L. Zhang, G. R. Hood, J. R. Ott, and S. P. Egan, unpublished data). Here, we define B. treatae populations from sympatric Qg (east), sympatric Qv (east), and allopatric Qv (west) as three lineages that are genomically dis-tinctive. We then compared the fitness of offspring produced among these three lineages to isolate the role of host association from intertree variability while controlling for factors such as geography in influencing divergence (Funk 1998).

Life Cycle and Sample Collection
Belonocnema treatae exhibits cyclical parthenogenesis (Lund et al. 1998). The sexual generation emerges from multichambered root galls in the spring, immediately mates, and oviposits into the veins on the underside of newly flushed leaves to initiate gall induction (Zhang et al. 2017;Hood et al. 2019). Single-chambered leaf galls develop during summer and mature between October and December, when asexual females emerge and oviposit into rootlets to complete the life cycle. For this study, newly emerged sexual B. treatae were collected in the spring of 2017 and 2019 from three allopatric and four sympatric Qv sites as well as three sympatric Qg sites (table S1, fig. S1A; tables S1-S6, figs. S1-S3 are available online). Specifically, unmated females were sampled by isolating individual root galls, which contain either all-female or all-male broods within fruit fly culture vials. To synchronize mating of adults sampled from Qv, which typically emerge earlier than Qg, half of each collection of root galls from each Qv site were temporarily stored at 47C to delay development. Additionally, we used the natural variation in leaf flush among the Qv and Qg trees employed in this experiment to match B. treatae emergence of host-associated populations ).

Experimental Methods
To evaluate hybrid inviability between Qv-and Qg-associated lineages, we estimated and compared the fitness surrogate (survival) of hybrid and resident crosses developing in each host environment. Resident crosses represent matings between individuals that feed on the same native host plant, whereas hybrid crosses employ males and females from different host plants in both directions (e.g., ♀ Qv #♂ Q g , ♀ Q g #♂ Qv ). For each replicate, we placed two males with three to four females inside a 30#46-cm mesh enclosure that was secured around a live oak branch containing newly flushed leaves of individual plants in a common-garden environment. In 2017, we mixed males and females from many different galls in a cross to increase the genetic diversity. In 2019, to increase the likelihood of mating, we placed the pairings into a 2-L plastic container for 48 h prior to the enclosure. To evaluate the effect of individual trees on wasp survival, both residents and hybrid crosses were deployed onto the same individual trees. The number of replicates (enclosures) per tree depended on the availability of newly flushed leaves. In total, 305 crosses were deployed across 66 individual trees (table S2; for details of potted sapling trees, see Hood et al. 2019). In the fall of 2017 and 2019, leaves with evidence of oviposition and/or galls were harvested from each replicate and then inspected by observers with no knowledge of cross details.

Survival Measurement
We evaluated survival per replicate as (1 2 the proportion of host plant immune response per oviposition event). Oviposition by B. treatae produces a permanent "scar" on lateral leaf veins, providing a record of oviposition attempts per leaf/replicate (Egan et al. 2007;Zhang et al. 2017). Oviposition may elicit a host immune response characterized by observable local tissue necrosis termed "hypersensitive response" (HR; Fernandes 1990Fernandes , 1998, which reduces the incidence of gall formation (Fernandes and Negreiros 2001), including for B. treatae (fig. S1B; Campbell 2010; Zhang et al. 2021b). Thus, the proportion of oviposition events exhibiting HR is a critical component of gall former fitness at this early life stage. In total, we examined 51,311 oviposition scars and 12,051 HRs from 10,840 leaves (table S2).

Test of Hybrid Inviability and the Effect of Individual Trees
To test whether hybrid crosses result in lower fitness compared with resident crosses, we used a generalized linear mixed model to compare "wasp survival following host immune response" out of the total number of oviposition scars per replicate modeled with a binomial probability distribution. The factors "host plant environment" (Qv or Qg), "cross type" (resident or hybrid), and "geographic context" (sympatric or allopatric) were considered as fixed effects, while "source population" was considered as a random effect. Geographic context did not explain a significant amount of variation and was excluded from the final model (see table S4, fig. S2).
To assess the effect of individual trees on hybrid inviability, we included individual tree identity (N p 1-66) as a random slope over cross types and tested the significance of the random slope in the model. Following Grossenbacher et al. (2015) and Martín-Robles et al. (2018), a significant random slope effect indicates that the direction and relative strength of hybrid viability varies among trees. To assess the significance of this random effect, we conducted likelihood ratio tests (function anova in R) among three models: (1) a model without individual tree (described above), (2) a model with individual tree as a "random intercept" (Gillies et al. 2006;Bolker et al. 2009), and (3) a model with individual tree as a "random slope" (i.e., "random slope over cross types"; Johnson 2014). All other modeled effects remained the same. Following Grossenbacher et al. (2015) and Martín-Robles et al. (2018), a significant effect of individual tree as a random intercept suggests that offspring fitness, regardless of cross type, varied by individual tree, while a significant effect of individual tree as a random slope over cross types suggests that differences in offspring fitness between cross types varied from tree to tree. Importantly, this latter test indicates that the direction and relative strength of hybrid viability varies among trees. We found no significant deviation in the dispersion of residuals in our model ( fig. S3) using the DHARMa package in R (Dunn and Smyth 1996). To test whether the significant effect of individual tree as a random slope over cross types depended on individual trees with larger sample sizes and thus weighted more in the model, we conducted a crossvalidation analysis by excluding one tree out of the model at a time. We then compared the Akaike information criterion (AIC) value between the model with individual trees as a random intercept and the model with individual trees as a random slope.
We also calculated the marginal and conditional R 2 values for the best fitted model to compare variance explained by fixed effects and by both fixed and random effects combined, using the function tab model in the R package sjPlot (Lüdecke 2020). Least square means for survivorship from the full model were then used to conduct multiple pairwise comparisons between cross types with Tukey's post hoc test in the emmeans package (Lenth et al. 2020).
Finally, when multiple replicates of both resident crosses and hybrid crosses were available on the same tree (table S3), we implemented a complementary procedure to test the effect of individual trees on hybrid inviability. Here, we conducted a generalized linear model with the response variable per replicate (total number of oviposition scars without HR) with a binomial distribution and cross type as a fixed effect. By doing this, we compared whether the direction of hybrid fitness relative to resident crosses was uniformly distributed across different trees. For this subset analysis, we examined 167 crosses on 25 individual trees.

Results
Survival of hybrid cross offspring did not differ from resident crosses when development proceeded on Qv. However, in the Qg environment, hybrid fitness was reduced when females were from the alternative host (Qv), reflecting a higher rate of host immune response than residents ( fig. 1; tables S5, S6; Tukey's test, resident vs. hybrid ♀ Qv #♂ Qg : t p 21:499, P ! :001).
In contrast, individual host tree explained a large proportion of variation in hybrid fitness (table 1; fig. 2). Including individual tree as a random slope over cross type significantly increased the model fit compared with model 1 (individual tree excluded, table 1; P ! :001) and model 2 (individual tree included as a random intercept, table 1; P ! :0001). Furthermore, our cross-validation results show that the DAIC between model 2 and model 3 ranged from 328 to 446 with a mean of 425, suggesting that individual trees with a large sample size did not influence the results of model comparison. Thus, the relative fitness of all cross types varied significantly across individual trees. Fixed effects models that included only host plant environment (Qv and Qg) and cross type (resident and hybrid) explained !5% of the variance in survival (table S5). Notably, when individual tree was considered as a random slope over cross types (full model), the total variance explained increased by 1450% (marginal R 2 p 0:045, conditional R 2 p 0:251). Despite the significant difference in immune response between 2017 and 2019 (table S5; P ! :0001), the direction of relative hybrid fitness on each of the three trees used in both 2017 and 2019 remained the same despite differences between individual trees. This observation is consistent with the individual tree environment influencing hybrid viabil-ity more than other environmental factors (e.g., betweenyear variation).
The regression coefficients for cross types (resident and hybrid) for each individual tree exhibited a nonuniform distribution (sign test for direction of slopes: P p :69). Among the 25 trees with multiple replicates of residents and hybrid crosses, fitness differed significantly between cross types in 16 cases (the probability of 16 false positives at a p :05 p 2#10 215 ). Hybrid inviability of Belonocnema treatae (lower fitness of hybrids than residents) was exhibited on nine trees, while hybrid vigor (higher fitness of hybrids than residents) was detected on seven trees.

Discussion
We found evidence of variability in the direction and strength of hybrid viability across individual hosts. On some trees hybrid crosses exhibited lower fitness than residents consistent with hybrid inviability, while on other trees hybrid crosses exhibited fitness equal to or higher than Figure 1: Fitness of resident, hybrid crosses (migrant female with resident male and resident female with migrant male; least squared means with 95% confidence interval) for each host plant environment based on the individual tree as a random slope. Different letters denote statistical significant at P ! :05.  residents ( fig. 2). The variation in hybrid fitness across individual trees, especially hybrid vigor, could limit RI between host-associated populations in sympatry. In contrast, little variation in hybrid fitness was explained by host association or geography. This latter result is surprising given the large degree of genomic differentiation among lineages on different host plants (G ST p 0:31) or between geographic regions (G ST p 0:81; Driscoe et al. 2019). The high variation in relative fitness of hybrids across individual trees could be due to differences in overall host plant defense (Egan and Ott 2007), which overrides differences between hybrid and residents. To address this a posteriori hypothesis, we examined the association between host immune response induced by residents on individual trees, a proxy for individual host defense, and the regression coefficient of cross types for each individual tree ( fig. 2), which estimates the fitness difference between resident and hybrid crosses (where negative values demonstrate lower hybrid fitness). We observed across all trees from both host species that the individual trees with the highest immune response to resident offspring tended to exhibit lower immune responses to hybrid offspring (Pearson's correlation, r p 20:58, t p 3:49, P p :002). This general association between resident and hybrid fitness on individual trees holds true when testing the correlation separately for each host plant species (Qv: r p 20:45, one-tailed t p 1:68, P p :060; Qg: r p 20:61, one-tailed t p 2:59, P p :0125). This result implies that those individual host plants that inflict higher degrees of immune response to locally adapted wasps inflict lower immune responses to hybrid offspring. One exception to the observation that individual host plant variation was the main driver of variation in hybrid fitness occurred when offspring developed on Qg. On this host, when females of the hybrid cross were from Qv, hybrid offspring exhibited lower average viability (i.e., induced higher host immune response) compared with offspring from residents, after accounting for the effect of individual trees. Interestingly, we found the reduction of fitness in one hybrid cross direction (♀ Qv #♂ Q g ) on Qg but not in the other cross direction (♀ Q g #♂ Qv ); such a pattern could also be a product of maternal effects rather than hybrid inviability, as the offspring produced by hybrid cross ♀ Qv #♂ Q g and ♀ Q g #♂ Qv should be genetically similar.
The subtle hybrid inviability observed between hostassociated lineages of Belonocnema treatae at the early life stage is unlikely to be caused by genetic incompatibility and is more consistent with ecological postzygotic isolation (Egan and Funk 2009). First, rather than observing hybrid inviability regardless of host plant environment as an expectation of genetic incompatibility, reduction in the fitness of hybrids of the same cross type (e.g., females from Qv mated with males from Qg) were found only in the Qg environment (table S4; Craig et al. 2007;Egan et al. 2011). Second, no evidence of hybrid inviability was found between allopatric same-host crosses between Qv-associated lineages despite significant genomic divergence of Qv-associated lineages (G ST p 0:81; Driscoe et al. 2019). Such a lack of genetic incompatibility despite the substantial amount of genetic divergence observed between B. treatae lineages is surprising given the observed genetic differentiation (Hartl and Clark 1997;Frankham et al. 2002). We tracked hybrid fitness only at this early stage of the insect life history, during avoidance of the host immune response (Campbell 2010). Thus, it is possible that as yet unidentified forms of hybrid inviability occur during later stages of the gall wasp's life cycle or in later generation hybrids (F 2 and beyond). Another possible explanation for the lack of genetic incompatibility is that haplodiploid organisms, such as B. treatae, have a unique association between the strength of genetic incompatibility and genetic differentiation (Koevoets and Beukeboom 2009;Beukeboom et al. 2015;Bendall et al. 2020).
While a few studies demonstrate that environmental context (mostly temporal variation) can dramatically change the course of speciation by reversing speciation or generating new species through hybridization (Seehausen et al. 1997;Grant and Grant 2008), our study provides one of the few examples in which RI is highly context dependent at a small spatial scale (among individual trees). Such context-dependent RI might be a common pattern during population divergence and speciation, particularly among plant-feeding insects (Mitter et al. 1988;Wiens et al. 2015), given the common observation of variation among individual plants in susceptibility to herbivorous insects (e.g., Gloss et al. 2016;Barker et al. 2018;Verçosa et al. 2019). Furthermore, considering that many speciation studies have documented partial RI yet no genetic differentiation between lineages experiencing divergent environments (Nosil 2012;Lackey et al. 2017), context-dependent RI could be an important but underexplored aspect constraining lineage divergence and speciation.

Data and Code Availability
Data are available from the Dryad Digital Repository (https://doi.org/10.5061/dryad.v41ns1rt6; Zhang et al. 2021a). Data and code are also provided in a zip file, available online. 1 Literature Cited