1. Introduction

Bioregionalization is the effort to objectively determine natural limits of distinct biotas at different spatial scales (Kreft and Jetz, 2010; Morrone, 2014; Smith et al., 2018). Global and continental bioregionalizations have been used as geographical templates to address several ecological and biogeographical issues, such as the amount of biotic interchange among distinct regions (Antonelli et al., 2018), climate change (Li and Zhang, 2017; Reygondeau, 2019), and global biodiversity conservation (Myers et al., 2000; Olson et al., 2001). Some important issues, however, can only be addressed at a regional level, where finding natural boundaries is key to conservation planning and management, as they lead to the conservation of unique and endemic biotas (Silva-Souza and Souza, 2020; Smith et al., 2018).

At a regional level, natural boundaries based on species composition data of taxon-specific groups are referred to as subregions (Moura et al., 2016; Saiter et al., 2016; Silva-Souza and Souza, 2020). These subregions emerge from turnover patterns that are especially complex in biologically rich settings such as Neotropical vegetation domains. Neotropical flora has been shown to be distributed along environmental gradients linked to productivity, such as climate and soil variation (Bueno et al., 2018; Esquivel-Muelbert et al., 2017; Neves et al., 2015; Silva and Souza, 2018a; Souza and Eisenlohr, 2020; Wittmann et al., 2017), that directly influence metabolic rates and performance, and indirectly influence suitability to natural enemies and competitors (Santiago et al., 2016). The correlation between community composition and environmental gradients also depends on the ability of species to overcome natural barriers and colonise habitats with suitable conditions (Pulliam, 2000). In addition, historical events (e.g., the Quaternary glacial-interglacial oscillations) may have influenced the spatial patterns of current flora by driving the expansion of species distribution boundaries, large-scale dispersal limitations, local to regional scale extinctions, and speciation events associated with patterns of habitat stability and refugia (Bueno et al., 2017; Carnaval and Moritz, 2008; Costa et al., 2017; Werneck et al., 2012). Over the last 10, 000 years, human activity in Neotropical biomes may have produced spatial structure in the floristic turnover that has thus far been thus been little studied (e.g., Lewis et al., 2015; Cantidio and Souza, 2019; Silva-Souza and Souza, 2020; Silva and Souza, 2018a). Current and past human activities may alter the spatial structure of plant communities through intensive permanent agriculture, industrial logging, fire, and fragmentation (Lewis et al., 2015). A long history of plant domestication also plays a role in the structure of modern tree communities to an important extent, altering relative abundance and richness of species in a persistent way (Levis et al., 2017).

The Cerrado vegetation in South America represents the most diverse tropical savanna in the world (Silva and Bates, 2002), and favours the study of patterns and processes shaping neotropical systems at a regional scale. This is because (ⅰ) it is distributed in a complex topography with high environmental heterogeneity (Furley, 1999; Sano et al., 2019), which is covered by distinct vegetation types (physiognomies), including seasonally dry forests, forested savannas and moist formations in addition to savannas (IBGE, 2019); (ⅱ) its territory was exposed to many climatic changes during the Pleistocene (Costa et al., 2017; Vidal et al., 2019); and (ⅲ) it has been subjected to high human pressure (Françoso et al., 2015), having most of its territory included in the list of global biodiversity hotspots (Myers et al., 2000; Silva and Bates, 2002). Previous bioregionalization schemes of the Cerrado have focused on the herbaceous layer (Amaral et al., 2017) or were restricted to the woody flora of the savanna physiognomies (Françoso et al., 2020; Ratter et al., 2003, 2011). These contributions represented important advances to our knowledge about the spatial patterns of biodiversity in the Cerrado. However, a number of issues remain to be addressed. Although the tree layer represents a smaller fraction of the Cerrado flora compared to the herbaceous layer (Amaral et al., 2017; Reflora, 2020), knowledge of the spatial patterns of tree species composition is needed because of the key role this group plays in determining ecosystem processes (Hättenschwiler et al., 2008), and habitat structure (Kissling et al., 2008). Furthermore, although the savanna constitutes the predominant vegetation type in the Cerrado (IBGE, 2012), a complete representation of the turnover patterns in the region, including forested areas, is necessary to investigate the spatial patterns of the tree biodiversity in the Cerrado as a whole. Although it is known that different vegetation types have compositional differences (Oliveira-Filho et al., 2021), it is possible that there are other compositional patterns, independent of the physiognomic differences, and thus far unknown. Finally, the Pantanal, a seasonally flooded depression in the upper Paraguay River catchment area (Cunha et al., 2007), has a considerable degree of floristic continuity with the Cerrado, with most of its tree flora originating from this vegetation complex (Junk et al., 2006; Por, 1995; Pott et al., 2011; Pott and Silva, 2015). Despite this, and to the best of our knowledge, no bioregionalization proposal has encompassed different vegetation types as well as the Pantanal flora until now.

Besides the mentioned complexity of vegetation types and floristic continuity in the Cerrado and the Pantanal biomes, the analysis of abundance data may yield macroecological patterns markedly different from those generated by occurrence data. Occurrence-based subregions tend to highlight rarity and endemism patterns (Cantidio and Souza, 2019), because rare species are emphasized in occurrence-based analyses due to equal weights they receive relative to common species (Legendre and Legendre, 2012). In contrast, abundance-based subregions are expected to portray dominance patterns in local communities that reflect demographic success resulting from adaptations related to life history, ecological history, and competitive processes besides habitat suitability (Nielsen et al., 2005).

