Assessing Across-Scale Optical Diversity and Productivity Relationships in Grasslands of the Italian Alps

Abstract

The linearity and scale-dependency of ecosystem biodiversity and productivity relationships (BPRs) have been under intense debate. In a changing climate, monitoring BPRs within and across different ecosystem types is crucial, and novel remote sensing tools such as the Sentinel-2 (S2) may be adopted to retrieve ecosystem diversity information and to investigate optical diversity and productivity patterns. But are the S2 spectral and spatial resolutions suitable to detect relationships between optical diversity and productivity? In this study, we implemented an integrated analysis of spatial patterns of grassland productivity and optical diversity using optical remote sensing and Eddy Covariance data. Across-scale optical diversity and ecosystem productivity patterns were analyzed for different grassland associations with a wide range of productivity. Using airborne optical data to simulate S2, we provided empirical evidence that the best optical proxies of ecosystem productivity were linearly correlated with optical diversity. Correlation analysis at increasing pixel sizes proved an evident scale-dependency of the relationships between optical diversity and productivity. The results indicate the strong potential of S2 for future large-scale assessment of across-ecosystem dynamics at upper levels of observation.

1. Introduction

As extensive changes in biodiversity are occurring globally, ecosystem studies have been focusing on investigating how species composition, distribution, and abundance affect ecosystem functioning. To improve our understanding of such interactions, a deeper insight into processes that regulate vegetation biodiversity and productivity is fundamental [1].

Plant primary production (expressed as above-ground biomass) and productivity (the rate of new biomass generation, estimated through biomass harvesting or expressed as Net Ecosystem Exchange, NEE; [2]) are key metrics of ecosystem functioning. Although positive correlations between biodiversity and productivity (BPRs) were often highlighted in grassland ecosystems [2,3,4,5], other studies showed negative relationships [6,7,8]. A humped-back model was recently proposed for high above-ground biomass rates based on the hypothesis that plant diversity peaks at intermediate productivity, while at low and high productivity biodiversity is lower, as only a few species are able to deal with stress and competition, respectively [8]. Uncertainties about the methodologies beneath BPR studies also exist. The analysis of the factors affecting BPRs is often carried out at small scale using a “within-site” approach to keep environmental conditions constant as much as possible among treatments [9] and is often based on manipulation of species richness. To this regard, Grace et al. [6] highlighted the limits of manipulations experiments and the need for more analyses focused on mature natural ecosystems.

The scale-dependency of BPRs is under intense debate [7,9,10]. Remote Sensing (RS) techniques provide strong potential for across-scale BPRs spatial patterns investigations: Since ground surveys are time-consuming and costly, larger areas can be analyzed using remotely-sensed proxies of diversity and productivity, leading to innovative insights in ecological applications. However, only sporadic experimental RS studies have been carried out to analyze the scale-dependency of optical diversity and productivity relationships. Airborne hyperspectral sensors are able to retrieve and combine detailed information on biochemistry and structure and represent a solid basis for across-scale investigations and for the validation of existing and upcoming satellite observations [2].

Solar radiation interactions with vegetation in different spectral domains provide different information about the observed ecosystems. The visible (VIS) spectral region is dependent on biochemistry and pigments absorptions, while the relatively high reflectance in the near-infrared (NIR) is linked to structure and internal leaf and canopy scattering [11]. By combining structural and biochemical information obtained from light reflected at different wavelengths, optical data can be adopted to map both vegetation productivity and biodiversity [2,12,13].

Canopy structure can be described by structural traits such as Leaf Area Index (LAI), aboveground biomass (AGB), and the proportion of photosynthetic and non-photosynthetic elements. In herbaceous plants, AGB is a spatially and temporally-dynamic trait retrievable from RS observations [14]. A commonly-adopted approach integrating ecosystem productivity and RS is based on the theory of light use efficiency (LUE) [15] which states that a relatively constant relationship exists between productivity (expressed as Gross Primary Production (GPP) or Net Ecosystem Exchange (NEE)) and absorbed radiation at the canopy level:

$$\mathrm{C}\mathrm{uptake}=\mathsf{\epsilon }\times (\mathrm{fAPAR}\times \mathrm{PAR})=\mathsf{\epsilon }\times \mathrm{APAR}$$

(1)

where ε is the photosynthetic light use efficiency expressing the carbon sequestration efficiency per amount of the absorbed solar energy (µmol CO₂·µmol⁻¹ APAR); PAR is the incident photosynthetically active radiation (µmol·m⁻²·s⁻¹); fAPAR is the fraction of PAR absorbed by the vegetation canopy (-), and APAR is the amount of PAR absorbed by the aboveground green biomass (µmol·m⁻²·s⁻¹).

Productivity in herbaceous ecosystems such as grasslands and crops is strongly related to the seasonal variation of aboveground green biomass and is mostly driven by APAR, while ε is often assumed as a constant in the models [16,17]. The fAPAR term was demonstrated to be directly linked to canopy chlorophyll content, and in the Monteith model can be estimated using spectral vegetation indices (SVIs) [18,19], which are mathematical combinations of vegetation surface reflectance at two or more bands [20]. Productivity is usually retrieved using SVIs, including in their formulation a NIR band characterized by a high reflectance due to leaf and canopy scattering, and a VIS band where absorption by the chlorophyll pigments is predominant (e.g., the Normalized Difference Vegetation Index (NDVI), and other Chlorophyll-related indices [20]). The adopted models include an incident PAR component [16,21,22], although some other studies [17,23] demonstrated the suitability of SVIs to estimate productivity without PAR information.

The European Space Agency (ESA) satellite Sentinel-2 (S2) recent mission has opened unprecedented opportunities for BPRs monitoring. The S2 twin satellites are equipped with a Multispectral Instrument (MSI) measuring reflected radiance in 13 spectral bands ranging from VIS and NIR to shortwave infrared (SWIR) with a high spatial resolution. The red-edge domain, which is on the borderline between chlorophyll absorption (in the red wavelengths) and leaf/canopy scattering in the NIR wavelengths [24], is well covered by S2 (B5, B6).

S2 allows the calculation of both “traditional” VIS-NIR SVIs (which are greenness-related, as chlorophyll absorption takes place mostly in the VIS wavelengths up to 700 nm [25]) and NIR indices based on NIR bands ≥740 nm. Such band combinations can provide important information related to canopy light scattering and canopy structure, and therefore canopy structural controls on photosynthesis [26,27].

However, can NIR-based canopy structure observations be an indicator of productivity? To explain this, we cite Ollinger [25], who highlighted the paradox that “the physiological activity of vegetation is often more strongly related to reflectance at wavelengths that are not used in photosynthesis than to those that are”. The nature of the “paradoxical” relationship between the NIR spectral response and ecosystem physiology is related to the complex scattering dynamics of leaf, stem, and canopy-level structural traits, which appear to co-vary with traits related with physiology. Such dynamics have only been partially explored, and more studies are needed to investigate the links between ecosystem physiology, canopy structure, and NIR spectral responses [25,28].

Assuming that leaf biochemistry and canopy structure variability are linked to different strategies to cope with different environmental resource limitations [25], spectral heterogeneity can track relative diversity, which incorporates species richness, biochemical properties, and canopy structure [13,29]. A new objective method to relate hyperspectral RS and ecosystem diversity was proposed [29], and significant correlations were observed between grassland ecosystem diversity (expressed as conventional species diversity indices such as richness and Shannon index) and “optical diversity”, an optically-derived metric of biodiversity defined as the coefficient of variation (CV) in spectral reflectance across space. Optical diversity proved to be a valid approach for the detection of grassland α-diversity (at a local scale, e.g., at the plot level), while larger-scale models have been applied to measure the species turnover among sites (β-diversity, sensu Whittaker, 1972 [30]). We refer to Rocchini et al. [31] for a review of the matter. Following the optical diversity approach, a significant positive linear correlation was observed between optical diversity and productivity, but only for very low productive grasslands (maximum productivity: 136 g·m⁻²; [2]), while studies covering a wide range of productivity are missing.

