Occurrence-based models reveal greater spatial details in tropical species richness

Species richness maps are essential for biodiversity assessment and conservation planning, but their accuracy is limited by uneven sampling effort and large data gaps, especially in tropical regions. Here, we combined the Uniform Sampling from Sampling Effort (USSE) framework with deep neural networks to generate occurrence-based richness maps for tropical vertebrates, invertebrates and plants. The USSE approach incorporates sampling intensity during model calibration and projects richness under a standardized high-sampling scenario, reducing the influence of spatial sampling bias. Model predictors included climatic, topographic and remotely sensed habitat-conditions, allowing richness predictions to reflect both environmental gradients and habitat alteration. Predictive performance was high across taxonomic groups, with spatially blocked cross validated R² values ranging from 0.60 to 0.95, and a value of 0.90 for the combined multi taxon model. Independent validation against 107 published field inventories yielded R² values between 0.55 and 0.87. For terrestrial vertebrates, predicted richness patterns were broadly consistent with International Union for Conservation of Nature polygon-based maps in areas of high ecological integrity, while revealing greater fine-scale spatial heterogeneity associated with habitat variation. Across groups, richness was generally highest in tropical forest regions. Similarity analyses indicated that invertebrate groups, particularly Lepidoptera, showed strong agreement with overall multi-taxon richness patterns. Our results show that USSE combined with deep learning provides a scalable framework for producing spatially detailed biodiversity maps that account for sampling bias and current habitat conditions, offering a complementary tool for conservation prioritization in tropical regions.

Species richness, the number of species present in a given area, is a fundamental dimension of biodiversity (Burgess et al., 2024) and a central focus in macroecology, biogeography, and conservation (Wiens and Donoghue, 2004; Rangel et al., 2018; Burgess et al., 2024). Mapping species richness is used in many biodiversity conservation approaches (Fleishman et al., 2006; Burgess et al., 2024), making accurate mapping essential for robust conservation studies and proposals. However, large biological sampling gaps and sampling bias hinder the accurate mapping of richness (Oliveira et al., 2016; Cardoso et al., 2024), resulting in some regions having little or no available biodiversity information (Bini et al., 2006; Brito, 2010; Oliveira et al., 2016). These challenges are particularly pronounced in the tropics (Deikumah et al., 2014; Oliveira et al., 2016). Tropical ecosystems are highly heterogeneous, ranging from evergreen and seasonally dry forests to savannas, shrublands and deserts, so that richness varies steeply over short distances. Most primary biodiversity information for these ecosystems is mobilized by infrastructures such as the Global Biodiversity Information Facility (GBIF), which aggregate occurrence records across taxonomic groups and institutions. These records provide unprecedented taxonomic and geographic breadth, and their georeferencing allows each occurrence to be associated with independent environmental data, such as climate surfaces and remotely sensed descriptors of habitat structure. Their accumulation across space is nonetheless uneven, concentrated near roads, research institutions and accessible areas (Deikumah et al., 2014; Oliveira et al., 2016).

Historically, data limitations have led to the widespread use of richness estimates based on distribution polygons, which are produced by specialists interpolating occurrence points using field knowledge and expertise (Hurlbert and Jetz, 2007; Anderson, 2023). Because these polygons are drawn by interpolation between known records rather than inferred from environmental conditions, their spatial grain is coarse and largely independent of the environmental heterogeneity contained within the delimited range (Mainali et al., 2020; Broekman et al., 2022; Marsh et al., 2022).

Summing up of these individual polygons yields richness maps used by databases such as the International Union for Conservation of Nature (IUCN) (IUCN, 2025). Yet this approach has notable limitations: polygon-based maps are usually derived from occurrence records, expert knowledge, and geographic interpolation, but they may incorporate subjective biases and generalize species presence across unsuitable or degraded habitats (Mainali et al., 2020; Broekman et al., 2022; Marsh et al., 2022). When multiple species maps are stacked to estimate richness, these individual uncertainties and commission errors can accumulate, thereby increasing uncertainty and potentially overestimating local richness at fine spatial scales (Mainali et al., 2020; Broekman et al., 2022; Marsh et al., 2022). Such maps also often fail to reflect current ecological conditions, as they do not account for habitat degradation or differences in richness between intact and impacted areas, an essential consideration for conservation strategies (Ferrier et al., 2010). This limitation is inherited rather than incidental, since neither the occurrence databases from which these maps are built nor the maps themselves carry information on the current condition of the habitats they cover, so the resulting richness estimates remain static with respect to deforestation and degradation (Mainali et al., 2020; Broekman et al., 2022; Marsh et al., 2022). Additionally, databases like the IUCN are taxonomically biased, focusing predominantly on vertebrates and largely excluding more diverse and less-studied groups such as plants and invertebrates (Seddon et al., 2005; McRae et al., 2017; Titley et al., 2017).