Here we used occurrence and abundance data with objective analytical methods to propose two bioregionalization schemes for tree species of the Cerrado and the Pantanal in South America. For this purpose, we analysed a tree flora database containing 894 sites, and included rare species as well as different physiognomies (hereafter vegetation types) found in the Cerrado and the Pantanal. We address five questions about the spatial patterns of the Cerrado-Pantanal tree flora and its determinants: (ⅰ) what are the main spatial patterns of the tree flora of the Cerrado and the Pantanal? (ⅱ) do species abundances in the Cerrado and the Pantanal form spatial patterns different from those based on species occurrence? (ⅲ) does the Pantanal form exclusive subregions or is its tree flora just a continuity of the Cerrado? (ⅳ) is the distribution of the Cerrado-Pantanal subregions related to environmental, human, and historical factors? (ⅴ) what is the relative importance of these factors in determining subregion distribution? We expect that areas with different topography and productivity, mediated by climate, soil, and vegetation types, are related to distinct floristic subregions. We also expect that the current spatial structure of the Cerrado tree flora reflects climatic and vegetation changes that occurred during the late Pleistocene. Finally, we expect to find distinct species subregions in areas with different intensities of Pre-Columbian and current human activities.

2. Materials and methods 2.1. Study area

The Cerrado covers about two million km² in central South America, mostly in central Brazil, while the Pantanal covers about 195, 000 km² between the southwestern Cerrado in Brazil and the Chaco domain in Paraguay and Bolivia (Fig. 1a). Together, the Cerrado and the Pantanal span at least 11 Brazilian states and are situated between the main Brazilian biomes (Amazonia, Atlantic Forest and Caatinga; Fig. 1b and c) with vast ecotones between them (Marques et al., 2020; Vieira et al., 2019). We used the official map of Brazilian biomes to delimit our study area in this study (IBGE, 2019). The Cerrado biome is recognised as the most diverse tropical savanna in the world and a global biodiversity hotspot that harbours more than 12, 000 plant species, 35% of which are endemics (BFG, 2015; Myers et al., 2000; Silva and Bates, 2002). The Pantanal biome, a seasonally flooded depression, comprises the world's largest and most important tropical wetlands (Ioris, 2013; Pott and Pott, 1994; Junk et al., 2006; Lima et al., 2020), and has been recognized as a World Natural Heritage Site. The Pantanal provides critical ecosystem services to humans, including buffering of hydrological changes, water purification, and wildlife habitat (Wantzen et al., 2008).

The Cerrado and the Pantanal contain a mosaic of distinct vegetation types that include seasonally dry forests, forested savannas, and moist formations such as semideciduous and evergreen riverine forests in addition to savannas (Bueno et al., 2018; Oliveira-Filho and Fontes, 2000; Pott and Silva, 2015, Figs. 1d and S1). The Cerrado and the Pantanal also comprise 12 ecoregions in which the Pantanal has been identified as a single unit as well as most of the Cerrado (Dinerstein et al., 2017; Olson et al., 2001; Figs. 1e and S2). Mean annual temperature varies from 15 to 28 ℃ in the Cerrado and from 22 to 27 ℃ in the Pantanal (Fick and Hijmans, 2017). In the Cerrado precipitation is more seasonal, with most of the region receiving 800–1500 mm of rainfall annually with a marked dry period from May to September (Fick and Hijmans, 2017; Franco et al., 2014; Fig. S3). There is a high diversity of soil types in the Cerrado, usually attributed to the water regions and the complex landscape in which ancient plateaus alternate with younger depressions (Furley, 1999; Sano et al., 2019). In most of the biome, soils are acidic, deep, well-drained and nutrient-poor, with low cation exchange capacity and low organic carbon content (Fig. S3). In the Pantanal, the soils are more alkaline and are poorly drained, with flooding occurring by overflowing rivers and local rainfall (Fig. S3; Junk et al., 2006).

2.2. Floristic and explanatory data

We compiled occurrence data on tree species composition from 1295 local assemblages. Trees were defined as freely standing woody plants that were classified in the Flora do Brasil project database (Reflora, 2020) as trees or shrubs, including palms but excluding tree ferns and subshrubs. We adopted this criterion because phenotypic plasticity is a common phenomenon in woody savanna species, which makes the distinction between trees and shrubs often arbitrary. For 999 of the 1295 local assemblages, abundance data were also available. These data integrate the Caaporã database of Brazilian and Amazonian domains (native Tupi word for beautiful forest, used by Cantidio and Souza, 2019; Silva-Souza and Souza, 2020; Silva and Souza, 2018a) and the forest inventories of the Minas Gerais and Tocantins states (Haidar et al., 2013a; Oliveira-Filho and Scolforo, 2008). See the 'Species and explanatory data' section in Figs. S4-S5, and Table S1 for a more detailed description of our data set. Methodological differences may bias the results due to variation in the proportion of the local species pools captured in each sampling point. To minimise this problem, we standardised the assemblages used in the analysis according to the range of species richness found in 1-ha samples. Under this criterion, we excluded assemblages with fewer species than the minimum value of species count found in 1-ha samples (15 species), and those with more species than the maximum species count found in 1-ha samples (186 species). A total of 401 local assemblages were excluded, leaving a total of 894 local assemblages for analyses (Fig. 1c, see Appendix B), 658 of which with abundance data. Of the 894 local assemblages retained for analysis, 65.3% were from savanna, 21% were from deciduous forests, and 13.7% were from humid forests (semideciduous and evergreen; Fig. A4). We used the updated taxonomic resources found in the Flora do Brasil 2020 project to harmonize family, genus, and species names (revised in December 2020). Further details on the species data set can be found in Appendices A and B.

