Habitat‐specific impacts of climate change in the Mata Atlântica biodiversity hotspot

Elucidate the potential impacts of climate changes on the distribution and conservation of the multiple habitats of the Mata Atlântica biodiversity hotspot, which are often treated as a unique entity in ecological studies.


| INTRODUC TI ON
The Mata Atlântica of South America is renowned worldwide for being one of the 36 biodiversity hotspots for conservation prioritization (Mittermeier, Turner, Larsen, Brooks, & Gascon, 2011;Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000;Williams et al., 2011). Less known facts, however, are that (a) the hotspot status is specifically referring to its core vegetation type, the rain forest, and that (b) the Mata Atlântica also houses a diverse and complex mosaic of vegetation types, with their occurrence and distribution determined by the harshest extremes of five key environmental factors ( Figure 1; Neves et al., 2017;Scarano, 2009). Thus, vegetation types are defined here as a plant assemblage and its associated environmental conditions (hereafter 'habitat'). Following Walter (1971), these factors can be classified into azonal (non-climatic) and zonal (climatic). The distribution of azonal habitats in the Mata Atlântica is determined by rocky substrates (rock outcrop dwarf-forests and savannas, henceforth rock outcrop habitats), salinity (white-sand woodlands, henceforth restinga woodlands) or waterlogged soils (tropical riverine forests, henceforth riverine forest), while the distribution of zonal habitats is determined by frost (montane and subtropical riverine forests, henceforth high elevation/latitude forest), drought stress (semideciduous forests) or high levels of rainfall (cloud and rain forest, henceforth rain forest).
In a seminal article, Scarano (2009) argued that environmentally marginal habitats in the Mata Atlântica comprise an impoverished subset of rain forest species that can tolerate the harshest extremes of their environmental conditions. A recent study, however, showed that all Mata Atlântica habitats are strikingly distinct both floristically and environmentally , suggesting that marginal habitats are not simply a nested subset of the more diverse Mata Atlântica rain forest. For conservation purposes, a pertinent takeaway message in Neves et al. (2017) is that a substantial portion of the plant diversity in the Mata Atlântica might be neglected if the spatial design for new protected areas is solely based upon studies that places these multiple habitats together (e.g., Zwiener et al., 2017).
The persistence of biodiversity through such global change will demand biogeographic shifts at all levels of biological organization (e.g. from populations to communities to functional groups, Bhatta, Grytnes, & Vetaas, 2018;Frainer et al., 2017;McLachlan, Hellmann, & Schwartz, 2007, respectively. See also Barnosky et al., 2017, for a recent review).
In the last decades, ecological niche modelling became a major tool to predict the impacts of climate changes on biodiversity, aiding conservation planning in future, dynamic scenarios (Peterson, 2001;Peterson, Egbert, Sánchez-Cordero, & Price, 2000;Peterson et al., 2002). With the development of novel learning machine algorithms (Guisan & Thuiller, 2005) and more accurate climate change predictions (Moss et al., 2010), we are now capable to reduce analytical F I G U R E 1 Distribution of Mata Atlântica habitats in South America (sensu Scarano, 2009) and main environmental factors (arrows) sorting species across these habitats (adapted from Neves et al., 2017). Ellipses indicate zonal habitats, and rectangles indicate azonal habitats uncertainties and provide the much-needed information to support conservation prioritization while accounting for global change scenarios (Elith et al., 2006). This is of particular relevance for biodiversity hotspots, where species are likely to be more susceptible due to its reduced population sizes caused by habitat fragmentation.
Our goal here is to elucidate the potential impacts of climate changes in Mata Atlântica habitats' distribution and conservation.
Because Mata Atlântica habitats occupy distinct climatic and geographic space, our hypothesis is that climate changes will severely impact all habitats, though to different degrees. In addition, because South America will experience increasing temperatures with reduced water availability (IPCC, 2013), we predict that future climate changes will have less severe impacts in restinga, rock outcrop habitats and semideciduous habitats, and more severe impacts in plant communities found at high elevation/latitude and in riverine and rain forests.

