Ecology of the Granular Spiny Frog Quasipaa verrucospinosa (Amphibia: Anura - Dicroglossidae) in Central Vietnam

Abstract

We conducted a large-scale assessment at 35 primary forest sites and 42 secondary forest sites in Bach Ma National Park, central Vietnam, using the detection/non-detection data for each site over multiple visits, to quantify the site proportions that were occupied by granular spiny frogs (Quasipaa verrucospinosa). We additionally investigated the effect of site covariates (primary versus secondary forests) and sample covariates (temperature, humidity, and precipitation) to examine the environmental needs that may be incorporated for conserving rain forest amphibians in Vietnam. From the best model among all candidate models, We estimated a site occupancy probability of 0.632 that was higher than the naïve occupancy estimate of 0.403 and a 57% increase over the proportion of sites at which frogs were actually observed. The primary forest variable was an important determinant of site occupancy, whereas occupancy was not associated with the variable of secondary forest. In a combined AIC model weight, the detection model p (temperature, humidity, precipitation) included 90.9% of the total weight, providing clear evidence that environmental conditions were important sample covariates in modeling detection probabilities of granular spiny frogs. Our results substantiate the importance of incorporating occupancy and detection probabilities into studies of habitat relationships and suggest that the primary forest factor associated with environmental conditions influence the occupancy of granular spiny frogs.

1. Introduction

Studies on site occupancy of various species and their spatial patterns have been conducted in recent years to inform and develop amphibian conservation programs [1, 2, 3, 4]. Data that were commonly adopted may include forest type, land ownership (e.g., state or private), soil class, and the proximity of a resource (e.g., streams and wetlands) to the sites occupied by a targeted species [4, 5, 6]. However, many studies of metapopulation dynamics that seek to understand for the factors that determine whether a species will exist at a location [7, 8, 9]. Large-scale monitoring programs for amphibian species [10, 11, 12] often relied on remotely sensed data (data on amphibians that has been gathered using biophysical variables derived from moderate resolution imaging spectroradiometers or from remote-sensing instruments on satellites) to depict spatial models in habitat occupancy [13, 14, 15]. Ignoring detectability may lead to biased estimations of site occupancy [16, 17, 18] and studies of habitat occupancy are often hampered by imperfect detectability for the species [1, 19, 20, 21, 22].

Previous studies suggest that it is preferable to use a sampling method that involves multiple visits to sites (or patches) during the appropriate season in which a species can be detected [3] and the proportion of sites that are occupied by the species is assessed in the face of imperfect detection [21, 22]. In these cases, sampling sites may represent separate habitat patches in a dynamic context of metapopulations or sampling units (quadrats) regularly visited as part of a large-scale monitoring program [6] because the presence or absence of a species from a collection permits inference to the entire region of interest. Detection and nondetection models using multiple visits to each site permit assessment of detection probabilities and interesting parameters, including determining the ratio of sites occupied by the target species.

Human activities such as timber harvest, land conversion to monoculture crops or other developments, and manipulation of natural waterways have heavily eroded tropical primary rain forests and secondary forests in central Vietnam [23]. These actions have created forest fragments while severely impacting biodiversity and associated natural interaction and processes. Similar to other landscapes in this region, habitats in Bach Ma National Park have been manipulated and transformed, creating metapopulations of species, including the granular spiny frog (Quasipaa verrucospinosa), in the entire park. This species has been listed in the IUCN Red List as a near threatened species due to environmental degradation, loss of forest and stream habitats, global climate change, and overexploitation for consumption [24]. However, little is known about many aspects of the population ecology of Q. verrucospinosa in Vietnam and large-scale studies of occupancy models for this species are nonexistent. Information on site occupancy and microhabitat use in granular spiny frogs in primary and secondary forests of Bach Ma National Park is lacking even though terrestrial habitats have been identified as important to conservation programs [23, 25, 26].

In this study, We estimated site occupancy for Q. verrucospinosa in Bach Ma National Park, central Vietnam to (1) compare occupancy and detection probabilities for two specific habitat types (primary and secondary forests); (2) to obtain an overall estimate of site occupancy for the entire national park; and (3) to determine the number of microhabitat use individuals. I examined the effects of site covariates, including temperature, humidity, and precipitation, on the occupancy and detection of frogs, and tested the hypothesis that differences among habitat types result in different levels of detection. The presence of a primary forest canopy that better regulates forest temperature and soil moisture is crucial in determining anuran survival, reproduction, movement models, and species richness [4, 27, 28, 29], and amphibian species may be particularly sensitive to the effects of habitat fragmentation [30, 31, 32]. We predicted a greater abundance of granular spiny frogs in primary forests (undisturbed habitat) than in secondary forests (disturbed habitat).