To investigate the main correlates of the floristic patterns, we obtained a set of 21 explanatory variables corresponding to environmental, human effects, and historical stability (Table S2). The variables were obtained for each local assemblage. The environmental variables included climatic, soil, fire, and topographic factors, as well as three variables representing the main vegetation types of the Cerrado and the Pantanal. Climatic variables were included mean diurnal range, temperature seasonality, precipitation seasonality, an aridity index, and annual potential evapotranspiration. Soil variables included silt and organic carbon content, bulk density, absolute depth to bedrock, cation exchange capacity, pH and the vertical distance to the nearest drainage, a proxy for soil water availability. Fire data included the fire frequency recorded from 1985 to 2020 (the sum of the values within a buffer of 314.16 km² around each locality). Topographic heterogeneity was represented by elevation and terrain roughness (the coefficient of elevation within a buffer of 314.16 km² around each locality). To represent vegetation variation, we summarized the vegetation types found in the Cerrado and in the Pantanal into three main categories according to Bueno et al. (2018): savanna (including savannas and Cerradão forested savanna), deciduous forest, and humid forest (evergreen/semideciduous forests including riparian forests). We included savanna, humid forest, and deciduous forest as three distinct categorical variables. We included the physiognomic variation as these three vegetation types, and as part of our environmental variables, because biotic interactions such as competition through shade are mediated by physiognomic variation. To estimate historical stability, we used the data on biome stability since the Last Glacial Maximum (21, 000 yr BP) provided by Costa et al. (2017). This variable consisted of counting how many times a change in climate was predicted to shift a grid cell between forest and non-forest states for every 1000 years back to 21, 000 yr BP. Human effect variables included the human footprint index, a proxy to current human pressures on the environment (Venter et al., 2016), and the distance to archaeological sites, indicating pre-Columbian indigenous activity (IPHAN, 2020). Although the continuous variables were obtained at different spatial resolutions, they were resampled to 2.5 arc-min (~5 km) to ensure uniformity. See SI for further details on explanatory data (Appendix A).

2.3. Bioregionalization procedure

The bioregionalization approach we used was proposed by Moura et al. (2016) and has been successfully applied in identifying and mapping plant subregions in other South America regions (Cantidio and Souza, 2019; Silva-Souza and Souza, 2020; Silva and Souza, 2018a). The approach is based on unconstrained community-level modelling of the compositional distance between assemblages (Ferrier et al., 2007; Ferrier and Guisan, 2006). We chose this method due to its independence from environmental factors and the possibility of identifying both abundance and occurrence patterns. All procedures were performed separately for the occurrence and the abundance data sets. Here, we summarize the bioregionalization procedure in three steps (Fig. 2). First, we produced a dissimilarity matrix using the Simpson index for the occurrence data set and the dBC-bal index for the abundance data set. The latter measures the turnover component of β-diversity contained in the commonly used Bray–Curtis index (Baselga, 2013). To correct the loss of sensitivity of the dissimilarity indices when compositional change approaches or equals one, we applied a step-across procedure to the dissimilarity matrices (Williamson, 1978).

In the second step, we used ordination and interpolation techniques to model compositional dissimilarity and map the turnover component of β-diversity over the entire extent of the Cerrado domain. We used non-metric multidimensional scaling ordination (NMDS) to represent the dissimilarity matrices into a few dimensions (Figs. S6–S8). We then used the NMDS-scores of each assemblage as input data in an ordinary kriging interpolation procedure (Wackernagel, 2003). We interpolated the NMDS-axes with a high resolution of 2.5 arc-min (~5 km, Figs. S9–S10).

Finally, we performed a partition of the stacked interpolated NMDS-axes of compositional variation using an unsupervised pixel-based classification based on Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA; Ball and Hall, 1967). We used ISODATA because it is a useful method for classifying continuous data surfaces. The ISODATA is similar to the k-means algorithm which has been widely used in biodiversity classifications (e.g. Françoso et al., 2020; Moura et al., 2016; Silva and Souza, 2018a). In both approaches, a predefined number of clusters must be provided. The ISODATA, however, has the advantage of evaluating the statistics of each cluster at each iteration of the algorithm, with pixels being reclassified into revised clusters (Abbas et al., 2016). Therefore, while in k-means the number of classes provided a priori determines the number of resulted clusters, in ISODATA it is only an initial estimate of the number of clusters and the algorithm can change this number by eliminating small clusters, merging, or splitting clusters towards a more natural division (Abbas et al., 2016; Memarsadeghi et al., 2007). We used the L-method to obtain an initial estimate of the best number of clusters (Salvador and Chan, 2004), which was 36 for occurrence-based data and four for abundance-based data as the best estimate solution for each data set (Fig. S11).

To ensure that the abundance-based subregions reflected real abundance patterns rather than occurrence patterns produced by a smaller number of samples (Silva-Souza and Souza, 2020), we applied the entire bioregionalization procedure to the abundance data set (658 local assemblages) converted into presence-absence and using the Simpson dissimilarity index. The abundance data, when subjected to the bioregionalization procedure in the form of presence and absence (Fig. S12), no longer generated the patterns found when the abundances were analysed. Rather, they approximated the patterns found for the entire occurrence data set. This shows that the results generated for the abundance data did not merely represent an artifact caused by a smaller number of points (N = 658) relative to the full occurrence data set (N = 894), but rather patterns produced by the spatial distribution of species dominances and rarities. We then regarded occurrence- and abundance-based subregions as distinct and complementary approaches in this study.

