Multi-temporal ecological niche modeling for bird conservation in the face of climate change scenarios in Caatinga, Brazil

Study area

Caatinga covers 844,453 km² and is the only domain restricted to Brazil. This Brazilian biome is part of a complex of forest vegetation types with unique characteristics throughout the Neotropics, named Seasonally Dry Tropical Forests (SDTFs; Pennington, Lewis & Ratter, 2006). This biome extends from 2°54′S to 17°21′S latitudes and occupies the Brazilian northeast and parts of Minas Gerais (Andrade-Lima, 1981; Paim & Franca-Rocha, 2009; Tabarelli & Da Silva, 2003). It is considered one of the most complex ecoregions worldwide (Nimer, 1972; Reis, 1976). The Caatinga is situated at the convergence zone of several highly unstable air masses. Also, this biome is marked by Brazil’s strongest insolation, high thermal averages (26–29 °C), low relative humidity percentages, and scarce and irregular rainfalls (annual 250–800 mm). Finally, it is a reasonably seasonal domain, with rains distributed through a short period of the year (2–3 months) and extensive periods of cyclic droughts (Nimer, 1972; Reis, 1976).

Herein, we assessed 127 protection areas (PAs) distributed throughout the Caatinga domain (Fig. 1). Of those, 38 are Strictly Protected Areas (SPA), 52 are Sustainable Use Areas (SUA), and 37 are Indigenous Lands (IL).

Target species and occurrence dataset

