Asymmetric, but opposing reductions in immigrant viability and fecundity promote reproductive isolation among host‐associated populations of an insect herbivore

Immigrant inviability can contribute to reproductive isolation (RI) during ecological speciation by reducing the survival of immigrants in non‐native environments. However, studies that assess the fitness consequence of immigrants moving from native to non‐native environments typically fail to explore the potential role of concomitant reductions in immigrant fecundity despite recent evidence suggesting its prominent role during local adaptation. Here, we evaluate the directionality and magnitude of both immigrant viability and fecundity to RI in a host‐specific gall‐forming wasp, Belonocnema treatae. Using reciprocal transplant experiments replicated across sites, we measure immigrant viability and fecundity by comparing differences in the incidence of gall formation (viability) and predicted the number of eggs per female (fecundity) between residents and immigrants in each of two host‐plant environments. Reduced immigrant viability was found in one host environment while reduced immigrant fecundity was found in the other. Such habitat‐dependent barriers resulted in asymmetric RI between populations. By surveying recent literature on local adaptation, we find that asymmetry in immigrant viability and fecundity are widespread across disparate taxa, which highlights the need to combine estimates of both common and overlooked barriers in cases of potential bi‐directional gene flow to create a more comprehensive view of the evolution of RI.

Divergent natural selection between environments can play a major role in initiating population divergence, reducing gene flow in the intermediate stages of local adaptation, and/or maintaining species-level differences in the later stages of speciation (Rundle and Nosil 2005;Via 2009;Nosil 2012;Ingley and Johnson 2016;Lackey and Boughman 2017). Many studies spanning a diversity of taxa have examined a single reproductive barrier between lineages inhabiting different environments, revealing the common effect of ecology in facilitating and/or maintaining speciation (Huber et al. 2007;Bolnick et al. 2009;Dopman et al. 2009). However, to more fully understand the relative contribution of ecologically based divergent selection as an isolating mechanism between lineages, a comprehensive inspection of multiple reproductive barriers within single study systems is required (Ramsey et al. 2003;Nosil 2007;Richards and Ortiz-Barrientos 2016;Lackey and Boughman 2017;Karrenberg et al. 2019). To accomplish this goal in a given system requires a detailed understanding of natural history (Futuyma 1998), coupled with an assessment of commonly measured forms of RI, such as habitat isolation (e.g., Egan et al. 2012a;Bolnick et al. 2009), sexual isolation (e.g., Egan et al. 2012b;Nosil et al. 2002), temporal isolation (Dopman et al. 2009;Hood et al. 2019), as well as more system-specific or rarely measured reproductive barriers, such as cryptic isolation (Nosil and Crespi 2006) or extrinsic postzygotic isolation (Egan and Funk 2009;McQuillan et al. 2018).
Immigrant inviability is a taxonomically widespread reproductive barrier that can contribute significantly to RI (Via et al. 2000;Nosil et al. 2005;Ingley and Johnson 2016;Richards and Ortiz-Barrientos 2016). Defined as a reduction in gene exchange between populations due to selection against immigrants in nonnative environments, immigrant inviability is especially important in initiating ecological speciation as it represents a direct product of divergent adaptation to different environmental conditions . While selection against immigrants could manifest as higher mortality and/or lower fecundity in nonnative environments, most studies have focused specifically on viability (Nosil 2012;Gosden et al. 2015;Ingley and Johnson 2016;Kaufmann et al. 2017), despite both empirical and theoretical evidence suggesting that fecundity may play an important role in isolation during local adaptation (Hereford et al. 2004;Hereford 2009;Porter and Benkman 2017). A greater reduction in fecundity than viability could occur when a divergent selection from different environments favors allocation of energy to maintain survival, but not reproductive effort (Porter and Benkman 2017). The few studies that have assessed the effect of both reduced immigrant viability and fecundity to reproductive isolation provide mixed evidence of the relative importance of these two reproductive barriers at the early versus late stages of speciation. For example, in host-associated ecotypes of Timema stick insects, reduced immigrant viability, but not reduced fecundity, contributes significantly to RI at the early stages of population divergence, while both barriers serve to isolate taxa at the later stages of speciation (Nosil and Sandoval 2008). On the other hand, the most recent survey of studies of local adaptation across diverse taxa found that reduced immigrant fecundity contributes to RI at a similar or stronger magnitude than reduced immigrant viability, particularly in the early stage of speciation (Porter and Benkman 2017). Clearly, more empirical studies that assess and compare the relative magnitudes of immigrant viability and fecundity as reproductive barriers in the same study system are needed to improve our general understanding of the roles these barriers play during the speciation process.
A second issue attendant to the exchange of immigrants is the directionality of gene flow between diverging populations. During the initial stages of divergent adaptation, barriers to gene flow tend to be asymmetrical (Lowry et al. 2008;Lackey and Boughman 2017). For example, the strength of reduced immigrant viability is often asymmetrical, with reduced immigrant viability often found in one environment, but not in the other (Räsänen and Hendry 2014;Gosden et al. 2015;Kaufmann et al. 2017). Asymmetries in the strength of other reproductive barriers also appear to be common during speciation (Lowry et al. 2008b;Lackey and Boughman 2017). The question of how populations maintain locally adapted traits and genetic divergence despite asymmetry in the strength of barriers to gene flow stresses the need for a comprehensive assessment of RI in both directions of gene exchange (Slatkin 1987(Slatkin , 1995Lackey and Boughman 2017). One possible mechanism is that multiple asymmetric, but opposing reproductive barriers evolve and cumulatively generate symmetrical reductions in gene exchange between populations.
In this study, we measure the magnitude and directionality of both immigrant viability and fecundity between host-associated populations of the gall-forming insect specialist, Belonocnema treatae. Populations of this species sampled from two sister species of oaks are genomically distinctive with low levels of immigration, but exhibit no measurable admixture in sympatry (Driscoe et al. 2019). Multiple prezygotic reproductive barriers between host-associated populations of B. treatae, including temporal isolation , habitat isolation (Egan et al. 2012b), and sexual isolation (Egan et al. 2012a), contribute to total RI which ranges from 0.81 to 0.99 across different population pairs, suggesting unmeasured barriers to complete RI. Additionally, Zhang et al. (2017) found that B. treatae populations from allopatric host plants have significantly lower rates of gall formation on non-native host plants, suggesting that divergent host use could exert strong selection on survival. On the basis of these results, we hypothesize that reductions in both immigrant viability and fecundity help explain the missing link between incomplete RI and distinct genotypic clustering among host-associated B. treatae populations.
To address these hypotheses, we first test whether immigrants have lower survival and fecundity on non-native host plants relative to residents among sympatric B. treatae populations. We then compare the strength of RI contributed by reduced immigrant viability and fecundity in each host environment. To accomplish this we conducted a reciprocal transplant experiment replicated across populations of B. treatae from two sister species of oaks, Quercus virginiana (Qv) and Q. geminata (Qg) in Florida, where the geographic ranges of these species overlap. Next, we compared the strength of reduced immigrant viability and fecundity in each host environment using pair-wise comparisons between populations to highlight how the two barriers can combine to generate RI. We then combined our estimates of RI with previous estimates of temporal, habitat, and sexual isolation among host-associated populations of B. treatae to assess the relative contribution of each barrier to total RI in this system (Ramsey et al. 2003). Finally, to place our results in a broader context, we used a recent compilation of studies documenting fecundity and viability of immigrants associated with locally adapted populations (Porter and Benkman 2017) to calculate the relative strength of each barrier, assess their directionality, and determine whether the barriers are asymmetric across different taxa.

STUDY SYSTEM
The gall-forming wasp, Belonocnema treatae Mayr (Hymenoptera: Cynipini: Cynipidae), is host-specific to a subset of live oak species (Fagaceae) in the subsection Virentes, including the sister taxa Quercus virginiana Mill (Qv) and Quercus geminata Small (Qg) distributed across the southeastern United States (Cavender-Bares et al. 2015;Fig. 1A). In regions of Florida where Qv and Qg co-occur, each host plant occupies different microhabitats (Cavender-Bares and Pahlich 2009): Qg inhabits drier, sandier, and higher pH soils than Qv, which is typically found in mesic, lower pH soils. Although each tree species specializes on slightly different soil types on average, there are many sites where both host plants co-occur within the "cruising range" of the wasps generating a blurry mosaic of host populations (Mayr 1947). The geographic distribution of B. treatae parallels the distribution of the host plant (Driscoe et al. 2019). In this study, we focus on documenting RI between host-associated populations of B. treatae in the region of sympatry of Qv and Qg in southcentral Florida (Fig. 1A).
Belonocnema treatae exhibits cyclical parthenogenesis wherein temporally and spatially segregated sexual and asexual generations alternate to complete an annual lifecycle, which is common among cynipid gall wasps (Stone et al. 2002). Sexual generation adult males and females emerge from multichambered root galls during the spring, coinciding with new leaf flush. Upon emergence, wasps mate immediately, females then oviposit into the newly flushed leaves, from which singlechambered galls arise that house the asexual generation (Lund et al. 1998;Figure 1B, 1C). Oviposition timed to match the onset of new leaf growth is critical for successful gall initiation (Zhang et al. 2017;Hood et al. 2019). Each ovipositor insertion by B. treatae generates a permanent scar on the lateral vein on the underside of each leaf, providing a record of oviposition attempts (Fig. 1D). Upon emergence in the fall, asexual females oviposit into roots inducing multi-chambered galls that give rise to the sexual generation to complete the life cycle (Lund et al. 1998).
Populations of B. treatae associated with sympatric Qv and Qg in Florida represent two highly diverged genetic clusters: average genetic distance (G ST ) between host-associated populations in sympatric areas is 0.371 (Driscoe et al. 2019). Moreover, this study documented limited migration among host plant species with no apparent admixture, suggesting migration between host plants is possible but that hybrid formation is rare or absent. In contrast, laboratory experiments have found little evidence of hybrid inviability (Zhang et al. unpublished data), indicating that prezygotic reproductive barriers are likely the main mechanisms maintaining reproductive isolation and genetic cohesion. Prior research has also documented multiple incomplete prezygotic reproductive barriers. First, leaf flush occurs 2−3 weeks earlier in Qv, promoting divergence in the emergence time of the sexual generation between host-associated populations and generating temporal RI . Second, sexual generation populations of B. treatae exhibit habitat preference for their native host plants, generating habitat isolation (Egan et al. 2012b, Egan et al. 2013. Third, assortative mating independent of host environments is present among populations with different host-associations, creating sexual isolation (Egan et al. 2012a, Egan et al. 2013. In total, these barriers combined to reduce gene flow by 81% to 99% between populations . Taken collectively, these results position the current study to examine whether reduced immigrant viability and/or fecundity can help explain the missing link between incomplete RI and distinct genotypic clustering among host-associated B. treatae populations. To test these hypotheses, we performed a set of reciprocal transplant experiments.

RECIPROCAL TRANSPLANT EXPERIMENTS
In early to mid-March, 2017, (prior to emergence of the sexual generation from root galls), at each of six sites (three Qg and Qv sites, respectively) distributed across southcentral Florida, we placed 32 Qv and 36 Qg 2-m tall, 5-year-old potted saplings of each host plant species under the canopy of old-growth live oaks that each possessed high densities of developing root galls (Fig. 1A, Table S1). To ensure that the potted plants would be available for oviposition by newly emerged B. treatae, we monitored sexual generation emergence at each site using yellow sticky traps placed under the host oaks beginning in early March thru late April and deployed experimental potted plants when emergence was detected. Naturally emerging B. treatae were allowed to oviposit into newly flushing leaves of the saplings over the course of three weeks at each site. Leaf flush and the emergence of sexual generation B. treatae are synchronized and occur 2-3 weeks later on Qg than on Qv. Therefore, we made use of the natural variation in the timing of leaf flush among individual potted saplings to match the sexual generation emergence dynamics from root galls (see Hood et al. 2019 for detailed methods of monitoring and matching the timing of leaf flush). As well, we manipulated leaf flush to ensure that the timing of tissue availability for Qv and Qg saplings at each location would be equivalent. Specifically, we advanced the phenology of 20 Qg saplings by placing them in a warm and humid greenhouse prior to the experiment and delayed leaf flush of 35 Qv saplings by placing them in the shade to reduce sunlight and lower ambient temperatures. Following exposure to oviposition, saplings were kept in an enclosed area at Archbold Biological Station until the end of April, and then transported to the greenhouse and placed in controlled conditions where galls were allowed to develop and mature (see below).
To guard against natural enemy mortality and/or unpredictable biotic or abiotic events in the field, we supplemented the field-based reciprocal transplant study with a similar, but more controlled, experiment enclosing a similar set of sapling trees within cages at Archbold Biological Station, Venus, FL, which is centered within the range of our experimental sites. This approach effectively eliminated natural enemy attack on developing B. treatae and precluded unwanted oviposition by naturally emerging B. treatae. We first placed five to eight saplings of each host plant species inside each of six Nytex screen cages (1.8 m × 3.66 m × 1.83 m), one for each of the six field sites described above. Informed by the timing of adult emergence at each site, we collected root galls that harbored the penultimate stage male and female B. treatae between mid-and late-March for Qv, and between mid-March and mid-April for Qg. Root galls from each site were isolated in 2-L plastic bottles that were housed outside under ambient environmental conditions under a shaded alcove at Archbold Biological Station. Adult emergence from galls was monitored daily and upon emergence, three to five males and three to five females were aspirated into individual Organdy mesh bags (30 cm x 46 cm) and placed on an oak branch with newly flushed leaves. Bags were removed after one week as adult life span of B. treatae ranges from one to seven days . A total of 309 replicates were deployed across 38 Qv and 30 Qg saplings (Table S2). All experimental saplings were maintained at Archbold Biological Station within the Nytex cages until transport to the greenhouse.
On May 1, 2017, all experimental saplings from the field and cage were placed in a greenhouse and maintained at 23°C under ambient light conditions and watered every other day. Aphid and other greenhouse pest activities on saplings were monitored every other day and physical removal (hand picking) and biological control (ladybug release) were conducted as needed. Greenhouse rearing conditions for all sapling trees were equivalent.

SAMPLE PROCESSING
In mid-November 2017, upon gall maturation, all leaves with oviposition scars and/or galls were harvested from each replicate prior to the emergence of asexual generation B. treatae from leaf galls for both the field and cage studies. All leaf galls from each replicate plant from the field experiment, and from each bagged branch corresponding to experimental replicates from each sapling in the cage study, were removed from leaves, placed inside a separate clear 0.5 L cup covered by a paper towel secured with a rubber band, and misted once daily to maintain humidity and moisture. Emergence of asexual gall wasps was monitored once every other day from mid-November to mid-January and all emerging wasps were preserved in 95% ethanol and stored at −20°C.

MEASURING SURVIVAL
We used two measures of immigrant and resident viability of B. treatae: host immune response to oviposition, and the proportion of oviposition events that resulted in gall establishment in native and non-native host plant environments. We define "immigrants" as eggs deposited into non-native host species and "residents" as eggs deposited into native host species. Oviposition into leaves is often followed by a hypersensitive-like immune response (HR) characterized by local tissue necrosis around the deposited egg (Fernandes 1990;Fernandes and Negreiros 2001). This response can reduce the incidence of gall formation (Fernandes and Negreiros 2001), including gall establishment in B. treatae (Campbell 2010 ; Fig. 1D). Thus, the proportion of oviposition events that exhibit hypersensitive response is a good indicator of one critical component of wasp performance (i.e., the ability to evade initial host plant defense; Zhang et al. 2017). In addition, we counted the number of galls, defined as threedimensional tumors-like structures, on each branch/experimental tree and used this information to estimate the proportion of oviposition events that resulted in gall establishment (a direct measurement of wasp survival at later life stages). We did not include the number of wasps produced per gall as a measurement of survival because the emergence rate of wasps in the laboratory is lower (< 20%) than natural populations when protected from natural enemies (∼60%) (Egan and Ott 2007;Hood and Ott 2009).

ESTIMATING FECUNDITY
To estimate potential fecundity of B. treatae as a function of the source population and rearing environment, we made use of the fact that asexual generation B. treatae are pro-ovigenic (i.e., females emerge as adults with the entire potential lifetime complement of eggs fully matured), and that body size is highly correlated with egg number (an estimate of potential fecundity) (Hood and Ott 2017). Potential fecundity was estimated using the relationship: potential fecundity = −225.4 + 423.9 * tibia length (R 2 = 0.85, P < 0.0001) established by Hood and Ott (2017), where hind tibia length is used as an index of body size, a common method for Hymenoptera (Honek 1993). Here, estimating fecundity using body size is independent of our estimation of viability (proportional gall formation) as all individual wasps used to measure body size formed and successfully emerged from a gall. Choosing independent measures of viability and fecundity are critical for comparison within a system or between other systems (Porter and Benkman 2017). The length of each tibia was digitally measured by using a LEICA M125 microscope at 100X magnification and LAS v4.4 software . In total, we measured tibia length and estimated fecundity for 991 individuals distributed across sites, host plants, and treatments (Table S3). We prioritized measuring individuals from the field treatment, where the rearing condition most closely approximated natural conditions. For sites where there were less than 100 individuals in the field treatment, we measured every individual and supplemented the sample with individuals from the cage treatment. For sites with more than 100 individuals in the field treatment, we sampled individuals evenly distributed across emergence dates. A comparison of field and cage reared wasps showed no significant differences in viability or fecundity metrics (see Supplemental Results and Fig. S2).

STATISTICAL ANALYSIS
To compare survival of immigrants and residents, generalized linear mixed models (GLMM) were conducted within the R package "glmmTMB" followed by multiple pairwise comparisons between residents and immigrants from different host environments via Tukey's post hoc test using the "emmeans" package (Brooks et al. 2017;Lenth et al. 2020). In all comparisons, we examined the a priori hypothesis that immigrants would exhibit lower fitness than residents due to reduced viability and/or fecundity, thus we performed one-tailed hypothesis testing. In our experimental design, we set up the field experiment and then, as a back-up, set up equivalent experiment in cages with bagged groups of males and females on a similar set of sapling trees from the same B. treatae populations. We first tested whether different rearing location (field or cage) affected wasp fitness (survival and fecundity). To accomplish this, we analyzed data from field and cage treatment separately. We found no significant changes in the qualitative results (see Supporting Information Results, Figs. S1 and S2). Therefore, we combined the field and cage data in the formal analysis by including the rearing location (field or cage) as a fixed effect to control its potential variation in the model for all analyses.
A separate GLMM was conducted for each component of survival (proportion of host immune response per tree and proportion of gall formation per tree) designated as the response variable with a binomial distribution. In each model, host environment (Qv or Qg), treatment (immigrant or resident), host environment × treatment, and rearing location (field or cage) were fixed effects, and sampling location was a random effect. We also tested for zero inflation under the assumed probability distribution using the function testZeroInflation in the R package "DHARMa". If significant, we included a zero-inflation function in the model.
To compare potential fecundity of immigrants and residents, a linear mixed model (LMM) was conducted using the R package "nlme" followed by multiple pairwise comparisons between residents and immigrants from the alternative host environments using Tukey's post-hoc test (Pinheiro et al. 2020). The LMM includes predicted number of eggs per individual as the response variable; host environment (Qv or Qg), treatment (immigrant or resident), host environment × treatment, and rearing location (field or cage) as fixed effects; and sampling locations as a random effect.
Lastly, to illustrate our results throughout, we calculated the least square means (lsm) of each response variable using the 'emmeans' package (Lenth 2020). All analyses were conducted with R v.3.6.1 (R Core Team 2020).

ESTIMATING INDIVIDUAL REPRODUCTIVE BARRIERS
We combined field and cage data to estimate the strength of RI in the form of reduced immigrant viability and fecundity for each of the six populations, as results suggest that analyzing the data sets separately does not change the conclusion of the study (see Supplemental Results, Fig. S1, S2). To calculate RI for both measures of viability and the single measure of fecundity, we used the method described by Sobel and Chen (2014)  (1) To estimate the strength of RI associated with reduced immigrant viability, we used the proportion of galls formed as the measure of fitness that estimates gall wasp survival at a later life stage than the degree of host immune response. To estimate the strength of RI associated with reduced immigrant fecundity, we used the predicted number of eggs contained per resident and immigrant female. Three Qv-and three Qg-associated populations colonizing both native and non-native host plants resulted in three populations of residents and three populations of immigrants, generating nine pairwise comparisons between different populations in each host environment (see Table S4). The strengths of reduced immigrant viability and fecundity in RI were estimated for each comparison. The average strength ± SE of RI was obtained from values across all pairwise comparison.

ESTIMATING TOTAL PREZYGOTIC RI
Total prezygotic RI between populations reared in each host environment was calculated by combining estimates of RI caused by reduced immigrant viability and immigrant fecundity with estimates of temporal isolation, habitat isolation, and sexual isolation from previous studies of B. treatae (Egan et al. 2012a;Egan et al. 2012b;Egan et al. 2013;Hood et al. 2019). To estimate total prezygotic RI, we calculated and then summed the relative contribution of each barrier after accounting for the cumulative RI that occurs before this barrier following Ramsey et al. (2003) using equation 2: where n is the total number of reproductive barriers, i is the i th reproductive barrier, RI i is the estimates of RI for i th reproductive barrier, and total RI n is the total RI summed over n reproductive barriers, where total RI n-1 is the total RI combined by the first n-1 reproductive barriers. Here, the order of each barrier was arranged according to the sequence in which they reduce gene flow (temporal isolation, habitat isolation, immigrant inviability, reduced immigrant fecundity, sexual isolation). We argue that reduced immigrant viability/fecundity occurs before sexual isolation because the immigrants are the eggs of B. treatae individuals oviposited into the non-native environment which occurs before the subsequent cycle of sexual generation mating between individuals from different host plants (sexual isolation). The strength of each reproductive barrier varied among population comparisons, therefore we estimated total prezygotic RI using the single lowest, the average, and the single largest pairwise estimate for each barrier following the methods of Hood et al. (2019).

ACROSS TAXA
To examine patterns of the strength, asymmetry, and directionality of reduced immigrant viability and fecundity across taxa, we extracted data from 16 studies spanning plants, arthropods, and cnidaria recently reviewed by Porter and Benkman (2017). While Porter and Benkman (2017) estimated RI by averaging the performance of residents and immigrants across different environments, we estimated the strength of RI within each environment based on our observations of habitat-specific and asymmetric barriers in B. treatae. Therefore, in each study, viability was estimated by individual survival rates, whereas fecundity was estimated by different measures of reproductive success (e.g., number of seeds, number of eggs). We used WebPlotDigitizer to extract data from figures in the literature when the raw data was not publicly available, which allowed us to retrieve 16 of the 18 studies from Porter and Benkman (2017). If survival was estimated throughout all life stages (e.g. seed germination, seedling, adult), viability was calculated as multiplicative survival across all stages. Additionally, if there were multiple population comparisons in the study, we calculated the mean across all comparisons. Lastly, if the study surveyed populations across an environmental gradient (e.g. altitude), we used data from the two most contrasting environments to capture the strongest estimates of selection against immigrants. The strength of reduced immigrant viability and fecundity in each environment was calculated following Sobel and Chen (2014) as described above. The magnitude and direction of asymmetry of each reproductive barrier were estimated for each study by calculating the differences in the strength of RI between two environments. To compare the degree to which a barrier was asymmetrical, we summarized observations of the absolute value of asymmetry, which ranges from 0 to 2, with 0 indicating equivalent RI in both directions for a specific barrier, and 2 indicating complete asymmetrical RI, where predicted gene flow is 1 in one direction while the probability of gene flow is 1 in the opposite direction. To determine whether reproductive barriers opposed or complemented each other, we used the sign of the degree of asymmetry, which ranges from −2 to 2. A value of 0 indicates symmetric gene flow in both directions whereas positive and negative values of asymmetry indicate stronger potential gene flow from population A to population B, and from population B to A, respectively. Opposing reproductive barriers were deemed present when the signs of the degree of asymmetry of reduced immigrant viability and fecundity were opposite (one positive and one negative value).

AND IMMIGRANT FECUNDITY TO RI
We observed reductions in immigrant viability across all nine pairwise comparisons in the Qg environment (sign-test: P < 0.01; Table S5). The average reduction in viability across pairwise comparisons translated into a predicted 60% decrease in gene flow from Qv-associated populations to Qg-associated populations (Fig. 4A, Table 1). In contrast, the strength of immigrant  inviability based on percent gall formation estimated for individuals moving from Qg to Qv was highly variable, and occasionally exceeded resident survival. Overall, predicted changes in potential gene flow from Qg-to Qv-associated populations ranged from −62% to 34% across population pairs, with an average of ∼0% change in gene flow across all pairwise comparisons in this direction (Fig. 4A, Table 1). The predicted fecundity of immigrants reared on Qg host plants varied across populations, with immigrants in some populations exhibiting higher fecundity than residents. Thus, changes in potential gene flow ranged from −11% to 14%, and averaged ∼0% from Qv-to Qg-associated populations (Fig. 4B, Table 1).
Conversely, lower immigrant fecundity was observed in the Qv environment across every population pair (sign-test: P < 0.01) resulting in a reduction in potential gene flow of 12.2 ± 5% (mean ± standard deviation) from Qg-associated B. treatae populations in the Qv environment (Fig. 4B, Table 1). Overall, the effect of immigrant inviability as a reproductive barrier to gene flow was approximately five times stronger than the effect of reduced immigrant fecundity.

POPULATIONS OF B. TREATAE
All prezygotic barriers associated with divergent host-use, including reduced immigrant viability and fecundity (this study), temporal isolation , habitat isolation (Egan et al. 2012b), and sexual isolation (Egan et al. 2012a), combined to significantly reduce predicted gene flow between host-associated populations of B. treatae, but RI varied by direction of migration and population comparison. For Qg-associated populations immigrating to Qv-associated populations, predicted gene flow was reduced by 78% to 100% with an average reduction of 95% (Table 2). Conversely, when Qv-associated populations immigrate to Qg-associated populations, the estimated reduction in gene flow ranged from 91% to 100%, with an average reduction of 99% (Table 2). Comparison of the absolute strength of each reproductive barriers shows that temporal isolation is the strongest barrier and is symmetrical with respect to the direction of gene flow (average = 0.82, range: 0.65 -0.95). Immigrant inviability is the second strongest barrier to gene flow of wasp populations associated with Qv immigrating onto Qg (average

ACROSS TAXA
The magnitude of asymmetry varied across the taxa examined, with the absolute value of asymmetry for immigrant inviability ranging from 0 to 1.5. In comparison, the absolute value of asymmetry for reduced immigrant fecundity ranged from 0 to 0.6. A magnitude of asymmetry above 0.5 was found in six of 17 studies of immigrant inviability, but only one study on reduced immigrant fecundity had a value of asymmetry above 0.5 ( Fig. 5A and B). Opposing reproductive barriers, where reduced immigrant viability and fecundity combine to form RI in both directions within a single system, were found in six of 17 studies (Fig. 5C).

Discussion
Reductions in gene flow between ecologically divergent populations as a result of decreased immigrant fitness in non-native environments may be a common and direct consequence of local adaptation promoting speciation (Via et al. 2000;Nosil et al. 2005;Ingley and Johnson 2016). While most studies have concentrated on the effect of reduced immigrant viability on RI, reduced immigrant fecundity is also a potentially important, but un- derstudied contributor to RI (Smith and Benkman 2007;Hereford 2009;Porter and Benkman 2017). This may be especially true for specialized herbivorous insects, where shifting to a different host plant might not affect survival, but instead could greatly reduce  individual growth thereby reducing fecundity (Bernays and Chapman 1994;Via et al. 2000;Nosil and Sandoval 2008). Our study presents an assessment of both reproductive barriers and their relative contribution to total prezygotic RI in a host-specific, gall-forming insect, B. treatae, on two sympatric host plants, Q. geminata and Q. virginiana. We found that both reduced immigrant viability and reduced fecundity diminish gene flow, but that both barriers were dependent on the host environment. Specifically, immigrants exhibited lower viability than residents on Qg, while immigrants had lower predicted fecundity than residents on Qv. Collectively, these two barriers combined to reduce gene flow in both directions of immigration, however, the strength of immigrant inviability was much stronger than decreased immigrant fecundity (0.60 versus 0.12). Therefore, the total strength of RI combined across barriers generated an overall greater reduction in predicted gene flow between populations when B. treatae from Qv-associated populations migrated to the Qg host plant environment.

Strength of RI for each of five reproductive barriers (Egan
To estimate and compare immigrant viability and fecundity within a single study, these measures of fitness variability must be independent. For the study reported herein, we deem our measures of immigrant fitness to be independent based on 1) the biology of the system and 2) the patterns observed in this study as well as the previous study of Egan et al. (2011). We measured viability by determining whether individual eggs and early instar larva avoided the hosts' immune defense. We then estimated fecundity using adult body size to calculate the predicted egg number. No adult can emerge without successfully forming a gall, therefore the adults we used to measure body size are all successful survivors in our viability measurement. Thus, by definition, viability and fecundity cannot be correlated. Moreover, if our measurement of fecundity was positively correlated with our measurement of viability, the direction of asymmetric RI caused by reduced immigrant viability and fecundity would be in the same direction. Instead, we found the opposite: viability differed in one environment only (Fig. 2), and fecundity differed in the alternative environment (Fig. 3). Finally, Egan et al. (2011) found that gall size, which is positively correlated with body size, is predominantly under stabilizing selection in the absence of parasitoids and other natural enemies. This mode of selection suggests that there is no intrinsic linear correlation between increased body size and increased survival in this system.

REDUCED IMMIGRANT VIABILITY AND FECUNDITY
A common explanation for the occurrence of reduced immigrant viability and fecundity is the fitness trade-off hypothesis: local adaptation to one environment comes at a cost of maladaptation in a different environment (Futuyma and Moreno 1988;Via et al. 2000;Nosil et al. 2005;Bono et al. 2017, but see Hereford 2009;Hall et al. 2010;Gompert et al. 2015). While this mechanism can hold true when symmetric reduced immigrant viability or reduced fecundity occurs, asymmetric reproductive barriers could be caused by other factors. For example, if one population simply outperforms the other population in either native or non-native environments (without exhibiting a fitness trade-off), this can lead to an asymmetric reproductive barrier (Hereford 2009;Räsänen and Hendry 2014;Gosden et al. 2015). Our study suggests that asymmetric immigrant inviability in B. treatae is driven by a combination of both mechanisms: individuals from both host plants formed fewer galls on non-native host plants than on native host plants, and individuals from native Q. geminata hosts have higher gall formation rates than individuals from native Qv in either host environment (Fig. 2B, Table S4). We also found that asymmetry in reduced immigrant fecundity may be driven by a combination of fitness trade-offs and differences in mean fitness of populations. Our results showed that individuals from the Qv host had lower predicted fecundity when developing on non-native host plants compared to native host plants, supporting the fitness trade-off hypothesis (Table S4). However, individuals from Qg tend to have lower fecundity on their native host plant species than individuals from Qv on their native host species, suggesting the differences in the mean fitness of populations inhabiting different host environments (Fig. 3, Table S4). These results are of interest as they contrast with the pattern we have observed in natural populations: wasps from Qg are significantly larger than wasps from Qv (Egan et al. 2012a). Differences in the body size of Qv-and Qg-associated B. treatae populations pattern may due to maladaptation to non-natal tree genotypes (Egan and Ott 2007), especially among Qg-associated gall wasp populations. Since asymmetric immigrant inviability has widely been reported across different study systems, more research is needed to address its underlying causes.

MULTIFARIOUS SELECTION ASSOCIATED WITH DIVERGENT HOST-USE
Multifarious selection in divergent environments has been proposed to play a critical role in ecological speciation (Rice and Hostert 1993;Nosil and Sandoval 2008;Nosil 2012). Treating different host plant species as different selective environments is a common practice in studying host-shift speciation events. Transplant experiments, while illuminating overall fitness differences between native and novel herbivore populations, provide few details about the specific mechanisms whereby different host plants exert selection on host-specific insects, since alternative host plants differ in many important aspects such as host-defense, nutrition, and natural enemy communities (Via 1990;Kawecki and Ebert 2004 Nosil and Crespi 2006) and physiological composition (Via et al. 2000;Nosil and Sandoval 2008) can drive strong selection on host-specific insects. However, few studies have investigated multi-dimensional selection on divergent host-use (but see Nosil and Sandoval 2008). In the B. treatae-live oak study system, Hood et al. (2019) showed that divergence in plant phenology results in strong selection on wasp phenology as mismatches in the timing of ephemeral leaf tissue required to initiate gall formation and adult wasp oviposition result in near-zero gall formation (Zhang et al. 2017). In this study, we found that the degree of host immune response parallels the proportion of galls formed, suggesting that the host plant immune response is also an important source of divergent selection on host-associated gall wasp populations. Moreover, we observed an overall reduction in body size of individuals when they developed on non-native host trees for both Qv-and Qg-associated populations. Such a reduction of body size suggests that differences in the chemical composition of host plant tissues and the wasp's physiological ability to confront plant defensive chemistry could be another strong source of selection by limiting individual growth. In total, multifarious selection as a result of divergent host-plant use (temporal isolation, reduced immigrant viability and reduced fecundity) combined to generate strong RI between populations, resulting in a predicted total reduction of 92% in gene flow from Qv-populations to Qg-populations and a reduction of 81% in gene flow from Qg-populations to Qv-populations (Table S5). Altogether, the prezygotic barriers associated with different host plant use limit predicted gene flow between B. treatae populations to 1% to 5% in both directions (Table 2).

VARIATION OF THE STRENGTH OF RI
Due to the homogenizing effects of gene flow even when gene exchange is low (Slatkin 1987), temporal and/or spatial variation of the strength of RI can constrain divergence by providing a steppingstone for individuals to migrate between populations experiencing differing environments (Arnold 2016). It is therefore important to understand mechanisms underlying variation in RI. The two major sources of variation observed in the present study stem from differences among populations and varying levels of selection within each micro-environment (individual host plant genotype). Our study showed that both the direction and strength of RI varied across population pairs (Fig. 4). Under the scenario of Qv-associated B. treatae populations colonizing Qg host plants, the positive values of RI observed for population pair comparisons suggests reduced gene flow from the Qv to the Qg environment. However, the strength of combined RI varied from as high as 79% for some population comparisons to as low as 39% in others. Additionally, we found a wide range of variation in wasp performance across individual trees of the same host species (see Fig. 2), which is consistent with previous findings that individual trees varied in suitability for B. treatae wasps (Egan and Ott 2007), resulting in a patchy distribution of gall formers across space. Further research will compare the fitness differences between immigrants and residents while feeding on the same individual tree to explore the contribution of variation in microenvironments to the total variation of RI.

CROSS TAXA COMPARISON
We compared our results to those of 16 published studies that measured immigrant viability and fecundity in both directions using pairs of locally adapted populations (Porter and Benkman 2017). The first general observation we draw is that reduced immigrant viability and fecundity are commonly asymmetric and typically operate to reduce gene flow in a single direction (Fig. 5A, 5B), consistent with previous reviews of barriers being heterogeneous between recently diverged lineages (Lowry et al. 2008). Interestingly, reduced immigrant fecundity appeared to be less asymmetric than reduced immigrant viability (Fig. 5). However, using a paired sample t-test to compare study systems with measures of asymmetry for both viability and fecundity found that this difference is not clear (t = 1.473, P = 0.161). Thus, further work is required to thoroughly assess whether the degree of asymmetry is stronger for reduced immigrant viability than fecundity. A second important observation is that reduced immigrant viability and fecundity often act in opposing directions, suggesting that each barrier is free to evolve independently between populations, but may also complement each other to reduce gene flow in both directions and promote speciation. Third, comparison of the relative strengths of immigrant fitness differences associated with viability and fecundity in our study of B. treatae produced results that are in opposition to the overall pattern observed in Porter and Benkman (2017), where the reduction in gene flow due to immigrant inviability was roughly five times greater than reduction due to reduced immigrant fecundity (∼60% versus ∼12%). Finally, given our observation of asymmetric barriers, it is noteworthy that analyzing reproductive barriers in each environment separately may change the evaluation of the total importance of each reproductive barrier (i.e., how you summarize, or average barriers matters). For example, we re-analyzed the same 16 studies from Porter and Benkman 2017 plus our study on B. treatae considering separate environments and directionality in gene exchange and found that reduced immigrant viability and fecundity had no significant differences with respect to contribution to the overall strength of RI (Fig.  S3). By considering directionality, our findings differ from the general findings from Porter and Benkman (2017), which considered each barrier across all directions and observed a stronger effect of reduced immigrant fecundity compared to immigrant inviability. The difference in the conclusion, however, could be a result of the differences in the degree of asymmetry between different reproductive barriers. We observed that reduced immigrant fecundity tended to be less asymmetric than reduced immigrant viability, thus, moderately strong reproductive barriers in both environments could be more effective at reducing overall gene flow than a strong but asymmetric reproductive barrier. The published literature on asymmetric barriers and their importance, or relative role, in reducing gene flow and promoting speciation is sparse, therefore leaving an open question for future work.

Conclusions
Our study has assessed the relative contribution of reduced immigrant viability and immigrant fecundity to RI between ecologically divergent populations and how asymmetric and opposing reproductive barriers can combine to reduce gene flow. Reduced immigrant fecundity is rarely considered in the context of population divergence and ecological speciation (Nosil 2012). Our study provides an important contribution by considering this barrier together with the more commonly studied immigrant inviability and by assessing how both of these barriers contribute to prezygotic isolation in a single system. Moreover, our findings and those of similar studies highlight the importance of incorporating directionality in assessing RI, as asymmetrical reductions in fitness can influence estimates of the magnitude and relative importance of RI during the evolutionary process of speciation.

AUTHOR CONTRIBUTIONS
L.Z., G.R.H., J.R.O., 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.
Associate Editor: A. Rice Handling Editor: T. Chapman

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Location of field sites in this study. Table S2. Number of individual sapling oaks deployed at each study site and the number of replicates (bags per tree) conducted with the caged transplant study. Table S3. Number of asexual generation individuals used to estimate potential fecundity Table S4. Results of pairwise comparison from linear models, one-tail test Table S5. Temporal isolation, reduced immigrant viability and fecundity measured between different host-associated populations of B. treatae on each host environment for specific geographic locations where multiple barriers were measured from the same site. Figure S1. Boxplot depicting viability of resident (gray) and immigrant (white) B. treatae developing in each host environment subdivided by field and cage experiment. Figure S2. Predicted fecundity (least square means with 95% confidence intervals) of resident (gray) and immigrant (white) B. treatae developing in each host environment subdivided by field and cage experiment. Figure S3. Comparison of the strength of each reproductive barrier when considering directionality and different environments.