Beta-Diversity Modeling and Mapping with LiDAR and Multispectral Sensors in a Semi-Evergreen Tropical Forest

Abstract

Tree beta-diversity denotes the variation in species composition at stand level, it is a key indicator of forest degradation, and is conjointly required with alpha-diversity for management decision making but has seldom been considered. Our aim was to map it in a continuous way with remote sensing technologies over a tropical landscape with different disturbance histories. We extracted a floristic gradient of dissimilarity through a non-metric multidimensional scaling ordination based on the ecological importance value of each species, which showed sensitivity to different land use history through significant differences in the gradient scores between the disturbances. After finding strong correlations between the floristic gradient and the rapidEye multispectral textures and LiDAR-derived variables, it was linearly regressed against them; variable selection was performed by fitting mixed-effect models. The redEdge band mean, the Canopy Height Model, and the infrared band variance explained 68% of its spatial variability, each coefficient with a relative importance of 49%, 32.5%, and 18.5% respectively. Our results confirmed the synergic use of LiDAR and multispectral sensors to map tree beta-diversity at stand level. This approach can be used, combined with ground data, to detect effects (either negative or positive) of management practices or natural disturbances on tree species composition.

1. Introduction

Mexico is the fourth country with the greatest biological wealth, in which approximately between 9% and 12% of the world’s discovered species live, despite representing approximately only 1.4% of the global territory [1]. Much of this biological wealth is found in its 65 million ha of temperate and tropical forests [2], and approximately 6.2 million ha are formally managed for timber production [3].

The management of forests for the production of any goods, particularly timber, undoubtedly impacts their species composition and, consequently, their biodiversity. Especially in the tropics, many efforts have been done to incentive a more whole use of biodiversity, but it’s yet to accomplish, considering the existence, for example, of more than 100 timber-yielding species recognized in the Mexican Yucatan peninsula [4], but the commercial interest focus on only six [5]. Even when the national strategy for biodiveristy conservation [6], through forest managment approaches at habitat scale [7,8], suggests keeping a mosaic representing the natural or historic variation with its corresponding tree-species mixing [9], it is hard to monitor such a task because of the relative complexity of describing the dissimilarity in tree species composition between the vegetation patches of the managed landscape.

Tree community composition is a key indicator of forest degradation [10,11], because it is one of the last recovering variables after a disturbance [12,13], and consequently it is an excelent measure of diversity loss. Crown-level assessment of species composition has been increasingly successful for boreal and temperate forests [14,15,16], but its applicability remains challenging for the most complex tropical forests. Alternatively, the stand-level approach provides parameters such as the position of a specific plot on a floristic gradient of tree species dissimilarities ( $\beta$ -diversity). It denotes the variation in species composition among sites, giving the possibility of identifying areas with high species turnover, and planing areas of conservation not exclusively based on $\alpha$ -diversity, but considering simultaneously the minimum overlap in species composition, and maximizing the number of species preserved in the whole region [17,18,19]). While $\alpha$ -diversity is helpful for identifying areas with high local species richness and is amply documented (even for the study area [20]), $\beta$ -diversity has seldom been studied and it is as important as $\alpha$ -diversity.

The tree $\beta$ -diversity of a region of interest can be summarized in a floristic gradient in which the more similar the calculated value, the more akin the tree species composition is. It is usually generated by an ordination method like the detrended correspondence analysis [19,21], or more often with ordination methods based on Bray–Curtis distances as the non-metric multidimensional scaling (nMDS) [22,23,24,25,26,27,28,29], or the isomap [30]. It helps to illustrate the mechanistic basis for the interpretation that is not apparent from high dimensional species composition data, and it is also advantageous because it preserve the general pattern in the original composition matrix, avoiding problems related to the large number of zero observations [26,31].

Compositional ordination has been widely used along with remote sensing methods because, contrary to field mapping, it proffers the potential to extrapolate the floristic gradient over space, from the plots to the surrounding area in a continuous way [22,28,32]. In some studies, it is just an intermediated step to ulterior purposes like analyze topics on vegetation ecology [33,34], predict the pertinence probability of each pixel to certain habitat types [27] or classify intact and degraded forest with multispectral imagery [10]. Some other studies have focused specifically on modeling compositional ordination, in order to propose a predictive map of the floristic gradient, which can directly give us an intuitive access to dissimilarity patterns, very difficult to discern from categorical maps [11]. Several pieces of research have been done in grasslands [21,27,29,30], and temperate [22,24,26] or tropical forested landscapes [19,23,25].