An alternative methodological approach involves deriving richness maps from Species Distribution Models (SDMs), which sum the predictions of individual species models that relate known distributions to environmental predictors (Cord et al., 2014; Pouteau et al., 2015; Schmitt et al., 2017). However, this method can accumulate errors from individual models, potentially reducing the accuracy of richness predictions (Cord et al., 2014; Henderson et al., 2014; Vale et al., 2017). Moreover, both stacked SDMs and direct richness modeling approaches remain sensitive to sampling bias and collection gaps, as under-representation of environmental conditions in poorly sampled areas can bias model outputs (Hortal and Lobo, 2011; Henderson et al., 2014; Oliveira et al., 2019a; Alves et al., 2020).

Addressing sampling bias and gaps is crucial for improving species richness mapping (García-Roselló et al., 2023; Oliveira et al., 2016, Oliveira et al., 2024). The Uniform Sampling from Sampling Effort (USSE) approach incorporates sampling intensity as a predictor during model calibration, allowing the model to estimate how much of the observed richness pattern is associated with environmental conditions and how much is associated with uneven sampling effort (Oliveira et al., 2024). During prediction, the sampling-effort layer is replaced by a uniform high-sampling surface, so that richness is estimated as if all regions had been sampled with comparable intensity (Oliveira et al., 2024). Thus, USSE is not a specific algorithm, but a sampling-standardization framework that can be implemented with different regression or machine-learning models and has demonstrated the ability to mitigate sampling bias by producing more accurate richness predictions (Oliveira et al., 2024). In this study, we apply the USSE approach to model species richness across the tropics, aiming to minimize sampling-related biases and incorporate predictors of habitat condition derived from remote sensing, allowing predictions to account for richness changes associated with habitat alteration, such as deforestation and degradation, for more accurate richness patterns. We broaden taxonomic coverage by modeling richness for tropical vertebrates (birds, mammals, amphibians, reptiles), a diverse set of invertebrate groups (bees, spiders, cockroaches and termites, flies, dragonflies, terrestrial gastropods, hemipterans, lepidopterans, orthopterans, trichopterans, wasps), and plants (selected angiosperm families and all pteridophytes). Because IUCN polygon-based maps generally represent potential species distributions and do not explicitly incorporate current habitat condition, this integration of climatic and remote-sensing predictors represents an important advance for estimating spatial variation in richness under present environmental conditions. We also expected that the use of these predictors would allow the models to capture finer-scale spatial variation in richness than polygon-based maps. Additionally, we analyze the similarity of richness patterns among these groups and compare our model results for terrestrial vertebrates with IUCN polygon-based richness maps, both in terms of spatial consistency and in their ability to represent fine-scale map detail. By addressing sampling bias, current habitat condition and taxonomic breadth within a integrated framework, rather than treating them as separate problems, our approach yields richness surfaces intended for direct use in conservation prioritization, impact assessment and land-use planning, which are typically conducted at the local and regional scales where polygon-based products are little informative. The framework is furthermore transferable to other occurrence derived biodiversity metrics and to any region where collection effort is uneven.

Species occurrence data collection and handling