of six
Mata Atlântica habitats, with each habitat being distributed across limited ranges of six environmental gradients: rain forest (warm and wet climates), high elevation/latitude forest (environments associated with seasonal cold), semideciduous forest (seasonal drought), restinga (salinity), rock outcrop habitats (seasonal fire and shallow soils) and riverine forests (seasonal soil waterlogging). In order to reflect these limiting environmental conditions in the analyses, we modelled the species of each habitat using distinct geographic delimitations, detailed below (see Figure S2). Spatial scope for species from high elevation/latitude and rain forests comprised the whole extent of the Mata Atlântica and the biogeographical Domains found in the neighbouring South American dry diagonal, namely Caatinga, Cerrado and Chaco. Because species from semideciduous forests are widely distributed across the dry diagonal, their spatial scope comprised the Mata Atlântica, dry diagonal Domains and the neighbouring lowland Amazon (warmer climates). Restingas, riverine and rock outcrop habitats are constrained within conditions that are primarily related to soil. Therefore, despite species from restingas, riverine and rock outcrop habitats having climatic suitability in other habitats (e.g. rain forests), these species are restricted to specific edaphic conditions (e.g. soil waterlogging in riverine forests). Thus, we modelled the potential distribution of these species within their edaphically suitable areas, which we established as the current distribution of restingas, riverine and rock outcrop habitats, respectively. We defined the distribution of the Mata Atlântica habitats, dry diagonal Domains and lowland Amazon in geographic space by creating polygons from a set of points. The 6,243 NeoTropTree sites (points) were previously classified into one of the South American biogeographic Domains and into one of the Mata Atlântica habitats where applicable. The size of each polygon was then estimated based on the distance between a given site and the other sites around it (wall-to-wall map).
Bioclimatic variables were obtained from WorldClim v.1.4 (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). Climatic layers were obtained at a 5-arcmin grain size (~10 km). This spatial resolution is particularly appropriate for this study because species checklists (sites) in NeoTropTree are defined by a single habitat, following the classification system proposed by Oliveira-Filho (2017), contained in a circular area with a 10 km diameter. NeoTropTree data were originally compiled from an extensive survey of published and unpublished (e.g. PhD theses) literature, particularly those on woody plant community surveys and floristic inventories. New species occurrence records obtained from both major herbaria and taxonomic monographs were then added to the checklists when they were collected within the 10 km diameter of the original NeoTropTree site and within the same habitat. The habitat delimitation was conducted using the package 'dismo' (Hijmans & Elith, 2015) in R Statistical Environment (R Development Core Team, 2011).

| Variable selection
Variable selection was very conservative in order to build understandable and ecologically meaningful models ( Figure 2). We followed a multiple-step variable selection routine, consisting of the following: (a) using variance inflation factors (VIF) to identify highly collinear variables, which were progressively excluded through a stepwise procedure. VIFs were computed using two methods: VIFcor (threshold = 0.5) and VIFstep (threshold = 10; see Marquardt, 1970, for method details). We then extracted bioclimatic values from presence points and (b) performed a principal components analyses (PCA) to visualize which variables were more effective in segregating the climatic space of each habitat relative to the climatic space of all other habitats. We also (c) performed PCAs for each habitat separately to assess which climatic variables showed higher correlations with the first three principal components (there was a negligible increase in constrained variance by adding a fourth component). Lastly, we (d) used Pearson's correlation to test whether all variables selected for a given habitat showed low correlation (cut-off = −0.5 < p < .5). We also legitimated the variable selection with literature review, which allowed us to select variables that better represented the climatic space occupied by the species of each habitat, while taking into account their ecological relevance (see Table S3).

