Occurrence data

We obtained a global occurrence data set for T. kirbyi from GBIF (Global Biodiversity Information Facility, www.gbif.org) and Biodiversidad Virtual (www.biodiversidadvirtual.org), a citizen science platform with over 2 million georeferenced records for the Iberian Peninsula. Additionally, we cross-checked these data with the Odonata Database of Africa [31, 32] to verify the information for that particular region. Occurrences were then filtered to remove records prior to 1970 and unverifiable data points in water far from the shore. To avoid autocorrelation, only one presence per pixel of the climatic variables at 1 × 1 km was retained. Although we assumed a possible bias towards better sampled areas, such as Europe, than to less sampled areas, such as Central Africa, we kept all the raster cells with occurrences to accurately follow the hypothesized colonization. The final data set included 1,717 occurrences (S2 Fig).

Models were evaluated with occurrences from Odonata Database of Africa, which were not used at all in the modelling procedure.

Species distribution modelling

We estimated the climate-based potential distribution of T. kirbyi using the biomod2 package [33] and six models: Generalized Linear Model (GLM), Generalized Additive Model (GAM), Artificial Neural Network (ANN), Classification Tree Analysis (CTA), Random Forest (RF), and Maximum Entropy (MaxEnt). The final ensemble model was based on the average of 60 individual models (i.e., 10 iterations × 6 models).