Many remote sensing covariates are recognized as valuable for biodiversity mapping [8,35,36,37,38]. Among them, multispectral textures have proven to be a tool of great worth [20,25,39], because they quantify the variability in reflectance values among pixels and this spectral heterogeneity is expected to be related to the heterogeneity of resource distribution or ecological variability, so that an increase in complexity allows more variability of habitat niches, as well as a higher number of species to coexist (spectral variation hypothesis) [18,40,41].

On the other hand, active sensors are widely used to derive forest structural features for many purposes [42] as biomass estimations [43], successional stages characterization [44], or biodiversity mapping [45]. LiDAR metrics like canopy height are also directly related to the complexity of vegetation structure. Having low values of canopy height will indicate a simple vegetation structure (generally associated with low $\alpha$ -diversity values and certain characteristic species mixing), whereas high values will indicate a more complex structure generally associated with higher species $\alpha$ -diversity and a different tree species mixing [46,47]. Active sensors have been used alone and complementing passive instruments, frequently outperforming either when used on its own [16,22].

We hypothesized that $\beta$ -diversity, summarized in a floristic gradient through an nMDS ordination, can be efficiently modeled with variables derived from LiDAR, and multispectral textures and that the common disturbance history areas will have less inner dissimilarity. It was supposed because as a result of the disturbance, a variation in tree species composition occurs and forest general attributes change [25], being some of the changing traits, like leaf structure, foliar chemistry or stand structure [21] trackable from remote sensing tools [23,48]. Our objectives were to map $\beta$ -diversity as a floristic gradient in a continuous way with remote sensing technologies over a tropical landscape, and to relate it to different disturbance histories.

2. Materials and Methods

2.1. Study Site

The site is located in the Ejido Felipe Carrillo Puerto (18 $°$ 53 ${}^{\prime }$ N–88 $°$ 14 ${}^{\prime }$ W) in central Quintana Roo, México [49]. The topography is karstic and varies from flat to undulating at less than 20 m of elevation above sea level [50]. It has a tropical warm climate with summer rain, its mean annual temperature is 26 °C, and has an annual average rainfall of 1200 mm [51]. Dominant species are Gymnanthes lucida, Manikara zapota, Bursera simaruba and Psidia psipula, with two or three canopy layers and a maximum canopy stature of 25 m.

The 900 ha study area consists of a semi-evergreen tropical forest mosaic after slash-and-burn agriculture (Ac: 8–23 years ago), selective logging (Fs: 43–53 years ago), and selective logging and forest fire (Fc: 21–28 years ago). We have sum up some features of stand structure (basal area, height, biomass and tree density) and $\alpha$ -diversity metrics (species richness and Shannon-Wiener index) in Table 1 in order to give a general context around the dissimilarities in species composition. The whole pipeline for data analysis is depicted in Figure 1.

2.2. Field Data Collection and Processing

Twenty eight sample units, located under a systematic sampling design, with a few modifications to represent all the vegetation patches present in the area, were surveyed. Each unit consisted of four 400 m ${}^2$ circular subunits, on an inverted Y array, similarly to Mexico’s national forest inventory and soils design [52]. Ground data were collected between March and August 2013. For all woody plants that had a DBH >10 cm, height, and diameter at breast height (DBH) were measured, and the taxonomic identity was determined. Botanical nomenclature followed The Plant List [53].

The ecological importance of the surveyed species was assessed through the percentage Importance Value index (IVI(%)) [54] at 400 m ${}^2$ level. Compositional ordination of the floristic data collected was performed using the non-metric multidimensional scaling (nMDS) technique computed over the log-transformed IVI(%) for the 112 subunits using the “metaMDS” function of the “vegan” package [55] in the R statistical software [56]. The log-transformation was applied to reduce the importance of very high values [31].

