Pervasive impacts of climate change on the woodiness and ecological generalism of dry forest plant assemblages

Climate emergency is a significant threat to biodiversity in the 21st century, but species will not be equally affected. In summing up the responses of different species at the local scale, we can assess changes in the species quantity and composition of biotic assemblages. We used more than 420K curated occurrence records of 3060 plant species to model current and future patterns of plant species distribution in one of the world's largest tropical dry forests—the Caatinga. While allowing different model extrapolation scenarios, we estimated potential changes in the species richness and composition of dryland plant assemblages in response to projected climate change, and assessed how the ecological generalism and woodiness of plant assemblages can be impacted by the climate crisis. More than 99% of plant assemblages were projected to lose species by 2060, with biotic homogenisation—the decrease in spatial beta diversity—forecasted in 40% of the Caatinga. The replacement of narrow‐range woody species by wide‐range non‐woody ones should impact at least 90% of Caatinga plant assemblages. The exacerbated species loss in the dryland plant assemblages was connected to the heterogenisation and homogenisation of biotic assemblages. Still, the magnitude of climate change impacts on ecological generalism and woodiness patterns of dryland plant assemblages differ according to the direction of the biotic change process. Synthesis. The future increase in aridity will change the patterns of woodiness and ecological generalism of tropical dry forest plant assemblages by decreasing vegetation diversity and complexity. The projected biotic changes in dryland plant assemblages indicate the erosion of ecosystem services linked to biomass productivity and carbon storage. We highlight the importance of long‐term conservation planning for maintaining tropical dry forests.


| INTRODUC TI ON
Climate change has altered the environmental conditions experienced by many species on Earth (Scheffers et al., 2016). If species tolerances do not encompass the novel climatic conditions, they may be forced to change their phenology or geographic range to track suitable climates (Lenoir & Svenning, 2015;Román-Palacios & Wiens, 2020). Spatial changes in the geographic range of species can alter the composition of species assemblages (Olden et al., 2004).
While certain species can colonise new sites in the future, most may not disperse quickly enough to avoid local extinctions, with the extinction risk being greatest for species with low vagility and narrow distribution (Urban, 2015). High local extinctions of narrow-range species and the potential colonisation of new sites by wide-range species can lead to the biotic homogenisation of species assemblages (Socolar et al., 2016), and the eventual loss of ecosystem functions provided by such species (Pecl et al., 2017). Because the climate emergency is a higher threat for tropical than temperate species (Román-Palacios & Wiens, 2020;Urban, 2015), long-term conservation planning will benefit from understanding how different tropical ecosystems are subject to biotic changes (Aguirre-Gutiérrez et al., 2020;Allen et al., 2017).
Among the impacts of climate change figure out the biotic homogenisation of plant assemblages in several ecosystems worldwide (Savage & Vellend, 2015), including drylands (Aguirre-Gutiérrez et al., 2020;Fauset et al., 2012;Vellend et al., 2017). It has been suggested that dryland plants already experience a high water deficit and are close to their climatic tolerances (Allen et al., 2017). One of the most extensive and floristically richest tropical dry forest in the world is found in north-eastern Brazil-the Caatinga-with 912,529 km 2 (Banda-R et al., 2016;Silva et al., 2017). Future climate projections indicate increases in aridity, with subsequent desertification of some areas within the Caatinga (Marengo et al., 2017).
Climate change should drive the range contraction of endemic Caatinga plant species, especially those with more specialised lifehistory attributes (Silva et al., 2019). Indeed, the colonisation and extinction rates of species assemblages may be affected by species geography and life-history attributes (Pearson et al., 2014). Narrowrange species (often rare and with restricted distributions) tend to be more sensitive to climate change, whereas wide-range species often exhibit broader climatic niches and thus high ecological generalism (Slatyer et al., 2013). Among flowering plants, the growth form reflects species ecophysiology (Pearcy & Ehleringer, 1984), as woody plant species likely exhibits limited adaptability to climate change due to their longer generation time and slower rates of climatic niche evolution relative to non-woody plants (Smith & Beaulieu, 2009).
However, woody species may exhibit higher resistance to drought by accessing deep water resources that non-woody species are unable to access (Maeght et al., 2013;Pivovaroff et al., 2016), potentially increasing their adaptability to climate change.
Herein, we apply ecological niche models (ENMs) to estimate current and future patterns of geographic distribution for Caatinga flowering plant species and assess the potential biotic changes in local plant assemblages in response to projected climate change. At the species level, we expect that narrow-ranged and woody plant species will be more impacted by climate change and exhibit more local extinction (i.e. species loss at the assemblage level) than wideranged and non-woody species. Consequently, the relative contribution of narrow-ranged and woody species to plant assemblage will decrease. That is, plant assemblages composed of woody and narrow-ranged species will be more susceptible to biotic homogenisation, whereas those assemblages dominated by non-woody and wide-ranged species should exhibit lower biotic homogenisation, or even be subject to biotic heterogenisation. Given that biotic homogenisation often underlies the systematic replacement of highdiversity biotas by low-diversity and more similar biotas (Filgueiras et al., 2021), we also expect a higher species loss in local assemblages dominated by woody and narrow-ranged species, while those assemblages dominated by non-woody and wide-ranged species might face less pronounced species loss or even species gain.