2. Materials and methods

2.1 Study sites

The fieldwork took place in Bach Ma National Park (15^o59’12″-16^o16’09”N, 107^o37’06″-107^o54’14″E, approximately 37,487 ha in size), central Vietnam (Figure 1). The study area is dominated by montane rain forests at elevations of 700–1400 m a.s.l. and cloud forests from about 1400 m up to summits at 1712 m [23, 33]. Seasonal monsoons and a tropical climate characterize this study area, with annual temperatures averaging 22.6 ± 0.26°C (ranging from 16.9 ± 0.39°C in January to 26.6 ± 0.36°C in June) and an annual mean precipitation of 3492.1 ± 228.7 mm. Most of the rain is concentrated in the main rainy season from September to December (monthly mean: 629.2 ± 44.1 mm) and the little rainy season (monthly mean: 149.2 ± 14.5 mm) from May to August, whereas the period from January to April is relatively dry, with monthly mean rainfall of 94.6 ± 12.8 mm.

Figure 1.

Location of survey sites in Map of Bach Ma National Park, Thua Thien-Hue Province, Central Vietnam, showing 35 sites in primary forests (●) and 42 sites in secondary forests (○), where granular spiny frogs were monitored during the breeding season of 2013.

2.2 Field sampling

The study area comprised primary forests (ca. 32.2%; canopy is not fragmented), secondary and restored forests (54.0%; fragmented canopy), and administrative areas (13.8%; plantations; [23, 34]). From September to December 2013, we conducted seven surveys for Q. verrucospinosa during the peak breeding season [35], but only in primary and secondary forests (i.e., distributed regions of this species). We set 77 sampling plots of 20 × 50 m (1000 m² each site), 35 in primary forests and 42 in secondary forests (Figure 1). We only selected sample sites that contained water bodies, either a part of a stream or marsh where Q. verrucospinosa are commonly active, and each site was located about 300 m apart to ensure independence among sites. No addition, removal, or alteration of plots was made during the entire study period. We considered the primary and secondary forest variables as site covariates to describe habitat occupancy.

Each site was visited every two weeks, and each site was visited and sampled at the same time. In the night, a team of two people walked slowly with a roughly equal pace along the plot, and visually searched for frogs using spotlights from 19:00 to 02:00 hours for 50 m. We searched for Q. verrucospinosa in the water where they were visible and reachable, on land up to 10 m away from the stream or the marsh, and on tree trunks and vegetation. After locating frogs, we collected them by hand [36]. We also adopted the auditory survey method [37, 38] of using calls to detect granular spiny frogs and to count hidden individuals at each site.

2.3 Data analyses and model selection

To determine site occupancy, it is necessary to simply record whether an individual is detected “1” or not detected “0” during each survey at each site when visited. We estimated the effect of the secondary forest variable (with strong disturbance) on occupancy and detection probability. Using field observations on forest canopy gathered prior to this study, we determined the level of the secondary forest variable at each site, and a covariate of secondary forest was defined as 1 if the site showed evidence of the fragmentation of canopy and 0 otherwise. At each time a site was visited, I also recorded air temperature (temp), relative humidity (humi), and precipitation (rain). These variables were considered sample covariates to estimate detection and presence probabilities, respectively. When an individual frog was detected, we recorded the habitat in which it was found as aquatic, terrestrial, or arboreal. We used these variables to estimate microhabitat use in Q. verrucospinosa.

Each site has its own detection history that can be represented by a mathematical equation (Appendix 1). For example, supposing 30 sites were each sampled four times within a season and the target species (Q. verrucospinosa) was detected at site 1 during the first and last survey occasion (1001). The site was occupied (ψ), the probability of detecting the species during the j^th survey was p_j, and the species was detected on the first and last surveys (p₁ and p₄) but not on the second and third surveys. We can write the probability of this detection history as follows: Pr(H_i = 1001) = ψp₁(1 – p₂)(1 – p₃)p₄.

Sites 2 and 30 represent a case where the target species was never detected (detection history = 0000). These sites could either be unoccupied, which mathematically is (1 – ψ), or they could be occupied but not detected. In this case, we can write the probability of this detection history as follows: ψ(1 – p₁)(1 – p₂)(1 – p₃)(1 – p₄) or ψ(1 – p_j)⁴. Thus, we can write the probability of detection history (0000) as follows: ψ(1 – p_j)⁴ + (1 – ψ).