This ordination technique is widely used for synthesizing the typical matrix of a large number of species by plots, into a reduced dimensional space (two metric variables in this case), that could be regressed against remotely sensing data. nMDS is a non-parametric technique based on Bray–Curtis distances that preserves the ranked dissimilarity in the compositional space [57,58], with precedent in both ecology and remote sensing [11,22,25,59]. The ordination results were evaluated using the stress function [55], which characterizes the disagreement between the rank-order of distances in the original multidimensional data, relative to the distances given by the ordination in the bidimensional space. The lower the stress, the closer the ordination to the original distribution is [58].

The normality of the floristic gradient was tested with the Shapiro–Wilk’s goodness-of-fit test using the function “saphiro.test”, and the differences in the ordination axes scores between the three vegetation patches were tested using a Kruskal-Wallis and a Wilcoxon test for paired comparisons using the functions “kruskal.test” and “pairwise.wilcos.test” all implemented in the R base package “stats” [56].

2.3. Remotely Sensed Data

2.3.1. Multispectral Imagery

We acquired the 1,649,010 scene of RapidEye satellite imagery in 3A processing level. It has a spatial resolution of 5 m, and a spectral resolution of five bands: blue (440–510 nm), green (520–590 nm), red (630–685 nm), redEdge (690–730 nm) and infrared (760–850 nm). The presence of the redEdge channel is a unique feature of RapidEye imagery that was found to be dependent on chlorophyll concentration, and leaf and canopy structure [60,61]. The 16 bit raw digital numbers were transformed to top of atmosphere reflectance after a radiance transformation [62]. With this procedure illumination variations were corrected [63] assuming a uniform atmosphere.

We also calculated the RedEdge Enhanced Vegetation Index (reEVI, Equation (1)) and the RedEdge Normalized Difference Vegetation Index (reNDVI, Equation (2)).

(1) $reEVI=G\frac{\mathrm{redEdge}-\mathrm{red}}{\mathrm{redEdge}+C_1\times \mathrm{red}-C_2\times \mathrm{blue}+L},$

(2) $reNDVI=\frac{\mathrm{redEdge}-\mathrm{red}}{\mathrm{redEdge}+\mathrm{redEdge}}.$

where G is the gain factor (2.5), $C_1$ and $C_2$ are the aerosol resistance weights (6 and 7.5 respectively) and L is the canopy background adjustment factor (1). In this index the blue band is used to stabilize aerosol influences in the red band, being less prone to saturate in high biomass scenarios. It is designed to extract canopy greenness, independent of the underlying soil and atmospheric spray variations [64,65].

The texture, or spatial arrangement of the reflectance values within the image [66] of the five original bands and the indexes, was used as a proxy for spectral heterogeneity, which has long been recognized as a key landscape characteristic with strong relevance for biodiversity [39,67,68]. Several surface metrics were extracted from each band. The Mean metric is considered an occurrence metric because it only quantify the tone and variation of the reflectance values within the window. While variance, homogeneity, contrast, dissimilarity, entropy, and angular second moment, are known as co-occurrence or Haralick’s texture metrics, because those consider all the spatial relations between groups of adjacent pixels.

The co-occurrence metrics derive from a gray-level co-occurrence matrix (GLCM), which is a tabulation of how often different combinations of gray levels occur at a specified distance and orientation in an image object [66]. It was computed in 64 gray levels, with an offset of three pixels (15 × 15 m window size considering the rapidEye spatial resolution) in order to approach the plot size without exceeding it. Four matrices for each one were calculated (directions: 0 $°$ , 45 $°$ , 90 $°$ , 135 $°$ ) and then averaged to obtain a single textural directionless metrics. Using the plots coordinates and size, we extracted the mean values of each layer by plot. The data was processed using the R statistical software and the “raster” [69] and “GLCM” [70] packages.

2.3.2. LiDAR

LiDAR data were collected on January of 2013, by a private contractor, CartoData, operating a Cessna T202 aircraft with an airborne laser scanner, RIEGL-QV-480 LiDAR, equipped with a NovAtel GPS/IMU. The system was operated at an average height of 396.2 m above ground level, a 30 $°$ field of view and a pulse repetition frequency of 200 kHz, for which the aircraft maintained a ground speed between 80 and 90 kph. Flights had an approximate overlap of 50% between adjacent flight lines, averaged more than five pulses per square meter and included up to five returns for each pulse.