| Species data
We compiled occurrence records of Caatinga flowering plants based on herbarium records obtained after the year 1970 and available at four online databases: (i) the Global Biodiversity Information Facility (GBIF, 2022), (ii) SpeciesLink (Canhos et al., 2022), (iii) Brazilian Biodiversity Information Facility Repository (SiBBr, 2022) and (iv) Chico Mendes Institute for Biodiversity Conservation (ICMBio, 2022). For GBIF, SpeciesLink and SiBBR databases, we filtered available records based on preserved specimens only.
In addition, we complemented our database by incorporating occurrence records from 332 literature sources (Appendix S1). All binomial names followed the Flora do Brasil (2022) taxonomy, with synonymies and invalid names revised whenever possible or else removed from the database.
We restricted the spatial coverage of the species occurrence dataset to the Neotropical region and recorded a total of 3,942,866 occurrences for 9260 species. Records that were duplicated, with georeferencing errors (records distant less than 1 km from municipality, state or country centroids, or located over water), or uncertain identification (defined above the species level) K E Y W O R D S beta-diversity, biotic homogenisation, climate emergency, drylands, ecological niche models, extrapolation, flowering plants, plant ecology were excluded. To reduce the potential effect of sampling bias and spatial autocorrelation in the occurrence dataset, we randomly filtered one occurrence record for each species within a radius of ~10 km (Kramer-Schadt et al., 2013) leading to 1,328,872 occurrences of 8720 species. Preliminary inspections indicated that many species occurred marginally in the Caatinga. Because our focus was on typical Caatinga plants, we kept only those species with at least 10% of their occurrences within the Caatinga, resulting in 437,745 occurrences of 4744 species. We excluded species with fewer than 15 occurrence records given the high dimensionality of our predictor data (see Section 2.2) and the large extent of our study area, which inevitably contribute to reduce species prevalence (van Proosdij et al., 2016). Our subsequent modelling procedures included 426,875 records of 3060 species belonging to 862 genera and 148 botanical families ( Figure S1).
We gathered data on the growth form of Caatinga flowering species using the Botanical Information and Ecology Network , the Plant Trait Database (Kattge et al., 2011) and the Brazilian Flora 2022(Flora do Brasil, 2022 complemented by pertinent literature (Fernandes et al., 2020;Moro et al., 2014;Zanne et al., 2014). For each species, we assigned one of seven growth form types: tree, shrub, palm tree, woody vine, herb, herbaceous vine or succulent. For species with multiple growth form categories, we assigned the growth form agreed upon by most sources (Engemann et al., 2016). We then grouped the species into two categories: (i) woody (trees, shrubs, palms and woody vines) and (ii) non-woody species (herbs, herbaceous vines and succulents). Growth form was assigned to 2656 angiosperm species (1484 woody and 1172 nonwoody species, see Data Availability Statement section), covering 86.8% of species in our database.

| Current and future projections
We used the 19 bioclimatic variables from the WorldClim database, version 2.1 (Fick & Hijmans, 2017) to represent the current climate.
Bioclimatic layers were downloaded at 5 arc-min (~10 km) and cropped to the extent of the Neotropical realm (our background).
We complemented the current climate conditions with 10 edaphic layers for the depth interval of 5-15 cm representing physical (content of clay, coarse fragments, sand and silk, and bulk density), chemical (cation exchange capacity, total nitrogen, soil organic carbon density, pH) and derived (organic carbon density) soil properties, all extracted from SoilGrids 2.0 (Poggio et al., 2021) at 5 arc-min of spatial resolution. To avoid problems with multicollinearity and reduce the dimensionality of predictor layers, we performed principal component analyses (PCA) separately for climatic and edaphic layers, and retained the first PCA axes that cumulatively explained at least 95% of the variation in the original data (six PCA climatic and six PCA edaphic axes). We used the PCA loading coefficients derived from climatic data to project the linear relationship between raw predictors and principal components onto new layers representing future climate scenarios. We assumed that the current edaphic conditions would remain unaltered along each future climate scenario.
We considered future climate projections for the periods 2041-2060 (hereafter 2060) and 2081-2100 (hereafter 2100) following the 6th Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2021). We used two shared socioeconomic pathways (SSP) for each period: SSP245 and SSP585 as a businessas-usual and non-mitigation scenario respectively. Because the selection of different generalised circulation models (GCM) is recognised as a source of uncertainty in projecting the future habitat suitability of species (Thuiller et al., 2019), we used five GCMs for each combination of the period and SSP, namely BCC-CSM2-MR, CNRM-CM6-1, IPSL-CM6A-LR, MIROC6 and MRI-ESM2-0.