To compile species distribution data, we searched online biological collection databases, including the Global Biodiversity Information Facility (GBIF) (https://doi.org/10.15468/dl.4wcg7z, https://doi.org/10.15468/dl.mcg2u9, https://doi.org/10.15468/dl.wvqpbn) and the Global Ant Biodiversity Informatics (GABI) (https://antmaps.org). GABI was included to complement the relatively low coverage of ant occurrence records in GBIF and to improve the representation of Formicidae in the database. In these databases, we applied three filters: 1 - records from tropical regions, between 23.5°N and 23.5°S; 2 - records based only on specimens deposited in biological collections; and 3 - records restricted to the following organism groups: vertebrates (birds, mammals, amphibians, and reptiles); invertebrates (bees, spiders, cockroaches and termites, flies, dragonflies, terrestrial gastropods, hemipterans, lepidopterans, orthopterans, trichopterans, and wasps); angiosperms (represented by all families), and pteridophytes (all fern families).

We screened occurrence records using a multi-step geographic validation and cleaning procedure. For records with coordinates, we first removed points located in the ocean and then verified whether each reported point fell within the locality described in the record, based on the available administrative information, including municipality, city, state, or country. Records with coordinates falling outside the expected locality were flagged as geographically inconsistent and were georeferenced again. Whenever possible, georeferencing was based on other GBIF records from the same locality whose coordinates had already been validated; otherwise, locality descriptions were georeferenced using OpenStreetMap (https://www.openstreetmap.org/). All newly assigned coordinates were then rechecked using the same locality-based validation procedure. Records without coordinates were georeferenced following the same workflow, and records for which no reliable geographic assignment could be obtained were excluded from the database.Species names were validated against the Catalogue of Life Checklist 2023 (https://www.catalogueoflife.org/), which was used as a single taxonomic backbone to ensure nomenclatural standardization across all taxonomic groups analyzed, including vertebrates, invertebrates, and plants. Although taxon-specific databases may provide finer taxonomic resolution for particular groups, the use of a common reference avoided applying heterogeneous nomenclatural criteria across groups and allowed a standardized comparison of richness patterns. We also excluded nomina dubia (names of doubtful application). After these steps, the initial database of 46,469,994 records was reduced to 14,904,219 occurrence records representing 400,443 species (Supplementary Table S1). Record filtering and georeferencing were performed using the GBIF cleaning data tool in BioDinamica (Oliveira et al., 2019a).

Modeling species richness

Because species richness is strongly influenced by sampling intensity (Dennis et al., 1999; Loiselle et al., 2007; Fründ et al., 2016; Oliveira et al., 2016), it is essential to account for this factor in modeling richness patterns. We applied the Uniform Sampling from Sampling Effort (USSE) method (Oliveira et al., 2024), which uses sampling effort as a predictor in regression or deep learning models to reduce bias by simulating uniformly high sampling during prediction. This approach has been shown to minimize the effects of bias and sampling gaps in spatial biodiversity modeling (Oliveira et al., 2024). After model training, sampling effort is replaced by a uniform surface representing maximum sampling density (using a 100 km kernel), allowing the models to estimate biodiversity patterns based on environmental and sampling variables while discounting sampling effects. We calculated species richness within a regular grid of equal-area hexagons (42.681 km²) and used these as dependent variables in spatial modeling with the USSE approach and various algorithms (see algorithm testing section). We tested multiple hexagon sizes and selected 42.681 km² because it consistently yielded the highest predictive capacity across taxonomic groups (Supplementary Table S2). Hexagons with fewer than 100 samples were excluded to avoid undersampling. We also tested alternative minimum-sample thresholds and retained 100 samples because it produced the best average predictive performance across groups and algorithms (Supplementary Table S3). Predictor variables included climate, topography, a Composite Vegetation Index (CVI), and sampling effort (described below).

Biodiversity predictors

Predictors were selected to represent three well-established classes of correlate of terrestrial richness, applicable across the taxonomic breadth of this study. Climate was included because temperature and water availability constrain physiological tolerances and primary productivity and are the strongest broad scale correlates of richness for plants, vertebrates and invertebrates alike (Hawkins et al., 2003; Willig et al., 2003; Vázquez‐Rivera and Currie, 2015; Barreto et al., 2021). Topography was included because elevation and slope generate local climatic and habitat heterogeneity that interpolated climate surfaces do not capture. Vegetation structure, summarized by the Composite Vegetation Index (CVI) (Oliveira and Soares-Filho, 2026), was included because canopy height and the seasonal dynamics of greenness determine resource availability, microclimate and niche space, thereby constraining how many species can persist locally (Boucher-Lalonde et al., 2014; Ferger et al., 2014; Hyman et al., 2019). Moreover, it is the only predictor in our set that registers habitat conversion and degradation. Sampling effort was included as a methodological rather than an ecological predictor, as described above.

All predictor layers were assembled on a common 1 km² raster grid. Each layer was then aggregated to the hexagonal modelling unit by averaging the values of all raster cells whose centroid fell within the hexagon. Occurrence records were assigned to hexagons by spatial overlay, each record being attributed to the hexagon containing its coordinates. Response and predictors therefore share the same spatial units. Temporal matching was defined by reference period rather than by date-to-date correspondence. WorldClim variables represent long term climatic normals, whereas the CVI was derived from a single year, 2020, because the signature of land use and degradation on vegetation structure is only detectable at a defined point in time and would be blurred by multiyear averaging. The capacity of the index to recover these conditions was validated against independent land-cover data (Oliveira and Soares-Filho, 2026). Although occurrence records accumulate over a long historical period, the majority were collected recently. In addition, restricting calibration to cells with Ecosystem Integrity above 0.9 further reduces the residual mismatch between the date of a record and the habitat condition observed at its location.

Climate

To characterize climate, we used 19 temperature and precipitation bioclimatic variables from WorldClim (Fick and Hijmans, 2017) at a 1 km² resolution. To reduce multicollinearity and prevent model overfitting, we performed a principal component analysis (PCA) on these variables, summarizing them into axes that capture the main environmental gradients. The first four PCA axes, selected according to Kaiser’s criterion (Jolliffe, 2002), explained 89.02% of the variation and were used as climate predictors (Supplementary Table S4). PCA raster maps were generated using the BioDinamica package (Oliveira et al., 2019a) from the Dinamica-EGO software (Soares-Filho et al., 2013), was used to perform all analyses.

Topography

Two topographic variables were used at a spatial resolution of 1 km²: 1) Elevation (Danielson et al., 2011), and derived from it using its original spatial resolution, 2) Slope.

Composite vegetation index