If a site is not surveyed at the j^th survey occasion (θ, the probability of miss detecting the species or missing observations), no information regarding the detection (or non-detection) of this species was collected from that site at that time. For example, the probability of detection histories (sites 3 and 4) could be expressed as: Pr(H₃ = 10x0) = ψp₁(1 – p₂) × θ × (1 – p₄) and Pr(H₄ = x0x1) = ψ × θ × (1 – p₂) × θ × p₄.

Finally, the mathematical equation of all detection histories are combined into model likelihood as follows: L(ψ, p/H₁,…, H₃₀) = Π(i = 1, 2, 3,…30)Pr(H_i). Maximum likelihood methods are incorporated in the program PRESENCE and this software was used to obtain estimates of occupancy and detectability for frogs at Bach Ma National Park.

We used the program PRESENCE and single-season occupancy models to assess occupancy and detection probabilities. This pattern assumes that sites or patches were closed to changes in occupancy between the first and last surveys of a given sampling season (i.e., no colonization or extinction events within the sampling season), and detection of the target species at a site is independent of detecting the species at other sites [3, 21]. We used the following parameters of interest for the present study: ψ is the probability of a site occupied by Q. verrucospinosa and p_j is the probability of detecting the species during the j^th survey given that it is present.

We used two essential models for the present study. The first model that assumes that occupancy and detection probabilities with respect to Q. verrucospinosa are constant across sites and surveys [denoted ѱ(.)p(.)]. The second model assumes constant occupancy among sites, but detection probabilities are allowed to vary among seven survey occasions [denoted ѱ(.)p(survey)]. Previous studies suggest that the detection probability ≥0.15 in the model ѱ(.)p(.) is needed for unbiased occupancy modeling [2, 39]. Our model ѱ(.)p(.) with a detection probability of 0.329 (SE = 0.035) is appropriate to pursue inference and to estimate the details of the parsimonious process of model selection. In the present study, detectability was either constant across all survey occasions and sites p(.), or varied in three possible ways: among seven survey occasions p(survey), or across sites according to weather conditions p(temp, humi, rain), or both p(survey, temp, humi, rain). We also estimated the value of the coefficient for the secondary forest variable with respect to its influence on occupancy probability and our candidate set contains 16 models without considering interactions between factors (Table 1).

Model	Np	AICc	∆AICc	w	–2 l	SF1	±SE
ѱ(PF), p(sur, temp, humi, rain)	11	322.33	0.00	0.242	300.33	—	—
ѱ(PF), p(temp, humi, rain)	4	322.45	0.12	0.228	314.45	—	—
ѱ(PF), p(temp, humi, rain, SF)	5	323.39	1.06	0.143	313.39	—	—
ѱ(SF), p(temp, humi, rain)	5	324.05	1.72	0.102	314.05	−0.549	0.727
ѱ(SF), p(sur, temp, humi, rain)	12	324.12	1.79	0.099	300.12	−0.401	0.781
ѱ(PF), p(sur, temp, humi, rain, SF)	12	324.19	1.86	0.095	300.19	—	—
ѱ(SF), p(temp, humi, rain, SF)	6	325.36	3.03	0.053	313.36	−0.221	1.141
ѱ(SF), p(sur, temp, humi, rain, SF)	13	326.06	3.73	0.038	300.06	−0.341	0.910
ѱ(SF), p(SF)	4	352.95	30.62	0.000	344.95	−0.219	1.562
ѱ(PF), p(SF)	3	354.83	32.50	0.000	348.83	—	—
ѱ(PF), p(sur, SF)	9	360.85	38.52	0.000	342.85	—	—
ѱ(SF), p(sur, SF)	10	362.83	40.50	0.000	342.83	−0.225	1.555
ѱ(SF), p(.)	3	371.27	48.94	0.000	365.27	−1.929	0.555
ѱ(SF), p(sur)	9	381.26	58.93	0.000	363.26	−1.928	0.554
ѱ(.), p(.)	2	383.14	60.81	0.000	379.14	—	—
ѱ(.), p(sur)	8	393.13	70.80	0.000	377.13	—	—

Table 1.

The summary of AIC model selection from 77 sites in the primary and secondary forests of Bach ma National Park. ∆AIC_c is the difference in AIC value for a particular model when compared with the top-ranked model; w: The AIC model weight; Np is the number of parameters; −2 l is twice the negative log-likelihood value; SF₁: The coefficient for the secondary forest variable with respect to its effect on occupancy probability; SE: Standard error; PF: Primary forest; SF: Secondary forest; temp: Temperature; humi: Humidity; rain: Rainfall; Sur: Survey.

