From Remote Sensing to Species Distribution Modelling: An Integrated Workflow to Monitor Spreading Species in Key Grassland Habitats

Remote sensing (RS) has been widely adopted as a tool to investigate several biotic and abiotic factors, directly and indirectly, related to biodiversity conservation. European grasslands are one of the most biodiverse habitats in Europe. Most of these habitats are subject to priority conservation measure, and several human-induced processes threaten them. The broad expansions of few dominant species are usually reported as drivers of biodiversity loss. In this context, using Sentinel-2 (S2) images, we investigate the distribution of one of the most spreading species in the Central Apennine: Brachypodium genuense. We performed a binary Random Forest (RF) classification of B. genuense using RS images and field-sampled presence/absence data. Then, we integrate the occurrences obtained from RS classification into species distribution models to identify the topographic drivers of B. genuense distribution in the study area. Lastly, the impact of B. genuense distribution in the Natura 2000 (N2k) habitats (Annex I of the European Habitat Directive) was assessed by overlay analysis. The RF classification process detected cover of B. genuense with an overall accuracy of 94.79%. The topographic species distribution model shows that the most relevant topographic variables that influence the distribution of B. genuense are slope, elevation, solar radiation, and topographic wet index (TWI) in order of importance. The overlay analysis shows that 74.04% of the B. genuense identified in the study area falls on the semi-natural dry grasslands. The study highlights the RS classification and the topographic species distribution model’s importance as an integrated workflow for mapping a broad-expansion species such as B. genuense. The coupled techniques presented in this work should apply to other plant communities with remotely recognizable characteristics for more effective management of N2k habitats.


Introduction
Grassland ecosystems are a key spot of biodiversity worldwide [1]. They are characterized by noticeable plants and animal diversity and have a crucial role in species conservation. For these reasons, most of these communities in Europe are considered priority habitats of community interest under the European Community Habitats Directive (92/43/EEC) [2]. However, changes in agro-pastoral practices and land-use evolution mean that grasslands are rapidly disappearing and are nowadays among Europe's most threatened ecosystems [3]. The land-use change in these ecosystems is mainly due to the Pastoral activities strongly affect the structure and composition of grassland habitats over time. During the past century, the evolution of land use from sheep breeding into cattle and equine pasturage has favored the spreading of competitive and not palatable species. Moreover, the study area is crucial for biodiversity conservation, as one of the most extensive upland plains in Europe. Indeed, despite its limited size (~ 9400 ha), it is extremely rich in species composition [40,41]. For these reasons, the area falls within Special Protection Areas (SPAs) [42] and in Sites of Community Interest (SCIs) [43].

Remote Sensing Dataset and Pre-Processing
Brachypodium genuense (DC.) Roem. et Schult. is a graminoid perennial herbaceous species and is one of the most widespread dominant herbaceous species in Europe. In the central Apennines, this species follows distinct phenological phases [33] exploited in this study to identify it through satellite images.
The distribution of B. genuense was assessed using multispectral images of S2 satellites (Copernicus Program) at 10 m per pixel [44] accessed by ESA's Open Access portal (https://scihub.copernicus.eu, accessed on). The satellite images were chosen considering the phenological stages of the B. genuense species in the central Apennines [33]. Three phenological phases were considered: the beginning of vegetative growth (Green up-June 2019), the flowering and fruiting (Maturity-July 2019), the beginning of leaf yellowing and complete drying (Senescence-September 2019) [45].
"Level 1C" (TOA-top-of-atmosphere reflectance) was chosen by filtering images with low cloud coverage (0-20%). Then, cloud masking [46] and the atmospheric Pastoral activities strongly affect the structure and composition of grassland habitats over time. During the past century, the evolution of land use from sheep breeding into cattle and equine pasturage has favored the spreading of competitive and not palatable species. Moreover, the study area is crucial for biodiversity conservation, as one of the most extensive upland plains in Europe. Indeed, despite its limited size (~9400 ha), it is extremely rich in species composition [40,41]. For these reasons, the area falls within Special Protection Areas (SPAs) [42] and in Sites of Community Interest (SCIs) [43].