After ground normalization, the clouds of points were clipped to 400 m ${}^2$ and a set of 95 LiDAR metrics were extracted from them. To reduce the noise within the near-ground cloud of returns we use a 1 m threshold as a minimum height above ground. We also calculated metrics by strata (1–10 m and above 10 m) in order to found metrics reliable to depict the canopy of vegetation patches at different heights (lower at slash and burn agriculture antecedent, and higher at extensive logging antecedent). Some of those are metrics based on height statistics, as mean, mode, maximum and minimum elevation, metrics of dispersion of return height (as variance, standard deviation), and statistical height distribution (as percentiles or L-moments). The others describe the density of the return cloud and are mainly counts of returns (number of first returns above 10 m, number of returns above mean), and ratios (Canopy Relief Ratio, All returns above 10 over total first returns). We also considered the metrics extracted from the smothered Canopy Height Model (CHM), the Digital Terrain Model (DTM) and cover metrics at different thresholds: 5 and 7 m for assessing the understory cover, and 10 and 15 m for the canopy cover. LiDAR data were processed using FUSION software [71], and the CHM were smoothed with a gaussian filter with a window size of 3 × 3, in order to eliminate spurious local maxima caused by tree branches with the “rLIDAR” [72] package in R.

2.4. Correlation Test

Spearman’s coefficients were calculated to evaluate the magnitude and significance of the correlations between the remotely sensed covariates and the nMDS axis-1 scores ( $\alpha =0.05$ ). We have chosen this association measure considering that our response variable was calculated with a rank based method (nMDS), and likewise, it was used to mitigate the outliers effect, and facilitate the application to non-normal data distributions [73].

2.5. Model Fitting

LiDAR and RapidEye metrics altogether were used as the predictor variables for estimating the spatial variability of the nMDS axis-1 scores, our proxy to $\beta$ -diversity, meaning the dissimilarity scaling in tree species composition. In order to map $\beta$ -diversity for the entire area of study, we fitted a linear regression model with the nMDS axis-1 scores of surveyed plots as the response variable, and the remote sensing metrics as covariates. Note that there is a regression coefficient associated with each metric, that once that is estimated can be used to predict $\beta$ -diversity.

The variable selection was based on a preliminary ranking of the explanatory variables through a mixed model (Equation (3)).

(3) $\underline{y}=\mu \underline{1}+{\underline{X}}_j{\beta }_j+\underline{u}+\underline{e},$

where $\underline{y}$ is the response vector nMDS-1, $\mu$ is an intercept, ${\underline{X}}_j$ is the jth column of X matrix wich contains all covariates, ${\beta }_j$ is a regression coefficient for the jth covariate, and $\underline{u}$ is a random effect, $\underline{u}\sim NM(0,{\sigma }_g^{2}G)$ , where G is $X{X}^{\prime }$ .

The linear mixed model includes as fixed effect the covariate to be tested and a random effect that takes into account that there exists other covariates that needs to be taken into account to avoid correlation issues, a strategy similar to this is used in genome wide association studies in genetics [74]. The model was fitted in R using the “BGLR” package [75]. We tested $H_0:{\beta }_j=0$ vs. $H_1:{\beta }_j\ne 0$ , $j=1,\dots ,p$ (number of covariates), and asses iteratively the significance of each variable without multicollinearity damage, having a number of covariables that exceeds sample size (n).

Then, a backward stepwise multiple linear regression of the top ten variables was performed to select the best explanatory variables. Linear dependencies were explored by computing the variables’ variance inflation factors (VIF). The “lmg” function in the “Relaimpo” package [76] was used to assess the relative importance of LiDAR and RapidEye derived variables. This function calculate the proportion of contribution of each variable makes to the $R^2$ , considering both, its direct effect and its effect when combined with the other variables in the regression equation. Metrics were normalized to sum to 100% of the explained variance, and 95%-bootstrap confidence intervals of the relative weights were calculated using the “boot.relimp” function [76] with 1000 bootstrap replications.

