The palm E. edulis has a wide geographic distribution, occurring in Brazil, Argentina and Paraguay (Henderson et al., 1995). Most of the native distribution range is in the Brazilian Atlantic Forest, where the species has a wide latitudinal distribution, but the species also occurs naturally in the Brazilian Cerrado (Henderson et al., 1995). Euterpe edulis is considered the most important non-timber forest product of the Atlantic Forest, as its palm heart has been extensively exploited for human consumption (Silva-Matos et al., 1999). This economic activity eliminates adults and reduces seed production, thus decreasing population density and growth (Freckleton et al., 2003; Souza and Prevedello, 2020). Also, this species has important ecological functions, as its fruits are considered a keystone resource consumed by a large variety of animals, especially birds and mammals (Galetti et al., 1999).

Density estimates of E. edulis were obtained from an extensive literature search on different databases (Web of Science, SciELO, Scopus, and Google Scholar), and from Melito et al. (2014) (see Souza and Prevedello, 2020 for details). We only included density estimates when precise coordinates were provided or could be obtained based on sampling sites description. We retained information from studies that measured the density (individuals/ha) separately for different ontogenetic stages of E. edulis. We did not use density estimates from studies that only sampled E. edulis individuals above a minimum size (e.g., diameter at breast height >10 cm). For studies reporting density values for different years in the same sampling site, we averaged density across the whole period, thus including a single density estimate per site in our analyses. Similarly, for studies with nearby sampling plots (<1000 m) in the same site, we calculated the average density across plots. Some studies showed density estimates for sites spread across a broader region, which were considered as independent estimates (e.g., Oliveira et al., 2014). With these procedures, we obtained density data for a total of 50 sites spread across most of the Atlantic Forest (Fig. S1a; see also Souza and Prevedello, 2020 for further details). The large geographical area covered by these population data includes large variation in latitude and climatic conditions, allowing a robust test of how variation in environmental suitability affects abundance.

From each study, we estimated E. edulis density for three ontogenetic stages separately: seedlings, adults and intermediate-class individuals. The seedling stage included stemless individuals with palmate leaves. The adult stage included reproductive individuals, i.e. individuals with flowers and/or fruits (Souza et al., 2018). The intermediate stage included sapling and immature individuals; finer subdivisions of this stage were not possible because different studies used different classification criteria. We also measured per capita seedling density (Comita et al., 2007) in each population, as the ratio between seedling density and adult density.

A previous study has shown the occurrence of density-dependence effects on E. edulis, with a lower probability of seedling growth and survival in higher densities of seedlings and adults, which in turn can promote the regulation of population growth (Silva-Matos et al., 1999). Therefore, we considered “density class” as a predictor variable in our analysis, to check a possible density-dependent effect on the abundance-suitability relationship. To do so, we classified each of the 50 E. edulis populations as either “low-” versus “high-density” populations, based on the median density observed across all populations (see “Data analysis”, below).

We also recorded the protection status of each population as a second predictor variable. We included protection status as a proxy for anthropogenic disturbance, in particular harvest occurrence, since E. edulis populations inside protected areas are less subjected to palm-heart exploitation (Souza and Prevedello, 2020). In addition, protection status may also be a proxy for forest cover and seed disperser diversity, since both factors tend to be higher inside protected areas than outside. In addition, a previous study showed that E. edulis adult density is higher inside protected areas (Souza and Prevedello, 2020). This study, however, has not tested for differences in the other ontogenetic stages and per capita seedling density, neither investigated abundance-suitability relationships.