Most studies focused on habitat characterization are conducted at a local scale, describing specific environmental features (Pappas et al., 2015; Lipsey et al., 2017; Nagy et al., 2019). Global approaches typically offer a qualitative characterization of habitats (Jung et al., 2020). Here, we used a Composite Vegetation Index (CVI) to describe quantitatively different types and structure of vegetation, also a measure of habitat changes due to deforestation, fire, and climate change. To build a quantitative index that directly qualifies vegetation types and structure, and human-induced changes, we used remote sensed data derived from Light Detection and Ranging (LIDAR) (Lang et al., 2023) and the Moderate-Resolution Imaging Spectroradiometer (MODIS) sensor. We used canopy height obtained by machine learning interpolation of LIDAR (Lang et al., 2023) data from 2020, as a direct measure of differences in vegetation structure and between natural and anthropogenic habitats. We also used the mean monthly variance per cell over the year 2020 for the Normalized Difference Vegetation Index (NDVI), which characterizes vegetation stress conditions, density, and seasonality, effectively distinguishing differences in seasonality, vegetation types, and anthropogenic land uses. For this purpose, we used monthly data obtained from NASA’s MYD09A1 (Vermote, 2015) product for the year 2020, at a spatial resolution of 1 km². To obtain the quantitative index, we performed a PCA analysis on these data. Vegetation types and structure was characterized by the first PCA axis significant according to Kaiser’s criterion (Jolliffe, 2002), representing, representing 70.61% of the observed variation in environmental variables (Supplementary Table S5).

To validate the ability of this index to characterize different vegetation types and anthropogenic changes, we tested whether the index would differentiate between the classes of the land-cover and land-use map of 2020, i.e., “The Global 2000–2020 Land Cover and Land Use Change” (Potapov et al., 2022), by using a generalized linear model (GLM). The CVI was the dependent variable, and the land-use classes were the predictor variables. The maps, originally with a spatial resolution of 30 meters, include as land-classes forests, cultivated areas, built-up areas, surface water, and perennial ice from 2000 to 2020. We resampled the maps to 1 km² spatial resolution using the Nearest Neighbor method. We chose the Nearest Neighbor method to preserve the original land-cover classification integrity so as to avoid artificial mixing of classes, given the fine-scale heterogeneity of the original 30m data. In the tests, the CVI showed distinct value ranges for different environments and land uses, as indicated by the GLM analysis, which yielded a pseudo-R² of 0.97 and a p-value < 0.0001. High values characterize forest vegetation, while lower values characterize savannahs and desert environments (Supplementary Figure 1). To rigorously assess pairwise differences between land-cover classes, we performed Tukey’s post hoc test. The results confirmed statistically significant distinctions (p < 0.0001) of CVI values across all class combinations (Supplementary Table S6), reinforcing the GLM results. Since the CVI was included as a predictor variable, and it shows distinct values for different land-uses, only samples located in areas with an Ecosystem Integrity (Hill et al., 2022) greater than 0.9 were employed for training to ensure that the biological diversity was correctly associated with the original habitat of the location. Because the CVI describes vegetation structure, The CVI is applied together with four climatic axes and two topographic variables, so that a natural savanna and a deforested rainforest cell occupy distinct positions in the predictor space even when their CVI values coincide, the former embedded in the climatic regime under which open vegetation is the natural state and the latter in the climatic regime of moist forest. The discriminatory capacity of the index across vegetation types and land uses is described in Oliveira and Soares-Filho (2026).

Sampling effort

We used sampling density to estimate the sampling effort for training the USSE method. The density was estimated using the kernel index, based on all the samples available in the species occurrence database. We used 100 km as the kernel radius. We tested alternative kernel radii for estimating sampling effort and selected the radius used in the analyses because it provided the best predictive performance in cross validation (Supplementary Table S7).

Algorithm testing for spatial modeling

