Regular articles
Special Issues



Caribbean Naturalist
    CANA Home
    Range and Scope
    Board of Editors
    Staff
    Editorial Workflow
    Publication Charges
    Subscriptions

Other EH Journals
    Northeastern Naturalist
    Southeastern Naturalist
    Urban Naturalist
    Eastern Paleontologist
    Eastern Biologist
    Journal of the North Atlantic

EH Natural History Home

Occupancy and Index of Abundance of Eleutherodactylus wightmanae and E. brittoni along Elevational Gradients in West-Central Puerto Rico
Kelen D. Monroe, Jaime A. Collazo, Krishna Pacifici, Brian J. Reich, Alberto R. Puente-Rolón, and Adam J. Terando

Caribbean Naturalist, No. 40

Download full-text pdf (Accessible only to subscribers. To subscribe click here.)

 

Site by Bennett Web & Design Co.
Caribbean Naturalist 1 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 22001177 CARIBBEAN NATURALIST No. 40N:1o–. 1480 Occupancy and Index of Abundance of Eleutherodactylus wightmanae and E. brittoni along Elevational Gradients in West-Central Puerto Rico Kelen D. Monroe1, Jaime A. Collazo2,*, Krishna Pacifici3, Brian J. Reich4, Alberto R. Puente-Rolón5, and Adam J. Terando6 Abstract - Populations of Eleutherodactylus species in Puerto Rico have declined in recent decades due to habitat loss and long-term climatic changes. The conservation of these habitat specialists requires an understanding of factors influencing their abundance and distribution, which at present is scant. We estimated occupancy probability and the probability of encountering ≥2 indivivuals of E. wightmanae (Melodious Coqui or Wightman’s Robber Frog) and E. brittoni (Grass Coqui), species with contrasting habitat affinities, using multi-season, multi-state occupancy models. These parameters also served as an index of abundance (non-presence, 1, and ≥2 individuals). We modedled parameters as a function of seasonal temperature and humidity, long-term average monthly precipitation, and habitat covariates measured at survey sites along 2 elevation gradients in the southern slopes of west-central Puerto Rico. We collected survey data using passive acoustic recorders during 3 seasonal periods between February and July 2015. Occupancy patterns of both species was unimodal, containing higher probabilities (e.g., ≥0.5) at elevations between 400 m and 700 m, where long-term monthly precipitation vaired between 120 mm and 160 mm. Chances of encountering ≥2 individuals increased with ground cover for E. brittoni, and decreased with increasing canopy cover for E. wightmanae. Seasonal temperature and relative humidity did not influence occupancy or the probability of encountering ≥2 individuals, likely because covariates varied within known tolerance levels for Eleutherodactylus. Our findings help reduce local extinction probability through management of habitat conditions that increase the likelihood of encountering ≥2 individuals. We also detailed an analytical framework suitable to test hypotheses aimed at predicting potential impacts from land use and climatic changes, and species responses to conservation actions. Introduction The long-term persistence of amphibians worldwide is threatened by land-use and land-cover changes, and by the projected rise in temperatures of between 2 °C and 4 °C due to anthropogenic climate change (IPCC 2014, Lips et al. 2005, 1North Carolina Cooperative Fish and Wildlife Research Unit, Department of Applied Ecology, North Carolina State University, Raleigh, NC 27695, USA. 2US Geological Survey, North Carolina Cooperative Fish and Wildlife Research Unit, North Carolina State University, Raleigh, NC 27695, USA. 3Department of Forestry and Environmental Resources, Program in Fisheries, Wildlife, and Conservation Biology, North Carolina State University, Raleigh, North Carolina 27695, USA. 4Department of Statistics, North Carolina State University, Raleigh, NC 27695, USA. 5University of Puerto Rico, Department of Biology, Mayaguez, PR 00681, USA. 6US Geological Survey, Southeast Climate Science Center, Raleigh, NC 27695, USA. *Corresponding author - jaime_collazo@ncsu.edu. Manuscript Editor: Robert Powell Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 2 McCallum 2007, Wake and Vredenburg 2008, Watling et al. 2009). These threats, not surprisingly, are also prevalent in and projected for the Caribbean. Puerto Rico is an example where historic land-use changes have been extensive, although forest cover has increased in recent decades (Pares-Ramos et al. 2008, Radeloff et al. 2015), and where climate-model projections suggest that drier and warmer conditions (4.6 °C and 9 °C, depending on emission scenario) will take hold in the future (Henareh Khalyani et al. 2016). Puerto Rico harbors 25 species of amphibians, 17 in the genus Eleutherodactylus. This taxonomic group exhibits high endemism (15 species), and most species are considered habitat and microclimatic specialists (Joglar 1998). Two of these species are federally listed as endangered (E. jasperi Drewry & Jones [Golden Coqui] and E. juanariveroi Ríos-Lopez & Thomas [Plains Coqui]) under the United States Endangered Species Act. The remaining are considered species at risk of endangerment due to habitat loss and threats mediated by climatic changes that suppress successful reproduction and exacerbate diseaseinduced mortality (Burrowes et al. 2004, 2008; Joglar 1998; Pounds and Crump 1994; Pounds et al. 2006). The conservation predicament of Eleutherodactylus spp. underscores the need for a greater understanding of factors influencing their distribution, abundance, and other vital rates (e.g., survival, reproduction). This understanding is essential to categorize habitat quality, and importantly, guide management actions to enhance local populations (Frishkoff et al. 2015, McDonald-Madden et al. 2011, Watling and Donnelly 2009). Available information is valuable, but a substantial amount of it either stems from anecdotal observations and associated habitat, or studies driven by other primary objectives (Drewy and Rand 1983, Joglar 1998, Rivero 1978). With the exception of Barker and Rios-Franceschi (2014) and Burrowes et al. (2004), it is difficult to determine how demographic parameters change as a function of varying habitat or climatic conditions. The former study estimated occupancy probability to assess the status of E. portoricensis Schmidt (Mountain Coqui) and created a suitability map that reflected explanatory geographic and climatic covariates to guide habitat conservation. The latter study reported the relationship between changing micro-habitat and climatic conditions and vulnerability of amphibians, including 8 species of Eleutherodactylus, to several diseases. We also note that these studies were conducted primarily in moist to wet forests (86% of Puerto Rico; Ewel and Whitmore 1973). For example, the lowest sampling location reported by Barker and Rios-Franceschi (2014) was 564 m in elevation. Drewry and Rand (1983) covered a greater range of elevations as they categorized calling sites for Eleutherodactylus, but as indicated above, species–habitat relationships were descriptive. In this study, we expand the aforementioned body of work in Puerto Rico in 2 ways. First, we estimated local occupancy probability of E. wightmanae Schmidt (Melodious Coqui or Wightman’s Robber Frog) and E. brittoni Schmidt (Grass Coqui) over 3 seasonal periods from February to July 2015. In addition, if the site was occupied, we estimated the probability of encountering ≥2 individuals at the survey site. We treat the 3 states or categories reported in this study (i.e., nonpresence, occupied [1 individual], and ≥2 individuals [couple to many]) as an index Caribbean Naturalist 3 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 of abundance. Occupancy is relevant to our understanding of distribution (Yackulic et al. 2015) and abundance is essential to assess the status of species and responses to conservation actions (Nichols and Williams 2006). Second, we obtained these parameter estimates along 2 elevational gradients that spanned dry-lowland to wethighland secondary forests on the southern mountain slopes of west-central Puerto Rico. We used passive acoustic monitoring devices to survey both species (Campos- Cerqueira and Aide 2016). Parameter estimates were obtained using multi-season, multi-state occupancy models, a framework that adjusts for imperfect detection (Nichols et al. 2007). We selected E. wightmanae and E. brittoni because they have contrasting habitat affinities (Drewry and Rand 1983, Joglar 1998). Accordingly, we expected occupancy and the probability of encountering ≥2 individuals to vary as a function of forest cover for E. wightmanae, a forest-dwelling specialist, and ground cover for E. brittoni, a herbaceous-dwelling specialist. In concert with evolutionary adaptations of Eleutherodactylus spp., we expected occupancy and probability of encountering ≥2 individuals of both species to be greater at higher elevations because theses areas are characterized by higher long-term average monthly precipitation, lower temperatures, and higher humidity (Burrowes et al. 2004, Drewry and Rand 1983, Joglar 1998). We also expected that occupancy and probability of encountering ≥2 individuals of both species would decline with warmer temperatures and lower humidity recorded in 2015, particularly for E. brittoni because microclimatic conditions in herbaceous habitats would likely be more variable than in forested habitats. We discuss the implications of this work in the context of improving predictive species distribution models and prospective adaptation strategies in Puerto Rico and elsewhere in the Caribbean. Field-Site Description We established acoustic survey stations within dry, moist, and wet subtropical forest life zones along 2 elevational gradients in the west-central mountains of Puerto Rico (Fig. 1; Ewel and Whitmore 1973). These zones created 3 sampling strata within which we allocated 28 survey stations equally between gradients. We established 2 more stations in the moist life zone (6 vs. 4/gradient) because data were also used to estimate occupancy among selected land-use classes within that life zone (Dowdy 2016). Ultimately, we used data from 26 stations because habitat in 2 stations was altered during the course of the study. Stations were 50 m from major roads to minimize edge effects, and ~500 m apart from each other to reduce spatial correlation. Survey stations were between 18°2'18.68"N, 66°38'34.36"W (low elevation) and 18°10'58.27"N, 66°40'29.91"W (high elevation) for the eastern gradient (townships of Ponce to Adjuntas), and between 18°5'12.54"N, 66°58'55.27"W (low elevation) and 18°8'37.56"N, 66°58'47.87"W (high elevation) for the western gradient (townships of Sabana Grande to Maricao). Elevations within these bounds varied between 77 m and 865 m above sea level. We pooled data to summarize vegetation and environmental variables because they were similar between the eastern and western gradients (ANCOVA: P > 0.05), but provide gradient-specific values in Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 4 Figure 1. Map of Puerto Rico depicting the general location of survey stations along elevational gradients on the southern mountain slopes of west-central Puerto Rico. The top panel depicts temperature (degrees Celsius/10), and the bottom panel the annual precipitation (mm) encompassed by survey stations based on WorldClim (Hijmans et al. 2005). See details provided in the text. Table 1. Average (± SE) and range for covariates measured at 26 stations in secondary forest along elevational gradients in west-central Puerto Rico. Average (± SE) of covariate values are also summarized per gradient. Elevation varied from 77 to 865 m. Temperature and humidity are seasonal averages (n = 3). Period of study was February–March (1), April–May (2), June–July (3), 2015. Ground cover was scaled as: 1–5 (1 = 0–20%, 2 = 21–40%, 3 = 41–60%, 4 = 61–80, 5 = 81–100%). Overall East West Covariate average (SE) Min–max average (SE) average (SE) Ground cover 2.89 (0.24) 1–5 2.90 (0.31) 2.87 (0.42) Canopy cover (%) 76.65 (4.85) 23–98 76.62 (6.10) 81.90 (8.31) Forest cover (%) 85.38 (3.50) 29–100 88.67 (3.52) 80.12 (7.09) T1 (°C) 20.02 (0.22) 17.76–23.18 20.13 (0.34) 19.85 (0.19) T2 (°C) 22.76 (0.28) 20.34–25.95 23.00 (0.42) 22.37 (0.25) T3 (°C) 23.81 (0.28) 21.61–27.39 24.25 (0.38) 23.12 (0.32) H1 (%) 95.59 (0.54) 86.49–99.94 95.52 (0.85) 95.71 (0.39) H2 (%) 90.37 (0.92) 80.52–97.56 89.48 (1.39) 91.80 (0.78) H3 (%) 91.10 (0.89) 78.41–98.55 90.06 (1.30) 92.77 (0.87) Caribbean Naturalist 5 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 Table 1. Average seasonal temperature and humidity varied from 17.16 °C to 27.40 °C and 78.42% to 99.95%, respectively (Table 1). Average temperature decreased with increasing elevation over the duration of the study (F = 23.48; df = 1, 24; P < 0.001). Conversely, average humidity increased with elevation (F = 8.29; df = 1, 24; P < 0.008). On average, temperature increased from February through July by nearly 3 °C, and humidity decreased by about 4%. Long-term average monthly precipitation (1950–2000) per survey station between February and July was 136.52 ± 5.42 mm (n = 156) and varied from 25 mm in February to 320 mm in May (World Clim; Hijmans et al. 2005). Coefficient of variation for precipitation seasonality across stations ranged from 46 to 58 (mode = 48; WorldClim [Hijmans et al. 2005]). Methods We conducted acoustic surveys during 3 primary periods and 3 secondary sampling occasions within each primary period, a scheme that followed Pollock’s (1982) robust design. Primary seasons were February–March, April–May, and June–July, representing periods of contrasting cumulative amounts of precipitation in southcentral Puerto Rico (Calvesbert 1970). Long-term average cumulative rainfall totals during these bi-monthly sampling periods at survey sites was 136.23 ± 7.89 mm for February–March, 368.30 ± 20.92 mm for April–May, and 314.57 ± 17.35 mm for June–July (World Clim; Hijmans et al. 2005). Seasonality was also evident when total monthly precipitation in 2015 was averaged for each primary sampling period as illustrated by the Maricao Hatchery Weather Station (station 665911): 162.31 mm for February–March, 273.02 mm for April–May, and 89.91 mm for June–July (Southeast Regional Climate Center, University of North Carolina at Chapel Hill, NC, USA). We used passive acoustic recording devices (Campos-Cerqueira and Aide 2016), surveying each station 3 times within primary seasons (21 days), with a time span of 4–5 weeks in between to avoid sampling toward the end or beginning of each primary period. Each secondary sampling occasion consisted of recording sessions of 1 minute every hour for 2 consecutive evenings from 1800 h to 0600 h. This schedule and frequency of recording made possible (1) placement, collection, and removal of multiple recorders from multiple locations throughout the study; and (2) determination if at least 1 individual (presence) of a species was represented in 24 recordings per secondary sampling occasion. We collected 68,131 recordings throughout the study. We used 3 types of vegetation cover to represent habitat quality: canopy, ground (horizontal), and contextual forest cover. These covariates provided a means to gauge conditions that promote food availability, shelter, and suitable microclimatic conditions for amphibians (Huang et al. 2014, Joglar 1998, Kearney et al. 2009, Whitfield et al. 2014). At the survey station level, we measured canopy and ground cover. Both covariates were measured within 5 m of the survey station. Canopy cover was measured using a densitometer and expressed as percent cover. We measured ground cover using a density pole 0.3 m wide, segmented into 5 sections 0.5 m tall each (ground level to 2.5 m above ground; Nudds 1977). Measures of ground cover were made from 5 m away from the pole, and represent the averCaribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 6 age of 4 readings (i.e., cardinal directions) using an ordinal scale: 1 = 1–20%, 2 = 21–40%, 3 = 41–60%, 4 = 61–80%, 5 = 8–100%. We also quantified percent forest cover within a 100-m radius of the survey station. This covariate represents the extent over which Eleutherodactylus spp. could potentially move based on translocation experiments of E. coqui Thomas (Common Coquí; Gonser and Woolbright 1995). It is also a measure of the forest cover surrounding the station, henceforth termed contextual cover. We obtained estimates of forest cover by creating habitat polygons within a 100-m radius of the station to distinguish forest cover from alternative cover types using satellite imagery from Google Earth (2015). Cover types were forest, low vegetation (e.g., shrubs), and other cover (impervious surfaces, bare soil, grass). We recorded relative humidity (H) and temperature (T), within-season weather covariates, because these parameters influence calling activity and habitat quality in Eleutherodactylus (Burrowes et al. 2004, Ospina et al. 2013). We obtained data using HOBO data loggers (Pro v2 Temp/RH/weatherproof housing; accuracy = 2.5%) at 18 of 26 sites on every sampling occasion (i.e., sampling covariates). Estimates for the remaining 8 sites were interpolated using regression analyses and data on humidity and temperature collected on 288 occasions and pertinent site covariates (e.g., elevation, east–west location). Although temperature and humidity thresholds (thermal limits) were not available for E. brittoni and E. wightmanae, we modeled the aforementioned weather covariates to gain insights on the functional relationships between covariates and model parameters. We also modeled occupancy by the long-term average monthly precipitation (Precip, mm) from 1950 to 2000 for each station (WorldClim; Hijmans et al. 2005). This covariate represented the climatic regime of surveyed sites. Although elevation and precipitation are correlated (r = 0.56), this association does not preclude the possibility that occupancy responds to one covariate (precipitation) differentially at varying levels of the other (elevation). We employed a multi-season, multi-state occupancy modeling framework (Nichols et al. 2007) to obtain estimates of occupancy (Psi), probability of encountering ≥2 individuals (R), detection (p), and delta (dlta) using the program PRESENCE (Hines 2006). Occupancy is defined as the probability that a site is occupied, given that individuals of Eleutherodactylus are available to be detected. Detection (p) is defined as the probability of detection in a survey station given that the station is occupied. Parameter R is the probability that ≥2 Eleutherodactylus were present, given that the survey station was occupied. Parameter delta is defined as the probability that state 2 (≥2 frogs calling) was true, given that the survey station was occupied. We developed a suite of candidate models to test our a priori hypotheses (Table 2). We first modelled the detection process, determining whether detection was constant (.), varied by season (S), and whether detection probability was influenced by temperature (T) and humidity (H) over time (sampling covariate). The latter environmental covariates may influence availability (i.e., calling activity; Ospina et al. 2013). The model with highest support (AICwt) was used to determine if occupancy (Psi) and the probability of encountering ≥2 individuals (R) were constant (.), varied by the initial season as compared to the other two (init), or was season-specific. We included a model evaluating the influence of the initial Caribbean Naturalist 7 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 Table 2. Hypotheses and a priori predictions about the influence of environmental and habitat covariates on local occupancy (Psi) and probability of encountering ≥2 individuals (R) of Eleutherodactylus wightmanae and E. brittoni during 2015 in secondary forest in the leeward side of mountains in west-central Puerto Rico. See the text for detailed descriptions of covariates. Covariates with asterisks (*) had a strong influence (Beta 95% CIs did not overlap zero) on occupancy or abundance categories. Species Hypothesis Covariate E. wightmanae E. brittoni Eleutherodactylus spp. exploit high elevation, Elevation (m) Occupancy (+), Occupancy (+), mesic environments (Joglar 1998) abundance (+) abundance (+) Occurrence of Eleutherodactylus reflects Long-term average Occupancy (+)* Occupancy (+)* adaptations to long-term climatic conditions monthly precipitation (mm) (Burrowes et al. 2004, Joglar 1998) Occupancy and abundance will be affected by Temperature (T), humidity (H) Occupancy: seasonal T (-), H (-); Occupancy: seasonal T (-), H (-); humidity and temperature in June–July as abundance: seasonal T (-), H (-) abundance: seasonal T (-), H (-) compared to February–March due to changes in microclimatic conditions (Huang et al. 2014). Survey site cover influences shelter, food, and Ground (GC), canopy (CC), Occupancy: GC (-), CC (+), Occupancy: GC (+), CC (-), microclimatic conditions (Huang et al. 2014, forest (FC) FC (+); FC (-); Tews et al 2004, Whitfield et al. 2014). abundance: GC (-), CC (+), abundance: GC (+)*,CC (-), FC (+)* FC (-) Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 8 season because long-term rainfall, which affects humidity and habitat quality, was substantially lower than during the second and third primary sampling seasons. Parameters were modeled by elevation (Elev), seasonal humidity (H) and temperature (T), and habitat. Habitat covariates were ground cover (GC), canopy cover (CC), and forest cover (FC). We also modeled occupancy by long-term average monthly precipitation (Precip), the interaction between elevation and precipitation (Elev*Precip), and a quadratic term for elevation to detect possible non-linear patterns in occupancy (ElevQ). We did not use interpolated humidity and temperature data to model occupancy. Instead, we filled missing station values (n = 5) with the average of available station data, standardizing the covariate to retain the same mean and variance (Cooch and White 2010). All continuous covariates were normalized (i.e., mean = 0, standard deviation = 1) before running analyses. Input data (i.e., encounter histories) for analyses were created manually because it required not only identification of species, but also determination of how many individuals were calling concurrently. Because this was a substantial effort, we used a subset of the data to conduct the analyses. We randomly selected 6 recordings from each secondary sampling period for a total of 1404 recordings. We inspected >100 spectrograms visually and by sound, ascertaining species identification and noting call duration, frequency, and bandwidth using the program ARBIMON (Fig. 2; Campos-Cerqueira and Aide 2016). This “standard” served to inspect all remaining spectrograms in the same fashion to determine species presence, and, if present, whether ≥2 individuals were calling. We explored the possibility of classifying parameter R as ≥3 frogs calling concurrently, but we opted not to do so because discerning those occasions from spectrograms was uncertain, leading to misclassifying data. Data were coded as 0 (non-presence), 1 (single individual), and 2 (≥2 individuals present), our indices of abundance. We used Akaike’s information criterion (AIC) to evaluate the support in the data for models in the candidate sets and the strength of each covariate’s effect on Eleutherodactylus species occupancy (Burnham and Anderson 2002). Models with a ΔAIC ≤ 2 were considered to have substantial support in the data. We report only models with ΔAIC ≤ 10 (Burnham and Anderson 2002), but the complete modelselection tables are listed in the supplemental information (see Supplemental File 1, available online at https://www.eaglehill.us/programs/journals/cana/supplementalfiles/ C165-Collazo-s1). The effect of covariates (i.e., b coefficient) on occupancy was considered to be strong if the 95% confidence interval (CI) did not overlap zero, and weak otherwise. Parameter estimates ± SE are reported with occupancy model results. Careful consideration of model assumptions was important for the interpretation of model results. Multi-season, multi-state occupancy models assume that: (1) sites are “closed” to changes in occupancy within primary sampling periods; (2) there are multiple site visits within secondary periods; (3) that there are no false detections; and (4) detection across sites are independent (MacKenzie et al. 2002). We believe we met assumptions for the following reasons. Each survey site was acoustically sampled 3 times (secondary periods) during each primary season (assumption 2). Caribbean Naturalist 9 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 Surveys within primary periods were completed over a short period of time (21 days) relative to life-history events (e.g., reproduction; assumption 1). Encounter histories were derived from acoustic recording devices using training data (>100 recordings evaluated by an experienced local herpetologist and co-author), which provided a means to detect false positives (assumption 3). Survey stations were separated by >500 m, greater than the longest known distance (100 m) from which frogs returned to their site of capture when translocated (Gonser and Woodbright 1995; assumption 4). The model defines parameter delta (dlta) as the probability of correctly identifying state 2, given it is in state 2 (MacKenzie et al. 2006). We assumed that was the case when call patterns (1 vs. ≥2 individuals) could be visually identified using program ARBIMON’s spectrogram (Fig. 2). We believed that 2 states (occupied and ≥2 individuals), as opposed to 3 states or more (e.g., occupied, Figure 2. Spectrogram created with program Arbimon to distinguish instances when 1 individual (top panel) or ≥2 individuals (bottom panel) of Eleutherodactylus wightmanae were calling. Acoustic data were collected in February–March, April–May, and June–July 2015. Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 10 2, and ≥ 3 individuals), produced more accurate results. However, because there is a chance that the spectrogram may not have registered instances of ≥2 individuals, thus not visible to us, we let PRESENCE estimate delta. Results E. wightmanae We recorded this species in 10 of 26 survey stations throughout the study. Variation in occupancy and levels of abundance were best explained by a model in which occupancy was influenced by the interaction between elevation (elev) and long-term average precipitation (precip), and the probability of encountering ≥2 individuals was influenced by station-level canopy cover (AICwgt = 0.79) (Table 3). Detection was lower during the initial primary seasonal period (0.31 ± 0.08) than in the subsequent 2 seasons (0.87 ± 0.06), but not influenced by humidity or temperature. This model lends support the predictive influence of precipitation and elevation on occupancy, but not on the probability of encountering ≥2 individuals, which was contrary to our prediction (Table 2). The interaction between elevation and long-term average monthly precipitation was negative and strong (b = -5.05 ± 2.50), yielding a unimodal occupancy pattern (Fig. 3). Higher occupancy rates (e.g., ≥ 0.5) were recorded at elevations between 400 m and 700 m and were associated with average monthly precipitation of 110 mm and 165 mm (Fig. 4). The probability of encountering ≥2 individuals was strongly and negatively influenced by percent canopy cover (b = -1.81 ± 0.78), and positively, but weakly, influenced by percent contextual forest cover (b = 2.08 ± 1.24). The influence of canopy cover on the probability of detecting ≥2 individuals was illustrated by predicting probability values (R) after removing the effect of forest cover (Fig. 5). The probability of correctly identifying a survey station as state 2 (≥2 individuals) at any survey station in the study was 0.96 ± 0.03. E. brittoni We recorded this species in 14 of 26 survey stations throughout the study. Variation in occupancy and levels of abundance were best explained by a model in which occupancy was influenced by a contrast between occupancy on the first season (init), which was drier when compared to the subsequent 2 seasons, and the interaction of elevation and long-term average monthly precipitation, and the probability of encountering ≥2 individuals was influenced by ground cover and contextual forest cover (AICwgt = 0.85) (Table 3). Detection varied as a function of humidity, with initial season rates varying from 0.25 to 0.28, and subsequent samplingseason rates varying from 0.66 to 0.70. This model lends support to the predictive influence of precipitation and elevation on occupancy, but not on the probability of encountering ≥2 individuals. The predictive influence of ground cover (+) and forest cover (-) matched predictions (Table 2). The influence of the interaction between elevation and average precipitation on occupancy was strong and negative (b = -7.24 ±2.97). Occupancy rates were higher (e.g., ≥ 0.5) at elevations between Caribbean Naturalist 11 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 Table 3. Model selection table for multi-season, multi-state occupancy models for Eleutherodactylus wightmanae and E. brittoni surveyed on southern mountain slopes of west-central Puerto Rico. Acoustic surveys were conducted in January–March, April–May, and June–July 2015. Model parameters were local occupancy (Psi), probability of encountering ≥2 individuals (R), probability of correctly classifying a survey site as harboring ≥2 individuals (dlta), and detection probability (p). Parameters were modeled as constant over time (.), season-specific (S), and initial season vs. remaining two (init). Covariates were temperature (T), humidity (H), long-term average monthly precipitation (Precip), elevation (Elev), and seasonal long-term average precipitation (Season+Precip or Season*Precip), and 3 habitat covariates (ground cover (GC), canopy cover (CC), and forest cover (FC). The full list of competing models is provided as supplemental information. Number of Species/ model AIC ΔAIC AIC wgt parameters -2*LogLike E. wightmanae Psi(Elev+Precip+E*P),R(CC+FC),dlta(.),p(init) 196.91 0.00 0.79 10 176.91 Psi(Elev),R(CC+FC),dlta(.),p(init) 201.48 4.57 0.08 8 185.48 Psi(Elev+ElevQ),R(CC+FC),dlta(.),p(init) 202.71 5.80 0.04 9 184.71 Psi(Precip),R(CC+FC),dlta(.),p(init) 203.53 6.62 0.03 8 187.53 Psi(Elev),R(CC),dlta(.),p(init) 205.59 8.68 0.01 7 191.59 Psi(Elev),R(FC),dlta(.),p(init) 206.67 9.76 0.01 7 192.67 Psi(Elev),R(GC),dlta(.),p(init) 206.86 9.95 0.005 7 192.86 Psi(Elev),R(CC+GC),dlta(.),p(init) 207.36 10.45 0.004 8 191.36 E. brittoni Psi(initSeason+Elev+Precip+E*P),R(FC+GC),dlta(H),p(init+H) 229.07 0.00 0.85 13 203.07 Psi(initSeason+Elev),R(FC+GC),dlta(H),p(init+H) 234.44 5.37 0.06 11 212.44 Psi(initSeason+Elev+ElevQ),R(FC+GC),dlta(H),p(init+H) 235.24 6.17 0.04 12 211.24 Psi(initSeason+Elev+Precip+P*E),R(.),dlta(H),p(init+H) 235.83 6.76 0.03 11 213.83 Psi(initSeason+Elev),R(FC),dlta(H),p(init+H) 238.28 9.21 0.01 10 218.28 Psi(initSeason+Elev),R(GC),dlta(H),p(init+H) 238.59 9.52 0.01 10 218.59 Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 12 500 m and 700 m (Fig. 3) and were associated with average monthly precipitation varying from 130 mm to 165 mm (Fig. 4). The influence of ground cover on the probability of encountering ≥2 individuals was positive and strong (b = 1.48 ± 0.74), whereas the influence of contextual forest cover was negative and weak (b = -2.35 ± 1.43). The influence of ground cover on the probability of encountering ≥2 individuals was illustrated by predicting probability values (R) after keeping forest cover constant (Fig. 5). On average, the probability of correctly identifying a station as state 2 (≥2 individuals) was 0.89 ± 0.06, and varied as a function of humidity Figure 3. Occupancy probability (± SE) for (A) Eleutherodactylus wightmanae and (B) E. brittoni as a function of elevation (m) in forested habitats on southern forested slopes in west-central Puerto Rico. Estimates were obtained from the best-supported model with a strong precipitation*elevation influence on occupancy. Caribbean Naturalist 13 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 Figure 4. Occupancy probability (± SE) for (A) Eleutherodactylus wightmanae and (B) E. brittoni as a function of long-term average monthly precipitation on southern forested slopes in west-central Puerto Rico. Estimates were obtained from the best-supported model with a strong precipitation*elevation influence on occupan cy. Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 14 levels per sampling occasion (b = 1.71 ± 0.79). Discussion We report occupancy probability and the probability of encountering ≥2 individuals at occupied survey sites of E. wightmanae and E. brittoni along elevational gradients on southern mountain slopes in west-central Puerto Rico. Our interest in these parameters was justified because they vary as a function of habitat conditions and serve as measures of success to gauge the value of habitat management actions (Nichols and Williams 2006, Yackulic et al. 2015). This work also integrated the use of passive acoustic devices, an approach that is non-invasive, efficient, and Figure 5. Predicted probability (dots) and 95% confidence intervals (CIs; open triangles) of encountering ≥2 individuals of (A) Eleutherodactylus wightmanae as a function of canopy cover and of (B) E. brittoni as a function of ground cover (scale 1–5) on southern mountain slopes of west-central Puerto Rico. Caribbean Naturalist 15 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 suited to generating reliable survey data for parameter estimation of rare, hard-toidentify species (Campos-Cerqueira and Aide 2016). Occupancy of both species was unimodal with respect to elevation. Highest probabilities were recorded between ~500 m and 700 m for E. brittoni and starting at ~400 m for E. wightmanae. These findings are in agreement with elevational ranges reported by Drewry and Rand (1983), but extend the elevational range inhabited by E. brittoni, previously believed to occur ≤700 m. Both species occurred at elevations up to 865 m, albeit occupancy probabilities were 0.10–0.15 at those highest elevations. As predicted, the probability of encountering ≥2 individuals of E. brittoni was positively influenced by ground cover, a finding that matches known habitat preferences (herbaceous; Drewry and Rand 1983, Joglar 1998). Canopy cover between 20% and 70% yielded the highest probability of encountering ≥2 individuals of E. wightmanae. A plausible explanation for this relationship is that partial cover promotes a heterogeneous understory, thereby enhancing food availability, shelter, and suitable micro-climatic conditions (Huang et al. 2014, Tews et al 2004, Whitfield et al. 2014). For instance, Drewry and Rand (1983) noted that E. wightmanae climbs low understory structures (to heights of 30 cm above the ground), perhaps expanding the number of types of calling sites where multiple individuals occur. The unimodal occupancy patterns suggested that habitat at intermediate elevation and precipitation amounts were most suitable for both species. Lower occupancy values recorded at lower elevation sites might reflect limits imposed by warmer and drier conditions. Admittedly, factors driving lower occupancy at higher elevations were less clear. It is not known if higher precipitation and lower temperatures impose physiological constraints, factors that do not preclude the occurrence of E. portoricensis at higher elevation (Barker and Rios-Franceschi 2014). It is also possible that factors not measured in this study, such as interspecific interactions, contributed to lower occupancy rates. Contrary to predictions, neither occupancy nor the probability of encountering ≥2 individuals of either species was influenced by seasonal weather covariates. This finding is likely due to the fact that the spread of temperatures and humidity values recorded during the study were within the long-term average monthly minima (13 °C) and maxima (32.5 °C) recorded over 50 years (Hijmans et al. 2005). Seasonal effects would likely require values closer to the species’ tolerance limits. Although thermal limits for E. wightmanae and E. brittoni are not known, recorded temperatures were within known minima and maxima for E. coqui (7.2–37.2 °C) and E. portoricensis (7.5–36.3 °C) (Christian et al. 1988). The conservation implications of this study are twofold. First, it provides a stronger foundation to implement habitat management actions. For example, lower local extinctions rates for species like E. wightmanae, listed as endangered in the IUCN Red List (Angulo 2008), could be achieved if actions target canopy cover, and thus abundance, at sites predicted to have high occupancy rates. Our index of abundance is a step forward in guiding conservation actions and evaluating species-specific responses to them, but we acknowledge that expanding the number of abundance categories or states (e.g., non-presence, 1, 2, 3–4, ≥5) should be Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 16 adopted as soon as advances in acoustical analyses make it possible. Second, we detailed a methodological framework suitable to test hypotheses about processes and mechanisms that might constrain abundance and distribution over a wide range of environmental conditions. One possibility is exemplified by long-term climatic changes, which in Puerto Rico have been shown to contribute to population declines and restricted distributions of Eleutherodactylus spp. (Burrowes et al. 2004, 2008; Pounds and Crump 1994; Pounds et al. 2006). Actions addressing conservation challenges posed by habitat and climatic changes are better informed with increased knowledge about species-habitat relationships, particularly if derived from dynamic (multi-season) analytical frameworks like the one used in this study (Frishkoff et al. 2015, Kearney and Porter 2009, McDonald-Madden et al. 2011, Watling and Donnelly 2008, Watling et al. 2009, Yackulic et al. 2015). Acknowledgments This work was supported by a grant from the US Geological Survey, Science Support Program in collaboration with the US Fish and Wildlife Service. We want to thank D. Hardgrove, M. Rivera, A. Irizarry, and L. Rodriguez for their assistance with field work, and J. G. Monroe for assistance with data management. We also thank M. Aide for his advice on managing and modelling acoustic data, J. Hines for statistical advice, and M. Campos- Cerqueira for comments that improved earlier versions of the manuscript. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US government. Literature Cited Angulo, A. 2008. Eleutherodactylus wightmanae. The IUCN Red List of Threatened Species 2008. Available online at http://dx.doi.org/10.2305/IUCN.UK.2008.RLTS. T57056A11575451.en. Accessed 23 July 2016. Barker, B.S., and A. Rios-Franceschi. 2014. Population declines of Mountain Coqui (Eleutherodactylus portoricensis) in the Cordillera Central of Puerto Rico. Herpetological Conservation and Biology 9:578–589. Burnham K.P., and D.R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Second Edition. Springer Science and Business Media. Fort Collins, CO, USA. Burrowes, P.A., R.L. Joglar, and D.E. Green. 2004. Potential causes for amphibian declines in Puerto Rico. Herpetologica 60(2):141–154. Burrowes, P.A., A.V. Longo, and C.A. Rodríguez. 2008. Potential fitness cost of Batrachochytrium dendrobatidis in Eleutherodactylus coqui, and comments on environmentrelated risk of infection. Herpetotropicos 4:51–57. Calvesbert, R.J. 1970. Climate of Puerto Rico and the US Virgin Islands. Climatography of the United States No. 60–52. US Department of Commerce, Environmental Science Services Administration, Environmental Data Service, Washington, DC, USA. Campos-Cerqueira, M., and T.M. Aide. 2016. Improving distribution data of threatened species by combining acoustic monitoring and occupancy modelling. Methods in Ecology and Evolution 7:1340–1348. Christian, K.A., F. Nunez, L. Clos, and L. Diaz. 1988. Thermal relations of some tropical frogs along an altitudinal gradient. Biotropica 20:236–239. Caribbean Naturalist 17 K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 Cooch, E., and G. White. 2010. Program MARK: A gentle introduction. Colorado State University, Fort Collins, CO, USA. Dowdy, K.E. 2016. Occupancy and Abundance of Eleutherodactylus frogs in coffee agroecosystems and along an elevational gradient in the mountains of southwestern Puerto Rico. M.Sc. Thesis. North Carolina State University, Raleigh, NC, USA. Drewry, G.E., and A.S. Rand. 1983. Characteristics of an acoustic community: Puerto Rican frogs of the genus Eleutherodactylus. Copeia 4:941–953. Ewel, J.J., and J.L. Whitmore. 1973. Ecological life zones of Puerto Rico and US Virgin Islands. USDA Forest Service research paper ITF-18. Institute of Tropical Forestry, Rio Piedras, PR, USA. December 1973. Frishkoff, L.O., E.A. Hadly, and G.C. Daily. 2015. Thermal niche predicts tolerance to habitat conversion in tropical amphibians and reptiles. Global Change Biology 21:3901–3916. Gonser, R.A., and L.L. Woolbright. 1995. Homing behavior of the Puerto Rican frog Eleutherodactylus coqui. Society for the Study of Amphibians and Reptiles 29: 481–484. Henareh Khalyani, A., W. Gould, E. Harmsen, A. Terando, M. Quinones, and J. Collazo. 2016. Climate change implications for tropical islands: Interpolating and interpreting statistically downscaled GCM projections for management and planning. Journal of Applied Meteorology and Climatology 55:265–282. Hijmans, R.J., S.E. Cameron, J.L. Parra, P.G. Jones, and A. Jarvis. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25:1965–1978. Hines, J.E. 2006. PRESENCE – Software to estimate patch occupancy and related parameters. USGS-PWRC. Available online at http://www.mbr-pwrc.usgs.gov/software/presence. html. Accessed 14 September 2015. Huang, S.-P., W.P. Porter, M.-C. Tu, and C.-R. Chiou. 2014. Forest cover reduces thermally suitable habitats and affects responses to a warmer climate predicted in a high-elevation lizard. Oecologia 175:25–35. Intergovernmental Panel on Climate Change (IPCC). 2014. Climate change 2014 synthesis report. Contribution of Working Groups I, II, and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Available online at https://www.ipcc. ch/report/ar5/syr/. Joglar, R.L. 1998. Los Coquíes de Puerto Rico: Su Historia Natural y Conservación. Editorial de la Universidad de Puerto Rico, San Juan, PR, USA. Kearney, M., and W. Porter. 2009. Mechanistic niche modelling: Combining physiological and spatial data to predict species’ ranges. Ecology letters12:334–350. Kearney, M., R. Shine, and W.P. Porter. 2009. The potential for behavioral thermoregulation to buffer ‘‘cold-blooded’’ animals against climate warming. PNAS 106(10):3835–3840. Lips, K.R., P.A. Burrowes, J.R. Mendelson III, and G. Parra-Olea. 2005. Amphibian population declines in Latin America: A Synthesis. Biotropica 37:222–226. MacKenzie, D., J.D. Nichols, G.B. Lachman, S. Droege., J.A. Royle, and C.A. Langtimm. 2002. Estimating site occupancy rates when detection probabilities are less than one. Ecology 83:2248–2255. MacKenzie, D.I., J.D. Nichols J.A. Royle, K.H. Pollock, L.L. Bailey, and J.E. Hines. 2006. Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence. Academic Press, Cambridge, MA, USA. 324 pp. McCallum, M.L. 2007. Amphibian decline or extinction? Current declines dwarf background extinction rate. Journal of Herpetology 41:483–491. Caribbean Naturalist K.D. Monroe, J.A. Collazo, K. Pacifici, B.J. Reich, A.R. Puente-Rolón, and A.J. Terando 2017 No. 40 18 McDonald-Madden, E., M.C. Runge, H.P. Possingham, and T.G. Martin. 2011. Optimal timing for managed relocation of species faced with climate change. Nature Climate Change (letters) 1:261–265. Nichols, J.D., and K.B. Williams. 2006. Monitoring for conservation. Trends in Ecology and Evolution 21:668–673. Nichols, J.D., J.E. Hines, D.I. MacKenzie, M.E. Seamans, and R.J. Gutierrez. 2007. Occupancy estimation and modeling with multiple states and state uncertainty. Ecology 88:1395–1400. Nudds, T.D. 1977. Quantifying the vegetative structure of wildlife cover. Wildlife Society Bulletin 5:113–117. Ospina, O.E., L.J. Villanueva-Rivera, C.J. Corrada-Bravo, and T.M. Aide. 2013. Variable response of anuran calling activity to daily precipitation and temperature: Implications for climate change. Ecosphere 4(4):47. Available online at http://dx.doi. org/10.1890/ ES12-00258.1. Parés-Ramos, I.K., W.A. Gould, and T.M. Aide. 2008. Agricultural abandonment, suburban growth, and forest expansion in Puerto Rico between 1991 and 2000. Ecology and Society, 13(2):1. Pollock, K. 1982. A capture–recapture sampling design robust to unequal catchability. The Journal of Wildlife Management 46:752–757. Pounds, J.A, and M.L. Crump. 1994. Amphibian declines and climate disturbance: The case of the Golden Toad and the Harlequin Frog. Conservation Biology 8:72–85. Pounds, J.A., M.R. Bustamante, L.A. Coloma, J.A. Consuegra, M.P.L Fogden, P.N. Foster, and B.E. Young. 2006. Widespread amphibian extinctions from epidemic disease driven by global warming. Nature 439:161–167. Radeloff, V.C, et al. 2015. The rise of novelty in ecosystems. Ecological Applications 25(8):2051–2068. Rivero, J.A. 1978. Los Anfibios y Reptiles de Puerto Rico (The Amphibians and Reptiles of Puerto Rico) (First Edition). Editorial de la Universidad de Puerto Rico, San Juan, PR, USA. Tews, J., U. Brose, V. Grimm, K. Tielbörger, M.C. Wichmann, M. Schwager, and F Jeltsch. 2004. Animal species diversity driven by habitat heterogeneity/diversity: The Importance of keystone structures. Journal of Biogeography 31(1):79–92. Wake, D.B., and V.T. Vredenburg. 2008. Are we in the midst of the sixth mass extinction? A view from the world of amphibians. Proceedings of the National Academy of Sciences 105(Supplement 1):11466–11473. Watling, J.I., and M.A. Donnelly. 2008. Species richness and composition of amphibians and reptiles in a fragmented forest landscape in northeastern Bolivia. Basic and Applied Ecology 95:523–532. Watling, J.I., K. Gerow, and M.A. Donnelly. 2009. Nested species subsets of amphibians and reptiles on Neotropical forest islands. Animal Conservation 125:467–476. Whitfield, S.M., K. Reider, S. Greenspan, and M.A. Donnelly. 2014. Litter dynamics regulate population densities in a declining terrestrial herpetofauna. Copeia 3:454–461. Yackulic, C.B., J.D. Nichols, J. Reid, and R. Der. 2015. To predict the niche, model colonization and extinction. Ecology 96(1):16–23.