We used the Akaike Information Criteria for small sample size (AIC_c), the differences in the Akaike Information Criteria for a particular model when compared to the top-ranked model (∆AIC_c), the AIC model weight (w), the number of parameters for each model (Np), and twice the negative log-likelihood value (−2 l), to establish the process of model selection [40]. All models with AIC differences of <2.0 have a substantial level of empirical support and should be considered when making statistical inferences or reporting parameter estimates of the best models, and the AIC weights summed to 1.0 for all of the members of the model set [40].

We followed a two-step process to detect site occupancy patterns [4, 6, 41, 42]. We first examined occupancy patterns as a function of site covariates employing the best detection patterns (as indicated by AIC_c weight). Burnham and Anderson [40] recommend using AIC_c for model ranking to help account for small sample sizes. Then, we modeled detection probability as a function of the sample covariates while keeping occupancy constant [i.e., ѱ(.)p(covariates)]. Finally, we examined the null hypothesis, the model ѱ(.)p(.), and one alternative hypothesis, the model ѱ(.)p(survey), using the χ²-test by the linear interpolation [6, 43] prior to inference for the next models.

We employed the equation, Ψ_mle = S_D/SP_mle, to assess a constant detection probability. S_D is the number of sites where Q. verrucospinosa was detected at least once, S is the total number of sites, and P_mle is the estimated probability of detecting the species at least once during a survey given it is present. P_mle = 1 – (1 – p_mle)⁷, where p_mle is the maximum likelihood estimate of a constant detection probability in a single survey of an occupied site, and mle represents the maximum likelihood estimate of the respective model parameters. Thus, the probability of detecting at least once Q. verrucospinosa after k surveys (k = 7) of the site will be p* = 1 – (1 – p)⁷, where p assumes that Q. verrucospinosa is detected imperfectly and gives the probability of detecting Q. verrucospinosa in a single survey of an occupied site. This is one minus the probability of Q. verrucospinosa being undetected in all k surveys [6]. For the probability of occupancy given that a species is not detected at a site (i.e., that Q. verrucospinosa was present at a site given it was never detected), we used the following equation: Ψ_condl = ѱ(1 – p_j)⁷]/[(1 – ѱ) + ѱ(1 – p_j)⁷], where ѱ is estimated occupancy probability, p_j is the detection probability estimates in the j^th survey. We obtained the standard errors for Ψ_mle and Ψ_condl by application of the delta method, where the variance–covariance matrix for ѱ and p given by the program PRESENCE.

To estimate the fit of the model to data, accounting for an over-dispersion in model selection, and to account for missing observations, we used a simple Pearson’s Chi-square statistic to test whether there was sufficient evidence of poor model fit [5]. We calculated the observed χ² goodness-of-fit statistic and estimated an over-dispersion parameter (ĉ) from 10,000 bootstrap samples using the program PRESENCE for the global model (the most general model in the model set or the model with the most parameters) from the candidate models (i.e., AIC weight > 0.001). Over-dispersion is common in ecological models, and adjusting the model selection criteria is recommended [5, 6]. In the absence of a global model, the weighted-average (w_i) of ĉ was used to portray a goodness-of-fit for all candidate models [40]. If the values of weighted ĉ are greater than one, it suggests that there is more variation in the observed data than expected by the model (over-dispersed occupancy models), while values less than one suggest less variation. If the values of weighted ĉ are equivalent to one, then the target model is an adequate description of the data [5, 6].

We analyzed the data using the program PRESENCE (USGS-Patuxent Wildlife Research Center, Maryland, USA). We used SPSS v.16.0 (SPSS Inc., Chicago, Illinois, USA) for Windows 10 to analyze the data of microhabitat use, and set the significance level at P ≤ 0.05 for all analyses. To test the number of individuals using microhabitats and among surveys, we used a one-way analysis of variance (ANOVA). We used the χ² tests to examine the significance level between the first model ѱ(.)p(.) and the second model ѱ(.)p(survey) through the seven surveys. We tested the possible effects of climatic factors (air temperature, relative humidity, and precipitation) on the detection of individuals using multiple regression analyses. All data are presented as mean ± 1 SE (unless otherwise noted).

3. Results

Quasipaa verrucospinosa was detected at least once at 31 of the 77 sites, yielding an overall naïve occupancy estimate of 0.403, clearly indicating that detection probabilities are less than one. There conceivably can be a number of locations where granular spiny frogs were present but simply never detected during the seven survey occasions. Our detection-corrected occupancy estimates by site in the primary and secondary forests of the national park ranged 0.143–0.714 (average naïve occupancy = 0.351 ± 0.032). As a general approach, two essential models [ψ(.)p(.) and ψ(.)p(survey)] need to consider before inferring next models. The first model assumes that occupancy and detection probabilities are constant across sites and surveys. The rate of sites occupied by Q. verrucospinosa from the constant model [ψ(.)p(.)] was 0.433 (SE = 0.061). The second model assumes constant occupancy among sites, but detection probabilities are allowed to vary among the seven surveys. The rate of sites occupied based on the second model of [ψ(.)p(survey)] was 0.432 (SE = 0.061).

