For this study, we focused on the secondary introduction event of the hornet in Portugal and used all of the available Portuguese presence data of V. velutina (8610 records of nests, from 2013 to 2018). Data was obtained from Bombeiros Voluntários de Viana do Castelo and from the online platform ‘STOPvespa’ (http://stopvespa.icnf.pt/), which is managed by the Instituto da Conservação da Natureza e das Florestas (ICNF) and aggregates all validated Portuguese records of V. velutina nests that were previously registered in the platform by citizens. To avoid spatial autocorrelation, we reduced the number of occurrence data points through the spatially rarefy occurrence data tool (pixel size resolution: 300 m) in SDMtoolbox (Brown 2014) in ArcGIS 10.4.1 (ESRI 2016); this resulted in a total of 7847 points (Suppl. material 2: Appendix A). To calculate the continuous distribution area of V. velutina (2013 to 2018) we drew a 5 km buffer around each nest. The workers of V. velutina probably forage less than 1000 m from their nest, and this buffer of 5 km corresponds to the estimated maximum homing ability of the species, since few workers have the ability to perform long trips (Poidatz et al. 2018). Moreover, and according to (Lioy et al. 2019), most new nests (>90%) in NW Italy were located within this radius to their nearest source of the previous year. All the contiguous buffers were aggregated to establish each year’s continuous distribution; all records outside the continuous area of the previous year were considered expansion nests. From these expansion nests, those located within the new continuous area were considered to result from diffusion dispersal and those found outside this new limit were considered outposts (i.e., an evidence of jump dispersal, through either self-mediated or human-mediated).

To identify the range expansion trajectory, we calculated the annual increment in the continuous area. The number of new outposts per year was counted and their contribution to the overall expansion was estimated by identifying those outposts that could have functioned as a source for other nests. To ascertain a possible origin for each expansion colony and outpost, we compared its distance to the nearest edge of the continuous area and to the nearest outpost of the previous year; all the records for which the difference between both distances was lower than 5 km (corresponding to 1419 records) were discarded, being considered of non-attributable origin.

The yearly expansion resulting from diffusion dispersal along the N-S and W-E axes was estimated by measuring the distance to the south and east between consecutive limits of the continuous distribution area. The number of new nests established exclusively to the south and east from the previous continuous limit was counted and we identified how many of these were outposts. Yearly, for each outpost, its distance to the nearest source of the previous year was measured. To test for an acceleration of both types of expansion, the slopes of the relationships between these distances and year was compared with zero.

Assuming that the same variables influencing distribution have the potential to promote its dispersal, we considered three climatic and eight land cover and anthropogenic variables (see Suppl. material 1: Table S1). Variables in this study were adapted from Bessa et al. (2016), with the following changes: i) we excluded NDVI and isothermality, ii) we used the distance to each specific land cover class instead of its percentage because using distances assures better performance for landscape features (Rainho and Palmeirim 2011), iii) we incorporated the classes “distance to forest” and “index of human influence” and iv) we included predictors related to the distance to linear structures (motorways and railways). To avoid collinearity, we inspected if there were highly correlated variables (r≥0.70; Dormann et al. 2013) (see Suppl. material 1: Table S2). For this analysis we used three datasets: i) bound records (the series of points defining a minimum convex polygon that represents the leading edge of the continuous invaded area for each particular year, presumably resulting from dispersal by diffusion), ii) all outposts (resulting from jump dispersal) and iii) > 18 km outposts, representing the subset of outposts located more than 18 km (distance travelled in flight mill experiments – pers. comm., Dr. D. Sauvard, INRA, France to the authors of Robinet et al. 2017) from the continuous area of the respective year, for which there is probably a higher contribution of human-mediated dispersal.

To assess which variables influence the dispersal of V. velutina for each of the three datasets we used generalized linear mixed models (GLMM) with the package ‘lme4’ (Bates et al. 2015) in R (Core Team 2019). We began by running full models with climatic, land cover and anthropogenic drivers simultaneously. As climatic variables are acknowledged to be the main factors influencing the species distribution across varying spatial scales (Pearson and Dawson 2003), we decided to run additional models with land cover and anthropogenic variables only, in an attempt to find other possible predictors at a regional scale. For each dataset, we set the dependent variable as the minimum distance of the records to the continuous area of the previous year (we discarded three records that were located less than 5 km from an outpost established in the previous year, as that could be an offshoot of that outpost). To detect collinearity between explanatory variables we used the Vifstep function in the usdm R package (Naimi et al. 2014) to calculate the variance inflation factor (VIF) and excluded the variables in models with a VIF value greater than the threshold (th=3). A variable “year” was included as a random effect to account for yearly climatic variations that may affect the dispersal of the hornet. We then selected the best model (using the Akaike Information Criterion – AIC) with the dredge R function, and generated average estimates of the effect of each variable using the model.avg R function (models with delta AIC values < 2) from the MuMIn package (Bartón 2009). The results were plotted using the package visreg (Breheny and Burchett 2017).

