A reaction-diffusion model was developed to simulate the local spread of T. erytreae since its arrival in Portugal for the whole continental area of the country. This type of modelling is commonly used for describing the spatial spread of invasive species (e.g. Shigesada and Kawasaki 1997). The model which incorporates the dispersal of the species and the population’s growth can be expressed using the following Fisher equation:

$\frac{\partial N}{\partial t}=D\left(\frac{{\partial }^{2}N}{\partial x^2}+\frac{{\partial }^{2}N}{\partial y^2}\right)+rN\left(1-\frac{N}{K}\right)$ eq.1

where N is the population density of T. erytreae (km^-2), dependent on time t and spatial location (x,y); D is the coefficient of diffusion (km² /year); r is the population growth rate (year^-1); and K is the carrying capacity (km^-2). The model exhibits a travelling wave with a constant spread rate (C; km/ year) defined by:

$C=2\sqrt{rD}$ eq.2

where C, K and r were all estimated beforehand (see text below), while D was assumed to be homogeneous across Portugal, based on eq. 2 (similarly to Robinet et al. 2017), considering C and r of the area infested by T. erytreae at the end of 2015.

The model described was applied over a grid of 5 km × 5 km resolution, from 2015 to 2021 with a yearly time step.

The population dynamics and spread parameters were obtained, based on previous research studies on the biology and ecology of the psyllid, so that model simulation could be then validated using the independent presence/absence of observed data.

The estimate of the local spread rate C (6 km/year) was based on the maximum flight capacity of T. erytreae determined by Van den Berg and Deacon (1988), i.e. 1.5 km, multiplied by 4, the estimated number of yearly generations that T. erytreae can have in Porto, where T. erytreae was first detected.

The carrying capacity K was estimated for each cell in the model, as the product of the average maximum capacity of T. erytreae per host tree, k_tree and the estimated number of citrus trees in the cell. k_tree was determined using data from Catling (1972), corresponding to the mean number of individuals per citrus tree (adults, nymphs and eggs), observed in a 3-year study, under favourable conditions, using the following formula:

k_tree = Adults + (0.289 × Nymphs + (0.289 × 0.95 × Eggs) eq. 3

where 0.289 is the nymphal survival rate under natural conditions (Caitling 1970) and 0.95 is the egg viability in optimum environmental conditions (Catling 1972; Van den Berg et al. 1991) to obtain k_tree = 2719 adults of T. erytreae per tree.

To estimate the spatial density of citrus trees in orchards throughout the country, we used the data from the 2019 agricultural census (INE 2021). This dataset does not include citrus plants in urban and peri-urban areas, as well as isolated citrus trees in rural landscapes. The dataset from INE (2021) provides the area of citrus orchards per county. The density of citrus trees in orchards was assumed, for simplicity, to be the same all over the territory, i.e. 400 citrus trees/ha, considering 5 m × 5 m per tree. Although citrus-tree spacing may vary between 5 m × 4 m or lower and 7 m × 5 m (Cavaco and Calouro 2005; Vacante and Gerson 2012), we considered for simplicity a median value of 5 m × 5 m. The citrus tree density for each 25 km² grid cell of the model was estimated, based on spatial data from the Land Use and Occupancy Mapping - COS2018, (available at https://www.dgterritorio.gov.pt/) and the data from the 2019 agriculture census (INE 2021), about the area of citrus orchards.

As data on the density of citrus trees in urban and peri-urban areas were not available from the 2019 agriculture census (INE 2021) and considering its possible influence on the dispersal of T. erytreae, we developed an innovative approach to estimate it. We used the spatial data of the COS2018 dataset, to classify the urban and peri-urban areas. Then, we divided these areas into three different classes: Vertical urban areas; Horizontal urban areas; and Discontinuous urban areas (see S1). For each class, the mean density of citrus trees was estimated using Google Street view imagery to survey the number of visible trees in randomly selected polygon areas extracted from the COS2018 spatial dataset throughout the country (Rousselet et al. 2013; Berland and Lange 2017). The estimations were made using the survey counts of citrus trees in each area, weighted against the sample area sizes. Only areas with good image quality were used. We surveyed at least 250 ha for each of the three urban classes considered to provide a confident estimation of citrus tree density. All the surveys were conducted by the first author.

To calculate the growth rate, r, for T. erytreae in Portugal, we used climatic modelled data collected for 30 years (Palma 2017). These data were collected for each centroid of a grid of 25 km × 25 km covering Portugal. The climatic variables considered were the daily mean temperature, daily maximum temperature and the daily minimum relative air humidity. The average daily climatic data were grouped into three periods per month. Each period had 10 days, except the last third of each month, which varied from 8 to 11 days, depending on the month. These periods are henceforth called as “10-day periods”, needed to calculate the T. erytreae survival rate, using the method developed by Catling (1969).

For each 10-day periods, the number of “viable days” (i.e. the number of days above the lower temperature threshold for development) was calculated. Temperatures were estimated, based on the average of the 30 years of the climatic data. A lower temperature threshold of 12.0 °C was considered, based on citrus tree growth not occurring below these temperature values (Webber and Batchelor 1943; Kumar 1977), as well as the inability of T. erytreae larva growth (Catling 1973). An upper temperature threshold was not used since we took into consideration the effect of high temperature and low humidity in the variable weather survival.

Weather Survival (WS%) was calculated for each 10-day period, using the mean Vapour Pressure Deficit (VPD) of the three days with the highest daily values of maximum temperature, using Saturated Vapour Pressure (SVP) (Murray 1967; Green and Catling 1971):

VPD = ((100 – RH)/100) × SVP eq. 4

SVP = 610.7 × 10^{7.5Tmax/(237.3+Tmax)} eq. 5

Weather Survival (WS%) = 0.0308 X₃² – 4.1825 X₃ + 137.7709 eq. 6

where X₃ is the mean value of the VPD in millibars, of the three days with the highest maximum temperatures during the 10-day period.

For the model calculations, we used Weather Mortality (WM%):

WM = 1-(WS/100) eq. 7

For each area, the number of possible yearly generations was then estimated, as the sum of life cycle progress rounded down from each yearly 10-day period from February until the end of September, the most important period of leaf flushing for citrus trees in Portugal (Paiva et al. 2020). Life Cycle Progress was calculated by dividing the average viable days and the estimated total life cycle duration in days for each 10-day period, calculated using the average temperature of the viable days (VD) and the life cycle duration (G), that is the expected total number of days to successfully complete the insect life cycle from egg to adult.

Life Cycle Progress = VD / G eq. 8

The life cycle duration G, in days, was calculated, based on the number of days needed to complete egg incubation (Idays) plus the number of days needed to complete nymphal growth (Ndays), based on the equations proposed by Catling (1969). To this period, we added 5 days of the pre-oviposition period, being the mean of the pre-oviposition time of the species (Van der Merwe 1923; Catling 1973; Van den Berg 1990) and another 10 days to reach the peak of oviposition (Catling 1973).

Idays = 4.9763 + 3.3443 × 0.8452^Tmed-20 eq. 9

Ndays = 16.7974 + 5.2726 × 0.7843^Tmed-20 eq. 10

G = Idays + Ndays + 15 eq. 11

The growth rate at a given time period t in the year, R_t was calculated as:

R_t = f x 0.5 × (1-WM) × (1- 0.8 WM) × (1- 0.3 WM) × (1- 0.15 WM) × (1- 0.075 WM) × 0.289 eq. 12

where f is the female fecundity, estimated by the average number of eggs per female. Due to the different fecundity estimates reported in literature, we tested two different values, i.e. 827 and 327 eggs per female (Moran and Blowers 1967; Catling 1969) in the model. The mean fecundity value was multiplied by egg viability rate, estimated as 0.95 in optimal environmental conditions (Catling 1972). We assumed a sex ratio of 1:1 (Van den Berg 1990). WM, 0.8 WM, 0.3 WM, 0.15 WM, 0.075 WM correspond to the weather mortality for the 1^st, 2^nd, 3^rd, 4^th and 5^th instars, respectively, since larval stages increasingly become more resilient to the adverse environmental conditions (Catling 1972). Finally, 0.289 is the natural survival rate from egg to adult in “perfect conditions”, with the presence of natural enemies, according to Catling (1970).

For the remaining generations, we considered a geometric growth, N_t = N_t-1R_t

The resulting R^G=, ∏R_t where G is the maximum potential generations in each cell per year, will reflect the total growth potential of the species in the area for one year (Pointeau et al. 2021). The intrinsic yearly growth rate of the cell is then calculated as follows:

r = ln(R^G)

We repeated this procedure for each cell of the model.

A cold limitation factor was further generated considering the reproduction of this species is limited by extremely cold temperatures. Cold and long winter periods can hinder population reproduction and limit the species’ capacity to stay in the region throughout the entire year. Thus, we consider that regions having on average less than 2 days of viable days (i.e. with a mean daily temperature above 12.0 °C), for the three consecutive winter months (December, January and February), were not suitable for T. erytreae to establish. The model found that, in these areas, the number of T. erytreae was always 0. The period of three months is based on the maximum longevity of 82 days recorded for the species during winter (Catling 1969).

Modelling long-distance dispersal is always very challenging since it is based on stochastic and relatively rare dispersal events. To model human-mediated long-distance dispersal, we calculated the number of long-distance dispersal events (NB) that occur in each year, which was randomly chosen using the following equation:

NB = 1 + e eq. 13

where e denotes an independent and identically distributed random variable with a Poisson distribution with mean λ = (ln((P/2642)+1/ln(2)) × 5), with P being the number of cells estimated to be infested at a given time t and 2642 is the total number of cells that can be infested in the model.

We used the Poisson distribution as it is a simple discrete distribution that is often used to model jump processes (Hooten and Wikle 2008). The mean λ was defined as a concave increasing function of P, with a minimum value of 2.8 when P = 0 and a maximum value of 6.0 when P = 2642. In this way, the number of simulated long-distance jumps increases with the area infested, taking relatively realistic values compared to the spread pattern observed.

For parameter testing, we made simulations using low-frequency jumps, with low λ = (mean λ) /2 and high frequency jumps with high λ = (mean λ) x 2.

For each long-distance dispersal event, the model randomly chooses a cell that is not yet infested (N < 1), with a growth rate r > 1 and a suitable human population density. The minimum human density threshold (H) allowing the arrival of a long-distance jump was set to 125 habitats/km², using the same threshold considered for the spread of the yellow-legged wasp in France (Robinet et al. 2017), which also delimits urban areas relatively well in Portugal, being mostly in the coastal areas of the country, represented as dark red areas in Fig. 1. A recent study showed that both the distance to roads and urbanisation intensity play an important role in spatial and temporal dynamics of the dispersal of the Asian citrus psyllid, D. citri in California (Bayles et al. 2022). We did not use road network data nor urbanisation intensity due to lack of reliable complete data available for Portugal. Instead, we used human population density (data from 2017 available at https://www.ine.pt/), as it was demonstrated to being a suitable proxy for road traffic in Portugal (Barata 2012).

Figure 1. 

Distribution of the human population density in Portugal (number of inhabitants/km²) in 2017. The human population threshold of 125 habitants/km² was used to characterise locations where long-range distance dispersal could occur in the model.

Cells infested by such long-distance jumps received an arbitrary set of 50 individuals.

A molecular study on the genetic diversity of T. erytreae populations in the Iberian Peninsula suggested the possibility of multiple introductions of the psyllid (Ruíz-Rivero et al. 2021). Two main genetic clusters were observed along the Portuguese coast. Based on these results, we defined two possible scenarios in the model: 1) one single introduction of T. erytreae in the north (Porto in 2014); 2) two introductions, the first in the north (Porto in 2014), and a second in the region of Lisbon (Lisbon in 2017), which has an international port and airport. For the second introduction scenario, we added 1000 individuals to a specific cell in the Lisbon area, where T. erytreae was first detected in the region (data provided by the Portuguese National Plant Protection Authority reports).

We combined estimates of local spread and growth with estimates of long-distance dispersal. We applied these models on a grid that covers Portugal with a spatial resolution of 5 km × 5 km. The simulation began in 2014 in one cell in the Porto region with 1000 individuals for the initial condition of the invasion, which was later discovered already spreading, in 2015.

Since the first detection of T. erytreae in Porto, in 2015, the Portuguese National Plant Protection Authority, has been monitoring the species spread at the parish level, with the deployment of a trapping protocol surrounding the infested areas within a 3 km buffer range, with yellow-sticky traps and additional, monitoring within 10 km range. In addition, they also included reports from citizens and stakeholders, after confirming their authenticity. Reports are publicly available whenever there are newly-infested areas (available at: https://www.dgav.pt/plantas/conteudo/sanidade-vegetal/inspecao-fitossanitaria/informacao-fitossanitaria/trioza-erytreae/). This complete dataset was provided to us by the Portuguese National Plant Protection Authority (DGAV). We compiled it yearly and used it to validate our model, as well as determining the dispersal capacity, from 2015 until the end of 2021. All the parameters used for the model were estimated independently from this presence/absence observed data, allowing for independent validation.

To validate the model and identify the dispersal scenarios that best fit the observed data, the species’ spread was simulated between 2014 and 2021 for different combinations of model parameters (Table 1).

For the simulations considering long-distance dispersal (stochastic sub-model), we ran 300 replicate simulations using randomly generated long-distance events.

We consider that a cell is infested when N ≥ 1 individual.

For the replication of simulations considering long-distance dispersal, we calculated the percentage of simulations that classified each cell as infested. For each cell, if 50% or more simulations predicted the infestation, then the model classifies those cells as infested. Thereafter, for each model, using the simulation data from 2021, we calculated the F1-Score performance criteria (Chinchor and Sundheim 1993), using the following equation:

$F1-\mathrm{Score}·=\frac{2x\left(\text{ precision }x\text{ recall }\right)}{\text{ precision }+\text{ recall }}=\left(\frac{TP}{TP+\frac{1}{2}(FP+FN)}\right)$ eq.14

where TP is the sum of true positives, FP is the sum of false positives and FN is the sum of false negatives. F1-score is the harmonic average of precision and recall.

We compared models with different parameters using the model’s performance criteria, F1-Score. We also calculated standard errors for each model’s F1-Score, using 2000 bootstrap simulations taken from the 300 original simulations of each stochastic model. These were used to obtain the p-values for testing the equality of the F1-Score for pairwise comparisons between two models. P-value was calculated by the inversion of the bootstrap normal confidence interval for the difference in means (Thulin 2021) using the following equation:

$P\text{-value }=2\left(1-2\left(1-{\Phi }^{-1}(F1a-F1b)/(\sqrt{SEa+SEb}\right.\right.$ eq.15

where Φ^-1 denotes the inverse cumulative distribution function of the standardised normal distribution. F1a and F1b are the F1-score of model a and model b, respectively and SEa and SEb are the Standard error values of model a and model b, respectively.

We also calculated the Area Under Cover (AUC) of the receiver operating characteristics plots (Fielding and Bell 1997) and the Youden index of each model (Youden 1950), but they led to the same conclusions as F1-Score.

The tested parameters were the two different female fecundity rate values (Low or High), the inclusion or not of long-distance dispersal events (Yes or No), the frequency of the long-distance dispersal events (Low, Medium, High), the inclusion of the estimated residential urban citrus trees (Yes or No) and the occurrence of a second introduction of T. erytreae in Lisbon (Yes or No) (Table 1).

We compared the estimated local spread rate value against the observed short-range dispersal, based on presence-absence data reports from DGAV. For each year, we calculated the mean least distance between all newly-infested parishes centroids and past infested parish centroids (DP). Infested parishes attributed to long-distance dispersal were removed from the short-range dispersal rate calculation. We identified such parishes, when their DP was higher than 30 km or was higher than the distance between its centroid and a long-distance dispersal parish centroid. Thirty km is an arbitrary distance value that is significantly higher than the yearly estimated flight capacity of T. erytreae and three times the 10 km radius used by the Portuguese Plant Protection Authority for the species monitoring. Finally, with the average value of DP from each year, we calculated the average short-range dispersal rate of T. erytreae in Portugal for each year and in total using all DP values independently of the year of infestation. Furthermore, to calculate the total dispersal capacity to the east and to the south of the country, the main directions the species could spread in Portugal, we used the distance between the infestation origin and the furthest parish towards the east and the south, as was done for the spread pattern of V. velutina (Verdasca et al. 2021). Additionally, we calculated the total infestation area of the infested parishes along the invasion years.

Model simulation running, validation and all statistical analysis were done with the statistical language R version 4.2.0 (R Core Team 2022). Modelling Data and R code are available at https://zenodo.org/record/7096566.

The reports of the Portuguese Plant Protection Authority denote a fast dispersal of T. erytreae in Portugal (Fig. 2). Between 2015 and 2021, the African citrus psyllid was able to spread mostly southwards, along the coastal area of Portugal, covering a maximum distance of about 461 km, between Porto and western Algarve and a cumulative area of about 14,239 km² (Fig. 2). This corresponds to an average of about 65.9 km/year and 2034 km²/year. The dispersal towards the east was only 100 km (14 km/year). However, removing long-distance dispersal events, the observed mean dispersal rate of T. erytreae in Portugal was 7.8 ± 0.3 km/year (Table 2). The estimated short-range dispersal capacity used in our model simulation was 6 km/year, which turned out to be very close to the observed data (ranging from 5.6 to 10.4 km/year) (Table 2).

Figure 2. 

Spatio-temporal representation of Portugal’s invasion by Trioza erytreae, between 2015 and 2021. Elaborated, based on data from the published by the Portuguese Plant Protection Authority.

Our estimates of the number of yearly T. erytreae generations in Portugal varied from 3 to 4 generations per year. The estimated grow rate (r) of T. erytreae in the Portuguese territory was found to be higher along the coast area (Fig. 3). A different female fecundity rate had a major impact on the growth rates estimated, especially in the interior central and southern regions (Fig. 3). The model included a cold limiting factor, portraying areas whose winter was deemed as too extreme for T. erytreae survival, where the growth rate was 0. The cold limited areas are all located in the northern interior part of the country (Fig. 3), where most areas are mountainous and the climate is colder, especially in the winter.

Figure 3. 

The estimated growth index of Trioza erytreae in Portugal, considering two fecundity levels: 327 eggs per female (left); 827 eggs per female (right).

According to the most recent agriculture census of Portugal, there is a total surface of 21,681 ha of citrus orchards in continental Portugal, 74% of which located in the south, in the Algarve Region (INE 2021). We estimated a total of 11,993,645 citrus trees in orchards and 7,427 trees in urban and peri-urban areas (Fig. 4). The estimated citrus-trees density in Vertical, Horizontal and Discontinuous urban areas were 0.37, 3.2 and 5.14 trees per hectare, respectively (see Suppl. material 1: table S1).

Figure 4. 

The spatial representations of the estimated citrus trees density (number of trees/km²) in Portugal. The left map represents the estimation of citrus trees density in Portugal, based on the area of citrus orchards reported in the last agricultural census (INE 2021), while the right map was obtained using both the data from the agricultural census (INE 2021) and our estimates of the number of citrus trees in urban and peri-urban areas, based on Google Street imagery.

The model simulations that included long-distance dispersal fit the observed data better than those that did not (Table 3; see Suppl. material 1: tables S2 for p-values). This is shown by the significantly higher F1-Score between every model with the same combination of parameters besides the long-distance dispersal, independently of the frequency considered, i.e. low, medium or high (e.g. model 6 F1-Score = 0.733 vs. model 30 F1-Score = 0.803, p-value < 0.001; Table 3, see Suppl. material 1: table S2.1). The difference is even greater if the second introduction scenario is not considered (e.g. model 5 F1-Score = 0.583 vs. model 29 F1-Score = 0.801, p-value < 0.001; Table 3, see Suppl. material 1: table S2.1).

Different frequencies of long-distance dispersal events (low, medium and high) were not consistent in model improvement in all parameter combinations (Table 3, see Suppl. material 1: table S2.2).

The inclusion of the estimated urban and peri-urban citrus trees significantly increased the model performance, with significantly higher F1-score values in every model combination (e.g. model 30 F1-Score = 0.803 vs. model 32 F1-Score = 0.686, p-value < 0.001; Table 3, see Suppl. material 1: table S2.3).

The scenario of considering a second introduction was beneficial only for the simulations not using long-distance dispersal, when compared with similar models (e.g. model 6 F1-Score = 0.733 vs. model 5 F1-Score 0.583, p-value < 0.001; Table 3, see Suppl. material 1: table S2.4). The same was not true for model simulations considering long-distance dispersal. The parameter was sometimes not significant (e.g. model 9 F1-Score = 0.791 vs. model 10 F1-Score = 0.790, p-value = 0.97), sometimes beneficial (e.g. model 19 F1-Score = 0.622 vs. model 20 F1-Score = 0.643, p-value = 0.011) and sometimes negative towards model performance (e.g. model 13 F1-Score = 0.804 vs. model 14 F1-Score = 0.795, p-value = 0.006; Table 3, see Suppl. material 1: table S2.4).

Finally, model simulations that considered high fecundity (827 eggs per female) performed better those with low fecundity (327 eggs per female) in 6 out of 8 parameter combinations. In two cases, changing fecundity did not significantly affect model performance (e.g. model 22 F1-Score = 0.800 vs. model 18 F1-Score = 0.794, p-value = 0.126; Table 3, see Suppl. material 1: table S2.5).

Overall, the model simulations with the highest performance were 13, 21, 22 and 30, showing no significant differences (see Suppl. material 1: table S2.6). All these models used long-distance dispersal, high female fecundity and urban citrus trees, but differed in the long-distance dispersal frequency and in considering the second introduction of T. erytreae.

Altogether and considering the temporal evolution between 2015 and 2021, our best model simulations showed a high concordance between the observed and predicted distribution of T. erytreae over the seven years after its detection in Portugal (Fig. 5).

Figure 5. 

Spatio-temporal dynamics of Trioza erytreae spread in Portugal from 2015 up to 2021, predicted using 300 replicate simulations of model simulation #30. The colour gradient represents the probability of each grid cell to be infested according to the model, calculated as the relative number of simulations that predicted the area’s infestation. The border of the observed distribution area of T. erytreae is represented by a black line (elaborated, based on data obtained from DGAV reports) to allow visual assessment of model performance.

Human-mediated dispersal is a well-known documented phenomenon, recognised as a key issue in invasion science (Ricciardi et al. 2017; Bullock et al. 2018; Gippet et al. 2019). Human activities leading to insect dispersal can be divided into three pathways: Contamination, hitchhiking and harvesting (Pergl et al. 2017; Gippet et al. 2019). For the spread of T. erytreae, we believe the major pathways behind human-mediated dispersal of the species would be hitchhiking, as suggested for the invasion of D. citri in southern California (Bayles et al. 2017) and the invasion of the yellow-legged wasp, Vespa velutina in Portugal (Verdasca et al. 2021). In this pathway, adults’ psyllids would be accidentally attached to a vehicle vector, from where they may be transported further away from their flight capacity, increasing the potential dispersal capacity of the species. This dispersal pattern coincides with higher dispersal along the coastline, where the human population is denser. It also reflects the distribution of north-south highways along the coast. Additionally, the movement of infested citrus plants is another possible pathway for the dispersal of T. erytreae in Portugal.

Since its detection, in 2015, the African citrus psyllid was able to spread mostly southwards, along the coastal area of Portugal, covering a maximum distance of about 461 km, between Porto and western Algarve, in seven years, corresponding to an average of about 66 km/year. This dispersal rate is about 4 to 8 times higher than the values reported for other Hemiptera, such as the hemlock woolly adelgid, Adelges tsugae Annand (Adelgidae) (8–13 km/year) and the beech scale, Cryptococcus fagisuga Lindinger (Eriococcidae) (14–15 km/year) (Liebhold and Tobin 2008). However, without considering long-distance dispersion events, the observed mean dispersal rate of T. erytreae in Portugal was 7.8 km/year, similar to the spread rate of the other insect cases (Liebhold and Tobin 2008).

These results highlight that the spread of T. erytreae corresponds to a combination of short-range and occasional long-range dispersal events. This dispersal pattern, called “stratified dispersal” is commonly observed in the spread of invasive insect species (Liebhold and Tobin 2008). For T. erytreae, the inclusion of long-distance dispersal events greatly improved model performance in predicting the observed data (Table 3). In the current study, the large difference between the mean global dispersal rate (66 km/year) and the mean diffusion dispersal rate, estimated excluding the long-distance dispersal events (7.8 km/year), greatly highlights the importance of long-distance dispersal events for the species spread, which often are anthropogenic (Liebhold and Tobin 2008).

Likewise, the predicted and observed spread of T. erytreae along the coastal area of Portugal is also related with the high population density in the area, mostly between Porto and Setubal regions (Fig. 1), where long-distance dispersal events were concentrated. The role of human mediated dispersal was also reported in the yellow-legged hornet’s rapid expansion along the coast of Portugal, attributed to the density of motorways (Verdasca et al. 2021). Nevertheless, motorways density and vehicle traffic are correlated with population density (Barata 2012). Although human-mediated movement of insect life stages is usually the dominant modality of long-distance dispersal, other mechanisms may also be involved, such as the wind, which has not been considered here. For example, Antolínez et al. (2022) showed that the dispersal of the Asian citrus psyllid, D. citri may be influenced by wind speed. On the contrary, wind direction was not found to be a significant factor in an experimental trial on D. citri dispersal conducted by Lewis-Rosenblumet et al. (2015). Nevertheless, Bayles et al. (2017) suggested that the observed spread pattern of the psyllid in California could be related with the prevailing wind direction, but without supporting a definitive conclusion. Future studies should investigate the possibility of assisted dispersal of psyllids in the upper wind, as an additional long-distance dispersal mechanism.

Our model provided contrasting results regarding the hypothesis of additional introductions of T. erytreae in Portugal during its invasion, as suggested by Ruíz-Rivero et al. (2021), based on the genetic diversity of T. erytreae populations. When not using long-distance dispersal in the model, including a second introduction, improved model performance towards predicting T. erytreae spread in Portugal (Table 3). Yet, when coupled with the long-distance dispersal parameter, the effect of a second introduction on model performance was inconsistent. This was shown by the best model simulations (13, 20, 21 and 30), that either used or did not use the second introduction parameter. This is likely due to long-distance dispersal diluting the importance of the introduction in the model since both have a similar impact on the model. In fact, the approach, which was used to include in the model an additional introduction of T. erytreae, was basically based on an input of individuals (1,000) in a defined cell of Lisbon area, which is not much different from a non-random long-distance dispersal event. This outcome further reveals that secondary introductions of invasive species might be frequently misled with long-distance dispersals. Nevertheless, multiple introductions are more than just an addition of individuals to the founding invasive population. They may have an important role in incrementing genetic diversity of the invasive population, compensating the low genetic variability associated with founder effects (Handley et al. 2011). If this increment of genetic diversity is related with new adaptive traits, such as higher fecundity and/or survival rates, then additional introductions are expected to influence the dispersal dynamics of the invasive population. This scenario could be considered in the model, by changing the parameters fecundity and survival. In this respect, it is interesting to note that the higher fecundity value tested was associated with a higher performance of the model.

Urban areas often facilitate the introduction, establishment and spread of non-native species, in biological invasions (Cadotte et al. 2017; Gaertner et al. 2017; Hui et al. 2017; Padayachee et al. 2017). Their green areas, including ornamental trees, public gardens, parks and backyard gardens, may function as stepping stones for non-natives species to disperse and invade agroecosystems (Hui et al. 2017). In Portugal, citrus trees are one of the most common plant species present in urban and peri-urban landscapes, used as ornamental plants in street trees, public gardens, parks, as well as food plants in backyard gardens. Using an innovative method, based on Google Street View imagery, we estimated the spatial distribution of those citrus trees outside orchards and its density according to Urban Areas typology (see Suppl. material 1: table S1). Our results showed that these citrus trees played a very important role in the dispersal of T. erytreae throughout the country (Table 3, see Suppl. material 1: table S2.3). For large areas of the observed distribution of the psyllid in 2021, where citrus orchards are almost non-existent (Fig. 2), the major source of host plants are the citrus trees in the urban and peri-urban areas (Fig. 4). This corroborates the previous claim that ornamental host species can contribute to connecting fragmented citrus-producing lands, as well as act as reservoir areas for the psyllid, especially in the framework of management actions in citrus-producing lands (Van den Berg et al. 1991).

Similarly, a recent study in California (Bayles et al. 2022), where citrus trees are also common ornamental and food plants in urban and peri-urban areas, also pointed out the importance that these trees played in the invasion dynamics of the Asian citrus psyllid, including its spill-over between urban and agricultural habitats (commercial citrus orchards).

We found no recent information in literature regarding female fecundity and no data available from Portugal on this parameter. For the modelling scenarios, we used two values provided in old literature, one considering a high fecundity of 827 eggs/female (Catling 1973) and an alternative one estimating a lower fecundity of 327 eggs/female (Moran and Blowers 1967) (Table 3, see Suppl. material 1: table S2.5). Our modelling results showed significant differences according to this parameter, with the simulations using the higher fecundity performing significantly better than those with the lower fecundity. This outcome evidences the relevance to retrieve this type of basic biological data for the understanding of population invasion dynamics. Regrettably, these biological data are sometimes scarce and with low sampling power. Additionally, biological information may change with populations and differ in the invaded range. As these data are essential for modelling potential spread over several generations, we recommend that efforts should be spent on collecting such biological information on the invaded range of the species.

Globally, our model was able to predict most of the spatio-temporal dynamics of T. erytreae spread quite well, except the recent invasion of the south-western area in 2021, in the coast of Alentejo and west coast of Algarve, for which the model predicted a low colonisation probability of 17% by the psyllid (Fig. 5). This low probability associated to a long-distance dispersal event, into the referred region, is explained by the low human population density in the region. However, a high seasonal touristic flow from the north occurring in this coastal area during Spring and Summer periods was not considered in our model. This large movement of people, including many residents from T. erytreae infested areas, may favour hitchhiking mechanisms of dispersal (Gippet et al. 2019). Nevertheless, the observed presence of the psyllid in the area consisted mainly of small colonies or damage in isolated trees or small groups of trees in backyards and gardens (Amílcar Duarte, University of Algarve and Celestino Soares, DRAPALG pers. com., 2021). Furthermore, even if the model predicts low probability of invasion there, the probability is above 0, so it does not predict absence (Fig. 5). This infestation results from relatively rare and stochastic events, which are difficult to predict with a high probability.

The fast spread of T. erytreae in Portugal occurred despite the efforts carried out by the Portuguese Plant Protection Authority to contain its dispersal and eradicate it. The measures implemented included the delimitation of demarcated areas, being the infested areas plus a buffer zone surrounding the infested areas with trap placement and active monitoring, along with various other control measures within the infested areas. Such control measures have not been accounted for in our spread model simulations and this may explain the reason why some areas in southern and inner Portugal, with relatively moderate predicted infestation probabilities, were not invaded by T. erytreae (Fig. 5). The apparent failure of the model in some parts of the country may result from a certain level of success of the implemented control measures.