To perform the USSE approach, which is a methodological framework for standardizing sampling effort and can therefore be implemented using different regression or machine-learning algorithms, we tested the models Deep Neural Network (DNN) (Schmidhuber, 2015), Generalized Linear Model (GLM) (Nelder and Wedderburn, 1972), Generalized Additive Model (GAM) (Hastie and Tibshirani, 1990), Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996), Multivariate Adaptive Regression Splines (MARS) (Friedman, 1991), Random Forest (RF) (Tin Kam Ho, 1998) and Support Vector Machines (SVM) (Cortes and Vapnik, 1995). We opted to test deep learning methods because our goal was to achieve the best possible predictions. Therefore, we tested even complex models like DNNs, as they were expected to provide superior prediction accuracy, which was indeed observed (Chang, 2023). The set of algorithms was chosen to span the methods frequently applied to continuous biodiversity responses. Rather than to be exhaustive, we aimed at assessing methods along a gradient of flexibility. This gradient comprises a parametric linear baseline (GLM), a penalized linear model performing variable selection under collinearity (LASSO), two approaches that relax linearity while remaining interpretable (GAM and MARS), two machine learning regressors with contrasting inductive biases, namely tree ensembles (RF) and a kernel-based regression (SVM), and finally a deep neural network capable of approximating high order interactions between predictors (DNN). Other approaches traditionally used to map richness, such as nonparametric estimators and stacked SDMs, were not tested because they are conditioned on the survey intensity or subjected to accumulated errors from the summed surface of species, as discussed above, whereas modelling richness directly with sampling effort as a predictor addresses both. As described above, predictors are represented in all algorithms as axes of correlation based PCA rather than as raw variables. This procedure removes multicollinearity and places inputs on a comparable scale, as required by the gradient based optimization of the DNN and by the SVM kernel. Elevation, slope and sampling effort were rescaled to the same range. Richness was modelled untransformed, as its distribution is approximately symmetric after excluding hexagons with fewer than 100 records, which conform to the Gaussian structure used in the GLM and GAM. In the GLM, a Gaussian adjustment was used, as the biodiversity index data are continuous, representing numerical values that are expected to follow distributions close to normal; in the GAM, a spline function was applied to each predictor variable to capture nonlinear relationships between them, and a Gaussian distribution was used because the biodiversity index data are continuous and expected to follow a distribution close to the normal; in the LASSO algorithm, cross-validation was used to estimate the optimal value of lambda for regularization; in the MARS, the maximum degree of allowed interactions was set to 3, with a maximum of 100 optimized models and 50 nodes to avoid model overfitting, using a penalty of 2. The RF algorithm used 300 random trees, and the SVM was configured for an epsilon-regression function. The regressive DNN employed five layers of 64 neurons, with layers 1 and 3 incorporating kernel regularization (=0.01), interspersed with 20% dropout layers to reduce overfitting, followed by another 64-neuron layer and a final layer with a single neuron. Backpropagation was used to adjust the weights efficiently, minimizing error during training. Internally, 20% of the data was used for validation during training. We tested different structures for optimizing the DNN in prediction, including:

A layer of 32 neurons, with 20% dropout, followed by a layer of 64 neurons and one with a single neuron.

A layer of 64 neurons, with 20% dropout, followed by a layer of 64 neurons and one with a single neuron.

Two layers of 64 neurons, one with kernel regularization (λ=0.01) and another with 20% dropout, followed by a layer with a single neuron.

Five layers of 64 neurons, with kernel regularization (λ=0.01) in layers 1 and 3, interspersed with 20% dropout layers, followed by a layer of 64 neurons and one with a single neuron.

A layer of 128 neurons with kernel regularization (λ=0.01), a layer of 128 neurons with 20% dropout, followed by two layers of 64 neurons (one with kernel regularization (λ=0.01) and another with 20% dropout), followed by a layer of 32 neurons and one with a single neuron.

A layer of 128 neurons with kernel regularization (λ=0.01), a layer of 128 neurons with 20% dropout, followed by five layers of 64 neurons, with kernel regularization (λ=0.01) in layers 1 and 3, interspersed with 20% dropout layers, followed by a layer of 64 neurons and one with a single neuron.

Model selection and validation