Model predictive ability was assessed via a ten-fold cross-validation, with 100 repetitions. Random samples of the full dataset were iteratively split into ten training and validation subsets to test the fit of predicted versus observed values. The model was qualified over the Pearson’s correlation coefficient between observed and predicted values (r), coefficient of determination ( $R^2$ ), the root mean squared error (RMSE), and the Mean Absolute Error (MAE). Normality of residuals were checked graphically and with the Kolmogorov-Smirnov test, and its randomness were assessed through a residual vs. fitted values plot. The final model was applied in a cell-by-cell basis to the 3 km ${}^2$ study area, using the function “predict” available in the “raster” package [69]. For modeling it, the LiDAR selected variables were gridded from the cloud points at the same extent and spatial resolution of the rapidEye imagery layers. All analyses were performed in the R statistical analysis environment [56].

3. Results

3.1. Species Compositional Ordination

The nMDS scores had not a normal distribution (Shapiro-Wilk test, p > 0.0001). There were significant differences in the nMDS axis-1 scores among patches of vegetation with different disturbance antecedents (Kruskal-Wallis test, p > 0.0001, Wilcoxon test, p < 0.0001), but the axis-2 showed no evidence for significant differences among patches (Wilcoxon test, p > 0.05), so there is not a clear ecological effect to report, and we considered the nMDS axis-1 alone as the representation of the floristic gradient, and henceforth we will be referring to it when writing $\beta$ -diversity.

The plots located in the vegetation patch with slash and burn agriculture antecedent have the highest scores of nMDS, whereas the lowest values were observed for the vegetation patch with the antecedent of extensive logging without forest fire disturbance (Figure 2). There is an overlapped area between plots with extensive logging, with and without forest fire disturbance, coherent with the shared land use antecedent and the IVI% of the most dominant species (Gymnanthes lucida). Even in this case, there is a significative difference between means scores of the gradient.

3.2. Correlation Analysis

Spearman coefficients of the textural metrics achieved higher correlations than the coefficients of the LiDAR metrics (Table 2). Mean and variance of the red, green, reNDVI and redEdge channels were all significantly positively correlated ( $r_s$ > 0.70, p < 0.0001) with $\beta$ -diversity values. The higher correlation of the LiDAR variables was found with the Canopy Relief Ratio ( $r_s$ = −0.616, p < 0.0001). Then appears one metric of statistical height distribution (Height elevation L3 moment; $r_s$ = 0.537, p < 0.0001), one of height dispersion (Height skewness; $r_s$ = 0.537, p < 0.0001), and a metric of density (All returns above mean over total of first returns × 100; $r_s$ = 0.537, p < 0.0001). The cover metric with 7m threshold also achieved negative correlation ( $r_s$ = −0.51, p < 0.0001), closely followed by a metric derived from the CHM ( $r_s$ = −0.517, p < 0.0001). The other metrics of height distribution such as the n-percentiles also showed statistically significant correlations ( $r_s$ > 0.428, p < 0.0001).

3.3. Predictive Model

The final model explained 68.5% (RMSE = 0.20, MAE = 0.14) of the spatial variability of the $\beta$ -diversity (Equation (4)).

(4) $\widehat{nMDS}=-0.7276-0.02068\times Canopyheightmodel+14.0307\times b4mean-0.0009207\times b5variance$

The three most highly significant predictors (p < 0.01) were the redEdge band means (RapidEye derived), the CHM (LiDAR-derived), and the infrared band variance (RapidEye derived) (VIF < 1.5). Each coefficient represented in average 49.04, 32.47 and 18.48% of the explained variability, respectively (Figure 3).

The residuals are distributed normally (Kolmogorov-Smirnov test, p = 0.675), and spread randomly around the 0 line indicating homogeneity of error variance, no residual is visibly away form the random pattern of residuals (Figure 4). The correlation between the observed and predicted values was 82.3% (Figure 5).

3.4. Composition Mapping

The $\beta$ -diversity values ranged from −1.13 to 6.94 but highly concentrated between −1 and 2 (Figure 6). The few values higher than 2 were found along the road located at the southeast corner of the study area, followed by the south vegetation patch which has the most recent disturbance history (slash and burn agriculture). While the lowest values were observed in the areas that just experienced extensive logging. This patch seems to share space with the extensive logging plus forest fire patch, but in general, the latter area showed higher predicted values (Figure 7).

4. Discussion

