Effects of variation in forest fragment habitat on black howler monkey demography in the unprotected landscape around Palenque National Park, Mexico

Study site

Our study was conducted in 34 forest fragments located within a ~10 km radius of PNP, Mexico (17°28′N, 92°03′W; Fig. 1). PNP encompasses 1,771 ha, ~1,000 ha of which are primary and secondary tall evergreen and tropical rainforest. PNP contains a stable population of numerous black howler groups (Van Belle et al., 2012; Klass, Van Belle & Estrada, 2020). The study fragments (Table 1; Table SM1) ranged in size from 0.2 to 36.2 ha (mean = 11.1 ha ± 9.8 SD; median size = 7.5 ha), for a total area of 377 ha. The fragments ranged in elevation from 34 to 380 m and were located 0.3–9.7 km (mean = 5.3 km ± 2.7) from the core area of PNP.

Fragments around PNP were found in pastureland matrix and on steep slopes (Fig. 2). Many streams formed the basis of narrow riparian corridors of remnant primary forest (Fig. 2C). Fragments were mainly composed of secondary growth forest and remnant primary forest, and were subject to disturbance (e.g., cattle grazing) of varying intensity. Human activities in the matrix consisted mainly of cattle pastures and agriculture (Figs. 1 and 2), including oil palm plantations; some fragments were adjacent to urban areas. Interviews of local people indicated that anthropogenic deforestation and fragmentation began about a century ago (A. Estrada, 2015, unpublished data). Between 1984 and 2001 ~33% of forest cover was lost in and around PNP and the number of fragments increased (Estrada et al., 2002).

Data collection

Fragment (explanatory) variables

We analyzed the effect of fragment-scale variables (hereafter referred to as fragment or explanatory variables; Tables 1 and 2) of different types, including dimensions, forest physical structure, and forest composition, and variables quantifying landscape connectivity in the vicinity (10–1,000s of m) of each fragment. We collected data on fragment area, shape, isolation from nearest fragment and from PNP, and number of treelines extending from the fragment for all 34 study fragments, by examination of high-resolution (1 m²) satellite images of the study area in Google Earth© (2018), combined with on-the-ground validation using handheld GPS data during surveys of each fragment and its perimeter (Tables 1 and 2).

We calculated an area-based isolation measure (Bender, Tischendorf & Fahrig, 2003), the proportion of forest habitat present within a pre-defined buffer zone around each fragment, excluding the area of the fragment itself. The size of the buffer was equivalent to a circle with a radius of the study species’ median dispersal distance. When dispersal distance has not been directly measured, it can be estimated based on the species’ home range size (Bowman et al., 2002). We used the mean home range size (9.3 ha) of five groups studied long term in PNP (Table 2), in the following formula (Bowman et al., 2002):

(a) dispersal distance = 7 × (diameter of home range size)

which gave an estimated median dispersal distance, and buffer radius, of 2,400 m.

In ArcGIS (ESRI, 2018), we calculated the proportion of forest cover in 2017 using the University of Maryland’s open access Global Forest Change Dataset (Hansen et al., 2013). We defined forest as pixels (30 m²) with ≥70% canopy cover. To determine this cut-off point, we first calculated the range and distribution of canopy cover values of random points in the study fragments and PNP at two time-steps, 2000 and 2017, to identify the canopy cover values of known howler monkey habitat. The first quartile value in the distribution of % canopy cover represents 25% of the data; in other words, 75% of the pixels sampled from the study fragments and PNP had canopy cover values greater than the first quartile value. In 2000, the first quartile was 73% canopy cover, and in 2017 it was 67.75%. We used the average of these values, 70% canopy cover, as the minimum to define pixels as “forest”; all other pixels were defined as “matrix”. We then calculated the proportion of forest in the buffer area around each fragment.