Remote Sensing Dataset and Pre-Processing
Brachypodium genuense (DC.) Roem. et Schult. is a graminoid perennial herbaceous species and is one of the most widespread dominant herbaceous species in Europe. In the central Apennines, this species follows distinct phenological phases [33] exploited in this study to identify it through satellite images.
The distribution of B. genuense was assessed using multispectral images of S2 satellites (Copernicus Program) at 10 m per pixel [44] accessed by ESA's Open Access portal (https: //scihub.copernicus.eu, accessed on 24 January 2021). The satellite images were chosen considering the phenological stages of the B. genuense species in the central Apennines [33]. Three phenological phases were considered: the beginning of vegetative growth (Green up-June 2019), the flowering and fruiting (Maturity-July 2019), the beginning of leaf yellowing and complete drying (Senescence-September 2019) [45].

Data Collection
Presence/absence data for classification analysis were collected in the field integrating with some absences remotely gathered (e.g., Catonica, et al. [53], Congedo et al. [54], Copernicus [55], Google Earth [56]). A random stratified sampling based on elevational gradient and land-cover was used. The ground truth data were collected close to the exact date when the imagery used for the classification was acquired (2019-07). Moreover, the samplings were made in homogeneous areas to avoid mixed pixel signal [51]. A total of 80 sample points (presences/absences) ( Figure 1) were collected within the study area. The presence sampled points were collected in 40 homogeneous patches of 100 m 2 , where B. genuense cover is higher than 90%. The same number (40) of absence points was also collected as follows: 30 absences were directly observed in the field (absence of the target species in the patches), while the rest were acquired from remotely sensed layers (see above) in sites where the target species cannot occur (e.g., forest canopies, bare soils, screes, urbanized).

Classification Analysis by Machine Learning
To capture the extents of B. genuense in the study area, we performed a binary classification using presence/absence data [57]. Three machine learning classifiers were tested: Random Forest (RF) [58], Support Vector Machine (SVM) [59], and Classification and Regression Trees (CART) [60]. The "superclass" function within the RStoolbox package [61] in R environment [62] was used to perform all the classification analyses. To assess the model accuracies, we first varied the number of input variables (with/without ancillary data). At the same time, we modulated the mtry parameter for RF, the C (Cost) parameter for SVM (svmRadialCost), and maxdepth for CART, based on the work of Li et al. [63].

Accuracy Assessment
We assessed overall classification performances by cross-validation and the Kappa through the "superclass" function of the RStoolbox package [61]. We used an iterative k-folds (k = 10) cross-validation approach; for each run, the dataset was split into 80% train and 20% test [61].
To estimate the area where the classification model can be reliably applied, we used the method proposed by Meyer and Pebesma [64], estimating the area of applicability (AOA). The dissimilarity index (DI), based on the minimum distance from the training data in the multidimensional predictive space, was calculated to derive the AOA. The predictors are weighted by their respective importance in the model [65]. The AOA is derived by applying a threshold based on the DI of the training data. The DI of the training data is calculated concerning the cross-validation strategy used for training the model. The CAST package [66] was used to estimate the AOA in the R environment. This function estimates the DI and the derived AOA of spatial prediction models by considering the distance of new data (i.e., a raster stack of spatial predictors used in the models) in the predictor variable space to the data used for model training.
The spatial association of the classified B. genuense patches were tested using the local indicators of spatial association (LISA) analysis. This analysis provides for each observation an indication of the extent of significant spatial clustering, and the sum of LISAs for all observations can be considered as a global indicator of spatial association [67].