For the bioclimatic variables, we used the 19 variables of the WorldClim 2.1 data set (http://www.worldclim.org) at 1 × 1 km cell size. This data set provides information on seasonality trends, averages, and extreme values of temperature and precipitation from 1970 to 2000, and these bioclimatic variables are widely used to predict the potential distribution of species by reflecting temperature and precipitation patterns and their variation, factors that are overall ecologically meaningful for species [34–36]. We visually inspected the variables and eliminated those with anomalous patterns in the study area, which were mostly due to the difficulty of getting reliable interpolated precipitation values across areas with low precipitation, such as the Sahara Desert or the western Mediterranean [bio 8, bio 9, bio 15, bio 18 and bio 19, see 37]. To establish a set of uncorrelated climatic variables, we intersected the remaining variables with 10,000 randomly selected points in the extent of the study area, performed a correlation analysis, and removed one of the variables from each pair with a Pearson correlation value > 0.7 (S3 Fig), leaving the one perceived to have higher relative ecological importance and predictor effectiveness. The final data set included Mean Diurnal Range (bio2), Isothermality (bio3), Maximum Temperature of Warmest Month (bio5), Minimum Temperature of Coldest Month (bio6), Precipitation of Driest Month (bio14), and Precipitation of Wettest Quarter (bio16).

The construction of background and pseudoabsences was based on a previous simple environmental coverage model with only occurrences, performed with the range between the maximum and minimum values of each bioclimatic variable [38]. Areas that had all their values within the maximum and minimum range of each variable were considered habitable for T. kirbyi and were used to establish the background points, whereas areas that did not fulfil at least two of these variables were used to establish pseudoabsence points [39, 40]. During the development of this preliminary model, we observed which bioclimatic variables discriminated the greatest amount of area between habitable and non-habitable areas (i.e., areas outside the maximum and minimum range) and considered them in the final selection of variables (as those that barely discriminate between habitable and non-habitable areas provide less information in the subsequent distribution model).

Presence, pseudoabsence, and background data were split into 75% training and 25% testing data sets to generate an extrinsic Area Under the receiver operating characteristic Curve (AUC) evaluation for the final models that was independent of the intrinsic AUC evaluation of each individual model generated by biomod2. A final ensemble model was obtained by a weighted averaging of the six models in each replica. Weights for each AUC model were calculated, using only those with an AUC > 0.7 [41]. Finally, the ensemble models were evaluated through the extrinsic AUC test with the 25% testing data sets. The continuous ensemble model was binarized to a presence/absence model using the cut-off threshold (0.619) that maximized the True Skill Statistic (TSS). A post hoc evaluation of the final model on Africa was carried out using the occurrences of T. kirbyi included in the Odonata Database of Africa [31], which were not used in the training or testing data sets used for model generation. For this purpose, we inspected how many points (and their suitability) were either inside or outside the presence/absence threshold suitability value. We selected a number of pseudoabsences equal to the number of occurrences in Africa (828) to perform an AUC test and repeated the process 100 times to obtain an average AUC value to evaluate the ensemble model for Africa.

Variable selection and species distribution models were performed in the R 3.5.0 and RStudio 1.1.453 environments. ArcGIS for Desktop 10.3 (www.esri.com) was used to generate background and pseudoabsence points and model maps.

Distribution and dispersal analyses

To assess the environmental characteristics of the distribution area of T. kirbyi in the western Mediterranean, we converted the continuous suitability values to a presence/absence raster using the cut-off threshold value that maximized TSS (0.619). At each level, we generated 10,000 random points to run ANOVA tests and generate boxplots using presence/absence as a dependent variable. We did not perform this analysis nor the subsequent ones at a global scale due to the high environmental variability of the global study area (e.g., the desert areas of the Sahara can bias the models).

We also studied which variables were most important for the suitability of T. kirbyi in the colonized area. The relative importance of each variable and algorithm was calculated by the biomod2 package. We also performed an independent conditional Random Forest analysis (2,500 trees, 4 variables) with the 20,000 presences and absences previously generated, followed by a simple regression analysis between the resulting most important variables and suitability to assess the significance effect (R²).

To assess the role of summer temperature peaks on dispersal speed and spread, we used the annual mean temperature change with respect to the data from 1880 to 1920, i.e., the annual thermal anomaly, using data at the latitude of the Iberian Peninsula [42]. We calculated annual and cumulative increases in the distribution area of the species (as minimum convex polygons) and compared them against the annual and cumulative thermal anomalies for each year since 2010 (the scarcity of observations between 2007 and 2009 precluded the inclusion of these years) using the Mann–Whitney–Wilcoxon U test. We use both the area (raw variable) and the transformed (linearized) area to perform these analyses.

Finally, we calculated the maximum observed dispersal distance, defined as the absolute distance between the highest latitude point(s) of consecutive years, since the northernmost point(s) can be interpreted as the colonization avant-garde of this dragonfly. We chose from 1 to 5 points for each pair of consecutive years to calculate the dispersal distance mean, standard deviation and, finally, the total dispersal distance for 2010–2022.

Potential distribution and dispersal of T. kirbyi

Trithemis kirbyi thrives in three biogeographical realms, the Palearctic, Afrotropical and Indomalayan [43], covering a wide distribution in Africa and Eurasia. Its affinity for warm and dry climates makes the western shores of the Mediterranean Sea an ideal location for its northward expansion into hitherto vacant areas that experience an increase in temperature (Fig 3) and dryness. Rapidly rising temperatures, with never-before-seen peaks, may permit this northward expansion to occur at an unprecedented rate (Fig 5).

Although T. kirbyi shows a wide but scattered distribution, the model shows areas of maximum suitability connecting the known populations. The Indian area, considered as the original distribution area of the species, is connected to Africa through the coast of Pakistan, Iran, and the Arabian Peninsula. From this last region, the species’ distribution area branches towards southern Africa through the Rift Valley, and towards West (tropical) Africa through an almost continuous high suitability belt. According to the model, the North African coast, the Iberian Peninsula, and the remaining Mediterranean coastal areas are also highly suitable for the species. Although the Sinai Peninsula and coastal areas of the eastern Arabian Peninsula are connected to high suitability areas in the southern Arabian Peninsula, T. kirbyi has not been recorded in the region from the Sinai Peninsula to Tunisia except for an isolated locality in Egypt [20]. Thus, an African-European expansion through the Tibesti and Hoggar mountains, as implied by Boudot and Kalkman (20), or through the west coast of Africa, as suggested by our suitability model and the presence of the species in scattered localities in Mauritania, seems more plausible than a colonization from the east (e.g., Anatolian Peninsula), as supported by the lack of records from southeastern Europe. A west African expansion into Europe would also fit with the documented western to eastern colonization of the Mediterranean islands [24, 25].

Trithemis kirbyi has a high dispersal potential and can disperse an average of 430 km per year. A mark-recapture study of dragonflies that has been carried out by the Biodiversity Monitoring Group of the Universidad Complutense de Madrid in a peri-urban area of Madrid (Spain) since 2013 supports the long range dispersal capacity of the species. In 2015, several individuals of T. kirbyi were marked in this long-term study, representing the first record of this dragonfly in the central Iberian Peninsula [44]. One individual that was first captured and marked in Madrid city (40.4467N 3.7254W) on July 7 was observed approximately 38 km away in La Pedriza (40.7721N 3.8687W) on July 28 (Victor Salvador, personal communication). Its high dispersal capacity is also supported by the fact that T. kirbyi colonized the Balearic Islands from the Iberian Peninsula, demonstrating that it can cross great distances over the sea: T. kirbyi reached the island of Ibiza from the Iberian Peninsula (127 km) and subsequently Mallorca, where it was observed in a locality 130 km from the nearest occurrence in Ibiza. The species crossed similar distances (100–150 km) during its colonization of the Sicilian Channel islands from the mainland or among the islands [24, 27], as happened in Sardinia in 2003 [26].

Colonization of an exotic species driven by climate change

Trithemis kirbyi is an exotic species in Europe and North Africa [45] that has naturally colonized the Iberian Peninsula through Africa, and is currently expanding into northern European territories. Species occupy areas that fulfil their abiotic (e.g., minimum temperatures that influence the activity of the species) and biotic (presence of water bodies to complete their life cycle) requirements, and that can be reached [the "BAM diagram", 46, 47]. Human-facilitated dispersal combined with climatic niche shifts [48] are responsible for the introduction and expansion of most invasive alien species. In the case of T. kirbyi, however, there is not an a priori niche shift: the species occupies the same climatic and environmental niche reported for the African populations (such as Mediterranean climate areas as the south of South Africa or North Africa). It also differs from other invasive alien species by being a strong disperser, capable of migrating to areas far from well-established populations. Climate change is generating ideal conditions for T. kirbyi in southern Europe, and its dispersal pattern reveals its northern expansion in relatively well-defined waves that coincide with thermal anomaly peaks (Fig 5). This is an example of how previously unsuitable areas have become potentially colonizable or invadable [49], and how climate change can favour biological invasions [50], in addition to natural colonization following successful dispersal events. Further studies could investigate this niche conservatism of the T. kirbyi to reinforce our hypotheses.

Case studies of dispersion and establishment in insects are scarce and rely on accurate monitoring of biodiversity [51–54]. Despite the importance of some insect species as agricultural pests or disease vectors, their monitoring is generally not a high priority for governmental administrations. Citizen science platforms such as Biodiversidad Virtual, Observation (https://observation.org), and iNaturalist (https://www.inaturalist.org), however, have proven highly useful for gathering up-to-date, expert-validated information on species introduced into a new territory, and can serve as monitoring systems that can be used by both academic researchers carrying out scientific studies and policymakers responsible for implementing health and environmental measures [55–58].

Conservation implications

Trithemis kirbyi has been able to colonize areas in the Iberian Peninsula that have experienced a high increase in temperature and aridity. This region represents one of the most arid in the southwestern Mediterranean, and models predict it will be severely affected by climate change [59]. Although the species is thermophilic and a generalist, which facilitate its colonization of warmer and drier areas with respect to its abiotic requirements, increased aridity implies a lower availability of water masses, a resource this species needs to complete its life cycle. Changes in water mass availability in increasingly arid environments will impact the future distribution of the species, as is being currently observed with its increasing scarcity in the Sahara Desert.

Rising temperatures caused by climate change have been identified as a factor leading to the expansion of some insect distributions [e.g., 60]. Trithemis annulata is another example of a dragonfly expanding into Europe following temperature increases in the 1990s [20], much like our model of T. kirbyi [61]. The impacts of these two exotic species of Trithemis on native communities have not been studied so far, thus we do not know if they could displace native species or if they are simply new additions to existing plastic and generalist communities. So far, there is no evidence to suggest that T. annulata and T. kirbyi have become invasive species in the western Mediterranean, but they certainly have become naturalized. Again, citizen science platforms can play a crucial role in detecting the potential impact of these species on native fauna almost in real time [62]. In fact, during the writing of this manuscript, T. kirbyi was also found in Belgium during a high temperature episode, reinforcing the idea of the importance of this type of platforms and data, as well as evaluating and validating our model (suitability of this area of Belgium = 0.6) [63].

Final remarks

The steady and marked increase in temperature since the 1980s, combined with extremely intense heat waves in the last decade [28, 64], has turned a hitherto inhospitable ‘cold’ region into a new colonizable area with ideal conditions for alien species from southern latitudes. In the last decade, several North African arthropods have been detected in southern areas of the Iberian Peninsula, including Xylocopa pubescens Spinola, 1838 [65], Azanus ubaldus (Stoll, 1782) [66], and Cheilomenes propinqua (Mulsant, 1850) [67]. Although these species colonized the area naturally, similar to T. kirbyi, many exotic species have been introduced into Europe from areas much farther away owing to the international trade in goods and inefficient quarantines [68]. We demonstrate natural colonization by an insect that has, so far, shown no signs of harm to natural autochthonous ecosystems; however, our study shows that the road is paved for large-scale invasions that will cost money for the direct control measures they should trigger (e.g., measures against agricultural pests and vector-borne diseases). They will also cause large economic losses due to a decline in biodiversity complexity and associated ecosystem services (e.g., pollination or quality in freshwater ecosystems).