| Ecological niche models
We initially established the species accessible area to restrict predictions of habitat suitability to regions considered reachable by a species within the time frame of projected climate change.
For each species, the accessible area was defined by buffer with a width size equal to the maximum nearest neighbour distance among pairs of occurrences . Within the accessible area of each species, we computed pseudoabsences using the same number of observed presences to maintain the presenceabsence ratio of 1:1 ( Barbet-Massin et al., 2012;Liu et al., 2019).
Pseudoabsences were allocated following the environmentally constrained method, based on the lowest suitable region predicted by a climate envelope (Engler et al., 2004). The models produced with this method have higher discriminatory and explanatory capacities than models when pseudoabsences are allocated randomly in space (Lobo & Tognelli, 2011;VanDerWal et al., 2009).
Since the choice of the statistical method or algorithm can affect predictions from an ENM depending on the initial modelling conditions (Qiao et al., 2015), we computed an ensemble of projections for each species to minimise uncertainty around the ENM method (Araújo & New, 2007). For species with at least 30 non-overlapping presence records, the ensemble included projections for six methods, all computed with linear predictor terms unless otherwise specified: Climate envelope (BIOCLIM), Gower Environmental Distance (DOMAIN), Generalised Linear Models (using linear and quadratic terms), Generalised Additive Models (using smooth terms with three dimensions), Maximum Entropy (using 10,000 background points and default features based on MaxNet package; Phillips et al., 2017) and Random Forests (with the mtry parameter automatically tuned by growing 1000 trees through tuneRF function in randomForest package; Breiman, 2001;Liaw & Wiener, 2002).
Species holding from 15 to 29 presence records were considered rare and modelled under the ensemble of small models (ESM) approach (Breiner et al., 2015). More specifically, we computed bivariate models using all pairwise combinations of our 12 predictors (six PCA climatic and six PCA edaphic axes). The bivariate models were computed with linear predictor terms using four statistical methods or algorithms: Generalised Linear Models, Generalised Additive Models (using smooth terms with two dimensions) and Gradient Boosting Models (using learning rate of 0.1 and 100 trees). For each method and rare species, we obtained the ESM by averaging the habitat suitability of bivariate models weighted by their respective model Somers' D [D = 2 × (AUC − 0.5)] (Breiner et al., 2015). The ESMs computed for the four above-mentioned methods were then used to build an ensemble of projections for each rare species.
In projecting ENMs to either novel regions or time periods, it is possible to estimate habitat suitability for environmental conditions outside the range represented by the training data (Elith et al., 2010). To evaluate model extrapolation for each species projection, we used the mobility-oriented parity (MOP) approach (Owens et al., 2013). Briefly, the MOP analysis involves the calculation of the Euclidean distance between environmental conditions of the training data and conditions in each projected pixel, which is further normalised to 1 and subtracted from 1 to reflect environmental similarity (Owens et al., 2013). For each projected pixel, the MOP analysis allows for restricting the portion of the training data used to compute the Euclidean distance (Owens et al., 2013).
Within species accessible area, we computed environmental similarity between each projected pixel and the nearest 10% of training data points (Montti et al., 2021), and then filtered habitat suitability estimates for projected pixels showing very high (MOP values ≥0.9), high (MOP ≥0.8) and moderate (MOP ≥0.7) environmental similarity with the training data. To minimise issues with unlimited dispersal, we restricted all projections to the respective calibration area defined for each species.
Cross-validation strategies based on geographical partition can improve error estimates when extrapolating (Roberts et al., 2017), but depending on the configuration of spatial partitions (e.g. checkerboard blocks), they can also reduce the number of occurrences and constraint predictor ranges used for model training. Since our projections were constrained to regions holding high similarity with training data, we calibrated all models for the baseline period (current climate) using fivefold cross-validation, with 80% of randomly selected samples used for model training and the remaining 20% used for testing in each iteration. To evaluate model performance, we measured the similarity between predictions and observations using the Sørensen similarity index, which is independent of species prevalence (Leroy et al., 2018). It is necessary to made the species habitat suitability binary according to some threshold value to compute the Sørensen index. We used the species suitability value that maximised the Sørensen index for each algorithm at the baseline period Leroy et al., 2018). For the current climate, and for each combination of GCM, SSP and year, the ensemble model for each species was computed as the average weighted habitat suitability across algorithms, with weights given by the Sørensen index calculated for each algorithm. Similarly, we also computed the average weighted binarisation threshold across algorithms, with weights equal to the Sørensen index (baseline period only). For each future scenario (year and SSP combination), we extracted the average habitat suitability across the GCMs. Then, we binarised the resulting ensemble projection into presence-absence maps using the average weighted binarisation threshold. We used the standard deviation of habitat suitability across the GCMs as future models' uncertainty.
We applied an occurrence-based restriction to keep only the patches of suitable habitat considered reachable by a species (Mendes et al., 2020) to minimise overprediction issues associated with presence-absence maps derived from ENMs. Patches of current suitable habitats are assumed to be reachable by the species if they overlap with a presence record or are within an edge-edge distance threshold of an occupied suitable patch (Mendes et al., 2020). This distance threshold was defined as the maximum nearest neighbour distance among pairs of occurrences of the respective species (the 'OBR' approach in Mendes et al., 2020). Computations were performed in R 4.1.0 (R Core Team, 2021) using the package ENMTML  for non-rare species and the package flexsdm (Velazco et al., 2022) for rare species.

