The Dieng Mountains are located in Central Java Province, Indonesia (109°32'–109°56'E and 7°04'–7°13'S). The study was conducted in the remaining 90–175 km² of natural forest, which mostly covers the northern part of the mountains (Nijman and Van Balen 1998, Setiawan et al. 2012). The altitudinal range of 250–1700 m above sea level is covered by forest blocks from lowland to the submontane tropical rainforest. The area has a hilly terrain and complex riverine. The forested area consists of a mixture of natural and plantation forests dissected by a large number of secondary roads and most relatively flat areas amongst the forest patches were built as settlements and croplands. Most natural forest patches are secondary forests due to logging in the past and current non-timber extractions by local people (shade-grown coffee planting, honey collecting, firewood collecting etc.) (Setiawan et al. 2012). However, the other forest patches are topographically inaccessible. Although diverse flora and fauna inhabit the forest, especially for all primates in Java (Nijman and Van Balen 1998), this area received less protection for biodiversity conservation. All forested areas are administratively managed by Perum Perhutani and Perusahaan Perkebunan Negara (an Indonesian state-owned forestry and plantation enterprise), which are mainly responsible for forest production. Some forest blocks are protection forest, which are not converted to plantations to maintain soil fertility and prevent landslides (Nijman and Van Balen 1998). Plantation forests such as pine, teak, agathis, rubber and tea plantations are generated in the surrounding the area, mostly adjacent to the natural forest. Moreover, many natural features, such as rivers and waterfalls, have been developed as tourism parks.

This study used ecological niche modelling (Franklin 2010) to estimate the total habitat potentially occupied by Javan gibbons in the Dieng Mountains. The maximum entropy algorithm (Phillips et al. 2004, Phillips and Dudık 2008) was generated for the modelling in the MaxEnt 3.4.4 programme (Phillips et al. 2020). The model has been generated in a previous study (Widyastuti et al. 2020), but it was re-run in the present study with improved environmental variables. The total of nine independent environmental variables were involved as predictors (Table 1). The model used the 70% (214) of recorded occurrence points during the field survey in 2018 as a training sample (see Widyastuti et al. 2020). The model was run under 10 replicates, a subsample replicated run type and used the 10 percentile training presence as a threshold rule (Young et al. 2011). The value from area under the curve (AUC) of the receiver operating characteristic (ROC) was used to evaluate the model performance (Peterson et al. 2011).

The index produced from the model ranged from 0 to 1, with 1 indicating the highest probability of Javan gibbon presence. The map with the probability range was then delineated into suitable and unsuitable habitats, based on the 10% training presence logistic threshold (Young et al. 2011, Fitzgerald et al. 2018), which is 0.27. The suitable habitat was expected to be the area where Javan gibbons potentially occur. The suitable habitat was then divided equally into two types, low suitable and high suitable. The 91 remaining occurrence points were then used for validation.

