Describing the spatial spread in invasive species and understanding the underlying processes across different organisms have been long-standing research goals across many disciplines ranging from biology to mathematics and epidemiology (Fisher, 1937; Skellam, 1951; Andow et al., 1990; Okubo and Levin, 2001; Hastings et al., 2005; Alford, 2021). Among the stages of an invasion process, spread corresponds to invading organisms moving across the landscape and colonizing new habitats. Spread is a critical phase in the invasion process from a biological point of view (e.g., Andow et al., 1990) and control and management perspective (Epanchin-Niell and Hastings, 2010; Epanchin-Niell and Wilen, 2012; Robertson et al., 2020). Based on the early works of Fisher (1937) and Skellam (1951), most work on spread modeling has followed a mechanistic bottom-up approach (Hordijk and Broennimann, 2012), using integrodifference or reaction–diffusion models combined with a population growth model (van den Bosch et al., 1990; Kot et al., 1996; Soubeyrand and Roques, 2014). This approach has the advantage of providing an explicit representation of the different processes involved in the invasion process (Shigesada and Kawasaki, 1997; Hastings et al., 2005; Soubeyrand and Roques, 2014). It allows the estimation of spread descriptors by fitting models with mathematically explicit functional forms and parameters (e.g., Abboud et al., 2019; Bruzzone and Utgés, 2022). However, spread is a complex process, and the assumptions of these models could be too simplistic, making it difficult to obtain robust reconstructions of the invasion dynamics and the rate of spread through space and time (Hastings et al., 2005; Hordijk and Broennimann, 2012). In addition, most mechanistic bottom-up invasion spread models are usually defined using a deterministic framework (Kot et al., 1996; Lewis, 2000; Lewis and Pacala, 2000; Soubeyrand and Roques, 2014), ignoring the importance of extrinsic or intrinsic stochastic factors (Hastings et al., 2005). Extrinsic stochasticity sources may include the host’s actual spatial distribution, the climate effect on insect phenology, the presence of natural barriers (Carrasco et al., 2010). Intrinsic sources of stochasticity may include the exact form of the dispersal kernel or stochastic population dynamics (Lewis, 2000). On the other hand, empirical approaches to the description of invasion processes and invasion spread emphasize the role of spatial heterogeneity, temporal variability, and ecological interactions (Hastings et al., 2005). These more versatile top-down approaches depend on few or fewer assumptions, facilitatating their implementation (Hordijk and Broennimann, 2012). Thus, top-down approaches may provide spread descriptors without assuming specific mechanisms or explicit mathematical models (Hordijk and Broennimann, 2012), which may be argued to be their major disadvantage.