| Spatial patterns of beta diversity, woodiness and ecological generalism
Multiple strategies currently exist to model spatial biodiversity patterns (Pollock et al., 2020). Two of the commonest ones are 'assemble first, predict later', in which species are initially sampled into assemblages and then modelled in relation to environment, and 'predict first, assemble later' when individual species are modelled as a function of environmental conditions and then stacked across geographical units such as grid cells (Ferrier & Guisan, 2006).
While each approach has its strengths and limitations, the 'predict first, assemble later' emerge from the view that species respond individually to environmental change, but it does not explicitly consider the role of biotic interactions in shaping assemblage composition (Ferrier & Guisan, 2006;Pollock et al., 2020).
Consequently, one potential drawback of stacking binary maps is the trend of overestimating species richness across geographical units (Dubuis et al., 2011;Zurell et al., 2020). We assumed this issue to be similar across current and future projections, and used the 'predict first, assemble later' strategy to assess changes in spatial biodiversity patterns of Caatinga plants.
We defined our geographical units using an equal-area projection grid cell of 10 × 10 km of spatial resolution across the entire Caatinga. We overlaid this grid cell system with the binary maps of each species ensemble projections to build species presenceabsence matrices for the current time and each future year-SSP scenario combination (2060 SSP245, 2060 SSP585, 2100 SSP245, and 2100 SSP585). We only considered species presence in a grid cell if they occupied at least 50% of the cell area. The spatial beta diversity for each grid cell was measured by the Sørensen-based multiple-site dissimilarity index, β SOR (Baselga, 2010), computed for the cell set formed by the focal cell and its immediately adjacent neighbour cells. Because the size of the cell set can affect the β SOR value (Baselga et al., 2012), we applied a subsampling procedure to randomly select four neighbour cells around each focal cell 100 times to compute the average β SOR across iterations.
We used the species growth form to compute the proportion of woody species in each assemblage (WoodyProp) to assess potential changes in the assemblage-level patterns of woodiness and ecological generalism. To measure the assemblage-level ecological generalism, we initially classified as narrow-range plant species whose range size distribution was below the first quartile (~100 mil km 2 ) of current projected distribution (vi) WideRatio differ between assemblages subject to biotic homogenisation (∆β SOR < 0) or heterogenisation (∆β SOR > 0). Linear relationships between changes in species richness (∆S) and spatial beta diversity (∆β SOR ), and changes in the relative contribution of woody (WoodyRatio) and wide-range (WideRatio) species were verified through a modified t-test (Dutilleul, 1993) to spatially correct the degrees of freedom of correlation coefficients.
Computations were performed in R using the packages SpatialPack (Osorio et al., 2014).