The fixed-point count was used owing to its efficiency and suitability for surveying gibbons, based on its call (Brockelman and Srikosamatara 1993, O'Brien et al. 2004, Whittaker 2005, Buckley et al. 2006, Cheyne et al. 2008, Cheyne et al. 2016, Hamard et al. 2010, Kabir et al. 2021). At each sampling site, three listening posts (LPs) were located on an either linear or triangular formation on the ridge or at a higher elevation to cover a wider listening area free from disturbing noise (Brockelman and Ali 1987, Cheyne et al. 2008, Höing et al. 2013). The posts were positioned 300–500 m apart to maximise the listening area and to allow detection of groups from at least two LPs (Brockelman and Ali 1987, Kidney et al. 2016, Hankinson et al. 2021). The outmost LPs of the adjacent sampling sites were positioned ≥ 2 km away from each other to avoid counting the same gibbon group at more than one site. This study used a 1-km buffer around each LP as effective listening radius (ELR) to define the effective listening area (ELA) of each sampling site (Brockelman and Ali 1987, Cheyne et al. 2008, Hamard et al. 2010) (Fig. 1). According to the identified habitat suitability type (Table 2), four sampling sites were placed in the area where the percentage of high suitable habitat was relatively high (≥ 19%), namely high suitable sites and four sampling sites were placed in the area where the percentage of high suitable habitat was relatively low (≤ 6%), namely low suitable sites (Fig. 1).

Figure 1.  

Study area and sampling sites (1 km buffer from each listening post) within the Javan gibbon habitat suitability type in the Dieng Mountains, Central Java, Indonesia. (Sampling sites: 1: Sikesod; 2: Tombo; 3: Sawangan Ronggo; 4: Sawahan; 5: Salakan; 6: Tinalum; 7: Kalipaingan; 8: Linggo Asri).

The auditory sampling in each site was sequentially conducted from 5 October to 20 November 2021. Two observers sat at each LP in one sampling site to record the gibbon calls simultaneously. At each post in each site, each gibbon call heard between 06:00–10:00 h for four consecutive days were recorded, including the compass bearing, start and end times and number of female great calls (Brockelman and Srikosamatara 1993, Cheyne et al. 2016). The data collection was started at 06:00 h, when female gibbons in Dieng most frequently call (Geissmann and Nijman 2001). Data were not collected during rainy or stormy days and when the rain had continued for 2 h before the start of sampling, as rain and strong winds can affect negatively gibbon vocalisation (Brockelman and Srikosamatara 1993, Cheyne et al. 2007).

Each gibbon call recorded per day, based on compass bearing, was plotted into the map using ET GeoWizard in ArcMap 10.5. Imaginary lines represented the compass bearing of each call recorded. The lines were then imported to Google Earth Pro for group identification. Only calls that consist the female great call (indicating a gibbon group) were included in the analyses to avoid counting solitary gibbons. First, the bearing lines from two or more LPs were paired when the singing timeframe and number of great calls were closely similar. If the paired lines produced an intersection or triangulation, the intersected point was then identified as the estimated singing location. Finally, gibbons at two or more singing locations were identified as the same group if the locations were within 500 m of each other. Otherwise, they were identified as different groups if they sang at the same time or too close in time to be from the same group (Brockelman 2019, Brockelman and Ali 1987, Brockelman and Srikosamatara 1993, Cheyne et al. 2016, Hankinson et al. 2021, O'Brien et al. 2004) or if confirmed with other evidence (e.g. identification by sighting or fragmented by an unforested barrier).

All recorded Javan gibbon calls, including non-triangulated calls, which follow these conditions were used in the group density calculation: 1) consisted of great calls (female great call) as a representative of the family group, 2) had a compass bearing to the identified singing group location and estimated to be that particular group and 3) within an ELA (Nijman and Van Balen 1998, Geissmann and Nijman 2001). The gibbon group density data were calculated per array using the standard triangulation calculation (Brockelman and Ali 1987, Nijman 2004, Cheyne et al. 2008, Hamard et al. 2010) and an ASCR programme (Stevenson et al. 2021). However, on the basis of cases reported in previous studies (Hankinson et al. 2022), we used two alternative calculation methods for the ASCR calculation, which are explained in the Acoustic spatial capture-recapture subsection . As a result, this study compared the gibbon group densities obtained from three methods using paired t tests to analyse the differences between the methods (Hankinson et al. 2021).

In the field, visual encounters were also recorded to estimate the group size. The number of individuals seen, age class, time and coordinates were recorded. Considering that the number of groups encountered in 2021 was not robust, the group size used for estimating the population size was calculated from the visual encounters in the survey in August–December 2018, which was conducted in a longer period and obtained a larger sample (Table 4). In addition, several groups monitored in 2019 and 2020 were also summarised to compare the group size data during the period 2018–2021.

Individual densities were calculated per site by multiplying the group density with the overall group size. Individual densities were averaged for each habitat type. The population size in each habitat type was calculated by multiplying the average individual density in each habitat type by the total area for the corresponding habitat type. The population sizes in the low and high suitable habitats were then summed to obtain the overall population size.

Habitat characteristics were measured in 10 plots of 10 x 10-m vegetation plots randomly placed within the ELA of the eight sampling sites (Hamard et al. 2010, Hankinson et al. 2021). The data on the plots were measured at the same period with the population survey. The following variables for the trees with a DBH of ≥ 10 cm were recorded: 1) DBH (cm), the perimeter measured using a measuring tape and then converted to diameter; 2) tree height (m), measured using a Nikon Laser Rangefinder; 3) crown area (m²), calculated using following formula: \(A=π(\frac{N-S width}{2})x(\frac{E-W width}{2})\), 4) canopy cover (%), measured at each corner and in the middle of the plot using visual estimation by the same observer for all plots; and 5) tree density (number of trees per km²) (Hamard et al. 2010, Kim et al. 2011, Hankinson et al. 2021). The elevation for each plot was also recorded from the GPS receiver.