Conversion of Plant Associations into NATURA 2000 Habitats
The plant associations in the study area were extrapolated from "La vegetazione di Campo Imperatore (Gran Sasso d'Italia)" [41]. The plant associations were transformed in N2k habitats according to the Interpretation Manual of European Union Habitats-EUR28 [2] and the Interpretation Manual of Italian Habitats [68]. The complex vegetation systems like those of the sinkholes that occur in a small area cannot be reported as one habitat type; thus, they were named "Mosaic". To evaluate the spread of B. genuense on N2k habitats, we calculated the covered area of B. genuense obtained by classification analysis. All spatial processes and geographic analyses were managed through QGIS 3.10.11 [69]. The raster obtained from the classification process has been polygonized and intersected with another vector layer of the N2k habitats [41]. This post-modelling analysis was performed to investigate the spatial relationship between the B. genuense and the N2k habitats.
To assess how each habitat was affected by the spread of B. genuense, the patches' cumulative contribution in each habitat and the differences between "realized" cover with "expected" cover were calculated. The realized cover is the real cover of the target species in each habitat. The expected cover is an equal repartition of B. genuense cover proportional to each habitat area; thus, the target species results equally distributed among the habitats. The "expected" cover in each habitat (Expected Cover (i) ) was calculated by using the following formula: AreaH (i) is the total area of each habitat, TotAreaBra is the total area covered by B. genuense, and TotArea is the entire study area.
We also investigate how the distribution of B. genuense is associated with the ecotones among habitats using the "moving windows" approach. The moving windows approach is widely used to detect the boundaries or ecotonal areas; indeed, a peak in standard deviation is associated with the transition zone [70][71][72]. Each pixel that shows a standard deviation higher than zero was identified as ecotones. The presence/absence data of B. genuense detected in these ecotonal areas were compared with each habitats' core areas. To test if B. genuense is associated with habitats' ecotones, three different buffer areas were used: 3 × 3, 5 × 5, and 7 × 7 [51].

Topographic SDMs
Species distribution models (SDMs) were used to investigate the topographic drivers of the investigated species. They are numerical tools that combine observations of species occurrence with environmental estimates [24].
The topographic variables considered were elevation, eastness, northness, slope, topographic wetness index (TWI), roughness, terrain ruggedness index (TRI), topographic position index (TPI), and solar radiation. All these variables were extrapolated from the digital elevation model (DEM) at 10 m per pixel, free available for the Abruzzo Region. Autocorrelation among variables was tested using Pearson correlation coefficient, excluding variables with Pearson R > 0.75 [73,74]. The occurrences data were extrapolated from the RS classification. Moran's test was performed in R studio using the usdm package [75] to assess the spatial autocorrelation among occurrences. We reduce the high spatial autocorrelation among occurrences filtering out 95% of data based on topographic similarity. Occurrences were classified into 100 classes based on topography similarity. Within each class was sampled the 5% of points, obtaining a total of 884 presence points.
Three sets of 1000 pseudo absences were randomly sampled in all areas where the classification does not identify the B. genuense occurrences. For the species distribution model were used two different algorithms: Generalized Linear Models (GLMs; type Remote Sens. 2021, 13,1904 6 of 16 = "quadratic", interaction level = 2) and Generalized Boosting Model, also known as Boosted Regression Trees (BRT; number of trees = 10000, interaction depth = 3, crossvalidation folds = 10). The choice of these techniques permitted us to explore responses from different classes of models, ranging from more classical statistical techniques (GLMs) to machine learning-oriented approaches (BRT) [25,76]. GLMs are based on parametric linear functions [73,77]. BRT combines the regression-tree and boosting algorithms to optimize predictive performance from an ensemble of trees sequentially fitted focusing on residuals from the previous iterations [78]; this technique generally results in high discrimination performance fit of accurate function [23]. The SDMs for the investigated species was performed using BIOMOD2 [79] and ECOSPAT [80].
All the methods used in this work are summarized in the following flowchart ( Figure 2).
autocorrelation among occurrences filtering out 95% of data based on topographic similarity. Occurrences were classified into 100 classes based on topography similarity. Within each class was sampled the 5% of points, obtaining a total of 884 presence points. Three sets of 1000 pseudo absences were randomly sampled in all areas where the classification does not identify the B. genuense occurrences. For the species distribution model were used two different algorithms: Generalized Linear Models (GLMs; type = "quadratic", interaction level = 2) and Generalized Boosting Model, also known as Boosted Regression Trees (BRT; number of trees = 10000, interaction depth = 3, cross-validation folds = 10). The choice of these techniques permitted us to explore responses from different classes of models, ranging from more classical statistical techniques (GLMs) to machine learning-oriented approaches (BRT) [25,76]. GLMs are based on parametric linear functions [73,77]. BRT combines the regression-tree and boosting algorithms to optimize predictive performance from an ensemble of trees sequentially fitted focusing on residuals from the previous iterations [78]; this technique generally results in high discrimination performance fit of accurate function [23]. The SDMs for the investigated species was performed using BIOMOD2 [79] and ECOSPAT [80].
All the methods used in this work are summarized in the following flowchart ( Figure 2).

