Geographic potential of the world’s largest hornet, Vespa mandarinia Smith (Hymenoptera: Vespidae), worldwide and particularly in North America

Occurrence and environmental data

We downloaded occurrence data for V. mandarinia from the Global Biodiversity Information Facility database (GBIF.org, 2020; https://www.gbif.org/). We kept records from the species’ native range (Fig. 1) separate from non-native occurrences facilitated by human introduction. We cleaned occurrences from the native distribution following Cobos et al. (2018) by removing duplicates and records with inconsistent georeferencing (coordinates outside country limits, on the sea, or missing, as recommended in the literature of data cleaning; Chapman, 2005). To avoid model fitting influences of spatial autocorrelation and overdominance of specific regions due to sampling bias, we thinned these records spatially in two ways: by geographic distance and by density of records per country (Fig. 2). In the first case (distance-based thinning; Anderson, 2012), we excluded occurrences that were <50 km away from another occurrence record. In the second thinning approach (country-density thinning), we accounted for potential differences in survey and reporting activities among countries and randomly reduced numbers of occurrences in countries with the densest sampling, namely Japan, Taiwan, and South Korea (from 30, 6, and 5, to 6, 2, and 2 occurrences, respectively), to match the approximate reference density of India, Nepal, and China. We used the package ellipsenm (Cobos et al., 2020; available at https://github.com/marlonecobos/ellipsenm) in R 3.6.2 (R Core Team, 2019) to clean and thin the data. We retained 172 occurrence records for V. mandarinia after initial data cleaning, 49 records after the distance-based thinning approach, and 18 records after the country-density thinning approach (Fig. 1). We then treated both data sets independently in all subsequent analysis steps.

For environmental predictors, we used bioclimatic variables at 10′ resolution (~18 km at the Equator) from the MERRAclim database (Vega, Pertierra & Olalla-Tárraga, 2018). We excluded four variables because they are known to contain spatial artifacts as a result of combining temperature and humidity information (Escobar et al., 2014): mean temperature of most humid quarter, mean temperature of least humid quarter, specific humidity mean of warmest quarter, and specific humidity mean of coldest quarter. The 15 variables remaining were masked to an area for model calibration (M, see Ecological niche modeling).

These 15 variables were further processed into two subsets, each used separately in all subsequent analyses. One set consisted of submitting the variables to a principal components analysis (PCA) to reduce dimensionality and multicollinearity. The other set was created with the raw environmental variables that had a Pearson’s correlation coefficient ≤0.85 measured in the calibration area, choosing the most biologically relevant or interpretable variables based on our knowledge of AGH natural history (Simões et al., 2020). As a result of this selection process, we obtained six raw variables: isothermality (BIO3), maximum temperature of warmest month (BIO5), minimum temperature of coldest month (BIO6), temperature annual range (BIO7), specific humidity of most humid month (BIO13), and specific humidity of least humid month (BIO14). The PCA was performed with the raw variables masked to the M area, and principal components for the entire world were obtained by transforming raw variables (at world extent) using the scaling and rotations from the PCA obtained for M. Projecting the PCA results from the accessible area to the entire world ensures that values are comparable and prevents further problems when models are projected. For further analyses, we kept the first four PC axes, as they explained 97.9% of the cumulative variance (Fig. S1). All analyses were done in R; specifically, raster processing was completed using the packages raster (Hijmans et al., 2020), rgeos (Bivand et al., 2020b), and rgdal (Bivand et al., 2020a); PCA was performed using the ntbox package (Osorio-Olvera et al., 2020).

Ecological niche modeling

To identify a calibration area (ostensibly equivalent to M; Owens et al., 2013) for our models, we considered a region contained within a buffer of 500 km around the known occurrence records after the 50 km thinning process (Fig. 1). This buffer distance was selected considering that little is known about the species’ dispersal ability (Matsuura & Sakagami, 1973; APHIS, 2020), as no information exists about queens’ dispersal capacities; however, we consider that given the species’ body size, long-distance dispersal events should not be discarded. We used all pixels in M (15,411) as the background across which to calibrate the models.

Given uncertainty deriving from specific treatments of occurrence records and environmental predictors in ecological niche modeling (Alkishe et al., 2020), we calibrated models via four distinct schemes: (1) using raw variables and distance-based thinned occurrences, (2) using PCs and distance-based thinned occurrences, (3) using raw environmental variables and country-density thinned occurrences, and (4) using PCs and country-density thinned occurrences (Fig. 2). For each scheme, we calibrated models five times, each time randomly selecting 50% of the occurrences for calibrating models (random k-fold evaluation, where k = 5), and using the remaining records for testing (Cobos et al., 2019a).

Each process of model calibration consisted of creating and evaluating candidate models using Maxent (Phillips, Anderson & Schapire, 2006; Phillips et al., 2017). Since many choices are needed to parameterize Maxent (Merow, Smith & Silander, 2013), we chose distinct parameter settings: 10 regularization multiplier values (0.10, 0.25, 0.50, 0.75, 1, 2, 3, 4, 5, 6), eight feature classes (lq, lp, lqp, qp, q, lqpt, lqpth, lqph, where l is linear, q is quadratic, p is product, t is threshold, and h is hinge), and all combinations of more than two predictor variables (Cobos et al., 2019b; Tables S1 and S2). This resulted in a total of 4,560 models using raw variables and 880 using PCs that were tested, in tandem with the two subsets of occurrence data described above. We assessed model performance using partial ROC (for statistical significance; Peterson, Papeş & Soberón, 2008), omission rates (E = 5%, for predictive ability; Anderson, Lew & Peterson, 2003), and Akaike Information Criterion corrected for small sample sizes (AICc, for model complexity; Warren & Seifert, 2011). We selected models with delta AICc ≤2 (Cobos et al., 2019a) from those that were statistically significant and had omission rates below 5%.

After model calibration, we created models with the selected parameter values, using all occurrences after the corresponding thinning process, with 10 bootstrap replicates, cloglog output (Phillips et al., 2017), and model transfers using three types of extrapolation (free extrapolation, extrapolation with clamping, no extrapolation; Owens et al., 2013). Not all calibration processes identified models that met all three criteria of model selection; we did not consider those models in further analyses (Fig. 2; Table 1). As a final evaluation step, we binarized all model replicates given a modified least presence (5% of omission) and then tested whether each replicate of the selected models was able to anticipate the known invasive records of the species in British Columbia, Canada and Washington, USA. For each scheme, using only those model replicates that met the selection criteria and correctly predicted independent occurrences (known invaded localities in North America; independent testing), we created a consensus per sample and two types of final consensus: (1) a median of the medians obtained for each parameterization (continuous), and (2) the sum of all suitable areas derived from binarizing each replicate using a modified least presence (5% omission) threshold (this represents the number of coincidences; Pearson et al., 2007; Fig. 2). Our exploration of different modeling pipelines allowed us to highlight how different methodologies produce different results.

As we transferred models to the entire world, we used the mobility-oriented parity metric (MOP; Owens et al., 2013) to detect areas where strict or combinational extrapolation risks could be expected, given the presence of non-analogous conditions with respect to the environments manifested across the calibration area. The areas where extrapolation risks were detected using MOP were deleted from our binary results (suitable areas) to avoid potentially problematic interpretations based on extrapolative situations. Model calibration, production of selected models with replicates, and MOP analyses were done in R using the package kuenm (Cobos et al., 2019a); raster processing and independent testing of models were done using the package raster.

Dispersal simulations

We used the binary outputs from the final consensus models (suitable and unsuitable areas, without areas of strict extrapolation) to simulate invasion dynamics of the AGH. All simulations were started from the Pacific Northwest, from sites already known to be occupied by the AGH. The simulations were performed using the cellular automaton dynamic model included in the bam R package (Osorio-Olvera & Soberón, 2020; available at https://github.com/luismurao/bam). Under this discrete model, given an occupied area at time t, two layers of information are needed to obtain the occupied area at time t + 1: (i) the binary layer of suitability for the species, and (ii) a connectivity matrix determined by the species’ ability to reach neighboring cells in one time unit (known as “Moore’s neighborhood”; Gray et al., 2003, that defines patches that are connected by dispersal). At each step, each of the suitable cells can be either occupied or not by the species. If a cell is occupied, adjacent cells can be visited by the species, and if suitable, they become occupied. This method is similar to the one implemented in the MigClim R package (Engler, Hordijk & Guisan, 2012), but uses a simpler dispersal kernel and parameterization.

With each of the final consensus models for V. mandarinia, we performed a set of simulations in which we explored different degrees of connectivity (1, 2, 4, 8, 10, and 12 neighbor cells pixels per unit time) and different suitability thresholds (10 equidistant levels from 3% to 10% of the presence points to explore variability in sensitivity to the amount of area classified as suitable) to create the binary maps. Since no information is available about dispersal capacities of queens of V. mandarinia, all simulations were done with 200 arbitrary time steps that ensure reaching a steady state. In the end, we visualized the simulation results by summing the occupied distribution layers obtained from each set of simulations. A value of 100 in these final layers means that the species reached that cell in 100% of the simulations, whereas a value of 0 means that the species never reached that cell. Further details regarding the simulation processes can be found in the Supplemental Information.

Honey production and native bee richness in North America

To explore potential ecological and economic impacts of the invasion of the AGH in North America, we explored annual, state-level production of honey (for Mexico, Untied States, and Canada) as well as species richness of bumble bees (Bombus Latreille) and stingless bees of the genus Melipona Illiger in Mexico and the United States. We extracted data on 2016 honey production (in US dollars) for the United States from the U.S. Department of Agriculture (USDA; available at https://quickstats.nass.usda.gov/#4A0314DA-F3E5-3B06-ADD1-CA8032FBD937), from the Instituto Nacional de Estadística, Geografía e Informática (INEGI) for Mexico (https://atlasnacionaldelasabejasmx.github.io/atlas/cap5.html), and from the government of Canada website (https://www.agr.gc.ca/eng/horticulture/horticulture-sector-reports/statistical-overview-of-the-canadian-honey-and-bee-industry-2018/?id=1571143699779) for Canada. For native species richness, we obtained a list of species of bumble bees and stingless bees of the genus Melipona that occur in Mexico and the United States from Discover Life (https://www.discoverlife.org/) and downloaded their occurrence data from GBIF.org (2020). We chose these bee taxa as likely targets of AGH because the species in these groups are of similar body size and behavior to the typical prey of these hornets: they are social insects that form annual or perennial colonies that can have a few hundreds to as many as 10,000 individuals (Cueva del Castillo, Sanabria-Urbán & Serrano-Meneses, 2015; Viana et al., 2015), and store honey and pollen inside their nests (Michener, 2000). To summarize species richness of these two genera, we created a presence-absence matrix (PAM; Arita et al., 2008) for North America, based on geographic coordinates of occurrence data, with a pixel size of one degree. The PAM was created in R with the package biosurvey (Nuñez-Penichet et al., 2020; available at https://github.com/claununez/biosurvey).

To assure transparency and reproducibility of our work, we include an Overview, Data, Model, Assessment, and Prediction protocol (ODMAP; Zurell et al., 2020) in our Supplemental Materials. This metadata summary provides a detailed key to the steps of our analyses. The data and R code used in this research are openly available at http://hdl.handle.net/1808/30602 and https://github.com/townpeterson/vespa repositories, respectively.

Model calibration

The number of models that met the selection criteria was considerably smaller than the total number of models tested (Table 1, see Supplemental Information for more details of the calibration results). The calibration schemes including raw variables had fewer models selected than those using PCs (11, 19, 6, 15 models selected for raw/distance-thinned, PC/distance-thinned, raw/country-density, and PC/country-density, respectively). Not all replicates of selected models anticipated the V. mandarinia invaded areas in North America successfully, so we kept only those that predicted all known invasive records. The number of replicates retained varied among distinct calibration schemes and types of extrapolation used (Table 1).

Ecological niche model predictions

In our models, areas predicted as suitable for the AGH varied among calibration schemes, in both extension and geographic pattern (Fig. 3; Figs. S2–S4). The differences are conspicuous between the two types of thinning approaches, which resulted in models created with different numbers of occurrence records. Models with country-density thinning (18 records) resulted in broad predicted suitable areas worldwide, with areas of higher values of suitability concentrated in tropical regions (Fig. 3; Figs. S2–S4). In contrast, models created with the greater number of occurrences (49 records) from the geographic distance thinning predicted more patches of suitable areas across large extensions of Southeast Asia, Europe, West Africa, Central America, northern South America, and the Pacific Northwest and southeastern United States (Fig. 3; Figs. S2–S4). In the calibration area, the areas detected with high levels of suitability were larger in the scheme with geographic distance thinned occurrences and the raw variables and smaller in the predictions obtained with the country-density thinned occurrences and the PCs as environmental predictors (Fig. 3). In all schemes, the two northernmost occurrence points of this species in China were accorded relatively low levels of suitability (Fig. 3). Predicted suitable areas for this hornet worldwide were also different among types of extrapolation considered in this study, especially as regards its distribution size rather than location (Figs. S2–S4).

In North America, across multiple model calibration schemes, our various models agreed in predicting suitable areas for AGH in the Pacific region of southwestern Canada, the Pacific Northwest, the southeastern United States, and from central Mexico south to southernmost Panama (Fig. 4). Our model calibration schemes also agreed in identifying the Rocky Mountains and Great Plains as unsuitable for this species (Fig. 4).

The proportion of area identified as suitable varied among the data thinning schemes. In the case of models created with raw variables, the proportion was 0.171 and 0.164 for spatially thinned and country-density thinned records, respectively. When PCs were used, suitable proportions were 0.248 and 0.239, for spatially thinned and country-density thinned records, respectively (Table S3).

Extrapolation risks in model projections

The pattern of areas detected with risk of extrapolation was similar worldwide between thinning methods, but different between raw variables and PCs (Fig. 5; Fig. S8). Most tropical areas predicted as suitable were identified as regions with high extrapolation risk (Fig. S8). For raw variables, the areas with extrapolation risk in North America included most of Canada and Alaska, whereas for PCs areas with extrapolation risk included large portions of Mexico and, the central-southwestern United States, as well as the islands north of Hudson Bay in Canada (Fig. 5).

Simulations of potential invasion

The simulations of potential sequences of colonization and dispersal of AGH in North America, starting from the known invaded localities, showed agreement among calibration schemes in predicting an invasion across the Pacific Northwest from southernmost Alaska to southernmost California in the United States (Fig. 6). In contrast, we found that the dispersal distance required to invade all the way to the East Coast of North America varied among calibration schemes. In the schemes using raw variables, the route of invasion to reach the East Coast goes from the Pacific Northwest down to California and Mexico, and then up the East Coast of North America. A dispersal distance of 10 cells (where each cell represents ~18 km) was enough to reach the East Coast (see Figs. 6A and 6C). For the scheme using the 50 km spatially-thinned occurrences and PCs, the invasion follows a more direct route from the Pacific Northwest to the East Coast that goes through the United States, and the required dispersal distance to reach the East Coast was only 4 cells (Fig. 6B). Finally, in the case of country-density thinned occurrences and PCs, the invasion goes from the Pacific Northwest through Canada to the Atlantic, and then down the East Coast to the United States. A distance of 8 cells was needed to make this invasion route possible (Fig. 6D).

Honey production and native bee richness

The areas in North America that our models identified as highly suitable for AGH overlapped broadly with the states where honey production is highest. This overlap was particularly noticeable in southern Mexico and in some states of the Pacific Northwest and eastern US (Fig. 7). We found a similar pattern with the species richness of Bombus and Melipona.