As we verified that one anthropogenic predictor (distance to motorways; see Results) was influential on hornet jump dispersal we decided to further explore the data. First, we inspected if the outposts’ distance to motorways was random, i.e., we tested whether motorways may be acting as drivers of jump dispersal. To accomplish this, we generated a twin random point for each outpost, located at the same Euclidean distance to the continuous distribution area as the outpost, and compared their distance to motorways with a paired samples Wilcoxon test. Second, for both data sets of outposts we ran another GLMM model, but this time with the distance to the entire road network to inspect the relative importance of each road category in hornet jump dispersal.

To generate a risk map of V. velutina dispersal and identify regions most at risk of imminent invasion, we combined information from suitable areas (regions with rainy winters and pleasant summers, mainly located along the Atlantic coast: Verdasca et al., unpublished data) with the geographical information of the significant dispersal predictors of a model with climate, land cover and anthropogenic variables (see Suppl. material 1: Table S4). These predictors were combined according to their estimates to produce a dispersal map. To define the risk areas around motorways, we analyzed the pattern of the number of outposts as a function of distance to motorways. As the number of new nests established alongside the motorways decreased linearly with distance up to 17 km from the highway (after that there was no apparent relation with distance – see Results), we calculated this relation to estimate the width of the areas that contained 50% and 75% of the outposts.

The number of outposts varied across the different years from a minimum of 4 in 2016 to a maximum of 46 in 2015. Such outposts had a very high importance for the expansion of the hornet. Indeed, the number of new expansion nests that were located near the outposts established in the previous year was higher than the number of new nests found near the previous continuous limit in all years except 2017 (Table 1).

In the first three years (2014–2016), the mean distance of new nests to the nearest outpost was lower than the distance to the continuous area (Suppl. material 1: Fig. S1). The reverse scenario occurred in 2017 and 2018, when almost all outposts were established in Mediterranean-climate regions. As the core distribution area expanded to the south and east, some outposts that gave rise to new nests nearby were engulfed into the continuous distribution area (i.e., coalescent colony model; Fig. 2).

Outposts established southwards were over 3 times more frequent than those established eastwards (Table 2). There was a decrease in the number of successfully established nests since 2016, especially southwards (Table 2). The slope of the relationship between time and dispersal distance to the south and east was not significantly different from zero, for both types of dispersal (diffusion or jump dispersal) (Suppl. material 1: Table S3). Most outposts (90%) were located more than 18 km from the continuous area of the previous year.

Models with both climatic and land cover variables explained more variability of the dispersal patterns of V velutina than models solely with climatic or land cover variables (Suppl. material 1: Tables S6–S8). A climatic variable – precipitation of the wettest month – was the single variable selected in all the models. Temperature annual range was the only additional climatic variable identified as a driver of diffusion dispersal (Table 3, Suppl. material 1: Fig. S2, see also Suppl. material 1: Table S4 for a model with climatic and land cover variables). In the models with land cover variables only distance to shrub land (plus natural meadows) was identified as influential on diffusion dispersal; however, for jump dispersal, distance to motorways was the only significant predictor (Table 3, Suppl. material 1: Fig. S2). Distance to the entire road network (instead of distance to motorways) had no effect upon either dispersal pattern (Suppl. material 1: Table S5). In both datasets (all outposts and the subset located more than 18 km from the continuous area of the previous year), outposts were significantly closer to motorways than expected by chance (all outposts: 103 pairs compared, V = 1820, p-value = 0.0024 (also see Fig. 4); >18 km outposts: 93 pairs compared, V = 1306, p-value = 0.0004). We found that the number of outposts decreased in a linear function (y = -0.6103x + 10.61) with distance to motorways in a 17 km-width strip along the motorways (Fig. 4 – dark grey bars). Further away, there is no apparent relationship with distance. In fact, 50% of the new outposts were located within a 6 km wide buffer zone alongside motorways and 75% up to 12 km from these linear structures.