Classification Output and Accuracy Assessment
All the classification models used showed improved performances by adding ancillary data to the individual spectral bands deriving from S2 (Supplementary Materials, Table S1). RF showed the best performance for mapping the target species compared to SVM and CART (Supplementary Materials, Table S1); therefore, the RF model was chosen for the subsequent analyses.
The coverage map of B. genuense (Figure 3), obtained through the RF classifier, showed a large spread of this species within the study site. The classification algorithm (RF) detected 776.98 ha of B. genuense area out of the 9402.84 ha total study area (8.26% of the total area).

Classification Output and Accuracy Assessment
All the classification models used showed improved performances by adding ancillary data to the individual spectral bands deriving from S2 (Supplementary Materials, Table  S1). RF showed the best performance for mapping the target species compared to SVM and CART (Supplementary Materials, Table S1); therefore, the RF model was chosen for the subsequent analyses.
The coverage map of B. genuense (Figure 3), obtained through the RF classifier, showed a large spread of this species within the study site. The classification algorithm (RF) detected 776.98 ha of B. genuense area out of the 9402.84 ha total study area (8.26% of the total area).  Before performing the classification process, we estimated the spatial areas for which the classification model should provide reliable predictions through the AOA calculation. Forecasts outside the AOA (Supplementary Materials, Figure S2) were strongly limited in the study area's north-western sector. In the presence-absence context, the best result obtained (RF) on the Campo Imperatore upland plain provided a value of overall accuracy of 94.79% with a Kappa of 88.0% (Table 1). More details about the classification results were reported in Supplementary Materials (Table S3). In general, the RF classifier showed relatively consistent overall performance at the study site. The spatial association analysis (LISA) on patches of B. genuense shows that most patches have a high clustering value. The general model in the study area is "clustered." Only in the central and western parts of the study area, some patches are more dispersed (see Supplementary Materials, Figure S4).

Conversion of Plant Associations into Natura 2000 Habitats
As observed by Figure 4, the dominant vegetation cover in the study area is represented by high altitude grasslands and meadows (N2k habitats 6170, 6210*, 6230*). Before performing the classification process, we estimated the spatial areas for which the classification model should provide reliable predictions through the AOA calculation. Forecasts outside the AOA (Supplementary Materials, Figure S2) were strongly limited in the study area's north-western sector. In the presence-absence context, the best result obtained (RF) on the Campo Imperatore upland plain provided a value of overall accuracy of 94.79% with a Kappa of 88.0% (Table 1). More details about the classification results were reported in Supplementary Materials (Table S3). In general, the RF classifier showed relatively consistent overall performance at the study site. The spatial association analysis (LISA) on patches of B. genuense shows that most patches have a high clustering value. The general model in the study area is "clustered". Only in the central and western parts of the study area, some patches are more dispersed (see Supplementary Materials, Figure S4).