In this study, we selected 19 bird taxa from Caatinga that can be more sensitive to climate change, all the bird taxa included in the analysis were based in the following criteria: restricted distribution (endemic, near-endemic) and/or threatened status by The International Union for Conservation of Nature red list - IUCN (see Table S1 and Table S2 in Supporting Information).We gathered distribution data for each taxon based on literature reports, online databases Species Link (http://splink.cria.org.br/), Xeno-canto (http://www.xeno-canto.org/), Wikiaves (http://www.wikiaves.com.br/), Vertnet (http://www.vertnet.org/index.html), Ebird (http://www.ebird.org/), specimens held in Museu Paraense Emílio Goeldi, and personal data (PVC, GSRG, MPDS). All the occurrence records were adequately checked, undergoing a strict analysis process by bird experts. Therefore, we obtained a total of 1,523 unique occurrences (Table S1). We followed the taxonomy adopted by the Brazilian Committee of Ornithological Records (Pacheco et al., 2021). We controlled sampling bias on georeferenced data by removing duplicate records and leaving a single record per pixel. We used a thinning technique using the package spThin (Aiello-Lammens et al. , 2015) to reduce autocorrelation in occurrence data. Afterward, we used Moran’s I and variograms that minimize spatial autocorrelation to define the thinning distance (Andrade, Velazco & De Marco, 2020).

Environmental data

For the current scenario, we used 19 environmental variables available at WorldClim version 2.0 (http://www.worldclim.org/current; Hijmans et al., 2005) at a scale of ∼5 km². To avoid collinearity between the climate variables, we performed a principal component analysis (PCA) to reduce the dimensions, thus using as predictors of the response functions the axes that explained 95% of the total variance (De Marco & Nóbrega, 2018). To maintain the dimensionality of climate data over time, the coefficients obtained from the PCA performed with present climate data were used to compute scores for future climate data (Sillero & Barbosa, 2021). For future scenarios, we used the same climatic variables from the general circulation models (Atmosphere-Ocean General Circulation-AOGCMs) available for 2070. We performed a cluster analysis to select the AOGCMs proposed by Varela, Lima-Ribeiro & Terribile (2015) because it aims to maximize uncertainties between climate models. The selected AOGCMs were: CCSM4 (CC), HadGEM2-AO (HD), IPSL-CM5A-LR (IP), MIROC-ESM (#) (MR), MRI-CGCM3 (MG). We used two representative concentration pathways (RCP4.5 and 8.5) to assess the effect of climate change on target species distribution.

Data partition and modeling proceedings

We used the “leave-one-out” (LOO) method proposed by Pearson et al. (2007) to predict the distribution of species with less than 20 occurrence records so that we could partition the occurrences and evaluate whether the produced distributions were reliable. The LOO method is recommended when modeling species with a small number of occurrences since it indicates if they are not within the modeled distribution area (Lima-Ribeiro & Diniz-Filho, 2012). For our species the lowest number of occurrence records was six, in the literature it is estimated that the minimum number of occurrence records should be three to five in order to have reliable predictions (Almeida, Côrtes & De Marco, 2010; Lima-Ribeiro & Diniz-Filho, 2012). Thus, we jackknifed the observations to produce occurrences subsets with n-1 occurrences for our data input. After, we used these subsets to predict the distribution of the species. The occurrence record that was left out was used to evaluate the goodness-of-fit of the models. For distribution predictions to be considered reliable, p values must be less than 0.05, indicating that there would be no relevant sampling bias in the occurrence data set used to model species distributions. Conversely, p-values greater than 0.05 indicate an unreliable outcome.

For species with more than 20 occurrences, models were generally evaluated by cross-validation, considering the independent data set. However, it is necessary to have some caution when maximizing the independence of training and test subsets (Roberts et al., 2017). Therefore, we have partitioned the occurrence data similarly to the checkerboard partition method (Muscarella et al., 2014; Valavi et al., 2019; Velazco et al., 2019), using the ENMTML R package (Andrade, Velazco & De Marco, 2020) available on GitHub (https://github.com/andrefaa/ENMTML). Data were geographically partitioned based on grids of different sizes. We then chose the grid size that maximized the independence and environmental similarity between the subsets. In this approach, both subsets are used first to adjust the model and then later to evaluate the model.

Many algorithms are used to predict species distribution based on different statistical approaches and data inputs. Algorithm performance varies depending on the modeling condition (Norberg et al. , 2019; Qiao, Soberon & Peterson, 2015) and is one of the primary sources of model uncertainty (Thuiller et al., 2019). Thus, multiple algorithms allow us to identify the most suitable species and report model uncertainty (Norberg et al. , 2019; Thuiller et al., 2019). We used three algorithms (Maxent, SVM, and random forest) encompassing a range of statistical techniques for modeling species distributions. Maxent, SVM, and Random Forest were fitted using the R software v.3.5.1 (R Core Team, 2018), with the packages maxent v.0.1.2 (Phillips, 2017), kernlab v.0.9-25 (Karatzoglou et al., 2004) and randomForest v.4.6-12 (Liaw & Wiener, 2002), respectively. We evaluated models using the Jaccard Index (Jaccard, 1908). This index calculates similarity between the predictions and observations of the partitioned data, varies between 0 and 1 (Leroy et al., 2018). Index value closer to 1 means greater correspondence between predictions and observations, consequently, lower number of false positives and negatives, and better evaluated models (Leroy et al., 2018).

We made an ensemble forecast procedure to obtain the final ENM for each species (Araujo & New, 2007). The ensemble was obtained by calculating the arithmetic average of the suitability predicted by the best algorithms for each species. Therefore, the final model is the mean of the algorithms whose performances were greater than or equal to the algorithms’ average Jaccard value. For predicting the future distribution of birds, after performing the ensemble forecast to construct a single consensus model of the algorithms for a given GCMs, a new average of suitability values among the five GCMs was estimated to obtain a single future projection for each emission scenario. Subsequently, suitability maps (current and designed models) were transformed into binary maps based on the threshold values calculated from the Jaccard index (Pearson et al., 2007; Marco Jr & Siqueira, 2009). Finally, we performed pre- and post-processing modeling procedures with the ENMTML R package (Andrade, Velazco & De Marco, 2020; Mendes et al., 2020).

Stacked species distribution models

We performed stacked-species distribution models (S-SDM) of all species to obtain species richness maps in the different scenarios in both current and future scenarios. This method has already been shown to be effective in several different situations (Distler et al., 2015; Guisan & Theurillat, 2000; Wisz et al., 2008). First, we stacked the current and future areas to obtain stability areas for each species and thus obtained a common (stable) area for each species in all predictions. Subsequently, we built a richness map with all stable areas by stacking stable regions of all species. For all stackings, we used the raster calculator tool of QuantumGis 2.8.

Effectiveness of protected areas and conservation status of species

We used the shapefiles provided by the Ministry of Environment’s Cadastro Nacional de Unidades de Conservação do Ministério do Meio Ambiente (national registry of PAs; MMA, 2021) that contains Strictly Protection Areas (SPA), Sustainable Use Areas (SUA), and Indigenous Lands (IL) to represent the PAs. In total, the final dataset consisted of 127 PAs. We performed the analysis considering three types of PAs varying according to their protection levels: (1) Only SPA; (2) SPA + SUA; and (3) SPA + SUA + IL (Fagundes, Vogt & De Marco, 2016). The method employed in the gap analysis admits that species with restricted distribution must have all their occurrence area within PAs as they are more susceptible to extinction (Purvis et al., 2000).

However, widely distributed taxa should have at least 10% of their occurrence extension protected. Therefore, species with distribution areas covering less than 1,000 km² must have 100% of their distribution protected, while those with more than 250,000 km² must have at least 10% of their geographic distribution under protection. For species with intermediate distributional range sizes, a calculation is performed through interpolation using a logarithmic transformation, following the methods proposed by Rodrigues et al. (2003).

Finally, considering current and future scenarios, we classified species, regarding the targed of protection, as (1) protected (P) when the target percentage (≥90%) of species distribution size was within PAs; (2) partially protected (PP) when only one portion of the target percentage (<90% ≥ 20%) lays within PAs; (3) gaps (G) when only a percentage (<20%) of the target was within PAs, and (4) not protected (NP) when a 0 part of the target percentage was within the PAs (Frederico, Zuanon & De Marco, 2018; Velazco et al., 2022).

We evaluated the PA effectiveness using a null model approach: the ability of PAs to retain higher richness than would be expected by chance. For this, we compared the number of species within each PA with the estimated number of predicted species found within the PA according to a null model. This null model randomly allocated PAs within the Caatinga, maintaining their size, shape, and orientation (see Lemes & Loyola, 2013; Ferro et al., 2014; Ribeiro et al., 2016). In each run, we calculated the average value of species richness based on the cells encompassed by each PA for both current and future scenarios. PAs were effective if their observed species richness was greater than or equal to the null species richness obtained from randomizations in at least 95% of runs (i.e., PAs with p < 0.05).

We also identified the representativeness degree of the species within the PAs in current and future scenarios (Araújo et al., 2011). We calculated the representativeness degree within the PAs network as the mean percentage overlap (MPO) for each species. MPO corresponds to the mean percentage of overlap between the PAs in the Caatinga biome and the given species occurring within the PAs. First, we calculated the spatial overlap of each cell with the PAs polygons. Then, we used null models to test whether the MPO of each species was statistically significant, considering the range size of the species. The MPO value observed for a given species was compared to the MPO values obtained from 1,000 random species with an equivalent size interval. This meant the same number of grid cells modeled for the species but extracted randomly within all study extent. Such a procedure allows us to identify whether the representativeness of a given species within the PAs considered, along with their MPO value, was significantly higher or lower than expected at random, considering a significance level of p < 0.05 (Araújo et al., 2011).

All five criteria of the IUCN Red List (criterian A–E) must be taken into account to assess the conservation status of species. However, in some cases, classification can also be done if at least one of the criteria is considered (UNEP-WCMC, 2016). One may adopt this strategy when detailed information about species is unavailable, making it challenging to meet all IUCN criteria for assessing their conservation status. Thus, we calculated the area occupied by the species given current and future scenarios to determine the distribution dynamics of adequate areas for the studied species and assess their extinction risk. We assigned to each climate scenario a threat category from the IUCN and Natural Resources (IUCN, 2019; Thuiller et al., 2005; Akçakaya et al., 2006; Rodríguez et al., 2015). Following the IUCN Red List criterion A3(c), we use the following thresholds for each threat category: Extinct (EX): 100% projected area loss within a maximum of 100 years; Critically Endangered (CR): ≥ 80% area loss; Endangered (EN): area loss <80% and ≥ 50%; Vulnerable (VU): area loss <50% and ≥ 30%; Near Threatened (NT): loss <30%. Despite being a simplistic approach and taking into account only the effects of climate change, it can provide us with an overview of the threats for each species individually, thus being vital for decision-makers in conserving these species. We performed all analyses in R software v4.0.3 (R Core Team, 2020) using the raster package (Hijmans, 2015).

Species distribution models

We built SDMs for 19 bird species from the Caatinga domain. The values for model evaluation (Jaccard) varied between 0.63 (Herpsilochmus sellowi) and 0.96 (Anodorhynchus leari, see Table S1). Sixteen out of the 19 studied species will lose suitable occurrence areas under RCP4.5, and 17 will lose areas under RCP 8.5. The number of highly vulnerable species, which lost more than 40% of their original ranges, was six and 11 in RCP 4.5 and RCP 8.5 scenarios, respectively (Table S2). In RCP 4.5, Anodorhynchus leari (62%) had the most significant area loss. Pyrrhura griseipectus faced a reduction of 91% of its original distribution in RCP8.5 (see Fig. S1 for maps of each species). On the other hand, our results showed that suitable areas would expand for three species in RCP4.5: Nyctidromus hirundinaceus, Penelope jacucaca, and Xiphocolaptes falcirostris, with gains of 0.2%, 3%, and 15% in their geographical ranges, respectively. Penelope jacucaca, and Xiphocolaptes falcirostris, in turn, expanded their ranges by 2% and 14%, respectively, in RCP8.5 (see Fig. S1).

Species richness and stable areas

According to the models for the current scenario, the regions with the highest taxonomic richness are located in a small part of the Caatinga where few PAs exist. Analyzing taxonomic representativeness within the PAs in future scenarios, we observed a loss of areas with high species richness (Fig. 2). Areas of high richness harbor between 13 and 16 species in both present and future predictions. However, we notice a significant decrease in species numbers in the richest areas. This was especially true in the RCP 8.5 projection, where areas became smaller and more sparsely distributed (Fig. 2). A decrease in richness from the present to the two future scenarios is also observed for the other species richness ranges (Fig. 2).

In climatically stable areas, in the overlapping ranges of the species, we observed that there are minimal ranges, for the maximum range of overlap (13–16 species), in the two future scenarios, especially in RCP 8.5, where we observed only 42 km² of range (two cells) (Fig. 3). In the stable area in the RCP4.5 projection, we detected the presence of few PAs: Chapada do Araripe Environmental Protection Area (APA), APA Serras e Brejos do Capibaribe, APA Ambiente das Onças, APA, Boqueirão da Onça, APA Marimbus-Iraquara, Parque Nacional do Catimbau, Estação Ecológica Raso da Catarina, Parque Estadual do Morro do Chapéu, Floresta Nacional de Contendas do Sincorá, Terra Indígena Xukuru (TI), TI Pipipã, TI Kambiwá, TI Brejo do Burgo, TI Pankararé. There was no PA in the stable area with a high overlapping range of species for RCP8.5. We also verified that the range size of overlapping species considerably decreased from RCP4.5 to RCP8.5. In the following range of overlap (9–13 species) are all PAs encompassed in the maximum range of RCP4.5 with the addition of the Serra das Confusões National Park.

PA effectiveness and conservation status of species

In the current scenario, the gap analysis including only SPAs showed that no species was considered protected and one specie is as unprotected (Table 1). Therefore, SPA areas alone cannot effectively protect the bird species analyzed here. Considering the SPA+SUA, no species could be regarded as protected. Seventy-eight percent of species were partially protected, 15% corresponded to gaps, and 5% were not protected. Finally, in the broader protection scenario (SPA+SUA+IL), only one species was considered protected; 73% were partially protected, 15% remained within the gap category, and 5% were not protected.

In future scenario RCP4.5, where only SPA areas were considered, no species met the criteria of the protected status (Table 1). Likewise, the gap analysis, including SPAs, shows that no species met the requirements to be protected (Table 1), while 5% of species were not protected. Considering the SPA+SUA scenario, no species is considered protected, 73% are partially protected, 21% are gap, and 5% fall into the non-protected status. In the broader protection scenario (SPA+SUA+IL), only one species was protected; 68% were partially protected, 21% remained as gaps, and 10% were not protected. In future scenario RCP8.5, where only SPA areas were considered, no species reached the protected status (Table 1), and 10% of species were not protected. Considering the SPA+SUA scenario, none of the species was protected; 73% were partially protected, 21% were gaps, and 5% were not protected. Lastly, in the broader protection scenario (SPA+SUA+IL), we did not find any protected species; 73% are partially protected, 21% fell into the gap category, and 5% were not protected.

The species A. leari, Lepidocolaptes wagler, and P. griseipectus in all scenarios and all categories of PAs were considered as gap species. In contrast, Sclerurus cearensis is considered a gap only in future scenarios regardless of the categories of PAs. Rhopornis ardesiacus in current and future, and all PAs categories met the criteria of a not protected species. Interestingly, these four species are currently classified as Endangered (EN) by the IUCN.

Interestingly, ineffective PAs in the current period remained ineffective in 2070. Under future climatic conditions, for RCP 4.5, we remain with <9% of PAs considered ineffective, while in RCP 8.5, this value decreased to 6% (Table S3).

The MPO between the distribution of bird species and the network of protected areas (PAs) in the Caatinga revealed inadequate protection. In the current scenario, the average percentage of overlap was 0.11% (ranging from 0.02 to 0.15%) (Table S3). In the future, the average percentage of overlap was 0.11% (ranging from 0.04 to 0.13%) for RCP 4.5 (Table S3) and 0.09% (ranging from 0.01 to 0.14%) for RCP8.5 (Table S3). The MPO relationship evidenced a more precarious conservation scenario for species within and below the random range. The protection of existing PAs was respectively equal to or worse than expected by a random distribution of the existing PAs (Table S3). Within the scenario of poor protection of Caatinga bird species, none had an MPO value above the expected by chance in the present and the future, representing 52.63%, 21.05%, and 21.05% of the species presented MPO not significantly different from the expected by chance in the current and future scenarios, respectively (Table S3). In even less favorable conditions, 47.36%, 78.94%, and 78.94% of bird species in the Caatinga were less representative than expected in current and future scenarios (RCP 4.5 and RCP 8.5), respectively (Table S3).

According to IUCN criterion A3(c), for the RCP 4.5 scenario, 42% of Caatinga bird species and 68% of species are projected to be threatened in RCP 8.5 by 2080 (Table S2). Our results also indicated that for RCP 4.5, four species were classified as EN and four as VU. For RCP 8.5, there will be three CR, four EN, and six VU species. A. leari, P. griseipectus, and S. cearenses had the most significant impact on their threat status, especially the RCP 8.5 scenario, classified as CR species. A summary of the impacts of future climate change on the conservation status of each species is shown in Table S2.

Protecting areas in the current and the future scenarios

Richness in all investigated scenarios and stable areas maps suggest essential areas for the biome’s southern, central and eastern regions. The southern portion of Caatinga is not currently protected. However, it is an area of high species richness and has been identified as a high-priority biodiversity conservation spot by the Brazilian Ministry of the Environment (MMA, 2021). Furthermore, this region is home to sites of interest for establishing new areas of integral protection due to extensive natural remnants (MMA, 2021; Antongiovanni, Venticinque & Fonseca, 2018; Antongiovanni et al., 2020). The south of the Caatinga is fundamental for the conservation of cacti species (Carvalho et al., 2022) and one of the few areas in eastern Brazil capable of safeguarding different Cactaceae species in climate change scenarios in the next 60 years (Simoes et al., 2019). The east and central parts are also crucial for conserving endemic plants from the Caatinga dry forest (Silva et al., 2019), corroborating the pattern found in the present study. In our predictions, we found that the richest site is located in the central region of the domain. However, the natural vegetation of this area is not only currently highly fragmented (Antongiovanni, Venticinque & Fonseca, 2018) but also is very susceptible to desertification (Salazar, Nobre & Oyama, 2007) due to the high temperatures, decreased rainfall, soil degradation, and their combination (Darkoh, 1998; Geist & Lambin, 2004; Sivakumar, 2007).

Although protected areas provide an essential service in protecting species, their current configuration in the Caatinga is inefficient for conserving the bird species studied herein, even considering all protection categories (PI+US+IL). Most protected areas we assessed tend to protect areas of low species richness, and just a few PAs protected high species numbers. The role of PAs in protecting species has been discussed regarding several taxa in different domains (Lemes, Melo & Loyola, 2014; Carvalho et al., 2017; Oliveira et al., 2017). For instance, PAs do not preserve restricted/rare Odonata species in the Cerrado (Nóbrega & De Marco, 2011). PAs in the Amazon do not fully protect fish biodiversity in watercourses. Prieto-Torres et al (2018), when evaluating the PAs of several neotropical dry forests, found that the current network covers less than 15% of the distribution ranges of 80% of bird species. The authors also detected that it would be necessary to double the area enclosed within the PAs so that their biodiversity could be more efficiently sheltered. The ineffectiveness of PAs in Caatinga might be directly related to their small number and size.

Another issue that has been discussed is which types of PAs are more effective for biodiversity conservation: more restricted use sites (e.g., National Parks, Biological Reserves, Ecological Stations) or sustainable use units (e.g., National Forests, Extractive Reserves; Locke & Dearden, 2005; Sims, 2010; Ferraro et al., 2013; Carranza et al., 2014). Teixeira et al. (2021) found that fully protected PAs reduce deforestation within their limits. This author also emphasized the ineffectiveness of sustainable use of PAs, especially the ones classified as “Áreas de Proteção Ambiental” (APAs). The same issue has already been reported in other domains such as the Cerrado (Carranza et al., 2014; Françoso et al., 2015; Paiva, Brites & Machado, 2015) and the Amazon (Nepstad et al., 2006; Adeney, Christensen Jr & Pimm, 2009; Soares-Filho et al., 2010; Nolte et al., 2013).

PAs cover 9% of the Caatinga region (Brasil, 2020); only 1.8% are fully protected areas, while 7.2% consist of sustainable use areas, where the direct use of natural resources is permitted by law (Brasil, 2020). According to our gap analysis, it is worth noting that Sustainable Use PAs are of great importance since they partially protect our target species. Thus, it becomes essential to reassess the protection category of these PAs, given that sustainable use areas are not quite efficient for containing human disturbances, as we previously mentioned (Antongiovanni, Venticinque & Fonseca, 2018; Antongiovanni et al., 2020). Most of the Caatinga PAs face problems related to administrative management, such as the absence of land tenure regularization and the lack of a management plan. Around 92.5% of the Caatinga protected areas do not only have a management plan (Brasil, 2020) but suffer from both the lack of human resources and financial support (Drummond, Franco & Ninis, 2009). Along with deforestation, PAs also regularly face synergetic, agricultural, cattle grazing, and burning activities (Feliciano et al., 2003; CSR/IBAMA, 2014).

Therefore, it is crucial to evaluate the role of protected areas in adequately protecting Caatinga’s biodiversity. It is equally important to note that, in addition to climatic conditions, the long-term viability of the endemic bird species of Caatinga will also depend on the persistence of vegetation remnants (which still exist). That is, preserving those remnants still within PAs and creating ecological corridors between the external fragments to make it possible for the species to move across the landscape. In this regard, conservation efforts should be directed towards expanding current PAs, possibly changing the protection category of PAs, and implementing priority areas for conservation proposed by the MMA. In addition, the creation of ecological dispersion corridors within stable areas to avoid population isolation and increase their viability under climate change should be considered (Beier & Noss, 1998).

It is also worth highlighting that all these measures would help us achieve the biodiversity conservation goals proposed in 2013 by the National Biodiversity Commission (Conabio, 2013). One of these proposals was to improve the Brazilian biodiversity conservation status by creating new PAs and other protected areas to protect at least 30% of the Amazon and 17% of the other terrestrial biomes, including the Caatinga (Brasil, Ministério do Meio Ambiente, 2013). However, only 9% of the Caatinga territory is covered by protected areas, and none of our studied species are under the protection of the current PAs. Lastly, we also emphasize the need for environmental education initiatives to increase awareness of the locals about the importance of both biodiversity and PAs.