| Environmental niche modelling
Models were calculated in three independent cross-validation runs with 30% of data kept to evaluate the model and two evaluation methods (true skill statistic, TSS, and area under the receiver operating characteristic, AUC) for every algorithm available in 'biomod2' r Package (Thuiller, Georges, Engler, Georges, & Thuiller, 2014; generalized linear models, generalized additive models, boosted regression trees, classification tree analysis, artificial neural networks, Bioclim, flexible discriminant analysis, multiple adaptive regression splines, random forest and MaxEnt). We only kept ensemble models with TSS higher than 0.7. We generated 1,000 pseudoabsences through different background areas for each habitat, since they have distinct spatial scopes in our analyses (see section 2.1 and Figure   S2 for more details). A caveat to this approach is the recommendations of Barbet-Massin, Jiguet, Albert, and Thuiller (2012) regarding the use of lower pseudoabsences in some algorithms. Nonetheless, here we followed Thuiller (2014), which points out that the main advantage of biomod2 lies in the capability to compare and combine multiple algorithms using the same set of initial data and parameterization. We controlled for spatial autocorrelation in our models using a generalized least squares framework (Zuur, Ieno, Walker, Saveliev, & Smith, 2009), which consists in modelling alpha diversity as a function of predicting variables using different spatial correlation structures (exponential, gaussian, spherical, linear and rational quadratics) and then selecting the best model (highest delta AIC relative to the null model; i.e. no spatial autocorrelation). We then built a raster with cell sizes as values weighted by presence probabilities to provide a more conservative measure of the potential area occupied by each habitat (Figure 2). Models were projected to CMIP5 data (Coupled Model Intercomparison Project Phase 5; downscaled at 5-arc-minute spatial resolution) using all General Circulation Models available in WorldClim v.1.4 (Hijmans et al., 2005)  into an alpha diversity raster for each habitat and weighted by the F I G U R E 2 Methods summary. Environmental niche models (ENM) were projected to 11 Atmosphere Ocean General Circulation Models (AOGCM) and four representative concentration pathways (RCP) to 2050 and 2070. To calculate potential occupied area, the presence probability rasters were multiplied by cell area rasters, generating a weighted area raster, following two approaches: (i) considering a presenceabsence map with a threshold = 0.5, that is, each cell with presence probability >0.5 sum 100 km 2 (grid cell size) of the total potential area. (ii) Considering that cells could be partially occupied, that is, occupancy models that are either gradually fading (a) or abruptly changing (b) ecotones maximum number of species before generating habitat suitability maps ( Figure 2). Ensemble models were generated for each habitat by first summing their diagnostic species distribution maps and then dividing the resulting map by the number of diagnostic species in a given habitat. This generated a final suitability map (ranging from zero to one) for each habitat.

| Potential area and conservation status
Our models showed that in current climatic conditions, the existing network of protected areas is more effective in protecting the potential distribution of azonal habitats (17.4% compared to only 9.0% in zonal habitats, Figure 3). Semideciduous forest is the least protected habitat, with only 7.1% of its potential area (537,640.29 km 2 ) occurring within protected areas (39,320.39 km 2 ). Amongst azonal habitats, riverine forest was the least protected, with only 8.5% of its potential area (91,492.64 km 2 ) occurring within protected areas (7,816 km 2 ). On average, 13.8% of the potential distributions of marginal habitats are found within protected areas, which is higher than the potential distribution of rain forest occurring within protected areas (10.2%; 41,203.47 km 2 ).
From current conditions to the worst climate change scenario, the high elevation/latitude forest was the least affected, with 5.2% of potential area shrinkage, followed by restinga (7.6%), and rain forest (12.2%). In contrast, future scenarios for semideciduous, riverine and rock outcrop habitats were worrisome. Potential area shrinkage in future climatic scenarios can be as high as 50.4% in semideciduous forest, 58.6% in riverine forest and 66% in rock outcrop habitats.
This loss of climatically suitable areas across all habitats is also reflected in their levels of protection. From current to worst scenario, restinga woodlands are predicted to lose climatic suitability in 6.6% of its currently protected area, followed by high elevation/latitude (8.0%), rain (13.6%) and riverine (55.3%) forests. The current network of protected areas in rock outcrop habitats is predicted to undergo the most severe impacts of climate change, with 60.1% of shrinkage in areas of climatic suitability for species of rock outcrop habitats in these protected areas. Conversely, shrinkage in areas of climatic suitability for species of semideciduous forest (50.4%) will mainly occur outside protected areas (19.0% of protected area loss).