All habitat characteristics and group density data were tested for normality using the Shapiro-Wilktest. After the normality was confirmed for all data, the differences in habitat characteristics between the sites were evaluated using one-way analysis of variance (ANOVA) test. The correlation between the gibbon group density and each habitat characteristic was examined using the Pearson correlation test. All statistical analyses were generated in R studio version 2022.02.3 then delineated into suitable and unsuitable habitats, based on the 10% training presence logistic threshold (Young et al. 2011, Fitzgerald et al. 2018), which is 0.27. The suitable habitat was expected to be the area where Javan gibbons potentially occur. The suitable habitat was then divided equally into two types, low suitable and high suitable. The 91 remaining occurrence points were then used for validation.

The Maxent model produced the probability index which ranged from 0 to 0.96. The range of the index was then reclassified into three habitat types, which were unsuitable, low suitable and high suitable habitats (Table 2). The model revealed that the overall potential habitat (total of low suitable and high suitable habitats) for Javan gibbons in the Dieng Mountains was 100.91 km², where the low suitable habitat was larger than the high suitable habitat (Table 2). The Maxent model showed high model performance, based on an AUC value of 0.946 ± 0.007. Of the testing points, 90% were in the suitable habitat, mostly in the high suitable habitat (51.65%), which indicated the appropriate model classification.

The total survey effort covered 33.4 km² of Javan gibbon habitat (on the basis of the ELR of 1 km from each LP) across eight sites during the 32-day survey. A total of 529 call events were recorded from 24 LPs and 52 groups of Javan gibbons within all ELAs were identified. The number of groups identified in each site ranged from 2 to 10. The detection probability of Javan gibbons [p(1), the probability of calls produced on any given day] during the 3-day survey in each array, ranged from 0.38 to 0.74. The correction factors [p(m), the proportion of gibbons expected to sing at an area during the sampling period] for the 4-day survey in each site ranged from 0.85 to 0.99.

The group densities for each sampling site that were calculated by triangulation, ASCR before division and ASCR after division were within the ranges of 1.68–7.45, 0.7–5.9 and 0.46–2.68 (in groups/km²), respectively (Table 3), which showed normal distributions (Shapiro-Wilk test: P = 0.5358, P = 0.1867 and P = 0.6106, respectively). The group density calculated by triangulation and ASCR before division were not significantly different (paired t test: t = 0.85, df = 7, P = 0.422). By contrast, the density calculated by ASCR after division was significantly lower than those calculated using the other calculation methods (paired t test with triangulation: t = 5.34, df = 7, P = 0.00107; paired t test with ASCR before division: t = 4.86, df = 7, P = 0.00183). The mean density in each habitat type (between low and high suitable sites) were significantly different for all calculation methods (ANOVA test: P = 0.0195 for triangulation, P = 0.00297 for ASCR before division and P = 0.00604 for ASCR after division). Considering that the traditional triangulation method is the most common method to calculate gibbon density and the ASCR is a newly-developed statistical programme, the group densities calculated using ASCR before division were used for further analyses. The mean density was 2.15 groups/km² in the low suitable habitat and 5.55 groups/km² in the high suitable habitat (Table 3).

During the fieldwork in 2021, more than nine groups were encountered, but only nine groups from five sites could be observed and counted accurately. However, during the previous work in 2018, more extensive encounters led to the identification and confirmation of 20 groups from nine sites. Finally, we used the group size from 2018 to calculate the population size in this study, considering the larger sample size and coverage of the sampling sites. The number of individuals per group ranged from 2–7 and the mean group size was 3.95 individuals per group. A number of infants and juveniles were observed during the period 2018–2021 (Table 4).

The total suitable habitat for Javan gibbons was calculated to be 100.91 km², of which 83.41 km² was low suitable habitat and 17.50 km² was high suitable habitat. The number of Javan gibbons was estimated to be 708 in the low suitable habitat, 384 in the high suitable habitat and 1092 in the whole area of the Dieng Mountains (Table 5).