Conversion of Plant Associations into Natura 2000 Habitats
As observed by Figure 4, the dominant vegetation cover in the study area is represented by high altitude grasslands and meadows (N2k habitats 6170, 6210*, 6230*).  Table S5), the plant associations, overlapped with B. genuense cover map, were divided into five habitats. The 6170, 6210* and 6230* are grassland habitats, the 8210 are screes habitats, and the 4060 represent shrubby habitats.

Overlay Analysis Output
Following the overlay analysis (Figure 4), B. genuense classified in the study area (8.26% of the total area) was found in five different types of habitats ( Table 2). The differences between the "realized" and "expected" cover of B. genuense in each habitat show as the most invaded habitat is the 6210* (Semi-natural dry grasslands and scrubland facies on calcareous substrates-Festuco-Brometalia-* important orchid sites) with 255 ha more than expected. Moreover, 74.04% of the total B. genuense distribution was detected within this habitat. On the other hand, the less affected is the Alpine and subalpine calcareous grasslands (6170), with 3.45% B. genuense cover. The other habitats (6230*, Mosaic and 8120) have slight differences between "expected" and "realized" cover. These habitats have shown a cumulative contribution of 5.17% in 6230* (Species-rich Nardus grasslands, on siliceous substrates in mountain areas and submountain areas, in Continental Europe), 1.90% in "Mosaic" (habitats 6170/6210*/6230*), and 0.48% in 8120 (Calcareous rocky slopes with chasmophytic vegetation).  Table S5), the plant associations, overlapped with B. genuense cover map, were divided into five habitats. The 6170, 6210* and 6230* are grassland habitats, the 8210 are screes habitats, and the 4060 represent shrubby habitats.

Overlay Analysis Output
Following the overlay analysis (Figure 4), B. genuense classified in the study area (8.26% of the total area) was found in five different types of habitats ( Table 2). The differences between the "realized" and "expected" cover of B. genuense in each habitat show as the most invaded habitat is the 6210* (Semi-natural dry grasslands and scrubland facies on calcareous substrates-Festuco-Brometalia-* important orchid sites) with 255 ha more than expected. Moreover, 74.04% of the total B. genuense distribution was detected within this habitat. On the other hand, the less affected is the Alpine and subalpine calcareous grasslands (6170), with 3.45% B. genuense cover. The other habitats (6230*, Mosaic and 8120) have slight differences between "expected" and "realized" cover. These habitats have shown a cumulative contribution of 5.17% in 6230* (Species-rich Nardus grasslands, on siliceous substrates in mountain areas and submountain areas, in Continental Europe), 1.90% in "Mosaic" (habitats 6170/6210*/6230*), and 0.48% in 8120 (Calcareous rocky slopes with chasmophytic vegetation). The results obtained using the moving windows techniques identify the ecotonal area among habitats. Compared to the B. genuense cover between the core and the ecotonal areas, a slight difference was observed (≈1.6%). All three different moving windows size (3 × 3, 5 × 5, 7 × 7) give similar results, and in the core area, the percentage of B. genuense ranges from 7.1% to 7.4%, while in the ecotonal areas ranges from 8.7% to 8.9%.

Topographic Model Output
The ensemble model to investigate the topographic driver of B. genuense distribution shows high accuracy with sensitivity higher than 95 and~85 of specificity (Table 3). Six topographic variables with a Pearson correlation coefficient <0.75 were selected [73,74]: elevation, eastness, northness, slope, topographic wetness index (TWI), and solar radiation ( Figure 5). The main topographic drivers are elevation, slope, solar radiation, and TWI in order of importance. Indeed, B. genuense shows a strong preference for a slight slope area with an elevation between 1500 and 2000 m a.s.l. Moreover, the target species prefers humid areas, as highlighted by TWI, which is the most commonly topographic index used to describe a cell's tendency to accumulate water [81]. Moreover, this species shows a slight preference for areas with high solar radiation. On the other hands, the distribution seems not affected by aspect; indeed, both eastness and northness show flat response curves.