| Distribution of azonal habitats
Riverine forests, which are mainly found in Central Brazil, and rock outcrop habitats, which are mainly found in the transition between Mata Atlântica and Cerrado, are predicted to lose higher levels of climatic suitability in lower latitudes (see Figures S4 and S5). Restinga is predicted to lose lower amounts of suitable climatic space relative to the other Mata Atlântica habitats, suggesting higher climatic stability across coastal white-sand environments in eastern South America (see Figures S6 and S7).

| Distribution of zonal habitats
Our results showed a substantial degree of overlap in the climatic spaces occupied by species from high elevation/latitude, semideciduous and rain forests (for decoupled maps check Figures S8, S9 and S10). This suggests that the abrupt contours that are currently used for delimiting the distribution of these three habitats might be too simplistic (Figure 4). Under current climatic conditions, for instance, our models showed that for 3.7% of the geographic space covered F I G U R E 3 Potential area (in km 2 ) of Mata Atlântica habitats (total area in a and b; protected area in c and d) through scenarios of increase in CO 2 concentration  Figure 5).

| Climatically stable areas
Areas in southeastern Brazil showed a high probability of climatic stability for species of all three zonal habitats ( Figure 6). These potential refugia occur mainly in Rio de Janeiro and São Paulo states.
Potential refugia for species of high elevation/latitude forest are also  Table S1 and Section 2). Habitats were plotted using a red-green-blue colour scheme. The brightest shades of red, green and blue represent the highest probability of occurrence of species from semideciduous, rain and high elevation/ latitude forests, respectively. If a grid cell is potentially occupied by two habitats, its colour will represent an intermediate palette between the colours for these two habitats. However, these climatically stable restinga woodlands are mostly found outside existing protected areas (only 19.9% within protected areas; Figure 6).

| D ISCUSS I ON
Here, we showed that both core and marginal habitats of the Mata Atlântica will be severely impacted by human-induced climate change, though to different, uneven degrees. For instance, considering variation from current conditions to the most pessimistic scenario of climate change in our models (RCP8.5/2070), rain forest is likely to be more climatically stable relative to semideciduous, riverine and rock outcrop habitats, but more impacted than high elevation/latitude forest and restinga woodlands.