We collected data on stem density, tree genus (from which we calculated the Shannon’s diversity and Simpson’s evenness indices; Tables 1 and 2), tree height and diameter at breast height (DBH), and the number and DBH of stumps in 32 of 34 fragments using a modified Gentry Protocol (Gentry, 1982; Dias et al., 2014). Within each fragment, we recorded tree genus identity and DBH for all trees with DBH ≥ 10 cm along 10 transects, 50 m long and 2 m wide (total transect area per fragment = 1,000 m²). We recorded tree height for all trees with DBH ≥ 10 cm using a laser rangefinder (Larjavaara & Muller-Landau, 2013), and noted the DBH of all stumps with DBH ≥ 10 cm. Most trees were identified to the level of genus; the proportion of sampled trees per fragment which we were unable to identify ranged from 0 to 0.44 (mean = 0.14 ± 0.12). A total of 98 genera from 38 families were identified in the study fragments (mean = 16 ± 8.04 genera/fragment, range = 1–30; see Table SM2 for full list of genera).

Statistical analyses

Forest cover change and landscape connectivity

In ArcGIS, we calculated forest cover (see “Fragment (explanatory) variables” for details on the calculation of “forest” vs. “matrix”) for an area of 926 km² (92,605 ha) around PNP for two time steps, 2000 and 2017, using the Hansen et al. (2013) dataset. We subtracted the area of forest cover found in 2017 from that found in 2000, to arrive at the amount of forest cover lost in that 17-year period, and the remaining proportion of the landscape with forest cover in 2017.

We ran the 2017 forest cover layer through FRAGSTATS (McGarigal, Cushman & Ene, 2012) in R using the ConnCompLabel() and PatchStat() functions in the SDMTools (ver. 1.1–221) package (VanDerWal et al., 2014), to determine whether some fragments clustered together with greater connectivity in the landscape. The program used an eight-neighbor rule and returned a total of >6,000 fragments in the study area, using the 30 m² pixel size and definition of forest = >70% canopy cover of the input layer.

We identified our study fragments among the >6,000 defined by FRAGSTATS and found that 11/34 were embedded in what the program identified as one very large patch of forest including PNP and extending to the east and west of the park (Fig. 3). Although higher-resolution Google Earth (2018©) images and on-the-ground validation identified these fragments as individual patches of forest, these results indicated greater connectivity among these fragments, compared to the remaining 23. We used one-tailed Mann–Whitney U tests to compare the fragment variables between the two sets of fragments (11 “connected” and 23 “unconnected”; Fig. 3).

Forest loss and connectivity

From 2000 to 2017, forest cover in the study area decreased from 46,395 ha (50% forested) to 35,582 ha to (38.4% forested). Put differently, 23.3% of the forest cover present in 2000 was lost by 2017 (Fig. 3). Our study area is clearly composed of a northern section in which most forest cover had already been cleared by 2000, and a southern section, extending to the east, west, and south of PNP, where larger swathes of forest cover remain. While this area still contains more forest than the northern section, the forest loss analysis shows that most of the deforestation between 2000 and 2017 occurred in this southern area, making this previously large and relatively intact stretch of forest smaller and increasingly fragmented (Fig. 3).

The 11 fragments FRAGSTATS identified as more connected are located within or adjacent to this southern area. When comparing them to the remaining 23 fragments (Fig. 3), we found that 5/15 fragment variables differed significantly between the two groups: maximum tree height was significantly larger in the connected group (Mann–Whitney Tests: U = 48, z = −2.50, p = 0.0062); Simpson’s evenness index was significantly higher in the connected group (U = 51, z = −2.38, p = 0.0087); Shannon’s index was significantly higher in the connected group (U = 28.5, z = −3.29, p = 0.001); distance to PNP was significantly shorter in the connected group (U = 50, z = 2.80, p = 0.0026); and the proportion of forest cover in the fragment buffer areas was significantly larger in the connected group (U = 15, z = −4.09, p < 0.00001).

Effects of fragment variables on black howler population size, density and composition

The Spearman correlation heatmap (Fig. 4; see Table SM3 for all correlation coefficients and p-values) showed that mean group size was significantly positively correlated with four other response variables: population size, population density, the mean #AM/grp, and the mean #AF/grp. Given this redundancy, we excluded mean group size from the RDA analysis.

