Population densities of Hylobates agilis in forests with different disturbance histories in Ulu Muda Forest Reserve, Malaysia.

Small ape habitat throughout Malaysia is rapidly being lost, degraded, and fragmented, and the effects of these changes on the abundance on this taxon are currently unknown. This study assessed the group density of Hylobates agilis in virgin forest, previously logged forest (1960s-1990s), and recently logged forest (2015-2017) of the Ulu Muda Forest Reserve (UMFR), Kedah, Malaysia. We conducted fixed-point active acoustic triangulation at nine survey areas to estimate group density. We used vegetation "speed plots" and satellite imagery to quantify habitat characteristics and used model selection to identify ecological predictors of group density variation. The estimated group density of H. agilis in UMFR was 4.03 ± 0.14 groups km-2 , with an estimated total of 2927 ± 102 groups in areas below 450 m a.s.l. in UMFR. Group density did not differ significantly among habitat types. The best ecological predictors for group density were canopy cover and proportion of deforested area. Areas with recent deforestation were associated with relatively high group densities, suggesting compression of the populations persisting in these habitat types. The consistently high group densities detected in all forest types emphasizes the importance of degraded forest as habitat for H. agilis. Because of the threats to small apes in Malaysia, and the uncertain status of most populations, we recommend a nationwide population census and regular monitoring to inform conservation planning and implementation. Most urgently, we call for immediate and permanent protection of UMFR and other forests in the Greater Ulu Muda landscape to protect the globally significant populations of H. agilis, as well as other charismatic and threatened megafauna, birds, and flora in the area.

corrected Akaike Information Criterion

National Institute of Advanced Industrial Science and Technology

analysis of variance

above sea level

American Society of Primatologists

Advanced Spaceborne Thermal Emission and Reflection Radiometer

confidence interval

coefficient of variation

Diameter at Breast Height

Digital Elevation Model

Global Positioning System

International Union for Conservation of Nature

Linear Model

Listening Post

Ministry of Economy, Trade and Industry, Japan

National Aeronautics and Space Administration

Passive Acoustic Monitoring

Department of Wildlife and National Parks Peninsular Malaysia (Jabatan Perlindungan Hidupan Liar dan Taman Negara Semenanjung Malaysia)

Permanent Reserved Forest

survey area

standard deviation

University of Maryland

Ulu Muda Forest Reserve

INTRODUCTION

Small apes (Hylobatidae) comprise 20 extant species, all of which are threatened and most of which have declining populations (Fleagle, 2013; International Union for Conservation of Nature, 2020; Rawson et al., 2011). Despite being broadly distributed across South and Southeast Asia and southern China, the habitat requirements of this taxon are fairly specific (Fleagle, 2013; Lappan & Whittaker, 2009). Small apes are vulnerable to deforestation, habitat degradation, and fragmentation as they are exclusively arboreal and live in the middle and upper canopy layers of the rainforest (Brockelman & Ali, 1987; Gittins, 1983; Whittaker, 2009), and they have not been reported to persist in monoculture plantations (Choudhury, 2006; Zhang et al., 2010). As true brachiators, small apes rely on continuous and dense canopy cover for locomotion (Cannon & Leighton, 1994), and their movements are limited by even small gaps in the forest canopy (Cheyne et al., 2013). For example, where small roads bisect the forest, Hylobates lar and Hylobates pileatus reduce the use of areas near the roadside and are less likely to cross roads than predicted by a null movement model (Asensio et al., 2021). As a result, the home ranges of roadside groups are partly delineated by roads, even though the roads are narrow and relatively permeable to gibbon movement (Asensio et al., 2021).

Hylobates agilis has a spatially disjunct distribution, with populations on Sumatra and on mainland Asia. The distribution range of H. agilis on mainland Asia is restricted to a narrow band stretching from southernmost Thailand (Yala Province) to the area between the Muda River and Perak River in three northern states (Kelantan, Kedah, and Perak) in Peninsular Malaysia (Gittins, 1978; Marshall, 1981; Mootnick, 2006). H. agilis lives in groups averaging ca. four individuals (Lee et al., 2015), usually with a pair‐bonded adult male and female, and up to four offspring. They defend territories of ca. 30 ha (Gittins, 1980) and produce loud, high‐pitched vocal duets, usually in the morning hours. H. agilis is listed as Endangered on the IUCN Red List, with a declining population, primarily due to habitat loss from anthropogenic disturbance, with poaching as a secondary threat in some areas (Geissmann et al., 2020). The situation for the population on mainland Asia is particularly dire, with much of its former range having been deforested, degraded, or fragmented in the last few decades, resulting in local extinctions of many populations (Lappan, in prep).

Given rapid anthropogenic loss, modification, and degradation of small ape habitat throughout Malaysia (Lappan & Ruppert, 2019), detailed knowledge about small ape abundance and occurrence, and their ability to persist in disturbed and undisturbed environments is critical for the development of an effective national conservation plan. In Malaysia, most studies on small apes took place in the 1970s and 1980s (e.g., Chivers, 1974, 1977, 1980; Johns, 1985, 1986). Though these provide useful baseline data on group density, group composition, and behavioral ecology, little information about the distribution and abundance of Malaysian taxa has been published in the past three decades.