In the top-down approach, models based on network theory have successfully been applied to describe the spread phenomenon. Disease propagation is a special case of a spreading phenomenon studied using networks (e.g., Newman, 2002; Sattenspiel and Lloyd, 2009; Bjørnstad, 2018). In this case, nodes in the network may be either people or sites (e.g., cities), and links represent contacts or connections (routes). In these spatial networks of sites or cities, the routes are defined a priori (e.g., highways), and spread occurs through these links. However, a slightly different approach exists to analyzing invasive species spread dynamics. This approach describes the geometrical properties of species invasion patterns by using graphs such as the minimal spanning tree (MST) (Labra et al., 2005). But also describes the invasion process by successive invaded sites or sampling areas where an invasive species has been recorded. These invaded sites are represented in the two-dimensional geographical space as a set of nodes, with the Euclidean distance describing the geographical separation among them (Figure 1A). These points in space may be joined in many different ways by lines (or edges), forming a graph. The resulting graph may be defined as connected if there is an edge between any pair of nodes, as shown in Figure 1B, which includes several circuits or loops (Labra et al., 2005). It has been shown that for complex spatial patterns, the dominant pattern of connectedness may be described by a minimal spanning tree (MST), which is a graph that connects all the nodes with no circuits or loops. It only considers those connections that minimize the total length across the graph (Figure 1C; Prim, 1957). MST has a long history of application in biology and ecology to describe various spatial patterns (Dussert et al., 1986, 1987; Cantwell and Forman, 1993; Lockwood et al., 1993; Jones et al., 1996; Keitt et al., 1997; Wallet and Dussert, 1997; Bunn et al., 2000, Urban and Keitt, 2001; Labra et al., 2005). However, to fully describe the invasion process, the MST should be rooted at the first recorded invasion site and directed to reflect the temporal sequence of site invasion. Optimum branching (OB), minimum cost arborescence (MCA, Figure 1D), or directed minimum spanning tree (DMST) are all names that have been used to describe same problem. Given a directed rooted weighted graph, what is the minimum cost (sum of weights) graph with a unique path from the root to any other node? (Chu and Liu, 1965; Edmonds, 1967; Tarjan, 1977). In the context of the spread of an invasive species, if the original point of introduction and successive dated points where the species was detected are known, and if we assume that spread follows the shortest path, then we can use the minimum cost arborescence to reconstruct the spread of the species in the landscape. The resulting MCA is a directed weighted graph with weights equal to the Euclidean distance between nodes (sites), the in-degree for all nodes is always one, as it is a branching process where new sites are assumed to originate from a single previously invaded site. However, the out-degree can show some variability depending on the topology of the resulting network, as any given site may be linked to one or more newly invaded sites in the following time step. Hordijk and Broennimann (2012) provide an excellent description of this approach; however, despite its many advantages, this approach has been scarcely explored in studies on the spread of invasive species. Among the advantages of this approach are that it facilitates the estimation of dispersal rate through the calculation of spread velocity and acceleration and allows the detection of long-distance jumps. In addition, metrics show changes in response to spatiotemporal heterogeneity, and hence may be valuable to identify different environmental drivers that may be forcing or influencing the invasive spread. Here we show how this approach may be used to describe an ongoing biological invasion.

Exotic insect pests colonizing new ecosystems have caused multiple detrimental effects around the world, having a direct impact on agricultural and forestry production (e.g., Bradshaw et al., 2016), ecosystems services (e.g., Clark et al., 2010), human health (Lounibos, 2002), and cultural values (Manachini, 2015). Among the most important insect pests currently attacking fruit production worldwide is Drosophila suzukii. This species originated in South East Asia, but in 2008, it was detected simultaneously in the USA, Italy, and Spain (Rota-Stabelli et al., 2013). Since then, it has been detected in several regions like Hawaii, North America, Europe, Central America, South America, and Africa (Rota-Stabelli et al., 2013; CABI, 2022). Adults of this species are 2–3 mm long with red eyes, pale brown or yellowish-brown thorax, and black transverse stripes on the abdomen (CABI, 2022). Drosophila suzukii shows a winter morph with increased cold tolerance and larger wings (Shearer et al., 2016). In 2017, D. suzukii was detected for the first time in Chile (Rojas et al., 2019; Devotto, 2020), and in the following years, the species spread through the country. In Latin America, D. suzukii has at least 64 hosts (Garcia et al., 2022). Berries are among the most important fruit host of D. suzukii. In 2020 Chile had ~22,000 ha planted with these species (Garcia et al., 2022). In blueberries (Buzzetti, 2020), reported damage ranges between 1 and 1.5 tonnes/ha, equivalent to 4,000 US$/ha. The same author reported average direct damage per year between 1.2 and 2.7 tonnes/ha with losses of 5,000 to 17,550 US$/ha in cherries. After the first detection, the Agricultural and Livestock Service (SAG in Spanish) implemented an intensive monitoring program that continues until today, where daily georeferenced captures are registered and reported (SAG, 2021). Currently, this is an ongoing invasion; hence the estimation of the spread rate is a fundamental metric with significant economic and management implications. In this study, we apply MCA methods to reconstruct and understand the invasion dynamics of D. suzukii in Chile since the first detection in 2017. We used this graph-based approach to address the following research questions: (1) what is the magnitude of spread? (2) How does the rate of spread vary when different periods and geographic regions are compared? (3) What is the relationship between the spatial and temporal variation and the spread rate with the species biology? Our results will contribute to understanding this species’s invasion biology and could provide valuable tools for developing appropriate control measures.