Figure 3.

Risk of dispersal of V. velutina in Portugal evidencing the buffers alongside motorways where dispersal is likely to be mostly human-mediated. The unsuitable area for the species (Verdasca et al. unpublished data) is depicted in a pale yellow. Almost all the isolated suitable areas located in the south of the country are also at risk as they are connected by motorways to other suitable regions.

The invasion of V. velutina is occurring at a slower pace in the northwest of the Iberian Peninsula (spread rate of approximately 45 km/year to the south and 20 km/year to the east) than in other temperate macroclimate regions (e.g., France), but faster than in other Supramediterranean climates (e.g., Italy). In the first few years of the invasion the number of new established nests was much higher near outposts than near the continuous distribution area, an indication that jump dispersal played an important role in the acceleration of the invasion process. Besides climate (namely, precipitation of the wettest month and the annual range of temperature), we found the distance to shrub lands to be influential in the dispersal of V. velutina. This finding adds new information to a previous study which also showed that land-use (namely, percentage of agricultural fields) has an important role in the expansion of this species at regional scales (Bessa et al. 2016). We also revealed that one anthropogenic driver (motorways) was important for the jump dispersal events of this flying insect, highlighting the role of these linear infrastructures in accelerating the natural invasion dynamics of V. velutina and the need to reinforce early monitoring programs in a 6 km wide buffer around motorways.

Figure 4.

Distribution V. velutina outposts (dark gray) and random points (light grey) according to the different classes of distance to motorways (km). The estimated linear function found up to 17 km is “No. nests = - 0.6103 * (Distance to motorways) + 10.61”.