The aim of this study was to assess the population status of H. agilis in the Ulu Muda Forest Reserve (UMFR) in Kedah, Malaysia, part of the largest remaining forested area within the H. agilis distribution range on mainland Asia. We conducted fixed‐point active acoustic triangulation at nine areas in UMFR to estimate group densities and the total number of H. agilis groups in lowland forest (<450 m a.s.l.) in UMFR. Because of the reported sensitivity of small apes to habitat disturbance, we predicted that group density would be higher in virgin forest than in previously logged and recently selectively logged forest habitats in the area. We also aimed to determine what habitat characteristics best predict H. agilis group densities at UMFR.

METHODS

Study area

The greater Ulu Muda landscape consists of several forest reserves (Table 1 and Figure 1) of which the UMFR, where we conducted the study, is the largest. The Ulu Muda Forest Reserve (UMFR, 100°51′30″−101°7′9″ E, 5°43′57″ N−6°16′23″ N) covers an area of ca. 1050 km2, consisting mainly of lowland dipterocarp forest with elevation ranging from 97 to 1256 m a.s.l. (Woo & Perumal, 2019), with an area of 725.6 km2 below 450 m a.s.l., 358.8 km2 between 450 and 850 m a.s.l., and 87.9 km2 above 850 m a.s.l. (Figure 3). The recorded daily minimum and maximum temperatures at the study area range between 15.2–25.8°C and 23.5–38.6°C, respectively. UMFR receives 2160–3000 mm/year of rainfall (Tan et al., 2019). Despite its recognition as a critical water catchment area for the states of Penang, Kedah, and Perlis in Malaysia and its legal status as “Permanent Reserved Forest” (PRF), most of the UMFR is classified as Production Forest within the jurisdiction of the Forestry Department, meaning that concessions can be issued to contractors to carry out selective logging, in which trees of designated species above certain sizes (e.g., 55 cm diameter at breast height [DBH] for non‐dipterocarps, 65 cm DBH for dipterocarps) are felled in a controlled manner and carried along a skid track constructed using a bulldozer (Saiful & Latiff, 2014). Land with this status has also been cleared to make way for timber tree plantations or subjected to excision, where the PRF status is removed and the area becomes open to conversion for other land use (Law, 2020).

Table 1

Forest reserves in the greater Ulu Muda landscape (Suksuwan, 2008).

Name of forest reserve	Area (km2)
Ulu Muda (UMFR)	1050.6
Pedu	153.0
Padang Terap	127.9
Bukit Keramat	102.3
Chabar Besar	88.3
Bukit Saiong	81.9
Chabar Kecil	11.8
Ulu Muda (extension)	13.6
Total	1629.3

Figure 1

Forested area in the Greater Ulu Muda landscape (main map) and position of Ulu Muda within Peninsular Malaysia (small map bottom right).

Figure 3

Elevation zones within UMFR. The area for each elevation category is: yellow: <450 m = 725.6 km2; light green: 450–850 m = 358.8 km2; white: >850 m = 87.88 km2).

With the help of trained field assistants (Supporting Information: Table S1), we conducted active acoustic surveys using fixed‐point triangulation following Brockelman and Ali (1987) and Brockelman and Srikosamatara (1993) in nine survey areas (SA; Table 2) from April 2018 to February 2019, covering a total SA of 33.4 km2. Trained field assistants with experience as principal investigators on other small ape population surveys helped to train new field assistants for the first few field trips, before the latter became experienced enough to lead others. We placed three SAs in each of three forest types (virgin, previously logged, and recently logged). All SAs were placed in forests below 450 m a.s.l. due to limitations in terms of logistics and budget. Forests were classified into habitat types using Landsat images (Hansen et al., 2013) and historical information from local forest guides. The habitat types were: (1) virgin forest, approximately 5 km northeast from the Earth Lodge, a small ecotourism lodge for tourists and researchers (6°06′54″ °N, 100°57′50″ E), (2) previously logged forest, located around the upstream area of Lake Muda, which was selectively logged between the 1960s and 1990s (Suksuwan, 2008), and (3) recently logged forest, located around the Ulu Legong Hot Spring (5°49′00″ N, 100°56′08″ E), where the forest is adjacent to villages and plantations, and was selectively logged between 2015 and 2017 (Table 2 and Figure 2). We selected the SA based on accessibility (e.g., presence of trails, topography), with the criterion that SAs were placed a minimum distance of 2 km apart to prevent pseudo‐replication (i.e., sampling the same H. agilis group(s) twice).

Table 2

Survey areas with forest type, location, and elevation range.

Survey area	Forest type	GPS coordinates	Elevation range (m)	Mean ± 95% CI elevation (m)
Sg. Labua	Previously logged	6°6′52″ N, 100°58′9″ E	222–301	262.3 ± 19.4
Tualang Trail		6°5′48″ N, 100°57′31″ E	190–276	249.4 ± 23.8
Lubok Petai		6°7′2″ N, 100°56′44″ E	194–279	248.6 ± 20.2
Camp Ninja	Virgin	6°8′24″ N, 100°59′22″ E	151–205	181.5 ± 13.7
Sg. Batu Hitam		6°7′22″ N, 100°59′16″ E	200–325	257.1 ± 28.3
Sira Siput		6°7′58″ N, 100°58′18″ E	165–321	234.7 ± 41.7
Air Legong	Recently logged	5°49′38″ N, 100°56′15″ E	290–362	340 ± 14.8
Bukit Gadung		5°47′42″ N, 100°56′5″ E	355–427	398.8 ± 21.0
Lata Gading		5°48′56″ N, 100°57′4″ E	250–334	302.7 ± 21.9

Abbreviation: CI, confidence interval; GPS, global positioning system.

Figure 2

