Effects of landscape modification on species richness patterns of fruit‐feeding butterflies in Brazilian Atlantic Forest

To assess the distributional patterns of fruit‐feeding butterfly species richness in Atlantic Forest (AF) based on stack species distribution models (SSDM); to evaluate the relative contribution of climate and landscape in the patterns of butterfly species richness; and to recommend conservation guidelines for AF regions based on the obtained results.


| INTRODUC TI ON
Extensive human-induced activities have been altering the dynamics and functioning of ecosystems over the last decades, affecting both biodiversity and societal wellbeing Hooper et al., 2012). Recent findings indicated that anthropogenic climate disruption and modification of natural habitats are the main drivers of species loss (Dirzo et al., 2014). Therefore, understanding biodiversity distribution patterns is an important step towards identifying priority areas and proposing more effective conservation policies to decrease the current biodiversity crisis (Margules & Pressey, 2000;Myers, Mittermeier, Mittermeier, Fonseca, & Kent, 2000).
Mapping areas of high species richness is one of the proxies used for the establishment of priority areas for conservation (Ceballos & Erhlich, 2006;Jenkins, Pimm, & Joppa, 2013). The success of this approach is based on the premise that conserving the richer areas necessarily implies in preserving more species.
However, conservation networks built under this criterion may not cover the species with small-range distributions (Nobrega & de Marco, 2011;Veach, Di Minin, Pouzols, & Moilanen, 2017). Therefore, both species richness and endemism are required to estimate the conservation value of a given area (Faith, Reid, & Hunter, 2004;Kier et al., 2009). Advanced macroecological techniques are increasingly able to fill the gaps in biodiversity knowledge (Hortal et al., 2015;Nobrega & de Marco, 2011;Siqueira, Durigan, de Marco, & Peterson, 2009). Thus, the use of species richness and endemism allied to the niche modelling techniques may represent an effective approach in conservation protocols (Costa, Nogueira, Machado, & Colli, 2010).
Climate and the barriers of dispersion are considered the main structuring factors of communities on a regional scale, from a hierarchical perspective of importance (Cornell & Harrison, 2014;Willig, Kaufman, & Stevens, 2003). Based on a species' physiological responses, one can identify the combination of climatic variables (i.e. the bioclimate envelope) that determine its distribution and then predict its potential distribution at broad-spatial scales (Guisan & Zimmermann, 2000;Pearson & Dawson, 2003). Conversely, at narrow-spatial scales, land use processes can be more important in determining species' distributions (Pearson & Dawson, 2003). A framework considering this hierarchical role of broad and narrow processes on the moulding of macroecological patterns of biodiversity has already been proposed, but has been rarely tested (Guisan & Rahbek, 2011). Therefore, tests involving the effects of narrow-scale processes on macroecological patterns are important, especially in highly human-impacted areas.
The Atlantic Forest (hereafter AF) domain occurs along a long latitudinal gradient of the Brazilian coast, extending west inland to areas of Paraguay and Argentina (Morellato & Haddad, 2000; Aguiar, Ribeiro, Metzger, & Peres, 2010). The AF suffered for centuries from the overexploitation of its natural resources and human land occupation (Laurance, 2009), causing large-scale landscape modification and forest fragmentation (Ribeiro, Metzger, Martensen, Ponzoni, & Hirota, 2009). These threats, coupled with high rates of species endemism, make the AF a biodiversity hotspot and one of the three biomes that are most vulnerable to global change (Bellard et al., 2014;Myers et al., 2000). However, the effect of landscape modification on AF biodiversity in a macroecological perspective still requires investigation.
Historical climate and geological processes are traditionally considered the drivers of the high rates of endemism and diversity in the AF (Sobral-Souza & Lima-Ribeiro, 2017). One of the most well known hypotheses, for example, links diversification rates to the regions that remained stable through the historical climatic fluctuations and served as refuges for forest species (Carnaval & Moritz, 2008;Graham, Moritz, & Williams, 2006). On the other hand, recent processes, such as the wide-scale conversion of native forest into human-modified landscapes, have been leading the AF to alarming rates of habitat loss (Tabarelli et al., 2010). The effectiveness of conservation efforts depends on a clear understanding of both (a) the contribution of climate to species distribution, and (b) species' responses to human-modified landscapes. Notwithstanding, some variables related to landscape ecology, such as forest cover, land use and habitat change, are generally neglected in macroecological modelling processes (Carnaval & Moritz, 2008;Carnaval et al., 2014;Metzger, 2001).
Here, we modelled the distribution of fruit-feeding butterfly species using both climate variables and landscape metrics, to understand the patterns of species richness distribution in the AF.
Tropical butterflies are considered good models for ecological studies due to their sensitivity to environmental changes (Bonebrake, Ponisio, Boggs, & Ehrlich, 2010). More specifically, the guild of fruit-feeding butterflies, whose adults primarily feed on rooting fruits, are employed as model organisms for environmental monitoring and conservation (Barlow, Overal, Araujo, Gardner, & Peres, 2007;Freitas et al., 2014;Santos, Marini-Filho, Freitas, & Uehara-Prado, 2016). Given the dependence on fruit and carrion resources (Horner-Devine, Daily, Ehrlich, & Boggs, 2003), the availability of natural forest habitats is essential for the maintenance of fruit-feeding butterfly species. Therefore, we hypothesized that the regions in which the large Atlantic Forest remnants are located would exhibit higher fruit-feeding butterfly species richness. Based on the above scenario, here we addressed the following questions: (a) How is the pattern of fruit-feeding butterfly richness distribution in the AF? (b) Are the patterns of species richness distribution predicted by climate conditions and landscape metrics congruent? Finally, we discuss conservation issues based on the obtained results.