Between 2017 and 2021, D. suzukii spread rapidly throughout the country. In four years, the species has colonized a ~1,000 km long strip in the central valley of Chile, ranging from 32° to 42° Lat S (Figure 2). The estimated MCA (Figure 3, Supplementary material 1) shows an initial phase with long paths and connections without a clear direction pattern (Figure 3). After this phase, a clearer north–east propagation pattern emerges (Figure 3). The median dispersal rate for the entire period was 8.8 meters/day with a 95% CI of (7.4–10.6) (Figure 4A). The coefficient of diffusion showed a median value of 19.6 meters²/day with a 95% CI (13.6–27.9) (Figure 4B).

A temporal look at the propagation pattern showed two clear peaks of spread during the autumn and winter of 2017 and 2018, where 75% of the total area colonized until 2021 in Chile was reached. During the autumn–winter of the first two years, dispersal rates reached a maximum of ~110 meters/day, and the diffusion coefficient reached a maximum of 3,000 meters²/day (Figures 5A,B). Both values are several times higher than the median values for the entire period. In a similar way, spread shows several acceleration–deceleration phases, also mainly in the autumn–winter season (Figure 5C). After the first two years, the dispersal rate stabilizes around the median values. The dispersal rates follow a log-normal distribution, typical of processes that involve multiple scales, as can be seen when we compare the empirical probabilities with the theoretical probabilities from the best-fitting log-normal (Figure 5D, see Discussion).

Drosophila suzukii has spread rapidly along Chile, occupying most of the region with berries plantations in less than four years. This high spread capacity is a well-known trait of this species. Hauser (2011) and Calabria et al. (2012) have described the fast dispersal of D. suzukii in the USA and Europe, respectively. In particular, Calabria et al. (2012) pointed out a dispersal rate of 1,400 km/year, probably combining active and passive dispersal. In Chile, the species did not reach such an extremely high dispersal rate, but 1,000 km in 3 years is still a significantly high spread capacity. As it was suggested by Calabria et al. (2012), this high dispersal rate is the consequence of a combination of active (diffusion) and passive (human transport) spread, which is called stratified dispersal (Shigesada and Kawasaki, 1997). In our case, observed dispersal rates follow almost perfectly a log-normal distribution, a well-known heavy-tailed distribution. This result suggests that the occurrence of extreme, long-distance dispersal events is “more common” than that expected under diffusion-only or Gaussian spread dynamics, in line with the observation of Calabria et al. (2012).

The general estimation of 8.8 meters/day is close to the results obtained experimentally by Vacas et al. (2019). These authors, performing mark-release-recapture experiments, show that most individuals of the species were recaptured at a distance below 10 meters after 24 h. However, several studies have reported that D. suzukii is capable of long flights. Tait et al. (2018), also performing mark-recapture experiments, reported high variability in the dispersal rate of this species. Depending on elevation, these authors reported several individuals with almost no movements after several weeks of sampling, whereas other individuals show dispersal rates as high as 9,000 meters/month (~300 meters/day). In the same vein, Wong et al. (2018), described individuals moving over one km in a single flying event, which also supports the occurrence of long-distance dispersal events.

The expansion D. suzukii continues today, most likely through human-mediated transport, to the arid-semiarid region of Chile. This arid region (including the Atacama desert) in the north and the Andes mountains in the east act as soft barriers for the dispersal of this pest, although D. suzukii arrived in Chile, probably from Argentina (Garcia et al., 2022). However, this region contains much less vegetation and crops, and dispersal will probably be slower.