The aim of this study was to use airborne optical RS data to investigate within-site spatial patterns of ecosystem diversity and functions in highly-diverse mature grassland canopies of the Italian Alps, considering different scenarios at decreasing spatial and spectral resolutions (using resampling techniques) to highlight the potential of S2 to capture ecosystem diversity and productivity patterns.

In this context, a number of research questions are addressed in this paper:

(i)

can we identify optical diversity and ecosystem productivity relationships (for a wide range of productivity) using remotely-sensed proxies?

(ii)

are such relationships scale-dependent?

(iii)

are the S2 data spatial and spectral resolutions suitable to detect such relationships?

Firstly, simulating the S2 bands, we compared the performance of VIS-NIR SVIs and the performance of solely NIR-based (≥740 nm bands) SVIs in estimating the NEE of five different grasslands, and we analyzed the relative contributions of different spectral domains for predicting the productivity of these ecosystems to highlight both structural and biochemical controls.

Secondly, we calculated the coefficient of variation (CV) of reflectance across space (which was demonstrated to be a proxy of ecosystem biodiversity [29]), and thus we explored across-scale optical diversity and productivity patterns for the investigated grasslands using both hyperspectral and S2 simulated multi-spectral bands. Finally, we investigated the S2 data spatial and spectral resolutions suitability for detecting optical diversity and productivity spatial patterns.

2. Materials and Methods

2.1. Study Area

The study area covered a grassland plateau located at Viote del Monte Bondone (Trentino province, Italy; 46°00′N 11°02′E; 1480–1550 m a.s.l.) where grasslands are managed mostly as meadows, with low mineral fertilization and one cut per year, usually in mid-July. The grasslands show a wide range of LAI (varying from 0.4–7 m²·m⁻²), productivity (31–735 g·m⁻² green dry biomass), and species composition ([32]; Figure 1). The vegetation of the area is very heterogeneous, and three main vegetation associations (one of them including two variants) can be found. The Sieverso-Nardetum strictae association (species richness: 71; average aboveground biomass: 236 g·m⁻², [32]) is characterized by short canopies, acid soils and is lightly-managed [33]. The Scorzonero Aristatae-Agrostidetum tenuis (species richness: 88; average aboveground biomass: 384 g·m⁻² [32]) grows on calcareous soils and includes some much more productive variants with species which can be typically found at much lower altitudes (e.g., Arrenatherum elatius, Dactylis glomerata). The plateau, in its Eastern part, consists of small peatland associations of the Caricion fuscae and Caricion davallianae characterized by very high productivity, with values of green aboveground biomass up to 735 g·m⁻² [34].

2.2. Eddy Covariance Measurements

Continuous Eddy Covariance (EC) measurements of CO₂, water vapor, and sensible heat fluxes were performed at five EC towers (Figure 1), one of which is a permanent FLUXNET EC site (IT_MBo, Tower A in Figure 1). The other four temporarily installed EC towers were located within a maximum distance of one kilometer from the IT_MBo tower and were operational three days before the flight campaign, and maintained for two weeks after the flight campaign (8–9 July 2011). All EC systems consisted of the same open-path infrared gas analyzer (Li-7500, Li-COR Inc., Lincoln, NE, USA) and a 3-D ultrasonic anemometer (Gill R3, Gill Instruments Ltd., Lymington, UK or YOUNG Model 81,000 V, R. M. Young Company, Traverse City, MI, USA), mounted at a height from 1.30–2.50 m. A preliminary footprint analysis demonstrated that in most cases for all the five sites, in unstable conditions, 70% of the total CO₂ flux originated within 30 m from the EC tower. Incoming total and diffuse PAR (BF3H, Delta-T Devices Ltd., Cambridge, UK) were measured at the permanent IT_MBo tower (Tower A; Figure 1) and assumed (in cloud-free conditions) to be representative of all the sites. Radiation variables were recorded at one-minute intervals and averaged over 30 minutes. Half-hourly measurements of NEE (µmol CO₂·m⁻²·s⁻¹) and photosynthetically active radiation (PAR, µmol·m⁻²·s⁻¹) were averaged and matched to the time of hyperspectral data flight acquisitions (Figure 1).

