Parsing propagule pressure: Number, not size, of introductions drives colonization success in a novel environment

Abstract Predicting whether individuals will colonize a novel habitat is of fundamental ecological interest and is crucial to conservation efforts. A consistently supported predictor of colonization success is the number of individuals introduced, also called propagule pressure. Propagule pressure increases with the number of introductions and the number of individuals per introduction (the size of the introduction), but it is unresolved which process is a stronger driver of colonization success. Furthermore, their relative importance may depend upon the environment, with multiple introductions potentially enhancing colonization of fluctuating environments. To evaluate the relative importance of the number and size of introductions and its dependence upon environmental variability, we paired demographic simulations with a microcosm experiment. Using Tribolium flour beetles as a model system, we introduced a fixed number of individuals into replicated novel habitats of stable or fluctuating quality, varying the number of introductions through time and size of each introduction. We evaluated establishment probability and the size of extant populations through seven generations. We found that establishment probability generally increased with more, smaller introductions, but was not affected by biologically realistic fluctuations in environmental quality. Population size was not significantly affected by environmental variability in the simulations, but populations in the microcosms grew larger in a stable environment, especially with more introduction events. In general, the microcosm experiment yielded higher establishment probability and larger populations than the demographic simulations. We suggest that genetic mechanisms likely underlie these differences and thus deserve more attention in efforts to parse propagule pressure. Our results highlight the importance of preventing further introductions of undesirable species to invaded sites and suggest conservation efforts should focus on increasing the number of introductions or reintroductions of desirable species rather than increasing the size of those introduction events into harsh environments.

tive importance of the number and size of introductions and its dependence upon environmental variability, we paired demographic simulations with a microcosm experiment. Using Tribolium flour beetles as a model system, we introduced a fixed number of individuals into replicated novel habitats of stable or fluctuating quality, varying the number of introductions through time and size of each introduction. We evaluated establishment probability and the size of extant populations through seven generations. We found that establishment probability generally increased with more, smaller introductions, but was not affected by biologically realistic fluctuations in environmental quality. Population size was not significantly affected by environmental variability in the simulations, but populations in the microcosms grew larger in a stable environment, especially with more introduction events. In general, the microcosm experiment yielded higher establishment probability and larger populations than the demographic simulations. We suggest that genetic mechanisms likely underlie these differences and thus deserve more attention in efforts to parse propagule pressure.
Our results highlight the importance of preventing further introductions of undesirable species to invaded sites and suggest conservation efforts should focus on increasing the number of introductions or reintroductions of desirable species rather than increasing the size of those introduction events into harsh environments.

K E Y W O R D S
biocontrol, conservation, invasion, microcosm, population dynamics, propagule pressure, reintroduction, simulation, stochasticity