The RDA model with 15 fragment variables accounted for over two-thirds of the variation in the demography of the black howler population in forest fragments (R² = 0.69). This model was globally significant (p = 0.003, F = 1.967, df = 15), and RDA axis 1 was also independently significant (p = 0.018, F = 20.93, df = 1), although axis 2 was not (p = 0.223, F = 10.63, df = 1; Fig. 5). Of the individual explanatory variables, the marginal effects of three were statistically significant (fragment area, Simpson’s evenness index, and the number of stumps) and two were close to the significance value (stem density and Shannon’s index; Table 3). The variables most strongly driving axes 1 and 2, which together accounted for 44% of the explained variation in the demographic variables (Fig. 5), were (in descending order) minimum tree height, Shannon’s index, minimum tree DBH, stem density, fragment area, and Simpson’s evenness index (axis 1); and area, the isolation distance to the nearest fragment, Simpson’s evenness index, the number of stumps, and the mean DBH of stumps (axis 2; Table 4).

Geographically close fragments did not cluster together in the triplot, indicating that proximity in the landscape did not predict similarity in black howler demography (Figs. 3 and 5). Similarly, those fragments located in the more connected southern area around PNP did not cluster together in the triplot; aside from the Shannon’s and Simpson’s indices, the main variables that drove the variation in fragment demography were not those that differed among the more vs. less connected fragments.

The proportion of AF and average #AF/grp variables clustered together (and with the fragment population size), correlating negatively with both axes (Fig. 5; Table 4). These response variables were most closely positively associated with fragment area, the number of treelines extending from the fragment, Shannon’s index, stem density, the proportion of forest cover in the buffer, and the shape index, although area had the strongest effect. These variables showed the strongest negative correlations with Simpson’s evenness, the minimum tree DBH and height, and the mean stump DBH (Fig. 5; Table 4).

Variables related to adult males also grouped together (and with population density), correlating negatively with axis 1 and positively with axis 2 (Fig. 5; Table 4). These variables were most closely positively associated with stem density, the distance from PNP, the Shannon index, maximum tree DBH, mean stump DBH, and the isolation distance to the nearest fragment (particularly density and the proportion of AM in the fragment; Fig. 5; Table 4); of these, Shannon’s index, stem density, and isolation from the nearest fragment had the strongest effect. The AM variables were negatively correlated with minimum tree height and DBH and the number of stumps. Population density was negatively correlated with area, the proportion of forest cover in the buffer, minimum tree height, and the number of stumps (Fig. 5; Table 4).

Several significant pairwise relationships emerged from the Spearman correlation heatmap (Fig. 4) that further clarify the patterns found in the RDA analysis. Some discrepancies in the heatmap and RDA results are expected, as RDA performs a joint analysis of the two datasets (i.e., demographic and fragment variables) to assesses the degree of co-variation between them; in other words, RDA identifies patterns that occur in both datasets, and ignores patterns that are unique to only one (Ulhaq et al., 2013).

Among the six response variables, only the mean #AF/grp was significantly and positively correlated with population size. Among the explanatory variables, the proportion of forest cover in the buffer was significantly positively correlated with maximum tree height and Simpson’s evenness index, and negatively correlated with the isolation distance from PNP. These relationships further supported the results of the analysis comparing the more connected fragments to the less connected ones, as defined by FRAGSTATS. Stem density was significantly negatively correlated with minimum tree height and DBH and positively correlated with Shannon’s index, and also (unlike the relationship seen in the RDA) positively correlated with the number of stumps. Similarly, the strong positive correlation between the Shannon’s and Simpson’s evenness indices seen in the heatmap was not apparent in the RDA.

The heatmap showed several significant relationships among individual explanatory and response variables. Fragment area correlated positively with population size and negatively with population density. Density also correlated negatively with the number of treelines extending from the fragment, and population size correlated negatively with the isolation distance from the nearest fragment. The proportion of AM and mean #AM/grp both correlated positively with stem density, while the mean #AM/grp also correlated negatively with minimum tree height.

Adult females