The ordinated biplot (Figure 2) demonstrated less dissimilarity in tree species composition among plots that belong to a common vegetation patch, observable in the clustering effect by disturbance history. It was efficiently extended into the modeled geographic domain (Figure 7), where we see less dissimilarity in $\beta$ -diversity values between pixels belonging to the same disturbance history, as a confirmation of the reliability of the $\beta$ -diversity map (summarized in a floristic gradient) for studying in a continuous way, how dissimilar is the tree species composition between two pixels, or even between stands with different forest management, land use history or successional state, as other studies had found [77].

The discretization of the plots according to the negative or positive nature of the $\beta$ -diversity scores has been reported before, also relating biomass mean values to certain classes [10]. The previous study and other related [11,45] had found a clear link between community taxonomic diversity and carbon stock through patches with human disturbances. Our findings support that relation, because the greatest differences in biomass and beta diversity were found between the same areas. In other words, the scores higher than 0.5 were associated with plots in vegetation patches with extensive logging antecedent (higher biomass, Table 1), while lesser and negative scores, were related to the plots with slash and burn agriculture history (lower biomass, Table 1).

Even when a large number of remote sensing variables achieved statistically significant correlations, as suggests the assumption of the variation in the tree species composition will be related to variations in canopy traits such as reflectance spectra or canopy form [26,48], the relative superiority of the RapidEye derived metrics were evident. Other research that compares the performance of multisensor metrics for predicting compositional ordination [26] and taxonomic indices of biodiversity [20] also found the textural variables of multispectral imagery more useful. Those metrics by themselves has proved to be an invaluable tool for the study of biodiversity, as other studies had been reported in temperate [14,25] and tropical landscapes [39,78].

Between the textures, the infrared band has been traditionally considered relevant for vegetation classification and monitoring, because the light waves of this region of the electromagnetic spectrum are reflected in order to avoid the overheating of the leafs [79]. A large variance in this wavelength may indicate low uniformity in the exposed foliar tissue, in the vigor of the vegetation, or in the species mixing. The redEdge band were also found significant, its mean values was found to be dependent on chlorophyll concentration and the structure of the leaf and canopy [60,61].

The LiDAR strongest correlation (Canopy relief ratio) is a quantitative descriptor of the relative shape of the canopy surface inherited from topography studies previous to the LiDAR blooming, defined as the mean height minus the minimum height divided by the maximum height minus the minimum height [80]. This metric reflects if the canopy is mostly in the upper or lower portions of the height range, and has also been considered as important variable on other models of community composition [28] and biodiversity [20,45].

It is relevant to notice that the model selected variables are not exactly the highest correlated variables with the $\beta$ -diversity, as another study reported [23], partially because of the high collinearity among the $\beta$ -diversity strongest related variables, evident in the rise of the variance inflation factor of the models including them. Actually, even with the relatively lower correlation of the LiDAR metrics, the addition of the CHM variable to the model improved its predictive power, representing almost a third of the total explained variability. Height variables had been considered before as key measures for evaluating variations in community composition [11], and in this study, the canopy height was especially momentous for distinguishing the dissimilarity in tree species composition between the shorter areas with agriculture disturbance history, and the taller areas with selective logging disturbance history (Table 1).

Most studies mapping $\beta$ -diversity as a floristic gradient through ordination techniques have found similar results. Many considered LiDAR descriptors, whether in joint with hyperspectral data (random forest modeling, temperate forest, $R^2$ = 0.71) [22], or just LiDAR metrics (linear modeling, tropical rainforest, $R^2$ = 0.71) [23], and others had focused on multispectral textures with LiDAR (random forest modeling, temperate forest, $R^2$ = 0.61) [24], and without LiDAR (linear modeling, tropical forest, $R^2$ = 0.60–0.74 [25], kriging-regression, tropical forest, $R^2$ = 0.71) [19], or even with imaging spectroscopy alone (linear modeling, temperate forest, $R^2$ = 0.67) [26]. The major amount of studies of community composition modeling without using LiDAR metrics are addressed to grasslands. Those also concentrate the major variation in calibration; even though some studies considering hyperspectral data presented similar results (partial least squares modeling, $R^2$ = 0.71) [21,27], or (partial least squares modeling, $R^2$ = 0.74) [30]. Others showed much less fitting (random forest modeling, $R^2$ = 0.36) [28], or a great improvement (partial least squares modeling, $R^2$ = 0.92) [29].