Discussion
The tested workflow allowed to map the B. genuense and to identify the main topographic drivers that affect its distribution in the study area. Moreover, the obtained distribution has been used to assess the potential spread of the target species on N2k habitats.
The classification result allowed us to distinguish the species with overall accuracy values > 90%. We agree with other studies on widely spread grasslands and invasive species where the accuracy values ranged from 70 to 90% [82][83][84][85]. RF proved to be the best classification methods among those tested, confirming its ability to map vegetation using various data types [86,87]. Moreover, RF demonstrates its lower sensitivity than other machine learning classifiers to sample quality training and overfitting [88]. Our results suggest that a binary classifier approach (presence/absence) well performs with small training samples, according to recent literature [63,89,90].
In the study area, some plants species, such as Bromopsis erecta and Festuca circummediterranea (with phenology and/or environmental needs similar to the target species), could alter the measure of the spectral signature of Brachypodium genuense, contributing to the spectral mixing [91] and overall classification error. Despite this, the high accuracy values of the RF classification may be due to the decrease in abundance of early and mid-flowering species, as the invasion of B. genuense decreases the heterogeneity of the vegetation from the phenological point of view [92]. Another factor that could explain the classifiers' high performance may be due to the physiological peculiarities (silica-rich and hairy leaves) [93] of the target species. These characteristics allow the overabundance of B. genuense compared to other species as it is not palatable by domestic herbivores (sheep and cattle) in poor stock conditions [94].
It is well-known that the images' spatial resolution must be carefully chosen concerning the spatial scale of the analyzed object [95]. Notwithstanding, the S2 data with a 10 m per pixel spatial resolution could map the target species. Through a careful ground-based survey campaign simultaneously as the acquisition of images by satellite sensors, we have shown how maps deriving from S2 can be a powerful source of information, as confirmed