From the initial propagule found in the north of Portugal (near the coast), self-mediated dispersal has been occurring faster towards the south than towards the east. The western Iberian Peninsula encompasses different bioclimatic belts (Mesotemperate, Supramediterranean, Mesomediterranean and Thermomediterranean, Rivas-Martínez et al. 2017) that spread more along the North-South axis near the Atlantic coast than along the West-East axis. In Portugal, V. velutina is therefore faced with two transition zones, differing in extent, from temperate to Mediterranean climates, as its expansion is predominantly occurring along the Atlantic coast. To better predict the risk of invasion in the short run, it is important to disentangle the spread rate across different climatic gradients. For instance, the average rate of V. velutina expansion here identified (45 km/year to south and 20 km/year to east) is different from the one recently estimated for Portugal (37.4 ± 13.2 km/year, but considering all the directions, and therefore an intermediate value) (Carvalho et al. 2020). We believe that it is important to refine the estimates, as the spread rate of this invasive hornet is clearly not uniform across Portugal, and this same invasion process may occur in other temperate-Mediterranean transition zones. We acknowledge that since our records are reported by citizens, a bias in V. velutina detection may be occurring due to more identifications in areas with higher population density (i.e., along the coast). However, this is a very mediatic species in Portugal, and most people are aware of this and its impacts on beekeeping, agriculture, and public health. Despite the lower population density in the eastern part of the country, there are still some important cities, numerous small villages and, even more important, more beekeepers in these rural regions. Given that we used the outermost records for each year (and not the density of records) for most of our estimates, we think that the major patterns detected, such as the differential expansion along the North-South and West-East axes, are barely affected by differences in human density. We did not find a difference in the distance of establishment of new nests between the two directions; however, there was a substantial difference in the number of new nests, as those established southwards were twice the number of those established eastwards. The climatic transitions are more abrupt towards the east, where the hornet is now facing Mediterranean climatic conditions (i.e., drier, higher range of temperatures), which may explain why the species has more difficulty in establishing new nests in this direction. This may be the reason for the decrease in the numbers of established nests since 2016 and supports the importance of climate for the expansion rate of this hornet. In fact, the expansion area is not increasing exponentially as would be expected if diffusion was occurring equally across all directions. In France, the species spread rapidly toward the northeast and not so much to the south (Robinet et al. 2019). In Portugal, the rate of expansion was lower in 2018, potentially due to the major and uncontrolled wildfires that occurred in 2017 precisely over the distribution limit of V. velutina in that year. Indeed, by December 2019, the spread rate towards the south was again near 50 km/year (see http://stopvespa.icnf.pt/ by ICNF; also check the invaded area in Portugal by May 2021 in Fig. 2). Extrapolating to other temperate-Mediterranean transition zones, such as the Italian and Balkan peninsulas, the rate of expansion and invasion pattern may be similar.

As in other countries, V. velutina in Portugal is dispersing by both diffusion and jump dispersal. This same pattern was noticed in France (Robinet et al. 2017) and Italy (Bertolino et al. 2016; Lioy et al. 2019), as well as for other social Hymenoptera (like the Argentine ant; Suarez et al. 2001). The frequency and distance of jump-dispersal events are thought to be stochastic, and therefore difficult to predict. For species that spread through stratified diffusion, the distance and rate at which new foci are created through jump dispersal may be more important than the rate of spread through diffusion from established foci (Suarez et al. 2001). Here, the successful long-distance dispersal events played an important role in the expansion of V. velutina in almost every year after its establishment, as the number of new nests was higher near outposts than the boundary of the continuous distribution area (Suppl. material 1: Fig. S1). These results support a coalescent colony growth model, similar to prior tudies that found outposts to accelerate range expansion (Shigesada et al. 1995). Previous studies from Italy (Bertolino et al. 2016; Lioy et al. 2019) found the dispersal of V. velutina to be hindered by high mountain ranges (above 700 m), and therefore argued that this may be one of the main reasons for the low spread rate in Italy (18 km/year: Bertolino et al. 2016) compared to France (78 km/year: Robinet et al. 2017). Spread rates similar to those in Italy were registered in Korea (10–20 km/ year), although the low spread rate there may be due to competition with six other hornet species (Choi et al. 2012). Nonetheless, the constant spread rate observed so far in Portugal may begin to decrease southwards when the species reaches Mediterranean climates. Although all outposts were located within the suitable area for the species, V. velutina is not (yet) in geographical equilibrium, since, according to our former work, there are still suitable areas to the south and east (Verdasca et al. unpublished data). This apparent limit on the establishment of outposts corroborates our estimates about the adequate areas for the species; however, as the species reaches its estimated limits, it is important to assess how robust they are, as the colonization of adjacent areas, or even the adaptation to novel environmental conditions is possible. If these limits hold, this means that jump dispersal will be the only dispersal mechanism allowing the species to reach the isolated suitable areas in the south of the country.

The dispersion of V. velutina is affected by precipitation and temperature gradients, a result that is similar to those of other studies that modeled the hornet’s bioclimatic niche (Villemant et al. 2011; Verdasca et al. unpublished data). Besides precipitation of the wettest month and the annual range of temperature, distance to vegetated, but treeless landscapes (covered by shrubs and natural meadows) seem to favor diffusion dispersal. As shrub lands provide a wide variety of nesting sites and food resources for wild pollinators (Chaplin-Kramer et al. 2011), hornets will probably need to fly over longer distances until such pollinator suitable habitats can be reached. In regions that are climatically suitable, the presence of shrubs may thus reduce hornet dispersal. However, shrub land cover is probably not related to the large difference between the eastward and southward rate of expansion, as this habitat is regularly found across the suitable area for the species in Portugal.

Precipitation in the wettest month, and motorways, were the only factors identified as drivers of jump dispersal, but the role of motorways in the dispersal of the hornet was only detected when the climatic predictors were not included in the models. This is in line with the scale dependencies outlined by Pearson and Dawson (2003) – different processes are more important at different scales i.e., at a continental scale, climate can be considered the dominant factor, whilst at more local scales factors including topography and land-cover type become increasingly important. Further down the hierarchy, if conditions at higher levels are satisfied, factors including biotic interactions and microclimate may become significant (for details on hierarchical modeling framework see Pearson and Dawson (2003)).

The fact that motorways were important predictors of outposts is an indication that they may have resulted from human-mediated dispersal. Yet, as motorways are heavily used by people, a potential bias in the detection of nests near these human infrastructures may have occurred. However, most motorways pass through remote places with low population density, and people cannot stop their cars over vast extensions. Therefore, it is unlikely that nest reports come from people using the motorways. Jump dispersal events were predicted by motorways, but not by all roads, railways, or the index of human influence, a variable highly correlated with human population density (e.g., cities). This is an indication that the establishment of outposts is probably mediated through the movement of vehicles and goods, such as wood products and bark or man-made goods (e.g., ceramic pottery associated with garden trade), which in Portugal occurs mostly through the motorways. These products provide suitable refuges for hibernating inseminated V. velutina queens (Marris et al. 2011); indeed this was the most probable route of incursion of V. velutina in Europe, on pots imported from coastal China, near Shanghai. However, it is also plausible that a dispersing gyne may simply land on a car that then travels a good distance away and starts a nest there. As V. velutina gynes can fly over long distances and generate stochastic patterns of spread similar to those resulting from human-mediated dispersal (Robinet et al. 2017), it is probable that some of the records may have originated from self-mediated dispersal. However, the overall detected effect of motorways regardless of the distance group considered, together with the decreasing trend in nest abundance as one proceeds away from the motorways, are difficult to explain through self-dispersal alone.

The association of nest establishment with motorways was only found for outposts (50% and 75% of them established within a 6 km and 12 km wide buffer zone alongside motorways, respectively), and not for records that originated from diffusion dispersal. Our findings corroborate a previous study in Italy (Porporato et al. 2014) where a high number of the observations and captures of V. velutina in bait-traps were recorded near highways, emphasizing that freight traffic can contribute to the transport of this species far from the invasion front. The fact that Bessa and collaborators (2016) did not find any relation with the road network could be due to the very restricted region that was used in their study (roughly 10% of the area that we used here). In France, Robinet et al. (2017) used human population density as a proxy for trade to test jump dispersal, not taking into account the road network, and concluded that the rapid spread of the hornet may not be necessarily mediated by humans. So, it is possible that long-distance dispersal events that occurred in France may also have contributed to unintentional introductions via motorways.

Despite it being extremely difficult to provide evidence for early introductions, other social insects have also probably been transported accidentally by humans over long distances since the establishment of long-distance trade routes (Bertelsmeier 2021). For example, New Zealand had no social wasp species prior to human colonization, but over the last century has been invaded by several species of social wasps (Lester and Beggs 2019). Indeed, Vespula germanica and Vespula vulgaris, both native from Eurasia, have become widespread throughout the New Zealand causing major impacts to native biodiversity (Lester et al. 2014). In Argentina, where V. germanica is also invasive, the observed stratified geographical expansion pattern (which frequently exceeds 30 km per year, although faster to south) does not match the observed queen dispersal abilities (only a few hundred meters naturally to find nest sites), suggesting that human-aided transport of hibernating queens is the central driver of the current distribution of these wasps in the country (Masciocchi and Corley 2012). At more local scales, the anthropogenic influence on the spread of invasive insects was also demonstrated. For example, the distance to railroad tracks influenced the spread of the invasive termite species Reticulitermes flavipes (Perdereau et al. 2019).

Identifying pathways that facilitate the dispersal of invasive species is essential for informing efforts to contain invasions (Suarez et. al. 2001). To be successful, every invasive species control program must consider the probability of detecting the species and the cost of the process. In this work, we showed that 50% of the presumed new nests resulting from human-mediated long-distance dispersal established within a 6 km wide buffer along motorways. To raise this proportion to 75%, the buffer must be increased to 12 km, which represents an almost 70% increase in the area to be surveyed. Based on results here, effective measures to contain V. velutina invasions should include early monitoring programs in a buffer of 5 km (the maximum homing ability that few hornet workers can reach: Poidatz et al. 2018) around the continuous distribution area of the previous year, and 6 km (ideally 12 km) around motorways. If the climatic conditions are met, the vicinity of the main roads is susceptible to be colonized faster through human-mediated transport. Even in highly fragmented habitats, the main roads can connect isolated suitable areas. For instance, in the regions at risk in southern Portugal, the area to be surveyed can be limited only to climatically favorable regions that are reachable by highway. This is particularly relevant in southern Portugal where the isolated fragments of suitable landscape are economically very important for beekeeping activities. The early detection and control of nascent populations in these areas may be a good way to manage its spread, rather than focusing efforts on established invasion fronts. Local outreach activities, especially those targeted to transportation companies, should also be prioritized to prevent the European motorway network from becoming an invasion route for the hornet to new countries. However, different types of cargo do not carry the same risk of being infested (as different species may differ in their commodity associations). Therefore, focused biosecurity policies for V. velutina, are needed, particularly targeted to the interception of wooden products’ transportation and man-made goods associated with garden trade, due to the potential of these commodities to shelter hibernating queens. It is also important to promote control actions on ports of species entry, namely harbors along the coast.