We performed an indicator species analysis to identify the species that best characterized each of the occurrence- and abundance-based subregions (De Cáceres and Legendre, 2009). We also extended our bioregionalization to larger spatial scales. For this aim, we conducted separated hierarchical clustering analyses on the species occurrence and abundance data sets aggregated by subregions to investigate the relationships between the subregions. We used the Simpson or the dBC-bal dissimilarities and the UPGMA linkage method. The UPGMA linkage method was chosen through cophenetic correlation (cc = 0.78; Daru et al., 2020) contrasting the UPGMA with other widely used hierarchical clustering algorithms, and to make the results more comparable with the literature, as many other authors have used this method (e.g., Cantidio and Souza, 2019; Dapporto et al., 2014; Kreft and Jetz, 2010; Souza and Eisenlohr, 2020). Following Holt et al. (2013), the clusters corresponding to at least 90% of the variation in ecological distances were considered as significant. We refer to the patterns at broader spatial scales resulting from the UPGMA classification of subregions as floristic zones. The regionalization steps were performed using software R 4.0.2 (R-Core-Team, 2020) and ArcGis 10.5 (ESRI, 2017). See SI (Appendix A) for further details on the bioregionalization procedure and R packages and functions we used.

2.4. Assessing the drivers of floristic patterns

To reduce the dimensionality and collinearity in our explanatory database, we performed a factor analysis of mixed data (FAMD) on the environmental variables, including the vegetation types as dummy variables. FAMD is an extension of principal component analysis (PCA) in which it is possible to analyse a data set containing both quantitative and qualitative variables (Pagès, 2014). We conducted two separate FAMD analyses, one with environmental variable values obtained from the 894 coordinates of occurrence data and another with environmental values from a subset of these coordinates corresponding to 658 coordinates of abundance data. Following Kaiser's criterion, from each FAMD we retained the first seven principal dimensions (eigenvalues ≥1) for subsequent analysis. Together, these dimensions accounted for more than 75% of the total variance in both occurrence and abundance data sets (Fig. S13). In total, 10 variables were used to investigate the subregion correlates: (ⅰ) the seven FAMD dimensions of environmental variables (Figs. S14–S17), (ⅱ) the two human-effect variables (human footprint index and distance to archaeological sites), and (ⅲ) the historical stability.

We used Multinomial Logistic Regression to investigate the environmental, historical, and human correlates of subregions. The dependent variable was a categorical variable formed by the list of local communities numbered according to the corresponding subregions. We constructed two separate Multinomial Logistic Regression models, one for the occurrence data set and other for the abundance data set. We associated the value of the nearest explanatory variable for each single community. Initially the models contained all 10 explanatory variables for the occurrence and the abundance data sets. We then used a stepwise procedure to select for the occurrence and abundance data sets the smallest set of explanatory variables (Burnham and Anderson, 2002) that minimized the Akaike information criterion (AIC). We considered that environmental, human or historical factors are relevant in explaining the occurrence- and abundance-based subregions distribution if at least one variable linked to these factors was retained in the final models. To control for spatial autocorrelation in the data sets, we used Moran's eigenvector maps (MEMs) following the protocol by Bauman et al. (2018). The minimum spanning tree (MST) was the spatial weighting matrix (SWM) chosen for the occurrence (MST linear), and abundance (MST binary) data sets, because it was the one with which the highest R-squared was reached among a number of candidate matrices (Bauman et al., 2018). Positive and significant (p < 0.05) MEM eigenfunctions were included as explanatory variables alongside the variables selected in the occurrence and abundance regression models (see Appendix A). We evaluated the relative importance of the predictors selected in the best models and spatial structure measured by MEMs through partition of the explained deviance. For this purpose, we grouped the explanatory variables selected in the best models into different predictor sets corresponding to the environmental, historical, and human effects, and to the spatial structure. The occurrence or abundance regression model with the lowest AIC had its deviance partitioned to obtain the unique and shared contribution of the different predictor sets in explaining the Cerrado-Pantanal tree subregions (Moura et al., 2016). To better explore the differences in the explanatory variables included in the final occurrence and abundance models between subregions, we used Kruskal–Wallis and Wilcoxon-rank-sum tests. The analyses were performed in R 4.0.2 (R-Core-Team, 2020). See SI (Appendix A) for further descriptions.

3. Results 3.1. Bioregionalization

The NMDS based on the occurrence of 1714 tree species showed a clear separation between the composition of savannas on the one hand, and humid and deciduous forests on the other (Fig. 3). Based on the L method and the ISODATA algorithm, we found 18 floristic subregions for this data set (Figs. 4ab and S18). The compositional variation based on species occurrence showed a complex pattern, with high interdigitation among the subregions and some discontinuity cases that produced disjunct patches within different subregions. The largest occurrence-based subregions were Paranaíba Cerrados (subregion 16; 212, 766.3 km²), Upper-Paraná Cerrados (15; 171, 184.5 km²) and Mid-Araguaia Cerrados (10; 160, 069.5 km²; Table S3). The subregions with the greatest number of indicator species are Western Minas Gerais Cerrados (subregion 17; n = 152 indicator species), Mato Grosso Parecis Cerrado (8; n = 102) and Lower-middle Parnaíba Cerrados (1, n = 61). Two main occurrence-based subregions were identified in the Pantanal (11- Eastern Pantanal and 12- Western Pantanal). Based on the UPGMA, the 18 occurrence-based subregions were grouped into nine floristic zones that explained at least 90% of the variation in ecological distances (Fig. 5a), with a large zone covering the central and southern portion of the Cerrado and minor floristic zones towards the edges of the domain (Fig. 5b).