| Model extrapolation and uncertainty propagation in stacked species assemblages
Each ensemble projection is accompanied by one uncertainty map informing the standard deviation of the species habitat suitability across the five GCMs herein investigated. To obtain the average standard deviation of habitat suitability in each assemblage (grid cell), we first averaged the variances (i.e. the squared deviations) for species in each cell, and then square rooted of it. We used the average standard deviation (AvgSD) of species habitat suitability as an aggregated uncertainty metric for each future year-SSP scenario combination (2060 SSP245, 2060 SSP585, 2100 SSP245 and 2100 SSP585). We further inspected if the aggregated uncertainty correlated with changes in species richness (∆S), spatial beta diversity (∆β SOR ) and relative contribution of woody (WoodyRatio) and widerange (WideRatio) species using a modified t-test (Dutilleul, 1993).
Computations were performed in R using the package SpatialPack (Osorio et al., 2014).
We observed a predominance of assemblages with a higher proportion of non-woody and wide-ranged species in the northern and middle-west regions of the Caatinga, whereas assemblages with relatively more woody and narrow-ranged species occurred in the southern and north-eastern Caatinga (Figure 4). Under the business-as-usual scenario, 93.1% of plant assemblages will experience a reduction in the proportion of woody species (WoodyRatio < 1), slightly less than in the non-mitigation scenario (93.9% of assemblages). The relative contribution of wide-range species will increase (WideRatio >1) in most assemblages in both SSP scenarios (63.88% in the SSP245 and 65.4% in the SSP585, Figure 4b,c). In all SSP scenarios investigated, the increase in spatial beta diversity of plant assemblages was directly related to species loss (Figure 5a,b), with changes in the relative contribution of wide-ranged species not linked to the changes in proportion of woody species (Figure 5c,d).
We found similar results for species-level projections targeting regions with at least 70% or 80% of environmental similarity with the training data, but results differed more substantially when regions with at least 90% of environmental similarity were considered ( Figure 6). By confining the projection space to regions that displayed at least 90% environmental similarity, we found a marked reduction in species range size (Figures S5 and S6). At the assemblage level, the propagation of such effects led to the decrease in species richness ( Figure S19) and increase in spatial beta diversity ( Figure S20) Figures S3 and S4). Overall, the uncertainty across different GCM was higher in species-rich regions ( Figure S15) subject to biotic heterogenisation ( Figure S16), and with the presence of relatively more woody species ( Figure S18).

