Assessing Suitable Habitat and Functional Connectivity for Margay ( Leopardus wiedii ) in Central Western Forests in Michoacán,México

Introduction

Habitat loss and fragmentation threaten biodiversity by reducing species’ distribution and affecting natural populations’ dispersal and connectivity (Beier et al., 2007; Imbernon et al., 2005). In the face of these threats, several strategies have been implemented to identify areas that can sustain viable populations of endangered species and plan effective conservation strategies. Habitat suitability analysis is a methodology used to evaluate and predict the quality and availability of suitable habitats for a particular species or group of species (Wang et al., 2022). This type of analysis is essential in biodiversity conservation as it assists in identifying areas with the potential to sustain viable populations of endangered species, and improve effective conservation strategies. Habitat suitability analysis is a methodology used to predict and assess the quality and availability of suitable habitats for any particular species (Hirzel & Le Lay, 2008). This kind of analysis is essential for biodiversity conservation as it helps to identify areas where viable populations of endangered species may survive and to elaborate effective conservation strategies (Muhammed et al., 2022).

A detailed understanding of the habitat characteristics required by different species may facilitate the design and zoning of nature reserves and biological corridors that maintain the continuity of contiguous stands of suitable habitat, facilitating movement and genetic flow between populations. Maintaining habitat continuity at a landscape level is a strategy to mitigate some of the worst detrimental effects of fragmentation (Gilbert-Norton et al., 2010; Resasco, 2019). Landscape connectivity implies the capacity of a landscape to maintain vegetation structures that allow the movement and dispersal of species across its habitats. Connectivity has two main approaches: structural connectivity, which describes the configuration of habitat patches in the landscape, and functional connectivity; defined as a behavioral response of organisms to landscape structural and physical configuration (Adriaensen et al., 2003; Tischendorf & Fahrig, 2000). Worldwide, avoiding habitat fragmentation and establishing and maintaining connectivity among remaining natural habitats and protected areas has been used in the conservation of large carnivores, including felids, which in most species require large habitat areas and adequate prey abundance and shelter due to their territorial behavior (Sunquist & Sunquist, 2002; Wallace et al., 2010).

For example, in a continental level diagnostic for jaguar (Panthera onca) conservation, two ecoregions in Michoacán: the Jalisco dry forests and Sierra Madre del Sur pine-oak forests were identified as potential biological corridors (Rabinowitz & Zeller, 2010). The forest cover of these regions is assumed to allow the suspected dispersal and gene flow for the jaguar population occurring on the Mexican Pacific slope, estimated as the largest in the country (Charre-Medellín et al., 2015, 2021; De la Torre et al., 2017; Rabinowitz & Zeller, 2010; Rodríguez-Soto et al., 2013). Other endangered or threatened felines such as ocelot (Leopardus pardalis), margay (Leopardus wiedii), and jaguarundi (Herpailurus yagouaourundi) also inhabit the tropical and temperate forests of the state (Charre-Medellín et al., 2015; Monterrubio-Rico et al., 2012, 2017). However, the status of the regional populations is unknown due to the lack of systematic feline surveys (Charre-Medellín, 2009; 2015; Monterrubio-Rico et al., 2012, 2017).

The margay is a particularly vulnerable species, listed by the IUCN as “Near Threatened” globally, as margay populations are in decline through its range due to conversion of forest habitats to agricultural areas and human infrastructure interference. In addition, it’s specialized aboral lifestyle presents a higher risk of extinction. The margays apparent ecological displacement by ocelots due to competition is also significant, which are sympatric in many areas (Oliveira et al., 2010, 2015). Although the species presents a large extent of occurrence at continental level, its area of occupancy is considerably smaller. In addition, the margay is one of the least studied felines as camera trap surveys are more likely dedicated to assess jaguars and ocelots (Brodie, 2009; de la Torre et al., 2017; Nagy-Reis et al., 2020; Tobler & Powell, 2013).

In Mexico, the margay is listed as endangered by the Mexican official protected species law known as the NOM-059-SEMARNAT-2010 (SEMARNAT, 2010). Although the knowledge of this species has expanded over the last decade, population assessments providing updated estimates of abundance or current distribution within its broad range in the country are scarce. In the central western state of Michoacán, recent wildlife assessments based on camera trap surveys revealed the presence of this species in tropical and temperate forests. Records were obtained mainly for the Sierra Madre del Sur pine-oak forests, the Jalisco dry forests, and to a lesser extent, in the Trans-Mexican volcanic belt pine-oak forests and the Balsas dry forest ecoregions (Charre-Medellín et al., 2015; Monterrubio-Rico et al., 2023).