All data on the habitat characteristics followed a normal distribution (Table 6). Significant differences amongst sites were found only for elevation (ANOVA test: P = 0.04), whereas all vegetation characteristics did not significantly differ amongst sites (Table 6). No significant correlations were found between density and all habitat characteristics tested in this study (Table 7).

Although the Dieng Mountains are a heterogenous and unprotected landscape, these areas hold a large portion of the Javan gibbon population (Nijman 2004). After a decade since the last study that updated the population size, this study found a significantly higher density, group size and population size of Javan gibbons in the Dieng Mountains than any previous studies (Table 8). The mean density in the low suitable site (2.15 groups/km²) demonstrates an ideal range of gibbon density (2–5 groups/km²) (Brockelman and Srikosamatara 1993) and is still within the same range as those reported in previous studies. However, the mean density in the high suitable sites (5.55 groups/km²) was the highest amongst those from any surveys in Dieng and demonstrates a high density range (5–6 groups/km²) (Brockelman and Srikosamatara 1993) (Table 8).

The higher Javan gibbon density in the present study than in the most recent study indicates two possibilities. First, it indicates that the Javan gibbon population in Dieng has been increasing. Alternatively, the population has not been increasing, but the higher density in the recent study resulted from the difference in survey technique. Compared with the most recent study (Setiawan et al. 2012), which used the visual line-transect method, this study used a fixed-point count method, based on auditory sampling. As Javan gibbons are fully arboreal and fast, but smoothly moving primates, they can flee before being detected under visual observation inside the dense forest; thus, some gibbons might have been missed in the count in the previous study (Brockelman 2019). Our survey using auditory sampling could probably detect more gibbon groups than the previous survey, but we cannot conclude that the previous study underestimated the number of groups on the basis of this reason alone. Regarding the difference in survey method, overestimation was more unlikely in the recent study for two reasons. First, we carefully used a conservative principle to count the number of gibbon groups within the sampling area. We only separated two or more singing locations as different groups if they met the set criteria. Second, we heard and recorded high numbers of gibbon calls in the sampling sites that showed very high group densities, such as Salakan and Kalipaingan and recorded very low numbers of calls in the sites that showed very low densities, such as Sikesod, for which 0.7 group/km² was recorded (Table 3).

Without any other robust supportive evidence, we cannot conclude that the Javan gibbon density has not been increasing. On the other hand, this study shows a significantly higher group size than any other previous study of Javan gibbons in Dieng (Geissmann and Nijman 2001, Nijman 2004, Setiawan et al. 2012), which was 3.95 individuals per group (n = 20 groups) in 2018 (Table 4). In all studies, the group size was commonly determined from visual encounter; therefore, it is the best parameter for comparison between studies. The group size in the present study is an increase from that in the previous study (Setiawan et al. 2012), which reported that the mean group size in several sites in Dieng was only 2.61 individuals (n = 31 groups). An earlier study (Geissmann and Nijman 2001) also reported a lower group size from the forest near Linggo Asri in Dieng, which was 3.5 individuals (n = 15). In addition, regeneration was observed. A total of nine groups were found with offspring in 2018, of which seven had juveniles, one had an infant and one had both a juvenile and an infant. Several numbers of offspring were also found in several groups encountered in 2019–2021 (Table 4). Thus, this indicates that the first possibility is more likely. The Javan gibbon population in Dieng has been increasing or has remained at least stable during the last decade.

The population increase suggests that the reproduction rate of Javan gibbons in the Dieng Mountains during the last decade was higher than their mortality rate. The larger group size and number of offspring found in this survey prove the high reproduction rate. Healthy reproduction is supported by enough food sources. Therefore, the high reproduction rate indicates that the food tree species for Javan gibbons in the Dieng Mountains were plentiful, although more work is needed to confirm this. Furthermore, the availability of habitat space for the territory of new groups is also important for the population growth of Javan gibbons. On the basis of the Javan gibbon home range size of 12–36 ha (Iskandar 2007, Malone 2007, Kim et al. 2011, Maya 2013), the identified suitable habitat of 100 km² (equal to 10,000 ha) could roughly support 400 gibbon groups (with the assumption of a 25-ha home range size), but the estimate in this study was only 276 groups. This means that the Dieng Mountains have provided enough habitat space for new gibbon groups in the last decade and can still support the population growth in the future.