Based on additional analyses of this population (Klass, Van Belle & Estrada, 2020) and of black-and-gold howler monkeys (Oklander & Corach, 2013) indicating that howler monkey females may be more philopatric in fragmented landscapes, we hypothesized that in our study fragments, AF become more philopatric than in continuous forests because of increased competition for resources (Gros, Poethke & Hovestadt, 2009). We therefore expected variables relating to AF to positively correlate with variables indicative of habitat quality, where we defined habitat quality as increasing with increasing similarity to primary, undisturbed rainforest.

The main patterns that emerged for AF can be summarized thus: there were more AF/group and a higher proportion of AF in the population in larger fragments with better connectivity, characteristics indicative of abundant secondary growth (e.g., high stem density, smaller minimum tree size), and those with more diverse vegetation, but lower Simpson’s evenness indices (indicative of a tree composition dominated by only one or a few genera; Morris et al., 2014). Adult females were also more abundant in fragments with smaller mean stump size, which may be an indication of less human disturbance, as larger stumps were often the result of the harvesting of large, mature trees by local people.

These results are somewhat surprising, not least because the RDA shows that while there were more AF in larger fragments and in fragments with more diverse tree genera, there was no direct correlation between area and tree diversity, and a negative correlation between area and the evenness of the forest composition; we would expect larger fragments, with fewer edge effects, to have more diverse vegetation, with a more even Simpson’s index (Santos et al., 2008; Dias et al., 2014). However, fragment size does not always positively correlate with food availability (Dias et al., 2014), and it is possible that the ability of black howlers to maintain dietary breadth in fragments with different levels of tree species richness and diversity (Dias et al., 2014) minimizes the importance of tree composition for black howler females. Productive secondary growth in fragments may provide howlers with abundant food resources (Arroyo-Rodriguez & Dias, 2010), albeit different from those found in undisturbed primary rainforest.

Another indication that secondary growth food resources may be important for howlers in fragments is the positive, although minor, association we found between the shape index and many of the response variables. A high shape index indicates a complex fragment shape with more forest edge, which is associated with changes to vegetation relative to primary forest that could be considered indicative of low-quality habitat with fewer resources (Ewers & Didham, 2007; Arroyo-Rodriguez & Dias, 2010). However, fragments with more complex shapes and more edges may also have an abundance of new, secondary growth, providing howlers with an important source of young leaves to feed on (Cristobal-Azkarate et al., 2005).

We found that fragment isolation negatively affected the abundance of AF, which may stem from females avoiding costly, long-distance dispersal through the matrix. Additional studies have shown that habitat fragmentation may reduce female dispersal frequency (Oklander & Corach, 2013) or distance traveled (Minhos Rodrigues, 2012) in arboreal primates.

In continuous forest populations of ursine howler monkeys (A. arctoidea), the costs of dispersal in terms of reproductive success and mortality risks were much greater for females than for males, and females would thus prefer to remain, and have their daughters remain, in the natal group and territory when conditions allowed (e.g., <4 AF in the group; Pope, 2000). When obligated to disperse, females encountered resistance to immigration into existing groups from resident females, thereby increasing dispersal costs. Evidence from PNP indicates that female black howler dispersal patterns in continuous forests may be similar to ursine howlers (Van Belle et al., 2012; S. Van Belle, 2019, unpublished data). Thus, a possible explanation for the patterns found in this study is that black howler females in fragments are more likely to remain in the natal group, establish a new group in the natal fragment, or disperse short distances to nearby fragments, rather than disperse long distances across the matrix. This is because the benefits of remaining in the natal group or fragment, in terms of access to resources and reproductive success, are even greater than in continuous forests, as the potential costs of dispersal outside of the natal fragment, in terms of mortality risk and energetic costs, are likely even higher in a fragmented landscape compared to continuous forest (Bowler & Benton, 2005).

Avoidance of longer-distance dispersal by females may only be possible for those in larger, less-isolated fragments, that can support multiple groups or are near easily accessible additional fragments; thus over time populations in these fragments might show higher proportions of AF and larger populations overall, as we found in this study. This strategy may not be available to females born in smaller, more isolated fragments, as establishing new groups in these small areas would be difficult; additionally, our results showed that populations in smaller, isolated fragments have fewer AF, not more.

