Modeling ocelot (Leopardus pardalis) distribution in the southern limits in Brazil

ABSTRACT Ocelot is a Neotropical wildcat with a broad geographic distribution. The southernmost region of its distribution is an ecotone zone between biomes Atlantic Forest and Pampa in Brazil. Ocelot is rarely seen in this region, especially in the Pampa, being one of the most endangered wildcats in this region. The objective of this study is to predict the potential distribution of ocelots in the southern limit of it known distribution and to measure the percentage of protected areas coverage. We collected 24 records of the species in the state of Rio Grande do Sul and we modeled the potential distribution using seven distinct algorithms tested using two evaluation metrics (TSS and AUC). The variables with the highest percentage contribution to the models were: BIO8, BIO19, BIO7 (related to temperature and precipitation) and NDVI. Additionally, we observed that of the areas considered most suitable for the species (8,388 km2), only (39 km2 or 0.46%) were included in the protected areas and these were primarily located in the Atlantic Forest biome. By contrast, the most suitable areas for the species in the Pampa biome were not covered by protected areas. We consider this to be a serious issue considering the species is endangered in the state.


Introduction
The ocelot Leopardus pardalis (Linnaeus 1758) is a wildcat with a broad geographic distribution, occurring from the southernmost part of Arizona and isolated populations in Texas in the United States into northeast of Argentina and the southernmost parts of Brazil, in the state of Rio Grande do Sul (Paviolo et al. 2015). It has being found in tropical and subtropical forests, temperate forests, semitropical scrub and semi desert scrub (Martínez-Calderas et al. 2011). In the southern limits of the species distribution, a study carried out on Argentinean Atlantic Forest, estimated habitat models and abundance for the species (Cruz et al. 2019). However, the southernmost region of species distribution, an ecotone zone between biomes Atlantic Forest and Pampa that is located in Rio Grande do Sul state, in Brazil, there is no studies related to their habitat suitability, a region that has been highly disturbed by anthropogenic activities.
Populations that live in the limit of the species distribution area are among the most susceptible to a variable array of impacts, including being the frontline of the increasing of climate-related environmental change (Pearson et al. 2009;Hetem et al. 2014). As a carnivore, L. pardalis plays an important ecological role in the community (Fonseca & Robinson 1990), influencing not only the abundance of their prey populations, but also the dynamics of the plant community and, consequently, the vegetation diversity (Terborgh 1992). Therefore, they may be more directly affected by habitat fragmentation and habitat loss since they need to use large areas and may exacerbate the effects of climate change (Hetem et al. 2014). Thus, the population trends at the edge of the distribution could range between two extreme scenarios: populations could become completely lost, resulting in isolation of the species and a latitudinal displacement of a species range, or the species could find the changing climate beneficial, resulting in a simple expansion of the range into newly favorable regions (Hampe & Petit 2005). In this context, the conservation of natural environments must be maintained and an example of this is the Conservation Units (UCs), which serve to protect biological diversity and associated genetic resources (MMA 2002;Oliveira et al. 2010), which also can serve as protection for those populations living within the border areas of the species known distribution.
Historical data and current studies indicated that the southernmost occurrences of ocelot could extend south of Rio Grande do Sul state, into the Pampas biome and Uruguay (Ihering 1892;Araujo 1897;Indrusiak & Eizirik 2003;Andrade-Núñez & Aide 2010;Gonzáles & Lanfranco 2010;Queirolo 2016;Peters et al. 2017). Nonetheless, there is no information whether the Pampa is a suitable biome for the species to maintain resident populations. In addition, these two Brazilian biomes cited are threatened due to anthropogenic impacts (Ribeiro et al. 2009;MapBiomas Project 2019). The Atlantic Forest has been deforested down to a 12% of its original distribution (Ribeiro et al. 2009) and the Pampa lost 21% of its' native grasslands between 1985 and 2017 (MapBiomas Project 2019). Additionally, the Pampa is the least represented biome in the Brazilian National System of Conservation Units (SNUC), representing only 0.4% of Brazilian protected area (MMA 2000). Protected areas coverage is important for ocelots, since its populations are threatened mainly by loss and fragmentation of the natural habitats (Paviolo et al. 2015). Although the species is classified as Least Concern (LC) at the global level (IUCN 2018), it is regionaly categorized as Vulnerable to extinction (VU) in Rio Grande do Sul (FZB 2014). Taking all these together, it is clear that the species needs a special attention in this ecotone zone that Rio Grande do Sul state is located.
Species distribution models (SDMs) are important conservation tools, as they can serve as the basis for the acquisition of knowledge about rare species or populations, such as ocelot in Rio Grande do Sul. The use of SDMs that correlate observations of species occurrences with environmental variables is one of the most common approaches to estimating potential distributions of species (Guisan & Thuiller 2005;Dormann et al. 2012). In the last decade, this method has been used to predict suitable areas for wild cats (Wilting et al. 2010a;Rodríguez-Soto et al. 2011;Abade et al. 2014;Peterson et al. 2014;Cuyckens et al. 2016), helping researchers better understand the habitat requirements and thus the most important areas for these species (Rushton et al. 2004).
Improved the understanding of ocelot distribution in southern limit could help to identify ecological aspects and protected areas coverage for the edge populations the species. In this context, in this study we explored different potential species distribution model evaluation algorithms along with different evaluation metrics, to test the suitability of ocelot habitat in the two representative biomes at their southern limit of distribution and to compare the areas of suitability within Conservation Units (UCs) in Rio Grande do Sul. In addition, the results allow us to direct future conservation efforts for ocelot populations in poorly studied areas, as well as to base protected area Management Plans on the identified distribution areas. Our hypothesis is that the species has distinct areas of suitability in both biomes and our prediction is that the areas of highest suitability will be located in the forested areas of the Atlantic Forest once the species is typically related to dense vegetation (Harveson et al. 2004;Di Bitetti et al. 2008) and other studies report their densities and occupancies correlate with forest cover (Paviolo et al. 2015;Bolze 2019;Cruz et al. 2019;Wang et al. 2019).