K E Y W O R D S
Atlantic Forest, butterflies, conservation, diversity, macroecology, species richness 2 | ME THODS

| Study area
The Atlantic Forest (AF) roughly comprises two major vegetation types: the coastal Atlantic Rainforest and the Atlantic Semideciduous Forest (Morellato & Haddad, 2000). Based on a non-objective analysis of three biological groups distribution (birds, primates and swallowtail butterflies), Silva and Casteleti (2003) (Figure 1). This scheme was followed by Ribeiro et al. (2009), and several other subsequent authors since then. Currently, 70% of the Brazilian population are settled within the AF domain and the natural forest cover is confined to approximately 11% of its original extent (Ribeiro et al., 2009;Tabarelli et al., 2010). In this study, we followed the Atlantic Forest delimitation proposed by Muylaert et al. (2018).

| Species database
We used a large dataset of fruit-feeding butterfly communities from the Atlantic Forest, which consists of a compilation of published and unpublished records from 1949 to 2018. It comprises 7,062 occurrence records of 279 species, distributed across 122 well-sampled local communities. The association of historical species records with landscape metrics in niche modelling analysis is potentially problematic. As the landscape changes in time, their current metrics may not represent an equivalent niche from the time of the species records. However, most species records in our dataset pertain to local communities sampled during the last 20 years, which endorses synchrony with the landscape metrics ( Figure S1.1). More details about the dataset can be found in Santos et al. (2018). We modelled the distribution of the 146 fruit-feeding butterfly species represented by at least 10 occurrence points. The remaining 133 species (with <10 occurrence points) were considered here as endemic or rare species and were counted only in the grid cells of occurrence and then summed to the final maps of richness.