However, forest loss is increasing in the state of Michoacán and fragmentation is reducing the size and connectivity of forest remnants (CONAFOR, 2020). The margay is vulnerable to the fragmentation and the overall reduction of forested areas in Michoacán state (Carvajal-Villarreal et al., 2012; Charre-Medellín et al., 2015, 2021; Rodríguez-Soto et al., 2013).

Although some studies on margay provide regional, national, or continental distribution estimations, these are inadequate to evaluate a species status in scenarios of heavy regional forest fragmentation, which is occurring at a fast pace, making updated and improved assessments necessary. To complement information on the species’ potential habitat availability, models should incorporate the species’ likely mobility and connectivity among potential populations affected by habitat encroachment (Espinosa et al., 2017; Martínez-Calderas et al., 2016; Monroy-Vilchis et al., 2019). Habitat suitability models represent an adequate approach to indirectly evaluate a species’ status, as the models correlate the environmental variables and the presence of the species to determine favorable habitat conditions for the species in the geographic space (Soberón et al., 2017). These models are the basis for subsequently evaluating functional connectivity by facilitating the identification of suitable habitat fragments for the species of interest (Alexander et al., 2016; Almasieh et al., 2016; Correa-Ayram et al., 2014; Rodríguez-Soto et al., 2013). As in the case of the jaguar, the forest fragments of Michoacan may be critical to maintaining margay populations and their mobility. Large forest stands may likely harbor populations with estimated densities varying from 17 to 26 per 100 km² (Cuellar et al., 2006). The species home ranges vary from 4.0 - 21.8 km² (Carvajal-Villarreal et al., 2012; Kasper et al., 2016).

Functional connectivity of felids has been evaluated using “Least Cost Paths” (LCP) and circuit theory; both operate from a similar logic of landscape resistance. This resistance is understood as a differential risk for an individual by the limits imposed on crossing an unsuitable or hostile area or habitat or the energy spent on moving through a specific element of the landscape (de la Torre et al., 2017; Espinosa et al., 2017; Ávila-Coria et al., 2015). “Cost values” are weighted to resistance; if the resistance or difficulty is high, the cost value will be high (Zeller et al., 2012). The LCP method estimates the cost of potential routes for an organism to cross from one habitat patch to another and selects the one with the least resistance. In contrast, the circuit theory method identifies the paths with higher connectivity probability. This method weights all the potential routes and displays them as connectivity probability values. Circuits are conceived as network nodes connected as electrical components that conduct current flows (McRae et al., 2008).

A few connectivity studies have assessed Felidae in México, including for jaguars the circuit theory (de la Torre et al., 2017) or least-cost paths (Duenas-López et al., 2015; Salazar et al., 2017; Van den Berg, 2021). Connectivity analysis with “least-cost paths” have been also implemented to evaluate ocelot (Leopardus pardalis) in Oaxaca, México (Van den Berg, 2021), puma (Puma concolor) in central region of Mexico (Aguascalientes) (González-Saucedo, 2011), as well as for bobcat (Lynx rufus) in the central region of México (Michoacan) (Correa-Ayram et al., 2014).

Considering that the margay is the most affected species by forest loss and fragmentation in the central location of the study region, it is essential to evaluate the species’ habitat availability and functional connectivity. Such information assists in designing complementary areas for protection or evaluating the degree of geographic isolation of established protected areas based on the specific requirements of the threatened species. In Mexico, establishing protected areas to improve and complement the protection of biodiversity and environmental services by state or municipal governments is mandatory (Diario Oficial de la Federación, 2014). Temperate and tropical forest remnants in the region may be assumed to function as a connective link for margay populations in western and central Mexico. Therefore, evaluating the degree to which the landscape maintains suitable conditions to provide functional connectivity for the margay is necessary. The aims of this study were: i) to determine the distribution of suitable habitat conditions for the margay in the state of Michoacán, México, and ii) to evaluate the study region’s stability for the species’ distribution and the habitat’s functional connectivity.

