Eight collections of D. carinulata were made from across the introduced range of the species. Four collections were from well-established original introduction sites in the north of the range and represent the range core. Another four collections were from the newly established southern edge of the range (Fig. 2, Table S1). One edge population (H in Fig. 2 and Table S1) was collected slightly behind the probable edge in 2017 in order to avoid other Diorhabda species that were moving northward and overlapping the D. carinulata range (Ozsoy et al., 2019). Two populations (one core and one edge) were collected Fall 2017 and the first lab generations were put into reproductive diapause to reduce the number of lab generations before the start of the experiment. All other populations were collected in Summer 2018 and cultured in the lab for one generation to standardize maternal environment effects prior to starting experiments. All insects were reared in growth chambers under reproductive conditions (D. W. Bean, Wang, et al., 2007), with a 16h/8h light/dark cycle and 25°C/20°C day/night temperatures and were fed fresh tamarisk as needed.
Life-history
Newly emerged adult females of the second lab generation were weighed before feeding and reared individually thereafter in 0.24 L plastic containers with mesh lids. Three days after eclosion, each female was paired with a male from the same population of about the same age and allowed to mate for 24 hours before the male was removed. Presence of eggs in each container was assessed daily. All eggs that were laid on the first day of oviposition were counted to provide a measure of early fecundity and the age at first reproduction in days was recorded. Individuals that had not oviposited within seven days of emergence were recorded as non-layers (Lewis et al., 2003). We collected mass and egg count data from 130 core and 140 edge female D. carinulata, and age at first reproduction from 104 core and 118 edge females.
Dispersal
Dispersal ability was measured for only male D. carinulata since, in the field, males have been observed dispersing first and using pheromones to attract mixed-sex aggregations of reproductive adults (Cossé et al., 2005). After emergence as adults, males were randomly assigned to density and mating treatments. Males assigned to the mated treatment were paired with a female from the same population for 24 hours. The males were thereafter reared in 0.24 L plastic containers with mesh lids in groups of five (high density) or alone (low density). All males in each high-density container were of the same mating treatment. All containers received the same surplus amount of fresh tamarisk, regardless of how many beetles were in the container. Male beetles were between 6 and 23 days after eclosion and were weighed on the day of the dispersal trial.
We assessed dispersal of male beetles using tethered flight mills (reviewed in Minter et al., 2018), similar in design to Maes et al. (2014) (Appendix S1). Each beetle was given one hour to take any number of flights on a flight mill. Each individual was between 6 and 23 days post-eclosion for its flight trial. Data from each trial was converted into four dispersal elements: occurrence of at least one flight, number of flights, total flight distance, and average flight speed (Appendix S2). Each dispersal element has different biological relevance (Stevens et al., 2013; Tung et al., 2017). Occurrence of flight and number of flights measure the probability and frequency of movement from the local patch. Total distance and average speed measure how far individuals disperse, after initiation. Spatial sorting may act on any one or combination of these elements. We collected dispersal data from 279 core males and 311 edge males, with 65 to 81 males in each density-mating treatment combination and at least 15 from each population. Average speed was calculated for 231 core and 266 edge males that took at least one flight, with 56 to 71 males in each density-mating treatment combination.
Statistical Analyses
The three life-history traits of female body mass, fecundity over 24 hours, and age at first reproduction were analyzed both individually and with a multivariate analysis of variance (MANOVA). In the MANOVA, only individuals that produced eggs during the experiment were included, since they were the only complete observations. The three life-history traits were the response variables, and range (core or edge), population, and eclosion date were fixed effects.
Mass of female beetles at adult emergence was analyzed individually with a linear mixed model, with range as a fixed effect and population of origin as a random effect.
Since some females (17%) did not lay eggs within 10 days, the data on the numbers of eggs were split into two datasets, one including the number of eggs for laying individuals, and the other including the binary response (laying or non-laying) for all individuals, to assess whether probability of laying and fecundity differed between core and edge populations. To account for overdispersion, a negative-binomial mixed model was fit to the count data (excluding non-laying individuals) using the glmmTMB package version 1.0.2.1 (Brooks et al., 2017) with range, mass, and age at first reproduction as fixed effects and population as a random effect. A logistic mixed model was fit with the glmmTMB package to the binary dataset with range and mass as fixed effects and population as a random effect. Age at first reproduction could not be included as a covariate due to convergence issues.
Age at first reproduction was analyzed for individuals that reproduced during the experiment using a Conway-Maxwell Poisson mixed model which accounts for under-dispersion in this dataset (Brooks et al., 2019), with range and mass as fixed effects and population as a random effect. We complemented this analysis with a Kaplan-Meier survival analysis, which accounts for censoring of individuals that either did not reproduce during the experiment or died before reproducing. Since covariates and random effects cannot be added in Kaplan-Meier survival analyses, the analysis was run twice, once with population as the predictor to visualize the spread among populations and again with range as the predictor to estimate the total effect of range.
Mass of male beetles at the time of the dispersal trial was analyzed with a linear mixed model, with range, rearing density (low or high), mating status (unmated or mated), and interactions between those factors as fixed effects, age at time of weighing as a fixed covariate, and population as a random effect.
Each of the four dispersal elements (occurrence of flight, number of flights, total distance, and average speed) were analyzed separately. We chose not to do multivariate analyses on the dispersal traits since all but one variable was highly zero-inflated and skewed, which violates assumptions of multivariate tests such as MANOVA, and univariate tests could better incorporate sampling design and environmental covariates during the trials, which greatly improve model fit. For each dispersal element, the same factors were included. Range, density, mating status, and all interactions were fixed effects. Mass at the time of dispersal trial, age, mill friction (Appendix S1), and temperature were fixed covariates. Population and trial date were random effects.
The occurrence of at least one flight was analyzed with a binomial model with the packages lme4 version 1.1-26 (Bates et al., 2015) and lmerTest version 3.1-3 (Kuznetsova et al., 2017). The number of flights during the 1-hour trail was analyzed using both a negative binomial and a Poisson mixed model with the glmmTMB package. QQ plots of model residuals and residual vs. fitted plots from the DHARMa package version 0.3.3.0 were used to assess model fit (Hartig, 2020).The negative binomial mixed model best met assumptions of normality of residuals and was chosen as the final model. Total distance was a count of revolutions of the flight mill and thus a discrete variable. Six models were fit for total distance: a linear mixed model on log-transformed distance, a negative binomial mixed model, a generalized Poisson model (commonly used for highly right skewed data with a high frequency of low counts (Brooks et al., 2019; Joe & Zhu, 2005)), and a zero-inflated version of each of those, using the glmmTMB package. Based on the same residual diagnostics as above, the zero-inflated generalized Poisson model best met assumptions of normality of residuals and was chosen as the final model. Average speed during the 1-hr trial included only trials in which the beetle made at least one flight and was analyzed with a linear mixed model using lme4 and lmerTest.
In all dispersal models, the three-way interaction was dropped from the model if it was not statistically significant. Post hoc comparison of means was done with the emmeans package version 1.5.4 (Lenth, 2020). All analyses were done in R version 3.6.2 (R Core Team, 2018).