Study area
The study area was focus in the southern limit of ocelot distribution, more specifically in the ecotone zone between Atlantic Forest and Pampa biomes and adjacent areas. The landscape composes a mosaic composed by a country matrix framed by riparian or hillside forest formations, the Atlantic plateau integrates the north of Rio Grande do Sul to the limit with the state of Santa Catarina. It is formed by an altitude steppe interspersed by mixed ombrophilous forest (northeast), deciduous seasonal forest (northwest and south portion) and dense ombrophilous forest (east slope) distributed up to 1,300 m in relation to sea level (Ab'Sáber 2003;IBGE 2004). The pampean steppe covers the southern half of Rio Grande do Sul and Uruguay. It includes both vast coastal country plains and mountainous regions thickened by seasonal semideciduous forest, which may slightly exceed 500 m (IBGE 2004) in elevation. Both areas have a humid mesothermal climate, with average temperatures ranging from 22°C during summer to 10°C during winter (Peel et al. 2007).

Points of occurence of the species
Occurrence data records were obtained through two distinct data sources: (i) primary data obtained from camera traps, direct visualization and identification of trails obtained by field expeditions headed by the authors between 2017 and 2019 and (ii) additionally, we considered records from bibliographic references (Jardim et al. 2005;Pires & Cademartori 2012;Queirolo 2016), whose records where obtained from interviews with local residents of the region. All species-related records were stored in a data sheet and the coordinates converted to geographic coordinates in decimal degrees and Datum WGS84.

Environmental predictors
We used 19 bioclimatic variables from WorldClim version 2.0 (http://www.worldclim.org), given that climate variables are among the most important factors influencing species distribution (Grinnell 1917;Guisan & Zimmermann 2000). In addition to bioclimates, we used altitude obtained from SRTM -Shuttle Radar Topography Mission (http://www2.jpl.nasa.gov/srtm), NDVI (Normalized Difference Vegetation Index) from MODIS -Moderate Resolution Imaging Spectroradiometer (https://search.earthdata.nasa.gov/ search), Distance to Roads and Distance to Water Resources, the latter were produced by the authors using the Euclidean Distance Tool of ArcGIS 10.3.5 software (ESRI 2019), the analysis was performed with maps of Water Resources and Highways downloaded from the site of 'Fundação Estadual de Proteção Ambiental Henrique Luiz Roessler' (http://www.fepam. rs.gov.br/biblioteca/geo/bases_geo.asp). All variables were standardized and processed using ArcGIS 10.3.5 software (ESRI 2019).
Collinearity between variables causes instability in parameter estimation in regression models (Dormann et al. 2013). To detect and remove highly correlated variables we used the VIF (Variation Inflation Factor) test using the 'vifstep' function of the 'usdm' package (Naimi & Araújo 2016) in the R environmental (R Core Team 2018). Here we consider higher VIF values >10 as adequate (Dormann et al. 2013;Naimi & Araújo 2016).

Potential distribution modeling
The search for the best suitability model should consider different modeling methods, since the existence of a single model explaining the distribution of a species is rarely ecologically reliable (Reichert 1997;Wintle et al. 2003), since the combinations of different approaches will always be more robust when faced with uncertainty (Araújo & New 2007;Gallien et al. 2012;Komac et al. 2016). To model the suitability of L. pardalis in Rio Grande do Sul, we used seven different algorithms: (i) Generalized Additive Models -GAM (Hastie & Tibshirani 1990); (ii) Boosted Regression Trees -BRT (Friedman 2001); (iii) Random Forests -RF (Breiman 2001); (iv) Generalized Linear Models -GLM (Mccullagh & Nelder 1989); (v) Domain (Carpenter et al. 1993); (vi) Support Vector Machine -SVM (Vapnik 1995) and (vii) Maxlike (Royle et al. 2012). For each model we used 70% of the points as training data and 30% of the occurrence points as test data, the replicates were elaborated with the bootstrapping method, for each algorithm 10 replicates were performed, in which the average of these samples was used as the result of each algorithm.
After the models were evaluated using true skill statistics, or TSS (Allouche et al. 2006) and the area under the ROC curve, the AUC (Fielding & Bell 1997). TSS is a threshold dependent accuracy measure, while AUC is a highly effective measure of the performance of ordinal scoring models and a threshold-independent precision measure (Thuiller et al. 2009;Lobo et al. 2008). TSS values ranged from −1 to +1, with −1 being predictions considered inaccurate and +1 considered accurate, TSS values >0.6 are considered excellent. AUC values range from 0 to 1, with 0 for predictions considered inaccurate and 1 considered systematically accurate model predictions, AUC values >0.7 are considered sufficient.
From the best algorithm, we calculate the relative importance of the different predictor variables based on the training datasets. Subsequently, we estimated the response curve of the suitability of the species in relation to the predictor variables that presented the highest relative importance. Finally, we make a raster map with a weighted average of all forecasts, specify which evaluation statistics (TSS or AUC) were used as weight in the weighted average procedure.

Results
We collected 24 records of Leopardus pardalis species in Rio Grande do Sul (Table 1), of which 15 occur in the Atlantic Forest biome and nine in the Pampa biome ( Figure 1). The results of the collinearity test resulted in 11 variables from the initial 23 (Supplementary Material, S1).
We modeled the potential distribution of the ocelot using seven different algorithms (Supplementary Material, S2), all of which performed reasonably to well, with TSS values ranging from 0.5 to 0.77 and AUC values from 0.73 to 0.92 ( Table 2). The potential distribution map of the species was constructed based on the mean AUC values and the mean TSS values (Figure 2). These results provide a satisfactory prediction for a potential distribution for L. pardalis. Among the variables of the best model, Random Forest algorithm, those that had a higher percentage of contribution (> 40%) to the ocelot potential model were Mean Temperature of Wettest Quarter (BIO8), Precipitation of Coldest Quarter (BIO19), Temperature Annual Range (BIO7) and NDVI -Normalized Difference Vegetation Index (Supplementary  Material, S3). These variables influenced the model in different ways (Figure 3). Of the total predicted area for species in the state (336,792 km 2 ), only (1,308 km 2 ) were in UCs. Regarding the areas considered most suitable for the species (8,388 km 2 ), only (39 km 2 ) were included in Conservation Units (Table 3).

Discussion
This study has helped to fill in some gaps about the knowledge of L. pardalis biology in the extreme south of Brazil, generating the prediction map of the species and observing how environmental variables influence its spatial distribution.
All the algorithms performed reasonably well, but the Random Forest (AUC: 0.92 and TSS: 0.77) performed better than the other. This method is a supervised learning algorithm; it consists of a large number of single decision trees that operate as an ensemble. Each model (tree) in the random forest is uncorrelated, but they operate as a committee, in which the most voted class spit out by each individual tree becomes the final model prediction of random forest, this results in a more accurate and stable prediction. Therefore, this algorithm could be more powerful than the traditional regression-based algorithms (Williams et al. 2009;Lei et al. 2011).
Among the variables of greater relative importance, those related to temperature indicated a negative tendency in relation to temperature. We can infer that, in the study region, L. pardalis seems to decrease its suitability with temperatures higher than 18°C. The suitability of the species in the region seems to be positively related to precipitation over 300 mm. Additionally, one of the variables that contributed most to the model was NDVI, it is a vegetation index that indicates the closer to 1, the denser the vegetation (Weier & Herring 2000). Our results indicated highest suitability among NDVI values of 0.8 and 1, indicating that this species is positively related to the dense vegetation areas in Rio Grande do Sul. When we look at the suitability map we observe highest suitability in the less degraded areas. These areas often have higher altitudes, being places of more difficult accessibility to human activities (such as agricultural production and road construction). According to Oberosler et al. (2017), human disturbance plays an important role in the probability of detecting rare species, indicating that the activity patterns of some mammals are totally dependent on the level of disturbance generated by humans, translating into a process of avoidance. This fact is applied to L. pardalis distribution in Rio Grande do Sul, as disturbed areas are avoided.
The regions with the best suitability for ocelot were estimated to occur mainly in the northeastern areas of the state, where part of the Atlantic Forest biome is located, corroborating our prediction, the species is typically related to dense forested vegetation     (2019), in the other extension of southern limit of species distribution, in Argentina. The authors observed ocelot has suitable habitat areas that correspond to native forest and that half of these suitable areas have some protection status. Our results also highlight the importance of protected areas for the conservation of the species in the southernmost limit of the Atlantic Forest. These forested protected areas are important because the species is characteristically related to dense vegetation (Harveson et al. 2004;Di Bitetti et al. 2008;Paviolo et al. 2015;Bolze 2019;Cruz et al. 2019;Wang et al. 2019). Unfortunately, areas of high suitability for the species in the Pampa biome (in the municipalities of Quevedos, Itaara and São Martinho da Serra) are not within protected areas or protected areas, which could be problematic for the individuals from these areas. Regional population estimates are restricted to the study by Kasper et al. (2015) and Bolze (2019), held at Turvo State Park, in the far north of Rio Grande do Sul. The authors assume that in a scenario of virtual total isolation of the park, the population of L. pardalis could reach about 0.11% of the national population, which is estimated to be approximately 40,000 mature individuals (Oliveira et al. 2013). However, local reality shows that park is connected to the Green Corridor of Misiones, which makes it possible to assume a representative population of up to 7.5% compared to the national estimate. Such a finding points to the importance of intact and connected areas (either by the protection of the UCs system or the difficulty of human access) to the regional viability of the species, a fact present along the Atlantic Forest and absent along the Pampa.
From the above considerations, it is possible to assume that the north-northeast of Rio Grande do Sul gathers attributes to be considered important to the regional population of L. pardalis. Dense protected forests, difficult to reach and connected to each other, could be considered as a regional population source area, supporting the occurrence of a viable population and allowing dispersions toward rural and fragmented areas located at the southern end of their range.
The opposite situation is verified for the interior of the Pampa biome, which was predicted to have very low to medium suitability for the presence of L. pardalis. The greater availability of rural environments, low forest cover, the flat relief and land favorable to human activities, as well as the smaller representativeness of protected areas reflects the rarity of records of the species toward the south of Rio Grande do Sul. The recent quotation by Andrade-Núñez and Aide (2010) for the Uruguayan Pampa still needs confirmation, however, the historical occurrence pointed out by Queirolo (2016) confirms the past occurrence of the species in the neighboring country. Already the isolated record recently confirmed by Peters et al. (2017) at the southern tip of the Brazilian Pampa corroborates the existence of individuals on the southern tip of the Pampa. The absence of connected UCs, the low availability of forest formations and the increasing human occupation along the Pampa making difficult to establish L. pardalis population in the Pampa, so that the few occurrences can be considered as a result of occasional dispersions. Such conditions imply the great importance in the establishment of protected and connected areas along the Pampa for the maintenance of a locally viable populations and avoidances of isolated presences.
Finally, the ecotone zone between the two biomes were predicted to medium to high suitability for the species (especially in the municipalities of Quevedos, Itaára and São Martinho da Serra), this region has less geographical representation and is dominated by certain anthropogenic activities, however, it is characterized by an irregular pattern of agricultural land and human settlements. Despite the absence of UCs, the presence of riparian forest corridors allows the dispersal of species between the Atlantic plateau and the Pampean plains, which also tends to favor the full population establishment of L. pardalis in the region.
This study contributes relevant information to early conservation strategies for L. pardalis in Rio Grande do Sul, identifying new records and describing areas of greater population suitability along its southern distribution area. Thus, it is emphasized the need to carry out studies of greater temporal scope to evaluate abundance, density and spatial ecology of L. pardalis at the southern limit of it population and dispersal range.