Invasive hornets on the road: motorway-driven dispersal must be considered in management plans of Vespa velutina

Understanding the mechanisms that potentiate the dispersion of an invasive species is essential to anticipate its arrival into new regions and to develop adequate management actions to minimize damage to biodiversity and society. One of the most successful invaders in Europe, the yellow-legged hornet (Vespa velutina), is dispersing through self-diffusion and jump dispersal. Using information on species occurrence in Portugal from 2013 to 2018, this study aimed to understand the range expansion trajectory of V. velutina and to identify the role of climate, landscape and anthropogenic variables on the two mechanisms of spread. We found that in Portugal the invasion is proceeding faster southwards (45 km/year) along the Atlantic coast than eastwards (20 km/ year) where the climatic suitability gradient is more compressed, with jump dispersal playing an important role in this difference and in the acceleration of the invasion process. Dispersal by diffusion was best explained by the annual range of temperature and precipitation of the wettest month, with distance to shrub land also having an important role. Additionally, jump dispersal appeared to be facilitated by motorways, hinting at the role of human-mediated dispersal. Indeed, the number of nests that resulted from this dispersive mechanism were significantly closer to motorways than expected by chance. To prevent the dispersal of V. velutina into Mediterranean regions, and in addition to a special attention to the advancing front, early monitoring programs should also target a buffer zone on both sides of motorways, and at freight shipping hubs.


Introduction
Invasive species can have important environmental and socioeconomic impacts. Knowing the dispersal routes of such species is crucial to anticipate their arrival and define adequate management practices in a timely fashion. Invasiveness (a dynamic property of the species) and invasibility (a property of a location that can change with anthropogenic disturbance, seasons, climate change), two key components of biological invasions, are thought to be primarily determined by species' dispersal ability and habitat suitability, respectively (Brooks 2007). Invasiveness is thus mediated by life strategies (Pysek and Richardson 2007), whereas invasibility is related to local conditions at the site, habitat or landscape levels (Vicente et al. 2010). For a risk assessment of the vulnerability of a site to an invasive species it is therefore important to have detailed information on both the species' ecological niche (i.e., static information, such as climatic tolerance) and movement ability, which will regulate if and when suitable areas away from the site where the invasion started will be reached (i.e., dynamic information). Including such information on distribution models of invasive species can help to distinguish suitable habitat that is, or can be potentially occupied, from suitable habitat that is inaccessible (Miller and Holloway 2015).
As different environmental conditions and landscape heterogeneity may accelerate or hamper the invasive process (Hastings et al. 2005), it is also important to identify the patterns of range expansion. Three different trajectories of range expansion versus time can be considered: linear (Andow et al. 1990), biphasic (with an initial shallow slope followed by a steep linear slope), and accelerating with time, quickly reaching the saturation phase (Shigesada et al. 1995). The three expansion patterns occur through either self-mediated dispersal from an initial location (diffusion; Lockwood et al. 2007) or jump dispersal into regions relatively far from the core distribution area without colonizing the regions in between, leading to the establishment of nascent colonies (Suarez et al. 2001), "outposts" hereafter. If outposts establish in environmentally suitable areas the species can continue the expansion process from there, accelerating it. When a species spreads through both processes (natural diffusion and jump dispersal), stratified diffusion occurs (Hengeveld 1989;Suarez et al. 2001). In such a process, the initial range expansion occurs mainly by diffusion, but as the area of the founder population expands, new nests created by long-distance events accelerate range expansion in later phases (Shigesada et al. 1995; see also some insect examples in Andow et al. (1990)).
Insects are the dominant group among non-native terrestrial invertebrates in Europe (Roques et al. 2009). Social insects, in particular social Hymenoptera, are likely to become successful invaders due to their excellent dispersal abilities, high reproductive rates in an annual life cycle, broad diets and habitat ranges, colony initiation by a single inseminated queen, and to their close association with human transportation with relatively low probability of detection (Moller 1996;Beggs et al. 2011). These characteristics favor the invasiveness of eusocial insects, giving them a plasticity of responses that allow their survival and establishment in new environments. Impacts of invasive social insects include changes of ecosystem functions, competitive displacement of native species, hybridization with native species, threats to human or animal health through stings or the transmission of pathogens (Schneider et al. 2004;Lester and Beggs 2019). One example of an invasive eusocial insect that is currently spreading in Europe is the yellow-legged hornet, (Vespa velutina Lepeletier, 1836), an aggressive predator of honeybees and wild pollinators, that is putting honey and agricultural production at risk (Monceau et al. 2014;Verdasca et al. 2021). The high dispersal ability of this invasive combined with the lack of true competitors and with the availability of food resources in Europe, has been favoring its expansion in this continent. After the accidental introduction of V. velutina in France in 2004, probably by a single female originating from temperate south-eastern China (Arca et al. 2015), the invasion spread to other European countries -Spain, Portugal, Belgium, Italy, Germany, United Kingdom, Netherlands and Luxembourg, being the hornet now considered an invasive alien species of concern in the European Union (European Commission 2016). In 2011, a new invasion was detected in the north of Portugal (Grosso-Silva and Maia 2012) that then expanded southward into the center of Portugal (Carvalho et al. 2020) across a climate gradient between Temperate and Mediterranean bioclimates (Rivas-Martínez et al. 2017;Sayre et al. 2020) and northward into Galicia (Spain) (Rodríguez-Flores et al. 2019). The current European distribution of this species, mostly concentrated along the Atlantic coast, aligns with the climatic preferences predicted by Villemant et al. (2011). As a central place forager, workers of V. velutina optimally explore foraging areas 500 m to 800 m from their nests, although the maximum homing ability (the maximal distance an individual is able to travel on its way home) of the hornet is much greater (5000 m; Poidatz et al. 2018).
The spread of V. velutina in Europe has been considered a stratified diffusion process, including a mixture of natural diffusion and jump dispersal events (Bertolino et al. 2016;Robinet et al. 2017;Lioy et al. 2019). Jump dispersal may occur by two different processes: self-dispersal of gynes (queens of the next generation), which are able to fly long distances on their own (18 km/day in flight mill experiments -pers. comm., Dr. D. Sauvard, INRA, France to the authors of Robinet et al. 2017), or through human-mediated dispersal (Robinet et al. 2019). Due to the difficulty in disentangling these processes, only occurrence points found more than 78 km (yearly spread rate in France) from their nearest neighbors were considered to result from long distance dispersal events by Robinet et al. (2017). However, Bertolino et al. (2016), working in Northwest Italy, considered a much shorter distance as the limit over which human-mediated dispersal was considered the most likely explanation for jump dispersal events (yearly spread rate in Italy = 18.3 ± 3.3 km). Besides the different methodological approaches (including the scale of analysis), the differences between the two studies are attributed by Bertolino et al. (2016) to the Italian mountainous territory, when compared with the largely flat north western and central France, while Robinet et. al. (2017) justify the high spread rate in France by the ability of founder females of flying long distances in flight mill experiments.
Precipitation and temperature are thought to be the strongest predictors of the invasive range of V. velutina (Villemant et al. 2011), with land-use also playing an important role at a regional scale (Bessa et al. 2016). This hornet is now colonizing the Mediterranean peninsulas, and therefore being exposed for the first time to a Mediterranean climate. As the species spreads along the transitional temperate-Mediterranean climate regions, it is important to assess the extent of suitable area at a finer resolution, as well as to understand how the invasion process is unfolding and the role of humanmediated dispersal. A recent study highlighted the need to identify the colonization pathways and plan management approaches to halt the spread of V. velutina in Portugal (Carvalho et al. 2020). At a regional scale, the identification of a limited number of key variables explaining the dispersal and establishment success of V. velutina can facilitate the creation of effective preventive and control measures. In this study, we incorporated land cover and anthropogenic drivers to predict the risk of invasion by V. velutina within the transitional temperate-Mediterranean zones and infer colonization pathways. More specifically, we (i) assessed the roles of diffusion and jump dispersal on V. velutina expansion into Mediterranean-type climates and (ii) identified which environmental attributes are most influential on the direction and speed of its dispersal.