The NMDS of abundances also showed a clear separation between the composition of savannas on the one hand, and humid and deciduous forests on the other (Fig. 3). Based on the L method and the ISODATA algorithm, we found four subregions based on the abundances of 1387 species in 658 assemblages (Figs. 4cd and S19). The Northern and Central Cerrados (subregions 2 and 3) were the largest subregions, covering 762098.2 and 1004707.3 km², respectively (Table S3). The subregions with the greatest number of indicator species are Southern Cerrados (subregion 1; n = 332 indicator species) and Central Cerrados (3; n = 117). Based on the UPGMA, abundance-based subregions were grouped into two broad floristic zones with a north-south division that explained at least 90% of the variation in ecological distances (Fig. 5c and d). The bioregionalizations proposed in this study based either on species occurrences or species abundances showed that the three main vegetation types existing in the Cerrado and the Pantanal are geographically subdivided according to the species composition (Figs. 3, S20 and S21).

3.2. Environmental, historical, and human drivers

The best multinomial logistic model related to the affiliation of assemblages to different occurrence-based subregions included all variables entered in the full model (the seven environmental dimensions, human footprint, archaeological sites distance, historical stability, besides MEM occurrence variables). Most of the variability in occurrence-based subregions was explained by the shared contribution (mixed fraction) between environment and the spatial structure (51.5%) followed by the pure spatial structure fraction (36.2%; Fig. S22a). The best model for abundance-based subregions included the first six FAMD environmental dimensions, biome historical stability, and distance to archaeological sites, besides MEM variables. Like the occurrence-based subregions, most of the variability in abundance-based subregions was explained by the shared contribution between environment and the spatial structure (59.9%) followed by the pure spatial structure fraction (36.5%; Fig. S22b).

The Kruskal–Wallis test and the multiple pairwise-comparison between groups (Wilcoxon test) revealed significant differences between subregions for all variables included in the occurrence and abundance models (p < 0.05, Figs. 6 and 7), with an exception for the biome stability variable, which showed no significant differences between the four abundance-based subregions (Fig. 7).

4. Discussion 4.1. Regional patterns

Our bioregionalizations differed from those of previous studies based on biophysical, multi-taxa, and ecoregions (Dinerstein et al., 2017; Sano et al., 2019; Silva and Abdon, 1998). These approaches are often presented as proxies for compositional bioregionalizations (Smith et al., 2018). However, our results suggest that this premise is not accurate for the Cerrado and the Pantanal, and that floristic patterns may be obscured in multi-taxa and biophysical regionalization schemes. For example, both the Cerrado and the Pantanal form nearly homogeneous units in the World Wildlife Fund (WWF) ecoregion map (Fig. 1e), which contrasts with the internal divisions we found in both occurrence- and abundance-based bioregionalizations.

Our occurrence-based patterns point to high compositional turnover in space, resulting in a finer and more complex division than previous compositional regionalizations carried out within the Cerrado boundaries (Pantanal not included) for both herbaceous (Amaral et al., 2017) and savanna-only woody (Françoso et al., 2020) plants. The higher turnover reflected in our maps probably resulted from the larger data set we used, which included vegetation types other than savannas as well as more species. Françoso et al. (2020) found subdivisions (their districts) in the northern and north-western (Rondônia state) parts of the Cerrado that may correspond in part to our occurrence broad-scale UPGMA groupings. Their broad subdivisions of the central and southern parts of the Cerrado, however, partly resemble our abundance-based southernmost subregion 1 and central subregion 3. The occurrence floristic zones we found point to a large central floristic zone of typical and endemic species of the Cerrado and minor floristic zones probably representing marginal enclaves of neighbouring biomes towards the periphery of the Cerrado (Françoso et al., 2016; Vieira et al., 2019). Although the correspondence between a bioregionalization of the herbaceous and woody floras may not be expected to be extensive, the phytogeographic regions found by Amaral et al. (2017) for the Cerrado herbaceous flora also found a north-south divide, with a distinct flora in the northwestern portion of the region, and the periphery of the Cerrado with close floristic affinities with neighbouring biomes. The community-level interpolation and partition techniques we used identified floristic groups regardless of any environmental variable. This allowed us to capture a high degree of spatial interdigitation and several floristic enclaves in our occurrence-based subregions map (Silva-Souza and Souza, 2020). We consider as floristic enclaves the small group discontinuities occurring as islands within or bordering other groups (for example, the discontinuities of subregion 18 occurring within subregion 6, and the discontinuities of subregion 3 bordering subregions 4 and 2). Although it may not be visually clean, floristic enclaves and spatial interdigitation are common in biologically complex regions and have been found in several neotropical vegetation maps (Cantidio and Souza, 2019; Marques et al., 2020; Silva-Souza and Souza, 2020; Silva and Souza, 2018a), and seem to be particularly evident in the Cerrado (Miranda et al., 2018). The occurrence of floristic enclaves has been associated with ecological and historical factors, such as island-like environmental conditions and the frequent expansion/retraction events of forests and arid environments throughout the Pleistocene (Lima, 2021; Silveira et al., 2019). In the Cerrado, these processes may be particularly strong due to the complex topography associated with a mosaic of edaphic conditions, and the extensive ecotones with neighbouring biomes (Caatinga, Atlantic Forest, Amazonia and Chaco), which have advanced and retreated at different intensities over Pleistocene cycles (Costa et al., 2017), giving ample opportunity for the formation of remnant island species groups with similar environmental niches (Silva and Souza, 2018b).