| DISCUSS ION
Climate change will drastically alter plant biodiversity in one of the world's largest seasonal tropical dry forests, the Caatinga.
Our projections show almost 90% of Caatinga plant species will lose suitable areas in any future scenario of climate change herein considered, particularly narrow-ranged species. At the assemblagelevel, more than 93% of plant assemblages in the Caatinga will lose species by 2060. Biotic homogenisation is expected in about 40% of Caatinga plant assemblages, particularly in species-poor regions currently dominated by non-woody and wide-ranged species.
Overall, our projections point out to Caatinga plant assemblages largely composed by non-woody and widespread plant species in the future.
One of the most important findings in this study is the predominant species loss across >99% of Caatinga plant assemblages by 2060, which can erode ecosystem services in Caatinga by 2060. In drylands worldwide, the role of plant richness in productivity stability is as important as that of climate and edaphic conditions (García-Palacios et al., 2018). Under high aridity, species-rich assemblages are more important for ecosystem stability, whereas functionally distinct species minimise variation in the temporal delivery of ecosystem services at low aridity (García-Palacios et al., 2018). Climate change is expected to increase aridity in Caatinga, particularly in the Aridity may favour the establishment of wide-range plant species (Gallagher, 2016). Since most wide-range plant species in Caatinga are non-woody (Fernandes et al., 2020), the projected increase in aridity in this region will lead to structural changes in vegetation complexity. With higher aridity, dryland ecosystems face a vegetation decline phase due to the reduction of leaf area and canopy cover (Berdugo et al., 2020). Aridification can also promote compositional change (Aguirre-Gutiérrez et al., 2020;Ulrich et al., 2014) and reduce the beta diversity of dryland plant assemblages (Yan et al., 2020).

F I G U R E 3
Assemblage-level metrics across regions subject to different levels of projected biotic change by 2060. Plots on top indicate group comparissons across current patterns of species richness (a), woodniness (b) and ecological generalism (c), while plots at the bottom indicate the respective changes in each pattern (d-f). Each box denotes the median (horizontal line), the 25th and 75th percentiles, the 95% confidence intervals (vertical line) and outliers (dots). Small capital letters denote the results of the Kruskal-Wallis tests for the difference in medians across assemblages subject to different levels of biotic homogenisation (p = 0.05, using Bonferroni correction). Woody and wide ratios above 1 indicate the future increase in the assemblage-level proportion of woody and wide-range species. The computation of all metrics was based on species projections holding at least 80% of environmental similarity with training data. See Table S1 and Figures S10 and S11 for results on complementary projections.