Considering the background, it is evident that tropical landscapes have relatively few pieces of research [19,23,25]. And as far as we know, there are no record of sensors synergic use for mapping $\beta$ -diversity in the tropics, being a specially important tool for monitoring Latin American tropics recovery, a region dominated by second-growth forest [81]. Even when as a result of the synergy, we derived many covariates, we have attached to the Ocamm’s Razor principle of parsimony because, beyond the computationally demanding issues, for forestry applications like this, each argument is a remotely sensed covariate that may function as tacit reference to a process involved in the phenomenon under investigation, in our case, the variation in tree species composition. In this study, just three remote sensing variables were enough to map almost the same variability represented by the models of the literature reviewed, even when many of the suggested models required dozens of covariates [21,22,24,27,28,29,30].

Besides that, most of the studies suggesting few-variables models are focused on just one sensor, allowing them to use efficiently the stepwise variable selection method [19,23,25]. However, in genomic analysis or cases like this, the number of predictors (p) is larger than the number of records available (n), and if we think the prediction problem from the regression perspective, the parameters are not likelihood identifiable neither from the classical nor the Bayesian point of view in statistics [82]. So, we didn’t find another remote sensing research using the pre-selection variable strategy that we used in this study. None adjusted one linear mixed model for each covariate for assessing its individual significance and rank them. That was meaningful because each model includes as fixed effect the covariate to be tested, and a random effect that takes into account that there exists other covariates, in order to avoid issues of colinearity.

Finally, we want to highlight that one of the main reasons for considering forest structure and $\alpha$ -biodiversity (as species richness or biodiversity indexes) instead of dissimilarity in community composition, is the relative complexity of the latter, but nowadays it is recognized the potential of this metric, assessed as simple as the first axis of nMDS ordination, like a sensitive and effective tool for monitoring recovery on tropical managed landscapes [11,26,59]. One of the research topics that ought to consider $\beta$ -diversity standardly, is the tropical forest restoration. There are studies promoting passive restoration as an effective or better alternative to active restoration in the tropics [83], and at global scale [84]. Their conclusions are based primarily on ecological metrics of abundance and species richness, and forest structure that take far less time to achieve restoration success than dissimilarity on specific composition, particularly for canopy trees due to their slow turnover time [12,13,85,86]. Otherwise, most of the research that considered community composition encourage active restoration in order to assist more complete and faster recovery of old-growth biodiversity [87,88], highlighting that differences are long lasting. But even highly disturbed forests could, with appropriate management, provide important opportunities for better conserving and using tropical forest [89]. Aside from them, other research considering species richness, biomass, and composition dissimilarity, had found passive restoration fast and effective [90]. More studies that take into account dissimilarity in tree species composition over succession are needed.

5. Conclusions

Our findings support compositional ordination as a reliable tool to represent the general dissimilarity pattern of the whole species assemblage in a single metric of $\beta$ -diversity (nMDS axis-1 scores), with significative differences among the means of the scores corresponding to vegetation patches with different disturbance histories.

The hypothesis were confirmed. The floristic gradient (nMDS-1) variability is strongly correlated with some LiDAR metrics and RapidEye textures. It could be mapped for the whole study area when knowing the values of CHM, the mean reflectance of redEdge channel, and the variance of the infrared one. Besides, the areas with common disturbance histories have lesser inner dissimilarity.

This strategy for assessing $\beta$ -diversity, could be an important tool to monitor and understand forest succession, as to advise forest management decisions makers, with spatially explicit information.

Author Contributions

Conceptualization, A.d.P.O.-F., J.R.V.-L., G.Á.-P., J.I.V.-H., H.M.d.l.S.-P.; Statistical analyses, A.d.P.O.-F., J.L.H.-S., P.P.-R.; Original draft, A.d.P.O.-F.; Review and editing, J.R.V.-L., G.Á.-P., J.I.V.-H., J.L.H.-S., H.M.d.l.S.-P.

Funding

Project 12-IJ-11242306-033: Reducing Greenhouse Gas Emissions and Improving Forest Management in Mexico.

Acknowledgments