The distribution of Cerrado-Pantanal subregions differed markedly between the occurrence- and abundance-based data sets. Occurrence-based patterns depict commonness and rarity but do not necessarily reflect species abundances. Our occurrence-based subregions were more finely divided than those based on abundances and are indicative of rarity and endemism patterns, which are influential in the detection of turnover patterns (Cantidio and Souza, 2019; Silva-Souza and Souza, 2020). Some of our occurrence-based subregions may reflect the presence of species that occur with limited extension in the Cerrado and the Pantanal, either because these species occupy non-optimal habitats in the Cerrado-Pantanal and are maintained through a regular flux of migrants (mass effect process), because they originate from recent speciation, because they are involved in the process of colonization from neighbouring biomes, or because they are remnant populations of historical retractions of the limits between biomes (Bueno et al., 2017; Costa et al., 2017; Pulliam, 2000; Sexton et al., 2009; Vellend, 2010, 2016). The fine divisions disappeared in the abundance-based subregions, which is indicative of domination by an oligarchy of species that are abundant over vast areas (Bridgewater et al., 2004; ter Steege et al., 2013) and may reflect demographic success resulting from adaptations related to life history, site history, and competitive processes besides habitat suitability (Nielsen et al., 2005). For example, the Pantanal biome boundaries showed unique occurrence-based subregions (i.e., subregions 11 and 12, which were not shared by the Cerrado) but did not show regionally restricted abundance-based subregions. This means that the Pantanal vegetation is dominated by the same dominant species found in the Cerrado, and from this point of view the Pantanal may be viewed as forming a floristic continuity with the Cerrado. At the same time, our occurrence-based results highlight that the Pantanal has sets of rare species that do not occur in the Cerrado. Indeed, it has been recognised that the Pantanal has very few endemic species due to its recent geological origin but is inhabited by several species considered rare or endangered in South America (Junk et al., 2006; Pott et al., 2011), probably species from the dry seasonal forests, Chaco, Amazonia and Atlantic Forest (Pott et al., 2011).

4.2. The drivers of floristic subregions

Our results indicated that the distribution of the Cerrado-Pantanal subregions emerged from multiple processes. Environmental, historical, and human factors were retained in both the best occurrence- and abundance-based models, while the Kruskal–Wallis tests identified differences between subregions for variables related to these factors. Moreover, deviance partitioning analysis indicated shared effects of the variables linked to the environmental factors and spatial structure as well as purely spatial effects on the distribution of subregions.

Besides confirming the importance of environmental factors such as climatic, soil, and fire gradients on the structure of the Cerrado flora (Abrahão et al., 2019; Bueno et al., 2018; Françoso et al., 2020; Neves et al., 2015), our results may help to build more precise hypotheses about which sets of variables should be more important in which specific regions of the Cerrado and the Pantanal. This possibility may be particularly relevant regarding the Pantanal, whose floristic organization in space has been less studied. In this sense, different subregions of the Cerrado and the Pantanal seem to experience distinct disturbance and stress regimes. The Pantanal (occurrence-based subregions 11 and 12), as well as the northern portion of the Cerrado (abundance-based subregions 2 and 4) and its subregions bordering the Caatinga seasonally dry biome (occurrence-based subregions 1 and 18) suffer from more intense droughts according to the distribution of the aridity index (Fig. 6), evapotranspiration, and precipitation seasonality (Fig. 7). Fire frequency broadly followed this pattern, being more intense in subregions in the northern Cerrado (occurrence-based subregions 1 and 7), in the Pantanal plain (occurrence-based subregion 12) and at the boundary between the states of Minas Gerais and Bahia (occurrence-based subregion 18). Each subregion had a unique set of conditions. The occurrence-based subregion 5, for instance, experiences a less arid climate for the Cerrado and the Pantanal yet has a relatively high fire frequency. Specific combinations of climate, soil conditions, fire, and historical and human effects make each subregion we identified a potential metacommunity with differentiated selective pressures for functional traits and which are expected to display particular population and community dynamics. Previous work in the Cerrado indicated effects of soil variation on the plant species composition (Abrahão et al., 2019; Bueno et al., 2018). Our work confirmed these effects. For example, the two subregions that characterized the Pantanal flora with occurrence data are associated with soils that were richer in organic matter and less well-drained than the Cerrado soils. We thus expect that the Cerrado vegetation dynamics framework proposed by Bueno et al. (2018), in which humidity, fire, and soil fertility drive the proportions of different vegetation types and species composition, should work differently in each of our subregions.

We found that the main Cerrado and Pantanal vegetation types were associated with different environmental variables, and harbor distinct groups of species, as found previously (Bueno et al., 2018; Felfili and Silva Junior, 1992). However, as found for the Amazon (Silva-Souza and Souza, 2020), the Caatinga (Silva and Souza, 2018a), the Atlantic Forest (Cantidio and Souza, 2019), and for the herbaceous layer of the Cerrado (Amaral et al., 2017), the match between vegetation types and floristic subregions was incomplete, and subregions were not formed by exclusive physiognomies (savanna, humid or deciduous forest; Fig. 6 Dim-5 and 7; Fig. 7 Dim-5; Figs. S20 and S21). This means that the same vegetation types occurred in distinct compositional subregions, as has been found for other biogeographical regions (Amaral et al., 2017; Françoso et al., 2020). Vegetation types have been shown to differ in species composition (Bueno et al., 2018; Santos et al., 2012), besides being proxies for ecosystem productivity and environmental stress variation (Dexter et al., 2018). Our results indicate, however, that there are floristic turnover patterns beyond physiognomy variation in the Cerrado and the Pantanal. The patterns we found based on species composition at both the subregion and zone spatial scales were remarkably different from the complex and fine-scale pattern of vegetation types (Fig. 1d), and this strengthens the conclusion of previous workers that there are relevant patterns of species composition beyond those associated with vegetation types (Françoso et al., 2020; Ratter et al., 2011).