F I G U R E 4
Patterns of woodiness and ecological generalism of plant assemblages in the Caatinga Per assemblage proportion of woody and wide-ranged plant species for the (a) current time and expected relative changes for the future scenarios (b) business-as-usual, SSP245-2060, and (c) non-mitigation, SSP585-2060. Woody and wide-range ratios above 1 indicate the future increase in the assemblagelevel proportion of woody and wide-range species. The computation of all metrics was based on species projections holding at least 80% of environmental similarity with training data. See Figure S12 for results on complementary projections. of each species and applied spatial restrictions to remove (potentially) unreachable patches of suitable habitats and minimise overprediction issues related to ENMs (Mendes et al., 2020). Moreover, evidence available on niche conservatism in plants indicates that the upper limits of plant thermal tolerance are highly conserved, particularly in tropical species (Araújo et al., 2013). Finally, we limited model extrapolation issues by constraining estimates of habitat suitability to environmental conditions similar to those in the training data. While extrapolation constraints can avoid unrealistic projections (Owens et al., 2013), the application of highly conservative constraints can affect the geographical space available for each species. Since wide-ranged species have 'more space' to be constrained relative to narrow-range species, extrapolation masks led to substantial decrease of range size for wide-ranged species. In decreasing the filling of the presence-absence matrix, these highly conservative extrapolation constraints increased beta diversity by restricting the sharing of wide-range species across assemblages (Socolar et al., 2016). Thus, we advise caution in interpreting results derived from projections a single (and potentially highly conservative) extrapolation constraint until further investigations shed light on how species biology interacts with methodological choices.
The impacts of climate change are often expected to be less severe in mountainous regions (Randin et al., 2009). Although elevational gradients can allow species to track more suitable climates over time, the spatial configuration of mountainous areas can limit elevational range shifts (Chen et al., 2011), particularly for woody F I G U R E 5 Change in species richness, spatial beta diversity, woodiness and ecological generalism in Caatinga plant assemblages by 2060. Relationship between (a-b) changes in species richness (∆S) and spatial beta diversity (∆β SOR ) and (c-d) changes relative contribution of woody (WoodRatio) and wide-range (WideRatio) species in the scenarios (a,c) business-as-usual, SSP245, and (b,d) non-mitigation SSP585. Relationship between change in the relative contribution of woody (WoodRatio) and wide-range (WideRatio) species. Symbol colours follow the species assemblage representation in Figure 4a. Pearson correlations at the top of each panel were based on spatially corrected degrees of freedom. The computation of all metrics was based on species projections holding at least 80% of environmental similarity with training data. See Figures S13 and S14 for results on complementary projections. species (Zu et al., 2021). In the Caatinga, the four most relevant highlands-where many narrow-range woody species concentrateare disconnected from each other and located in transitional zones in the south (Chapada Diamantina), east (Planalto da Borborema) and central-northwest (Chapada do Araripe and Serra da Ibiapaba), which impose additional dispersal constraints on woody species. Threats to woody species and their ecosystem functions can be considerably greater than those we project. Vertebrates and invertebrates interacting with woody plants-herbivores, seed dispersers and pollinators-may not be able to couple with the sudden changes in the availability and composition of plant resources (Oliveira et al., 2019(Oliveira et al., , 2021, ultimately scaling up the potential for disruption of biotic interactions (Andrade, Alvarado, et al., 2020).
We have shown how climate change can jeopardise Caatinga plant biodiversity, but much of this region is also affected by chronic disturbances (Antongiovanni et al., 2020) that can operate synergistically with climate change, and intensify the impacts of biodiversity loss on ecosystem functions (Frishkoff et al., 2016). Caatinga already lost half its original cover, and more than 90% of the remaining fragments have less than 500 ha (Antongiovanni et al., 2020). This scenario may not prevent range expansion for some non-woody, self-pollinated and wind-dispersed species, but our projections F I G U R E 6 Impact of extrapolation constraints on species projections and aggregated biodiversity metrics for the current time. (a) Species range size (km 2 ) measured across the Caatinga extent using equal-area projections. (b) Log 10 species richness. (c) Spatial beta diversity (β SOR ). All metrics were computed using species-level projections for regions without extrapolation constraints, and with at least 70%, 80% and 90% (columns) of environmental similarity with the training data. The red dashed line indicates the 1:1 relationship. See Figures S5, S6, S19 and S20 for results on complementary projections.
show an extensive reduction in suitable areas for most woody and non-woody species. Interestingly, the western half of the Caatinga, where biotic changes are most pronounced, also concentrates the most conserved dry forest remnants (Antongiovanni et al., 2020).
This coincident spatial configuration imposes both a challenge and an opportunity to expand the protected area network to assure the connectivity and long-term persistence of Caatinga plant assemblages. Caatinga still figures with around 1% of the original extension covered by strictly protected areas (Oliveira et al., 2017), far beyond the more recent thresholds established under the post-2020 Global Biodiversity Frame of conserving 30% of the Earth's land by 2030 (Xu et al., 2021). As a member of the Convention on Biological Diversity and the sole holder of the Caatinga, Brazil will have a crucial role in conserving the most extensive tropical dry forest in South America.

ACK N O WLE D G E M ENTS
We are grateful to E. G. Moura-Jr and the two anonymous

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no competing interests.

PEER R E V I E W
The peer review history for this article is available at https:// w w w.web of scien ce.com/api/g atew ay/wos/p e er-revie w/ 10.1111/1365-2745.14139.

DATA AVA I L A B I L I T Y S TAT E M E N T
The R-script and raw dataset supporting the results of this work are available at Zenodo Digital Repository: https://doi.org/10.5281/ zenodo.7982973 (Moura et al., 2023).

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information can be found online in the Supporting Information section at the end of this article.