2.3. The Hyperspectral Flight Campaign and Imagery Processing

A hyperspectral flight campaign was organized in July 2011. Three flights (Figure 1) were carried out in cloud-free conditions on 8–9 July 2011 at different times of the day using an AisaEAGLE hyperspectral sensor (VIS-NIR, 130 bands; SPECIM, Finland). Hyperspectral images (0.9 m spatial resolution; 400–1000 nm; spectral resolution: 4.3–4.8 nm FWHM) were ortho-rectified, atmospherically corrected, and mosaicked using the ATCOR-4 package of the PARGE software (ReSe applications, Zurich, Switzerland).

According to the preliminary EC footprint analysis, a circular area with a radius of 30 m (region of interest: ROI) was considered, and data were exported and averaged using ENVI 4.8. In order to calculate the SVIs, S2 bands (

Table 1. Specification of the Multi-Spectral Instrument (MSI) of the Sentinel-2 (S2) satellite system. The bands which were simulated for spectral vegetation indices (SVIs) and for the S2 CV calculations are shown in bold.

Band Number	Central Wavelength (nm)	Bandwidth (nm)	Spatial Resolution (m)
B1	443	20	60
B2	490	65	10
B3	560	35	10
B4	665	30	10
B5	705	15	20
B6	740	15	20
B7	783	20	20
B8	842	115	10
B8a	865	20	20
B9	945	20	60
B10	1380	30	60
B11	1610	90	20
B12	2190	180	20

; as in Clevers et al. [35]) were considered, and a simulation was obtained by calculating the average reflectance over the bandwidth of the respective MSI bands following Peng et al. [36]. Both traditional VIS-NIR and NIR-based (NIDI; Normalized Infrared Difference Index) SVIs were investigated (

Table 2. Spectral vegetation indices (SVIs) presented in this study: normalized difference vegetation index (NDVI), green normalized difference vegetation index (NDVIg), red-edge normalized difference vegetation index (NDVIre), MERIS terrestrial chlorophyll index (MTCI), red-edge chlorophyll index (CIre), green chlorophyll index (CIg), and enhanced vegetation index (EVI). The list includes two near-infrared difference indices (NIDI1, NIDI2) using NIR bands ≥740 nm. R refers to reflectance at a specific S2 simulated band (nm).

SVI	Formulation	Reference
NDVI	(R865-R665)/(R865+R665)	[37]
NDVIg	(R865-R560)/(R865+R560)	[38]
NDVIre	(R865-R705)/(R865+R705)	[38]
MTCI	(R865-R705)/(R705+R665)	[39]
CIre	(R865/R705)-1	[40]
CIg	(R865/R561)-1	[40]
NIDI1	(R865-R740)/(R865+R740)	Proposed in this study
NIDI2	(R865-R783)/(R865+R783)	[26]
EVI	2.5*(R865-R665)/(1+R865+6*R665-7.5*R490)	[41]

).

2.4. Models for NEE Estimation

In order to estimate NEE, three different approaches were used: (i) linear regression between NEE and SVIs (model 1); (ii) linear regression between NEE and a product of SVIs and incoming PAR (model 2); and (iii) partial least squares regression (PLSR) using the full set of AisaEAGLE hyperspectral reflectance spectra (model 3). The first two models utilized S2 simulated spectral bands, while in the third approach the full set of hyperspectral AisaEAGLE reflectance spectra (130 bands between 400 and 1000 nm) was used simultaneously to estimate NEE.

