Borneo's amphibians face an uncertain future due to high levels of forest degradation and a paucity of data for effective conservation management. Several studies identified strong species–habitat relationships in pristine and conventionally logged forests. However, these studies did not account for detectability or habitat associations in sustainably managed forest reserves. Here, we determined detectability and species habitat relationships in stream amphibians within the Deramakot forest reserve, a reduced impact logging concession in the Malaysian state of Sabah, northern Borneo. We analyzed data for 10 stream species collected along 32 stream transects. An occupancy modeling framework was used to determine the climatological, temporal, and environmental covariates associated with detection and occupancy probabilities. We identified high variability in detection probability between species, including significant associations with moon phase (six species), time since sunset (five species), humidity (five species), rainfall (four species), and temperature (three species). Stream slope and volume provided by far the best predictors of occurrence, with significant positive or negative associations with the occupancy of six species each. These associations were more similar to those found in pristine compared with conventionally logged habitats. The highly variable detectability associations within our amphibian community suggest a level of temporal separation in regard to activity and breeding phenology in these species. This stresses the importance of accounting for detection probability via surveying sites across varying climatic/temporal conditions to obtain a representative sample of amphibian communities in pristine and disturbed tropical forests.
The global amphibian extinction crisis is well documented (Stuart et al., 2004). Although this decline is poorly documented in Southeast Asia, high land conversion rates are likely to severely impact understudied and diverse rainforest amphibians (Rowley et al., 2009; Bickford et al., 2010). Borneo boasts among the highest amphibian biodiversity and endemism in the region (Goutte et al., 2017; Inger et al., 2017). However, the ecology of many species is poorly understood, and the environmental factors impacting species occurrence have seldom been systematically analyzed. Several authors (Keller et al., 2009; Konopik et al., 2015; Goutte et al., 2017) addressed this discrepancy via determining stream anuran habitat associations in northern Borneo. Their research found that environmental heterogeneity, specifically variations in stream volume and slope, provided the best predictors for anuran occurrence in pristine habitats (Keller et al., 2009; Goutte et al., 2017). However, in oil palm and conventionally logged forests, only stream slope was identified as a predictor of anuran occurrence (Konopik et al., 2015). Because amphibian species–habitat relationships can be dependent on the degree of disturbance (Ernst and Rödel, 2005), an understanding of these relationships across disturbance gradients is essential for effective conservation management.
Borneo has some of the highest deforestation rates in the world (Gaveau et al., 2014), with logging concessions covering the majority of Borneo's remaining forests (Edwards et al., 2014; Struebig et al., 2015). Forest disturbance and biodiversity losses in these reserves are entirely dependent on the management regime, with reduced impact logging (RIL) methods using reduced harvesting rates (<30 m3 timber per ha) and best practice harvesting techniques (e.g., directional felling, reduced skid trail construction, limitations on extractable timber; Pinard et al., 1995). As such, RIL results in lower biodiversity loss (Bicknell et al., 2014) and 50% less damage to remaining forests than conventional logging techniques (Sist et al., 1998). Malaysian and Indonesian commitments to increase sustainably certified forest reserves has led to a sharp increase in RIL use throughout Borneo (Sabah Forestry Department, 2009; Ellis et al., 2019). Unfortunately, no studies have focused on Bornean amphibians within sustainably managed RIL reserves.
In addition to the lack of studies in RIL sites, previous studies investigating amphibian–habitat associations on Borneo did not account for detectability, potentially undermining the general applicability of their findings. Underestimated nondetection can increase negative and positive bias in habitat covariate associations with occupancy (Gu and Swihart, 2004). Furthermore, studies in the United States (Homyack et al., 2016; Guzy et al., 2019) and Brazil (Ribeiro et al., 2018) identified highly variable detection probabilities between amphibians in the same community. The climatological factors associated with detectability in these communities were linked to calling, breeding behavior and phenological cues, suggesting some form of temporal activity separation between co-occurring species. As a result, accounting for the various detection probabilities of amphibian community members is essential for determining the covariates associated with breeding phenology and for determining site occupancy where a species is imperfectly detected (Mackenzie et al., 2002).
In this study, we therefore assessed 1) which covariates contain the strongest associations with species-specific detection probabilities and 2) whether stream anuran species–habitat associations within an RIL forest were similar to those found in pristine forest (Keller et al., 2009; Goutte et al., 2017), conventionally logged forest, and oil palm plantation (Konopik et al., 2015) streams.
Materials and Methods
We conducted this study within the Deramakot Forest Reserve (5°14′28″N, 117°19′36″E, WGS 84), central Sabah, Malaysian Borneo (Fig. 1). The 550-km2 reserve comprises predominantly hilly, lowland dipterocarp forests (50–350 m a.s.l.). Throughout Deramakot, RIL techniques are used in accordance with the Forest Stewardship Council guidelines (e.g., <30 m3 timber per ha, directional felling, reduced skid trail construction, pre- and postharvest planning). In addition, RIL is conducted only within a 291–770-ha area per annum (0.52–1.4% total reserve area). Deramakot was the first Southeast Asian forest certified by the Forest Stewardship Council in 1997 and has been credited for its sustainable forest management (Lagan et al., 2007). We selected 32 streams sites of varying structural dynamics, spaced 150 m–12 km apart, for stream amphibian sampling to provide a representative sample of Deramakot's variable stream habitats (Fig. 1).
The 32 streams were surveyed between March and August in 2017, 2018, or both years (n = 9 in 2017, n = 19 in 2018, n = 4 in both years). Four streams were surveyed in both years to improve the estimates of occurrence due to possible fluctuations in amphibian populations between years. All streams were surveyed on three to eight occasions, with at least 5 d between each survey. We conducted standardized visual and acoustic transect sampling to determine community composition. This method provides cost-effective, repeatable, quantitative data for amphibians, while maintaining a low impact on the study organisms (Rödel and Ernst, 2004). Each stream contained one transect located either along the stream center (stream width <2.5 m) or along the streambank (>2.5 m), so at least half of the transect area is occupied by streambank regardless of width. All stream transects were 100 m long and 5 m wide and were divided into 10 adjacent 10 × 5-m plots (plots per Keller et al.  and Goutte et al. ). These plots were searched for amphibians from the forest floor to a height of 2 m. All surveys were conducted between 1830 and 2300 h, the time period that encompasses the peak calling/activity times for the majority of Borneo's stream anurans (Konopik et al., 2015; Inger et al., 2017). On each survey, two persons searched a transect for 30 min each (3 min per person per plot). All amphibians encountered visually and acoustically within the 5-m transect width were recorded and identified to species following Inger et al. (2017).
Habitat and Detectability Covariates.—
Covariates assumed to affect detectability in Bornean amphibians (Inger et al., 2017) and that have direct associations with amphibian reproductive cues, activity, and calling behavior were recorded (Mackenzie et al., 2002; Allentoft and O'Brien, 2010; Grant et al., 2012; Homyack et al., 2016). Time since sunset (TSS) (0–4.11 h) and current moon phase (0–100% of lunar disc visible) were recorded at the start of each transect visit. We collected daily rainfall, temperature, and humidity data from a weather station located within 400 m–13 km of all transects (Fig. 1). Habitat covariates known to affect Bornean stream amphibian occurrence were collected: covariates associated with egg laying and tadpole habitats (stream volume, speed, slope, and siltation; Inger, 1966; Inger et al., 2017); covariates associated with sheltering/foraging sites for adult amphibians (substrate size, leaf litter cover, and canopy closure; Keller et al., 2009; Goutte et al., 2017); and covariates associated with anuran calling sites (understory density and vegetation height density; Goutte et al., 2017). All covariates were measured every 10 m along transects (recorded from the start and end of each plot). Because adjacent plots contained similar habitat covariate values, we included a random effect covariate for each transect (see Analysis). We measured stream volume by multiplying the plot length (10 m) by the stream's wet width and depth at each 10-m interval. We measured stream speed (meters per second) by timing how long a flotation device traveled along a 1-m stream section. Speed and volume were measured during each survey, and the average of all values was taken. All subsequent covariates were collected once. We measured stream slope as the height difference along the stream every 10 m. We ran a 10-m tape from the start to the end of a plot, with a spirit level used at the plot end to determine the height of the slope; we then measured this height by using a 2-m measuring rod. We visually assessed siltation within a quadrat located at the bottom of the stream center in four categories: 0–25, 26–50, 51–75, and 76–100% siltation cover. Average substrate sizes from five categories: <0.5, 0.5–1.5, 1.5–5, 5–50, and 50–200 cm, were visually assessed within two quadrats (1 × 1 m) along the left/right side (bank substrate) and center (stream bottom substrate). We assessed leaf litter cover within each of these quadrats in four categories: 0–25, 26–50, 51–75, and 76–100%. These categorical covariates were then converted to a continuous scale. We measured canopy cover by using images taken with a Nikon Coolpix S33, held 2 m from the ground pointing north facing the canopy. Images were taken on the left, center, and right edge of a plot. Images were converted to black and white by using Habitat-Net (Abrams et al., 2019). Next, we divided the number of black pixels by total pixels in the images to generate canopy coverage estimates. Understory density measurements were taken in two categories: understory density and vegetation height density. We measured understory density by using an image taken with the above-mentioned camera from the transect center of a 1 × 1.5-m red sheet held 2.5 m from the left transect edge. Because images taken on the right side in larger rivers would fall within the stream, no right-side images were taken. Images were cropped to show only the extent of the red sheet and converted to black and white manually. We estimated the percentage of flysheet covered by vegetation. We recorded vegetation height density by using a modification following Keller et al. (2009), with vegetation of different heights (0–30, 31–100, and 101–200 cm) counted within each quadrat in a plot.
We restricted our analysis to species occupying >16 plots (5% of sites occupied) to ensure that sufficient numbers of records and sites were available to distinguish which covariates impacted detection and occupancy and to prevent zero-inflated nonconvergence (Mackenzie et al., 2005). To determine species' habitat associations and detectability, single-species, single-season occupancy models were created using the package unmarked version 0.12-3 (Fiske and Chandler, 2011). This method estimates occupancy where species may be detected imperfectly; a variable number of site visits is conducted and allows occupancy (psi) and detectability (p) to be modeled as a function of covariates (Mackenzie et al., 2002, 2017). Because of the multiyear sampling scheme, we determined the effect of year on site occupancy by “stacking” the data, representing the four sites sampled in 2017 and 2018 twice. Before model selection, all continuous covariates were scaled to have a mean of 0 and variance of 1. All covariates were then tested for collinearity by using Spearman's rank correlation in the package Hmisc version 4.2-0, with substantially correlated covariates (coefficient > 0.7) excluded from subsequent analysis (Dormann et al., 2013). Because the vegetation height density covariate is associated with the density of available perches for calling, we only included it within models for species that exhibit perching and calling behavior. We conducted model selection via a two-part process in the package unmarked version 0.12-3 (Fiske and Chandler, 2011). Initially, we determined the optimum detection model for each species by using all possible additive combinations of the detection covariates, including a global model (all covariates) and a null model (no covariates). Subsequently, all possible additive combinations of the habitat covariates were included in the occupancy models, including a global model and a null model. No interaction covariate combinations were included. We defined covariate combinations resulting in the optimum model via Akaike's Information Criterion model selection in unmarked (Burnham and Anderson, 2004). Because our dataset consists of 320 plots nested within 32 transects (one transect = 10 plots), we included a random effect for each transect in the final models to account for the nested spatial effect between plots in the same transect. Because unmarked lacks this function, we used a Bayesian framework in JAGS (Plummer, 2003) via R using the package rjags version 4.3.0 (Plummer, 2014) to incorporate a random effect into the optimum occupancy model for each species (Appendix 1). We ran three parallel Markov chains with 50,000 iterations of which we discarded the first 10,000 as burn-in and applied a thinning rate of 20. We assessed convergence using the Gelman–Rubin statistic (Gelman, 2004). We conducted all analysis in R 3.5.2 (R Core Team, 2019).
Across the 32 stream transects, we conducted 146 transect surveys (mean surveys per transect: 4.6; SD: 1.9). We recorded 24 amphibian species from six families, with only 10 species recorded in 16 or more plots (Table 1). Because all other species were recorded in no more than five plots (1.56% naïve occupancy), we restricted analysis to these 10 taxa. Because Limnonectes kuhlii and Chalcorana raniceps/megaloensa are species complexes, we herein refer to them as Limnonectes cf. kuhlii and Chalcorana cf. raniceps, respectively (Inger et al., 2017).
Spearman's rank correlation tests for habitat and detection covariates identified no collinearity; therefore, all covariates previously listed were included in the occupancy models. Optimum occupancy models for each species yielded a range of detection and occupancy covariates (Table 2). Moon phase is most represented in optimum detectability models (six species), followed by humidity and TSS (five species), followed by temperature and rainfall (four species). Higher moon phases significantly improved detection probability for four (Leptobrachella parva, Alcalus baluensis, Limnonectes cf. kuhlii, and Leptobrachium abbotti) of six species (two significant negative associations: Chalcorana cf. raniceps and Limnonectes leporinus) (Fig. 2). Rainfall exhibited strong, significant negative associations with detectability in the optimum detection models of four species (Leptobrachella parva, Chalcorana cf. raniceps, Meristogenys orphnocnemis, and Staurois guttatus) (Fig. 2). Humidity, TSS, and temperature exhibited significant, yet weak, positive and negative associations with detectability (Fig. 2). In one species (Limnonectes ingeri), covariates did not explain detection probability better than the null model (no covariate effects).
Although a range of habitat covariates were selected in the optimum occupancy models, four covariates—stream slope (eight species), stream volume (eight species), canopy closure (five species), and substrate size (four species)—are most represented. All species except Chalcorana cf. raniceps exhibited a strong significant association with either stream slope or stream volume (six species each). Four species (Limnonectes ingeri, Leptobrachium abbotti, Pulchrana picturata, and Limnonectes leporinus) exhibited significant negative associations with stream slope, whereas two species (Leptobrachella parva and Alcalus baluensis) exhibited the opposite associations (Fig. 3). Four species exhibited significant positive associations with stream volume (Meristogenys orphnocnemis, Pulchrana picturata, Limnonectes leporinus, and Staurois guttatus), compared with significant negative associations in two species (Limnonectes cf. kuhlii and Leptobrachium abbotti). Vegetation height density covariates were included in the model selection process for six perch-calling species: Alcalus baluensis, Chalcorana cf. raniceps, Leptobrachella parva, Meristogenys orphnocnemis, Staurois guttatus, and Pulchrana picturata (Tables 1, 2). However, only two strata categories (0–30 and 31–100 cm) were ever selected within optimum occupancy models (0–30 cm, three species; 31–100 cm, two species). Within optimum models, year was only significantly (negatively) associated with the occupancy probability of Limnonectes leporinus (Table 2).
We found different species-specific associations between detection and habitat covariates in the stream amphibians of this study. Although our habitat associations mostly conform to those of previous studies (Keller et al., 2009; Konopik et al., 2015; Goutte et al., 2017), we found highly variable detectability associations between species. Four species exhibited significant positive associations with detection probability and higher moon phases. Several studies have identified either positive or negative associations with moon phase and amphibian activity/detectability (Grant et al., 2012; Underhill and Höbel, 2018). These associations are assumed to be driven by predator avoidance, foraging activity, and breeding/calling behavior (Grant et al., 2012). Reproductive behavior associated with moon phase may be the result of increased visibility of conspecifics or correspond to the synchronization of temporally separated explosive breeding events. The two species with the strongest positive associations with moon phase (Leptobrachella parva and Leptobrachium abbotti) are both explosive breeding megophryids with high calling intensity (Inger et al., 2017). However, during our study no cloud cover data were recorded, potentially undermining visibility provided by higher moon phases. Regardless, given the number and strength of significant associations with this covariate, our observations suggest that moon phase could be an important cue for breeding in these stream species.
Increased rainfall is generally linked to increased activity and abundance in amphibians (Allentoft and O'Brien, 2010); however, within our study, four species exhibited strong, significant negative associations with this covariate and detectability. All our study species breed within stream networks; thus, water availability, a key factor in increased amphibian breeding activity, is likely not a confounding issue. Each species also has specific egg deposition site habitat requirements for successful egg and tadpole development (Inger 1966; Inger et al., 1986, 2017). High rainfall shifts both the acoustic qualities of the surrounding soundscape (Lengagne and Slater, 2002) and the waterbody dynamics (e.g., stream size, speed, aquatic microhabitats) over short time periods (Townsend, 1989). Reduced activity after heavy rains could be due to avoiding suboptimal calling time, thereby reducing the risk of egg clutches and young tadpoles drifting into suboptimal habitats. Reduced calling activity and breeding behavior after heavy rainfall have been recorded in a variety of stream anurans in West Africa (Rödel and Bangoura, 2004), Japan (Fukuyama and Kusano, 1992), and Brazil (Hatano et al., 2002). Similarities in responses to rainfall could therefore be a universal trait in stream frog guilds that use sensitive riparian tadpole/egg laying habitats, particularly in streams where desiccation risk/water persistence is not a confounding issue.
Increased humidity resulted in significant positive (three species) and negative (two species) associations with detectability, whereas temperature exhibited two significant positive and one negative association. Humidity and temperature were found to positively and negatively affect amphibian detectability in other amphibian communities (MacKenzie et al., 2002; Homyack et al., 2016; Guzy et al., 2019). These studies suggested temperature and humidity are cues for amphibian breeding phenology, conforming to previous research of anuran phenological cycles (While and Uller, 2014; Ficetola and Maiorano, 2016; Canavero et al., 2019; Chmura et al., 2019). Given the converse associations with these covariates, it is possible that this represents separation in phenological cycles between our study species.
Previous authors suggested that Bornean stream amphibian calling activity peaks between 1800 and 2130 h (Konopik et al., 2015; Goutte et al., 2017; Inger et al., 2017). In our study, two species (Limnonectes leporinus and Staurois guttatus) confirmed this peak time, with decreasing detection probability after 1830 h (TSS), whereas three species (Leptobrachella parva, Leptobrachium abbotti, and Alcalus baluensis) exhibited increasing detection at later times. Within anurans, the acoustic niche hypothesis suggests that call partitioning must occur spatially, temporally, or with call dynamics for multiple species to occupy the same acoustic environment (Krause, 1993). Several studies of call partitioning in African and Neotropical anuran communities support this theory but found that calls were predominately separated spatially and by frequency dynamics (Sinsch et al., 2012; Villanueva-Rivera, 2014; Protázio et al., 2015), with few species in the same community exhibiting temporal separation (Lüddecke et al., 2000; Gottsberger and Gruber, 2004). Our results support this dynamic, as detectability between species is positively and negatively associated with TSS and other covariates (moon phase, humidity, and temperature). In addition, our species varied spatially at the macroscales (habitat associations) and microscales (perch heights) (Table 1). This suggests that although temporal separation does occur in three of our study species, this separation is dwarfed by spatial separation (habitat associations) and by variation in detectability associations with other covariates.
Although our detection results provide new insights into these species' behavior, habitat associations indicated that stream slope, volume, canopy closure, and substrate size were the best predictors of stream amphibian occurrence. All 10 stream amphibians are classified as stream breeding, with calling behavior and tadpole stages tied to stream microhabitats, which is supported by our model results (Inger, 1966; Inger et al., 1986). Species exhibiting positive associations with stream slope had tadpoles associated with mobile, loose-substrate habitats (Inger, 1966; Inger et al., 1986) commonly found in these streams (Winemiller et al., 2008). Those associated with large flat streams (Pulchrana picturata and Limnonectes leporinus) exhibited tadpole stages requiring side pools with heavy aquatic leaf litter (Inger, 1966; Inger et al., 1986, 2017). Three species (Limnonectes cf. kuhlii, Meristogenys orphnocnemis, and Staurois guttatus) exhibited positive associations with substrate size. These species also have tadpoles with adaptations such as gastric suckers and long tails (Meristogenys; Shimada et al., 2007) and dorsoventral compression (Limnonectes cf. kuhlii; Inger, 1966) for living in rocky, turbulent waters, or they use torrent-adjacent pools (Staurois; Inger et al., 2017). Canopy closure was present within five optimum species models and exhibited significant positive and negative associations with two species, Alcalus baluensis and Limnonectes ingeri, respectively. Limnonectes ingeri is primarily a secondary forest species, abundant in disturbed habitats, whereas Alcalus baluensis is detected in primary/undisturbed forests (Inger et al., 2017). These results conform to similar Bornean anuran studies that linked canopy cover to forest disturbance and subsequently amphibian occurrence/abundance (Konopik et al., 2015). In addition, research in other regions found both positive and negative anuran responses to canopy cover, dependent on species disturbance tolerance/breeding ecology (Werner et al., 2007; Provete et al., 2014).
Although these covariate associations predominantly conformed to the results of Keller et al. (2009), Konopik et al. (2015), and Goutte et al. (2017), there were minor exceptions. Leptobrachella parva, Alcalus baluensis, and Leptobrachium abbotti exhibited the opposite associations with slope and stream size from what was observed in these studies. The cause of this disparity could be due to the inclusion of detection probability, which may have provided greater resolution for determining their habitat associations. Increased likelihood of nondetection can lead to significant bias in the estimation of habitat associations (Gu and Swihart, 2004). Although this may explain the disparity in associations between these species, all other species–habitat relationships appear to be more similar to primary forests (stream slope and volume; Keller et al., 2009; Goutte et al., 2017) than those in conventionally logged and palm oil plantation streams (stream slope only; Konopik et al., 2015). Amphibian species–habitat relationships in other areas can exhibit significant variation, often dependent on the degree of site level disturbance, the species' resilience, and the spatial scale of sampling (Ernst and Rödel, 2005; Stoddard and Hayes, 2005; Ernst et al., 2012). The level of disturbance within our RIL site exceeds that of pristine habitats but is far below that of conventionally logged sites (Pereira et al., 2002) or palm oil plantations (Konopik et al., 2015). Given the proliferation of logging reserves throughout Sabah, RIL areas may provide landscape-level refuge for these species outside of the comparatively small protected area network (Gaveau et al., 2014).
Our results demonstrate that detectability is highly variable between our study species. Because of their co-occurrence, we speculate that some form of separation exists between their breeding cycles and activity times. Furthermore, stream amphibian habitat associations within this RIL reserve are similar to those of primary forest sites (Keller et al., 2009; Goutte et al., 2017). As such, the protection of heterogeneous stream habitats within pristine and logged forests is paramount for the conservation of Borneo's stream anuran diversity.
We thank V. Vitalis, who served as the research assistant for the entirety of this project, providing invaluable assistance. We are indebted to the Sabah Forestry Department, particularly J. Kissing and P. Lagan who facilitated access to the study site, field accommodation, and a wealth of data. The Auckland Zoo, Columbus Zoo, Museum für Naturkunde, Berlin and the Leibniz Institute for Zoo and Wildlife Research, Berlin funded parts of this research. Sabah Biodiversity Council granted the access licence to conduct this research (permit JKM/MBS.1000-2/2 JLD.7 (63)). The manuscript profited from the constructive critique of two anonymous reviewers.
Description of the Global Single-Species Occupancy Model (including All Occupancy and Detection Covariates) including a Random Effect in JAGS
zj = true occupancy state (0 or 1) of the species of interest (of the 10 stream amphibians of this study) at plot j
ψj = respective occupancy probability
α = intercept of the logit-linear predictor of occupancy probability
β1–β11 = coefficients for environmental covariates.
In our model, yjk are the observations (0 or 1) of the species of interest at site j on occasion k; transect-specific random effects are drawn from a normal distribution with study-wide means and SDs; jk are the respective detection probabilities; α.p is the intercept of the logit-linear predictor of detection probability; and β.r1–β.r5 is the effect of the five detection covariates on detection probability.