Monitoring plant diversity to support agri-environmental schemes: evaluating statistical models informed by satellite and local factors in Southern European mountain pastoral systems

  Citation: Monteiro, A.T.; Alves, P.; Carvalho-Santos, C.; Lucas, R.; Cunha, M.; Marques da Costa, E.; Fava, F. Monitoring Plant Diversity to Support Agri-Environmental Schemes: Evaluating Statistical Models Informed by Satellite and Local Factors in Southern European Mountain Pastoral Systems. Diversity 2022,14, 8. https://doi.org/ 10.3390/d14010008 Academic Editors: Ingrid Belˇcáková and Michael Wink Received: 9 October 2021 Accepted: 21 December 2021 Published: 24 December 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). diversity Article Monitoring Plant Diversity to Support Agri-Environmental Schemes: Evaluating Statistical Models Informed by Satellite and Local Factors in Southern European Mountain Pastoral Systems Antonio T. Monteiro 1,2,* , Paulo Alves 3, Claudia Carvalho-Santos 4, Richard Lucas 5, Mario Cunha 6,7 , Eduarda Marques da Costa 1and Francesco Fava 8 1 Centro de Estudos Geográficos (CEG), Laboratório Associado TERRA, Instituto de Geografia e Ordenamento do Território (IGOT), Universidade de Lisboa, Rua Edmée Marques, 1600-276 Lisboa, Portugal; [email protected] 2 Istituto di Geoscienze e Georisorse, Consiglio Nazionale delle Ricerche (CNR), Via Moruzzi 2, 56124 Pisa, Italy 3Floradata-Biodiversidade, Ambiente e Recursos Naturais, Lda., Rua do Campinho 32, 4000-151 Porto, Portugal; [email protected] 4Centre of Molecular and Environmental Biology (CBMA), University of Minho, 4710-057 Braga, Portugal; [email protected] 5Department of Geography and Earth Sciences, University of Aberystwyth, Penglais, Aberystwyth SY23 3FL, UK; [email protected] 6Faculty of Sciences, University of Porto, Rua do Campo Alegre, s/n, 4169-007 Porto, Portugal; [email protected] 7Science INESC TEC and Geo-Space Sciences Research Centre (CICGE), Institute for Systems and Computer Engineering Technology, Universidade do Porto, Campus da Faculdade de Engenharia, 4200-465 Porto, Portugal 8 Department of Environmental Science and Policy (ESP), Universitàdegli Studi di Milano, 20133 Milan, Italy; [email protected] *Correspondence: [email protected]; Tel.:+351-210-443-000 Abstract: The spatial monitoring of plant diversity in the endangered species-rich grasslands of European mountain pastoral systems is an important step for fairer and more efficient Agri-Environmental policy schemes supporting conservation. This study assessed the underlying support for a spatially explicit monitoring of plant species richness at parcel level (policy making scale) in Southern European mountain grasslands, with statistical models informed by Sentinel-2 satellite and environmental factors. Twenty-four grassland parcels were surveyed for species richness in the Peneda-Gerês National Park, northern Portugal. Using a multi-model inference approach, three competing hypotheses guided by the species-scaling theoretical framework were established: species–area (P1), species–energy (P2) and species–spectral heterogeneity (P3), each representing a candidate spatial pathway to predict species richness. To evaluate the statistical support of each spatial pathway, generalized linear models were fitted and model selection based on Akaike information criterion (AIC) was conducted. Later, the performance of the most supported spatial pathway(s) was assessed using a leave-one-out cross validation. A model guided by the species–energy hypothesis (P2) was the most parsimonious spatial pathway to monitor plant species richness in mountain grassland parcels (P2, AICc = 137.6, ∆ AIC = 0.0, wi = 0.97). Species–area and species–spectral heterogeneity pathways (P1 and P3) were less statistically supported ( ∆ AICc values in the range 5.7–10.0). The underlying support of the species–energy spatial pathway was based on Sentinel-2 satellite data, namely on the near-infrared (NIR) green ratio in the spring season (NIR/Green spring ) and on its ratio of change between spring and summer (NIR/Green change ). Both predictor variables related negatively to species richness. Grassland parcels with lower values of near-infrared (NIR) green ratio and lower seasonal amplitude presented higher species richness records. The leave-one-out cross validation indicated a moderate performance of the species–energy spatial pathway in predicting species richness in the grassland parcels covered by the dataset (R 2 = 0.44, RMSE = 4.3 species, MAE = 3.5 species). Overall, a species–energy framework based on Sentinel 2 data resulted in a promising spatial pathway for the monitoring of species richness in mountain grassland parcels and Diversity 2022,14, 8. https://doi.org/10.3390/d14010008 https://www.mdpi.com/journal/diversity Diversity 2022,14, 8 2 of 16 for informing decision making on Agri-Environmental policy schemes. The near-infrared (NIR) green ratio and its change in time seems a relevant variable to deliver predictions for plant species richness and further research should be conducted on that. Keywords: biodiversity conservation; species richness; policy monitoring; generalized linear modeling; remote sensing; Sentinel-2 satellite 1. Introduction Agri-Environmental policy schemes are the most powerful tools to support conservation measures in the European agricultural landscapes (e.g., Common Agriculture Policy, CAP) and to ensure that they are supporting biodiversity and the best possible agricultural production is essential [ 1 ]. Among the different strategies that can be followed to pursue this goal, the spatially explicit monitoring of biodiversity represents one of the most relevant. Among the systems that would benefit most from a spatially explicit development of policy are the European mountain pastoral systems and their endangered species-rich grasslands [ 2 ]. Here, grasslands are target systems of Agri-Environmental policy schemes affected by agricultural land use changes and abandonment [ 3 , 4 ] that still need reliable spatial indicators to understand the status of biodiversity and, consequently, the progress made in policy targets [ 5 ]. The achievements of Agri-Environmental policy schemes for biodiversity in grasslands remain only partially known (e.g., in the cases where some field inventories are performed) or just assumed. Allowing the spatial monitoring of plant diversity at policy making levels (e.g., parcel level, unit of payment of conservation schemes to farmers) is then fundamental to efficiently allocate subsidies and delineate local protection schemes for species-rich sites. There has been consistent research on spatially explicit indicators for grassland plant diversity in the last decade, mostly driven by the developments in remote sensing and spatial data science among other factors (e.g., data accessibility, spatial resolution) [6,7] . Remote sensing driven approaches have shown significant relationships between spectral reflectance, spectral vegetation indices and plant diversity metrics [ 8 – 15 ]. For instance, Fava and colleagues [ 8 ], with field spectrometric measurements, found that the best predictive indicators were based on chlorophyll/nitrogen sensible bands in the blue and red spectral regions. Mockel et al. [ 10 ], with airborne hyperspectral data, identified a negative correlation between the reflectance in the near-infrared spectral region and species diversity. This study also suggested that species diversity increased as the above-ground biomass decreased and that spectral response can be a useful indicator for species richness monitoring from space. Interestingly, both studies suggested spectral bands related to plant and structural properties of grasslands as indicators for plant diversity, a relationship theoretically described as species–energy. Other recent efforts with satellite optical data described substantial support for spectral heterogeneity and its temporal variation [ 11 , 13 , 16 ] as species richness indicators. Spectral diversity in the green band of IKONOS multispectral satellite data may provide reasonable estimates for plant species richness [ 16 ]. Fauvel [ 13 ] emphasized the exploitation of spectral heterogeneity in time-series for more accurate prediction of grassland plant diversity, namely time-series of red and infrared bands of Sentinel-2 satellite. The spectral heterogeneity pathway, however, seems to be scale dependent [ 12 , 16 ]. In terms of predictive capacity, studies at higher spectral and spatial resolutions tended to better explain plant diversity data (66 to 72%; [ 8 , 9 ]) than those conducted at lower resolutions [ 13 , 16 ]. Despite the advances, the literature suggests that plant diversity may not be directly measurable by satellite remote sensing [ 17 ]. A widely usable stand-alone model to the European level may not exist [ 12 ]. Hence, calibration of spatially explicit models at local, landscape or regional scales seems necessary, particularly for policy making and for improvement of model accuracy. Addressing the geographical unbalance distribution of knowledge across European mountains would also be valuable (e.g., knowledge in the Diversity 2022,14, 8 3 of 16 Southern European mountains is very reduced). The integration of satellite data with other local spatial environmental factors (e.g., management, topography) is a strategy that may improve the capacity of statistical spatial models to monitor grassland plant diversity [ 13 ]. Remote sensing data, therefore, seems part of a pathway for more accurate spatial monitoring of biodiversity. The links between relevant remote sensing indicators and theoretical species-scaling frameworks (e.g., species–energy, species–heterogeneity) suggests that these frameworks could guide further research in this context. Species–scaling relationships (e.g., species–area [ 18 ], species–energy [ 19 , 20 ], species–disturbance [ 21 ] and species–heterogeneity [ 22 ]) are consolidated attempts to upscale plant diversity. They could serve as exploratory spatial pathways that would avoid an aimless search for indicators by remote sensing community and simultaneously facilitate the dialogue and convergence of spatial biodiversity knowledge existing beyond the remote sensing sector [ 1 ]. In fact, in literature, we can find consolidated in situ knowledge on spatially explicit environmental drivers of grassland plant diversity that has always been analyzed in relation to theoretical species–scaling relationships [ 23 – 25 ]. That could lead to the formulation of more hybrid spatial models. Mechanistic [ 26 ] or statistical modeling approaches [ 8 , 9 , 11 ] can be followed towards a spatially explicit monitoring of grassland plant diversity as a function of satellite and local spatial environmental factors. Statistical modeling represents a bridge that has been commonly used to link field observations and predictor variables [ 17 ]. Within the statistical approach, generalized linear models represent a modeling technique that is frequently considered [ 9 , 27 ]. In addition, in regions where little is known about the links between satellite remote sensing data and biodiversity and, consequently, multiple hypotheses are plausible (e.g., species–area, species–energy), multi-model inference (MMI) can be of help to evaluate the relative strength of evidence in the data for each one [28]. This study addresses the spatially explicit monitoring (hereafter called spatial monitoring) of plant diversity in mountain grasslands in view of more effective Agri-Environmental policy schemes devoted to their protection and management (e.g., CAP). Specifically, it conducts a preliminary assessment with the main goal of identifying a spatial pathway(s) for the monitoring of plant diversity at parcel level in Southern European mountains grasslands. These are understudied pastoral systems in this regard, where EU AgriEnvironmental funding is being used to safeguard grasslands through the payment of subsidies that compensate farmers by maintaining grasslands and adopting sustainable practices. The reference unit (policy making scale) for farmers compensation is the parcel size and sustainable practices are still mainly assessed by questionnaires during the annual payment application. Specifically, we made use of an MMI framework that related species richness collected in semi-natural grassland parcels (n= 24) to a set of three competing hypotheses based on the species–area, species–energy and species–spectral heterogeneity scaling theories. Each hypothesis represented a potential spatial pathway to monitor plant diversity. Specifically, it was assumed that a spatial model guided by the species–energy relationship could be a particularly effective pathway since it was the most responsive in semi-natural grasslands or in conditions similar to our study area [ 8 ]. In addition, spatial models that combine remote sensing and environmental factors may perform better than single remote sensing approaches because they benefit from the merging of well-studied spatial environmental driving factors of mountain grassland plant diversity. Two main questions are therefore addressed: (i) What is the statistical support for a spatial monitoring of plant diversity at parcel level, a policy making scale, in Southern European mountain grasslands? (ii) If support exists, what is the most supported pathway(s) and how do variables composing the best performing pathway(s) relate to plant diversity? The grasslands of the Peneda-Gerês mountain range, a biosphere reserve area and EU Natura 2000 site in the north of Portugal, were chosen to answer these questions. Diversity 2022,14, 8 4 of 16 2. Materials and Methods 2.1. Study Area The study area is located in the Peneda-Gerês mountain range (41 ◦ 42 0 25.2 00 N, 8◦1001.200 W , Figure 1), one of the best representatives of mountain pastoral systems in Portugal. The climate is temperate sub-Mediterranean with mean annual precipitation of 1300 mm and mean annual temperature of 13 ◦ C. Between cold and warm season, there is a large amplitude in mean values ( − 10 to 35 ◦ C). Steep slopes are a main feature (>20 ◦ ) and elevation ranges from 145 to 1253 m above sea level (a.s.l). In the study area, grasslands are mainly permanent hay meadows (habitat 6510) [ 29 ] and represent the main cover type of the traditional man-made farmland mosaics. Grasslands have shaped the landscape and livelihoods during centuries but, due to land abandonment and outmigration since the 1970s, they are now viewed as threatened natural and cultural heritage sites [ 29 ]. To safeguard this pastoral system and its grasslands, different EU policy schemes are acting in the study area, namely on behalf the Common Agriculture Policy (CAP). These policy schemes subsidize pastoral communities to maintain grasslands and sustainable management practices that would also indirectly safeguard biodiversity and ecosystem services. The management of grasslands is extensive, consisting mainly of one cut by the end of June– beginning of July for hay production. After cutting, grasslands can be grazed by livestock. Some grasslands are also irrigated by allowing water to circulate on it throughout the year. This management regime can vary with proximity to settlements, those closely located being more intensively managed. The overall conditions described makes Peneda-Gerês a good representative of Southern European mountain pastoral systems. Diversity 2022, 13, x FOR PEER REVIEW 4 of 17 (ii) If support exists, what is the most supported pathway(s) and how do variables composing the best performing pathway(s) relate to plant diversity? The grasslands of the Peneda−Gerês mountain range, a biosphere reserve area and EU Natura 2000 site in the north of Portugal, were chosen to answer these questions. 2. Materials and Methods 2.1. Study Area The study area is located in the Peneda−Gerês mountain range (41°42′25.2″ N, 8°10′1.2″ W, Figure 1), one of the best representatives of mountain pastoral systems in Portugal. The climate is temperate sub−Mediterranean with mean annual precipitation of 1300 mm and mean annual temperature of 13 °C. Between cold and warm season, there is a large amplitude in mean values (−10 to 35 °C). Steep slopes are a main feature (>20°) and elevation ranges from 145 to 1253 m above sea level (a.s.l). In the study area, grasslands are mainly permanent hay meadows (habitat 6510) [29] and represent the main cover type of the traditional man−made farmland mosaics. Grasslands have shaped the landscape and livelihoods during centuries but, due to land abandonment and outmigration since the 1970s, they are now viewed as threatened natural and cultural heritage sites [29]. To safeguard this pastoral system and its grasslands, different EU policy schemes are acting in the study area, namely on behalf the Common Agriculture Policy (CAP). These policy schemes subsidize pastoral communities to maintain grasslands and sustainable management practices that would also indirectly safeguard biodiversity and ecosystem services. The management of grasslands is extensive, consisting mainly of one cut by the end of June–beginning of July for hay production. After cutting, grasslands can be grazed by livestock. Some grasslands are also irrigated by allowing water to circulate on it throughout the year. This management regime can vary with proximity to settlements, those closely located being more intensively managed. The overall conditions described makes Peneda−Gerês a good representative of Southern European mountain pastoral systems. Figure 1. The study area located in the Portuguese Peneda−Gerês National Park, part of the transboundary biosphere reserve Gerês−Xurez. The upper left panel shows the study area captured by Figure 1. The study area located in the Portuguese Peneda-Gerês National Park, part of the transboundary biosphere reserve Gerês-Xurez. The upper left panel shows the study area captured by the Sentinel-2 satellite sensor in April 2016 (red-green-blue composite) and the grassland parcels surveyed for plant species richness (in red, n= 24). The lower panels present grassland mosaics in the study area and the traditional hay cutting performed by local pastoral communities in the summer season. Diversity 2022,14, 8 5 of 16 2.2. Field Data A vegetation assessment was conducted in 24 grassland parcels to collect data for plant species richness (Figures 1and S1). Surveys were performed between the mid and end of June 2016 before the traditional annual hay cut (see Table S1). In 2016, the traditional hay cut was particularly late due to the extreme precipitation levels registered in May (400 mm, +300%) and June (100 mm, +75%) in the region compared to historical mean values (1971–2000) [ 30 ]. The sampling size was defined based on the operational and budget capacity available for this field assessment. Grassland parcel was the sample unit because it is the unit used by CAP agri-environment scheme to compensate farmers by their sustainable practices in grasslands. The selection of parcels was made to cover the range of edaphic and topographic conditions shaping grasslands (Figure 2). The sampling strategy was opportunistic, and each record contained all vascular plant species present in a 100 m 2 sampling area established inside of each parcel. For the vegetation recording, the scale used was adopted from the “Handbook for Surveillance and Monitoring of European Habitats” [ 31 ], which is more detailed than the traditional Braun-Blanquet scale used in phytosociological studies. The scale has intervals of 5%, with everything less than 5% being treated at intervals of 1% to <1%. The surveyed parcels were mesic meadows (n= 8), meso-hydrophilic meadows (n= 10), meso-xerophytic meadows (n= 3) and hygrophilic meadows (n= 3). Regarding the topographical position, parcels were mainly bottom valley meadows (n= 13) and hillside meadows (n= 11). During the survey, the geolocation of the plot was registered using a GPS and the area of each parcel hand digitized combining field information and the very high resolution satellite imagery archive of Google Earth [32]. Diversity 2022, 13, x FOR PEER REVIEW 5 of 17 the Sentinel−2 satellite sensor in April 2016 (red−green−blue composite) and the grassland parcels surveyed for plant species richness (in red, n = 24). The lower panels present grassland mosaics in the study area and the traditional hay cutting performed by local pastoral communities in the summer season. 2.2. Field Data A vegetation assessment was conducted in 24 grassland parcels to collect data for plant species richness (Figures 1 and S1). Surveys were performed between the mid and end of June 2016 before the traditional annual hay cut (see Table S1). In 2016, the traditional hay cut was particularly late due to the extreme precipitation levels registered in May (400 mm, +300%) and June (100 mm, +75%) in the region compared to historical mean values (1971–2000) [30]. The sampling size was defined based on the operational and budget capacity available for this field assessment. Grassland parcel was the sample unit because it is the unit used by CAP agri−environment scheme to compensate farmers by their sustainable practices in grasslands. The selection of parcels was made to cover the range of edaphic and topographic conditions shaping grasslands (Figure 2). The sampling strategy was opportunistic, and each record contained all vascular plant species present in a 100 m2 sampling area established inside of each parcel. For the vegetation recording, the scale used was adopted from the “Handbook for Surveillance and Monitoring of European Habitats” [31], which is more detailed than the traditional Braun−Blanquet scale used in phytosociological studies. The scale has intervals of 5%, with everything less than 5% being treated at intervals of 1% to <1%. The surveyed parcels were mesic meadows (n = 8), meso−hydrophilic meadows (n = 10), meso−xerophytic meadows (n = 3) and hygrophilic meadows (n = 3). Regarding the topographical position, parcels were mainly bottom valley meadows (n = 13) and hillside meadows (n = 11). During the survey, the geolocation of the plot was registered using a GPS and the area of each parcel hand digitized combining field information and the very high resolution satellite imagery archive of Google Earth [32]. Figure 2. Examples of the different types of grasslands surveyed in the study area located in the Portuguese Peneda−Gerês National Park (photos were taken during the field survey mid−end of Figure 2. Examples of the different types of grasslands surveyed in the study area located in the Portuguese Peneda-Gerês National Park (photos were taken during the field survey mid-end of June 2016). Photo ( a ) shows the mesic grasslands, ( b ) meso-hydrophilic grasslands, ( c ) meso-xerophytic and (d) hygrophilic grasslands. 2.3. Remote Sensing and Spatial Environmental Data Collection Very high-resolution imagery (VHR), Sentinel-2 satellite imagery and a digital elevation model (DEM) were used as data sources to evaluate the set of competing spatial Diversity 2022,14, 8 6 of 16 pathways established to monitor plant diversity in grassland parcels. The VHR imagery was accessed through the imagery archive of Google Earth [ 32 ]. The satellite imagery was obtained from Copernicus Sentinel-2 satellite at Level-1 C in the spring (30 April) and summer season (9 July) of 2016. The images capture the study area before and after the field survey in the grassland parcels. The Sentinel-2 images were composed of 10 spectral bands: four at 10 m spatial resolution (blue b2: 497 nm; green b3: 560 nm; red b4: 664 nm ; NIR b8: 842 nm) and six at 20 m spatial resolution (red-edge b5: 704 nm, b6: 740 nm , b7: 782 nm; NIR narrow b8a: 865 nm; short-wave infrared SWIR b11: 1614 nm, b12: 2202 nm ). The scenes were obtained from the Sentinel Hub and corrected atmospherically with the Sen2Cor algorithm to Bottom-Of-Atmosphere (BOA) reflectance [ 33 ]. For this study, only the green, red and NIR spectral bands were selected. The DEM information was obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) at 30 m spatial resolution. These data were resampled to 10 m resolution using the nearest neighbor approach [ 34 ]. Satellite and DEM data were then stacked, and all data projected in UTM-WGS84—Zone 29 N at 10 m spatial resolution. Operations were performed in ENVI-IDL 5.0 (Exelis Visual Information Solutions, Boulder, CO, USA). 2.4. Competing Spatial Pathways To evaluate the support for a spatial statistical model-based monitoring of plant diversity at parcel level in mountain grasslands (policy making scale), a MMI framework was established. MMI is an approach that allows comparing and ranking a set of competing hypotheses according to their capacity to explain variation in data and to identify model selection uncertainty [ 35 ]. This approach was followed given the lack of previous studies in the study region in this regard. MMI requires a priori definition of competing hypotheses, here guided by the species-scaling theoretical framework and described as spatial pathways. Accordingly, three main spatial pathways were defined: species–area (Pathway 1, P1), species–energy (Pathway 2, P2) and species–spectral heterogeneity (Pathway 3, P3). Each pathway refers to a candidate spatial model composed by predictor variables derived from imagery (VHR, satellite) and/or spatial environmental data defined according to literature. The pathways, rationale and respective predictor variables are given hereafter. 2.4.1. Species–Area (P1) The species–area pathway relies on a primary question: can a spatial model based on grassland parcel area be used to predict plant species richness? Theoretically, it emerges from the empirical scaling relationship between species richness and sampling area (A) [ 36 ] and from the observed increase in species richness with size of grasslands (ha) [ 24 ]. In addition, it considers the fact that parcel size is the unit of compensation of agri-environment schemes to farmers and that policy authorities have georeferenced spatial data on parcel area for each farmer [ 37 ]. The parcel area of each surveyed parcel (n= 24) was used as the predictor variable to evaluate this hypothesis. Parcel area was derived from VHR imagery of Google Earth [ 32 ] by hand digitizing and log transformed (ln) for analysis. Figure 3a illustrates the features of the predictor using the predictor in a sampled parcel. 2.4.2. Species–Energy (P2) The species–energy pathway is based on the theoretical variation of species richness with the availability of environmental energy [ 19 , 20 ] as well as on the developments in remote sensing and ecological research that seem to support species–energy as a pathway for the spatially explicit monitoring of plant diversity. Remote sensing developments show a link between reflectance and plant diversity based on wavebands sensitive to plant pigment content (400–700 nm) and structural properties of grasslands ( 700–1300 nm) , namely a negative correlation with reflectance in the NIR and green regions, and a decline of species diversity with the increase in above-ground biomass [ 8 , 10 ]. The NIR region is a main source of information of grassland structural properties from space [ 38 ]. Ecological Diversity 2022,14, 8 7 of 16 research developments suggest that topographic parameters can enhance the modeling of functional attributes of grasslands [ 39 ]. Specifically, it emphasizes a link between topographic parameters (e.g., slope, aspect) and grassland plant diversity by their effect on light interception and farming management strategies [ 23 , 40 ]. To assess the strength of this hypothesis, four predictor variables were selected: near-infrared green ratio in the spring season (NIR/Green spring , Equation (1)) and its ratio of change between spring and summer (NIR/Green change , Equation (2)), slope and aspect. Each predictor variable represents the mean values inside the sample parcel (n= 24). The Sentinel-2 and DEM data were the source data for the satellite and topographic variables, respectively. Specifically, the NIR/Greenchange represented the monthly ratio of change between spring (April) and summer (July) Sentinel-2 scenes. R842 corresponded to band 8 and R560 to band 3 of Sentinel-2, and time = 2 months. Aspect is a circular variable, and it was converted into continuous data (from − 1 to 1) for analysis by applying the cosine (aspect) and sine (aspect) transformation. Cosine and sine aspect emphasize the north-south and east-west exposures, respectively. Figure 3b illustrates the features of NIR/Green ratio predictor using the predictor in springtime over a sampled parcel. NIR/Greenspring =R842 nm R560 nm (1) NIR/Greenchange =(NIR/Green spring −NIR/Green summer) time(months)(2) Diversity 2022, 13, x FOR PEER REVIEW 7 of 17 Figure 3. (a–c) Illustration of the features of predictors using the sampled parcel geometry, the NIR/Green ratio in the spring and the standard deviation of NIR band in the summer (NIR SD summer) in a sampled parcel. 2.4.2. Species–Energy (P2) The species–energy pathway is based on the theoretical variation of species richness with the availability of environmental energy [19,20] as well as on the developments in remote sensing and ecological research that seem to support species–energy as a pathway for the spatially explicit monitoring of plant diversity. Remote sensing developments show a link between reflectance and plant diversity based on wavebands sensitive to plant pigment content (400–700 nm) and structural properties of grasslands (700–1300 nm), namely a negative correlation with reflectance in the NIR and green regions, and a decline of species diversity with the increase in above−ground biomass [8,10]. The NIR region is a main source of information of grassland structural properties from space [38]. Ecological research developments suggest that topographic parameters can enhance the modeling of functional attributes of grasslands [39]. Specifically, it emphasizes a link between topographic parameters (e.g., slope, aspect) and grassland plant diversity by their effect on light interception and farming management strategies [23,40]. To assess the strength of this hypothesis, four predictor variables were selected: near−infrared green ratio in the spring season (NIR/Greenspring, Equation (1)) and its ratio of change between spring and summer (NIR/Greenchange, Equation (2)), slope and aspect. Each predictor variable represents the mean values inside the sample parcel (n = 24). The Sentinel−2 and DEM data were the source data for the satellite and topographic variables, respectively. Specifically, the NIR/Greenchange represented the monthly ratio of change between spring (April) and summer (July) Sentinel−2 scenes. R842 corresponded to band 8 and R560 to band 3 of Sentinel−2, and time = 2 months. Aspect is a circular variable, and it was converted into continuous data (from −1 to 1) for analysis by applying the cosine (aspect) and sine (aspect) transformation. Cosine and sine aspect emphasize the north−south and east−west exposures, respectively. Figure 3b illustrates the features of NIR/Green ratio predictor using the predictor in springtime over a sampled parcel. NIR/Green=R842 nm R560 nm (1) NIR/Green=(NIR/Green spring − NIR/Green summer) time(months) (2) 2.4.3. Species–Spectral Heterogeneity (P3) The species–spectral heterogeneity pathway is increasingly suggested as valuable approach to achieve spatially explicit representations for biodiversity [41]. It refers to the spatial variability (heterogeneity) in the remotely sensed signal as a proxy of species diversity. It assumes that greater spectral heterogeneity may be a proxy for greater local Figure 3. ( a – c ) Illustration of the features of predictors using the sampled parcel geometry, the NIR/Green ratio in the spring and the standard deviation of NIR band in the summer (NIR SD summer) in a sampled parcel. 2.4.3. Species–Spectral Heterogeneity (P3) The species–spectral heterogeneity pathway is increasingly suggested as valuable approach to achieve spatially explicit representations for biodiversity [ 41 ]. It refers to the spatial variability (heterogeneity) in the remotely sensed signal as a proxy of species diversity. It assumes that greater spectral heterogeneity may be a proxy for greater local species diversity [ 41 ]. The species–spectral heterogeneity relationship [ 9 ] relies on the species heterogeneity theory that assumes an increase in species richness with the number of ecological niches and that species coexistence is facilitated in more heterogeneous environments where different taxa can capitalize on different environmental conditions [ 19 ]. Recent efforts with satellite optical data described substantial support for spectral heterogeneity and its temporal variation [11,13,16]. Spectral diversity in the green band of multispectral satellite data was a consistent predictor for species richness [ 16 ]. Spectral heterogeneity in the red and infrared bands of Sentinel-2 satellite can be relevant explanatory variables for the prediction of grassland biodiversity, particularly when considered from a time-series Diversity 2022,14, 8 8 of 16 perspective [ 11 , 13 ]. Hence, to evaluate the support of this hypothesis, we used the standard deviation (SD) of the reflectance in the NIR and Red bands of Sentinel-2 in the spring and summer scenes in each sample parcel (Equation (3)) as predictor variables. SD is one commonly used diversity metric in the spectral heterogeneity research [ 41 , 42 ]. X ijλ is the reflectance for the spectral band λ and pixel (i, j), X is the mean reflectance in the sample parcel, n is the number of pixels in the sample parcel. NIR SD and Red SD (Equation (3)) for each sample parcel were obtained using the statistics function present in ENVI-IDL and the sample parcel polygons as zonal reference for statistics. Figure 3c illustrates NIR SD using the predictor in the springtime over a sampled parcel. NIRSDspring =s∑(Xijλ−X)2 n−1(3) 2.5. Model Selection and Accuracy Assessment To evaluate the support of the three candidate spatial pathways established for the monitoring of plant richness in mountain grasslands at parcel level (policy making scale), statistical generalized linear models (GLM) were used. GLMs with a Poisson distribution of errors were fitted and model selection based on Akaike information criterion corrected for small sample sizes (AICc) was conducted. The dependent variable was plant species richness registered in the grassland parcels. In order to avoid model overfitting by including more predictor variables than allowed by the sample size (n= 24), the maximum number of predictor variables included in each competing spatial pathway corresponded to one tenth of the number of observations used or two predictors [ 43 ]. Models with AICc difference ( ∆ AICc) < 2.0 were considered models receiving support to explain species richness [ 35 ]. Before fitting GLM models, spearman rank correlations were performed. They evaluated collinearity among predictor variables (no strong correlations were found (r > 0.6, Table S2 ) and helped to identify the final set of predictors included in the species–energy and species–spectral heterogeneity models, where the number of candidate predictors exceed the allowed (n= 4). The final candidate variables included in each of these were those achieving the highest correlation coefficient with species richness (Table S2). After model selection, multicollinearity in the most supported spatial pathway(s) was evaluated by estimating variation inflation factors (VIF) and by checking for overdispersion using the test of Cameron and Trivedi [ 44 ]. To quantify the contribution of each predictor to the most supported model(s), relative importance scores were estimated using the dominance analysis methodology [45]. The performance of the most supported spatial pathway(s) in predicting plant species richness was assessed using the leave-one-out cross-validation approach [ 46 ] since the species richness dataset is small (24 samples) and an independent test dataset would not be truly representative of the overall dataset [ 47 ]. In this approach, the model is fitted and used to predict one hold-out observation at a time. This procedure is repeated until all observations are used for validation. It provides unbiased estimates for in-sample predictive accuracies, though does not allow assessing the model performance for predictions outside the spatial range covered by the dataset [ 46 ]. The metrics cross validated R 2 , root mean squared error (RMSE) and mean absolute error (MAE) were used as performance metrics, and they assess the agreement between observed and predicted absolute values of plant species richness in grassland parcels. All operations were performed in the R statistical environment [ 48 ], using the standard “glm” function for the generalized linear models, MUMIN package for model selection, HMISC for correlation analysis, Domir for the relative importance scores, Caret for validation, AER and CAR packages for dispersion test and VIF factors, respectively. Diversity 2022,14, 8 9 of 16 3. Results The Statistical Support for Spatial Monitoring of Species Richness in Grassland Parcels A total of 104 species were registered in the 24 grassland parcels surveyed in this study. Mean species richness amounted to 19 species, with values spanning from 13 to 33 species (Figure 3and Table S1). The most abundant species were Agrostis capillaris,Holcus lanatus and Arrhenatherum bulbosum. Mean parcel area was 2379 m 2 , with size ranging from 493 to 10,900 m 2 . Summary statistics for the response and predictor variables are given in Figure 4. The results of generalized linear modeling and AICc model selection indicated the spatial pathway guided by the species–energy scaling relationship as the most parsimonious to monitor plant species richness in mountain grassland parcels (P2, AICc = 137.6, ∆AIC = 0.0 , wi = 0.97, Table 1). The spatial pathways species–area (P1, AICc = 147.4, ∆ AIC = 9.7, wi = 0.01 ) and species–spectral heterogeneity (P3, AICc = 145.9, ∆ AIC = 8.28, wi = 0.02) were less statistically supported according to the ∆ AICc values in the range 5.7–9.0. Results highlight a support for a species–energy spatial pathway based on information present in satellite data, namely on the near-infrared green ratio in spring obtained from Sentinel-2 satellite scene (NIR/Green spring ) and on its ratio of change between spring and summer scenes (NIR/Green change ). Both predictor variables related negatively to species richness, indicating that grassland parcels with lower near-infrared green ratio values and amplitude presented higher species richness records (Figure 5). When assessing the performance of the species–energy spatial pathway in predicting plant species richness in the grassland parcels covered by the dataset with the leave-one-out cross-validation technique, results indicated a moderate performance (R 2 = 0.44, RMSE = 4.3 species, MAE = 3.5 species). The VIF values of the most supported model were lower than 2.0, indicating low collinearity among variables (Tables 1and 2). The overdispersion ratio was significantly lower than 1.0 (0.64, z = 2.03, p-value = 0.97), indicating that our Poisson model fits well to the data and data are not overdispersed. In terms of variable importance, the standardized general dominance statistic highlighted that NIR/Green spring has a larger contribution (54.8%) in the model explanation than (NIR/Greenchange, 45.2%, Figure 6). Table 1. Results of the multi-model selection approach with the three different spatial pathways evaluated to monitor plant species richness in the grassland parcels of Peneda-Gerês, Portugal (n= 24) . Model selection was performed using generalized linear models with a Poisson distribution of errors and the Akaike information criteria corrected for small sample sizes (AICc) to identify and rank the pathways. Pathways with ∆ AICc ≤ 2.0 were considered models receiving strong support. aVariable definitions: Species richness is the plant species richness registered in the grassland parcels during the field survey; NIR/Green spring is the near infrared-green simple ratio estimated from the Sentinel-2 satellite in spring; NIR/Green change represents the monthly ratio of change in the near infrared-green simple ratio between spring and summer scene Sentinel-2 satellite scenes. Spring and Summer scenes refers to Sentinel-2 imagery collected on 30 April and 9 July 2016, respectively. Loglik, log-likelihood. AICc, Akaike information criterion value. ∆ AIC are the differences between the AICc values of the best spatial pathway and the other spatial pathways considered. wi are the Akaike weights and indicate the probability that a particular spatial pathway was best among those considered. VIFs are the variation inflation factors for the best spatial model pathway that indicate no multicollinearity affecting the model. Rank Model Structure aIntercept Loglik df AICc ∆AICc wi Spatial Pathway 1Species richness~NIR/Green spring + NIR/Greenchange 3.47 −65.2 3 137.6 0.0 0.97 species–energy (P2) 2Species richness~NIRSDspring + RedSDspring 2.9 −69.4 3 145.9 8.28 0.02 species–spectral heterogeneity (P3) 3Species richness~Parcel area (ln) 3.1 −71.4 2 147.4 9.77 0.01 species–area (P1) Null model 3.1 −71.4 1 145.1 7.42 0.02 VIFs NIR/Greenspring (1.60); NIR/Greenchange (1.57) Diversity 2022,14, 8 16 of 16 38. Rossini, M.; Cogliati, S.; Meroni, M.; Migliavacca, M.; Galvagno, M.; Busetto, L.; Cremonese, E.; Julitta, T.; Siniscalco, C.; Morra di Cella, U.; et al. Remote sensing-based estimation of gross primary production in a subalpine grassland. Biogeosciences 2012 ,9, 2565–2584. [CrossRef] 39. Magiera, A.; Feilhauer, H.; Waldhardt, R.; Wiesmair, M.; Otte, A. Mapping Plant Functional Groups in Subalpine Grassland of the Greater Caucasus. Mt. Res. Dev. 2018,38, 63–72. [CrossRef] 40. Monteiro, A.T.; Fava, F.; Gonçalves, J.; Huete, A.; Gusmeroli, F.; Parolo, G.; Spano, D.; Bocchi, S. Landscape context determinants to plant diversity in the permanent meadows of Southern European Alps. Biodivers. Conserv. 2013,22, 937–958. [CrossRef] 41. Rocchini, D.; Balkenhol, N.; Carter, G.A.; Foody, G.M.; Gillespie, T.W.; He, K.S.; Kark, S.; Levin, N.; Lucas, K.; Luoto, M.; et al. Remotely sensed spectral heterogeneity as a proxy of species diversity: Recent advances and open challenges. Ecol. Inform. 2010 , 5, 318–329. [CrossRef] 42. Imran, H.A.; Gianelle, D.; Scotton, M.; Rocchini, D.; Dalponte, M.; Macolino, S.; Sakowska, K.; Pornaro, C.; Vescovo, L. Potential and Limitations of Grasslands α -Diversity Prediction Using Fine-Scale Hyperspectral Imagery. Remote Sens. 2021 ,13, 2649. [CrossRef] 43. Fieberg, J.; Johnson, D.H. MMI: Multimodel inference or models with management implications? J. Wildl. Manag. 2015 ,79, 708–718. [CrossRef] 44. Cameron, A.C.; Trivedi, P.K. Regression-based tests for overdispersion in the Poisson model. J. Econom. 1990 ,46, 347–364. [CrossRef] 45. Azen, R.; Budescu, D.V. The dominance analysis approach for comparing predictors in multiple regression. Psychol. Methods 2003 , 8, 129–148. [CrossRef] 46. Wenger, S.J.; Olden, J.D. Assessing transferability of ecological models: An underappreciated aspect of statistical validation. Methods Ecol. Evol. 2012,3, 260–267. [CrossRef] 47. Meyer, H.; Lehnert, L.W.; Wang, Y.; Reudenbach, C.; Nauss, T.; Bendix, J. From local spectral measurements to maps of vegetation cover and biomass on the Qinghai-Tibet-Plateau: Do we need hyperspectral information? Int. J. Appl. Earth Obs. Geoinf. 2017 ,55, 21–31. [CrossRef] 48. Team, R. RStudio: Integrated Development for R; PBC: Boston, MA, USA, 2020. 49. Hughes, A. Disturbance and Diversity: An Ecological Chicken and Egg Problem. Nat. Educ. Knowl. 2010,3, 48. 50. Primdahl, J.; Vesterager, J.P.; Finn, J.A.; Vlahos, G.; Kristensen, L.; Vejre, H. Current use of impact models for agri-environment schemes and potential for improvements of policy design and assessment. J. Environ. Manag. 2010 ,91, 1245–1254. [CrossRef] [PubMed] 51. Šumrada, T.; Vreš, B.; ˇ Celik, T.; Šilc, U.; Rac, I.; Udovˇc, A.; Erjavec, E. Are result-based schemes a superior approach to the conservation of High Nature Value grasslands? Evidence from Slovenia. Land Use Policy 2021,111, 105749. [CrossRef] 52. Gillet, F.; Mauchamp, L.; Badot, P.-M.; Mouly, A. Recent changes in mountain grasslands: A vegetation resampling study. Ecol. Evol. 2016,6, 2333–2345. [CrossRef] 53. Rodríguez, A.; de Lamo, X.; Sebastià, M.-T. Interactions between global change components drive plant species richness patterns within communities in mountain grasslands independently of topography. J. Veg. Sci. 2018,29, 1029–1039. [CrossRef] 54. Rapinel, S.; Mony, C.; Lecoq, L.; Clément, B.; Thomas, A.; Hubert-Moy, L. Evaluation of Sentinel-2 time-series for mapping floodplain grassland plant communities. Remote Sens. Environ. 2019,223, 115–129. [CrossRef] 55. Sakowska, K.; MacArthur, A.; Gianelle, D.; Dalponte, M.; Alberti, G.; Gioli, B.; Miglietta, F.; Pitacco, A.; Meggio, F.; Fava, F.; et al. Assessing Across-Scale Optical Diversity and Productivity Relationships in Grasslands of the Italian Alps. Remote Sens. 2019 ,11, 614. [CrossRef] 56. Rossi, C.; Kneubühler, M.; Schütz, M.; Schaepman, M.E.; Haller, R.M.; Risch, A.C. Spatial resolution, spectral metrics and biomass are key aspects in estimating plant species richness from spectral diversity in species-rich grasslands. Remote Sens. Ecol. Conserv. 2021,2056–3485. [CrossRef] 57. Sun, H.; Hu, J.; Wang, J.; Zhou, J.; Lv, L.; Nie, J. RSPD: A Novel Remote Sensing Index of Plant Biodiversity Combining Spectral Variation Hypothesis and Productivity Hypothesis. Remote Sens. 2021,13, 3007. [CrossRef]