The estimated occupancy probability is very similar in both models, 0.433 and 0.432 from the first and secondary models, respectively. When estimating occupancy probabilities including only the two models [(ψ(.)p(.) and ψ(.)p(survey)], both models gave essentially the same results, and both are about 8% larger than the naïve occupancy estimate. The model-averaged estimate of occupancy probability between primary and secondary forest habitat categories was 0.433 (SE = 0.061). When examining the results in which parameter estimates only have the two models [ѱ(.)p(.) and ѱ(.)p(survey)], a difference of 9.99 ∆AIC_c units between these two models [with the AIC weight value of 0.993 in the model ѱ(.)p(.)] shows that the model ѱ(.)p(.) is the “best” model. However, the second model [ѱ(.)p(survey)] still has a reasonable relative level of support (the AIC model weight value of 0.007) and there is further evidence of this second model to pursue inference. We examined a likelihood proportion of the null hypothesis of detection probability being constant and the alternative hypothesis that detection probability differs among the seven survey occasions. The test statistic for this is 379.14–377.13 = 2.01 (Table 1), compared to the χ² distribution with 8–2 = 6 degrees of freedom, by the linear interpolation, resulting in a significant level of P = 0.933. Thus, there is insufficient evidence to reject the null hypothesis in this study.

Testing the global model (the model with the most parameters) from the candidate set (Table 1), the model ѱ(secondary forest)p(survey, temperature, humidity, precipitation, secondary forest), does not show any evidence of over-dispersion (weighted ĉ = 0.436), indicating insufficient evidence of the poor model fit using 10,000 bootstrap iterations. As a result, the adjustment has been made to the model selection procedure (AIC) and parameter assessments to estimate the details of this parsimonious process of model selection. Detectability varied among surveys and possibly among sites with previously disturbed and undisturbed histories (Figure 2).

Figure 2.

Estimating the average pattern of detectability across surveys and among sites with different disturbance histories of granular spiny frogs. Undisturbed habitat (−--○---) and disturbed habitat (─●─).

Our candidate set contained 16 models without considering interactions between factors due to limitation of software and complexity of models (Table 1). There was no single model that was demonstrably better than the others. As a general rule, the six top models are separated by less than 2.0 AIC units, which means that these models have substantial support and should be considered when reporting parameter estimates or making inferences (Table 1). The AIC model weight (w) was distributed across a number of models, indicating that a number of models may be reasonable for our collected data. In terms of model weights, the p(temperature, humidity, precipitation) models have 90.9% of the total, providing clear evidence that weather condition is an important factor in terms of accurately modeling detection probabilities. In terms of comparing hypotheses, the hypothesis that the detection probability varied among weather conditions, therefore, has much greater support than the hypothesis that it was constant. Many of the top-ranked models also contained the factor “survey” for detection probabilities, providing evidence that the survey occasions differed in their ability to find Q. verrucospinosa in the sites; a combined model weight for p(survey) models is 43.6% of the total. There was substantially less support for the hypothesis that the level of the secondary forest variable affected detection probabilities for Q. verrucospinosa, with a combined model weight of 23.8% (Table 1).

The primary forest variable ranked first among the set of models that accounted for differences in the survey, temperature, humidity, and precipitation to explain occupancy, and detection probabilities were approximately 2.5 times more likely than the next best model (evidence ratio [Akaike weight of top model/Aikaike weight of second best model] = 2.45). A model including temperature, humidity, and precipitation from primary forest sites ranked secondly among the set of models to explain the probability of occupancy and detectability were about 2.3 times more likely than the next competing model from secondary forest sites (evidence ratio [0.228/0.102] = 2.25). Estimating detection probabilities for each sampling covariate on each survey occasion is given in Table 2.

Sampling covariates	Survey occasions
	1	2	3	4	5	6	7
Temperature	0.283	0.259	0.306	0.241	0.176	0.257	0.254
Humidity	0.274	0.249	0.297	0.228	0.161	0.247	0.243
Precipitation	0.341	0.347	0.371	0.299	0.242	0.292	0.299

Table 2.

Estimating detection probabilities for each survey in Quasipaa verrucospinosa from the AIC occupancy model selection for sampling-specific covariates.

In terms of occupancy probability, based upon rankings and AIC model weights, the results are somewhat conclusive about the effect of secondary forest sites (29.2%) on the ѱ(primary forest) model. The combined weight for the ѱ(primary forest) models was 70.8%, and the ѱ(secondary forest) models was 20.1% (Table 1). The coefficient value for the secondary forest variable with respect to its effect on occupancy probability, the eight AIC selection models showing the negative SF values (all values of â₂ < 0), indicating certain evidence that the probability of occupancy is higher at the primary forest sites than at the secondary forest sites (Table 1). From the top-ranked model with ∆AIC_c < 2.0 units, the model ѱ(secondary forest)p(survey, temperature, humidity, precipitation) on the logit scale produces the following equation for estimating occupancy: Logit (ѱ_i) = 1 × â₁ + â₂ × SF_i = 1 × 0.608 + (−0.401) × SF_i.

For a primary forest site (where the secondary forest variable = 0, according to the value of SF_i = 0), which gives the odds of occupancy of e^0.608 = 1.837 (:1) and a probability of occupancy of 1.837/(1 + 1.837) = 0.647. The odds ratio for a secondary forest site being occupied (value â₂ = −0.401) by Q. verrucospinosa is e^–0.401 = 0.669. Thus, the odds of occupancy at a secondary forest site is 0.669 × 1.837 = 1.231 (:1) or a probability of occupancy of 1.231/(1 + 1.231) = 0.552. I also estimated a confidence interval for the influence of the secondary forest variable on site occupancy based upon the logit scale, an approximate two-sided 95% confidence interval is −0.401 ± 2 × 0.781 = (−1.963, 1.161), giving an interval of (e^–1.963, e^1.161) = (0.140, 3.193) for the odds ratio.

In terms of the overall estimate of site occupancy based upon the top-ranked model ѱ(secondary forest)p(survey, temperature, humidity, precipitation), an average from the estimated occupancy probabilities for the primary forest sites (35 sites) and the secondary forest sites (42 sites), an overall estimate based on the influence of the secondary forest variable was {(35 × 0.647 + 42 × 0.552)/(35 + 42)} = 0.595, with an SE value of 0.114. This is approximately 48% larger than the naïve occupancy estimate (the fraction of sites where Q. verrucospinosa was detected) of 0.403. However, this is about 9% smaller than the occupancy estimate in the “best” model ѱ(primary forest)p(survey, temperature, humidity, precipitation) of 0.632 (SE = 0.078). Clearly, accounting for detection probability has increased the estimated level of occupancy as expected (we discuss in detail below why the overall level of occupancy is larger than the naïve estimate). Based upon Bayes’ Theorem, we also estimated the probability of a site being occupied, given the granular spiny frog Q. verrucospinosa was not detected there in any of the seven survey occasions (Ψ_condl), from the “best” model [ѱ(.)p(.)], Ψ_condl = 0.433 × (1–0.329)⁷]/[1–0.433 × {1 – (1–0.329)⁷}] = 0.044, with an estimated SE value of 0.021. The value of p* in this model where the probability of detection is constant, the probability of detecting Q. verrucospinosa at least once after k surveys of the site was p* = {1 – (1–0.329)⁷} = 0.939.