Methods

The state of Michoacán is located in the central-western region of México between the external coordinates 20°23′40″, 17°54′54″ N, and 100°03′47″, 103°44′17″ W (Figure 1). The state presents high environmental heterogeneity due to complex topography and altitude range from 0 to 3,840 masl, with mountain ranges aligned to the coast, operating as an orographic barrier that maintains the humidity on the Pacific slope of the mountain ranges and creating dryer interior slope conditions, where dry and semi-dry climates prevail (Comisión Forestal del Estado de Michoacán, 2014; INEGI, 2014). At higher elevations, the environmental heterogeneity allows the presence of temperate forests composed of pine-oak, oak, fir, and cloud forests, in addition to tropical deciduous forests (Mas et al., 2017). All these landscape forest types constitute a complex land cover matrix where agricultural expansion reduces forest availability (Mas et al., 2017).OPEN IN VIEWER

Presence records were obtained from fieldwork in Michoacán from 2007-2019. Additional records were retrieved for Mexico from national and international databases: CONABIO, Global Biodiversity Information Facility, GBIF (https://www.gbif.org/) and SpeciesLink (https://splink.cria.org.br). Additionally, records documented in scientific articles were incorporated (Almazán-Catalán et al., 2013; Aranda & Valenzuela-Galván, 2015; Araya-Gamboa & Salom-Pérez, 2015; Botello et al., 2009; Cinta-Magallón et al., 2012). Once collected, the records were analyzed to exclude those with errors outside the historical distribution range or in current urban areas and filtered at 0.0083° (≈1 km²) to exclude redundancy.

The MaxEnt 3.3.3 software was selected to generate a prediction of the margay potential distribution as it presents the following characteristics: i) only uses presence data, ii) allows to show their habitats suitability spatially explicit, iii) quantifies the importance of each variable using a Jackknife test, and iv) allows the use of variables with continuous and categorical values (Elith et al., 2006; Qiao et al., 2015; Wei et al., 2018). The parameters used in this software normalize the contribution of each variable by using the maximum entropy approach (Elith et al., 2006; Feng et al., 2019; Phillips et al., 2006, 2017). The 19 WorldClim bioclimatic variables in ASCII format for model construction with a spatial resolution of 0.0083° (≈1 km2) were utilized (https://www.worldclim.org/bioclim). This set of variables represents annual trends, seasonality, and maximum and minimum temperature and precipitation values (Hijmans et al., 2005). These variables have been included in species distribution models, as these reflect the species tolerance to temperature and precipitation and adequately describe species distribution over extended areas (Dueñas-Lopez et al., 2015; Espinosa et al., 2017; Martínez-Calderas et al., 2016).

Parameters for software calibration in MaxEnt were 500 iterations and a convergence limit of 0.00001 during modeling, with logistic output format to create a response curve and calculate the contribution of the variables using the Jackknife test to identify the environmental variables with a significant contribution to the model creation (Moratelli et al., 2011). We used 70% of the presence records to generate the model and 30% to validate it. To evaluate the model performance, the Area Under the Curve (AUC) of the “Partial Receiver Operating Characteristic” (Partial-ROC) was used through the NICHE TOOLBOX platform (Osorio-Olvera et al., 2016; Peterson et al., 2008). The AUC measures the capacity of the model to discriminate between sites suitable for the margay presence and those where it is unsuitable (Elith et al., 2006). The partial ROC values oscillate from 0 to 2. Values greater than 1 indicate better performance than random (Osorio-Olvera et al., 2016). Therefore, we evaluated only the proportional areas on which the models made predictions, considering errors of omission of <10% and a differential weighting of errors of omission (false negatives) and commission (false positives) (Peterson et al., 2008). We generated a binary map (presence/absence) using the 10^th percentile presence threshold, widely used in SDM (Escalante et al., 2013; Liu et al., 2013).

To determine forest stands suitability, we employed the IDRISI.17.0 software, using Mas et al.'s vegetation and land use model of high definition (2017). We classify forests in primary condition as the more suitable in the following order: cloud and oak forest, mixed pine-oak forest, tropical dry deciduous forest, and tropical semideciduous forest. All constitute vegetation types that the scientific literature as margay habitat and where recent records were obtained in the field for Michoacán (Charre-Medellín et al., 2015). The size of forest stands is considered of adequate size if it includes a minimum area of 11 km², which is the average home range size (Carvajal-Villareal et al., 2012; Crawshaw, 1995; Kasper et al., 2016; Konecy, 1989; Oliveira et al., 2010). Subsequently, the selected stands were outlined according to the potential distribution model.

With the “Integral Index of Connectivity (IIC),” we measured the connectivity among forest patches. The IIC is a binary index that indicates whether a forest fragment maintains a connection to another, presenting values of 0 or 1. The 0 value indicates that forest fragments are not connected, whereas one suggests that they are connected (Decout et al., 2012; Pascual-Hortal & Saura, 2006). To measure this index, we used a distance threshold of 23 km; longer than the most extensive home range recorded: 21 km². The index was obtained with the software CONEFOR 2.6 (Saura & Torné, 2009) to select the most critical patches for connectivity.

The terrain resistance gradient for margay mobility in Michoacán was integrated considering the following species requirements. Home range (Carvajal-Villareal et al., 2012; Crawshaw, 1995; Kasper et al., 2016; Konecy, 1989; Oliveira et al., 2010), life history (Cuéllar et al., 2006; Oliveira, 1998; Pérez-Irineo et al., 2017; Pérez-Irineo & Santos-Moreno, 2016; Sunquist & Sunquist, 2002), and mobility features described in the literature. The variables that generated the resistance gradient included land use, human density/km², and vegetation vigor, for which a “Normalized Difference Vegetation Index” (NDVI) was used. The NDVI model corresponds to the “Eros” satellite image with the Moderate Resolution Imaging Spectroradiometer (eMODIS), dated March 11 to 20, 2019, corresponding to the dry season in Michoacán. These variables were weighted, giving 50% of the final resistance value to land use, as the vegetation is described as a fundamental variable to felids movement (Wikramanayake et al., 2004) and were proven successful in multiple studies of Felidae connectivity (Ashrafzadeh et al., 2020; Castilho et al., 2015; Correa-Ayram et al., 2014; Dickson et al., 2013; Dueñas-Lopez et al., 2015; Rabinowitz & Zeller, 2010) and assigned with the significant percentage of contribution to landscape resistance (Greenspan et al., 2021). Human density is considered to affect felid movement (Correa-Ayram et al., 2014; Dickson et al., 2013; Greenspan et al., 2021; Rabinowitz & Zeller, 2010) and is described as a significant limiting factor for mobility in felids as jaguars (Dueñas-Lopez et al., 2015). Therefore, we assigned a 30% resistance value to human density/ km². In connectivity studies with felids, quality variables are associated with vegetation prevalence like cover (Greenspan et al., 2021) or forest cover percentage (Rabinowitz & Zeller, 2010); we assigned 20% to NDVI.

The high-definition land cover model generated for the Michoacán state by Mas et al. (2017) was employed to assess the forest types and vegetation. Human density was obtained from the index of settlement and margination level database for 2020 from CONAPO (Appendix 1). The least-cost paths (LCP) criteria were implemented by using the Linkage Mapper 2 .0.0, an ArcGIS extension. These paths were outlined for all the important connectivity areas. The connectivity probability was obtained with the “circuitscape extension” for ArcGIS, and the selected calculation method was pairwise. The protected natural areas location polygons were obtained from Benzaury-Creel et al. (2009), and CONANP (2018).

Results

For the margay potential distribution model, 202 unique locations (coordinates) from México, including 21 locations obtained recently for Michoacán state during the 2007-2020 period,were used. The species potential distribution model exhibited adequate performance with an AUC of 0.96 and a “partial-roc” of 1.85 p<0.00. The crucial variables in the model generation were precipitation in the driest month (60.2%), isothermality (7.4%), and precipitation seasonality (6.6%) (Appendix 2).

The model generated for the state indicates suitable climatic conditions for the potential presence of margay in 82% of the state’s area if we consider all vegetation types and land uses. However, primary tropical or temperate forests are considered suitable habitats, constituting 50% of the state area or 13,061 km² of tropical and 15,576 temperate forests (Table 1). Fragmentation of native forests is broad across the state, but 116 forest stands were identified as containing suitable conditions for the margay, including seven (6% of suitable distribution) corresponding to protected areas. Three extend partially over the Biosphere Preserve Zicuiran-Infiernillo (RBZI by its Spanish acronym). Of these three, the one with the largest protected surface has a size of 2,197 km², 935.7 km² including under protection; another has an area of 3,647.7 km², and 7.8 km² under protection, and the last one that overlapped with the RBZI has an area of 6075.8 km², with 13.2 protected. Only one patch overlapped with the Tancitaro Mountain National Park. It has an area of 17,357 km², with 110.4 km² protected. Another forest stand overlaps with the Lagunas Costeras y serranías Patrimonial Reserve. This patch has an area of 83.5 km², with 28.4 km² protected. Another patch, extending over the El Jorullo Volcano Patrimonial Reserve, has an area of 11.3 km², with 7.7 km² protected. The last one extends over the Urban Ecological Park of Uruapan; it has an area of 65.3 km², with0.5 km² protected. Of the 11,597.9 km² we propose as primary habitat for the margay, only 1,103 km² (9.5%) are under protection.

OPEN IN VIEWER

Table 1.

Land Use Area in Distribution Model of Margay in Central Western Forests in Michoacán, México.

Vegetation type	Surface in Michoacán km2	Presence km2	No presence km2	Presence (%)	No presence (%)
Agriculture	15,849	11,678.50	4,170.40	73.69	26.31
Human settlements	858	585.5	272.5	68.24	31.76
Induced grasslands	6,165	5,536.30	629	89.80	10.20
Temperate forests	15,851	15,576.30	274.9	98.26	1.73
Tropical forests	17,442	13,061.70	4,380.40	74.89	25.11
Bodies of water	832.7	563.3	269.4	67.65	32.35
Mangroves	15.4	13.8	1.6	89.61	10.39
Popal-tular	117.7	112	5.7	95.16	4.84
No vegetation	66.3	61.4	4.9	92.61	7.39

In terms of ecoregions, suitable conditions were identified in decreasing area extent as follows: the Jalisco dry forest, the Sierra Madre del Sur pine-oak forests, and the Trans-Mexican volcanic belt pine-oak forests (Figure 2). In terms of forest categories, the potential distribution included 98.2% of the remaining temperate forests and 74.8% of the tropical forests (Table 1).OPEN IN VIEWER

The patches averaged 99.9 km² and ranged from 11.3 to 2,414.7 km² (the largest continuous stand of natural vegetation on the estate). The number of fragments and their size differed among forest categories (Table 1). From the 116 suitable patches, 21 were selected as the most important for connectivity according to the IIC analysis (Figure 3). These averaged 434.4 km², differing in size and forest category. In nine stands, the predominant forest category was tropical dry deciduous forest (Mean: 776.9: range = 63.9 km² to 2,414.7 km²), six in pine-oak forest (Mean: 94.2 km²: range = 64.5 km 2 to 206.2 km²), 5 in pine forest (296.3 km 2 average area ranging from 104.7 km² to 602.1 km²), and one in medium deciduous forest (83.5 km²).OPEN IN VIEWER

The circuit theory model map indicates how the observed margay areas of activity correlate with suitable habitat and connectivity probability (Figure 4), and the “least-cost paths” model identified 115 potential connectivity routes, which correspond to 476.9 km, averaging 4 km in length, with the longest path measuring 28.2 km (Figure 5). Evaluating the most critical suitable patches, the circuit theory model indicates high connectivity between the Jalisco dry forests and Trans-Mexican volcanic belt pine-oak forests. The lowest connectivity values were observed in the Balsas dry forests. With the “the least cost paths” methodology, 20 paths were designed totaling 598.3 km, with an average of 29.9 km, and the most extended patch measured 989.1 km long. A high overlap existed in the regions between high “probability of connectivity” and “low-cost paths” (Figure 6).OPEN IN VIEWER

OPEN IN VIEWER

Figure 5.

Least-Cost Paths.

OPEN IN VIEWER

Figure 6.

Probability of Connectivity and Least-Cost Paths for Most Important Patches.

Discussion

The results represent the first assessment in central western México to define suitable primary habitat available and the adequate area for the margay potential use, even though studies implying carnivore connectivity were performed previously (Balbuena-Serrano et al., 2021). Only two studies have assessed feline connectivity in central-western Mexico, including one for bobcats (Correa-Ayram et al., 2014) and another for ocelots (Avila-Coria et al., 2019). Nevertheless, we consider that addressing the habitat and connectivity requirements of the margay is a high priority due to the threat level at which it is listed.

The potential distribution model predicts a high presence probability in the state of Michoacán, and the potential distribution with suitable climatic conditions is significant for the Trans-Mexican volcanic belt pine-oak forests, Jalisco dry forests, and Sierra Madre del Sur pine-oak forests ecoregions, which in addition present also a high proportion of temperate (oak, pine-oak) and tropical forests (semideciduous and deciduous forests) respectively. We recommend paying particular attention to these areas with a high probability of presence for margay because they include areas that could function as habitat for margay. However, these areas face major land use change threats (Mas et al., 2009). Therefore, the presence of the species in these areas could be compromised.

The model exhibits a low presence probability in the Balsas and Bajío dry forests. This result can be explained by the most critical variable in the construction of the model which is the precipitation in the driest month, indicating that in these regions the precipitation levels in the dry season are not enough for margay requirements (Oliveira, 1998). Most relevant variables in model construction differ from those found by other authors (Charre-Medellín, 2009, 2017; Espinosa et al., 2017; Martínez-Calderas et al., 2016). This may result from the difference in the number, dispersion, and localization of records used, the scale at which the model was made, and the variables used to generate the model (Phillips et al., 2006).

Conservation actions should be conducted based on reliable, detailed, and spatially explicit species requirements (Rabinowitz & Zeller, 2010). Our analysis is the first based on specific margin requirements. The margay is probably the most susceptible felid to fragmentation in the state and country along with the jaguar (Oliveira et al., 2015; Quigley et al., 2017; SEMARNAT, 2010), considered vulnerable to human perturbation (Oliveira et al., 2015). Given the current situation in Michoacán regarding the change in land use, where native forests are being lost rapidly while the area of crops is increasing, coupled with the fact that the margay has arboreal habits (Mas et al., 2017), fragmentation of its habitat could threaten its subsistence in the state shortly. Furthermore, the margay needs to hunt to persist as a mesopredator. Margay’s diet includes small mammals such as rodents, birds, reptiles, amphibians, and invertebrates (Bianchi et al., 2011; Cinta-Magallón et al., 2012). Preserving Margay’s habitat involves habitat and prey as well.

Habitat loss and fragmentation affect felids dispersion (Abouelezz et al., 2018; Borah et al., 2016; Elliot et al., 2014), maintaining functional connectivity between isolated habitat patches could maintain gene flow and contribute to metapopulation dynamics (Stockwell et al., 2003). Functional connectivity is a helpful tool in felid conservation since it allows the recognition of areas with functional connectivity that present a greater richness of felids than those without (Gil-Fernandez et al., 2017).

The areas the model indicates as connective must be considered a priority for conservation for decision-makers (McRae et al., 2008). We identified enough habitat patches to maintain multiple margay populations in the state and the potential routes that margays could be used to move between patches. This connectivity would allow the establishment of a connectivity network between the protected natural areas of the central-western region of the country, whose margay populations could be critical in the species adaptability due in this region, the records of felids such as jaguarundi (Herpailurus yagouaroundi), ocelot (Leopardus pardalis), and margay correspond to maximum temperature and driest conditions in which the presence of these felid populations has been recorded in México (Charre-Medellín et al., 2015). Also, Michoacán connects the populations that are the limit of the margay’s distribution to the northern with the central and southern populations in México.

Our study observed high connectivity values in Coahuayana, Aquila, Tumbiscatío, Los Reyes, Uruapan, and Coeneo municipalities. Nevertheless, these regions are considered the main deforestation hotspots in Michoacán due to the advance in avocado cultivation in the Trans-Mexican volcanic belt pine-oak forests ecoregion and the advance of induced grasslands in the Jalisco dry forests and Sierra Madre del Sur pine-oak forests ecoregions (Mas et al., 2009). For this reason, the urgency of applying actions for conservation stands out.

Implications for Conservation

Natural protected areas are essential for conservation, but they alone cannot ensure sufficient habitat protection or effective conservation management for all species. These results provide a diagnosis of the current habitat availability based on margay requirements, highlighting relevant areas for population conservation and vegetation patches that promote connectivity.