A framework for near-real time monitoring of diversity patterns based on indirect remote sensing, with an application in the Brazilian Atlantic rainforest

a. Biological datasets and maps of species richness and phylogenetic diversity

To illustrate our biodiversity monitoring model and evaluate its sensitivity to the geographical range of the group, we first applied it to nine groups that differ in their spatial distribution (five groups of plants, and four groups of animals; Table 1). The plants included a group within the Bromelioideae subfamily (Aguirre-Santoro, Stevenson & Michelangeli, 2016; Aguirre-Santoro, 2017; Brown et al., 2020), three groups in the family Melastomataceae [in the Miconieae tribe (Goldenberg et al., 2008; Michelangeli et al., 2008; Caddah, 2013; Reginato & Michelangeli, 2016; Brown et al., 2020), in the genus Bertolonia (Bacci et al., 2020; Paz et al., 2021), and in the tribe Cambessedesieae (Bochorny et al., 2019; Paz et al., 2021)], and a group including Fridericia and allies (Kaehler, Michelangeli & Lohmann, 2019; Paz et al., 2021), in the plant family Bignoniaceae. The animals included the clearwing butterflies that compose the tribe Ithomiini (Brown et al., 2020), treefrogs of the genus Boana (Vasconcelos et al., 2014; Brown et al., 2020), horned frogs of the genus Proceratophrys (Brown et al., 2020), and a group of birds in the subfamily Thraupinae (Burns et al., 2014). These groups and their species differ in their distributions: some have small individual species distributions and small total geographical spread (like the Bromelioideae plants and the Boana frogs), others have broad individual species distributions and broad group-wide spatial coverage (like the Thraupinae), while others have discordant patterns (like the the Fridericia and allies, with small individual ranges but a broad group coverage; see Table 1, sorted from broader geographic distribution to smallest).

For each group, a map of species richness was downloaded from the Dryad repository of Paz et al. (2021); here, we provide a brief description of the methods used in that publication. Those maps were originally based on individual species distribution maps developed at a 10 km resolution (Paz et al., 2021), and created through the superimposition of alpha hulls representing the ranges of each species. Alpha hulls are a generalization of convex hulls, allowing for more shape flexibility, with lines that are allowed to curve and match the point distributions more closely (Meyer, Diniz-Filho & Lohmann, 2017). Once alpha hulls representing the ranges of each species in group were obtained, they were stacked and summed to obtain a species richness map, which was used in the present analysis (Paz et al., 2021). As any outline-based maps, alpha hulls may overpredict individual species ranges; however, they have been shown to do so less often than other commonly used techniques, such as minimum convex. Although they may lead to overprediction of species richness, this issue is more prevalent at smaller resolutions (Graham & Hijmans, 2006; Hurlbert & Jetz, 2007).

For each group, a map of phylogenetic diversity, as measured with Faith’s PD (Faith, 1992), was also downloaded from Paz et al. (2021). To obtain that map of the spatial distribution of phylogenetic diversity, the species richness map and community composition information of each group were previously combined with phylogenetic information describing the evolutionary relationships among all species in the group, using the software Biodiverse (Laffan, Lubarsky & Rosauer, 2010) as described in Paz et al. (2021).

b. Environmental descriptors