| Climate and landscape variables
First, we built species-specific distribution models (SDMs) for landscape and climate variables separately. We separated these models to infer differently scaled processes driving species distributions (narrow and broad-scale ecological filters processes, respectively).
Climate is known to be a driver of species distributions at broad-scale (Soberón & Nakamura, 2009) and landscape modification through habitat fragmentation is associated with major biodiversity changes, including species loss and local persistence at narrow-scale (Haddad et al., 2015;Thompson, Rayfield, & Gonzalez, 2017). Although we consider the hierarchical effects of climate and landscape, we constructed our SDM models with both variables in the cell-size resolution of 1 km 2 . We chose this resolution based on the differential occupation of butterfly species on ecological gradients in response to microclimatic variations. In addition, the Atlantic Forest terrain is subject to abrupt changes in altitude and slope due to its highly heterogeneous topography. Thus, we maintained the same cell-size resolution for landscape variables to catch the climatic variation within these different slope categories and directly compare with climate effects.
To build climate-based models, we used 19 bioclimatic variables (WorldClim version 2, http://world clim.org/version2, Fick & Hijmans, 2017). We used principal component analysis (PCA) to convert all climate variables into linear uncorrelated variables. For landscape-based models, we used five variables: forest cover percentage, patch size, edge distance, slope and functional connectivity-assuming 180 m as the limit of dispersal (following Ribeiro, Batista, Prado, Brown, & Freitas, 2012). These variables were selected based on the premise that most fruit-feeding butterfly species are adapted to forested habitats and depend on fruit resources. All these variables were obtained from the Spatial Ecology and Conservation Laboratory at São Paulo State University (LEEC-UNESP) (M. C. Ribeiro pers. comm.). To remove collinearity between landscape variables, we applied the same PCA procedure employed for climate variables to the landscape variables. We selected the set of PCA-derived variables for both climate and landscape based on Kaiser-Guttman rule (retention of axes with eigenvalues higher than 1) as proposed in de de Marco and Nóbrega (2018). Therefore, the first four axes of the climate and landscape ordinations were selected as variables, representing more than 90% of the cumulative proportion of variance in both cases (Table 1, see Figure S1.2). The PCA analyses were conducted in "RStoolbox" R-package (Leutner, Horning, & Schwalb-Willmann, 2018).

| Predicting fruit-feeding butterfly richness in Atlantic Forest
We used ecological niche modelling (ENM) techniques to predict geographical patterns of fruit-feeding butterfly species richness in the AF. ENMs infer relationships between environmental variables and species occurrence points to predict the suitability of habitats in unknown sites (Peterson et al., 2011). The ENM procedure consists of using three information sources: (a) known occurrences of a given species, (b) environmental layers and (c) mathematical algorithms (Franklin, 2009;Peterson et al., 2011).
To evaluate the resulting models, the species occurrence points were partitioned into two distinct subsets, 75% and 25% for model training and evaluation, respectively. As these subsets belong to the same data set, we repeated this procedure 10 times using the k-folding technique (k = 2) to reduce the collinearity. In total, we generated 40 predictions for climate-based models and another 40 for landscape-based models separately for each species (10 randomizations × 4 algorithms). To transform the model outputs into binary maps we calculated threshold values using maximum sensitivity and specificity. These thresholds maximize the correctness of presences and absences and have been shown to be effective in predicting occurrences from presence-only models (Liu, Newell, & White, 2016).
After defining the thresholds, we conducted the ensemble forecasting technique to obtain the prediction map of each species distribution (Araújo & New, 2007). The maps were generated based on climate and landscape separately. We produced the binary maps belonging to each algorithm (replicates) using the respective threshold values previously calculated and then summed the maps of the same algorithm and between the algorithms. To evaluate each of generated models, we estimated the values of True Skilled Statistic (TSS). The TSS values range from −1 to 1, where negative or near-zero values indicate that predictions do not differ from randomly generated models, whereas models with values closer to 1 indicate that observed and modelled distributions are in near-perfect agreement (Allouche, Tsoar, & Kadmon, 2006).
Finally, to create the fruit-feeding butterfly richness map we inferred the lowest presence threshold (LPT) for each species prediction map (Pearson, Raxworthy, Nakamura, & Peterson, 2007) F I G U R E 1 Distribution of the Atlantic Forest domain (following Muylaert et al., 2018), its biogeographical sub-regions (adapted from Silva & Casteleti, 2003), and the forest remnants (Ribeiro et al., 2009) to transform each continuous frequencies into binary maps (0 for absence and 1 for presence). Overlapping all species binary maps, we obtained the predicted number of species occurrence (richness) of fruit-feeding butterflies in each cell of the AF delimitation. Guisan and Rahbek (2011) described this method as stack species distribution models (SSDM). All modelling steps described above were followed separately to build a climate-based and a landscape-based map, thus resulting in one map of climate-richness and another map of landscape richness (Figure 2a,b).

| Predicting conservation issues in AF butterfly richness
To evaluate the contribution of climate conditions and landscape metrics in richness patterns and propose conservation efforts, we conducted an EcoLand analysis (adapted from Ferro e . We extracted the richness values from climate-based models and landscape-based maps and generated a scatter plot with the Xaxis representing climate-based richness and the Y-axis representing TA B L E 1 The first four axes of both Principal Component Analysis (PCA) of climate (19 variables extracted from Fick & Hijmans, 2017) and landscape (five variables extracted from the Spatial Ecology and Conservation Laboratory, Ribeiro pers. comm.), and their respective scores. (SD = Standard deviation)  (Figure 2c). The EcoLand map was used for designing conservation plans according to the richness scenarios.

| Climate-based models
According to the climate-based model, areas with high-predicted species richness have up to 162 species, and the areas with lowest richness have as few as 10 fruit-feeding species. Sites with high species richness are mostly concentrated in the south-eastern portion of AF, which include the Serra do Mar and Serra da Mantiqueira complexes, denominated as "Serra do Mar sub-region" by Silva and Casteleti (2003) (Figure 1). Within this sub-region, these areas comprise a narrow southern coastal strip that extends north-eastwards until the Doce river basin (Figure 3). The areas with lowest species richness are located within the São Francisco river basin, which is labelled the "São Francisco transitional sub-region" in Silva and Casteleti (2003)

| Landscape-based models
The landscape-based model predicted up to 190 fruit-feeding but-

| Predicting conservation issues in AF butterfly richness
According to the EcoLand models, contiguous areas with high species richness predicted by both climate and landscape-based models are found in the complex of Serra do Mar and Serra da Mantiqueira mountain ranges. There are also small patches of high climate and landscape richness in southern AF (Araucaria subregion), and in some north-eastern areas in Bahia, Pernambuco and interior sub-regions (red areas in Figure 5). Medium values of species richness predicted by both climate and landscape models were the most frequent cells in EcoLand analysis (light green areas in Figure 5). These areas are concentrated in one large patch in the mid-western AF (Interior sub-region) and two smaller patches in the north-eastern AF (Bahia and Pernambuco sub-regions).
Our models showed that most areas in the AF have high species richness predicted by the landscape-based models but not for cli-

mate-based models, being the second most frequent cells in the
EcoLand analysis (orange areas in Figure 5). These areas include the Araucaria (Paraná and Santa Catarina states) and the Bahia sub-regions. Areas with high richness according to climate-based models but low richness according to landscape-based models are mostly found adjacent to the high richness areas predicted by both models. These are probably the areas where butterfly richness was most affected by landscape modification and forest fragmentation. These areas are spread through practically all biogeographical sub-regions (dark green area in Figure 5), except the São Francisco sub-region. Finally, areas with low species richness predicted by at least one of our models (light blue) were found in the north-western region of AF, predominantly in the biogeographical sub-region of São Francisco basin.

| D ISCUSS I ON
Here, we modelled the species richness of Atlantic Forest fruit-feeding butterflies using both climate variables and landscape metrics.
Our findings demonstrate that the Serra do Mar sub-region is characterized by continuous high species richness areas, independently predicted by both climate and landscape models. Considering the combined effect of climate and landscape models in the EcoLand map, the high richness areas coincided with the location of the continuous forest remnants (Figures 1 and 5). Moreover, most of these species-rich spots are close to areas whose landscape-predicted richness values are lower than expected according to the climatepredicted richness (dark green areas in Figure 5). These findings clearly indicate the effects of landscape modification in the species richness patterns of fruit-feeding butterflies and the importance of the forest integrity for the maintenance of high butterfly richness indexes in the AF.
Forest cover is critical to the maintenance of specific microclimate conditions, the importance of which increases in seasonal areas (Horner-Devine et al., 2003). Some clades of fruit-feeding species are strictly associated with forest-shaded habitats and forage in specific forest vertical strata (DeVries, Murray, & Lande, 1997;Santos, Iserhard, Carreira, & Freitas, 2017). Fleshy fruits and organic matter, whose accumulation depends on the presence of forest cover, are important food sources for the adults of fruit-feeding species. Moreover, forest habitats represent good sources of host plants for their immature stages (Horner-Devine et al., 2003;Koh, 2007). The presence of forest fragments (source areas) and putative forest corridors are also critical for promoting the connectivity and species dispersion between areas, improving the persistence of populations of many forest species (Brown & Freitas, 2002;Lees & Peres, 2008).
The regions containing high species richness predicted by both climate and landscape models are within the Serra do Mar, Bahia and Pernambuco sub-regions (Figures 1 and 5). These sub-regions coincide with three out of the four regions recognized as butterfly endemism centres in the AF, according to Brown's (1979) original propositions based on the distribution of geographic races (butterfly subspecies). Moreover, most records of endemic or rare species in the database lie within the areas with high species richness (see F I G U R E 2 Framework demonstrating the steps of EcoLand analysis. The first step consisted of modelling the distribution of each butterfly species according to climate variables, using different algorithms. The ensemble method overlaps the maps from all algorithms and constructs a final consensus map for each species. From threshold values, the species suitability is converted in presence of species at each cell. These presences are summed to generate the distribution map (a). The same procedure is performed to landscape modelling, using landscape variables (b). Finally, the species richness predicted by the landscape in each cell was plotted against richness predicted by climate. The relation between the richness values was categorized. In the EcoLand map, these cell categories were identified in geographical space, thus allowing the understanding of the relative contribution of each variable in the species richness distribution (c) Figure S1.3). As the occurrence of these species corresponds to the precise record location (in other words, not modelled), they increase the richness values only in that specific cell. Therefore, the restricted distribution of these species might explain the overlap between endemism centres and areas of high species richness.
The SSDM techniques are expected to present upper values for species richness due to the uncertainty about the dispersal capacity of each taxon and energy availability, but principally because they do not consider other community constraints such as the biotic interactions (Guisan & Rahbek, 2011). For example, our landscape-based models predicted that at least 20 fruit-feeding butterfly species are present at any AF locality. Therefore, one could assume that AF localities with <20 species are possibly under-sampled or are naturally poor in fruit-feeding species, such as the high altitude grasslands and the pampas (e.g. Carvalho, Piovesan, & Morais, 2015;Lemes, Carvalho, Ribeiro, & Morais, 2015;Morais, Lemes, & Ritter, 2012). However, our model is predicting the maximum species richness value based on landscape niche properties, without considering that communities may also be subject to constraints of other assembly processes, as the case of biological interactions. The non-incorporation of other assembly processes does not invalidate our general patterns, but requires attention when interpreting them at the community level. As our models supposedly overpredict the richness values, it means they are not necessarily the same as the observed species richness.

| IMPLI C ATI ON S FOR CON S ERVATI ON
High human disturbance and agricultural expansion are continuously replacing natural habitats with more simplified and homogeneous landscapes (Hobbs, Hallet, Ehrlich, & Mooney, 2011). Habitat change may lead to generalist species replacing specialist species (i.e. biotic homogenization), with consequences for ecosystem functioning (Clavel, Julliard, & Devictor, 2011;Filgueiras et al., 2019;Olden, Poff, Douglas, Douglas, & Fausch, 2004). According to the "insurance hypothesis", a higher diversity of functional groups provides greater stability and resilience of ecological processes, thus buffering ecosystem functioning against environmental changes (Tscharntke et al., 2012;Yachi & Loreau, 1999). Decreasing species richness induced by landscape modification would reduce functionality, which in turn F I G U R E 3 Distribution of fruit-feeding butterfly species richness in Atlantic Forest according to the climate model. Warmer colours indicate a higher number of predicted species F I G U R E 4 Distribution of fruit-feeding butterfly richness in Atlantic Forest according to the landscape model. The lighter colour regions represent the landscapes with a higher number of predicted species compromise the buffering effect of communities (Gámez-Virués et al., 2015). However, in ecosystems with high species richness, species with low and high productivity are expected to coexist, thus not necessarily affecting the ecosystem productivity due to the asynchrony of their responses (Yachi & Loreau, 1999).
Following the premise that higher species richness reduces the probability of compromising the ecosystem functioning, we can use our EcoLand results as a proxy for proposing different conservation strategies for the Atlantic Forest. For instance, areas containing high species richness predicted by both climate and landscape models would be the core places for conservation, as the case of the Serra do Mar sub-region ( Figure 6). Indeed, this sub-region harbours the largest Atlantic continuous forest remnant, which covers high numbers of species richness and unique species ( Figure S1.3). On the other hand, it also contains the largest Brazilian urban areas, resulting in high pressure for the adjacent natural habitats (CNUC, 2018;Rezende et al., 2018).
In the regions where landscape structure predicted richness values below than expected by climate, the possibility of habitat restoration should be prioritized (dark green areas in Figures 5 and   6). Long-term restoration plans can expand suitable habitats for the maintenance of high local butterfly diversity (Sant'Anna, Ribeiro, Garcia, & Freitas, 2014) and, consequently, other biological groups.
However, forest restoration would be less effective and highly expensive in areas that harbour dense urban centres in the AF domain.
Therefore, a careful evaluation of any restoration plan is necessary, taking into account the resilience of the system, with the aim of optimizing strategies in terms of efficiency and cost (Tambosi, Martensen, Ribeiro, & Metzger, 2014). Nevertheless, we propose here the use of habitat restoration to increase local and regional species richness only; the reestablishment of other ecological functions and processes at narrow scales requires further studies before one could take decisions about the best restoration strategy for a given area.
In a scenario of climate change, some landscapes currently suitable for maintaining a high species richness may display a more suitable climate for harbouring a high species richness in the future. The Araucaria and São Francisco sub-regions (orange and yellow cells in Figures 5 and 6) are two examples that could fit the above scenario.
Conversely, in the Bahia sub-region, forests have been devastated in the last decades and now persisting as a plethora of small fragments (Ribeiro et al., 2009). Therefore, in this region, the implementation of ecological corridors would improve the connectivity between the fragments and of these with other small hotspots of richness (Dunwiddie et al., 2009).
Finally, in those regions where both landscape and climate-based models predicted low species richness, the perspectives of richness maintenance through restoration or conservation practices are F I G U R E 5 Fruit-feeding butterfly richness distribution in Atlantic Forest according to the EcoLand model, which takes into account the climate and landscape variables. The colour scheme indicates the different combinations between climate and landscape richness predictions F I G U R E 6 EcoLand colour scheme indicating different combinations between climate and landscape richness predictions and the recommended conservation and study approaches equally low. In addition to the inauspicious climate for supporting a high local or regional species richness, the landscapes are so devastated that the conditions would hinder the dispersal and maintenance of viable populations of many forest-related biological groups.
In short, these localities represent expensive restoration efforts, although they are locally important for conserving the most basic ecosystem services, such as pollination, water flow regulation and the reduction of greenhouse gases emission (Tambosi et al., 2014).

| CON CLUS IONS
Here, we aimed to unravel the relative contribution of climate and landscape variables to the distribution patterns of fruit-feeding butterfly species richness in the Atlantic Forest. Our results highlight the complementarity of different processes in shaping the patterns of species richness and the negative impact of human land use on natural habitats as biodiversity sources. Future studies addressing other biodiversity aspects (e.g. functional and phylogenetic diversity and beta diversity) are necessary to unravel the effect of humaninduced processes on the AF biodiversity and trace out the best conservation strategies in this vanishing biome.

ACK N OWLED G EM ENTS
We thank the following institutions for financial support:

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.