An overall estimate of microhabitat use in Q. verrucospinosa showed that the number of granular spiny frogs using the terrestrial habitat (121 individuals, 56.0%) was larger than the aquatic habitat (82 individuals, 38.0%) or the arboreal habitat (13 individuals, 6.0%). The number of individuals was significantly different among three habitat types (F_2,20 = 101.58, P < 0.001; Figure 3). In total, we found 216 individuals during the seven surveys. The number of individuals was found among seven survey occasions were not significantly different (F_6,75 = 0.94, P = 0.472). Multiple regression results for possible effects of air temperature, relative humidity, and precipitation on the detection of individuals were significant among surveys (R² = 0.139, F_3,251 = 27.92, P < 0.001).

Figure 3.

Microhabitat use in granular spiny frogs from Bach Ma National Park, Central Vietnam. Aquatic (░), terrestrial (□), and arboreal (■) habitats.

4. Discussion

Our results indicate that Q. verrucospinosa occupancy in the tropical forests of Bach Ma National Park is not associated with the secondary forest variable. Excluding the two last models, the six models include the secondary forest covariates (including both occupancy and detection probability) with the values of ∆AIC_c > 30.6 units and all AIC model weights are equal zero. These findings were similar to those of previous studies that Q. verrucospinosa frogs were mainly found in primary forests in central Vietnam [24, 33, 34, 35, 44, 45]. In fact, we sampled a relatively broad gradient of forest types (42 sites were classified as secondary forest and 35 sites were classified as primary forest). However, Q. verrucospinosa frogs were only found at eight sites with eight respective individuals in a total of 42 sampling sites in secondary forests. Thus, we speculate that air temperature, relative humidity, and abundant precipitation during our sample season may have lessened forest type effects, because weather conditions and survey factors were important covariates for occupancy and detectability of Q. verrucospinosa in Bach Ma National Park tropical forests. The presence of a forest canopy that regulates air temperature and forest soil moisture appears more critical in determining survival of amphibians and movement models (e.g., [46, 47]). Previous studies show that some anuran species (especially juveniles) have indicated a preference for habitat types with forested canopies compared with fragmented forests or open-vegetation types [48, 49].