The PLSR method, which is commonly adopted to analyze hyperspectral datasets [7], iteratively transforms predictor (reflectance data) and response (NEE) variables into latent vectors and generates band-wise calibration factors used to create a predictive linear model [19,42].

Model 1 was tested for each single flight on a separate basis (as radiation was assumed to be constant at all five towers), while models 2 and 3 (both including PAR) were tested for all the three flights together using the PAR data averaged for the same time period used for EC flux data (Figure 1).

Pearson’s correlation analysis was used to test the significance of the relationships between SVIs and NEE and SVIs*PAR and NEE. The PLSR analyses were carried out by means of the PLS package of the R software environment [43]. The optimal number of components in the PLSR analysis was determined by minimizing the leave-one-out cross-validated root mean square error of predictions (RMSEP_CV). The variable importance projection (VIP) statistic was computed in order to evaluate the relative importance of reflectance at different wavelengths in the PLSR model. VIP values >1 indicated high importance to the PLSR model [19].

The model’s coefficients were obtained by fitting each model against the measured NEE. The main goodness-of-fit statistics (coefficient of determination—R², percentage root mean square error (PRMSE)) were computed to compare the performances of the different models. All the statistical analyses were performed by means of the R soſtware (version 3.0.3, https://www.r-project.org/).

2.5. Ecosystem Function and Diversity Relationships: Spatial Dynamics

To investigate ecosystem productivity and diversity relationships across spatial scales, thirty 60 × 60 m ROIs located on the grassland plateau (Figure 1) were considered (five ROIs corresponding to the EC towers and the other 25 haphazardly chosen ROIs) in grassland homogeneous areas, i.e., where the disturbance due to the presence of rocks, roads, and small bushes was negligible (<1%). The average coefficient of variation in reflectance across space (CV, defined as “optical diversity”) for all the wavelengths between 430–900 nm was calculated (at increasing pixel sizes of 0.9, 2, 5, 10, 20, 30 m, obtained by spatial resampling) as in Wang et al. [2]:

$$\mathrm{CV}=({\mathsf{\Sigma }}_{\mathsf{\lambda }=430}^{900}(\frac{\mathrm{std}({\rho }_{\lambda })}{\mathrm{mean}({\rho }_{\lambda })}))/\mathrm{number}\mathrm{of}\mathrm{bands}$$

(2)

where ${\rho }_{\lambda }$ is the reflectance at wavelength λ, std(${\rho }_{\lambda }$) and mean (${\rho }_{\lambda }$) are the standard deviation and mean value of reflectance at wavelength λ across all pixels in one 60 × 60 m ROI, respectively. CV calculations were performed for both full-spectrum hyperspectral bands (AisaEAGLE) and S2 simulated multispectral bands (bands highlighted in bold, Table 1). CV is representing the variability of the spectral information content of pixels included in a given ROI (Figure 1). Instead of original reflectance values as in Wang et al. [2], the Continuum Removal (CR) reflectance was considered, adopting a spectral transformation technique (using the CR function of the R package “prospectr”; [44]) to enhance individual spectral features reducing noise from the sensor, atmosphere, soil background, topographic variation, and differences in albedo [13].

For all thirty ROIs, the average values of the investigated SVIs calculated based on S2 simulated bands was finally calculated and the across-scale correlations between these SVIs and CVs (calculated from both full-spectrum AisaEAGLE bands and S2 simulated bands) were assessed at increasing pixel sizes. To verify the possible collinearity between the aforementioned variables, a collinearity test was run using the NumXL Microsoft Excel add-in (Spider Financial Corp., IL, USA). Accordingly, maps coupling the optical diversity (CV) and MERIS terrestrial chlorophyll index (MTCI; Table 2) were obtained to visualize optical diversity and productivity patterns for the different investigated ecosystems.

3. Results

3.1. NEE Estimations

The linear regression analysis showed a large range of variation of R² values of the correlations between the investigated SVIs and NEE (

Table 3. Summary of the statistics (R²—coefficient of determination, PRMSE—percentage root mean square error) of the linear regression between measured Net Ecosystem Exchange (NEE) and spectral vegetation indices (SVIs; model 1) for single flights (F1, F2 and F3) and the product of SVIs and photosynthetically active radiation for all flights together (F1F2F3*: SVIs*PAR; model 2).

SVI	Model 1						Model 2
	F1		F2		F3		F1F2F3*
	R2	PRMSE (%)	R2	PRMSE (%)	R2	PRMSE (%)	R2	PRMSE (%)
NDVI	0.26	13.21	0.25	37.68	0.89	6.84	0.77	22.07
NDVIg	0.67	8.85	0.70	23.90	0.72	10.89	0.80	20.61
NDVIre	0.64	9.18	0.41	33.36	0.60	12.91	0.80	20.47
MTCI	0.75	7.68	0.52	30.20	0.79	9.40	0.85	17.89
CIre	0.65	9.02	0.42	33.23	0.61	12.81	0.82	19.34
CIg	0.70	8.37	0.66	25.28	0.73	10.52	0.81	20.17
NIDI1	0.73	7.93	0.52	30.22	0.75	10.19	0.90	14.54
NIDI2	0.00	15.33	0.20	39.05	0.00	20.43	0.75	23.14
EVI	0.84	6.21	0.25	37.65	0.74	10.51	0.80	20.84

). Traditional VIS-NIR SVIs showed highly variable performances when considering flights on a separate basis (model 1: R² up to 0.89). Red-edge based SVIs showed R² values ranging from 0.41–0.79. For model 1, NIDI indices showed extremely variable R² values (up to 0.75).

Considering the model 2, a strong increase of the predictive performance was observed. The R² of the model fed with NIDI1 reached 0.90 and PRMSE 14.54 % (Table 3), and red-edge SVIs R² values were very high (MTCI was the second-best index and reached an R² value of 0.85).

The PLSR (model 3) demonstrated that 95% of the variation in NEE could be explained by algorithms derived from AisaEAGLE spectra (7 latent components; PRMSE = 10.25%), with a heterogeneous statistical contribution of the various wavelengths (Figure 2). The most relevant region of the spectrum was at the very end of the chlorophyll absorption wavelengths (VIP peak at 741.9 nm; VIP value 3.14). The adjacent bands also showed high VIP values: For wavelengths <741.9nm, the VIP values were still higher than 2 for bands 737.3, 732.6, 727.9, and 723.3 nm. Also, considering the NIR shoulder bands >741.9 nm (towards the NIR-shoulder region), band 746.6 and 751 showed VIP values higher than 2. VIP values were higher than 1 up to 800 nm, indicating a relevant contribution of the NIR spectral domain.

3.2. Ecosystem Optical Diversity and Productivity

Ecosystem optical diversity showed linear positive correlations with the best NEE predictors (Figure 3), indicating evident optical diversity and productivity relationships. The collinearity test showed that there was no collinearity between SVIs and CVs as at the Variance Inflation Factor (VIF) values range was included between 1 and 2.61. The highest correlations between SVIs and CVs (CV calculated from full spectrum) were observed for MTCI (R² = 0.61) and NIDI1 (R² = 0.51), which were, in turn, the SVIs which were best correlated with NEE. Quite surprisingly, the widely-used SVIs NDVI and EVI, although they were good proxies of NEE (Table 3), showed very low correlations with CV. Across-scale observations indicated a clear scale-dependency of optical diversity and ecosystem productivity relationships, with an evident across-scale hump-shaped pattern of R² values for most SVIs. The best results were mostly obtained with 10 meters pixels (Figure 3), while the 20 meters spatial resolution (corresponding to the S2 resolution of the bands needed to calculate all the investigated SVIs) was still able to provide important information on ecosystem diversity and productivity relationships (MTCI R² = 0.54). Conversely, when the 30-meter spatial resolution was considered, a noticeable decrease in the coefficient of determination was observed (MTCI R² = 0.40).

When optical diversity (CV) was calculated from S2 simulated bands (Figure 3), a similar across-scale pattern was observed, with only slightly lower R² values compared to the ones obtained when CV was calculated from full-range hyperspectral data.

The MTCI and CV maps including the average optical diversity and the MTCI productivity proxy values within the 30 60 × 60 m ROIs enabled the visualization of the ecosystem diversity-productivity patterns for the investigated grassland ecosystems (Figure 4).

4. Discussion

4.1. Chlorophyll and Structural Controls on Ecosystem Function

Analyzing the relative contributions of different spectral domains for predicting productivity of different grassland ecosystems, both structural and biochemical controls were highlighted (Figure 2). Numerous studies demonstrated how the red-edge spectral domain [45] can be used for monitoring chlorophyll-related canopy traits (i.e., fAPAR and LAI) and ecosystem productivity of herbaceous ecosystems [16,17,18,19,24,47,48].

Band combinations investigated by these authors included at least one band <740 nm, in the region of the spectrum which is sensitive to chlorophyll absorption.

Indices calculated with NIR band combinations are expected to be influenced mostly by canopy structure [25,28,48]. As canopy chlorophyll is usually calculated as the product of foliar chlorophyll and LAI [36], the correlations between canopy chlorophyll content and NIR-based indices may be indirect and mostly driven by factors related to canopy structure and architecture. As demonstrated by Ollinger [25], variables such as cellular leaf anatomy, leaf clumping, and leaf angle distribution (which are contributing to the whole canopy architecture and, in turn, to scattering and reflectance in the NIR region) co-vary with traits strongly related to photosynthesis. The underlying mechanisms of such co-variations have only been partially explored and need more research.

Within model 1, the R² values between some SVIs and NEE showed variations among the three different flights (e.g., NDVI performance was much better in flight 3 in respect to the other two flights). Such variations very likely result from uncertainties which are typically associated with both EC flux [49,50,51] and spectral measurements [52,53], and from the constrained number of tower observations per flight available in the study. However, distinctive results were shown by both model 1 and model 2. The good performance of NIR-based indices is pointing out a substantial canopy structural control on photosynthesis. The correlations between NIR-based indices and productivity were often stronger than the ones observed with traditional VIS-NIR SVIs, as shown by Vescovo et al. [26] and Matthes et al. [27]. But how strong are structural controls on productivity? Can we expect relevant structural controls even in ecosystems where leaf chlorophyll content is strongly varying? According to the results of Sakowska et al. [54], who compared the photosynthetic activity of a soybean chlorophyll deficient mutant (with approximately 80% less chlorophyll than green varieties) and a wildtype, the leaf and canopy structure still appeared to be the key drivers for light absorption and photosynthesis dynamics.

4.2. Optical Diversity and Productivity Across-Scale Observations

In contrast with what was observed by Fraser et al. [8] at the global scale, our within-site observations on optical diversity and ecosystem productivity seemed to indicate positive linear BPRs for grasslands including a wide biomass range. More RS-based studies following the optical diversity approach (investigating the link between biodiversity indices and optical diversity, and combining multiple within-site observations, at upper scales) are recommended to clarify the possible impact of different sampling approaches and observation scales on the obtained results.

The optical diversity-productivity relationships were shown to be strongly scale-dependent as observed using a modeling approach for detecting BPRs by Grace et al. [1]. In our study, SVIs which were best correlated with optical diversity (MTCI and NIDI1; model 2 in Table 3) were, in turn, also the best proxies of productivity.

Spectral diversity incorporates not only ecosystem characteristics detectable at very small scale, such as species richness and leaf biochemical properties, but also canopy structure, which is expected to be detectable at relatively larger scales. Accordingly, Wang et al. [29] suggested that the relevant patch size (defined as a relatively homogenous spatial unit) for canopy structure is expected to be larger than the one for leaf biochemical traits. Our results on NIR-based indices confirm the importance of structural controls on diversity and productivity relationships, but on another hand suggest comparable biochemistry and structure optimal patch sizes, as the optical diversity vs. productivity R² values at increasing pixel sizes of both VIS-NIR and NIR-based indices (Figure 3) showed a similar hump-shaped trend.

NDVI and EVI (Figure 3) showed, in most cases, extremely low coefficient of determination values. These results are in strong disagreement with the observations of Wang et al. [2] for NDVI, and this is probably due to the strong saturation effect [26,34], which is detectable for medium-high values of productivity.

4.3. Interpretation of the Diversity Measures Detectable with Optical Sampling

Across-scale studies are vital to verify that the spatial resolution of remote observations are matching with the grain size of field data, considering that the “optimal” pixel size is varying among biomes and communities [29] and that actual spatial heterogeneity information may be hidden at the sub-pixel level [12]. Our study, showing that the optimal pixel size for detecting optical diversity and productivity relationships was mostly 10 meters, points out the potential of adopting coarser spatial data to analyze diversity and productivity spatial patterns in grassland ecosystems and is opening the way to further research on the nature of the observed diversity. The observed across-scale ecosystem diversity and productivity linear relationships may be linked to different measures of diversity. While the observed relationships at high spatial resolution (0.9 m) may be associated to grassland α-diversity at the plot level (as the pixel size is more comparable with the plant size), for larger-scale observations the diversity-productivity relationships may be related to other measures of diversity, such as β-diversity [29]. In the investigated grasslands ecosystems, β-diversity may be related to small-scale spatial variability of ecological factors (e.g., soil properties, slope) and varying disturbance effects (e.g., land use intensity).

4.4. The Potential of Sentinel-2 Optical Data to Analyze Ecosystem BPRs at the Global Scale

The most striking result of our study is that, when the optical diversity (CV) was calculated from S2 simulated data (thus both spectral and spatial resolutions corresponded to the S2 sensor), across-scale patterns of diversity and productivity still could be observed, with only slightly lower R² values compared to the CV calculated from full-range hyperspectral data (Figure 3).

Such results are opening exciting perspectives on the possible use of S2 to monitor grassland BPRs spatial patterns starting from within-site observations (focusing on the links between spectra and biodiversity indices) and scaling up to the landscape, regional, and global levels of observation. S2 is expected to help ecologists to extensively explore BPRs patterns linked to higher levels of diversity (such as β-diversity) for grassland ecosystem types where single plant size is small. Using remote sensing imagery, selected optical-based proxies of plant traits are expected to significantly improve the combined analysis of taxonomic and functional diversity [31]. Also, further studies are needed to assess the suitability of coarser resolution data to monitor diversity and productivity patterns in other ecosystem and vegetation types (e.g., forests, shrublands) where the higher single plant size is expected to facilitate the detection of lower-level diversity.

5. Conclusions

In this study, we adopted an optical-based approach to estimate productivity and optical diversity-productivity relationships. Our results demonstrated both structural and biochemical controls on grassland productivity. We provided empirical evidence that the best optical proxies of ecosystem productivity were also the best correlated with optical diversity. Our study assessed the linearity of optical diversity and productivity relationships for a wide range of productivity and highlighted their strong scale-dependency.

Using the optical diversity approach, we proved that S2 multispectral data are able to provide unprecedented datasets to extensively explore optical diversity and productivity patterns, providing a first basis for using satellite observations for assessing BPRs within and across different ecosystem types, scaling up to the landscape, regional, and global levels for novel developments in ecology studies.