All algorithms were compared using the mean squared error (MSE) through the random partitioning of five datasets, with 80% for training and 20% for testing in complete dataset (all groups). We selected MSE over alternative methods such as the area under the receiver operating characteristic (ROC) curve (AUC) or true skill statistic (TSS) because our model deals with continuous values (biodiversity indices), necessitating a validation metric that not only quantifies predictive accuracy but also captures the variation across the space and evaluates the correlation between modelled and observed values. Cross-validation was implemented using two distinct spatial partitioning strategies: a completely randomized hexagon distribution and a spatially blocked approach. In the latter, hexagons were aggregated into uniform, continent-specific blocks within each partition to account for spatial autocorrelation. The randomized partitions were utilized exclusively for model parameter sensitivity testing, whereas the spatially blocked approach was used to validate the final models. Additionally, we compiled 107 fauna and flora inventory papers (Appendix S2), independent of the model training dataset, for external validation. However, independent validation was not feasible for all taxonomic groups, and even the validated groups presented a limited number of samples. This independent validation was assessed using the coefficient of determination (R 2) from linear regressions between the observed inventory values and the model predictions. The BioDinamica (Oliveira et al., 2019b) package, from the Dinamica-EGO (Soares-Filho et al., 2013) software, was employed to perform all the analyses (All data and models are publicly available in the Zenodo repository at https://doi.org/10.5281/zenodo.20496882. The species richness models with the highest predictive capacity were those based on DNN, using the most complex structure tested (configuration 6), which returned the lowest error under both randomized and spatially blocked cross validation (Supplementary Table S8).

To test the sensitivity of the models to variations in predictor variables, we conducted a Monte Carlo simulation with 100 iterations, in which the best-performing models were run with random perturbations applied to their predictor variables. In each iteration, the variable values were resampled cell by cell from a normal distribution centered on the original value and with a variance equivalent to 30% of their mean value, simulating uncertainties and fluctuations in input data. This approach allowed us to assess the robustness and reliability of the model predictions under variations in environmental conditions and sampling data, ensuring predictive outputs remained coherent even under controlled perturbation scenarios. To characterize how predicted richness responds to habitat condition across the full range of the index, including the values typical of converted and degraded landscapes, we computed partial dependence curves for the CVI, for each taxonomic group and for all groups. Curves were built from a regular sequence of values spanning the observed range of the index. Each CVI was applied to all map grid cells while we kept the observed values for all other species richness predictors, we then average the predicted species richness for CVI value, for each average predicted richness value related to a level of vegetation structure, as indicated by corresponded a CVI value (Supplementary Figure 2).

Validation and comparison with IUCN polygon maps

We also evaluated, for each organism group, the ability of the models to spatially predict variation in richness. To do so, we contrasted observed and predicted richness using the fit (R²) of a linear regression. This procedure allows us to assess whether the proportion of richness, rather than the absolute number, tracks the observed variation. Because sampling effort influences this metric, we restricted the analysis to hexagons in the upper quartile of the number of samples and included the sampling-effort map as a control variable in the regression.

For terrestrial vertebrate groups, in addition to this evaluation, we fitted a regression between the models’ richness predictions and richness maps derived from IUCN polygons. In this case, we did not control for sampling effort because these maps are based on distribution polygons and are therefore likely less affected by sampling effort (Marsh et al., 2023). However, because IUCN maps do not incorporate richness changes associated with habitat change, as in the model proposed here, we restricted the regression to cells located in areas with Ecosystem Integrity (Hill et al., 2022) greater than 0.9, in order to make the comparison consistent. This restriction ensures that the two products are compared where they estimate the same quantity, since potential and realized richness converge only where habitat remains intact. In degraded areas the two are expected to diverge by construction, so that disagreement there could not be interpreted as model error.

Test of spatial variation in richness

The spatial resolution of a map depends not only on cell size but also on the scale at which the phenomenon is spatially represented. Thus, some maps better describe variation at broader scales, whereas others capture variation at finer scales. Because this characteristic is relevant when richness maps are used as inputs for conservation prioritization analyses, we assessed the spatial variation of richness values using the γ (gamma) index (Chilès and Delfiner, 2012). This index quantifies, at a defined scale (distance associated with the analyzed area), the degree of spatial variation in values, with low values indicating little variation and high values indicating greater heterogeneity. The index corresponds to the empirical semivariogram estimator, whose full definition and computation parameters are given in Supplementary Equation 1. It is expressed in squared units of the mapped variable, in this case species squared, and it is not dimensionless. Because γ scales with the magnitude and the variance of the mapped variable, values are directly comparable between maps of the same variable at the same lag distance, as done throughout this study, but not across variables of different magnitude nor across lag distances without reference to the underlying variance.

We compared richness maps generated by the USSE approach, based on deep neural networks, with richness maps derived from IUCN distribution polygons. Because IUCN provides richness maps only for terrestrial vertebrates (IUCN, 2025), the comparison was restricted to these groups. We evaluated three scales, with radii of 50, 200, and 500 km, in the calculation of the γ index. To spatially represent the areas with the greatest variation, we mapped the γ index for the IUCN vertebrate richness maps and for the USSE approach maps at the 50 km scale.

To quantify the relative similarity between the spatial patterns of richness of different taxonomic groups and identify which group best represents overall richness, we estimated the similarity between the richness maps modeled by the USSE using rank-based correlation (Spearman) calculated on a random sample of cells shared by all rasters and evaluated its significance after correcting for spatial autocorrelation (adjusting the effective degrees of freedom), minimizing the pseudo-replication bias associated with spatial dependence between cells (Clifford et al., 1989). Then, we converted the similarities into a (1 − ρ) distance matrix and applied hierarchical clustering to generate a dendrogram, in which the proximity between branches indicates greater agreement of the spatial gradients of richness between groups, allowing us to infer which group is the best proxy for overall richness and the degree of similarity between the richness patterns.The BioDinamica (Oliveira et al., 2019a) package was used to run all analysis.

The models demonstrated high predictive performance across all groups and in the combined analysis, with R² values ranging from 0.6 to 0.95 (Figure 1; Supplementary Figures 3–6). Predictive performance was lowest for Trichoptera (0.6) and Gastropoda (0.6), and highest for Aves (0.95), Mammals (0.92), Reptiles (0.93), and Angiosperms (0.87) (Supplementary Figures 3, S4). The independent validation using literature inventories confirmed the robust predictive performance of the models across the evaluated groups, with R 2 values ranging from 0.55 to 0.87. Performance was lowest for Lepidoptera (0.55) and Angiosperms (0.64), and highest for Aves (0.87), Araneae (0.87), Formicidae (0.84), and Odonata (0.77). Mammals and Amphibians also exhibited strong concordance, with R 2 values of 0.71 and 0.73, respectively. Model sensitivity to variation in predictor variables was low, as indicated by Monte Carlo tests showing changes between 0.26% and 1.31% (Supplementary Figure 7). For terrestrial vertebrates, richness maps generated using the USSE approach closely matched IUCN polygon-based maps in areas of high ecological integrity, with R² values of 0.63 for Birds, 0.80 for Reptiles, 0.84 for Mammals, and 0.85 for Amphibians (Supplementary Figure 8). Partial dependence curves for the CVI were smooth and continuous for all groups, with predicted richness increasing monotonically with the index across the lower half of the gradient and converging on a stable lower bound at low index values (Supplementary Figure 2).

To assess spatial heterogeneity, we used the γ (gamma) index, which quantifies the degree of spatial variation in richness values. USSE models consistently exhibited higher γ values than IUCN maps across all groups and spatial scales (Figure 2). The most pronounced differences were observed at the fine scale (50 km), where γ values for Reptiles and Mammals were more than ten times higher in USSE models than in IUCN maps, and twice as high for Birds and Amphibians. At the intermediate scale (200 km), differences were smaller (approximately 20%) for Birds and Amphibians, but larger for Mammals and Reptiles, with USSE γ values three to four times higher than IUCN (Figure 2). At the broadest scale (500 km), USSE γ values remained higher, ranging from 42% to 267% above those of IUCN (Figure 2). Additionally, areas with high γ values were more widely distributed in USSE maps, indicating greater spatial heterogeneity in richness, except in desert regions (Figure 3). Model predictions were strongly correlated with the Composite Vegetation Index (CVI), which highlighted differences in richness between native and anthropogenic areas; for example, the models captured reduced richness in the Amazon deforestation arc, southern African savannas, and deforested regions of the Malay Archipelago (Figure 4A). In addition, the models capture a higher level of spatial granularity, revealing variations that are absent in traditional polygon-based representations (Figures 4A–G).

Across all organism groups, the models indicated higher species richness in tropical forest areas, with richness decreasing in more open formations and reaching minimum values in deserts (Figure 1; Supplementary Figures 3–6). The Amazon forest exhibited the highest richness for all groups, while the African tropical forest and the Malay Archipelago displayed group-specific variations, for instance, Amphibians had smaller areas of high richness in these regions compared to other groups, and Mammals, Blattodea, and ants also showed reduced areas of peak richness in African forests. Madagascar exhibited substantial variation in richness patterns among groups, but consistently, the highest richness was concentrated in the island’s eastern portion. Among savanna regions, the South American Cerrado had the highest richness for all groups. Deserts generally had the lowest richness, with the Atacama Desert standing out as the least rich among deserts for most groups.

Overall, richness maps were highly similar across groups, with a maximum difference in correlation distance of 0.09 (Supplementary Figure 9). Amphibians and mammals exhibited the greatest similarity to each other and were closely aligned with the bird richness pattern. In contrast, reptiles’ patterns were more similar to those of bees than to other vertebrates. Ants’ richness patterns were closer to ferns than to other invertebrates, and trichopterans clustered with angiosperms, ants, and ferns (Supplementary Figure 9). Of all groups, the richness gradient of lepidopterans showed the highest rank correlation with the aggregate pattern across all organism groups, whereas spiders and coleopterans were more similar to each other and less similar to other invertebrates (Supplementary Figure 9).

The relatively high predictive performance of the models indicates that the set of predictor variables was adequate for estimating spatial richness patterns, supporting the use of climate and vegetation structure as predictors of biodiversity. Vegetation structure, represented here by the CVI, is to some extent a direct measure of vegetation itself and, therefore, was expected to show a strong relationship with plant diversity. Furthermore, the CVI proved informative for animal groups, which is consistent with evidence that habitat structure constrains resource availability, microclimates, and niches, thereby influencing which species and how many species can persist in a given location (Boucher-Lalonde et al., 2014; Hyman et al., 2019). Climate, in turn, in addition to acting directly on physiological limits and ecological processes associated with richness (Vázquez‐Rivera and Currie, 2015; Barreto et al., 2021), also influences vegetation structure (Ferger et al., 2014; Vázquez‐Rivera and Currie, 2015), which is quantified here in a direct way, thus, the combination of climate and habitat structure tends to capture complementary components of the environmental gradient, increasing the predictive capacity of biodiversity mapping. This performance is relevant to the applied objective of this study, which is to produce spatial estimates with higher precision and accuracy for conservation use, but it does not, by itself, provide a mechanistic explanation for the causes of richness patterns. Incorporating factors absent from our predictor set could refine these estimates further, notably biotic interactions, which impose local limits on coexistence, and lagged rather than contemporaneous responses to abiotic conditions, since productivity in one period may propagate to consumer richness only after an accumulation interval (Wang et al., 2018; Yates et al., 2018; Bouchet et al., 2019). The performance obtained here nonetheless indicates that climate, topography and vegetation structure, once sampling intensity is accounted for, already capture most of the spatial signal in richness across groups and scales, so that such additions would be expected to refine rather than to reshape the patterns reported. Sampling bias is only one of the data properties that condition richness estimates, which also depend on sample size and on detectability, both of which are most limiting for rare and cryptic species and therefore constrain how far standardizing sampling effort can improve predictive performance (Wang and Jackson, 2023). Even so, these results provide an empirical basis to guide future investigations into the determinants of richness variation in the tropics, a topic with extensive literature (Hawkins et al., 2003; Willig et al., 2003) but that still has gaps, especially when the aim is to reconcile different taxonomic groups and multiple spatial scales.

Polygon-based distribution maps are still widely used for their global coverage and standardization (Palminteri et al., 2011; Rangel et al., 2018; Anderson, 2023), despite ongoing criticism of their suitability for high-spatial resolution analyses (Mainali et al., 2020; Broekman et al., 2022; Marsh et al., 2022). These criticisms focus on commission errors and the use of broad boundaries, which can overestimate species presence when maps are stacked for community analyses, thereby affecting both richness estimates and conservation priority setting (Mainali et al., 2020; Broekman et al., 2022; Marsh et al., 2022). Such generalizations can distort assessments of representativeness and lead to overestimated protection levels, particularly at finer spatial resolutions. In this context, occurrence-based alternatives that incorporate habitat condition offer a basis for conservation prioritization that reflects current habitat conditions.

Vertebrates are commonly used as proxies for total biodiversity (Loyola et al., 2007; Monroy‐Gamboa et al., 2019; McRae et al., 2025), but reviews and assessments show that the effectiveness of such surrogates varies with geographic scale, metric, and taxonomic set, and is often limited when representing aggregated richness patterns (Forest, 2017; Oberprieler et al., 2019). This is particularly important because many studies aim to infer general biodiversity patterns to guide conservation priorities, not just to map vertebrates (McRae et al., 2025). Evidence from insects suggests that vertebrate-based priorities may overlook areas critical for the diversity of hyper diverse invertebrate groups (Diniz-Filho et al., 2010; Jenkins et al., 2013). Our results indicate that, among the groups analysed, invertebrates and particularly lepidopterans exhibit richness gradients more closely congruent with the aggregate multi taxon pattern than any vertebrate group does. Congruence of gradients is not by itself sufficient to establish a taxon as a surrogate for conservation planning, which depends on complementarity in the selection of priority areas, and we therefore present this result as an indication that the surrogacy potential of hyperdiverse invertebrate groups warrants formal evaluation. This finding underscores the importance of broadening the taxonomic scope of global products and increasing investment in surveys, taxonomy, and digitization of invertebrate collections to reduce knowledge gaps and improve macroecological inference. Recent global syntheses for bees and ants help contextualize these patterns and reinforce the need to interpret invertebrate richness maps in light of uneven sampling (Orr et al., 2021; Kass et al., 2022). For bees, our predictions are generally consistent with the regional variation described for tropical areas, but indicate proportionally higher richness in parts of the Malay Archipelago (Orr et al., 2021). This difference may reflect the strong sampling gaps reported for bees in tropical and developing regions, especially in Asia, Africa and Latin America, where available records may still underestimate regional richness (Orr et al., 2021). Thus, the Malay Archipelago may represent an area where environmental and habitat conditions support higher bee richness than is currently documented, although historical isolation and dispersal barriers may also influence realized richness. For ants, the agreement is more direct, as our models recover high-richness areas in tropical Africa and Southeast Asia, regions also identified as important centers or potential centers of ant diversity under increased sampling scenarios (Kass et al., 2022). Together, these comparisons suggest that divergences between observed and predicted richness for invertebrates may reflect not only model uncertainty, but also persistent sampling gaps and hidden diversity in poorly sampled tropical regions.

Previous studies using simulated data have shown that USSE can recover spatial patterns under various sampling scenarios by disentangling the biological signal from sampling effort (Oliveira et al., 2024). Here, we extend this evidence to a broad empirical context, demonstrating that the approach produces estimates compatible with established products for vertebrates and provides greater spatial detail, an essential feature for fine-scale spatial analyses in conservation. Moreover, the USSE approach is not limited to richness and can be generalized to other occurrence-derived biodiversity metrics, such as endemism, phylogenetic diversity, and phylogenetic endemism. The availability of these metrics within a unified framework can support prioritization analyses, impact assessments, and macroecological studies.

Therefore, the combination of the USSE approach with deep neural networks provides a methodological alternative to polygon-based maps by generating spatially explicit estimates directly derived from occurrence data, showing similarities with established products for vertebrates from the IUCN, nevertheless capturing spatial heterogeneity with greater granularity, and incorporating changes in richness due to current habitat alterations. This methodological advancement is particularly relevant for conservation applications, as prioritization decisions are often made at local and regional scales (Wiens and Bachelet, 2010; Ceauşu et al., 2015), where increased spatial detail ensures that estimates more accurately reflect the current state of biodiversity, hence enhancing the ability to identify areas at higher risk of richness loss.