Discussion
The tested workflow allowed to map the B. genuense and to identify the main topographic drivers that affect its distribution in the study area. Moreover, the obtained distribution has been used to assess the potential spread of the target species on N2k habitats.
The classification result allowed us to distinguish the species with overall accuracy values >90%. We agree with other studies on widely spread grasslands and invasive species where the accuracy values ranged from 70 to 90% [82][83][84][85]. RF proved to be the best classification methods among those tested, confirming its ability to map vegetation using various data types [86,87]. Moreover, RF demonstrates its lower sensitivity than other machine learning classifiers to sample quality training and overfitting [88]. Our results suggest that a binary classifier approach (presence/absence) well performs with small training samples, according to recent literature [63,89,90].
In the study area, some plants species, such as Bromopsis erecta and Festuca circummediterranea (with phenology and/or environmental needs similar to the target species), could alter the measure of the spectral signature of Brachypodium genuense, contributing to the spectral mixing [91] and overall classification error. Despite this, the high accuracy values of the RF classification may be due to the decrease in abundance of early and mid-flowering species, as the invasion of B. genuense decreases the heterogeneity of the vegetation from the phenological point of view [92]. Another factor that could explain the classifiers' high performance may be due to the physiological peculiarities (silica-rich and hairy leaves) [93] of the target species. These characteristics allow the overabundance of B. genuense compared to other species as it is not palatable by domestic herbivores (sheep and cattle) in poor stock conditions [94].
It is well-known that the images' spatial resolution must be carefully chosen concerning the spatial scale of the analyzed object [95]. Notwithstanding, the S2 data with a 10 m per pixel spatial resolution could map the target species. Through a careful ground-based survey campaign simultaneously as the acquisition of images by satellite sensors, we have shown how maps deriving from S2 can be a powerful source of information, as confirmed by Feilhauer et al. [96]. Thus, using the S2 images, the presented workflow can be ap-plied worldwide to map species and communities with distinctive characteristics remotely detectable, as demonstrated by [97,98].
The topographic model suggests that a high-resolution digital elevation model strongly improves plant species distribution knowledge [99]. The results obtained are consistent with the ecology of B. genuense; in fact, a preference for humid areas (high value of TWI) with a high solar radiation value suggests this species' competitiveness [32,100]. Topography can be considered a determinant factor of the distribution patterns in dry grasslands and controls soil moisture and pH. Thus, TWI is an excellent candidate to be a priority driver of local plant diversity patterns in grasslands [101]. These findings suggest the key role of hydrology in species distribution [102,103].
Moreover, the effect of topography on pH distribution could explain the strong effect of TWI in shaping B. genuense. Indeed, this species shows a small pH interval ranged from 6.17 to 7.2 [100]. However, the direct effect of solar radiation on local temperatures [104,105] may also affect the local vegetation patterns. Finally, the elevation and slope variables suggest that this species is limited by high elevation and high slope due to the combination of environmental constraints: decrease of soil depth, low temperature, and winter stress [33].
The overlay analysis on N2k habitats identify the most occupied habitat by B. genuense, and this result can be important to address management policy and conservation strategies for nature conservation. Concerning the relationship between species and N2k habitats, the most affected habitat for B. genuense spread is the semi-natural dry grasslands on calcareous substrates (habitat 6210*), which is also the most threatened habitat in Europe [5,106]. The strong diffusion of B. genuense on habitat 6210* is probably due to the variation in the livestock types (from sheep to cows/horses) and/or the reduction of extensive grazing in the study area. The increase of Brachypodium (B. pinnatum, B. genuense, B. rupestre) distribution cause a reduction in grasslands biodiversity [31,33,92,107,108]. As shown in a recent study [109], when B. genuense and/or B. rupestre becomes dominant (cover > 80%), a substantial reduction in biodiversity was observed. Moreover, as said above, many plant communities in the study area constitute habitats of interest for biodiversity conservation and hosting several plants of high naturalistic value that could be threatened by B. genuense spread.
Other habitats, such as 4060 and 6230*, have shown high B. genuense percentage probably associated with specific topographic characteristics and the limited area occupied in the study area. Indeed, the differences between "realized" and the "expected" cover are less than 20 ha. The 8120 and 6170 habitats have shown a low percentage of B. genuense cover; these results could be explained by the altitudinal distribution of these habitats that go up to 2000 m a.s.l., where the investigated species decrease its suitability.

Conclusions
The proposed workflow can be considered a practical approach to assess the distribution and spread of dominant plant species in grasslands habitats, especially in conservation perspectives. Binary classification coupled with SDMs is an affordable and reliable method for characterizing both distribution and topographic drivers of plant species such as B. genuense. Indeed, the occurrences data obtained from RS classification have proven to be crucial to reduce the sampling bias, and thus, to improve the SDMs calibration. The main results of SDMs have shown as elevation, slope, and TWI are the main distribution drivers of B. genuense.
The overlay analysis between B. genuense coverage and N2k habitats shows that the target species spread across all habitats with a strong preference for 6210*. The monitoring of potentially invasive species is fundamental for managing the N2k habitats, especially in the protected areas where biodiversity conservation is one of the priority goals. In this sense, further studies using high-resolution RS data and multitemporal diachronic analysis could strongly improve the knowledge on the spread of B. genuense. Moreover, the workflow presented in this study may be extended to larger areas and other plant communities for more effective habitats management.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/rs13101904/s1, Table S1: Performance of all tested algorithms by modulating hyperparameters. Figure S2: Results of the area of applicability (AOA) analysis. Table S3: Accuracy statistics. Figure S4: Maps derived from Local Indicators of Spatial Association (LISA) analysis of the classified B. genuense patches. Table S5: Plant associations recognized in the study area. Data Availability Statement: All the RS data are freely available. Presence/absence points used to calibrate the RF classification are available from the corresponding author on reasonable request.