The mismatch between vegetation types and floristic patterns is in part attributable to phenotypic plasticity among species occurring across different vegetation types and that may express distinct habits in different vegetation types (Silva et al., 2015, 2019). Extensive species overlap has been recognized among Cerrado vegetation types (Françoso et al., 2016). However, as the savannas in particular tend to show reduced floristic similarities with both deciduous and humid forests (Bueno et al., 2018; Felfili and Silva Junior, 1992; Hoffmann et al., 2009). This is clearly shown in our NMDS plots, with the savannas forming a compositionally more distant vegetation type in relation to humid and deciduous forests (Fig. 3). Some subregions were dominated by savanna or deciduous forest physiognomies. For example, subregion 18 is clearly dominated by deciduous forest (Fig. 6 Dim-7) with indicator species typical of this vegetation type (e.g., Machaerium scleroxylon, Dilodendron bipinnatum, Tabebuia roseoalba, and Aspidosperma pyrifolium; Haidar et al., 2013b). The disjunct patches of this subregion in the central Cerrado have been related to calcareous outcrop patches and basaltic spills that characterize high fertility soils typically associated with deciduous forests (Bueno et al., 2018; Haidar et al., 2013b). It is remarkable, however, that different compositional subregions were identified in different parts of each vegetation type across the Cerrado and the Pantanal. In the Cerrado, savannas, humid forests, and deciduous forests frequently occur in a spatially interdigitated pattern (Bueno et al., 2018; Miranda et al., 2018), which likely contributed to the complexity found in the occurrence-based subregion distributions. Regarding our abundance patterns, subregion 1 represents a group where semideciduous forests with a strong subtropical influence predominate. Subregion 4 represents a broad group of forests marked by climatic seasonality, including predominantly deciduous, semideciduous and evergreen seasonal forests. Savanna formations predominate in subregions 2 and 3, and fire-adapted species typical of savanna vegetation such as Andira vermifuga (Mart.), Benth, Byrsonima coccolobifolia Kunth, Caryocar brasiliense Cambess., and Eriotheca pubescens (Mart.) Schott & Endl (Simon and Pennington, 2012). are examples of indicator species found for these subregions.

The deviance partitioning showed only a discrete contribution of shared stability and spatial structure to occurrence and abundance patterns. Moreover, no abundance subregion showed differences with respect to historical stability. However, this does not mean that historical processes have not had a strong influence on current floristic patterns of the Cerrado. The influence of historical processes, both linked to Pleistocene changes in species distributions and longer-term speciation and extinction processes, may be reflected in the high fraction of pure spatial structure we found. The high contribution of spatial structure in explaining occurrence- and abundance-based subregions may also indicate that dispersal limitation plays an important role in differentiating large-scale patterns of the tree species in the Cerrado and the Pantanal (Neves et al., 2015). Several linages from adjacent biomes are thought to have colonized the Cerrado-Pantanal during the climatic and vegetation fluctuations of the Quaternary, while speciation events occurred in the more stable areas in the centre of the study region (Bueno et al., 2017; Vidal et al., 2019). Regardless of the source of Cerrado-Pantanal species (colonization events or speciation), many of them may not have had time and ability to disperse and expand their geographic distribution under the continuous climatic changes and vegetation fluctuations between forest and non-forest states during the Pleistocene (Costa et al., 2017; McGlone et al., 2016). According to McGlone et al. (2016), these rapid environmental changes across hundreds of kilometres increased the likelihood of species populations being stranded in less than optimal habitats, an effect they called the 'Pleistocene ratchet'. For many species, this effect can reduce population size and restrict ranges permanently, while gains are limited and temporary. The Pleistocene ratchet effect may be particularly strong in the Cerrado-Pantanal, which has a recent origin (Simon et al., 2009) and has experienced large areas of instability during the Quaternary period (Bueno et al., 2017; Costa et al., 2017; Werneck et al., 2012).

The spatial patterns in the studied assemblages of the Cerrado and the Pantanal may be produced to a lesser extent by human activities as suggested by the mixed spatial-human fraction. Studies have shown that the human influence on the structure of tropical plant communities has been taking place since pre-colonial times and occur in a variety of ways, from domestication and cultivation of preferred species to habitat alterations caused by land use and change in the fire regime (Levis et al., 2017; Lewis et al., 2015; Pivello, 2011; Stephens et al., 2019). The retention of distance to archaeological sites in our best model suggests that cumulative effects of pre-Columbian indigenous activities have left persistent marks on the spatial structure of tree biodiversity in the Cerrado and the Pantanal. The Central Cerrados (abundance-based subregion 3) suffered the greatest influence from pre-Columbian activities (Fig. 7). Indeed, fruit (e.g., Hymenaea stigonocarpa, Caryocar brasiliense, Eugenia dysenterica) and medicinal species widely used by traditional communities (e.g., Stryphnodendron adstringens and Dimorphandra mollis) were among the strongest indicator species of subregion 3, suggesting the cumulative effect of the domestication of these species by indigenous peoples and, possibly, also by traditional populations since the European colonization (Embrapa, 2006; Oliveira and Gondim, 2013; Vila Verde et al., 2003). The influence of historic indigenous activities has also been detected in Amazonia (Levis et al., 2017; Silva-Souza and Souza, 2020), probably through the domestication and marked dispersal of a subset of useful species, which often became hyperdominant in the area of influence of indigenous settlements (Levis et al., 2017). Over the last several decades, however, forests and natural areas in the Cerrado have undergone dramatic fragmentation and conversion due to livestock and agriculture, as well as hydroelectric reservoirs and the expansion of urban areas (Françoso et al., 2015). These factors explain the inclusion of both the current and historical human footprint into the occurrence-based model. Current human impacts on the Cerrado and the Pantanal are highest in the southern part of the Cerrado (Venter et al., 2016), where implementation of roads took place first, and largely comprise open vegetation formations (Françoso et al., 2015, 2020), which highlight the Paranaíba Cerrados (subregion 16) as having the highest human footprint index relative to all other subregions (Fig. 6).

4.3. Limitations and caveats