The Dieng Mountains are not protected under a conservation area and consist of heterogenous land-use and land-cover. This situation has led to various levels of disturbance that threaten gibbons and their habitats. However, this study suggests that the rate of mortality or disappearance due to illegal hunting of Javan gibbons in the Dieng Mountains has been lower than their natality rate in the last decade. This might be a result of the wide range of conservation efforts to reduce the threats at the grassroots level by various parties, mostly by a local non-government organisation (NGO; SwaraOwa), during the last decade. Several long-term conservation activities have been implemented, including livelihood development, environmental education, community awareness and ecological research (Setiawan et al. 2009, Herdiansyah and Setiyono 2019, Setiawan et al. 2020). In particular areas, such as near the Sokokembang forest, there is an indication that the local people have become more aware of the need to protect these endangered primates from illegal hunting. Furthermore, no case of Javan gibbon illegal hunting near Sokokembang has been reported during the last decade (Hendriati et al. 2020).

The IUCN Species Survival Commission Primate Specialist Group's Section on Small Apes recommended the recently-developed ASCR for gibbon population studies, as it is the most accurate way to analyse acoustic data (Cheyne, unpublished data; Hankinson et al. 2021). However, standard triangulation was still used in this study to calculate the gibbon group density because it has an important basic concept in density estimation using vocal count data (Brockelman and Ali 1987, Hankinson et al. 2021) and has been proven to be reliable in several surveys (Brockelman and Srikosamatara 1993, Nijman 2004, Cheyne et al. 2008, Höing et al. 2013, Gilhooly et al. 2015). Our results show that dividing the D from ASCR by the number of survey days, which was four and by the probability of gibbon calls in a day [p(1)] resulted in a significantly lower density than the original D before the divisions. Moreover, the densities obtained before the divisions were most likely to represent the number of calls heard and the groups identified during the survey. Furthermore, the densities obtained using ASCR before the divisions were within the same range as the densities obtained using the basic method, standard triangulation. For this reason, we used the densities from the ASCR before the divisions for the population size calculation and correlation analyses with habitat characteristics. The finding that the D of ASCR before division is more realistic than that of ASCR after division was similar to the findings of previous studies for lar gibbon (Hylobates lar) and siamang (Symphalangus syndactylus) in north Sumatra (Hankinson et al. 2021; personal communication, Hankinson 2022), but not for Thomas' langurs (Presbytis thomassi) (Hankinson et al. 2022). By contrast, in accordance with the instruction on the manual web-page of ASCR to divide the D by the number of survey days and p(1), the langur group densities appeared to be realistic and comparable with the results of a previous studies (Hankinson et al. 2022). This study indicates that the ASCR method is a useful modelling tool for calculating gibbon or primate density, based on calls, although the steps to divide the D by the number of survey days and p(1) may give different results in different species at different study sites. Therefore, in order to calculate the gibbon group density, we suggest to carefully consider the D obtained using the ASCR method, both before and after division, by comparing it with other Ds that were obtained using the other calculation method, such as the standard triangulation or with other Ds from previous studies.

The density was converted into population size using two separate extrapolations, based on the habitat suitability model for the low and high suitable habitats. This approach revealed that the gibbon densities in the two habitat types were significantly different. Thus, extrapolating by averaging the densities from all sites will result in bias estimation. The previous Javan gibbon population studies used separate extrapolations, based on altitudinal range (e.g.,] Kappeler (1984), Asquith et al. (1995), Nijman and Van Balen (1998)). However, the latest studies also support that elevation was one of the most important predictors of the presence of Javan gibbons, but there were other environmental variables which were also important, such as anthropogenic disturbances (Widyastuti et al. 2020). In the current condition in Dieng, not all lowland forests, which provide better support for higher gibbon densities than high-elevation forests, have good-quality vegetation and are free from anthropogenic disturbances. In our study, the group density in Sawahan, which is lowland forest, was much lower than those in other lowland sites. Therefore, the combination of multiple variables, rather than elevation alone, was taken into account in the extrapolation.

