The workflow uses occurrence data from the Global Biodiversity Information Facility (GBIF) (Fig. 1). GBIF is a large biodiversity data infrastructure that supports the FAIR Data Principles and is widely used for scientific purposes (Heberling et al. 2021; Ivanova and Shashkov 2021). The workflow uses an occurrence cube format, which aggregates occurrences along spatial (e.g. Extended Quarter Degree Grid Cell (EQDGC)), temporal (e.g. yearly) and taxonomic (e.g. species) dimensions (Desmet et al. 2023), making the data suitable for modelling and computational analysis (Oldoni et al. 2020). The cube is processed using the b3gbi R package (Dove 2026). This ensures standardised input data for further analysis and verifies that the data format is correct. To minimise the influence of sampling bias, multiple records of a species from the same site (grid cell, i.e. a single spatial unit on the Earth’s surface, for example, E033S02BB of EQDGC) within a given year are converted to binary presence (1) or absence (0) records for all computations (Boria et al. 2014).

The status of species (alien or native) constitutes an Essential Biodiversity Variable for monitoring biological invasions (Latombe et al. 2017). It is necessary to confirm that the species included in the analysis are alien within the region of interest before computing impact indicators. Several published datasets now facilitate the verification of species status. For example, Gómez-Suárez et al. (2025) developed a harmonised dataset of species status using the standardised SinAS workflow (Seebens et al. 2020), which integrates and harmonises global sources such as the Global Register of Introduced and Invasive Species (GRIIS) and the World Checklist of Vascular Plants (WCVP). In addition, reliable sources such as peer-reviewed literature and technical reports can be used to determine alien status.