Locations of the survey areas (colored circles) in the Ulu Muda Forest Reserve within the different forest types. The green polygon indicates the forest reserve boundaries (Kedah Forestry Department).

GROUP DENSITY

We estimated the group density (D) of H. agilis in UMFR following Brockelman and Ali (1987) as: where n is the total number of groups heard at each SA.  p(m) is the calling probability, which is the estimated proportion of groups expected to be detected during a survey period of m days, to account for the fact that all groups in a SA may not vocalize during the acoustic survey period (Thinh & Rawson, 2011). We estimated from the data recorded during the surveys using the spreadsheet package by Thinh and Rawson (2011).  is the area of the SA (in km2).

We used fixed‐point active acoustic triangulation following Brockelman and Ali (1987) and Brockelman and Srikosamatara (1993) to estimate the location of duetting H. agilis pairs and to estimate the number of groups detected in each SA. In each of the nine SAs, we established three listening posts (LP) spaced 300–400 m apart, resulting in a total of 27 LPs. All LPs were placed on prominent locations such as small hills or along ridges of 300–400 m a.s.l. to improve our ability to detect duets coming from all directions (Brockelman & Ali, 1987). We conducted surveys in each SA for three consecutive days from 7 a.m. to 11 a.m. (UTC + 08:00). These times were selected because Gittins and Raemakers (1980, p. 74) showed that the great calls of H. agilis in Sungai Dal, Malaysia peaked at around 7 a.m. (which is the earliest time of sunrise in UMFR), and our preliminary surveys confirmed that few great calls occurred in UMFR before 7 a.m. or after 11 a.m. During surveys of a SA, pairs of listeners at each LP recorded (1) the estimated compass bearing for all great calls/duets of a mated pair, (2) the start and end times for the calling bout, and (3) estimated distance (in m) between the singing pair and the LP. Only the great calls, which are sung by the adult female of the group and characterized by the long, high‐pitched notes that ascend before descending, followed by a coda (i.e., notes sung by the male at the end of the female's great call), were recorded, because they indicate the presence of a pair defending a territory (Haimoff & Gittins, 1985).

To determine the number of groups per SA, the locations of the groups whose duets were detected from two or more LPs were estimated using triangulation. Home ranges for Hylobates spp. average about 30 ha in size (Bartlett, 2011), so we followed Brockelman and Ali (1987) and Brockelman and Srikosamatara (1993) in assuming that duets that were produced at locations >500 m apart were considered as distinct groups, as 500 m is the approximate diameter of the taxon's home range. Other relatively recent studies of H. agilis have also used this assumption (Lee et al., 2015; O'Brien et al., 2004). Groups heard singing from locations that could not be triangulated (because they were not heard from two LPs simultaneously) were excluded if they were heard only once over the 3‐day survey period. In case a group was heard duetting from two LPs at the same estimated location but at different times or was heard calling more than once at the same location, it was also included in our density estimates if the estimated distance to the singing group from the LP was <1000 m (Thinh & Rawson, 2011). This modification was applied to the standard protocol because of the hilly topography as songs detected at one LP were in some cases not detectable from other LPs because land features obstructed the sound.

We estimated the maximum detection distance as 1000 m, as previous studies have shown that calls of Hylobates spp. can travel that far under good listening conditions (Brockelman & Ali, 1987), and groups that were further away from the LPs were unlikely to be detected (Vu et al., 2018) because their songs would generally be obscured by nearer songs or background noises. From our surveys, almost all songs that we could triangulate were estimated to be within a 1000 m (radius) from any of the LPs. We excluded songs that were triangulated to locations that appeared to be outside of the maximum detection distance (i.e., >1000 m from all LPs).

To estimate the area sampled for each SA more accurately, we used 2018 forest cover data from Hansen et al. (2013) and a digital elevation model (National Aeronautics and Space Administration/Ministry of Economy Trade and Industry/Advanced Industrial Science and Technology/Japan Spacesystems & U.S./Japan ASTER Science Team, 2018) to exclude the following areas from the effective SA: (1) areas that were acoustically obstructed from all three LPs – the sound of duets can be heard from up or down a slope (i.e., from a valley adjacent to a LP or from a higher point on the same slope), but not from beyond an intervening hill or ridge; therefore, groups that sing in a valley directly behind a hill or ridge from the perspective of a LP are unlikely to be heard from that LP (Brockelman & Ali, 1987), (2) water bodies (lakes and rivers); (3) riparian areas (areas < 50 m from rivers and lakes), because in UMFR these areas often have low canopy cover due to intermittent flooding; and (4) human settlements or areas that have been converted for agriculture (Thinh & Rawson, 2011).

To estimate the probability of acoustic detection (or calling probability ) during the 3‐day survey, we used the methods of Thinh and Rawson (2011), but with the modification that only groups that were detected from a location <600 m from the nearest LP were included in the calling probability estimation (Figure 4). Songs of groups with home ranges near an LP were likely to be detected whenever they were produced during the sampling period, as all of their songs should be produced at distances of <1000 m from the nearest LP (i.e., locations within the SA). However, groups with home ranges that are only partially within the SA may produce songs that are not detected from the LPs because they are produced >1000 m from the nearest LP. Including these groups in estimates of calling probability may therefore result in underestimation of the calling probability and overestimation of group density. To avoid this problem, Vu et al. (2018) suggested including only songs recorded from within 500–700 m of any LP, as most or all of the home range of the group singing from these locations should be within the SA.

Figure 4