| Potential area and conservation status
Through the scenarios, protected areas in riverine forest will have more stable climates across the southeastern portion of its current distribution, highlighting the importance of these areas for protecting viable population sizes of riverine species. Congruent with the results for riverine forest, our future scenario models showed that rock outcrop habitats will lose more climatically suitable areas in their lower latitudes, suggesting that southernmost sites may function as climatic refugia for this hyperdiverse habitat (Neves et al., 2018).
However, given the scattered spatial configuration of these rock outcrop sites, dispersal is likely to be very limited, which suggests that conservation strategies might need to consider new protected areas that connect these outcrop islands through the lowlands. In Upper-right groups represent grid cells where three (high elevation/latituderain-semideciduous, HRS) or two habitats overlap (semideciduous-rain, SR; high elevation/latitude-semideciduous, HS; high elevation/latitude-rain, HR). Habitats with climatic suitability ≥0.33 in a grid cell were considered present. Chord diagrams were made using circlize R package (Gu, Gu, Eils, Schlesner, & Brors, 2014) 1991). Therefore, here we stress that protected areas aiming to secure biodiversity of rock outcrop habitats should not be limited to rock outcrop areas. Rather, effective protected areas should function as ecological corridors connecting multiple rock outcrop sites through lowland environments.
Our models showed that while climate in restinga woodlands are expected to be more stable over time when compared to other habitats, this level of stability is highly variable within its distribution, with central and southern restingas being relatively more stable.
In addition to this uneven impact of climate change across restinga woodlands, coastal environments are also expected to be affected by erosion and sea level raise (EUROSION, 2004;IPCC, 2013). This suggests that conservation planning for restinga woodlands will require a high degree of complexity, with its effectivity depending on strategies that account for geomorphological variation changes associated with both climate and land use change. Restinga has suffered massive fragmentation due to high human occupation in coastal areas and a rapidly developing, unplanned tourism industry.
Amongst zonal habitats, semideciduous forest is predicted to be the most impacted, losing 64% of its current potential distribution under the most pessimistic scenario (RCP8.5/2070).
Moreover, while our models predict climatic stability for species of semideciduous forest in southeastern Brazil, there is a high degree of potential area shrinkage for species of semideciduous forests in northeastern Brazil (see Figure S10). These results therefore suggest that conservation strategies aiming to protect suitable climatic space for these northern species would have to consider corridors that could potentially link their current and future suitable climates. Conversely, high elevation/latitude and rain forests are relatively stable over time, indicating the need for tailor-made conservation strategies for each habitat of the Mata Atlântica.
Nonetheless, biodiversity in these forests is poorly and unequally captured by the current network of protected areas, especially in southern Brazil (Saraiva, dos Santos, Overbeck, Giehl & Jarenkow, 2018). Here, we suggest that accounting for climate change scenarios, in addition to multi-dimensional biodiversity assessments as in Saraiva et al. (2018), might improve current and future conservation strategies for these neglected high elevation/latitude and rain forests.

| Climate change and compositional complexity
Previous studies

| Climatic stability and protected areas
Biodiversity loss from climate change arises because species move to track suitable climate, and agricultural lands, urban development or transportation corridors may stop their movement (Hannah, Midgley, Hughes, & Bomhard, 2005;Heller & Zavaleta, 2009). Protected areas and biodiversity-friendly land uses lessen these barriers to movement (Urban, 2015), but the data needed to inform land use managers require insights from ecologists in which the movements of various species are modelled under multiple climate scenarios. In our models, climatically stable areas are mostly outside the existing protected areas (83.8%).
We, therefore, suggest that the areas identified as climatically stable in our analyses should be incorporated into systematic conservation planning and restoration projects to preserve Mata Atlântica habitats.
Altogether, these areas function as probable refugial areas and climatically stable corridors connecting unstable protected areas to currently protected refugial areas.

| CON CLUS IONS
Our study showed that 'lumping' the natural heterogeneity of the Mata Atlântica can bring great havoc for future conservation strategies, and highlighted three additional factors to be considered in conservation planning for this biodiversity hotspot: (a) we still have little understanding of how climate controls species distribution across the Mata Atlântica, and therefore, the future distribution of species from zonal habitats, namely high elevation/ latitude, semideciduous and rain forests, is highly uncertain. New conservation strategies will need to account for such uncertainty when estimating which areas in geographic space are more likely to protect species from a given habitat and which areas are likely to represent climatic overlaps that are suitable for species from two or more habitats. (b) The maintenance of habitat area through time will likely depend on major biogeographical shifts (see results for semideciduous forests). Thus, new conservation strategies will need to account for the climatic space that will likely facili-

ACK N OWLED G EM ENTS
We thank Marinez Siqueira, Gerhard Overbeck and Demétrio Guadagnin for their insightful considerations on an earlier version of this manuscript. We also thank two anonymous referees for their valuable contributions to this manuscript. D.M.N. was funded by a collaborative research grant from the US National Science Foundation (DEB-1556651) during the time this research was completed.