Potential distribution of and priority conservation areas for the Endangered Caatinga howler monkey Alouatta ululata in north-eastern Brazil

Article Potential distribution and priority conservation areas for the endangered Caatinga Howler Monkey Alouatta ululata, Northeastern Brazil 5 ROBÉRIO FREIRE FILHO and JORGE M. PALMEIRIM ROBÉRIO FREIRE FILHO (Corresponding author) Programa de Pós-Graduação em Biologia Animal, Universidade Federal de Pernambuco, Centro de Biociências, Departamento de Zoologia. Av. Prof. Moraes Rego, 1235, Cidade Universitária. Cep. 50670-420. Recife-Pernambuco. Brasil. E-mail 10 [email protected] JORGE M. PALMEIRIM Departamento de Biologia Animal, cE3c - Centro de Ecologia, Evolução e Alterações Ambientais. Faculdade de Ciências, Universidade de Lisboa, Lisboa, Portugal Abstract The Caatinga of NE Brazil is the largest and most diverse seasonally dry tropical forest in 15 the Americas and is home to numerous endemic species. It is declining alarmingly but only 1.2% is under full protection, so there is an urgent need to expand the protected network. The Caatinga howler monkey (Alouatta ululata) is an endangered species nearly endemic to the Caatinga, with a good potential as an umbrella species to protect much of its biodiversity. Using all available distribution data and own surveys we applied Maxent and Zonation spatial modelling to identify the range and 20 priority conservation areas for A. ululata, maximizing habitat quality and connectivity, while minimizing conservation constraints. The top 10% priority areas cover 34,400 Km2 and mostly coincide with good remnants of Caatinga. Only in the northern part of the species range priority areas are protected, so it is essential to create new protected areas in the center and south of the range. Maxent modelling shows that the species depends on good tree cover, but even inside protected areas 25 we observed recent deforestation, illustrating the urgency to improve management. Maxent also revealed that aridity limits the range of the species, so the ongoing aridification of the Caatinga is a menace to its future. In conclusion, the protection of the threatened A. ululata requires establishing new protected areas in priority regions and improving management in those already classified. This challenge is also an opportunity for the conservation of important biodiversity sharing the priority 30 areas of the species. Keywords Maxent, Primate, Species Distribution modelling, Zonation 2 Introduction Located in Northeastern Brazil, the Caatinga covers about 735,000 km2 (Leal et al., 2005), and is considered by many as one of the World’s major wilderness areas (Aguiar et al., 2002). It is 35 the largest and most diverse seasonally dry tropical forest in the Americas and harbors large numbers of endemic species (DRYFLOR). However, its natural vegetation has been declining at an alarming rate (Beuchle et al., 2015) due to land use intensification (Aguiar et al., 2002; Leal et al., 2005). There is an urgent need to take conservation measures to protect the Caatinga, where only 1.2% of the area is under full protection (DRYFLOR, 2016). Protection areas need to be greatly expanded, and 40 charismatic species, such as large primates, can facilitate this process (Ducarme et al, 2013), as they function as umbrella species for the conservation of valuable but more discreet biodiversity in their range. One of the species with greatest potential for this role is the Caatinga howler monkey (Alouatta ululata) a threatened large primate requiring vast areas of suitable habitat to maintain viable populations. 45 The Caatinga howler monkey has an Endangered status due to its small and declining population, consequences of loss of tree cover, habitat fragmentation and hunting (Oliveira & Kierulff, 2008). Most of the species’ range is in the Caatinga although it extends into the Cerrado (Oliveira & Kierulff, 2008), but its limits are poorly known, which is a major constraint for the planning of conservation measures (Oliveira & Kierulff, 2008). Species distribution modelling 50 (SDM) is a tool to map geographic distributions and study how environmental variables influence them (Miller, 2010). It is widely applied in conservation science and its models may support the selection of areas for conservation (Araújo et al., 2002). The most common approach in SDM is maximum-entropy modelling, often applied using Maxent software (Phillips et al., 2006). The identification of priority areas for conservation of species is a prime step to build 55 conservation plans (Pressey et al., 2007). The past few years have seen the development of computational tools to carry out this process of prioritization taking into consideration factors such as habitat quality, connectivity and conservation cost (Kukkala & Moilanen, 2012). C-Plan (Pressey et al., 2009), Marxan (Watts et al., 2009) and Zonation (Moilanen et al., 2005) are examples of approaches and packages developed for conservation prioritization. The general objective of all these 60 packages is quite similar, although the strategy to reach them varies. For example, Zonation, the package used in this study, prioritizes landscapes by iteratively removing the least valuable remaining areas while accounting for connectivity and generalized complementarity (Moilanen et al., 2011). The overarching objective of this project is to make spatially explicit analyses that contribute to the planning of measures to conserve A. ululata and the ecosystems to which it is associated. Our 65 3 specific objectives were (i) to carry out field surveys to collect information on the distribution of the species; (ii) to identify areas where further surveys are needed; (iii) to develop a model to generate a potential distribution species' map and understand the environmental determinants of this distribution; (iv) to identify the areas with most potential for the conservation of the species, (v) and to determine the degree of coverage of the priority areas by existing protection areas. Finally, (vii) we use our 70 results to make spatially explicit recommendations for actions needed to improve the conservation of A. ululata and of the many species that depend on the same Caatinga habitats. Study area The study area includes the known range of A. ululata, across the states of Maranhão, Piauí and Ceará (Fig.1). In the spatial analysis we included not only the area encompassing all the known locations 75 for the species, but also a 30-km wide buffer around it. The aim of this addition was to identify areas that may be suitable for the species but that are outside its currently known range. MethodsData Collection We compiled existing information on the distribution of A. ululata, mostly collected between 2004 and 2010 by the National Center for Research and Conservation of Brazilian Primates 80 (CPB/ICMBio). The number of direct observations by researchers is very low (20 records), so we also used reports obtained in CPB/ICMBio interviews (112 reports). Interview-based distribution analysis can complement direct monitoring data in the case of easily identifiable species (Anadón et al., 2010; Brittain et al., 2018). We also carried out surveys in two regions where information was scarce (August 2016 to May 2017) (Fig. 1). We interviewed 112 farmers and hunters that lived or 85 worked close to areas with natural vegetation. To minimize biased answers we hid that A. ululata was the focus of our interviews. We first asked about other mammals present in the region and only after that about A. alouata. We only interviewed persons comfortable with the the study and respected the interview ethical code of the British Sociological Association (https://www.britsoc.co.uk/media/24310/bsa_statement_of_ethical_practice.pdf). 90 Some records based on interviews were initially referenced with coordinates of the place of the interview, usually farmhouses. We replaced these coordinates by those of the nearest area with a natural environment, within a 2-km buffer around the original coordinates (approximately the distance at which the vocalization of the species can be heard). Points without natural environments inside the 2-km buffer were excluded. This procedure adds locational uncertainty, but Maxent 95 modelling can make useful predictions even when the occurrence data includes a moderate level of locational error (Graham et al., 2008). Our initial database included 184 occurrences, 52 from our surveys and 132 from CPB/ICMBio. However, to minimize problems associated with spatial biases 4 in sampling, we used spatial filtering (Kramer-schadt et al., 2013), reducing the number of occurrences in oversampled areas by using only one in a radius of 5 km. This filtering procedure 100 reduced the number of occurrences used in the modeling to 117 (20 direct observations and 70 reports from CPB/ICMBio, and one observation and 26 reports from our surveys). Modelling of potential distribution To identify variables influencing the distribution of A. ululata and generate a distribution map, we used Maxent software (Phillips et al., 2006). The choice of environmental variables was guided by 105 the species’ biology; A. ululata is arboreal, feeds on leaves, fruits and other plant parts, and lives in a semi-arid region influenced by strongly seasonal rainfall and high temperatures. We expected areas with higher precipitation to be more suitable during the most likely critical period of the year, the dry season. Furthermore, we hypothesized that tree cover and tree height would influence suitability, and that areas with rugged terrain would be more suitable since terrain ruggedness tends to be an obstacle 110 to habitat destruction and to provide some protection from hunting. With the help of a matrix of correlations between all candidate variables, we selected a set of six predictors (Table 1) that were not highly correlated (|r| < 0.70) (Rainho & Palmeirim, 2013) and that were biologically meaningful. Prior to running Maxent all layers were converted to the WGS 1984 geographical coordinate system and to a cell size of 30 arc seconds (about 1 km2), using IDRISI Selva (Eastman, 2012) and 115 QGIS 2.8 (QGIS Development Team, 2016). In Maxent we used the following settings: convergence threshold (10-5), maximum iterations (500), regularization multiplier (1), max number of background points (104), linear, quadratic, product and hinge features, random seed generation and 50 replicates. The resulting map is in logistic format with the probability of presence for each cell ranging between 0 and 1 (Phillips, 2008). To select a suitability threshold for the potential distribution map we used 120 the methodology described in Rainho and Palmeirim (2013), which allows the selection of the smallest area including most occurrences. The area selected was that corresponding to the Maxent suitability threshold 70%, which encompassed 83% of the occurrences (Supplementary Fig. S1); above this threshold the inclusion of more occurrences would force the addition of a disproportionally large area. 125 Prioritizing areas for conservation of Alouatta ululata We used the Zonation package (Moilanen et al., 2005) to prioritize areas for the conservation of A. ululata. Zonation generates a priority map that can be used to inform decision-making. It allows the consideration of cost efficiency in this prioritization, through the inclusion of a “cost layer”. We assumed that conservation cost efficiency is higher in areas with fewer anthropogenic constraints, and 130 5 thus generated a single “constraints” layer combining three thematic layers: human population density, proximity to roads, and anthropic areas (i.e. urban and farmland areas). These three layers, described in Table 1, were given equal weight in the generation of the “constraints” layer. Zonation assigned this layer negative weights and combined it with the map of potential distribution of A. ululata generated with Maxent. 135 In Zonation we selected “distribution smoothing” as the aggregation method. It considers fragmentation to be undesirable and thus retains areas that are well interconnected. The size of the smoothing kernel used was 6 km, a value based on the distances often crossed by various Alouatta species outside their usual home-ranges (Glander, 1992; Crockett, 1998). As a cell removal rule we used “core-area Zonation”, because it is appropriate when importance is given to core-areas, i.e. 140 locations with the highest suitability in terms of abundance or high probability of occurrence (Moilanen et al., 2017). To avoid losing valuable areas and to keep structural connectivity we selected the options “add edge points” and “edge removal” (Moilanen et al., 2017). To identify the priority regions lacking protection, we overlaid the Zonation map with the existing protected areas. Finally, we used a map of forest loss between 2000 and 2014 (Hansen et al., 145 2013) to identify recent deforestation in priority areas. Results The average test area under the curve (AUC) for the Maxent distribution model of A. ululata was 0.857 and the standard deviation 0.032, indicating a high efficiency distinguishing presence from random background locations (Fig. 2). The variables with greatest percent contributions to build the 150 potential distribution model were percentage tree cover, precipitation of driest quarter (bio17) and aridity index (Table 2). The jackknife analysis corroborates these results (Table 2). Probability of presence increases with tree cover, canopy height and terrain roughness (Fig. 3). In the case of the climatic variables (aridity, precipitation in driest quarter and precipitation seasonality) the highest probability of occurrence tends to be in the intermediate values (Fig. 3). 155 Results of the Zonation prioritization are shown in the map on Fig. 4A. For graphical clarity, we only show two levels of priority: high priority (the best 20%) and top priority (the best 10%). These priority areas can be separated in four regions, which are ecologically quite distinct: (1) Mangroves, in the Parnaíba river mouth; (2) Enclaves, encompassing the humid enclaves of Northwest Ceará and Northern Piauí; (3) Caatinga, a vast region fully inside the Caatinga biome; and 160 (4) Border, located along the border between Piauí and Maranhão (Fig. 4A). 6 About 21% of the most important areas for the conservation of the species are inside legally protected areas (Fig. 4A). However, the coverage of priority regions by these protected areas is very uneven; important parts of Mangrove and Enclaves have some level of legal protection, whereas Border and caatinga are virtually unprotected (Table 3). 165 Although the recent (2000-2014) loss of tree cover is scattered, it occurred throughout the high priority areas, including inside protected areas. In some regions of Northern Piauí (Enclaves), deforestation is more widespread and there are larger continuous deforested patches (Fig. 4B). Overall, about 3.8% of the best areas for the conservation of A. ululata has been deforested in the recent past (Table 3). 170 Discussion Potential distribution of Alouatta ululata and variables that influence the probability of occurrence of the species A visual analysis and the means number of the presence points in higher levels of probability of presence show a good correspondence between the distribution and density of the occurrences and 175 the Maxent map. However, there are suitable areas without occurrences, which may be due to local extinctions (e.g. due to hunting), or a lack of survey effort (Oliveira & Kierulff, 2008). There are also a few occurrences in areas with low Maxent suitability, which can be explained by the extensive loss of natural habitat in the region (Oliveira & Kierulff, 2008). An isolated population inhabits a small humid enclave in Acopiara, Ceará, separated from the 180 range of the species by over 100 km of unsuitable dry Caatinga (Fig. 1) (Oliveira et al, 2007). It is probably a remnant of a broader distribution of the species when humid forest dominated the region (Carmignotto et al., 2012), but it may also result from an ancient human introduction. The eastern limit of the range of A. ululata is clear and defined by the high aridity conditions in central Ceará. However, the western limit is ill defined and possibly explained too wet conditions 185 prevailing in Maranhão. Moreover, competition with Alouatta belzebul a close species occurring further west, may help shaping the western limit of the range of A. ululata. In fact, it is unclear if the two taxa are distinct species or just differentiated populations of the same species (Viana et al., 2015). The variables that entered the Maxent model are coherent with the biology of the species (Oliveira & Kierulff, 2008). Percentage tree cover was the most important of these variables and 190 suitability was very low up to approximately 50% tree cover, which is probably explained by the marked arboreal habits of the species (Oliveira & Kierulff, 2008). The folivorous-frugivorous diet of A. ululata may also contribute to make the species dependent on dense tree cover, as trees are the 7 main food source (Oliveira & Kierulff, 2008). Moreover, the results suggest that A. ululata prefers tall forest, even though they occupy a variety of woodland types (Oliveira & Kierulff, 2008). 195 The second and third most influential variables in the model are both climatic and related: aridity index and precipitation in the driest quarter (Bio17). The results show that regions with an extreme dry season are unsuitable habitats, suggesting vulnerability to the ongoing aridification of the Caatinga (Torres et al., 2017): it is predicted that at the end of the century temperatures may be up to 3.5-4.5ºC higher and rainfall 40-50% lower (PBMC, 2013). 200 Priority areas for conservation and current level of conservation The priority areas identified with Zonation (Fig. 4A) should conciliate a high habitat suitability, identified by Maxent, with good connectivity and low conservation constraints. This should make our results a good basis for conservation planning. However, it is important to be aware that both Maxent and Zonation modelling are affected by sources of error, such as inaccuracies in species occurrences 205 and environmental layers, that create uncertainty (Graham et al., 2008; Moilanen et al., 2017). New and better models should be generated as more information on species or better environmental layers become available, and management decisions may have to be adjusted to the new results. However, the inevitable uncertainty associated to models should not be an obstacle to carefully using them to plan conservation action; in the presence of rapid environmental change the risks of inaction are 210 probably greater than those of judiciously using models (Wiens et al., 2009). The identified priority areas encompass ecologically distinct regions. For example, the A. ululata population that inhabits the mangroves of the Parnaíba river mouth (Mangrove), has a unique ecology (e.g. a diet composed of mangrove plants). The groups that inhabit the Enclave region of Ceará live mostly in humid enclaves with open ombrophilous forest. Adaptations to these different 215 environments may have resulted in populations of A. ululata with distinct behaviors, ecologies, gene pools and even morphologies. It is clearly desirable to protect the various ecological contexts in which the species is present. Virtually all priority areas with a protected status are in the northern part of the species range. The south, including the region with most occurrences, is completely unprotected. This worrisome situation is in line with the scarcity of protected areas in the Caatinga (de Marques & 220 Peres, 2015; DRYFLOR, 2016). It is important to note that, even on paper, the level of conservation provided by most protected areas in the range of A. ululata is very slight. Of the nine relevant protected areas only two are "full protection" units (Ubajara and Sete Cidades National Parks) having nature preservation as a main objective and only allowing the indirect use of natural resources (MMA, 2016). All other conservation 225 units allow sustainable use of natural resources (MMA, 2016), which needs to be well managed to 8 avoid damaging habitats. However, a serious lack of human and financial resources results in highly insufficient management. Moreover, in this region awareness of the protected areas is extremely low (Drummond et al., 2009). Our analysis revealed that deforestation is ongoing throughout most of the priority areas, even 230 inside protected areas (Fig. 4B), in line with the general trend of tree cover loss in the Caatinga (Beuchle et al., 2015). The situation is only better in federally protected National Parks, where we did not identify recent deforestation. It is thus evident that protected areas currently make a very small contribution to the protection of A. ululata, and of the other natural values of the Caatinga. Conclusions and Conservation consequences 235 The results of the spatial analyses performed in this paper suggest a variety of lines of action to preserve A. ululata. These actions vary across the range of the species, as shown in Fig. 4C, and are discussed here. Maxent modelling indicates that the range of A. ululata includes poorly surveyed areas where the presence of the species has not been confirmed. Moreover, the Zonation analysis shows that some 240 of these areas have a high level of conservation priority. It is thus urgent to survey these areas and despite the recent surveys by CPB/IBAMA, the absence of information is still flagrant (Fig. 4C). The Maxent analysis indicates that good tree cover and levels of aridity lower than those prevailing in the Caatinga region are critical for the maintenance of A. ululata, suggesting that the species is affected by both deforestation and climate change. Unfortunately, tree cover is declining 245 (Beuchle et al., 2015), even in protected areas (Fig. 4C), and arid conditions are increasing (Torres et al., 2017). Our results show that the legal protection of the areas suitable for the protection of A. ululata is very uneven. The central and southern areas of the species’ range, which includes some of the potentially best areas for its protection, is not included in any protected area. It is important to fill this 250 gap by designating new state, federal and private protected areas, especially in the larger contiguous priority areas (Fig. 4C). Moreover, it is critical to improve the management of the existing protected areas. Knowledge on the biology of A. ululata is still insufficient to plan efficient management measures. Therefore, in addition to survey work to clarify the status of the species, it is important to research aspects of its ecology that are critical to conservation. 255 The priority areas identified by the Zonation models mostly coincide with good remnants of Caatinga, which are also important for a substantial part of the great biodiversity of this biome, threatened by habitat destruction but poorly covered by protected areas (DRYFLOR, 2016). Due to 9 its large body size, the conservation of viable populations of A. ululata requires the maintenance of large areas of well-preserved habitat, making it a good umbrella species, with potential to contribute 260 to the protection of the rich biodiversity that shares its threatened habitat. Acknowledgements Study funded by Rufford Foundation grant 19646-1. RFF holds a FACEPE scholarship (IBPG-12362.05/16). We are thankful to collaborating rural communities, and to Thieres Pinto, Gabriela Linhares, Andressa Fraga and Nádia Freitas for assistance in the field. Marcos Fialho and Leandro Jerusalinsky 265 (CPB-ICMBio) kindly provided presence points of the species. Author contributions RFF and JMP conceived and designed the study and wrote the article. RFF performed the field work and analysed the data. References 270 Aguiar, J., Lacher, T.E. & da Silva, J.M.C. (2002) The Caatinga. In Wilderness: Earth’s last wild places (eds R.A. Mittermeier, C.G. Mittermeier, P. Robles Gil, J. Pilgrim, G.A.B. da Fonseca, T. Brooks & W.R. Konstant), pp. 174-181. Cemex, México, MX. Araújo, M.B., Williams, P.H., & Fuller, R.J. (2002) Dynamics of extinction and the selection of nature reserves. Proceedings of the Royal Society of London B: Biological Sciences, 269, 1971275 1980. Anadón, J.D., Giménez, A. & Ballestar, R. (2010) Linking local ecological knowledge and habitat modelling to predict absolute species abundance on large scales. Biodiversity and Conservation, 19, 1443-1454. Balk, D.L., Deichmann, U., Yetman, G., Pozzi F., Hay, S.I. & Nelson, A. (2006) Determining 280 global population distribution: methods, applications and data. Advances in Parasitology, 62, 119156. Bean, W.T., Stafford, R. & Brashares, J.S. (2011) The effects of small sample size and sample bias on threshold selection and accuracy assessment of species distribution models. Ecography, 35, 250-258. 285 TABLE 2 Maxent results. Variables Percent contribution Permutation importance Training gain without Training gain with only Test gain without Test gain with only AUC without AUC with only Percent tree cover 55.411 51.928 0.789 0.502 0.714 0.475 0.814 0.751 Bio17 15.786 16.139 1.009 0.145 0.909 0.149 0.846 0.661 Aridity Index 10.516 17.884 1.043 0.363 0.949 0.359 0.847 0.727 Roughness Index 7.102 1.778 1.025 0.140 0.966 0.120 0.854 0.622 Bio15 5.179 7.624 1.025 0.041 0.907 0.042 0.845 0.578 Forest Canopy Height 6.004 4.644 1.031 0.204 0.985 0.169 0.849 0.663 TABLE 3 Protected area and recent forest loss in each of the four priority regions mentioned in the text. Inclusion in protected areas Forest loss (2000 to 2014) Percentage (%) Km2 Percentage (%) Km2 Mangrove 69 2,359.028 8 272.640 Enclave 50 9,814.051 4 776.058 Caatinga 10 63.619 2 926.061 Border 3 15,519.234 6.7 1,039.788 Total of four regions 21 27,755.932 3.8 3,216.066 415 420 FIG. 1 Location of study area, surveyed regions and sites where the species has been recorded. 425 18 FIG. 2 Probability of occurrence of the Caatinga howler monkey, as predicted by the Maxent model. 430 FIG. 2 Probability of occurrence of the Caatinga howler monkey, as predicted by the Maxent model. 19 435 440 445 FIG. 3 Curves showing the relationship of each environmental variable with the probability of presence of the species. 20 450 FIG. 4 (A) Priority areas for the conservation of A. ululata, limits of the regions referred to in the text (Mangrove, Enclave, Caatinga and Border), and existing protected areas. (B) Forest loss between 2000 and 2014 within priority areas. (C) Areas requiring survey work, new protected areas, and improved management. 21 FIG. 4 (A) Priority areas for the conservation of A. ululata, limits of the regions referred to in the text (Mangrove, Enclave, Caatinga and Border), and existing protected areas. (B) Forest loss between 2000 and 2014 within priority areas. (C) Areas requiring survey work, new protected areas, and improved management. 455 22 FIG. S1 Relationship between percentage of presence points and the presence probability percentage 460 to define the suitability threshold.