Similarities between the study fragment population and populations in PNP and other continuous forests in group size and mean number of AF (Klass, Van Belle & Estrada, 2020) may indicate that, when possible, females in fragments may be more commonly avoiding high dispersal costs by establishing new groups in the natal fragment or dispersing to nearby fragments, rather than remaining in the natal group more often than in continuous forest populations. The high per-fragment densities in our study fragments may provide some support for this scenario; these densities were often the result of multiple groups simultaneously occupying fragments smaller than some black howler group home range sizes observed in continuous forests (Gavazzi et al., 2008; De Guinea et al., 2019).

Adult males

Based on additional analyses of this population (Klass, Van Belle & Estrada, 2020) and on the hypothesis that males continue to disperse in the fragmented landscape to seek reproductive opportunities (Gros, Poethke & Hovestadt, 2009) and suffer higher mortality during dispersal (Banks et al., 2005), we predicted that the proportion of AM in the fragment and mean #AM/group in forest fragments would correlate positively with variables quantifying connectivity around fragments.

The main patterns that emerged for AM can be summarized thus: there were more AM/group and a higher proportion of AM in the population in fragments that were more isolated from nearby fragments (and to a lesser extent, from PNP), with characteristics indicative of abundant secondary tree growth (e.g., high stem density, smaller minimum tree size) but larger maximum tree DBH, with fewer, larger stumps, and more diverse vegetation, but lower Simpson’s evenness indices. Thus, contrary to our prediction, we found AM were more abundant in more isolated fragments. The abundance of AM was uncorrelated to fragment area, in contrast both to Dias et al. (2015) study on black howlers in fragments in Campeche, Mexico, in which groups were found to include more males in smaller fragments, and Puig-Lagunes et al. (2016), who found that numbers of all age-sex classes increased with fragment size in mantled howlers.

Pope (2000) found that ursine howler males in continuous forests most often dispersed very short distances, to neighboring groups, and when attempts to join or take over those groups were unsuccessful, they returned to the natal group, thus incurring generally low dispersal costs. Evidence from PNP indicates that black howler male dispersal patterns in continuous forests may be similar (S. Van Belle, 2019, unpublished data). However, AM in forest fragments are faced with a scarcity of nearby, familiar groups to take over or join, and as a result may be dispersing at similar rates but forced to travel longer distances during dispersal in the fragmented landscape as compared to continuous forest, to seek reproductive opportunities in more isolated fragments.

While engaged in longer-distance dispersal through the matrix, black howler males may be subject to increased mortality risks. Pozo-Montuy, Serio-Silva & Bonilla-Sanchez (2011) showed that while overall rates of emigration may be similar in continuous and fragmented habitat for black howlers in Balancán, Mexico, increased mortality for dispersers moving through the matrix in fragmented landscapes may result in reduced effective dispersal rates and gene flow. Hunting of black howlers occurs in our study area, and dogs are known to kill monkeys forced to move on the ground (Pozo-Montuy, Serio-Silva & Bonilla-Sanchez, 2011), which occurs in the matrix in our study site (Fig. 2D). The occurrence of increased male mortality in the fragmented landscape is supported by additional analyses of this population, which indicated that there were fewer AM in forest fragments relative to the population in PNP (Klass, Van Belle & Estrada, 2020).

The minimal effect of distance from PNP on the demography of fragment populations indicated that the mainland-island metapopulation model (Hanski, 1998) does not apply to the study population and the nearby large, stable population in PNP; that is, PNP does not act as a significant source population for the fragments (source-sink dynamics; Pulliam, 1988), including those that are closer to the park. However, it is possible that the larger, more productive fragments in the landscape, which are able to support multiple groups, may act as a source for small, sink populations in the more isolated, lower-quality fragments (Pulliam, 1988; Arroyo-Rodriguez, Mandujano & Benitez-Malvido, 2008). According to Banks et al. (2005), while high male dispersal mortality may bias the overall population sex ratio towards females, small populations maintained by immigration (i.e., sink populations) may have a relatively higher proportion of the dispersing sex. This pattern is in line with our observations of the study population: both fewer AM overall compared to PNP (Klass, Van Belle & Estrada, 2020) and higher proportions of AM and more AM/group in more isolated fragments.