Climate defined but not soil-restricted: the distribution of a Neotropical tree through space and time

Brosimum rubescens, a tree species with Neotropical distribution, can achieve local monodominance in Southern Amazonian forests. Understanding how and why this species varies across space and time is important because the monodominance of some species alters ecosystem complexity. Here we evaluated the fundamental ecological niche of B. rubescens by species distribution models (SDM), combining predictive environmental variables with occurrence points, and determined the temporal persistence and how the spatial distribution patterns of this species vary with different environmental predictive variables. To generate the SDMs, we incorporated predictive environmental variables as main components of climatic, hydric and edaphic variables. All algorithms showed higher performance in spatial predictions for hydric variables and for the combination of climatic, hydric and edaphic variables. We identified that the potential niches of B. rubescens seem to be defined by climatic fluctuations, with the edaphic conditions not limiting the presence of this species in the evaluated spatial scale. From the last glacial maximum to the present, this species seems to have increased its spatial amplitude; however, from the present to the future, predictions suggest that B. rubescens will experience a considerable loss of its range. Our findings showed independent and combined effects of different environmental variables, allowing us to identify which are limiting or facilitating the spatial distribution of B. rubescens. We corroborate the spatial persistence and geographical fidelity of the species’ distribution patterns over time.


Introduction
The global distribution of biodiversity follows a general pattern of spatial radiation from tropics to polar latitudes, and tropical forests are no exception to this pattern (Wiens and Donoghue 2004;Murphy and Bowman 2012). Within this spatial radiation, some tropical tree species dominate 50 to 100% of the total biomass forming the so-called monodominant forests (Hart et al. 1989). The monodominance of some species creates complex ecosystems, with no single or independent physical, chemical, or biological mechanism to explain the dominance over space and time (Torti et al. 2001;Peh et al. 2011;Marimon-Junior et al. 2019;ter Steege et al. 2019). The mechanisms promoting and maintaining the monodominant forests remain poorly understood and do not fit the spacio-temporal patterns of global biodiversity (Wiens and Donoghue 2004). One study suggested that these unique communities are associated with hydric, edaphic, and ecological traits occurring together (i.e. shade tolerance, slow decomposition rate, mast seeding, large seeds and poor dispersal, herbivore avoidance, ectomycorrhizal associations) under low exogenous disturbance over long time periods (Peh et al. 2011).
Neo-and paleo-tropical monodominant forests share some common edaphic characteristics, such as poor soil nutrient availability or an excess of certain elements (Marimon et al. 2001;Marimon-Junior 2007;Peh et al. 2011;Elias et al. 2018Elias et al. , 2019Marimon-Junior et al. 2019). Nevertheless, there is no clear evidence that only soil properties are the main determinants creating or maintaining the monodominance (Brito-Rocha et al. 2016;Elias et al. 2018;Marimon-Junior et al. 2019). Modeling the spatial distribution of monodominant species, using different environmental variables and different temporal projections, allows a better understanding of their combined and independent spacio-temporal effects and a better approximation to their Grinnellian niche (Soberón 2010;Hijmans et al. 2015;Alvarez et al. 2020).
Brosimum rubescens Taub. (Moraceae), the target species of the present study has, a Neotropical distribution and is associated with non-seasonally flooded tropical forests, with the dispersion center in the Amazon Basin and Guyana, with some records in Central America and the Caribbean (Berg 1972;Palacios 2005). In the Amazonia-Cerrado transition, B. rubescens can form monodominant forests of up to 5000 ha and can represent up to 85% of all trees >10 cm diameter, one of the highest values recorded for a Neotropical monodominant species (Marimon et al. 2001). In addition, the high density of B. rubescens affects species diversity, unbalancing the uniformity of their abundance by acting as a biological barrier (Marimon-Junior 2007;Marimon et al. 2008;Marimon-Junior et al. 2019).
Although the occurrence and monodominance of some Neotropical forests are not conditioned by restrictive climatic conditions, differences in the distribution of the canopy understory light and speciesspecific capacities for the use of this resource is one of the factors of most significant repercussion for structuring such forests (Marimon et al. 2008;Peh et al. 2011). B. rubescens show development closely related to the amount and quality of sun radiation received by the young individuals (specially seedlings), and to the morphological and anatomical characteristics that provide phenotypic plasticity and periodical growth (Marimon et al. 2001(Marimon et al. , 2008(Marimon et al. , 2012. Also, soils affect the demography of B. rubescens in different phases of development, and the soil where this species is found has been reported to be acid and nutrient poor, dystrophic (Quesada et al. 2010) and with Ca/Mg ration unbalanced (Marimon-Junior et al. 2019). Some studies have shown that B. rubescens invests in greater root biomass under more solar radiation as a response to seasonal hydric deficits (Butler et al. 2013), but the response of the species to climatic environmental dynamics is unknown (Marimon et al. 2008(Marimon et al. , 2012. The temporal and spatial maintenance of B. rubescens monodominant forests is linked to episodic recruitment in response to extreme climatic interference, probably associated with the impacts of prolonged drought disturbances related to El Niño Southern Oscillation phenomenon (ENSO) (Marimon et al. 2014a(Marimon et al. , 2020Marimon-Junior et al. 2019). Despite the morpho-anatomical and physiological attributes and environmental adaptations of B. rubescens, these drought events have caused a drastic reduction in the natural regeneration of the species, which may contribute to its population decline in the Amazonia-Cerrado transition (Marimon et al. 2012(Marimon et al. , 2020. Therefore, future climatic variations may affect the sensitive balance between mortality and recruitment, affecting the local and regional distribution of B. rubescens. Despite several studies about the monodominant forest of B. rubescens, the patterns of the spatial distribution of this species, using different sets of variables and its spatial trends in future climate change scenarios, are not yet understood. In this work, we aimed to: 1) determine how the spatial distribution patterns of B. rubescens vary with different environmental predictive variables, and 2) evaluate the temporal and spatial persistence of B. rubescens in the Neotropics. Knowing the spatial and temporal patterns of this species would facilitate the planning of management strategies for conservation of this species and the unique monodominant tropical forest ecosystems.

Study area
The results of the models were restricted to the northern portion of the South American continent ( Fig. 1 and Supplementary Material S1), but to generate the models, the Neotropical region was chosen (Leroy et al. 2018). B. rubescens is a canopy species belonging to Moraceae family (Berg 1972) that can reach 90 cm diameter at breast height with heights of up to 35 m (Laurance et al. 2004). In Brazil, the occurrence of this species is reported in forests in the Amazonia, Fig. 1 Geographic distribution of spatial records and Kernel density estimation (KDE) of Brosimum rubescens B. rubescens in tropical South America, with higher values of density (red) associated with zones with higher frequency and intensity of collection Cerrado, and Atlantic Forest biomes, with spatial interruptions coinciding with the Caatinga Biome and monodominant forests registered only in Amazonia-Cerrado transition zone (Marimon et al. 2001(Marimon et al. , 2014a(Marimon et al. , 2020. In Amazonia and Atlantic Forest biomes the soils are dystrophic but shallow for the former and humid and deep for the latter (Quesada et al. 2010;Joly et al. 2012). Cerrado soils are highly weathered and deep, acidic, dystrophic, and with high levels of iron and aluminum (Ribeiro and Walter 2008). The predominant climate types in Cerrado, Amazonia and Atlantic Forest biomes, especially in the Brazilian areas where B. rubescens was registered are Af, Am and Aw (Dubreuil et al. 2018).

OCCs: Obtaining and preprocessing
We obtained the occurrence data (OCCs: Fig. 2a) for B. rubescens from the virtual platforms GBIF (http:// www. gbif. org) and speciesLink (http:// splink. cria. org. br), complemented with data from the Laboratório de Ecologia Vegetal (Universidade do Estado de Mato Grosso, Nova Xavantina). We incorporate all OCCs data regardless of the sampling date. The initial database had 256 records that we filtered spatially to identify and eliminate data outside the Neotropics. Subsequently, we eliminated all repeated, null, or dubious taxonomic data. A common problem in data obtained from virtual platforms is the spatial biases associated with the frequency, intensity and heterogeneity of the sampling efforts (Pearson et al. 2007;de Oliveira et al. 2014;De Marco and Nóbrega 2018). Using Kernel density estimation (KDE) spatial heterogeneity (Silverman 2018), was observed in the OCCs of B. rubescens (Fig. 1). To obtain the KDE, in the QGIS software (version 3.14; QGIS.org 2020), we used the interpolation extension `Heatmap´ with Epanechnikov Kernel function and 10,000 km as bandwidth value. We evaluated the spatial autocorrelation (SAC) of the OCCs using Moran index (Fig. 2a). As the data presented SAC (values greater than |0.2|) and spatial heterogeneity, we applied the "spThin" package with 10 km as the cutoff value (Aiello-Lammens et al. 2015) from the R platform (R Core Team 2020).

Preparation of environmental predictors
To better approximate the fundamental niche of B. rubescens we used predictive environmental variables (Supplementary Material S2) that show an effect at different scales and magnitudes (Soberón 2010;Domisch et al. 2015). Environmental variables act with spatial connections that follow a spatial hierarchical dependence pattern, with climate acting on the global scale, geomorphology on the regional scale and hydrological variables on the local scale (Wiens and Donoghue 2004;Domisch et al. 2015). In this sense, despite the system emergence properties, evaluating each set of variables in blocks (climate, edaphic and hydric) may promote a better understanding of the effect of each set on the spatial distribution patterns of different species (Velazco et al. 2017;Alvarez et al. 2020). We established three subsets of environmental predictors: climatic, edaphic, and hydric variables. In the climatic subset, we used 19 predictor variables obtained from the WorldClim platform (http:// www. world clim. org, Hijmans et al. 2005). For the edaphic characterization, we used 17 variables (Supplementary Material S2 and Fig. 1) obtained from the SoilGrids project (http:// www. soild grids. org).
As described in different studies, life cycles, biomass dynamics, growth, stature, and species distribution are conditioned by the seasonality of climatic variables related to hydric stress (Bruijnzeel 2004;Feldpausch et al. 2011Feldpausch et al. , 2016. In this sense, we incorporated 23 variables in the hydric subset (Supplementary Material S2 and Fig. 1) obtained from different platforms: six variables from the potential evapotranspiration (http:// www. cgiar-csi. org), one from actual evapotranspiration and three from soil water stress (http:// www. cgiar-csi. org/ data), being six variables of relative humidity (https:// www. climo nd. org), three of cloud cover (http:// www. cgiar-csi. org/ data), three of water vapor pressure (http:// world clim. org/ versi on2) and one topographic wetness index variable (http:// envir em. github. io). Considering the relative temporal homogeneity and importance for the characterization and physical description of the landscape of the five Fig. 2 Methodological procedure to obtain SDMs of Brosimum rubescens B. rubescens using Climatic (C), Edaphic (E), Hydric (H), and all combined variables (CEH-1) in the present, Last Glacial Maximum (LGM), middle Holocene (MH) and future (F) scenarios. The main points are OCCs: Obtaining and pre-processing (a), Preparation of environmental predictors (b), SDM methods applied (c), and SDM configuration and evaluation (d) following variables: global relief model (http:// geoda ta. grid. unep. ch), elevation, slope-aspect (http:// srtm. csi. cgiar. org) and terrain roughness index (http:// envir em. github. io), we included them to each one of the three subsets described above.
To generate the species distribution models (SDMs) of past and future projections, we selected three scenarios; last glacial maximum (LGM: 22,000 years ago), middle Holocene (MH: 6000 years ago), and future (F: 2041-2060). We also used three models of global circulation for the selected periods: CCSM4, MPI-ESM-P, and MIROC-ESM, available in the WorldClim platform and selected because they generate robust temporal projections (Hijmans et al. 2005(Hijmans et al. , 2015Carnicer et al. 2020). Due to the increase in greenhouse gas emissions and trends expected for the future (Riahi et al. 2011; IPCC 2013), we considered the most extreme emission scenario [representative concentration paths (RCPs) = 8.5]. Due to the global climate context and intensified anthropic actions, we selected RCP 8.5 which, according to the Intergovernmental Panel on Climate Change in 2004, represents the most pessimistic and probable scenario for projections until the year 2100 (Schwalm et al. 2020).
We incorporated all the predictor variables with a resolution of 5 arc minutes (~10 km in Ecuador), and we delimited them geographically using a Neotropical cut mask. We standardized all environmental predictor variables (μ = 0; SD = |1|) to avoid different weights (Fig. 2). Subsequently, we performed a principal component analysis (PCA) to reduce environmental collinearity (De Marco and Nóbrega 2018): PCA1 (C: 24 climatic variables), PCA2 (E: 22 edaphic variables), PCA3 (H: 28 hydric variables), PCA4 (CEH-1: all variables together, 69). Before performing the PCAs, we projected the linear coefficients obtained for the current scenario in each past and future climate scenario. In each climate scenario, we used the first eight axes of the PCA that explained more than 95% of the original environmental variation of the original predictor variables as the new predictors for species distribution. We also obtained the independent contributions of each variable considered in the axes of all PCA (Supplementary Material S2 and Fig. 1).

SDM methods configuration and evaluation
We applied the SDMs that present greater predictive performances and flexibility with presence-only data (Fig. 2c), using the following algorithms: bioclim (BIO), Gower (GWR), maxent (ME), generalized linear models (GLM), random forest (RF) and support vector machine (SVM). To meet the requirements of the algorithms, we generated 10,000 background points randomly distributed within the study area, without spatially coinciding with known occurrence points. Then, we divided the data into 70% for test and 30% for training and validation. We generated the first models for the present scenario (PCA1, PCA2, PCA3, and PCA4) and later for the LGM, MH, and future scenarios. To obtain the potential distribution of B. rubescens, we used the most straightforward configurations of the models (Pearson et al. 2007;Varela et al. 2011), leaving the results to depend on the interactions between the algorithms and the input data (Fig. 2d).
We configured BIO and GWR algorithms by default using the "bioclim" (Hijmans and Graham 2006) and "gower" (Gower 1971) packages, respectively, and ME algorithm ("maxnet" R package: Phillips and Dudík 2008) in a predetermined way using linear characteristics. Where the mean and variance of the variables remained close to the observed values (Elith et al. 2011), the output was logistic type, with 1000 iterations and 10,000 background points. For GLM ("glmnet" R package: Hastie et al. 2016), the response data assume values of zero (absence) and one (presence), we used a logistic link function, considering all single predictor variables and without interaction between them. We optimized the estimation of the smoothing parameter using the proximal Newton method without selection of covariates backwards or forwards (Velazco et al. 2017). We adjusted the RF algorithm ("randomForest" R package: Liaw and Wiener 2002) automatically, used 500 trees in the adjustment step, and considered only those models with a minimum out-of-bag error, maintaining standard values in the means of the variables (Hastie et al. 2016;Velazco et al. 2017). Finally, we adjusted the SVM ("libsvm" R package: Chang and Lin 2011) algorithm in a predetermined manner with a linear kernel function type c-SVC with probabilistic output (Velazco et al. 2017).
We executed and built all SDMs using the R "dismo" package (Hijmans et al. 2015) with specific functions of each algorithm. With some exceptions, we configured the SVM models from the "ksvm" function of the "kernlab" package (Karatzoglou et al. 2004) and RF from the "randomForest" package (Liaw and Wiener 2002). Moreover, we defined the randomization process to 30 iterations for each algorithm (using a self-initialization analysis) to decrease the spatial structure between the training/test data (subsets of the same OCC points), to reduce possible correlation errors, and to reduce biased evaluations. Therefore, we generated 180 models for each set of predictive variables, totaling 1260 models.
The performance of the SDMs generated was evaluated by two families of statistical metrics: dependent and independent thresholds (Liu et al. 2005). In this sense, we applied the true skill statistic estimator (TSS) in which the response values vary between −1 to 1 where SDMs >0.5 are considered acceptable or adequate (Allouche et al. 2006). To complement, we applied the independent threshold area under the receiver-operator curve (AUC: Fielding and Bell 1997), which discriminates areas of omission from areas with known occurrence. Values range from 0 to 1, where the closer to 1 the better is the model's ability to represent reality (Elith 2000). To verify the differences in AUC and TSS values among the PCAs and the applied algorithms, we performed Kruskal-Wallis non-parametric tests since the data did not meet the assumption of normality.
We used as a threshold the lowest adequacy value within the incorporated OCCs [(lowest presence training (LPT)], allowing the generation of presence/ absence matrices (binary models) from the adequacy of continuous matrices projected for South America (environmental models). Applying the consensus of the weighted average we assembled all the suitability matrices obtained for the different temporal scenarios. In this way, to generate a final distribution matrix for B. rubescens in the temporal scenarios, we applied a consensus method, assembling all the predictions considering only those models with TSS values greater than the mean TSS value.

Model evaluation
All algorithms showed high performance and we registered statistical differences among them (Fig. 3). Moreover, all models showed TSS and AUC adequate values (Table 1), but the best performances were associated with the SVM, ME, and RF algorithms, with TSS values >0.5 and AUC > 0.8 (Table 1, Fig. 3). The models generated exclusively from hydric variables and combining all the variables performed better, regardless of the algorithms (Table 1). For both PCAs and algorithms, the null hypothesis was rejected (p value <0.001) (Fig. 3).

Present projections
The spatial complexity of the generated predictions was different depending on the built-in predictor variables. The suitability matrices generated from PCA1 (C: Climatic variables) and PCA3 (H: Hydric variables) showed wide continuity and spatial aggregation (Fig. 4). The matrices generated from the PCA3 (E: Edaphic variables) showed diffuse spatial distribution patterns of wide spatial amplitude (Fig. 4). When considering all the variables together (CEH-1), two situations can be observed: high connectivity in spatial distribution patterns and/or low spatial complexity concerning the other PCAs used. Regardless of the variables incorporated, all models showed centers of greater suitability in the central region of the Amazon Basin, the Atlantic coast, and the North region of South America (Fig. 4). As a particular case, CEH-1 showed predictions of suitability associated with the Central America and the Caribbean islands, although we have not incorporated occurrence records for these areas. Regardless of the set of built-in predictor variables, the predictions of environmental suitability were superimposed on the occurrence records obtained in the field.
When we combined all the variables to assess the independent contribution of each within a 200,000 km 2 hexagon, we observed in more detail the environmental structure (C: Climatic; E: Edaphic and H: Hydric) of the study area (CEH-2: Fig. 4 and Supplementary Material S3). In this sense, in the Amazon Basin and the Atlantic Coast, the balance between the three sets of variables makes it possible to recognize, geographically, the higher probabilities of environmental suitability and they coincide spatially with the reports of occurrence of B. rubescens.

Projections of temporal climate scenarios
The spatial distribution patterns of the models for different scenarios showed some overlap in the historical spatial distribution patterns of B. rubescens (Fig. 5).  Despite the "temporal distance" of the compared scenarios, the environmental suitability matrices (Fig. 5: Species model distribution) show high spatial overlap, with good connectivity, being more significant and evident in the environmental suitability centers associated with the Amazon Basin, North of South America and Atlantic coast. Although the matrices obtained from the projections have diffuse regions, both for the past (LGM and MH) and future (F) scenarios, the highest suitability values are superimposed on the present spatial distribution patterns (CEH-1). The projections were made using only climatic variables and not considered them as predictors of the landscape, such as topographic or edaphic variables.
While we did not incorporate occurrence records for Central America or the Caribbean islands, projections for both past and future scenarios showed environmental suitability associated with these regions (Fig. 5: Species model distribution). The suitability matrices for all temporal projection scenarios also show values associated with the Patagonian region, with wide spatial amplitude of predictions in the LGM and future scenarios. In the matrices of all projections two large areas can be seen, one between the Caribbean islands and Amazonia and the other between Amazonia and Atlantic Forest biomes (Fig. 5: Species model distribution). The pattern of spatial overlap between different scenarios can be highlighted in higher resolution in the binary matrices (Fig. 5: Binary models). The present, past and future (PPF) map is the result of overlapping all of them in a single matrix (Fig. 5: Binary models). The values associated with each binary matrix represent the area (km 2 ) occupied by B. rubescens within the study area. The central triangle shows losses and gains areas from one scenario to another and a gain of 932,659 km 2 can be noticed from the LGM to the CEH-1 (present). Also, it shows a substantial loss of 568,118 km 2 from the present to the future projection. However, this loss is not total: when comparing the surface calculated for the LGM with the amount of surface that was projected for the future, the result remains positive (Fig. 5: Binary models). Despite these variations in the total area values over time, in the PPF matrix, B. rubescens has a high and recurring presence associated with the coastal region of the South American Caribbean, Amazon Basin, and the Atlantic Forest, where we found highest adequacy values for all projections (Fig. 5: Binary models).

Discussion
The spatial distribution patterns predicted for B. rubescens, regardless of the variables, allowed a great spatial and historical description of the potential niches in the Neotropics (Fig. 5). All algorithms showed high performance, with greater responses for hydric  (Table 1 and Fig. 3). The spatial distribution patterns predicted using edaphic variables showed higher amplitude and spatial connectivity than those predicted in a similar way for climatic and hydric variables (Fig. 4). In the evaluated spatial scale, the potential niches of B. rubescens calculated from climatic variables presented smaller spatial amplitude (spatially restrictive) compared to those obtained from edaphic variables (not spatially restricted). Yet, historical patterns showed high spatial overlap and geographic fidelity, with increases in spatial amplitude from the past to the present and decrease from the present to the future (Fig. 5).

Brosimum rubescens and climate
The SDMs predictions obtained from climatic variables showed high performance, strong relationship, and dependence between climate and the spatial distributions of B. rubescens (Fig. 4), confirming results supported and found in multiple works of vegetation SDM with climatic variables (e.g. Peh et al. 2011;Velazco et al. 2017;Soong et al. 2020). Our models showed suitable areas highly connected and spatially continuous, however, with two spatial segregations, one in Amazonia and other associated with the Atlantic Forest (Fig. 5). Climatic conditions are very similar in the Central Amazonia and the Hileia Baiana region within the Atlantic Forest in Brazil (Andrade-Lima 1966;Gentry 1982;Sobral-Souza et al. 2015). Following the Köppen-Geiger (Beck et al. 2018) climate classification, these two regions share similar climatic conditions, corresponding to Af (equatorial), Am (monsoon), and Aw (savanna). Mainly Af and Am are characterized as tropical without (or short) dry season. In fact, the low prediction associated with the Amazonia-Cerrado transition (Aw type) corresponds to a tropical climate with more extended dry season (Beck et al. 2018).
Topographic variables, global relief model (GRM) and elevation were also selected within PCA1: these variables together with pluviometry, conditioned the temporal and spatial connectivity of hydrological networks (Marimon-Junior 2007;Domisch et al. 2015;Parreira et al. 2019) and species distribution (Elias et al. 2019). For some monodominant species, precipitation is important due to its impacts on water table fluctuations and flood pulses (Nunes Da Cunha and Junk 2001). However, there is no relationship between B. rubescens and soil moisture, but a set of other ecological conditions acting together on the mechanisms of monodominance origin and maintenance (Marimon-Junior et al. 2019). For B. rubescens, fruiting and seed dispersal processes are associated with the beginning of rainy season (Rivera et al. 2014) or end of dry season (Marimon and Felfili 2006). Also, ENSO-related extreme drought events affect negatively the recruitment and survival rates of B. rubescens (Laurance 2003;Phillips et al. 2010;Marimon et al. 2020), probably compromising the monodominance maintenance. This species is lightdemanding with episodic recruitment highly dependent on the temporal variations of light and temperature (Marimon et al. 2012). Although variables that could characterize this dynamic were not incorporated into our analysis, some variables such as those of cloud cover could act as proxies for this spatial variation of the environment (Clark and Clark 1994).

Brosimum rubescens and edaphic variables
Edaphic properties have been considered the best predictors of woody species spatial distribution (Velazco et al. 2017;Elias et al. 2018). The models obtained from edaphic variables enabled to correctly predicting the distribution of B. rubescens (Fig. 4). However, their performances were lower than those obtained with other variables, since the spatial distribution of tropical species is driven by the interaction of multiple variables in a set of ecological conditions acting together, in which soil attributes can predominate as a key determinant (e.g., soil texture, exchangeable cations) (Marimon-Junior et al. 2019;Elias et al. 2019). The spatial distribution patterns showed association with the central region of Amazonia and the Atlantic coast that have similar edaphic conditions, such as deep soils, predominantly Oxisols, and Argisols (Quesada et al. 2010;Marimon-Junior et al. 2019). The wide spatial amplitude obtained with edaphic variables suggests that the species could have a greater amplitude of niches, even so, it is subject to other environmental conditions acting together, such as climate and hydrology (i.e., soil texture, drainage/ depth of water table), that limit its presence in some places ( Fig. 4 and Supplementary Material S3).
The models resulting from edaphic variables show high spatial overlap with the areas of highest concentration of interchangeable cations (Supplementary Material S1: Cation exchange capacity). Thus, high concentrations of exchangeable cations in the top soil could act as a competitive advantage allowing the species to reach a wide range of niches (Elias et al. 2019) and, occasionally, reaching the monodominance with other ecological factors acting together (Marimon-Junior et al. 2019). However, some studies suggested that edaphic properties affect B. rubescens according to ontogenetic stage, thus the dependency relationship of this species with the edaphic properties increases in complexity (Brito-Rocha et al. 2016;Elias et al. 2018). In this sense, models generated from edaphic variables, for species with dissimilar edaphic dependencies in their development, could be conditioning the performance of the predictions.

Brosimum rubescens and hydric variables
The spatial distributions of the predictions of suitability generated from hydric variables showed high spatial connectivity, few diffuse areas, and spatial overlap with the Amazonas-Tocantins-Araguaia hydrographic networks and the East Atlantic and East-Northeast Atlantic Basins (Fig. 4). Water resources predictors are primarily responsible for defining hydraulic vegetation conditions (Borchert 1994;Anderegg et al. 2018;Bittencourt et al. 2020). Fluctuations in hydric properties act as a strong environmental filter, spatially segregating phytophysiognomies according to hydric and thermal stress tolerance thresholds (Borchert 1994;Bittencourt et al. 2020;Araújo et al. 2021). In the geographical extension (Leroy et al. 2018) delimited in our case, Neotropics, climate and edaphic properties may affect the predictions on macro-ecological spatial scales (regional or continental) and the hydric variables could be acting on smaller spatial scales (Domisch et al. 2015;Velazco et al. 2017). Since no relationship between soil moisture and B. rubescens monodominance has been recorded (Marimon-Junior et al. 2019), other soil hydraulic conditions, such as water table fluctuations, could be acting as one of the soil hydraulic factors influencing the distribution of the species. However, considering the difficulty of determining the water table depth on a large scale, further studies must be carried out to corroborate this hypothesis.
We observed that soil water stress (SWS) and topographic humidity index (TWI) were related to B. rubescens occurrence. Drought is an important factor related to the hyperdynamic of some tropical forests, affecting the structure, dynamics, and composition of tree species (Marimon et al. 2014b). The regeneration of this species may be conditioned by decreases in the soil moisture derived from droughts during and after ENSO events (Marimon et al. 2020). However, B. rubescens also produce larger amounts of seeds during drought periods (Marimon and Felfili 2006), and this strategy may be advantageous for maintaining its dominance (Rivera et al. 2014), specially due to the chance that some seeds survive the drought.
Brosimum rubescens, CEH-1 and CEH-2 Each set of variables can generate different predictions for a species, but the combination of all these variables in a single subset allows a better representation of ecosystems by improving the approximations of the fundamental niche (Soberón 2010; Velazco et al. 2017). When the different subsets of variables were combined in a single block (CEH-1: Fig. 4 and Supplementary Material S3), the result for B. rubescens showed a balance in the values of AUC and TSS (De Marco and Nóbrega 2018). When the variables are incorporated into subsets, they respond at their thresholds, but when they are complemented, they respond with spatial hierarchical dependence (Rivera et al. 2014;Domisch et al. 2015;Elias et al. 2018).
The environmental suitability matrices obtained from the CEH-1 scenario were adjusted to the spatial distribution patterns of occurrence reported for B. rubescens within the Neotropics (Fig. 4). The spatial distribution pattern of the Brosimeae tribe extends from Mexico to southern Brazil, with records in some Antilles islands (Cuba and Jamaica) and many species associated with the Amazon Basin (Berg 1972). According to this, B. rubescens are spatially discontinuous in the Neotropics, but the highest frequency of records is associated with the Amazon Basin in Brazil, Peru, Colombia, and Guianas. In the Brazilian Amazonia-Cerrado transition, B. rubescens form extensive monodominant forests (Marimon et al. 2001(Marimon et al. , 2020, and in Central Amazonia, this species had wide spatial connectivity in the environmental suitability. As evidenced in the CEH-2 map ( Fig. 4 and Supplementary Material S3), the regions where the balance between climatic, edaphic, and hydraulic factors were greater for this species and had a high spatial overlap with the matrices of environmental suitability obtained from the different algorithms. Other areas where the proportions of C, E, and H, were not in balance, showed the possible limiting factors for the absence. An example of this would be the Patagonia region where environmental suitability was predicted, but it would only be the edaphic variables (and climate in a smaller proportion) that would facilitate the presence of B. rubescens, as the hydric restrictions do not allow its development.

Brosimum rubescens and historical occurrence
In our results, we observed variation in the spatial amplitudes of the environmental suitability projected in different temporal scenarios (Fig. 5). Other studies suggested that decreases in population densities over time in conjunction with predominantly random spatial distribution pattern may be responding to densitydependent mortality, so when evaluating different dominant species of tropical forests, we observed that those with lower spatial densities tend to easily disappear (Felfili et al. 2000;Arieira et al. 2016). After recent drought events, B. rubescens showed variations in population density and demography, mainly in the regenerating stratum, but the robustness of the intraspecific spatial distribution pattern was maintained for adult trees, demonstrating resistance and resilience of the species to exogenous disturbances (Elias et al. 2018;Marimon et al. 2020). In this sense, the extreme climate events may not change the spatial distribution pattern of the Neotropical distribution of B. rubescens, as was observed on a smaller spatial scale (Elias et al. 2018). However, according to the silvigenetic cycle, B. rubescens may be emerging from a possible post-disturbance phase in the construction cycle (Hallé et al. 1978;Elias et al. 2018), suggesting that anthropic disturbances and clearing openings may have greater weight in their demographic and spatial distribution patterns.
It is unknown whether, at present, B. rubescens is in successional or climax stages, and if future climatic changes will affect the balance between mortality and recruitment of this species. Increases in mortality can decrease intraspecific competition and increase the spatial aggregation of this species (Elias et al. 2018). Such a pattern, on a larger spatial scale, was observed in our study in the transition from CURRENT-FUTURE environmental suitability projections (Fig. 5). Moreover, the reverse demographic pattern (decrease in mortality and increases in recruitment), could be a result of the diffuse and disaggregated spatial distribution pattern in the transition of the LGM-CURRENT scenarios (Fig. 5).
Following the BAM diagram (Soberón 2010), our models did not consider the bionomic factors, which may have masked some effects of unknown magnitude on the distribution patterns of B. rubescens in space and time (Rivera et al. 2014). For example, our model did not account for positive associations with species from different strata or local competitive exclusion (Elias et al. 2018) and, by extension, the representation of spatial dynamics and heterogeneity of environments where the species occurs. However, the strategies for spatial occupation and temporal maintenance of the occurrence areas are well characterized in the SDMs (Marimon et al. 2001;Soberón 2010;Elias et al. 2018).

Conclusion
Our SDMs results show the independent and combined effects of different environmental predictors for B. rubescens distribution, allowing the identification of variables that can increase or decrease the spatial amplitudes of its occurrence (Fig. 4). At the same time, we corroborate the spatial persistence and geographic fidelity of the spatial distribution patterns of the species over time (Fig. 5). The SDMs generated for B. rubescens from edaphic variables showed wide spatial amplitudes, greater than when considering just climatic or soil hydric variables, for which the predictions were spatially restrictive. Knowing the responses of species to environmental variables that condition their spatial and historical patterns are important in the global context of human actions and intensification of climate change.
Acknowledgments A special thanks for the contributions and suggestions of Maria Eduarda Maldaner (Dudita), as well as Angélica Faria de Resende and Naraiana Loureiro Benone, which helped to substantially improve this work. This study was financed in part by the Coordenação de Aperfeiçoamento