Nest occurrence data in Portugal
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).

Range expansion
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.

Environmental drivers
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.

Data analyses
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.

Availability of data and material
Due to privacy reasons, public data is available in Suppl. material 2: Appendix A of Supporting Information in a resolution of 5 × 5 km. The precise geo-localizations that support the findings of this work (which were used under license for the current study) are available from the authors upon reasonable request and after permission of the entities holding the data. Data requests can be addressed to the corresponding author, who will make them available jointly with the different entities that hold the data. Any further information can be obtained by request to the authors.

Results
From 2013 to 2018, the area occupied by V. velutina in Portugal experienced a 25-fold increase (from 845 km 2 to 20,561.26 km 2 ) in a linear manner without acceleration or deceleration (Fig. 1). Expansion was much faster along the north-south axis (45 km/year) than along the west-east axis (roughly 20 km/year), regardless of taking place in temperate (in 2014 and 2015) or Mediterranean climate regions (since 2016) ( Fig. 1 and Fig. 2).

Range expansion
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.

Environmental drivers
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  . 4); >18 km outposts: 93 pairs compared, V = 1306, pvalue = 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.

Major findings
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.

Range expansion
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.  (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.

Roles of diffusion and jump dispersal
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 studies 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.

Influential environmental attributes on the direction and speed of V. velutina dispersal
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)).

Motorways facilitate jump dispersal
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).

Implications for the management of Vespa velutina invasion process
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.  Table S1. Climate, land cover and anthropogenic variables with potential to affect the behaviour and establishment of Vespa velutina. Tables S2.
Correlation matrix of the climatic, land cover and anthropogenic drivers that have the potential to affect the behaviour and establishment of Vespa velutina in Europe: a) bound records; b) outposts; c) outposts >18 km. Table S3. Relation between dispersion distance (to south and east) vs time either for self mediated or jump dispersal by testing through t-test the significance of the slope of the regression when compared to zero (H0: the slope of the regression line is 0). Table S4. Effects of climatic, land cover and anthropogenic variables on the dispersal of V. velutina. Table S5. Effects of the different land cover and anthropogenic predictors on the dispersion of V. velutina, considering the predictor distance to the entire road network, instead of distance to motorways in bound records and both sets of outposts. Table S6. Set of best models with climatic and land variables according to the different datasets (1. bound records; 2. all outpost and 3. outpost 18 km). Table 7.
Set of best models with climatic variables according to the different datasets (1. bound records; 2. all outpost and 3. outpost 18 km). Table S8. Set of best models with land variables according to the different datasets (1. bound records; 2. all outpost and 3. outpost 18 km). Figure S1. Variation of the mean distance of the new records within a given year to the nearest potential source: continuous area or outpost. Error bars depict standard errors. Figure S2. Relation between the dispersion of Vespa velutina and significant climatic and land cover variables according to the different datasets: bound records, all outposts and outpost 18 km (see Table  3 of the main manuscript). Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.69.71352.suppl1