Our study did not find any correlation between the gibbon density in the Dieng Mountains and any of the habitat characteristics tested (Table 7). Although the density significantly differed between the habitat suitability types, it was not influenced by the vegetation structure or forest elevation. All vegetation characteristics did not show significant differences between the sites (Table 6). This indicates that the forest structures throughout the Dieng Mountains were similar and had no effect on the density variation. The same pattern of vegetation structure and no effect on gibbon density were also found in another study on western hoolock gibbons (Hoolock hoolock) in Namdapha National Park, India (Ray et al. 2015).

Other vegetation characteristics such the availability of food trees might have more significant influences on gibbon density. However, we could not examine this because of the lack of Javan gibbon food data and limitation in tree species identification. The tree species were identified with non-standardised local names by the experienced local guide; thus, the names given could differ between the sites. Therefore, the study of Javan gibbon dietary ecology that includes an analysis of vegetation composition in the Dieng Mountains is crucial for further research. Alternatively, as our study used multiple variables in the habitat type classification, anthropogenic disturbance or climatic factors could also be the major causes of the gibbon density variation. Statistical analyses that take into account multiple variable combinations to investigate the major causes of the density variation are recommended for further study.

This study found high gibbon group densities and estimated a large number of Javan gibbons inhabiting the Dieng Mountains. However, threats to Javan gibbons remain, the risk of population decline is still high and even local extinction is still possible. Long-term population monitoring should be strengthened along with other conservation programmes to ensure the long-term future of this endemic species. Long-term population monitoring is important to monitor the trend of the population size and also to detect the probability of local extinction in the early stage, which could be caused by hunting or deforestation (Vu et al. 2020).

Although our study indicates that the conservation efforts during last decade have led to the increase the Javan gibbon population in Dieng, more efforts are crucial to ensure the long-term future of this local population. The habitat degradation in some locations has been and will be a potential problem for Javan gibbons. For example, the forest surrounding Sawahan has been degraded because of the priority for unsustainable agroforestry, which has resulted in a low gibbon density. Furthermore, some forest patches in Kalipaingan and Linggo Asri hold high gibbon densities, as they are topographically inaccessible, but the remaining areas are prioritised for unsustainable agriculture. Conservation programmes, such as environmental education and community development for the villages surrounding these locations, should be strengthened to raise awareness and shift livelihoods to more sustainable ones. In addition, a new threat has arisen. In the past 5 years, the local people surrounding this Javan gibbon habitat, supported by the local government, developed many natural attractions which successfully attracted massive numbers of tourists. However, although this provides an alternative strategy to improve the local economy for the benefit of the local people, this poses a serious threat to Javan gibbons and their habitat, unless wise and careful management is employed. Gibbon watching as a nature-based tourism is a good tourism programme, which has been initiated in Petungkriyono (Supriatna et al. 2022). This provides an alternative sustainable tourism to this district and could be combined with sustainable or wildlife-friendly products of forest commodities to improve the local economy.

The results of this study strongly support that the Dieng Mountains are an important habitat for Javan gibbons and hold a significant proportion of the total population of gibbons in Java. Since the previous study, the forest in the Dieng Mountains has been suggested to be designated as a protected area (Nijman and Van Balen 1998), but this appears not to be a good solution for such a heterogenous landscape and complicated condition, as explained earlier. Many people live adjacent to the forest, who need to fulfil their economic well-being. Furthermore, a win-win solution should be attempted to allow the regulation of both interests, which are conserving the biodiversity in the Dieng Mountains, including Javan gibbons and improving the local economy. Since 2018, the local NGO (SwaraOwa) and local government have proposed a 5173.80 ha forest area in Petungkriyono for a collaborative forest management called Kawasan Ekosistem Esensial, or essential ecosystem area (EEA), to protect the habitat of Javan gibbons and the biodiversity in the Dieng Mountains. The EEA can be managed collaboratively by a forum of various stakeholders mainly to protect the biodiversity and ecosystem in the area and to manage the ecosystem sustainably for economic development (Sahide et al. 2020, Rosdiana et al. 2022). Although it only involves a small portion of the natural forest in the Dieng Mountains, it is a good decision to start with and can be expanded to other forest patches in the future.