Environmental data to train and project the near-real time models were obtained from remote sensing sources and used to create bioclimatic layers describing the conditions experienced at different time windows. For that, we collected Normalized Difference Vegetation Index (NDVI) data and temperature data from NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS; Vermote & Wolfe, 2015) instrument, and precipitation data from the Climate Hazards Infrared Precipitation with Stations (CHIRPS; Funk et al., 2015) dataset. The MODIS data products are available bi-weekly since 2002, at a minimum resolution of 500m and 1km respectively; CHIRPS data are available monthly between 1981 and the present, at a minimum resolution of 0.05 degrees (∼5.5 km). Using a custom-created python script for Google Earth Engine (GEE, Gorelick et al., 2017), we obtained all NDVI and temperature MODIS data and CHIRPS precipitation data representing conditions experienced in the Atlantic Forest from 2002 to the present, using the forest shapefile from Paz et al. (2021). In GEE, we summarized climate data as monthly means, maximum and minimum values and NDVI data as yearly medians at a 10 km resolution. All summary layers were downloaded into R for further processing and for model building. Based on those monthly data, sets of 19 bioclimatic layers representing summaries of temperature and precipitation (Table 2) were built using the bioclim function of the dismo package (Hijmans et al., 2011). To calibrate the machine learning model, we used built bioclimatic and NDVI layers for the 2002–2014 time period; for projections of the model, we used yearly layers from 2015 to 2019. The python script to summarize and download data from GEE, as well as the R scripts used to generate the bioclimatic layers, are available on the Github repository of this manuscript (https://github.com/andrepazv/NearRealTime). Before proceeding to modeling, we removed all highly redundant (correlated) variables (VIF > 10) using the vifstep function of the usdm R package (Naimi et al., 2014). This reduced the environmental dataset to ten variables that were used in all further analyses (Table 2).

c. The machine learning model and time projections

To predict the spatial distribution of species richness and phylogenetic diversity as a function of environmental change in near-real time, we independently modeled species richness and phylogenetic diversity with data describing temperature and precipitation. Those models with the best fit were then projected over time, yearly, based on each year’s environmental conditions.

We first trained an ensemble machine learning model to describe species richness, and an ensemble machine learning model to describe phylogenetic diversity, for each of the nine groups (Table 1). For that, we used bioclimatic layers built from MODIS and CHIRPS and NDVI estimates (Table 2) summarized for the 2002–2014 time period. To improve model performance, we performed an ensemble of four machine learning algorithms: Random Forests (rf from Liaw & Wiener, 2002), Neural Networks (nnet from Venables & Ripley, 2002), support vector machines (svmRadial from Karatzoglou et al., 2004), and generalized linear models. These algorithms allow the prediction of continuous variables (as a form of regression), but the first three make no assumption of linearity in the relationship between predictor and response variables, thus allowing for more complex and biologically relevant relationships. Seeking good transferability across time projections, we used a spatial partition method for model tuning and validation (Radosavljevic & Anderson, 2014): for each model, we partitioned the biodiversity data (pixels) into geographic regions, instead of randomly. To do that, for each biodiversity metric (species richness or phylogenetic diversity) and taxon (each one of the nine groups of organisms), we created four spatial partitions using an adapted version of the function get.block in the package ENMeval, partitioning the data into bins across latitude and longitude lines (provided with the full script on GitHub; Muscarella et al., 2014). For each one of the four machine learning algorithms, we employed k-fold cross validation, building k = 4 models, each leaving out one different spatial partition that was used for model validation. Models were trained and tuned using the caretList function in the caretEnsemble R package, with partitions passed through the trainControl argument (Deane-Mayer & Knowles, 2019). Individual tuned models were selected using the root mean square error (RMSE) metric. Using the settings determined as optimal in the validation stage, we then built a full model with all data (full model from hereon). Finally, an ensemble model was built for each organismal group and biodiversity metric, based on a linear combination of the models for all four algorithms, weighted according to their RMSE values. For that, we used the caretEnsemble R package and repeated 10-fold cross-validation (Deane-Mayer & Knowles, 2019). The script used for modeling is available in the GitHub repository for this manuscript.

These models allow for predictions of the response variable (here, species richness and phylogenetic diversity) to be transferred to a subsequent time period based on the same predictor variables (the bioclimatic variables + NDVI, Table 2, Fig. S4). To obtain an estimate of trends in biodiversity change over time, we projected the models in time over five different years, generating species richness and phylogenetic diversity projections for 2015, 2016, 2017, 2018 and 2019. For this, we applied the corresponding model to each year’s predictor variables. Using the model and the five projections over time (from 2015 to 2019), we stacked the maps together and computed the slope of a linear regression for each pixel to estimate trends in biodiversity change for each group and biodiversity metric. To obtain a spatial estimate of the expected mismatch between the biodiversity metrics, we computed the regression of phylogenetic diversity onto species richness for each year, and mapped the residuals of the regression for each year.

d. Ground truthing the modeling framework

We then used another dataset to ground-truth the model framework, contrasting the biodiversity change trends inferred by the model against trends documented in the field. For that, we used information available for two other groups of organisms—epiphytes and amphibians—for which survey data are available across multiple years and across a sufficient percentage of their ranges, using the ATLANTIC-epiphytes (Ramos et al., 2019) and ATLANTIC-amphibian datasets (Vancine et al., 2018). These datasets were made possible by a multi-institutional effort that gathered inventory data for different plant and animal groups, collected across several decades along the Atlantic Forest. Although geographical sampling along the years is sparse, these two datasets have survey information that spans most of the latitudinal extent of the Atlantic Forest and includes contemporary sampling that partly matches the available satellite information. The epiphyte dataset was composed of 89,270 observations from 1824 to 2018 (Ramos et al., 2019). That of the amphibians was composed of 17,619 records, with data from 1940 up to 2015 (Vancine et al., 2018). For both datasets, we created maps of species richness for both datasets at a resolution of 100 km. This resulted in maps of 188 and 137 pixels, respectively, with species richness per pixel varying from 1–365 epiphyte species and from 1–80 amphibian species. Because the ATLANTIC datasets do not include phylogenetic data, this analysis focused solely on species richness and not phylogenetic diversity.

For the purposes of ground-truthing the model, we limited our analyses to the years for which the datasets (i) documented the year of the observations, (ii) provided coordinates for the observation localities and (iii) provided a relatively broad sampling coverage of the region (i.e., sampled sites were not clustered geographically). This resulted in the use of data collected up to the year of 2015 for the amphibians, and up to the year of 2018 for the epiphytes. We then used the available inventory data from all but the last four years of observations to generate maps of species richness (up to the year of 2011 for the amphibians, and 2014 for the epiphytes), measured as the number of species present in each pixel in each year. Following the same modeling procedure applied previously, we combined those species richness maps with bioclimatic variables derived from remote sensing sources as described above, and for the same period, to train a predictive model. The model was then projected onto the four following (most recent) years, for which observational data are available for ground-truthing. As such, we obtained yearly model projections of amphibian species richness for the years of 2012, 2013, 2014, and 2015, and yearly model projections of epiphyte species richness for the years of 2015, 2016, 2017, and 2018. As done previously, we used these yearly model projections to create maps of inferred trends in species richness across time for each of the target groups, computing the slope of a linear regression for each group. We then compared these modelled trends with the observations recorded in the datasets, evaluating both directionality and magnitude of change. To evaluate if results were better than random predictions, we created 10,000 random predictions of gain and loss across the forest, and, for each one, computed the number of rightly predicted pixels. We used those results to create a null distribution of random predictions and compared the observed number of correct predictions with the null distribution. A script used to run the entire framework (in R and python) is available in the GitHub repository for this manuscript.

Model applicability and sensitivity to range size

The analysis of the first nine groups shows that prediction of biodiversity metrics from environmental variables describing temperature, precipitation and canopy cover during the years of 2002–2014 varied substantially in fit across datasets (species richness: R² values 0.18–0.83; phylogenetic diversity: 0.08–0.74; Table 3). In most cases, the fit of the ensemble model was better than that of any individual (tuned) model, with the fit of each individual model being similar across algorithms. In a few exceptions, however, the fit of the individual model was poor—in those cases, the fit of the ensemble model was even worse (Table S1).

Model fit increased with the total area of distribution of the group, when predicting both species richness and phylogenetic diversity (Fig. 1A). However, median individual species range does not appear to be related to model fit (Fig. 1B). Five out of the nine groups had well-fit models for both diversity measures (R² >  0.5): they included two plant groups (Miconieae and Fridericia & allies) and three animal groups (Thraupinae birds, Ithomiini butterflies and Boana frogs, Fig. 1, Table 3). These three groups are known for having the largest total geographical range in the Atlantic Forest (1,197,890 km², 1,282,733 km², 1,330,613 km², 1,147,465 km² and, 824,368 km² respectively), but they have variable individual species ranges (median individual species ranges 44,486 km², 85 km², 680,022 km², 321,331 km² and 240 km² respectively, Table 2). Models of all four remaining groups had an R² smaller than 0.49 for at least one of the two diversity measures (Table 3). All four groups had comparatively small total ranges in the forest, but variable individual ranges (Table 2).

In the five groups for which the models were able to predict species richness with good accuracy (Miconieae, the Fridericia & allies, Thraupinae birds, Ithomiini butterflies and, Boana frogs; species richness R² values of 0.71, 0.83, 0.72, 0.56 and, 0.53 respectively; Table 3), the most important predictor variables were all related to precipitation. In Miconieae and the birds, precipitation of the warmest quarter (bio 18) stood up (explaining 25% and 18% of the variation respectively, Table S2). For Fridericia and the Ithomiini butterflies, precipitation of the wettest month (bio 13) was most important variable (explaining 18% of the variation). Finally, for the Boana frogs, precipitation of the driest month (bio 14) was most important (22%, Table S2). Patterns of phylogenetic diversity in these groups were similarly well correlated with local environments (R² values of 0.74, 0.71, 0.68, 0.59 and, 0.52 respectively) and, again, the most important predictor variables were all related to precipitation. Precipitation of the warmest quarter (bio 18) was the main predictor of phylogenetic diversity in the Miconieae and the Boana frogs (25% and 22% respectively), whereas precipitation of the wettest month was the main driver for the other three groups (15% for the Fridericia, 21% for the Thraupinae and 17% for the Ithomiini, Table S2).

Projections in the five groups with good model fit show similarities and differences in the predicted trends of change in species richness and phylogenetic diversity from 2015 to 2019 (clusters of gain (red) and loss (blue) of predicted diversity across the forest; Fig. 2). In all groups, and for both diversity metrics, a trend toward diversity loss over time is identified in the central region of the forest, in the north of the São Paulo subtropical gap located between the coastal mountains (Amaral et al., 2018; Thom et al., 2020), and into the inland state of Minas Gerais. This is particularly strong in Fridericia & allies (Fig. 2, second row). On the other hand, a trend toward diversity gain (both in species richness and phylogenetic diversity) is detected along the coast, stretching northward from the state of Rio de Janeiro to Espírito Santo (Fig. 2). Models of all groups also detect trends of increasing diversity for both biodiversity metrics in the southern range of the forest (Fig. 2). Maps of the yearly predictions are available in the supplemental material (Figs. S1 and S2).

For all five groups, the predicted distribution of the residuals of the regression of phylogenetic diversity and species richness varies across space, and, to a lesser degree, over time (Fig. 3). In the Miconieae, the models suggest that more phylogenetic diversity is expected to accumulate given species richness (positive residuals expanding across time; Fig. 3, top row green colors). However, two areas of lower phylogenetic diversity are flagged over time: one to the northwest, where the deciduous forest is present (Peres et al., 2020), and a second one, along the coast, in the Serra do Mar region (Fig. 3, top row brown colors). In Fridericia & allies, the projections mostly suggest that phylogenetic diversity is expected to accumulate relative to species richness except for sites in the northwest (similar to Miconieae) and a large region in interior São Paulo (Fig. 3, second row, brown colors). The latter is the same region where both Miconiaeae and Fridericia & allies are inferred to undergo an increase in both diversity metrics over time. The animal groups (Thraupinae birds, Ithomiini butterflies and Boana frogs) differ from the two plant groups in the predicted pace of the changes; they appear less drastic over time (Fig. 3). The bird and frog models detect a large region accumulating more phylogenetic diversity than expected in the south of the forest, and suggest this trend to be mostly stable across time except for 2019 (Fig. 3, third and last row, green colors). These two groups also show a consistent prediction of lower phylogenetic diversity than expected in the north of the forest and a somewhat consistent trend in interior Sao Paulo (Fig. 3, third and last row, brown colors). For the butterflies, an opposite pattern emerges—with the south consistently predicted to have less butterfly phylogenetic diversity than expected, and the north having more phylogenetic diversity than expected (Fig. 3, fourth row, brown and green colors respectively).

Lessons from model testing

When the model framework is applied to inventory data that include multiple records per year (enabled by the ATLANTIC-amphibian and ATLANTIC-epiphytes datasets), this framework predicts the directionality of biodiversity change well. Although the predictive power of the amphibian and epiphyte models is lower than those observed for the groups analyzed in the first part of the study (R² of 0.39 for epiphytes, 0.21 for amphibians), the modelled trends of potential change in species richness are broadly congruent with those measured with empirical data, particularly regarding species loss (Fig. 4). A comparison of the modelled vs. observed trends in gain or loss of species richness indicates that the model framework correctly predicted the direction of change in 59% of the pixels (52/88 pixels, p-value = 0.08) in the epiphyte study, and 76% of the pixels (28/37 pixels p-value = 0.0008) in the amphibian group (Fig. 4). Even though richness was described using larger pixels in the ATLANTIC datasets relative to the nine other groups (100 km instead of 10 km), not all pixels contained enough information to compute observed trends; as such, fewer pixels are available for comparisons (88 for the epiphytes, 37 for the amphibians).

These modelled trends of change differ between the two groups and are spatially heterogeneous, with a tendency towards loss of species. A total of 46/88 pixels in the epiphytes model and 20/37 pixels in the amphibian model have a negative slope. The same tendency towards species loss is present in the observed data, although more strongly, especially in the epiphytes (79/88 pixels in the observed epiphyte data, and 26/37 pixels in the amphibian data; Fig. 4). While the predictive models of the epiphytes generally matched the observed trends of species loss, the magnitude of loss was generally underestimated; observed slopes ranged from -104 to 19, while modelled slopes varied from −75 to 28 (Fig. 4, bottom). In the amphibians, the patterns of biodiversity inferred by the model matched the direction of change reasonably well, with a similar estimated strength of change (observed slopes ranging from −5.62 to 5; modelled slopes spanning from −6.7 to 5.5, Fig. 4 top).