However, this median estimation can be misleading, as there is high variability in the dispersal rate of the insect depending on the seasons. Our results show two main peaks of dispersal in the first years after introduction. The dispersal rate during these peaks reaches 110 meters/day, more than ten times higher than the median value. These two peaks occurred during the first autumn–winter seasons. In southern Chile, D. suzukii individuals are mainly caught during the autumn–winter season (Drosoalert, 2020). The number of males and females caught during May–June is four to five times higher than the number of individuals caught during February–March (Drosoalert, 2020). In this season, fruit availability is lower, which can promote longer flight events (Little et al., 2020). Also, in Chile adults show a characteristic winter morph, with increased cold tolerance and larger wings which may be interpreted as an intrinsic factor associated with higher dispersal rates (Shearer et al., 2016). However, Tran et al. (2022) indicate that the winter morph does not show a higher-flying performance than the summer morph, despite its apparent morphological advantages. In this regard, our data is clear in showing a faster and accelerated spread in the autumn–winter seasons, but the determination and testing of the exact ecological mechanism behind this rapid spread will require new studies. After the initial fast propagation during the first two years, dispersal rates remained close to the median values, without new peaks, probably due to the species reaching most of the suitable habitats. In our MCA, this can be inferred by the lower median dispersal rates and few new nodes (colonized sites) after the second autumn–winter season (Figure 3). When we compare our results with those obtained in other species, we observe that the median dispersal rate of D. suzukii in Chile is lower than most of the results compiled by Evans (2016). However, dispersal rates during dispersal peaks can be comparable to those observed in some forest insects like Tomicus piniperda or Sirex noctilio (Evans, 2016).

On the other hand, extrinsic factors play an important role in modulating dispersal rates. For example, meteorological factors like temperature (Lantschner et al., 2014; Leach et al., 2019) or wind speed (Leitch et al., 2021) can be key for explaining the observed spread dynamics during an invasion. The high dispersal rates described by Tait et al. (2018) were observed from high to low elevations, which would suggest a major influence of winds in these events. This behavior has also been detected in other Drosophila species, where even longer wind-assisted jumps of ~12 km are possible (Leitch et al., 2021). In a more general analysis, Tait et al. (2020) showed that individual dispersion is a function of meteorological factors like temperature and humidity but also depends on the local diversity of alternative hosts. In this context, our results reinforce the complexities of the spread process and highlight the multiple difficulties of extracting detailed, ecologically meaningful information from real data. Intrinsic factors like population dynamics, phenology, or the emergence of seasonal morphs interact or depend on extrinsic forces like temperature, wind, or host availability, creating particular spread dynamics that can hardly be captured by one model, even in laboratory experiments (Melbourne and Hastings, 2009).

Our use of MCA allowed us to estimate several metrics of spread, successfully identifying both spatial and temporal variations. In this regard, the use of MCA, a rooted, directed weighted tree, shows several improvements over previously used MST, a minimal weighted tree method. In this regard, we show that MCA provides an efficient analytical process to describe observed invasions, which successfully identifies spatial and temporal heterogeneity in the observed rates of spread. By rooting the graph and requiring it to be directed, this approach integrates historical constraints. In addition, the use of minimal weights (similar to the MST) provides a set of parsimonious assumptions to describe the successive dispersal events from one time window to the next. While our data reflect a high-frequency standardized sampling effort, further research is needed to determine whether this method would still be successful with data captured at coarser temporal and spatial sampling grains. The complexities of real landscapes cannot be summarized in any model, but this study shows how an alternative top-down approach based on graph theory can facilitate the ecological analysis of the spread of an invasive species in a new territory.

Publicly available datasets were analyzed in this study. This data can be found at: https://www.sag.cl/ambitos-de-accion/estatus-fitosanitario-de-la-plaga-en-las-regiones-del-territorio-chileno.

SE designed the study, analyzed data, and wrote draft. CS designed the study, prepared the draft, and reviewed the final version. DL designed the study, prepared the draft, and reviewed the final draft. FL designed the study and reviewed the final version. All authors contributed to the article and approved the submitted version.

This study was supported by ANID PIA/BASAL FB0002 and Fondecyt 1211114 to SE, CS and DL and Fondecyt 1221153 to FL.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.