Acoustic survey at Lata Gading and triangulation of the Hylobates agilis groups. Only groups (numbered green circles 2 and 3) detected at or within 600 m from any listening post were included for calling probability calculation.

The weather during the surveys was consistently dry, with sporadic drizzles, mostly in the afternoon. This sporadic rainfall was unlikely to meaningfully affect the calling probability as groups were able to sing before and after the rain.

TOTAL NUMBER OF GROUPS

To estimate the total number of groups (G) in lowland forest areas (i.e., <450 m a.s.l.) in UMFR and to place 95% confidence intervals around that estimate, we used the software package designed by Thinh and Rawson (2011) for small ape population surveys. We used the same correction factor p(m) for all areas, assuming that calling probability does not vary across SA. To estimate the total area of habitat below 450 m a.s.l. in UMFR, we used forest cover data for 2018 from Hansen et al. (2013) and a digital elevation model (National Aeronautics and Space Administration/Ministry of Economy Trade and Industry/Advanced Industrial Science and Technology/Japan Spacesystems & U.S./Japan ASTER Science Team, 2018) to exclude areas that were clearly not habitat for this taxon, such as deforested areas (<25% forest cover at 5 m height; Hansen et al., 2013), water bodies, riparian areas (areas < 50 m from the reservoir and rivers, which tended to have low forest cover due to episodic inundation), and islands within the reservoir, to avoid overestimating the area of potential habitat (Thinh & Rawson, 2011). We also excluded areas >450 m a.s.l. as we did not survey SA with elevation >450 m a.s.l.

HABITAT CHARACTERISTICS

We compared group densities across the three forest types using one‐way analysis of variance (ANOVA). To compare habitat quality among forest types and identify the best ecological predictor(s) of H. agilis group density, we quantified indicators of habitat quality in two ways: (1) vegetation “speed plotting” and (2) satellite imagery. Following Hamard et al. (2010), we established ten 10 m × 10 m vegetation plots (VP; Supporting Information: Figure S1) around the three LPs in each SA, resulting in 90 VPs in total. We placed three VPs approximately 50 m from each LP, with one each to the north (0°), southeast (125°), and southwest (225°) of the LP (Supporting Information: Figure S1). We placed the final (tenth) plot for each SA in the approximate midpoint of the three LPs (Supporting Information: Figure S1). We recorded the following vegetation characteristics in each VP: (1) mean canopy cover percentage (scored visually using a GRS Densitometer™; estimated from three points within the VP; rounded to the nearest 5%), (2) total number of trees ≥10 cm DBH, (3) DBH of all trees ≥10 cm DBH, and (4) height of all trees ≥10 cm DBH. The DBH of each tree was consistently measured at a standard height of 130 cm above ground. To yield tree density (number of trees km−2), we divided the number of trees ≥10 cm DBH within a VP by 0.001. We conducted statistical tests using IBM SPSS Statistics Campus Edition V24.0, using α = 0.05. The VP variables were standardized into z‐scores before analysis. The elevation of each VP was recorded using DEM data from National Aeronautics and Space Administration/Ministry of Economy Trade and Industry/Advanced Industrial Science and Technology/Japan Spacesystems & U.S./Japan ASTER Science Team (2018).

To determine whether SA across the three forest types were similar in habitat characteristics, we used the Kruskal–Wallis test (Kruskal & Wallis, 1952) to compare vegetation variables (i.e., mean canopy cover, mean DBH, mean tree height, number of trees km−2, and elevation) from VPs across SA in different forest types.

We also measured percentage of forest cover in the year 2000 and proportion of forest cover lost between 2000 and 2018 for each SA using satellite imagery from the Hansen et al. (2013). The year 2000 cloud‐free composite LandSat image assigns a value of 0 to 100 to each pixel, which corresponds with the estimated percentage of forest cover for that pixel (Hansen et al., 2013). We calculated the mean percentage of forest cover in 2000 for each SA as: where α is the value for forest cover of each pixel, b is the number of pixels within the SA assigned each value of α.

To calculate the proportion of land deforested between the year 2000 and the study period (2018/19) for each SA, we measured the area (km2) deforested between 2000 and 2018 from the LandSat images (Hansen et al., 2013) and then divided this area by the total forested area (km2) in 2000 for the same SA. An area was considered deforested in the LandSat image when canopy cover was <25% at 5 m tree height (Hansen et al., 2013). As the variables mean forest cover percentage in the year 2000 and proportion of area deforested between 2000 and 2018 are proportional data, we transformed these variables to approximate a normal distribution using the arcsine‐square‐root transformation before statistical analysis.

To better understand the relationship between habitat characteristics and group densities, we built Linear Models (LM) with H. agilis group density as the outcome variable and different combinations of habitat variables as predictor variables. We then conducted model selection using the corrected Akaike Information Criterion (AICc; Akaike, 1974), an information criterion that includes a correction for small sample size (Hurvich & Tsai, 1989), to determine which model variable(s) best explain variation in group density, while prioritizing model simplicity (Wagenmakers & Farrell, 2004). Separate model selection procedures were conducted for (1) VP data: mean DBH, mean tree height, tree density and mean canopy cover and (2) remotely sensed data: mean percentage of forest cover in year 2000, forest lost between 2000 and 2018, and mean elevation of each SA (Hansen et al., 2013). Models were built for each possible combination of one, two, or three variables, and AICc and ΔAICc for each candidate model were compared. We then calculated Akaike weights for all candidate models in each category (vegetation plot data and remotely sensed data) to estimate the proportion of total information explained by each candidate model. We then built a third model with the combination of the best variables from both VP data and remotely sensed data to see if the model explains the total information better than each of them alone.