The workflow uses EICAT data as evidence of the impact of alien species. Impacts of alien species, based on EICAT assessments are classified into Minimal Concern (MC), Minor (MN), Moderate (MO), Major (MR) or Massive (MV), according to the severity of their negative effects on native species in the recipient community. Specifically, a species is assigned MC if there is no reduction in the performance of individuals of a native species; MN if the performance of individuals is reduced, but there is no decrease in native species population size; MO if the native species’ population is reduced; MR if the impact has led to local extinction, but is reversible; and MV if the impact has led to naturally irreversible extirpation or extinction. The impacts are also classified according to the 12 mechanisms described by the IUCN, which are competition, predation, hybridisation, transmission of disease, parasitism, poisoning/toxicity, biofouling or other direct physical disturbance, grazing/herbivory/browsing, chemical impacts on ecosystems, physical impacts on ecosystems, structural impacts on ecosystems and indirect impact through interactions with other species. Additional information, such as the location of the impact, is also reported (IUCN 2020a). EICAT data are available through the Global Invasive Species Database (GISD; https://www.iucngisd.org/gisd/). However, the EICAT assessments, hosted in GISD, are not yet fully open or FAIR, as they currently do not provide persistent identifiers and are not fully accessible for bulk or machine-readable extraction. For species without published assessments, users can conduct and add their own assessments within this workflow.

For computational purposes, the semi-quantitative (ordinal) EICAT categories must first be transformed into numerical scores. The transformation depends on the assumed encoding scale for the impact categories, such as a linear or exponential relationship and on how the MC category is interpreted. Specifically, MC may be treated either as zero impact, because no impact on native species has been recorded or as a non-zero category reflecting some degree of concern. We considered three encoding schemes, which can be used depending on the assumptions of the study. The first, MC = 0, MN = 1, MO = 2, MR = 3 and MV = 4, assumes a linear relationship and interprets MC as zero impact (Hagen and Kumschick 2018). This scheme is consistent with the interpretation of MC as the absence of detectable effects on native species. The second, MC = 1, MN = 2, MO = 3, MR = 4 and MV = 5, also assumes a linear relationship, but treats MC as a non-zero category to represent some degree of concern, despite no documented impact on native species (e.g. Jansen and Kumschick (2022)). The third, MC = 1, MN = 10, MO = 100, MR = 1000 and MV = 10,000, assumes an exponential relationship between categories and scores, while likewise treating MC as non-zero. This scheme gives disproportionate weight to higher-impact categories (Suppl. material 1: appendix 1, table SS1). There are many other possible schemes not considered in this framework, but these may be included in further versions of the software after careful methodological evaluation.

To produce impact indicators, we used a hierarchical aggregation framework with three levels: species, site and region (a whole study area, for example, South Africa). First, because species are often assigned multiple impact categories across study locations and impact mechanisms, we aggregated these records to obtain a single species-level impact score intended to represent the overall documented impacts of each alien species outside its native range. For example, Acacia dealbata has been classified as MR in the Drakensberg, South Africa, through structural impact, but as MN in the Biobío Region, Chile, through chemical impact (see Jansen and Kumschick (2022)). To reconcile such variation, we converted impact categories to scores and aggregated them within species using one of three approaches: Maximum, Mean or Sum of Maximum across Mechanisms (Table 1, Fig. 2).

We then used these species-level scores to estimate site-level impact. As multiple alien species can co-occur at a site, the impacts of all species recorded at a site were aggregated using one of five methods proposed by Boulesnane-Guengant et al. (2025): Precautionary, Precautionary Cumulative, Mean of Means, Mean Cumulative and Cumulative. These methods differ in how impacts are aggregated within and across species (Table 1, Fig. 2).

Finally, to estimate the impact of alien species at the scale of the study region, we aggregated site-level impact values across all sites within each year. The regional impact was, therefore, computed as the spatial sum of site-level impacts using the same five aggregation approaches as those used for the site-level indicator: Precautionary, Precautionary Cumulative, Mean of Means, Mean Cumulative and Cumulative (Table 1, Fig. 2).

To estimate the uncertainty associated with the impact indicators, we used the bootstrapping method. For each year in the occurrence cube, the procedure re-samples individual records (rows) from the cube (including impact scores) with replacement and then re-computes the indicator for each re-sampled dataset. The original cube is treated as the population and re-computing the indicator across multiple bootstrap samples yields a distribution of possible indicator values. The standard error and confidence interval are derived from this distribution (Suppl. material 1: appendix 2). We implemented the bootstrapping and calculated the confidence interval using the dubicube R package (Langeraert and Van Daele 2026). The uncertainty measure is currently available only for the regional impact indicator.

To estimate how sensitive the regional impact indicator is to individual species and to identify species that disproportionately influence the regional impact indicator, we used a leave-one-species-out cross-validation (LOSO-CV) technique. This technique evaluates the extent to which the regional impact indicator is driven by individual species by systematically recalculating the indicator after excluding one species at a time and comparing the resulting values to those obtained from the full occurrence cube. The percentage deviation between the indicator calculated after omitting a species and the full-data indicator provides a measure of the contribution of the left-out species to the regional impact indicator (Suppl. material 1: appendix 3). We implemented the LOSO-CV using the dubicube R package (Langeraert and Van Daele 2026).

R has become the primary programming language for ecological data analysis over the past decade (Lai et al. 2019; Gao et al. 2025). Therefore, we developed an R package, impIndicator (Yahaya et al. 2026), to facilitate the use of the workflow by ecologists and other stakeholders. The package implements the workflow described in sections: Occurrence data, Impact data, Impact indicators and Uncertainty and sensitivity (Fig. 1) and is fully open source. Its source code is available through the Biodiversity Building Blocks for Policy (B3) GitHub repository (https://github.com/b-cubed-eu/impIndicator). B3 is a European Union Horizon-funded project that aims to standardise access to biodiversity data for proactive policy-making (https://b-cubed.eu/). Each version release is automatically archived on Zenodo (https://doi.org/10.5281/zenodo.15052675). The package can be installed directly from the B3 R-universe repository (https://b-cubed-eu.r-universe.dev/builds):

We demonstrated the workflow using Acacia species in the continental regions of the Iberian Peninsula (Spain and Portugal combined), South Africa and California. We produced species-level, site-level and regional impact indicators for Acacia in these three regions, as well as uncertainty and sensitivity analyses for the regional impact indicators, using the impIndicator R package. We downloaded the GBIF occurrence cube for Acacia in Spain, Portugal, South Africa and California via the GBIF website (GBIF.org 2026). The downloaded cubes have an EQDGC Level 2 (0.25° by 0.25° grid cell) spatial resolution and yearly temporal resolution. We used the EQDGC because it has been widely used in African atlases and is well-suited to global analyses (Harrison et al. 1997; Larsen et al. 2009). In contrast, the European Environment Agency (EEA) reference grid is only available for Europe. The yearly interval is the basic format of the cube, which can be changed depending on the study. To test the sensitivity of the impact indicators to spatial and temporal resolution, we performed analyses at finer and broader spatial resolutions (0.125° by 0.125° and 0.5° by 0.5° grid cells, respectively) and finer and broader temporal resolutions (1/2-year and 2-year intervals) (Suppl. material 1: appendix 4). We used the EICAT assessment data from Jansen and Kumschick (2022), because these are currently more detailed than those available through the GISD. The data provide the Acacia species name, categories and the mechanism of the impact reported in each source. For all the analyses in this study, we used the default encoding scheme in impIndicator for the transformation of impact category to impact score; MC = 0, MN = 1, MO = 2, MR = 3 and MV = 4 because it is consistent with the original EICAT categories.

As an example, we first prepare acacia_cube_IBR, eicat_acacia and iberian_sf, representing the occurrence cube for Acacia species, EICAT data and the Iberian Peninsula’s shapefile, respectively (see Suppl. material 1: appendix 5). As the spatial object iberian_sf represents only the mainland of the Iberian Peninsula, the impact indicator is computed exclusively for the mainland. Consequently, only species occurrences located within the Iberian mainland area are included in the computation.

The impIndicator R package provides the function compute_species_indicator to compute the impact indicator per species using the Mean method (Table 1). The function returns a dataset containing the aggregation method used, the number of species, species names and the impact values for each species. The results demonstrate variability in the impacts of different species across regions. Overall, eleven unique Acacia species with impact evidence occur across the three regions: Acacia baileyana, Acacia cyclops, Acacia dealbata, Acacia decurrens, Acacia elata, Acacia longifolia, Acacia mearnsii, Acacia melanoxylon, Acacia pycnantha, Acacia salicina and Acacia saligna. Specifically, Acacia dealbata, Acacia mearnsii and Acacia baileyana have the highest impact values in the Iberian Peninsula, South Africa and California, respectively (Fig. 3). The output can be displayed in the console using the print function. The following shows an example of the code and output for the species-level impact indicator of Acacia in the Iberian Peninsula.

We used the compute_site_indicator function to compute the site-level impact indicator using the Mean Cumulative method. The function returns a dataset containing the aggregation method used, the number of sites, the number of species and the impact value for each site. We visualised the site-level impact indicator for 2014, 2019 and 2024 in each region, highlighting their spatial patterns of impact (Fig. 4). The following shows an example of the code and output for computing the site-level impact indicator for Acacia in the Iberian Peninsula.

We used the compute_regional_indicator function to compute the regional impact indicator using the Mean Cumulative method. Confidence intervals were estimated from the bootstrap replicates using the percentile method. The function returns a dataset containing the method used to compute the indicator, the number of sites, the number of species and the regional impact value together with its confidence interval. The resulting plot shows the regional impact indicator of Acacia for each region, indicating that regional impact values are higher in the Iberian Peninsula and South Africa than in California (Fig. 5). The following shows an example of the code and output for computing the regional impact indicator for Acacia in the Iberian Peninsula.

We applied the LOSO-CV Percentage Error analysis to the regional impact indicator for Acacia in the three regions using the Mean Cumulative method. The LOSO-CV Percentage Error quantifies how individual species contribute to the impact indicators across the three regions. For example, Acacia dealbata contributed approximately 55% to the regional impact indicator in the Iberian Peninsula in 2010 (Fig. 6). Specifically, Acacia dealbata, Acacia mearnsii and Acacia baileyana made the largest contributions to the regional impact indicator in the Iberian Peninsula, South Africa and California, respectively.

The impact values varied consistently with spatial and temporal resolution. Specifically, the range of site-level impact values across the three regions decreased as spatial resolution increased (Fig. 4, Suppl. material 1: appendix 4, figs S2, S5). In contrast, species-level and regional impact values increased consistently with increasing spatial resolution (Figs 3, 5, Suppl. material 1: appendix 1, figs S1, S3, S4, S6). This pattern arises because finer spatial resolution aggregates fewer species within individual sites (grid cells), resulting in lower site-level impact values. At the same time, finer resolution increases the number of occupied sites, leading to higher species-level and regional impact values.

The species-level and regional impact values decreased with coarser temporal resolution (Figs 3, 5, Suppl. material 1: appendix 1, figs S7, S9, S10, S11), whereas site-level impact values showed little variation across the three temporal resolutions (Fig. 4, Suppl. material 1: appendix 1, figs S8, S12). This is because finer temporal resolution partitions occurrence records into shorter time intervals, reducing the number of species and sites aggregated within each interval and consequently lowering species-level and regional impact values.

Rapid and reliable monitoring information is required to enable effective management of threats to biodiversity (UNEP 2022). In this study, we developed a reproducible workflow and an accompanying R package, impIndicator, to quantify species-level, site-level and regional impact indicators that track temporal and spatial changes in the impacts of alien species. Together, the workflow and software provide a practical framework for translating occurrence and impact data into indicators that provide information for evidence-based management and policy decisions related to biological invasions. By packaging the workflow into an accessible tool, the impIndicator R package facilitates the application of these indicators by researchers, conservation practitioners and policy stakeholders.

The values of the impact indicators generated by the workflow are driven by changes in the distribution of alien species within the study region. For example, the regional impact indicator increases when an alien species spreads to additional sites or when a new alien species is recorded in the study region (Fig. 5). Consequently, the impact indicators reflect both spatial expansion and the accumulation of impacts through time. As the indicators are sensitive to the spatial and temporal resolution of the underlying data, their absolute values should not be interpreted in isolation. Instead, they are most informative when compared across species, sites, regions or time using the same methodological settings. For example, the regional impact value for 2022 can only be interpreted as indicating a higher or lower impact relative to other years. Accordingly, the primary value of these indicators lies in supporting comparative assessments rather than providing absolute measures of impact. Nevertheless, additional factors, such as sampling completeness and species detectability, may influence the estimated indicator values. We discuss the implications of sampling completeness for the Acacia case study later in this section.

Species can be assigned their maximum recorded impact score, consistent with the approach adopted by the IUCN to support a precautionary framework for biodiversity conservation (Kumschick et al. 2024). However, this practice has been questioned by Cassini (2023), who argued that the maximum impact approach may not be appropriate for all purposes and departs from classical measures of central tendency, such as the mean, commonly used in ecological analyses. Moreover, the maximum method ignores the variation in magnitude, thereby obscuring the ranking of species for prioritisation (Volery et al. 2021). By offering flexible aggregation options, including the mean impact score and the sum of maximum impact across mechanisms, our workflow addresses these concerns by enabling users to select aggregation methods aligned with their specific research or management objectives. This flexibility also facilitates comparison with and integration of existing impact frameworks. For example, the impact risk maps described by Kumschick et al. (2025) are equivalent to our site-level impact indicators, while our workflow further extends this approach by allowing multiple options for aggregating species impact scores within sites (Boulesnane-Guengant et al. 2025).

As each indicator method employs a different aggregation strategy, the resulting indicators can produce substantially different estimates of impact. The choice of method should be guided by the objectives of the analysis and the assumptions underlying each approach. For example, the precautionary site-level indicator assigns the maximum impact score amongst species present at a site. As a result, when all species considered have been classified as causing high impacts (e.g. Major or Massive) at least once, the indicator may show little or no variation amongst sites because each site is assigned the same maximum score. In contrast, methods such as the Mean Cumulative indicator incorporate the contributions of all species and their average impacts, allowing greater variation in site-level impact values and providing a representation of cumulative invasion pressure (Boulesnane-Guengant et al. 2025).

Furthermore, the site-level impact indicator serves as a visual and analytical tool to represent the intensity of biological invasions across different parts of an area (e.g. Fig. 4). By enabling spatial comparisons, such as between provinces, states or conservation areas, it highlights hotspots and areas at risk of impact caused by alien species. These spatial data are useful for prioritising management actions, coordinating restoration projects and fostering cross-regional collaboration to address alien species impacts effectively (Potgieter et al. 2022).

The regional impact indicator offers a nuanced representation of temporal changes in alien species impacts at local, regional or global scales (e.g. Fig. 5). By tracking the increase and decrease of ecological threats over time, the indicator can provide insights into the dynamics of alien species impacts, helping assess whether current management practices are effective or need adjustment. This temporal perspective can also support more targeted resource allocation and proactive interventions to mitigate biodiversity loss and ecosystem degradation.

The LOSO-CV Percentage Error evaluates the sensitivity of the regional impact indicator to the removal of individual species, highlighting how strongly the regional impact depends on particular species and providing insight into the potential consequences of species control or removal. In contrast, the species-level impact indicator is intended for direct comparison amongst species, enabling the ranking and prioritisation of species, providing information for control strategies and advancing research on alien species’ ecological roles and adaptation patterns.

The Iberian Peninsula is one of the European hotspots for invasions by Acacia species (Botella et al. 2023). This prominence has spurred research on their management in recent years (Marchante et al. 2024), reflecting both ecological impacts and socioeconomic costs associated with these invasions. In Portugal, the significant environmental damage caused by Acacia species has led to strict regulatory measures: all Acacia species are currently subject to legal restrictions (Decreto-Lei no. 92/2019). Although Acacia species are widely distributed across mainland Portugal and the north-western and coastal regions of Spain, our results revealed substantial spatial variability in their impact across the Peninsula. Thus, a widespread distribution does not necessarily translate into uniformly high impact. This heterogeneity is especially relevant for management planning. Comprehensive control of all Acacia species across the entire Iberian Peninsula would require considerable financial and logistical resources. Therefore, species-level and site-level impact indicators provide a more strategic approach, enabling managers to prioritise the most impactful species and the most affected sites. Such targeted prioritisation can improve cost-effectiveness and enhance the overall efficiency of invasion management efforts in the region.

The California Invasive Plant Council (Cal-IPC) is a non-governmental organisation that works to protect California’s environment from invasive plants through science and policy. The Cal-IPC classifies species according to their impacts into five categories (High, Moderate, Limited, Alert and Watch), based on state-wide spread as a proxy for their ecological impact. Ten Acacia species have been listed in their inventory as having a significant current impact and/or a risk of becoming invasive (The Cal-IPC Inventory – California Invasive Plant Council 2026). Most of these species are classified as Watch due to their potential future impact, leading to the conclusion that Acacia species currently have relatively minor ecological impacts in California compared with other Mediterranean-climate regions (Rundel 2023). This aligns with our results, which showed lower indicator values in California than in the Iberian Peninsula and South Africa (Fig. 5).

Since the EICAT impact scores assigned to each species remain constant through time, the increase in impact values, for instance between 2018 and 2021 in South Africa (Figs 3b, 5), is mainly due to additional occurrence records, which may reflect either genuine spread to new sites or increased recording effort. However, the apparent rate of spread inferred from our indicators differs from that reported by Kotzé et al. (2025), who surveyed South African sites during 2007–2008 and 2016–2023 using observers in low-flying aircraft and estimated an increase of only 10.6% between the two survey periods. This discrepancy is most likely attributable to differences in sampling effort and data sources. In particular, the dataset used by Kotzé et al. (2025) is not included in the GBIF occurrence records analysed here and the increase in records between 2018 and 2021 is likely associated with the incorporation of iNaturalist observations into GBIF. Similar effects may also influence the impact indicators in the other study regions.

Apart from the sensitivity of the impact indicators to spatial and temporal resolution discussed above, we highlight that the indicators should be interpreted with caution, because observed trends may reflect changes in sampling effort in addition to genuine changes in species distributions. This is a common limitation of biodiversity occurrence datasets, arising from factors such as (i) delays between field observations and their publication; (ii) the exclusion of datasets that are not openly available; (iii) temporal changes in data availability as large datasets become publicly accessible (e.g. the incorporation of iNaturalist records) and (iv) non-systematic sampling. One way to reduce these limitations is to incorporate locally collected occurrence data, particularly those based on systematic or high-intensity sampling, even when they are not yet available through GBIF. Ideally, this should be complemented by more rapid and comprehensive integration of such datasets into GBIF.

Another limitation of this workflow is the current lack of EICAT impact assessments for most alien taxa, although ongoing initiatives are expanding the taxonomic coverage through coordinated assessment project, and few assessments are available in the scientific literature (e.g. Evans et al. (2016); Canavan et al. (2019); Volery et al. (2021); Jansen and Kumschick (2022)). At present, these assessments are neither fully open nor fully FAIR, requiring users to provide their own assessment data for many taxa and regions. The workflow also assumes that impact magnitudes recorded in one location are transferable to other locations (Hayes and Barry 2008), although this assumption may not hold in all ecological contexts, such as island ecosystems. Since most available impact assessments are not specific to the study region, the impact indicators presented here represent potential impacts, which may differ from the realised effects of alien species within the study region or at individual sites. As additional impact assessments become available, it will become possible to filter for impacts relevant to the region of interest or other ecological criteria. However, such geographically explicit impact data remain unavailable for most invaded regions.

Feedback will be collected via the issue tracker on the impIndicator GitHub repository (https://github.com/b-cubed-eu/impIndicator/issues). Future versions of the package will extend uncertainty calculation and visualisation to species-level and site-level impact indicators. Recent advances in deep learning species distribution models (SDMs) (Ryckewaert et al. 2025) provide an opportunity to extend the workflow to project the potential impact of alien species across space and time. Such models combine species occurrence records with environmental predictors to estimate species distributions, including in areas with limited sampling or undetected occurrences (Pearson et al. 2007; Algar et al. 2009). This approach could help reduce biases associated with uneven or unknown sampling effort, including delays between field observations and their subsequent publication through GBIF, thereby improving the robustness of the resulting impact indicators.