**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–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–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–9]. Large-scale monitoring programs for amphibian species [10–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–15]. Ignoring detectability may lead to biased estimations of site occupancy [16–18] and studies of habitat occupancy are often hampered by imperfect detectability for the species [1, 19–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 largescale 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–29], and amphibian species may be particularly sensitive to the effects of habitat fragmentation [30–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<sup>o</sup> <sup>59</sup>'12″-16o 16'09"N, 107o <sup>37</sup>'06″-107<sup>o</sup> <sup>54</sup>'14″E, approximately 37,487 ha in size), central Vietnam

*Ecology of the Granular Spiny Frog* Quasipaa verrucospinosa*… DOI: http://dx.doi.org/10.5772/intechopen.99656*

#### **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.*

(**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.

#### **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<sup>2</sup> 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 *pj*, and the species was detected on the first and last surveys (*p*<sup>1</sup> and *p*4) but not on the second and third surveys. We can write the probability of this detection history as follows: Pr (Hi = 1001) = ψ*p1*(1 – *p2*)(1 – *p3*)*p4*.

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 – *p1*)(1 – *p2*)(1 – *p3*)(1 – *p4*) or ψ(1 – *pj*) 4 . Thus, we can write the probability of detection history (0000) as follows: ψ(1 – *pj*) <sup>4</sup> + (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(H3 = 10x0) = ψ*p1*(1 – *p2*) θ (1 – *p4*) and Pr(H4 = x0x1) = ψ θ (1 – *p2*) θ *p4*.

Finally, the mathematical equation of all detection histories are combined into model likelihood as follows: *L*(ψ, *p*/H1, … , H30) = Π(*i* = 1, 2, 3, … 30)Pr(Hi).

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 *pj* 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**).

We used the Akaike Information Criteria for small sample size (AICc), the differences in the Akaike Information Criteria for a particular model when compared to the top-ranked model (ΔAICc), 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 AICc weight). Burnham and Anderson [40] recommend using AICc 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 χ<sup>2</sup> -test by the linear interpolation [6, 43] prior to inference for the next models.

We employed the equation, Ψmle = *SD*/*SP*mle, to assess a constant detection probability. *SD* 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) 7 , 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*) 7 , 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


#### **Table 1.**

*The summary of AIC model selection from 77 sites in the primary and secondary forests of Bach ma National Park. ΔAICc 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; SF1: 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.*

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 – *pj* ) 7 ]/[(1 – ѱ) + ѱ(1 – *pj* ) 7 ], where ѱ is estimated occupancy probability, *pj* 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 χ<sup>2</sup> 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 (*wi*) 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 (overdispersed 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 χ<sup>2</sup> 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).