To the people from the ejido Felipe Carrillo Puerto and the Civil Association Uyoolche for granting us permission and facilities to work on their property. Especially thanks to Filiberto Yam Buenfil and José Arreola Palacios. We also thanks to the Mexico’s National Forest Comission through the Gerency of the National System of Forest Monitoring for providing the RapidEye imagery, and APOF thanks to the National Council of Science and Technology for the masters grant. This research was financially supported by the Sustainable Landscapes Program of the Agency for International Development of the United States of America, through the USDA Forest Service International Programs Office and the Northern Research Station.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Flowchart representing the whole procedure. Log-transformed percentage Importance Value Indexes (IVI%) of the species at plot level were ordinated onto a 2-D non-metric multidimensional scaling ordination (nMDS). Differences between disturbance histories were tested for each axis-scores, and the axis sensitive to that was linearly modeled. As predicted variables, textures of the five bands and two vegetation indexes were calculated from RapidEye Multispectral imagery and 95 LiDAR-derived variables. Spearman’s correlations were calculated between the response variable and the predictor variables, and the variable selection was performed using iteratively a linear mixed model that includes as fixed effect the covariate to be tested, and a random effect that takes into account that there exists other covariates that needs to be considered. We finally rank the variables from the lower p-values, and apply a stepwise selection procedure to the top ten covariates. Linear dependencies were explored by computing the variables’ variance inflation factors (VIF). The relative importance of LiDAR and RapidEye derived variance was assessed. Model predictive ability was assessed via ten-fold cross-validation and the final model was used to map it for the study area extent, each disturbance history with a characteristic position in the compositional ordination space (β-diversity).

Figure 2. Non-metric Multidimensional Scaling (nMDS) ordination scatter plot, for the study plots by disturbance history. Stress value = 0.23.

Figure 3. Relative importance for nMDS-1 with 95% bootstrap confidence intervals, b4 = redEdge band, b5 = infrared band.

Figure 4. Residual vs. fitted values.

Figure 5. Cross-validated observed versus predicted values for all the repetitions.

Figure 6. Distribution of the nMDS-1 predicted values.

Figure 7. β-diversity predicted values, plots are gray dots, and disturbance histories are labeled. Map proyected in UTM16 N.

Table 1

Vegetation patches characterization by disturbance history.

Attribute	Agriculture	Selective Logging	S Logging + Fire
Tree density	$669\pm 16$ stems/ha	$988\pm 4$ stems/ha	$1306\pm 2.5$ stems/ha
Basal area	$10.7\pm 19$ m ${}^2$ /ha	$30.6\pm 6$ m ${}^2$ /ha	$32.6\pm 3.6$ m ${}^2$ /ha
80p Biomass	$164\pm 18.2$ kg/ind	$261.7\pm 9.1$ kg/ind	$215.5\pm 6.7$ kg/ind
Height	$9.5\pm 5.1$ m	$14.2\pm 1.5$ m	$13.8\pm 2$
Richness	$24\pm 10$ species	$37\pm 3.6$ species	$45\pm 3.8$ species
Shannon index	$2.5\pm 10$	$2.8\pm 2.8$	$3.1\pm 2$
IVI% dominant	Bursera simaruba	Gymnanthes lucida	Gymnanthes lucida
species	Psidia psipula	Manilkara zapota	Bursera simaruba

IVI% = Importance Value Index of the two most important species, 80p = 80th percentile.

Table 2

Spearman correlation coefficient ( $r_s$ ): Remote sensing predictors versus nMDS-1 ordination scores. All correlations were statistically significant (p < 0.0001).

Textural Variable	${\mathit{r}}_{\mathit{s}}$	LiDAR Variable	${\mathit{r}}_{\mathit{s}}$
b3 variance	0.780	Canopy relief ratio	−0.616
b3 mean	0.767	Height L3 moment	0.547
b2 variance	0.747	Height skewness	0.537
b2 mean	0.738	All above mean/first × 100	0.537
reNDVI mean	−0.736	Height L skewness	0.535
reNDVI variance	−0.735	Canopy height model	−0.517
b4 variance	0.721	Cover 7m threshold	−0.510
b4 mean	0.704	All returns above mean	0.461
b5 variance	−0.440	Height at 80th percentile	0.437

b2 = green band, b3 = red band, b4 = redEdge band, b5 = infrared band.