ETHICAL STATEMENT

This noninvasive research was permitted by the Department of Wildlife and National Parks Peninsular Malaysia (HQ‐00026‐15‐17 to SAMS) in compliance with the laws of Malaysia, following Universiti Sains Malaysia's Institutional Animal Care and Use principles, ASP Principles for Ethical Treatment of Non‐Human Primate, and was approved by the Appalachian State University's Institutional Animal Care and Use Committee (protocol 17‐17). Access to the study area (UMFR) was permitted by the Forestry Department of Kedah (PM KT 117/2018–125/2018).

RESULTS

Agile gibbon density

The estimated area surveyed across the nine SA in UMFR was 33.4 km2, and we detected 132 H. agilis groups in the SA. The estimated daily probability of singing for an average group was 0.72, and the estimated probability of a group singing at least once in 3 days, or was 0.98. These values resulted in a mean estimated group density in UMFR of 4.03 ± 0.14 groups km−2.

The group density estimates for individual SA ranged from 3.29 to 4.99 groups km−2 (Table 2). Group density did not significantly differ across habitat types (one‐way ANOVA: F 2,8 = 1.640, p = 0.270).

TOTAL NUMBER OF GROUPS IN LOWLAND FOREST IN UMFR

We estimated that 725.6 km2 of forested habitat within UMFR was <450 m a.s.l. Since gibbon density did not differ between virgin, previously logged, and recently logged forest, we used the overall mean group density to estimate that lowland (<450 m a.s.l.) forest in UMFR may support 2927 ± 102 H. agilis groups.

Although we surveyed only areas with elevation <450 m a.s.l., this elevation range covered ca. 62% of the UMFR landscape (Figure 3), while the remaining habitat is hill (approximately 31%) and lower montane forest (approximately 7%).

VEGETATION CHARACTERISTICS ACROSS SA

The mean tree DBH, canopy cover, and tree density (number of trees/km2) did not significantly differ across forest types (Table 3); however, mean tree height was higher in logged forests (recently logged and previously logged) than in virgin forest (Table 3), and previously logged forests had a lower mean elevation than recently logged forests (Table 3).

Table 3

Comparisons of vegetation “speed plot” variables across forest types.

	Mean (±SD) for each forest type
Variable	Virgin forest	Previous logged forest	Recently logged forest	p
Mean tree DBH (cm)	21.1 (±6.0)	24.3 (±7.2)	23.4 (±7.6)	0.253
Mean canopy cover (%)	67.7 (±29.8)	74.6 (±18.5)	75.6 (±20.8)	0.863
Mean tree density (number of trees km−2)	4700 (±2438)	5967 (±1921)	5700 (±2984)	0.076
Mean tree height (m)	12.7 (±4.2)	15.9 (± 3.2)	15.5 (±4.6)	0.006
Mean elevation (m)	224.4 (±51.9)	253.4 (± 29.4)	347.2 (±48.1)	0.000

Note: Bold values indicate variables significant at p (0.05).

Abbreviation: DBH, diameter at breast height.

ECOLOGICAL PREDICTORS OF GIBBON DENSITY

The best model for vegetation plot data included only the predictor variable mean percentage of canopy cover (Table 4a). This model predicted that for with every 1% increase in canopy cover, the group density increased by ca. 0.05 group km−2, and the 95% confidence intervals surrounding the estimate excluded zero. The second‐best model included only the predictor variable mean DBH and the third‐best model included mean percentage of canopy cover and mean DBH. These models had substantially less support than the best model, as indicated by the much lower Akaike weights, and the 95% confidence intervals for mean DBH in both models included zero, suggesting that these models are not reliable.

Table 4a

Three best models for H. agilis group density among models with 15 different combinations of predictor variables from vegetation “speed plot” data.

Predictors	AICc	ΔAICc	Akaike weight	Estimates	95% CI lower	95% CI upper
Canopya	17.07	0.00	0.48	0.05	0.01	0.08
DBHb	19.28	2.21	0.16	0.12	−0.36	0.60
Canopy + DBH	19.42	2.35	0.15	0.05	0.01	0.08
				0.02	−0.32	0.36

Note: Bold values indicate best predictor model.

Abbreviations: CI, confidence interval; DBH, diameter at breast height.

“Canopy,” mean percentage of canopy cover.

“DBH,” mean DBH.

For the remotely sensed data, the best model included mean percentage of forest cover in 2000 and proportion of forest lost between 2000 and 2018 (Table 4b). The 95% confidence intervals for both variables in this model do not overlap with zero, and the Akaike weight indicates a very high likelihood that this is the best of the candidate models (Table 4b). The second‐and third‐best models had substantially higher Akaike weight (lower Akaike weight indicate better models), and confidence intervals for all variables include zero, suggesting that these models are not reliable. Including the best predictor from the vegetation plots (mean canopy cover %) in the best model of satellite imagery did not lower the Akaike weight as compared to the previous models (Supporting Information: Table S3).

Table 4b

Three best models for group density among models including seven different combinations of predictor variables from remotely sensed data (Hansen et al., 2013; National Aeronautics and Space Administration/Ministry of Economy Trade and Industry/Advanced Industrial Science and Technology/Japan Spacesystems & U.S./Japan ASTER Science Team, 2018).