Wildlife occupancy relates only to site characteristics, whereas the probability of detecting a species during a single survey can vary with survey characteristics (e.g., temperature and precipitation) or site characteristics (e.g., habitat variables such as primary and secondary forests; [6]). An observed absence at a site occurs if either the species was truly absent, or the species was present at that site but not detected simply; while non-detection of a species does not mean that that species was truly absent unless the probability of detecting the species was 100%. That is the reason why previous occupancy studies of wildlife populations are often impeded by imperfect detectability [16, 21]. Thus, the rate of sites where a species of interest is detected will always be an underestimate with respect to the true occupancy level in the study region when detection is imperfect. Hence, inferences regarding the effects of site characteristics on habitat occupancy will be difficult or impossible to describe exactly [6, 16]. Our results from the best model ψ(primary forest)p(survey, temperature, humidity, rainfall) in the total candidate models were reliable, with the occupancy estimate of 0.632 (CI = 0.471–0.768) compared to the naïve occupancy estimate of 0.403. Although we did not consider colonization probabilities and local extinction factors, these two variables often influence parameter estimates in long-term monitoring programs of amphibians [21, 50, 51].

Moreover, our parameter estimates have satisfied normal assumptions of a model of single-species and single-season occupancy [6], including (1) the occupancy state of the sites does not change during the survey period, but can change between survey periods, (2) the detection of the target species in each survey occasion of a site is independent of detections within other survey occasions of the site, (3) occupancy probability is the same across sites or differences in habitat occupancy may be explained with site traits (covariates), (4) species detection probability at occupied sites is the same across all sites and surveys, or differences in detection probability can be explained with survey or site traits, and (5) the detection histories observed at each location are independent for a species of interest.

