We focus on four alien species that are established in Europe: the freshwater jellyfish (Craspedacusta sowerbii), the geranium bronze butterfly (Cacyreus marshalli), the floating primrose-willow (Ludwigia peploides) and the garden lupine (Lupinus polyphyllus). The levels of visual detectability for these species change considerably throughout the year. The freshwater jellyfish is presumed to be native to regions of Asia and has been introduced in most continents of the world. However, its alien distribution remains poorly known, largely because the most visible part of its life cycle involves small medusae that appear for only a few months each year (Marchessaux et al. 2021). The geranium bronze is a small butterfly native to southern Africa and currently invading parts of central and southern Europe. Like most insects, its adult (butterfly) stage has higher visibility due to increased mobility and conspicuous colours. The floating primrose-willow is an aquatic plant native to Oceania and the Americas, with invasive populations in countries of central and southern Europe. This species produces bright yellow flowers, which facilitate its identification amongst surrounding vegetation (Booy et al. 2015). Finally, the garden lupine is a plant native to western North America that is now widespread in many temperate regions of the world, including central and northern Europe. This species produces prominent flower spikes (often violet, but sometimes also pink or white) that greatly facilitate its detection, including via remote sensing (Wijesingha et al. 2020).

Following our previously described framework (Capinha et al. 2024), we collected observation records supported by photographs for each of these species from GBIF. These records were sourced globally (i.e. not limited to Europe) to maximise the range of environmental conditions sampled. To ensure temporal consistency with the environmental predictors (see below), we used georeferenced records from a 7-year period between 2016 and 2022 that included the complete date of observation (day, month and year). GBIF is a major online repository of biodiversity observation data, including from well-known and highly participated citizen-science projects (e.g. iNaturalist.org and observation.org), which typically provide supporting visual media. As the number of records obtained from GBIF for the freshwater jellyfish was low, we also searched for observation records of this species in the USGS Non-indigenous Aquatic Species portal (https://nas.er.usgs.gov) and a few additional sources (see link for full list in the Data Resources section).

We visually checked the photographs supporting each observation record of the four species and kept only those that clearly showed the medusa stage of the freshwater jellyfish, the butterfly stage of the geranium bronze and the flowering stages of both the floating primrose-willow and the garden lupine. Records with images suggesting that the specimens were under human-care (e.g. garden lupine in places showing garden-like features) were excluded. Likewise, we also excluded GBIF records where the observation date was the first day of the month and the observation time was ‘00:00:00’. These are typically records where only the month and year of observation are known and the first day of the month is assigned by default, i.e. the full date of the record may not be precise (Belitz et al. 2023). In total, we obtained 240 records for freshwater jellyfish medusae, 3,879 for geranium bronze butterflies and 1,688 and 10,345 records for floating primrose-willow and garden lupine flowers, respectively (Suppl. material 1: fig. S1).

We collected a time series of global-scale maps representing daily conditions of maximum, minimum and mean temperature, accumulated precipitation, wind speed and accumulated snow. These factors are expected to be drivers of the timing of occurrence of the species life stages of interest, according to previous research (Favilli and Manganelli 2006; Ludewig et al. 2022; Marchessaux et al. 2022). The data were collected from the Global Forecast System (GFS), a large-scale NOAA weather forecast model (https://www.ncei.noaa.gov), covering the period between 15 January 2015 (the first day the data is available) to 31 December 2022. GFS data are provided at 0.25° spatial resolution, for multiple hourly intervals and for each model run that takes place at 00, 06, 12 and 18 UTC daily. The daily-resolution maps were obtained by aggregating the first six-hour values provided by each model run. This means we used the forecast for 00:00–06:00, 06:00–12:00, 12:00–18:00 and 18:00–00:00 UTC, covering the full 24 hours of each day. The temporal immediacy of this period in relation to the timing of model runs results in highly accurate weather forecasts (NOAA 2022). Models employing these data achieve results comparable to those using climate re-analysis data - traditionally for ecological forecasting - like ERA5 (Capinha et al. 2024). Crucially, NOAA also offers real-time access to GFS data, including weather forecasts extending several days into the future, meeting our objective to deliver predictions in real-time and for short-term forecasting.

We implemented a set of procedures to minimise potential spatial and temporal biases in the observation data. Spatial bias refers to unequal numbers of records in distinct regions, which can lead to model responses being ‘dominated’ by the patterns occurring in oversampled regions. Temporal observation bias arises from unequal levels of recording effort within and across years, confounding the actual temporal signal of phenological events.

To address these biases, we followed the procedures we proposed earlier (Capinha et al. 2024). Specifically, and to minimise spatial bias, we kept only one record located in the same grid cell and having the same date of observation. Next, we downsampled observations in oversampled regions. For this purpose, we used a reference grid of 250 km × 250 km squares, for which we identified squares that were upper outliers in terms of record count (i.e. n > third quartile + 1.5 * interquartile range). For these areas, we randomly subsampled a number of records equal to the outlier threshold.

Our framework includes an optional procedure to minimise temporal bias, named ‘benchmark taxa approach’ (Capinha et al. 2024). This procedure uses taxa that maintain a consistent visual appearance throughout the year as indicators of variations in recorder activity (e.g. citizen scientists). We use pine species for this purpose because they are evergreen and have a fairly consistent appearance year-round. Since these taxa do not undergo significant seasonal changes, temporal variations in their recording frequency are likely due to fluctuations in recording effort rather than changes in the taxa themselves. In other words, we expect that the temporal patterns in recordings of these species serve as a proxy for general citizen-science activity levels (see further description of the rationale in Capinha et al. (2024)).

The temporal variation in the frequency of records for these taxa is related to variables expected to mediate levels of recording effort (e.g. days of the week, months of the year and weather conditions) by means of a statistical model such as a generalised linear model. Based on the relationships identified, the temporal biases in records of the phenomenon of interest can be minimised by a subsampling procedure, where records made in periods of higher levels of recording intensity receive a lower probability of being selected for model development. We demonstrated this approach previously and its application delivered similar performance to the models without using it. However, it is not clear if this outcome can be expected in the generality of phenomena. We therefore performed all the analyses using event observation data with this correction (described in Suppl. material: text S1) and without it. The results were similar in both approaches (see Results); therefore, the approach using the temporally corrected data is presented only in the Suppl. material 1.

We next characterised the meteorological conditions preceding each event record. We used a total of 67 features representing multiple features of temperature (e.g. maximum, minimum and mean values, growing degree days and cold accumulation), accumulated precipitation, accumulated snow and mean wind speed for distinct preceding periods ranging from days, weeks and months up to a year (see full list in Suppl. material 1: table S1).

Additionally, we assembled a second set of data aimed at representing the meteorological conditions that are generally available in the location of each of the events (i.e. the background environmental conditions). This was performed using ‘temporal pseudo-absences’ (Capinha et al. 2024), which correspond to records having the same geographical coordinates as event records, but with dates drawn at random within the temporal range of the event data. A total of 12 pseudo-absence records were generated from each event record and each was characterised using the same set of 67 environmental features. The number of pseudo-absence records per observation was determined subjectively, based on preliminary assessments, which indicated a good balance between the diversity of sampled conditions and the volume of data produced.

Prior to model fitting, we tested for the presence of multicollinearity amongst the predictors. For this purpose, we measured their variance inflation factor (‘VIF’) and excluded any predictor with a VIF value above 10 (Gareth et al. 2021). We then trained three machine-learning algorithms: random forests (RF), boosted regression trees (BRT) and generalised linear models with lasso regularisation (GLM-lasso), in distinguishing the conditions associated with the phenological stages of interest and those represented by the temporal pseudo-absences. These algorithms were selected because they are commonly used for ecological modelling and prediction and often rank amongst the best performing, including when transferred to new spatial settings (Zhang et al. 2019; Valavi et al. 2023).

The implementation of these models was performed in R (R Core Team 2024), using the ‘randomForest’ package for RF (Liaw and Wiener 2002), the ‘dismo’ package for BRT (Hijmans et al. 2017) and ‘glmnet’ for GLM-lasso (Hastie et al. 2021). We optimised several parameters within these functions to improve model fitting. Random forests used 2,000 individual trees (instead of the default 500) to increase the chances of the relatively large number of predictor variables and samples being adequately represented in the final ensemble. For the BRT models, the number of trees in each ensemble was automatically determined by the `gbm.step` function of the ‘dismo’ package, with a tree complexity of 3 (allowing for interactions amongst predictor variables) and a learning rate of 0.005. Additionally, for all three algorithms, we addressed the class imbalance in the data (i.e. one event observation record for every 12 pseudo-absence records), which could skew models towards over-predicting the dominant class (pseudo-absences). For GLM-lasso and BRTs, this was made through a weighting parameter (‘weights’ for GLM-lasso and `site.weights` for BRTs), assigning class-proportional weights to each sample. For Random Forests, the adjustment was made using the ‘sampsize’ parameter of the ‘randomForest’ function, ensuring that each individual tree used an equal number of event observation records and temporal pseudo-absences.

To evaluate the predictive performance of the models, we used a leave-one-year-out cross-validation procedure. This involved excluding the data from one year for model calibration and using it to assess the predictive ability of models trained on the remaining years. The procedure was iterated so that the data from each year served as an evaluation set. To measure the models’ performance, we used the Boyce Index, initially proposed for species distribution models (Hirzel et al. 2006). The Boyce Index measures the correlation between the frequency of observation records and predicted probability values. In the context of our work, a strong positive correlation implies that event observations are concentrated in periods predicted with higher probability values, whereas correlations near zero indicate predictions akin to those expected from a model that assigns probability values randomly throughout the year. We preferred this metric over discrimination-based metrics, like the commonly-used area under the receiver operating characteristic curve (Pontius and Parmentier 2014), because these latter metrics would involve assessing the misclassification of temporal pseudo-absences. However, in the context of our work, these pseudo-absences capture the environmental variation that is available, including periods that are potentially suitable for the occurrence of our focal phenological stages. This becomes particularly critical for models of phenological stages with extended seasons of occurrence, where temporal pseudo-absences during suitable periods are more likely. Consequently, this would likely result in an underestimation of the models’ true predictive performance due to the misclassification of pseudo-absences.

We performed the Boyce Index calculations using the ecospat.boyce function from the ecospat R package (Broennimann et al. 2015), employing the Pearson correlation coefficient as a metric. The measurements were performed for the three algorithms and for an ensemble of these predictions corresponding to their average.

Identifying regions with similar year-to-year phenological patterns (“phenoregions”; White et al. (2005)) allows defining areas where surveillance programmes for invasive species could be implemented in the same periods. To identify these regions, we used the phenology predictions from RF models at 5-day intervals for the 7-years period (from 2016 to 2022), across Europe. We used the predictions from this algorithm because it consistently performed well, ranking as the best or second-best across all species (see Results). To minimise the risk of extrapolation errors, i.e. making predictions for conditions far from those represented in the training data, we only made predictions for regions with Köppen climate classes (from Beck et al. (2018)) where each phenological stage had at least five observation records. Although simple, this method allows us to discriminate between regions with higher or lower propensities for model extrapolation. While this approach was effective for our purposes, future research could benefit from incorporating finer-scale assessments of extrapolation.

We applied a k-means algorithm to cluster regions based on the temporal variation in predicted values, using the ‘elbow’ method to determine the optimal number of clusters (Syakur et al. 2018). For each identified region, we then calculated the mean and standard deviation of inter-annual variation.

Overall, the predictions from models (Fig. 1, Suppl. material 1: figs S2, S3) achieved high performances. The predictions of imago-stage of the Geranium bronze attained a cross-model median correlation with the timing of observations of 0.95. GLM-Lasso achieved the highest median performance (r = 0.97), followed by RF (r = 0.96), multi-model ensemble (r = 0.95) and BRT (r = 0.94) (Fig. 2). Pairwise Kruskal-Wallis tests indicate non-significant differences in these performances (i.e. p ≥ 0.05; Suppl. material 1: table S2). Likewise, the flowering timings of the floating primrose-willow and the garden lupine were also well captured by the models, with median cross-model correlations of 0.94 and 0.91, respectively. For the floating primrose-willow, the multi-model ensembles, RF and GLM-Lasso models had the best performance (median r = 0.94), with BRT performing slightly the worst (median r = 0.88). However, as before, these differences were not statistically significant (Kruskal-Wallis p ≥ 0.05; Suppl. material 1: table S2). In contrast, for the garden lupine, RF displayed a statistically significant superiority (median r = 0.96; Kruskal-Wallis p ≤ 0.05; Suppl. material 1: table S2), over remaining algorithms. BRT was second best (median r = 0.92), followed by multi-model ensemble and GLM-Lasso (median r for both = 0.88), but with no statistically significant differences amongst these (Kruskal-Wallis p ≥ 0.05; Suppl. material 1: table S2). Predictive performance for the medusae stage of the freshwater jellyfish was the lowest, but still reasonable – with a median correlation of 0.73 across models (Fig. 2). Predictions from the ensemble approach demonstrated a substantially higher performance (median r = 0.82), followed by RF (median r = 0.79), GLM-Lasso (r = 0.64) and BRT (r = 0.42), but these differences were also not statistically significant (Suppl. material 1: table S2).

Figure 1.

Examples of daily predictions obtained from the modelling approach. These represent the probability of occurrence for each modelled phenological stage for four species (imago-stage of the Geranium bronze [a]; medusae of the freshwater jellyfish [b], flowering of the floating primrose-willow [c] and flowering of garden lupine [d]), for 1 July 2023. Predictions were obtained using random forests, the best performing algorithm, trained with observational data corrected for spatial bias.

Figure 2.

Results of Boyce Index corresponding to Pearson correlation values between predicted probabilities of event occurrence and the frequency of event observation records of the imago stage of the Geranium bronze butterfly (a), medusae of the freshwater jellyfish (b), the flowering of floating primrose-willow (c) and the garden lupine (d). The boxplots represent the variation of correlation values assessed for 7 years (2016 to 2022), using three modelling algorithms (boosted regression trees, BRT; generalised linear models with lasso regularisation, GLM-Lasso; random forest, RF) and an ensemble of previous algorithms (Ensemble), trained with observation data corrected for spatial bias.

Relevantly, models trained with data addressing both spatial and temporal biases showed similar predictive performances as those addressing only spatial biases (Suppl. material 1: fig. S4). Specifically, the predictions of imago-stage of the Geranium bronze retained the highest cross-model median correlation with the timing of observations (0.95), followed by the garden lupine (0.91), the floating primrose-willow (0.90) and the freshwater jellyfish (0.67). Performances for the Geranium bronze were highest for GLM-Lasso, RF and multi-model ensemble (median r = 0.95), followed by BRT (median r = 0.90). For the floating primrose-willow, the best predictions were achieved by multi-model ensembles and RF (median r = 0.91), followed by GLM-lasso (median r = 0.90) and BRT (median r = 0.87). For the garden lupine, the best performances were achieved by RF (median r = 0.97), followed by multi-model ensembles and BRT (median r = 0.90) and by GLM-Lasso (median r = 0.88). Finally, for the freshwater jellyfish, the best performances were achieved by RF and GLM-Lasso (median r = 0.68), followed by multi-model ensembles (median r = 0.66) and BRT (median r = 0.44). As for models using observation data corrected for spatial bias only, the only statistically significant difference amongst algorithms was verified for the garden lupine, with RF significantly outperforming remaining ones (Suppl. material 1: table S2).

We also assessed the performance of days-ahead forecasts across Europe over a 10-month period (Fig. 3; Suppl. material 1: fig. S5). For the freshwater jellyfish and the floating primrose-willow, the observation records used for model evaluation during this period were relatively limited (25 and 42, respectively), while they were more abundant for the Geranium bronze (1,187) and the garden lupine (126 records from a random sample of 1,000 processed records). Overall, the forecasts demonstrated a high correlation with the timing of observed life stages modelled, with values generally close to 0.8 or higher (Fig. 3; Suppl. material 1: fig. S5). Additionally, there was no apparent pattern of decreased model performance as the forecast horizon extended (e.g. from 1 day ahead to 9 days ahead).

Figure 3.

Boyce Index values for forecasts made 1, 3, 6 and 9 days in advance. These values correspond to the Pearson correlation coefficient between predicted probabilities of event occurrence from July 2023 to March 2024 and the frequency of event observations recorded during the same period. The values are reported for three modelling algorithms—boosted regression trees (BRT), generalised linear models with lasso regularisation (Lasso) and random forest (RF)—as well as an ensemble of these algorithms (Ensemble), all trained with observation data corrected for spatial bias.

Using the predictions of the random forest algorithm trained with observation data corrected for spatial bias, we identified five phenoregions for the Geranium bronze butterfly in Europe (Fig. 4a), with the peak occurrence ranging from early June to late November in the southernmost region (Fig. 4b). The remaining regions have peak occurrences between August and September, with northern latitudes having a shorter season of occurrence and the opposite being true for southern latitudes. Recorded observations of butterflies of this species are concentrated in the summer (Fig. 4c), especially in the mid- to higher latitude regions, in agreement with the predictions.

Figure 4.

Regional patterns of predicted and observed timings of the emergence of the imago stage in the Geranium bronze butterfly (panels a–c), the occurrence of medusae in the freshwater jellyfish (panels d–f) and the flowering phases of floating primrose-willow (panels g–i) and garden lupine (panels j–l). The maps display areas having similar phenological dynamics (‘phenoregions’), based on daily projections at 5-days intervals from 2016 to 2022. Grey areas represent Köppen climate classes where fewer than five observations of the species were made. These classes were not included in the analysis to minimise the risk of model extrapolation. Time series depict the inter-annual mean probabilities of occurrence of each event, along with their ranges (grey shading), for each region throughout the year. Histograms show the monthly frequency of effectively observed occurrences of each life stage within each phenoregion. Predictions and phenoregion delineations were made also for areas where the species have not yet been recorded, leading to observation records being absent from the histograms for certain regions.

Predictions of the timing of occurrence of medusae of the freshwater jellyfish were clustered into four phenoregions, peaking between late August and September. However, southern regions exhibit higher probabilities of occurrence over substantially longer periods and, conversely, shorter periods are predicted for northern regions. Medusae observations take place in the months of predicted peaks, except for the mid-latitude region covering most of Central Europe, where they concentrate in September and November — a period when the predicted probabilities are already declining.

For the floating primrose-willow, four regions were identified (Fig. 4g), with the southern regions showing earlier and more extended flowering periods (Fig. 4h). The timing of observational data matches well with the predictions for the region including most of Mediterranean Europe, although very few records are available for the other regions (Fig. 4i).

Predicted timings of flowering for the garden lupine were classified into four phenoregions (Fig. 4j), with probabilities peaking from early May to late June and occurring earlier in southern regions (Fig. 4k). Observational data also largely agree with the predicted patterns, especially for the two northernmost regions, which have most records (Fig. 4l).

Importantly, identified phenoregions and associated temporal patterns were strikingly similar to those obtained, based on models calibrated with temporally unbiased observation data (Suppl. material 1: fig. S6). This further suggests the robustness of predictions to temporally varying levels of recording effort.

The authors have declared that no competing interests exist.

No ethical statement was reported.

C.C. was supported by Portuguese National Funds through Fundação para a Ciência e a Tecnologia through support to CEG/IGOT Research Unit (UIDB/00295/2020 and UIDP/00295/2020) and by the EuropaBON project, funded by European Union’s Horizon 2020 Research and Innovation Programme under grant agreement No 101003553.

CC: Conceptualisation, Methodology, Validation, Formal analysis, Resources, Data Curation, Writing - Original draft, Project administration, Funding Acquisition. AC and AM: Writing - Review and Editing.

César Capinha https://orcid.org/0000-0002-0666-9755

Spatial weather data used for this work are publicly available online from NSF NCAR Research Data Archive (RDA) n (https://rda.ucar.edu/datasets/d084001/). Event observation data are available from USGS Non-indigenous Aquatic Species (https://nas.er.usgs.gov/queries/factsheet.aspx?SpeciesID=1068) and GBIF with DOIs: https://doi.org/10.15468/dl.3fve6q, https://doi.org/10.15468/dl.9drr85, https://doi.org/10.15468/dl.h5amhh, https://doi.org/10.15468/dl.krbguc, https://doi.org/10.15468/dl.2kbnxv, https://doi.org/10.15468/dl.ycjsyz, https://doi.org/10.15468/dl.7q6rke, https://doi.org/10.15468/dl.b9sze9 and https://doi.org/10.15468/dl.uh8apz. Additional data sources and R code are publicly available on Zenodo (https://doi.org/10.5281/zenodo.13847953).