We obtained occurrence records of E. edulis in South America from the literature (Henderson et al., 1995; Sales et al., 2021) and from online databases: GBIF (https://www.gbif.org/), REFLORA (http://reflora.jbrj.gov.br), Project JABOT (http://jabot.jbrj.gov.br), and Steere Herbarium (http://sweetgum.nybg.org/science). For modelling, we only kept occurrence records spaced at least 10 km apart (ten times the pixel size of the environmental layers – 1 km; see below) to maximize their independence. With these procedures, we obtained a total of 335 valid occurrence records for E. edulis (Fig. S1b).

To generate the background area over which the ecological niche models (ENMs) were built, we created a 200 km buffer around the minimum convex polygon encompassing all occurrence records of E. edulis (as in Souza and Prevedello, 2021). We obtained 19 bioclimatic variables and elevation from WorldClim 2.1 (Fick and Hijmans, 2017) at 30 s resolution (i.e., a cell with ∼1 × 1 km), representing historical (∼1970–2000) conditions. To select the bioclimatic variables for modelling, we extracted their values at 1,000 random points within the background area and portrayed them into a correlation matrix. We only kept six bioclimatic variables with low correlation (r < 0.65), in addition to elevation: temperature annual range, mean temperature of wettest quarter, annual precipitation, precipitation seasonality, precipitation of warmest quarter, and precipitation of coldest quarter.

We also considered the Normalized Difference Vegetation Index (NDVI) as an additional explanatory variable in the ENMs. This variable may be a proxy of biotic conditions and habitat suitability for E. edulis, as individuals of this species only grow and survive inside forest areas (Gatti et al., 2014), and their seeds are dispersed by many forest-dependent species (Sales et al., 2021). NDVI was obtained by averaging monthly data (2001–2010) from the MODIS-Terra MOD13C2 v006 dataset with resolution of 0.05 min (Didan, 2015). After averaging, we resampled the original resolution to 30 s through a bilinear interpolation. The ENMs including NDVI returned similar results as the models without NDVI, but with slightly lower coefficients of determination for explanatory models (Table S1). Therefore, we present in the main text only the results for the ENMs without NDVI.

We estimated environmental suitability for E. edulis using five algorithms that encompass different modelling assumptions and approaches: BIOCLIM (Booth et al., 2014), GLM (Guisan et al., 2002), MaxEnt (Phillips and Dudík, 2008), Random Forest (Breiman, 2001), and SVM (Drake et al., 2006). For BIOCLIM, we used as training data only the occurrence records of E. edulis. For GLM, SVM and Random Forest, we used as training data the occurrence records plus 10 pseudo-absence points for each record, totaling 3,350 pseudo-absences. Pseudo-absence points were randomly chosen from pixels located within the background area, but outside the environmentally suitable area as estimated from BIOCLIM (following Lobo and Tognelli, 2011). For MaxEnt, we used as training data the occurrence records plus 10,000 background points randomly chosen within the background area.

For each algorithm, we used 10-fold cross-validation, splitting data into 90% training and 10% test, replicated 10 times (Elith and Leathwick, 2009). For testing, we used the 10% remaining presences and pseudo-absences for all algorithms. We estimated the predictive performance of each model through the true skill statistic (TSS, Allouche et al., 2006) and the area under the receiver operating characteristic curve (AUC), separately for each algorithm and iteration. We only retained as valid models those with both TSS > 0.70 and AUC > 0.85 (as in Souza and Prevedello, 2021). We then used each valid model to produce a separate map of predicted environmental suitability across the Atlantic Forest.

To ensemble the continuous maps produced by each algorithm, we first rescaled pixel values to vary between 0 and 1 for standardization, as the original output scale varied among algorithms (i.e. minimum = 0 and maximum = 1 for BIOCLIM, MaxEnt and Random Forest, but not for GLM and SVM). We then calculated a TSS-weighted mean across the maps, obtaining a single final map of continuous suitability. Since the choice of the algorithm is the greatest source of uncertainty in ecological niche modelling (Watling et al., 2015), we used in the statistical analysis only the suitability values from the ensemble map. For each site with E. edulis density data, we extracted the environmental suitability value from the final ensemble map (Fig. S1b). Algorithms were implemented in R using the packages ‘dismo’ (Hijmans et al., 2017), ‘randomForest’ (Liaw and Wiener, 2002) and ‘raster’ (Hijmans, 2017).

We first tested the overall relationship between population density and environmental suitability across the 50 sites combined, using a linear regression model. We applied a logarithmic transformation to density values (i.e. ln(density+1), to normalize model residuals and thus attend the assumption of linear modelling. We focus the main text on the results for total population density only, i.e. the density of the three ontogenetic stages combined, as the results were similar for the three stages separately (Figs. S2 and S3). We used data from the 50 sites, including one site with zero density, sampled in a landscape where other populations of the species were recorded (Pires, 2006). A separate linear model was used to test the overall relationship between ln(per capita seedling density + 1) and environmental suitability. This model was fit using only 43 populations, as five populations had zero density for adults, and information on seedling density was unavailable for two populations, precluding calculation of per capita seedling density. For all analyses, we assumed that the 50 sites provided independent information, based on the analyses of Souza and Prevedello (2020), which showed that pseudo-replication is unlikely in this dataset, as similar results were obtained using simple or mixed-effect linear models.

To test whether the relationship between density and suitability differed between low- versus high-density populations, we first classified each of the 50 populations according to its total density. “Low-density” populations corresponded to populations with total density lower or equal to the median total density observed across all populations (1,844 individuals/ha), whereas “high-density” populations corresponded to populations with density higher than the median density. We used the median as the cutting value to obtain similar sample sizes for low- and high-density populations. We then built two separate linear models for each dependent variable (population density or per capita seedling density, both as ln(variable + 1)), one model for each density class (i.e. low- or high-density populations), including environmental suitability as the explanatory variable. Similar results were obtained using an alternative analytical approach, in which we included environmental suitability, density class (low- versus high-density) and their interaction as explanatory factors in a single model for each dependent variable (Table S2a,b). This alternative approach showed significant interaction terms, which indicate that the effects of environmental suitability on the dependent variables differ between low-versus high-density populations, the same result shown by the separate linear models (see Results).

Similarly, to test whether the relationship between density and suitability differed between protected and unprotected populations, we built two separate linear models for each dependent variable (population density or per capita seedling density, both as log(variable + 1)), one model for protected and another for unprotected populations, with environmental suitability as the explanatory variable. Similar results were also obtained when considering environmental suitability, protection status and their interaction as explanatory factors in a single model for each dependent variable (Table S2c,d).

Population density increased significantly with environmental suitability across the 50 sites (df = 48, t = 4.82, P <  0.001; R2 = 0.32; Fig. 1a). Similarly, per capita seedling density increased significantly with environmental suitability (df = 46, t = 3.50, P =  0.001; R2 = 0.21; Fig. 1b).

Fig. 1.

Overall relationships between environmental suitability and (a) population density and (b) per capita seedling density of Euterpe edulis across 50 populations in the Brazilian Atlantic Forest. Density values are represented as ln(individuals/ha).

(0.16MB).

Environmental suitability had different effects on population density in low- versus high-density populations. Environmental suitability had no significant effect on population density in high-density populations (df = 23, t = 0.31, P =  0.76), but a positive effect in low-density populations (df = 23, t = 4.29, P <  0.001; Fig. 2a). Similar patterns were observed for the three E. edulis ontogenetic stages analyzed separately (Fig. S2).

Fig. 2.

Effects of environmental suitability on (a) population density and (b) per capita seedling density in low-versus high-density populations of Euterpe edulis in the Brazilian Atlantic Forest. Low-density populations correspond to populations with density lower or equal to the median density observed across all populations; high-density populations had density higher than the median density. Density values are represented as ln(individuals/ha). Statistically significant relationships (P ≤  0.05) are highlighted in bold.

(0.41MB).

Similarly, environmental suitability effects on per capita seedling density were different in low- versus high-density populations. Environmental suitability had no significant effect on per capita seedling density in high-density populations (df = 21, t = 0.09, P =  0.93), but a positive effect in low-density populations (df = 23, t = 3.02, P =  0.006; Fig. 2b).

The inclusion of protection status (inside versus outside protected areas) increased the explanatory power of models for population density (R2 = 0.46) when compared to the model with environmental suitability only (R2 = 0.32; see above). Population density was on average 1.7 times higher inside (mean ± SD: 7,362 ± 13,885) than outside protected areas (4,312 ± 6,689 individuals/ha). Environmental suitability had different effects on population density inside and outside protected areas. Environmental suitability had no significant effect on population density inside protected areas (df = 29, t = 0.20, P =  0.84), but a positive effect outside protected areas (df = 17, t = 4.71, P <  0.001; Fig. 3a). Again, similar patterns were observed when different E. edulis ontogenetic stages were analyzed separately (Fig. S3).

Fig. 3.

Effects of environmental suitability on (a) population density and (b) per capita seedling density in populations of Euterpe edulis located inside and outside protected areas in the Brazilian Atlantic Forest. Density values are represented as ln(individuals/ha). Statistically significant relationships (P ≤  0.05) are highlighted in bold.

(0.39MB).

Finally, environmental suitability had positive effects on per capita seedling density both inside (df = 28, t = 2.11, P =  0.04) and outside protected areas (df = 16, t = 2.30, P =  0.04; Fig. 3b).