When considering the occurrence- or abundance-based subregions of the Cerrado and the Pantanal, it is important to be aware of some limitations in our work. Although we have compiled a comprehensive data set for mapping the Cerrado-Pantanal floristic patterns, some areas contain sampling gaps. This is the case for the far western Cerrado and the Pantanal. Well-spread sample points across the study area are particularly important for the modelling approach we use, since it is based on spatial autocorrelation of the tree communities. It is possible that regions, such as the Pantanal, have had their floristic complexity under-represented due to the existence of large areas without sampling points, considering the large number of rare species known for the area (Junk et al., 2006; Pott et al., 2011). The need for further sampling is more acute for abundance data because available data on species abundance are scarcer compared to occurrence data. The number and limits of subregions for under sampled areas were, thus, more speculative and should be seen as preliminary frameworks to be improved in the future. Another limitation emerges from the heterogeneity produced by different sampling methods contained in the databases we used. Different sampling methods may have added bias to our results by adding variation in the inclusion of different species or generating different dominance estimates due to differences in inclusion criteria and sampling areas (Beck et al., 2013). As we have argued elsewhere (Silva-Souza and Souza, 2020), sampling variation is difficult to avoid in large-scale studies and with different vegetation types. We chose to evaluate the largest amount of information we could amass concerning the tree biodiversity of the highly heterogeneous Cerrado and Pantanal, coupled with measures aimed at minimizing the effects of variation in sampling methods. These were the relativization of species abundances in each local assemblage, the exclusion of assemblages with species richness out of the range found in 1-hectare plots, and the use of the Simpson and dBC-bal dissimilarity indices, which are not influenced by richness differences between samples (Leprieur and Oikonomou, 2014).

Regarding the assessment of the drivers of the main subregions, possible limitations may be associated with the explanatory variables obtained and the statistical design. The time interval of the different environmental variables ranged from 1970 to 2020, whereas the human footprint variable ranged from 1993 to 2009. In principle, it is possible that this difference has added some degree of noise to the results. However, it is unlikely that this effect qualitatively altered our results, because both the pace of climate change and that of agricultural activities linked to fire change act on time scales far greater than the difference in years between the series we employed (IPCC, 2013; Lewis et al., 2015). We used deviance partitioning and rank sum tests as complementary approaches to disentangle the relationship between the subregions and possible correlates. However, the variation of some environmental variables overlaps extensively, and we could only include them as representatives of the FAMD dimensions. For example, fire is a recognized key factor in the functional and compositional structuring of savannas in general, and of the Cerrado in particular (Bueno et al., 2018; Hoffmann et al., 2009). In our data set, this influence appears as one among many others acting in concert and in a correlated manner on the floristic structure of the domain. This means that, with observational data, it is difficult to separate the effects of burning from the correlated effects of climatic, edaphic, and physiognomic factors. In addition, our fire data reflects recent and heavily anthropogenic time-series, and probably does not reflect the historical fire that shaped the internal compositional structure of the biome (Daru et al., 2016).

4.4. Conservation implications

Protected areas represent only 8.3% of the Cerrado, and most of them correspond to sustainable use categories that are known to be less effective in preventing deforestation (Françoso et al., 2015). Recently, the Cerrado and Pantanal have experienced a substantial increase in fire foci, and a 10-year ban on sugarcane cultivation in the Pantanal has been canceled (Lima et al., 2020). Understanding the natural divisions of biodiversity can lead to more effective conservation strategies (Marques et al., 2020), which become even more important in biodiversity hotspots like the Cerrado. Additionally, the identification of areas of conservation importance often needs the combined information from both occurrence and abundance models (Estrada and Arroyo, 2012; Nielsen et al., 2005). Finally, the relevance of environmental factors to the subregion patterns we found alert us to the profound impact on the spatial organization of the Cerrado-Pantanal tree flora in the context of global warming, which can cause the expansion of subregions linked to fire, rainfall seasonality, and human footprint, and the contraction or disappearance of subregions linked to the more humid and less disturbed habitats. The increase in frequency and severity of droughts may accelerate and considerably increase the effects caused by the Pleistocene ratchet, leading many species populations to collapse.

Acknowledgements

We thank Augusto C. Silva, Gabriela S. Faulhaber, Isabelle Cristina M. Coelho and Luiza S. Cantidio for their help with data entry in the caaporã database. We are grateful to Professor Fernando R. da Silva (UFSCar) who provided the R script for deviance partitioning. We thank José L.A. Silva for the constructive discussions on classification methods that helped in the implementation of the bioregionalization analyses. We are thankful to Carlos S.D. Fonseca, Gustavo B. Paterno, Marcelo F. Moro, Miriam P. Pinto, Vanessa L. Rezende, and Vinícius L. Dantas and three anonymous reviewers for useful comments on earlier versions of the manuscript. This research is registered in the Sistema Nacional de Gestão do Patrimônio Genético e do Conhecimento Tradicional Associado (SISGEN) under code AA38B00. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

Author contributions

A.F.S. conceived the research idea; M.G.P. and V.C.S. implemented the database and prepared some figures; R.F.H. provided part of the data; K.J.P.S.-S. analyzed the data and drafted the first version of the article; K.J.P.S.-S. and A.F.S. wrote the final version of the article; K.J.P.S.-S., A.F.S., and R.F.H. discussed the results and commented on the manuscript for approval of the final version.

Data availability

The data that support the findings and the analysis code are available on request from the corresponding author. The spreadsheet with the indicator species and the shapefiles of the maps of subregions and floristic zones of occurrence and abundance are openly available in Figshare at https://doi.org/10.6084/m9.figshare.19723603.

Declaration of competing interest

The authors declare no conflict of interest.

Appendice A and B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2022.09.006.