| INTRODUC TI ON
Colonization is the ecologically fundamental process of population establishment in an unoccupied location, and it underlies the past, present, and future distributions of species. Colonization occurs naturally but is increasingly prevalent due to anthropogenic influences (Cassey, Blackburn, Duncan, & Chown, 2005;Ricciardi, 2007;Sakai et al., 2001). Incipient populations often face environments that are entirely novel, which is especially likely in the case of anthropogenic colonization Ricciardi, 2007). Regardless of whether colonization events to novel habitats are natural (e.g., range expansion) or human-mediated (e.g., biological invasions, reintroductions of rare species, release of biological control agents), their successes or failures have significant implications for natural resource managers and society (Mack et al., 2000).
Propagule pressure is the total number of potentially reproductive individuals (e.g., adults, eggs, seeds, vegetative material) introduced to an area (Novak, 2007). It is often described in this broad sense, which belies its complexity. Two important components of propagule pressure are the number of introduction events-sometimes termed propagule number, and the average number of individuals per introduction event-sometimes termed propagule size (sensu Fauvergue et al., 2012). Here, we use "introduction regime" to refer to different combinations of the number of introduction events and the number of individuals introduced per event. The same total propagule pressure, N, is realized by different introduction regimes depending on how those N individuals are distributed in time or space (Haccou & Iwasa, 1996). An introduction regime of N individuals lies on a continuum bounded by maximizing the number of individuals introduced per event (all N individuals introduced in one event to the same location) and maximizing the number of introduction events (one individual introduced in each of N sequential events through time or to each of N unique locations).
Colonization success increases with total propagule pressure, but it is unclear whether the correlation is driven by the number of individuals introduced per event or the number of introduction events (Colautti et al., 2006;Lockwood et al., 2005;Simberloff, 2009). Historical data suggest that multiple introductions can facilitate population establishment, but conclusions from these studies are limited by the inability to control for the total number of individuals introduced (Blackburn & Duncan, 2001;Fauvergue et al., 2012;Grevstad, Coombs, & McEvoy, 2011;Hopper & Roush, 1993).
Models that hold the total propagule pressure constant agree that multiple, small introductions distributed across space will lead to greater establishment probability compared to a single, large introduction when Allee effects are weak (Grevstad, 1999;Haccou & Iwasa, 1996;Schreiber & Lloyd-Smith, 2009). However, both modeling and empirical approaches have generated conflicting views on how an introduction regime affects colonization success when introductions are distributed through time, a situation in which individuals from later introductions interact both demographically and genetically with individuals from previous introductions.
There is evidence from both models and experiments that colonization success can increase with more, smaller introduction events through time. Branching process models show that, in the long run, several small introductions will be more likely to successfully establish a population than a single large introduction (Haccou & Iwasa, 1996;Haccou & Vatutin, 2003). This finding was corroborated by simulations with no Allee effects (Drolet & Locke, 2016). In a Daphnia microcosm experiment, increasing introduction frequency (proportional to the number of introductions) positively affected population growth, but the number of individuals per introduction had no detectable effect (Drake, Baggenstos, & Lodge, 2005). Establishment probability was only affected by total propagule pressure and did not increase with increasing introduction frequency (Drake et al., 2005).
However, Drake et al. (2005) did not continue scheduled introductions if a population went extinct, denying those populations one of the main benefits of repeated introductions and creating variability in the total propagule pressure. In a more recent experiment, several small introductions through time led to a 65% increase in abundance of successfully colonizing invasive Pacific Oyster (Crassostrea gigas) compared to a single large introduction (Hedge, O'Connor, & Johnston, 2012).
There is also evidence from both models and experiments that colonization success can be greater with fewer larger introduction events than with multiple smaller introductions through time. For instance, in simulations of bird invasions by Cassey, Prowse, and Blackburn (2014), a single large introduction event always led to the greatest establishment probability. Drolet and Locke (2016) highlighted a key role of positive density dependence in favoring fewer larger introductions; when they included Allee effects in their simulations, colonization success was enhanced with fewer larger introductions-the reverse of their observed pattern without Allee effects. In a field experiment with the psyllid biocontrol agent, Arytainilla spartiophila, the number of individuals per introduction event was a better predictor of establishment success than the number of introduction events (Memmott, Craze, Harman, Syrett, & Fowler, 2005). In this case, however, the introduction regimes with the most individuals per event also had the highest total propagule pressure (Memmott et al., 2005). In an experiment that controlled total propagule pressure, a single, large introduction of the nonnative mysid, Hemimysis anomala, led to larger populations and greater survival probabilities compared to several, small introductions through time (Sinclair & Arnott, 2016).
Environmental stochasticity in the recipient environment may also affect which introduction regime is optimal for colonization.
Branching process models show that a more variable environment reduces the probability of population establishment for all introduction regimes (Haccou & Iwasa, 1996;Haccou & Vatutin, 2003), and simulations of introductions distributed in space suggest that greater environmental variability will magnify the benefit of multiple introductions (Grevstad, 1999). However, simulations of introductions through time by Cassey et al. (2014) did not support either of these outcomes-a single, large introduction was most likely to establish a population even with extreme levels of environmental stochasticity.

| General framework
In simulations and a microcosm experiment, we evaluated the outcome of introducing 20 total individuals to one of two environmental contexts (a stable or fluctuating novel environment), varying the number of introduction events used to distribute those individuals through time. This total propagule pressure was low enough to allow some population extinction within an observable timeframe, but high enough to be representative of documented introductions in the literature (Berggren, 2001;Drake et al., 2005;Grevstad, 1999;Simberloff, 1989Simberloff, , 2009Taylor, Jamieson, & Armstrong, 2005). The introduction regimes were: 20 individuals introduced in the first generation, 10 individuals introduced in each of the first two generations, five individuals introduced in each of the first four generations, and four individuals introduced in each of the first five generations.
To create the fluctuating environment, we imposed a magnitude of variability corresponding to environmental stochasticity in nature, which leads to frequent, mild-to-moderate perturbations in population growth rate due to external forces (Lande, 1993). Populations were tracked for seven generations following the initial introduction.
Establishment probability and the size of established populations were used as measures of colonization success. Populations were deemed "established" for any time step in which they were extant and population size was noted for all extant populations in every time step. Establishment probability and mean population size were assessed in generation 7. We also assessed establishment probability and population size 3 generations after the final introduction event (i.e., by assessing establishment and population size for the 20 × 1 regime at generation 3, 10 × 2 regime at generation 4, 4 × 5 regime at generation 6, and 5 × 4 regime at generation 7). By evaluating characteristics of these introduction scenarios at both an absolute time point (e.g., generation 7) and a relative time point (e.g.,

| Study system
Our simulations and microcosm experiment use Tribolium castaneum (red flour beetles) to model the life history of organisms with discrete, nonoverlapping generations (e.g., annual insects and fishes) following Melbourne and Hastings (2008). Individual simulated colonists were randomly sourced from a randomly mating, infinitely large population. Individual experimental colonists came from a thoroughly mixed source population maintained at 800 individuals in each of the four temporal blocks over which the experiment was replicated. The four source populations were themselves sourced from the "SF" laboratory colony which was collected in the wild and maintained at thousands of individuals for ~15 generations prior to use in this experiment (Hufbauer et al., 2015). To obtain colonists, beetles from the source population were mixed freely on a plate, selected from many sections of the plate, and introduced to subsets of the experiment populations in a random order (about 24 subsets per block). All migrant females were assumed to arrive mated from the source population.
The strong maternal effects exhibited by Tribolium flour beetles were reduced by rearing individuals from the source population on novel growth medium for one generation prior to using them as colonists (Hufbauer et al., 2015; Van Allen & Rudolf, 2013).

| Simulations
Population dynamics were simulated with a Negative Binomial-Binomial gamma (NBBg) model, which mechanistically describes change in population size in each discrete generation as a birthdeath process with cannibalism-induced density dependence and four distinct forms of stochasticity: environmental stochasticity, skewed sex ratio, demographic heterogeneity, and demographic stochasticity (Melbourne & Hastings, 2008). For each combination of the two introduction regimes and two recipient environment types (described in "General Framework" above), we simulated 500,000 replicate populations for seven generations. The NBBg model captures the stochastic and deterministic ecological dynamics of the Tribolium system well, but does not include evolutionary processes (Hufbauer et al., 2015;Melbourne & Hastings, 2008).
The core of the NBBg model is a Ricker function (Equation 1), where the mean size of population i at time t + 1 is a function of F mated(t,i) (the number of mated females in population i at time t), p (the probability of an individual being female, 0.5), R t,i (the densityindependent population growth rate for population i at time t), α (the egg cannibalism rate), and N t,i (the total size of population i at time t representing the sum of the residents and the migrants).
The expected equilibrium population size for the Ricker model, when N t+1 is expected to be equal to N t , is: Environmental stochasticity is treated as variability in the potential density-independent population growth rate arising from different unmeasured environmental conditions, and we model it Demographic heterogeneity is treated as variation in the expectation of the number of individuals in the next generation due to inherent differences in the egg-laying capacity of mated females (e.g., different sized females may lay different average numbers of eggs).
We model demographic heterogeneity (Equation 5) by drawing the expectation of population size for population i at time t + 1 from a gamma distribution with a mean of μ t+1,i (the expected mean size of population i at time t + 1 derived from the Ricker function which only incorporates environmental stochasticity and skewed sex ratio) and a shape parameter k D *F mated(t,i) (where a smaller value of k D indicates greater demographic heterogeneity, and the dependence of the shape parameter on F mated(t,i) indicates greater demographic heterogeneity with fewer mated females at time t).
We model variation arising from demographic stochasticity To estimate the key parameters for the simulation model (α, R 0 , k E , and k D ), we censused 125 Tribolium populations one generation after establishing them at various, known densities (between 5 and 200 individuals) on the novel growth medium following the rearing procedure described in "Microcosm Experiment" below. We combined Equations 1, 3, 4, 7, and 8 into a single hierarchical model with weakly regularizing gamma priors taken from Melbourne and Hastings (2008) (Equation 9) and fit the model to the population size We used the posterior distribution of the estimated environmental stochasticity parameter, k E , for simulations in a stable environment because the 125 populations were reared on a single novel growth medium mixture (mixture 5, Supporting information Table S1). We parameterized a fluctuating environment simulation by increasing the standard deviation of the density-independent per capita population growth rate by 10% compared to the stable environment for any given population at any given time step, a similar magnitude of variation as imposed by Cassey et al. (2014) (see Supporting information Methods S1). We used random sets

| Microcosm experiment
We    Table S1). We chose this range to mimic environmental stochasticity measured in nature (Saether & Engen, 2002) while remaining within the bounds of biologically realistic population growth (λ between 0.5 and 1.5) (Morris et al., 2008). We chose random sequences of growth media such that the expected geometric mean population growth rate for the population experiencing that sequence resembled expected growth of populations in the stable environment (λ expected = 1.2 ± 0.05).
We estimated the amount of environmental stochasticity that we achieved in the fluctuating environment as the difference in mean total stochasticity between populations in the fluctuating and stable environments (Saether & Engen, 2002). We assumed that total stochasticity was a combination of demographic and environmental stochasticity for populations in the fluctuating environment, and that demographic stochasticity was the sole contributor to total stochasticity for populations in the stable environment. Total stochasticity (demographic plus environmental) of each population that did not experience extinction (n = 667) was calculated as the variance of the natural logarithms of its population growth rates through seven generations following Saether and Engen (2002): where, for a particular population, s total is its total stochasticity, s demographic is its demographic stochasticity, s environmental is its environmental stochasticity (assumed to be 0 for populations in the stable environment), and λ t is its per capita population growth rate between generation t−1 and generation t (t = 1, 2, …, 7). We only calculated total stochasticity for populations that did not experience any extinction in order to capture the full temporal extent of environmental fluctuations and because extinctions would have an infinite effect on this measure of stochasticity.

| Statistical analyses
We evaluated how our environment treatment affected variability in population growth rate (total stochasticity from Equation 10) using a linear mixed effects model with environment (stable or fluctuating) as a fixed effect and block as a random intercept effect.
We used a mixed effects logistic regression with a logit link to predict the binary response of establishment, and a mixed effects Poisson regression with a log link to analyze population size. In both models, introduction regime, environment treatment, and their interaction were treated as fixed effects, and block was treated as a random intercept effect.
We assessed the effect of temporary extinctions on the establishment probability and mean population size by fitting the generalized linear mixed effects models described above to data from the multiple introduction regimes (i.e., not the 20 × 1 regime) and with additional predictor variables. To assess the effect of the presence of a temporary extinction, we included an additional Boolean predictor for whether a population went temporarily extinct or not. To assess the effect of total propagule pressure, we used the number of beetles introduced after the latest temporary extinction as a predictor because only introductions after the latest temporary extinction contribute to total propagule pressure.
Group-level significance of fixed effects was tested using likelihood ratio tests on nested models, and least-squares contrasts were used to compare levels of the fixed effects. All statistical analyses were performed in R, version 3.3.2 (R Core Team 2017). Generalized linear mixed models were fit using the lme4 package, version 1.1-12 (Bates, Machler, Bolker, & Walker, 2015) and pairwise comparisons were made using the lsmeans package, version 2.25 (Lenth, 2016).

| Data availability
The raw experiment data, simulated population trajectories, R code

| Note on egg contamination
Laboratory procedures after generation 3 resulted in occasional egg contamination between replicate populations of the same introduction regime/environmental variability treatment. To estimate the extent and magnitude of contamination, we examined trajectories of populations having only one individual and no additional introduction events, which should have deterministically gone extinct.
Of these 49 populations, 12 did not go extinct (24.5%) and instead persisted at a small size (min = 1, max = 3, mean = 1.9). For colonization success analyses, we manually edited population trajectories such that they went extinct in the generation after having only one individual.

| Simulations
Summary statistics for the posterior distributions of the four NBBg model parameters are given in Supporting information Table S2.
Posterior and prior distributions are shown in Supporting information Figure S1.
Our simulations showed introduction regimes with more introduction events were more likely to establish a population by the seventh generation (Figure 1). Introduction regimes with fewer introduction events were more likely to establish populations by three generations after the final introduction event. Simulated introductions into a stochastically fluctuating environment resulted in slightly lower population establishment for all introduction regimes (difference of ~1%), but did not favor any particular regime (not shown). The mean population sizes for each introduction regime/ environment combination were approximately equal in simulations.
Mean population sizes only incorporated extant populations, so they were slightly larger than the expectation for the equilibrium population size (Figure 2).
This value is near the median value of 0.055 measured in nature by Saether and Engen (2002) in a meta-analysis of 35 avian populations, indicating that we achieved biologically realistic fluctuations in population growth rate.
We found no evidence that the probability of establishment was affected by a main effect of environment (χ 2 = 0.72, df = 1, p = 0.40 for generation 7 assessment; χ 2 = 0.25, df = 1, p = 0.62 for assessment three generations after final introduction), nor by an interaction between environment and introduction (10) s total = s demographic + s environmental = var( log (λ t )) regime when establishment was assessed at generation 7 (χ 2 = 3.49, df = 3, p = 0.32). We detected a significant effect of an introduction regime/environment interaction when assessing establishment probability three generations after the final introduction event (χ 2 = 16.61, df = 3, p = 0.0008). There was strong support for an effect of introduction regime on establishment probability (χ 2 = 59.76, df = 3, p < 0.0001 for generation 7 assessment; χ 2 = 17.52, df = 3, p = 0.0006 for assessment three generations after final introduction). Pairwise comparisons of the different introduction regimes averaged across the environment treatments revealed that the 4 × 5 regime was the most likely to establish populations by generation 7, with a probability of about 0.98, whereas the 20 × 1 and 10 × 2 regimes were the least likely to establish populations, with a probability reduced to about 0.8 (Figure 1 left panel, Figure 3). A similar pattern emerged for establishment three generations after the final introduction, with the 4 × 5 regime being the most likely to establish a population with a probability of 99% (although not statistically distinguishable from the 20 × 1 or 5 × 4 regimes) and the 10 × 2 regime being the least likely to establish F I G U R E 1 Establishment probabilities assessed seven generations after the first introduction event (left panel) and three generations after the final introduction event (right panel) for each introduction regime. Because different introduction regimes took different numbers of generations to complete, assessments made three generations after the final introduction event were made in different absolute experimental generations (i.e., 20 × 1 assessment made at generation 3, 10 × 2 assessment made at generation 4, 4 × 5 assessment made at generation 6, and 5 × 4 assessment made at generation 7). Triangles represent results from simulations. Dot-whiskers represent estimates and 95% confidence intervals from the mixed effects logistic regression model fit to data from the microcosm experiment. Different letters over dot-whiskers represent introduction regime/environment combinations with significantly different establishment probabilities in the microcosm experiment and three generations after the final introduction event (right panel) for each introduction regime. Because different introduction regimes took different numbers of generations to complete, assessments made three generations after the final introduction event were made in different absolute experimental generations (i.e., 20 × 1 assessment made at generation 3, 10 × 2 assessment made at generation 4, 4 × 5 assessment made at generation 6, and 5 × 4 assessment made at generation 7). Triangles represent results from simulations. Dotwhiskers represent estimates and 95% confidence intervals from the mixed effects Poisson regression model fit to data from the microcosm experiment. Different letters over dot-whiskers represent introduction regimes with significantly different mean population sizes. The dashed line represents the theoretical equilibrium population size derived using Equation 2, which includes extinction and is therefore lower than the mean population size from the simulations For assessments made at generation 7 and 3 generations after the final introduction event, populations established via more introduction events were generally larger when averaged across the environment treatments. Extant populations in the stable environment were larger than those in the fluctuating environment when averaged across the introduction regimes. The interaction manifests as the benefit of a stable environment increases with more, smaller introduction events (Figures 2 and 4).  (Lande, 1993).
Population dynamics in demographic simulations were likely to have been affected by compensatory density dependence. Tribolium beetles in challenging environments experience negative density dependence arising from egg cannibalism, which is more common in harsh environments (Via, 1999). From our NBBg Ricker model fit, we estimated the egg cannibalism rate, α, to be relatively high compared to previous estimates on the more benign natal growth medium (our estimated mean α was 0.0087 while the mean of our prior taken from Melbourne and Hastings (2008) was 0.0037; see Supporting information Figure S1). This high cannibalism rate was incorporated as a key dynamic in our simulations and may have overwhelmed any effect of environmental variability. Furthermore, the strong negative density dependence in the simulations likely led to the convergence of the mean population sizes by generation 7 and 3 generations after final introduction events across all treatments. A similar convergence in population size was not observed in the microcosm, but oscillating population sizes were observed for all treatments (Figure 4) suggesting some role for compensatory density dependence in the experiment as well.
Density dependence is likely to interact with introduction regime to affect colonization success. For instance, fewer larger introductions are more likely to establish populations when positive density dependence is present (Drolet & Locke, 2016;Grevstad, 1999). Wittmann, Metzler, Gabriel, and Jeschke (2014) further show an effect of per capita population growth in systems with negative density dependence: when population growth rate is consistently greater than one in those systems, population establishment is faster with several smaller introductions, but when population growth rate is mixed (sometimes greater than one and sometimes less than one depending on population size), population establishment occurs fastest with fewer larger introductions. Strong negative density dependence led to an expected equilibrium population size below the total propagule pressure in this experiment (13.7 individuals; Equation 2; Figure 2). Thus, our simulated populations experienced the "mixed" scenario described by Wittmann et al. (2014), and predictably benefitted from fewer, larger introductions when standardizing our assessment of establishment probability to three generations after the final introduction event.
Population growth rate may also interact with introduction regime to affect colonization success. Cassey et al. (2014) suggested that their lower simulated mean population growth rates, where R was between 1.0 and 1.38, compared to that of Grevstad (1999), where R was equal to 2.0, explained why they found that a single large introduction always led to a greater establishment probability, while Grevstad (1999) found several small introductions to be more successful. However, our expected mean population growth rate was relatively low (R = 1.132; Supporting information Table S2) and we still found that several small introductions had the greatest colonization success at generation 7 (left panels of Figures 1 and 2). This the novel, harsh environment from standing variation is possible in this species within a similar timeframe as our experiment (Hufbauer et al., 2015;Szűcs, Vahsen, et al., 2017), and likely explains the greater establishment probability and population sizes in the microcosm compared to expectations derived from demographic simulations which do not include adaptation.
The rescue effect of multiple introductions can also act genetically by increasing the fitness of populations (Frankham, 2015;Hufbauer et al., 2015;Whiteley, Fitzpatrick, ChrisFunk, & Tallmon, 2015). The experimental immigrants all came from the same source population, so it is unlikely that gene flow from migrants united previously separated alleles into high-fitness genotypes (Novak, 2007). Thus, a more likely mechanism by which immigration increased mean population size beyond expectations was by relieving inbreeding depression or counteracting driftinduced allele loss. Small populations are more prone to experiencing increased homozygosity and inbreeding depression, which can reduce population growth rates and increase extinction risk (McCauley & Wade, 1981;O'Grady et al., 2006;Szűcs, Melbourne, Tuff, Weiss-Lehman, & Hufbauer, 2017). Even small amounts of gene flow can alleviate these effects, so the additional small introductions of mated individuals from the external source population were well-suited to bring about longer-term relief (Hufbauer et al., 2015;Slatkin, 1985). However, Cassey et al. (2014) found that simulated inbreeding depression was especially detrimental for several, small introductions through time, so other mechanisms are likely at play.
One such mechanism may be adaptation of the incipient population, which is affected by sustained immigration. Introductions to a harsh, novel habitat can result in adaptive evolution with the right amount of gene flow if additional immigrants to a declining population prevent extinction long enough to allow for adaptation to occur (Holt & Gomulkiewicz, 1997). The strong compensatory density dependence as a result of egg cannibalism may ultimately provide a pathway for adaptation, as cannibalism rates are genetically variable and can confer individual fitness advantage via reduced development time (Via, 1999). However, too much gene flow can lead to genetic swamping whereby the homogenizing effects of gene flow overpowers ongoing local adaptation (Lenormand, 2002). This may have been the case for populations in the 10 × 2 introduction regime, which had the highest average migration rate, lowest establishment probability, and lowest mean population size. Alternatively, negative density dependence may have reduced population fitness when migration rates were high, reducing population growth rates and hampering the spread of adaptive alleles (Holt & Gomulkiewicz, 1997). Although not significant, the lower establishment probability and mean population size in the 10 × 2 introduction regime warrants further work to assess whether they exemplify the yet-unseen scenario described by Blackburn, Lockwood, and Cassey (2015) in which maladaptive gene flow from multiple introductions hampers population establishment in a novel range. More broadly, further study is necessary to evaluate how immigration affected adaptation in this system, if at all (Boulding & Hay, 2001;Gomulkiewicz & Holt, 1995).

| CON CLUS ION
Our experimental results suggest that several, small introductions through time lead to greater colonization success in a novel habitat, and that introductions into a stable recipient environment lead to larger population sizes, but not greater establishment probability. Furthermore, introductions to a stable recipient environment are especially beneficial to populations established with more introduction events. These results defied our expectations derived from parallel simulations of a model that included demographic processes but not evolutionary ones, so we suspect a genetic mechanism might be at work. Genetic mechanisms are rarely incorporated when simulating the effect of introduction regime on colonization (but see Cassey et al., 2014), and our multigeneration microcosm is unique in bringing evolutionary processes to bear on parsing two key components of propagule pressure in an experimental setting. Sustained introduction efforts should also bring about concomitant benefits in the form of longer-term monitoring, increased data collection, and more opportunities for experimentation and adaptive management (Godefroid et al., 2011;Lockwood et al., 2005). Latimer, Truman Young, and Jenny Gremer for comments on the project and manuscript. We also thank two anonymous reviewers for their careful feedback, which greatly improved the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
MJK, RAH, and MFO designed the experiment. MJK, MFO, and RAH collected the data. MJK and BAM designed the simulations. MJK implemented the simulations. MJK, RAH, MFO, and BAM wrote the paper.