Predictors	AICc	ΔAICc	Akaike weight	Estimates	95% CI Lower	95% CI Upper
Cover in 2000a	−0.094	0.000	0.987	20.639	1.596	39.681
Forest loss 2000–2018b				8.214	2.924	13.504
Cover in 2000	9.916	10.010	0.007	20.620	−1.793	43.033
Mean elevation (m)				0.000	−0.009	0.008
Forest loss 2000–2018				8.225	1.613	14.836
Forest loss 2000–2018	10.999	11.093	0.004	4.910	−0.791	10.610

Note: Bold values indicate best predictor model.

Abbreviation: CI, confidence interval.

“Cover in 2000,” mean percentage of forest cover in 2000.

“Forest loss 2000–2018,” proportion of forest lost between 2000 and 2018.

DISCUSSION

Agile gibbon density and vegetation characteristics

The mean group density of H. agilis in UMF was consistently high, and did not differ across forest types with different logging histories, which could be attributed to (1) the high habitat quality across all habitat types and (2) the adaptability of H. agilis to disturbance caused by selective logging. While we placed three SAs in forest that had been subject to selective logging only a few years before the survey (2015–2017), vegetation characteristics in these SA did not significantly differ from those in other forest types. This likely reflects the small spatial extent of the logging activity within the recently logged forests that we sampled. The significant difference in elevation between habitat types may result from spatial autocorrelation, which is the tendency for a variable to have similar values when sampled from sites close to one another (Haining, 2001). The three SAs in each forest type were placed approximately 2 km apart, so they were generally in the same landscape. The higher mean tree height, canopy cover, bigger mean DBH, and greater tree density may also be attributed by competitor removal, in which some trees can thrive after other trees were removed by logging activity. Competition among neighbors is an inhibitor for tree growth (Coomes & Allen, 2007; Zhang et al., 2017). Other factors that we did not measure, such as differences in soil quality or plant species distributions, may also have caused tree heights to differ across forest types even before the logging.

H. agilis group density was consistently high across all forest types and comparable to the highest densities reported for H. agilis in other high‐quality habitats (see Table 5; Gittins & Raemakers, 1980; Lee et al., 2015; Nongkaew et al., 2018; O'Brien et al., 2004; Yanuar, 2009). Differences in H. agilis densities across these habitats may be explained in part by competition with sympatric Symphalangus syndactylus in some habitats (Table 5). While most small apes do not have overlapping distribution ranges, Symphalangus syndactylus are sympatric with H. agilis and H. lar across most of their distribution range. Symphalangus syndactylus may reduce competition by consuming more immature foliage (Palombit, 1997), and feeding on higher and larger trees than sympatric Hylobates species (Gittins & Raemakers, 1980). However, H. agilis and S. syndactylus are both highly frugivorous, with a preference for figs, and their substantial dietary overlap renders them ecological competitors (Chivers, 1974; Elder, 2013; Palombit, 1997). Moreover, because of their ability to adopt a more folivorous diet, S. syndactylus is thought to have an advantage in higher‐elevation forests, where fruits are not as abundant (Gittins & Raemakers, 1980). In habitats on Sumatra where both S. syndactylus and H. agilis occur, the group density of S. syndactylus follows a concave quadratic pattern, in which densities decrease as the elevation transitions from lowland dipterocarp (<450 m a.s.l.) to hill dipterocarp (450–900 m a.s.l.), then increase as in lower montane and montane forest (900–2400 m a.s.l.), whereas the reverse is true for H. agilis (O'Brien et al., 2004; Yanuar, 2009). However, at UMFR, S. syndactylus does not occur, and our study was conducted entirely in lowland forest, which may explain why H. agilis densities were consistently high. Similar conditions occur in the Harapan Rain Forest in Sumatra, from which high H. agilis densities are also reported (Table 5; Lee et al., 2015). High H. agilis densities have also been reported at Sungai Dal (Table 5; Gittins & Raemakers, 1980), where S. syndactylus is absent. Bala Forest in Thailand, which is relatively close to UMFR, includes the northern limit for the distribution of S. syndactylus, but S. syndactylus occurs at very low densities, and only in the southernmost part of the forest (Nongkaew et al., 2018). Thus, the relatively high H. agilis densities in Bala Forest (Table 5; Nongkaew et al., 2018) may also reflect low or absent competition with S. syndactylus in a mostly low‐elevation landscape.

Table 5

Hylobates agilis densities reported from sites in Peninsular Malaysia and Sumatra.

		Density (groups/km2 ± 95% CI)
Site	Siamang present?	0–450 m a.s.l.	450–850 m a.s.l.	>850 m a.s.l.	References
Ulu Muda	No	4.03 ± 0.14	‐	‐	This study
Sungai Dal	No	4			Gittins and Raemakers (1980)
Bukit Barisan Selatan	Yes	0.5	1.1	0.5	O'Brien et al. (2004)
Kerinci‐Seblat	Yes	1.7 ± 0.9	3.9 ± 1.0	3.1 ± 0.4	Yanuar (2009)
Harapan Rain Forest	No	2.6a	‐	‐	Lee et al. (2015)
Bala Forest	Yes	2.27 ± 0.18			Nongkaew et al. (2018)

Abbreviation: CI, confidence interval.

Lee et al. (2015) used multiple methods to estimate density. The value reported here is based on fixed‐radius point counts. The author also reports group density estimates with ±coefficient of variation (CV) instead of 95% CI, which is not included in this table.