A brief examination of the estimated detection probabilities clearly indicates why the overall level of occupancy is estimated to be 49% larger than the naïve occupancy estimate, the estimation based simply on the number of sites where Q. verrucospinosa frogs were detected during the seven surveys. There is clearly a reasonable level of survey variation and substantial differences among sampling covariates, including temperature, relative humidity, and precipitation (see results). Furthermore, the detection of a species at each site is indeed indicative of the presence but non-detection of a species is not equivalent to the absence, unless the detecting probability of the species was one, and a species can go undetected at a site or some sites even when present. Therefore, non-detection sites of a frog represent a case where the target species (Q. verrucospinosa) was never detected. These sites could either be unoccupied, which mathematically is (1 – ψ), or they could be occupied but we never detected the target species during k survey occasions, which mathematically is ψ(1 – p_j)^k. Both of these detectabilities have been included in maximum likelihood methods incorporated in the program PRESENCE to obtain estimates of occupancy and detectability. Although estimates in both essential models [(ψ(.)p(.) and ψ(.)p(survey)] are about 8% larger than the naïve occupancy estimate (suggesting that Q. verrucospinosa was never detected at one in every seven surveys), we believe that Q. verrucospinosa frogs can be more likely to occupy primary forest locations compared to secondary forest locations.

Parameter assessments and associated confidence intervals from pattern averaging indicated that the primary forest was an important determinant of Q. verrucospinosa occupancy in Bach Ma National Park. The appearance of both covariates (primary and secondary forests) in competing patterns is not a surprising result, but the weak relationship to habitat occupancy was unexpected. In terms of occupancy and detection probabilities, the negative values of the secondary forest variable indicated certain evidence with respect to its effect on occupancy and detectability of Q. verrucospinosa. The effects of forest and year factors on occupancy and detectability of tropical amphibians and habitat use emphasize the importance of conducting longer term researches for describing critical habitat relationships [52, 53]. Some populations of salamanders and frogs can fluctuate in the number of individuals (or even in number of species, genera, or orders) among breeding seasons [54, 55, 56], and with temperature and rainfall varying annually among forest categories, forest and year effects on occupancy are often common [4, 57].

Our results indicate that precipitation, temperature, and relative moisture were the most important sampling covariates for detection probabilities of Q. verrucospinosa. The importance of these environmental factors on amphibian breeding activity [58, 59], capture proportions [60], and calling of anuran species [61], and hence detection probability, has been well documented. In many cases, precipitation and temperature are expected to be good predictors of detectability for amphibian species. Although we did not analyze interactions between temperature, moisture, and precipitation on detectability in the present study, these variables often interact to affect amphibian timing [56, 62] and movement physiology [63]. A recent study indicated that detectability of anuran species is independently positively associated with temperature and precipitation, with temperature consistently having a greater effect [4]. The survey variable is also an important relative covariate for detection probabilities in this study, and detectability varied within the seven survey occasions and possibly among sites with different disturbance histories.

Missing observations are a common case in the model of the presence or absence of a species from a collection of sampling sites, and occur widely applied in wildlife and ecological studies [3, 4]. Missing observations may arise through a number of reasons, such as a vehicle or equipment breakdown or logistic difficulties in getting field personnel to all sites. Therefore, it may not be possible to survey all sites at all sampling occasions. These sampling inconsistencies can be accommodated using the proposed model likelihood. If a site is not surveyed at the j^th survey occasion, no information regarding the detection (or non-detection) of this species has been collected from that site at that time. Our observed data (77 sites with the seven surveys were conducted), including 17 missing observations, the percentage of missing data for the probability of Q. verrucospinosa presence was 3.15%. According to MacKenzie et al. (2002) on average, the standard error of ψ increased about 5% with 10% missing observations, and about 11% with 20% missing observations. Thus, our occupancy estimates and the bootstrap standard error estimates in the present study were reliable and accounted well for the loss of information for Q. verrucospinosa.

A common method of estimating over-dispersion is to use the observed chi-squared goodness of fit statistic for a global model (the most complex model with the greatest number of parameters), which should be estimated for lack of fit first [6]. According to previous studies estimating the occupancy and detection probability, it should be demonstrated that a fitted model adequately describes the observed data [64, 65]. Substantial lack of fit in the model may lead to inaccurate inferences, either in terms of bias or in terms of precision (e.g., reported standard errors are too small; [5]). Our global model does not indicate any evidence of lack of fit using 10,000 bootstrap samples, with an estimated over-dispersion parameter ĉ of 0.436. However, comparing our results with those assessing the fit of site-occupancy models given in MacKenzie and Bailey [2] and MacKenzie et al. (2006) showed that our models are suitable and that there is insufficient evidence of a poor model fit.

Amphibian species must choose terrestrial microhabitats that prevent loss of excessive water corresponding to each season and, thus, maintain hydration [66]. Moist environments are very important to the survival of juveniles [56, 67]. In this study, we detected frogs using all three habitat types during seven surveys in the main rainy season. In there, using terrestrial, aquatic, and arboreal habitats were 56%, 38%, and 6%, respectively. Although our sampling sites tended to remain moist throughout the sample period in the primary and secondary forests of Bach Ma National Park, there is evidence that in the summer (April to July), the rate of microhabitat use in the semi-aquatic model (about 60%) is larger than the terrestrial model (about 35%, B.V. Ngo, unpubl. data). This can explain why temperature, precipitation, and relative humidity were associated with detectability of frogs. We speculate that an overall moist forest environment (temperature and rainfall) coupled with species-specific behavioral adaptations (e.g., [56, 67]) allowed frogs to remain equally active across the range of precipitation events in our sampling period.

5. Conclusions

Based on the detection/non-detection data for each site over multiple visits for granular spiny frogs, from the best model among all candidate models, we estimated a site occupancy rate of 0.632 that was higher than the naïve occupancy estimate of 0.403 and a 57% increase over the rate of sites at which frogs were actually observed. The site variable of primary forest was an important determinant of site occupancy, whereas occupancy was not associated with the variable of secondary forest. The species detection model p(temperature, humidity, rainfall) included 90.9% of the total weight, providing clear evidence that environmental conditions were important sample covariates in modeling detection probabilities.

Acknowledgments

We thank the staff of A Luoi and A Pat Forest Ranger Stations, Bach Ma National Park, Border Stations 629 and 633, Phong Dien Nature Reserve, and Sao La Conservation Area, for support of our feldwork. We thank colleagues from Department of Life Sciences, National Cheng Kung University, Taiwan, for their support. The Ministry of Education and Training funded this study under grant number B2020-DHH-08. The authors also acknowledge the partial support of Hue University under the Core Research Program, grant number NCM.DHH.2018.10.

Site or Patch	Survey Occasion
	I	II	III	IV
Site 1	1	0	0	1
Site 2	0	0	0	0
Site 3	1	0	θ	0
Site 4	θ	0	θ	1
……	—	—	—	—
Site 30	0	0	0	0