Lee et al. (2015) found no significant relationship between H. agilis group density and forest quality or distance to the forest boundary at Harapan Rain Forest, Sumatra. Similarly, although the most recent logging activities in the southern part of UMFR (near Ulu Legong) were in 2017 and the village is 1–2 km away from the SA (Air Legong, Bukit Gadong and Lata Gading), group densities in this area were similar to group densities in the other two forest types. Our results, along with those of Lee et al. (2015) confirm the importance of degraded forests as potential habitat for H. agilis. It is important to note that Harapan Rain Forest and UMFR both contain large forested areas with a mosaic of different disturbance histories. The results, therefore, should not be extrapolated to smaller fragments or those not adjacent to less‐disturbed habitats.

High H. agilis densities in recently disturbed forests may indicate population compression. As forests become fragmented and degraded, surviving small ape groups are forced to live within a smaller area of remaining habitat, resulting in higher density despite reduced carrying capacity. For example, Cheyne et al. (2016) found that home ranges of Hylobates funereus, Hylobates abbotti, Hylobates muelleri, and Hylobates albibarbis in Kalimantan, Indonesia, were compressed in response to disturbance. Population compression may increase intraspecific competition with negative consequences for recruitment and long‐term persistence unless habitat quality is restored. For example, S. symphalangus in fire‐damaged habitat on Sumatra had lower infant and juvenile survivorship for at least five years post‐fire (O'Brien et al., 2003) and had smaller home ranges than S. symphalangus in adjacent areas of undamaged habitat for ≥12 years post‐fire (Lappan et al., 2021). If habitat quality remains poor, degraded areas may become population sinks for small apes, resulting in long‐term population declines or local extinctions (O'Brien et al., 2003). Therefore, while the high population densities found in recently logged forests in UMFR may indicate that the selective logging did not substantially affect habitat quality, it may also reflect population compression in these areas. In our sample, forests logged 30–60 years (i.e., three to six H. agilis generations) ago also supported high H. agilis populations, indicating that selectively logged forests can comprise important small ape habitats in the long term. However, our data set cannot shed light on the processes affecting H. agilis populations between logging and forest recovery. Specifically, because acoustic surveys can only be used to estimate small ape occurrence and group densities, but not group compositions, birth rates, or survivorship of immature animals, it is unclear whether recruitment of H. agilis in the first postlogging generation is adequate to support long‐term persistence of this taxon in logged areas of UMFR. Small apes are known to adapt and persist in some degraded forests, despite these forests having fewer large trees, including fruiting and sleeping trees, and less continuous canopies (Cheyne et al., 2013). However, without conservation management to increase habitat quality and connectivity between forest fragments, populations in low‐quality habitat may become locally extinct as adults die without recruitment (O'Brien & Kinnaird, 2011). Thus, our results confirm the importance of even recently logged forests as habitat for the taxon, and the importance of protecting large landscapes, such as the Ulu Muda landscape.

We were not able to directly assess the effects of food availability or season in this study due to resource limitations, but both of these factors influence small ape's ranging behavior and activity budgets (Inoue et al., 2016), which may affect the results of acoustic surveys by introducing variation in singing probability or home range sizes, both of which we assumed to be constant over time and space. We suggest that future studies should consider food availability, survey season, and other relevant variables, to develop more refined population estimates and to better understand ecological predictors of gibbon abundance.

Total number of groups in UMFR

We estimate that 2927 ± 102 H. agilis groups remain in the lowland (<450 m a.s.l.) forests of UMFR. Because our acoustic survey method did not provide information about H. agilis group sizes in UMFR, we did not produce an estimate of the number of individuals. Previous reports of H. agilis group sizes in stable populations from Sungai Dal, Perak (4.4 individuals/group, N = 7; Gittins & Raemakers, 1980) and Sumatra (ca. 4 individuals/group, N = 49; Lee et al., 2015) indicate that group sizes of ca. four individuals per group are common. However, there is substantial variation in group sizes across habitat types (ca. 3.5, N = 15 in “Low Secondary Forest”; ca. 4.4, N = 18 in “High Secondary Forest”; Lee et al., 2015), such that we cannot assume similar group sizes at UMFR. Similarly, we surveyed only areas with elevation <450 m a.s.l., an elevation range which constitutes more than half of the landscape (approximately 62%, Figure 3) in UMFR, while the remaining habitat is hill and lower montane forest (ca. 31% and ca. 7%, respectively). Several previous studies have reported lower densities of Hylobates spp. at higher elevations (e.g., Marshall, 2009; O'Brien et al., 2004; Yanuar, 2009), so we did not extrapolate to estimate H. agilis populations in higher‐elevation areas of UMFR. It is important to note that we assumed that our sample is representative of lowland areas of the UMFR more broadly in developing our estimate of the total population. However, our SA covered only about 5% of the total lowland forest area in UMFR, and the SA were not placed randomly in the landscape, although they were distributed across habitat types. Therefore, these numbers should be interpreted with some caution.

Ecological predictors

Mean percentage of canopy cover in the vegetation plots, and mean percentage of forest cover in 2000 and proportion of forest lost between 2000 and 2018 from the remotely sensed data were included in the best models for group density, which is consistent with the findings of many other studies indicating the importance of intact forest canopies for other small apes species (e.g., Brockelman & Ali, 1987; Fleagle, 2013; Gittins, 1983; Whittaker, 2009). The density of H. lar (but not S. syndactylus) also had a positive relationship with size (DBH) and height of trees (Hankinson et al., 2021). Conversely, in the present study, we did not detect an effect of DBH or tree height. It is noteworthy that in the model including mean percentage of forest cover in 2000 and proportion of forest lost between 2000 and 2018 as predictor variables, both had positive coefficients, indicating that forests with higher canopy cover in 2000 had greater H. agilis densities, but also that greater recent forest loss was associated with greater H. agilis density. These results are consistent with the interpretation that forest loss has both short‐term and longer‐term effects. In the long term, forest loss or degradation may reduce population densities in affected areas (as indicated by the positive effect of forest cover in 2000 on H. agilis densities), but recent forest losses may be associated with population compression (as indicated by the positive effect of forest loss 2000–2018 on H. agilis densities). To further assess the effects of forest loss on population abundance of small apes in Malaysia, future research should focus on long‐term population monitoring, behavioral ecology, and demography of H. agilis and other taxa in disturbed and fragmented habitats.

Our results show the usefulness of different types of habitat data for assessment of habitat quality for H. agilis. Remotely sensed data available from open access data sets such as the Hansen/UMD/Google/USGS/NASA data set (Hansen et al., 2013) do not require substantial investments in terms of funding, time, and labor for data acquisition, so it offers many advantages. Our results confirm the value of remotely‐sensed data for analyses of the relationship between ecological variables and H. agilis densities in UMFR. However, more detailed ground‐level vegetation assessment is necessary for ground‐truthing and to provide finer‐grained information, especially for studies involving behavior or diet selection.

Models including elevation (m) as a predictor variable explained the variation in group density poorly (Table 4b). The linear regression model also failed to detect a statistically significant effect of elevation on gibbon density across the SA (F 7,7 = 0.391, p = 0.0552; Supporting Information: Table S2), which could be explained by the limited range of elevations (0–450 m) covered by this study, given that H. agilis is known to thrive in these elevations (Lee et al., 2015).

Management and research recommendations

Our results suggest that UMFR contains a globally significant population of H. agilis, and support the interpretation that the Greater Ulu Muda landscape is likely to contain the largest remaining population of H. agilis on mainland Asia. As such, it is vital that this currently unprotected landscape receives full protection as a state or national park. Degraded forests in and surrounding UMFR should be included in conservation plans for full protection (i.e., as a national/state park) and restoration of the Greater Ulu Muda landscape, plus improving canopy connectivity between disturbed and pristine areas so that these valuable H. agilis habitats will be protected. UMFR also supports globally important populations of other threatened fauna; to date, 112 mammal species have been recorded here, including charismatic megafauna such as Asian elephants (Elephas maximus), clouded leopards (Neofelis nebulosa), leopards (Panthera pardus), Malayan tapirs (Tapirus indicus), and sun bears (Helarctos malayanus), as well as more than 300 species of birds, including all ten species of Malaysian hornbills (Bashir, 2014). The habitat also provides critical ecosystem services such as water filtration and supply, carbon storage, flood mitigation, oxygen production, and climate control that are crucial for economic and social sustainability across the region as well as biodiversity conservation (Suksuwan, 2008). Therefore, protection of this landscape is in line with national and global conservation and development priorities (Ministry of Natural Resources and Environment, 2016).

Considering the lack of recent (post‐1990) abundance and distribution data for small apes in Malaysia and the intense pressure on small ape habitats, we strongly recommend systematic nation‐wide surveys of potential small ape habitat, followed by consistent population monitoring of priority habitats. New technologies and methods, such as Passive Acoustic Monitoring (PAM), and new analytical methods, such as spatially explicit capture‐recapture modelling, may improve cost‐effectiveness and increase the coverage and accuracy of population surveys. Studies investigating variation in group density along elevational gradients in Peninsular Malaysia will be particularly useful, as hill dipterocarp forest (ca. 300–750 m a.s.l) is now the most extensive forest type (Forestry Department of Peninsular Malaysia, 2016) after the depletion of lowland forests due to extensive logging and habitat conversion from the 1960s through the 1980s (Law, 2020). We also suggest that future studies consider the abundance and distribution of important food trees to infer (seasonal) food availability in small ape habitats, variables which have been shown to influence group density and home range size. We further recommend studies of the behavioral ecology and population dynamics of small apes, especially in degraded and fragmented habitats, to support design of effective methods for restoration of important habitats and improvement of connectivity between forest fragments. Better understanding of the factors affecting group density is essential for the development of sound management practices and prioritization of conservation actions, such as construction of canopy bridges to re‐establish connectivity between forest fragments, translocation/reintroduction of rehabilitated small apes to viable but unoccupied habitat fragments, and habitat restoration projects.

AUTHOR CONTRIBUTIONS

Yi Heng Pang: data curation (lead); formal analysis (equal); funding acquisition (lead); investigation (lead); project administration (equal); resources (equal); writing – original draft (lead). Susan Lappan: conceptualization (lead); data curation (equal); formal analysis (lead); funding acquisition (supporting); investigation (supporting); methodology (supporting); supervision (equal); writing – review and editing (equal). Thad Q. Bartlett: conceptualization (equal); formal analysis (supporting); investigation (supporting); methodology (supporting); writing – review and editing (supporting). Shahrul Anuar Mohamad Sah: project administration (equal); resources (equal); writing – review and editing (supporting). Nik Fadzly N. Rosely: project administration (equal); supervision (supporting); writing – review and editing (supporting). Nadine Ruppert: conceptualization (supporting); funding acquisition (supporting); methodology (supporting); project administration (lead); resources (equal); supervision (lead); writing – review and editing (equal).

CONFLICT OF INTEREST

The authors declare no conflicts of interest. This study fully adhered to the legal and ethical requirements for wildlife research as stipulated by Malaysian law.

Supplementary information.

Click here for additional data file.