Demography and homing behavior in the poorly-known Philippine flat-headed frog Barbourula busuangensis (Anura: Bombinatoridae)

Study area and CMR monitoring

This study was conducted in the island of Busuanga, in the Calamian Archipelago, province of Palawan, Philippines (Figs. 2A–2B). Busuanga has an area of 915.16 km² and a maximum elevation of 620 m a.s.l. It is known for its rich and diverse vegetation, including, in order of abundance, early secondary growth forests, advanced secondary growth forest and old growth forests, land used for agricultural purposes, and mangroves (Supsup et al., 2023). The University of the Philippines at Los Baños and the Palawan Council for Sustainable Development approved the permit to conduct fieldwork in Busuanga (Wildlife Gratuitous Permit N° 2022-11).

We conducted field work over a total of 10 months, divided in three periods during 2022 and 2023, corresponding to April–July 2022, October–December 2022, and April–June 2023. We worked at two locations along rocky rivers surrounded by dense tropical rainforest (Fig. 1C). Our main study site, Malbato (Barangay Bintuan, 12°2′14″N, 120°5′49″E) (Fig. 2B, orange circle) is in the southern part of the island, likely near the vaguely described type locality of the species (Taylor & Noble, 1924). Here we established three transects adding up to 420 m (transect “A” = 100 m, transect “B” = 100 m and transect “C” = 220 m, Fig. 2C). Transects were marked with flagging tape every 10 m to identify specific locations where frogs were captured. Due to logistical reasons, we sampled in Malbato more frequently for a total of 40 nights throughout our study period, and on average surveys were conducted by two or three people. The second site, San Rafael, is in the north-western part of Busuanga, in the municipality of San Rafael (Barangay New Busuanga, 12°12′12″N, 119°55′9″E) (Fig. 2B). At this site, we monitored Barbourula busuangensis along a 320 m transect once a month during each field period, for a total of nine nights (Fig. 2C). Field parties at this site were always composed of four people such that in a way, the man-effort at both localities was balanced. The mean temperature during our work in Busuanga was 28 °C (range 26–31 °C) and the mean daily precipitation varied from 4.86–7.13 mm/day (data obtained from visualcrossing.com).

We monitored our transects along the river from 6:00 PM to approximately 3:00 AM, aiming to detect and capture every Barbourula busuangensis in the study area. For every captured individual, we measured body size as the snout-to-vent length (SVL) to the nearest 0.01 mm using a manual caliper. Individuals measuring more than 35 mm were marked with an 8 mm mini-Passive Internal Transponder (PIT) (ID 100A/1.4; Trovan, Douglas, Isle of Man, UK). To insert the transponder, we lifted the skin of the upper dorsal area of the frog with fine-tip tweezers and used surgical scissors to make a small incision. We placed the loaded canula (NDL-1xx/1.4; Trovan, Douglas, Isle of Man, UK) in the pistol grip implanter (IM-300C; Trovan, Douglas, Isle of Man, UK) at the incision point and inserted the microchip under the skin. Then, we used the blunt end of the tweezers to push the transponder posteriorly, close to the urostyle area. We did not seal the wound after insertion because experience with other tropical frogs proved that using this method yielded high transponder retention and no signs of infection (Burrowes, unpublished data). Each transponder has a unique digital code that can be detected with a portable reader (LID-573; Trovan, Douglas, Isle of Man, UK). Frogs with SVL < 35 mm were considered too small to be tagged with transponders and were marked with visual internal elastomers (VIE) instead. VIEs consist of a UV-fluorescent polymer that is injected subcutaneously using a 29–gauge syringe (Northwest Marine Technology Inc., Anacortes, WA, USA).

Age structure

To assess the age structure of the populations at both study sites, we classified individuals in six categories based on their SVL: Juveniles-I (individuals with SVL < 25 mm), Juveniles-II (SVL = 25–34 mm), Subadults (SVL = 35–49 mm), Adults-70 (SVL = 50–70 mm), Adults-90 (SVL = 71–90 mm) and Adults-100 (SVL > 90 mm). Because Barbourula busuangensis lacks obvious morphological features signaling sexual maturity, our criterion for determining adulthood was informed by the size of the smallest gravid female observed during our study (SVL = 51.2 mm); then, the younger age classes were subjectively estimated. Although skeletochronology (see for example Castanet & Smirina, 1990; Roček et al., 2016; Rahman et al., 2022) or dissection to observe gonad development (see for example López et al., 2017; Piazza et al., 2023) are more accurate methods to assess age and adulthood respectively, they require sacrificing individuals. Our work was limited in this respect, because we did not have permits to collect animals for the preservation. We applied different VIE colors (coded by fieldwork period, Table 1) to distinguish younger juveniles (Juveniles-I) from larger, and presumably older ones (Juveniles-II). These different color codes allowed us to recognize and track young cohorts through their growth until they reached subadult age class, at which we marked them with PITs (Peña-Jiménez & Burrowes, 2021).

Since Barbourula busuangensis also lacks external sexual dimorphism, we could only confirm the sex of gravid females when mature oocytes were visible through the ventral skin (Bosch et al., 2023). Consequently, we could not compare results by sex or determine sex biases in our sampling. To assess if the age structure was different between the southern (Malbato) and northern (San Rafael) populations, we conducted a Chi-square test of independence by sampling site.

Growth rates and age estimates

We estimated growth rates of post-metamorphic frogs as the increase in body size (SVL) divided by the time lapse between captures (mm/month) for both cohorts of juveniles, and for adult and subadult frogs marked with transponders. We applied a generalized lineal model (GLM) to explore the relationship between the SVL of adult and subadult frogs at first capture (explanatory variable), and their growth rate per month at subsequent recaptures (response variable). Because the response variable (monthly growth rate) adjusted to a gamma distribution, we implemented the family = gamma (link = “log”) arguments within the glm function in “R” statistical package to run the model, and qqplots and moments package to check the normality of residuals (R Core Team, 2020). We used average growth rates by age-class category to extrapolate the approximate age (in years) of Barbourula busuangensis of specific body sizes.

Population size estimates

We used MARK software (White & Burnham, 1999) to analyze the capture histories of all frogs marked during the study period. Specifically, we applied the POPAN formulation (Schwarz & Arnason, 1996) to estimate population size. This formulation is appropriate for the analysis of open populations where births, deaths, and movement of individuals in and out of the study area are likely to occur (Capellà-Marzo, Sánchez-Montes & Martínez-Solano, 2020; Fernández de Larrea et al., 2021). The POPAN models include four main parameters: (1) apparent survival of individuals (ϕ), (2) probability of capture (p), (3) rate of entrance of new individuals into the study area (pent), and (4) total population size, accounting for all the individuals present at the survey area during the study period (N). We tested a set of a priori biologically plausible models, with the parameters ϕ, p and pent as either constant (.) or time dependent (t) for both study sites. We analyzed the CMR datasets in two ways; first, to compare overall estimates between the north and south populations of Barbourula busuangensis, we analyzed each population separately by grouping individuals of all age classes over the three sampling periods at each site. Secondly, and only for Malbato, where the higher frequency of sampling occasions allowed for finer scale analysis, we divided individual capture histories by the three-field sampling periods, and included age-classes defined by body size as adults (>50 mm SVL) and sub-adults (35–49 mm SVL). In this detailed analysis at Malbato, ϕ, p and pent were modelled as either constant (.), dependent on age-class group (g), time (t) or their interaction (g*t). The likelihood of the tested models was assessed by Akaike’s information criterion corrected for small sample sizes (AICc), and estimates of population size (N) were obtained for each site based on weighted averaging among the AICc-ranked models. We used the software U-CARE to validate model assumptions via goodness of fit (GOF) tests (Choquet et al., 2009). The two most common sources of assumption violations are known as transience and trap-dependence effects, and both were examined using the tests 3.SR and 2.CT in UCARE, respectively. Transience effects caused by dispersing individuals violate the assumption that all frogs in the population have the same probability of being captured, while trap-dependence considers that marking has some effect (either positive or negative) on the probability of subsequently recapturing an individual (Pradel & Sanz-Aguilar, 2012).

Potential for homing behavior

We conducted displacement experiments at Malbato to test for the potential of homing behavior, defined as the probability of an individual to return to its original capture site after being displaced. Only marked frogs that had been recaptured at least once at the same location, were regarded as residents and considered for displacement experiments (Crump, 1986). We first moved 32 individuals 10 m away from their capture site; specifically, we translocated 10 adults and six subadults 10 m upstream, and nine adults and seven subadults 10 m downstream. We considered that an individual had exhibited homing behavior when recaptured within a 5 m radius of its original capture site. Once we had evidence that some individuals were returning to the original capture site, we extended the displacement distance to 30 m and 50 m during our second and third sampling periods. A total of 54 adults and 28 subadults were displaced at increasing distances (10, 30 and 50 m) upstream and downstream. During displacements, we collected frogs in transparent bags and carried them within a larger black bag. We did not run a control experiment to test if frog manipulation (i.e., capture and bagging) influenced the probability of homing. We used a multiple logistic regression model in R (R Core Team, 2020) to assess the effect of age class (adult or subadult), direction of displacement (upstream or downstream), and displacement distance (10, 30 or 50 m) on the probability of homing.

Permits to conduct this research in the field were obtained from Palawan Council for Sustainable Development (PCSD, Republic of the Philippines) and Ethics Competence B&C from Consejería de Medio Ambiente y Ordenación del Territorio de la Comunidad de Madrid, Spain. Ref. 10/096442.9/13.

Capture-mark-recapture histories

We registered a total of 1,251 capture events at Malbato, including all age classes. At this site, we captured and marked with transponders a total of 196 individuals of Barbourula busuangensis corresponding to adult (143) and subadult (53) age classes (SVL > 35 mm). Among them, 73 adults and 23 subadults (49%) were recaptured only once, 42 adults and seven subadults were recaptured twice, 13 adults and five subadults were recaptured three times, five adults and six subadults were recaptured four times, and 10 adults and 12 subadults were recaptured more than five times, including one subadult that was recaptured 14 times (Fig. 3A). At San Rafael, we recorded a total of 1,042 capture events across all age classes and marked 144 individuals with transponders corresponding to 112 adults and 32 subadults. Of these, 74 adults and 15 subadults (60%) were recaptured only once, 25 adults and nine subadults were recaptured twice, seven adults and six subadults were recaptured three times, and six adults and one subadult were recaptured four times, with a maximum of seven recaptures for a single subadult (Fig. 3B).

We captured more adults than subadults at both sites (Figs. 3A–3B), and adults comprised 73% and 78% of the recaptured individuals in Malbato and San Rafael, respectively. The total number of animals captured and recaptured increased with our sampling effort and recapture rates per sampling month (total number of individuals recaptured over total marked) approximated 50% of the population at both sites by the end of our study (Figs. 3C–3D). Thus, marking with transponders proved to be a safe way to track Barbourula busuangensis in its riverine habitat. We did not observe any signs of losing the transponders contrary to findings in frogs of the genus Litoria (Brannelly, Berger & Skerratt, 2014), or infection at the incision point despite not sealing the wound. Moreover, recaptured frogs had sealed the wound between 2–4 days of initial marking.

VIE proved to be a safe marking method to follow young froglets, too small to withstand transponders. We captured and marked with VIE 68 juveniles (<35 mm SVL) in Malbato and recaptured 44 (65%) either in the same or in subsequent field seasons (Table 1). Eight of the juveniles-II marked with blue VIE during the first sampling season were recaptured during the second season; six of these maintained a body size of SVL < 35 mm, and thus, were classified again as juveniles-II and re-marked with a pink VIE, whereas the other two were PIT-tagged as subadults because they had grown to SVL ≥ 35 mm (Table 1). Ten juveniles-II marked with pink VIEs during the second field season were recaptured during the third season, still at this age class, and re-marked with a purple VIE. In San Rafael, we marked 36 juveniles and recaptured only two of them during the first field season (Table 1).

Age structure

The proportion of individuals at each age class was independent of study site (X² (5, 450) = 0.61, p = 0.96), yielding a similar age class distribution for Malbato and San Rafael (Fig. 4). Adults of SVL between 50 and 90 mm were the dominant age-class at both sites. In contrast, we recorded very few adults larger than 90 mm in SVL (Adults-100) during our study, with only four in Malbato and one in San Rafael. We observed 23 individuals that we could accurately identify as adult females because they had visible eggs from the venter indicating that they were gravid. The size of these females ranged from SVL = 50–69 mm (N = 5), SVL = 70–79 mm (N = 18), and SVL = 80–90 mm (N = 4). Despite the greater sampling effort at Malbato, we found a similar number of subadults at both sites, 44 in Malbato and 33 in San Rafael (Fig. 4). We recorded almost the same number of Juvenile-II individuals for both sites (27 in Malbato and 28 in San Rafael), but almost twice as many Juvenile-I in Malbato than in San Rafael (Fig. 4).

Growth rates and age estimates

For cohorts of VIE-tagged juveniles that were recaptured in subsequent field seasons in Malbato we estimated average growth rates of 4.7 mm (range: 0.4–11 mm) in SVL over a period of 3–4 months (1.21–1.50 mm/month). Thus, if we consider the smallest froglet observed (SVL = 10.3 mm) a recent metamorph, we can infer that it takes a very young juvenile of Barbourula busuangensis approximately 17 months (roughly 1.5 years) to reach a SVL of 35.8 mm equivalent to our subadult age class (Table 2). For larger individuals marked with transponders, monthly growth rate decreased significantly with SVL (coefficient = −0.035, p < 0.0001) (Fig. 5, Appendix S1 for output table). Subadults grew at a mean growth rate of 1.56 mm/month (95% C.I [0.17–2.64]), while adults grew at a slower rate, mean of 0.60 mm/month (95% C.I. [0.04–1.92]) (Fig. 5).

Using the mean growth rates obtained per age class we extrapolated the age of metamorphosed individuals by body size (SVL). In this way, a young subadult of 36 mm and 1.5 years old (according to mean juvenile growth rates-see above), growing at an average rate of 1.56 mm/month, may take 12 months to grow to a SVL of 54 mm, thus, reaching Adult-70 category roughly at 2.5 years of age. Then, it would take this frog about 2.5 more years, growing at mean adult growth rate of 0.6 mm/month, to reach 72 mm in SVL (corresponding to a young Adult-90 category), and another 4–5 years to reach SVL in the range of 81–90 mm (advanced Adult-90), making the oldest individuals in the population over 11 years old (Table 2).

Population size estimates

Considering the data for all individuals over the entire study period, the best-ranked model at each sampling site was that where apparent survival (ϕ) was kept constant, and the probability of entrance of individuals (pent) was variable across time (Table 3). However, while the probability of capture (p) was kept constant in the best model for Malbato, it was variable across time for San Rafael (Table 3; Appendix S2 for all models). The model-averaged population size in Malbato (N = 268, 95% confidence interval CI [242–293]) was slightly higher than that for San Rafael (N = 232, 95% CI [191–273], Table 3). These estimates indicate that we captured 73% (196/268) and 62% (144/232) of the population at both sites respectively. The estimated survival rate was high throughout our study, with a mean daily survival rate of 0.99 for both sites (Appendix S3 and S4). The probability of capture of individuals across the sampling sessions was notably higher in San Rafael (0.33), compared to Malbato (0.09, Appendix S3 and S4). On the other hand, the rate of entrance of new individuals into the study area was low for both sites (Appendix S3 and S4).

When age class (adults and subadults) was included in the model and the three sampling periods were considered separately for Malbato, the best-ranked models were those where apparent survival (ϕ) was kept constant, the probability of capture (p) was dependent on the age class, and probability of entrance of individuals (pent) was variable across time and age classes (Table 4; see also Appendix S5 for all models). In the three sampling periods, a lower-ranked model showed age class dependent survival (Table 4). Daily survival rate in Malbato was concordant across all models, showing values ≥0.97 for all age classes in all sampling periods (Appendix S6, S7 and S8). The probability of capture, thus detectability, varied with age class, being on average higher for subadults (0.28) compared to that of adults (0.09). Finally, the mean probability of entrance of new individuals to the population was low across all sampling seasons and age classes (<0.001–0.265), but was generally larger for subadults than adults (Appendix S6, S7 and S8).

The abundance of Barbourula busuangensis estimated by this model varied between sampling periods for both age classes, but the model averaged estimates fall within the confidence intervals of the next sampling period (Table 4). While the average estimated abundance of subadults increased with time from the first to the last field season (N = 29–40), that of adults was highest in October–December 2022 (N = 185).

The GOF (goodness of fit) tests of the POPAN datasets for all individuals over the entire study period showed evidence of transience and trap-dependence effects in Malbato (Appendix S9). However, these effects were not consistently maintained when the data were analyzed by age-classes and field seasons. In this case, we only detected transience effects for subadults in two of the sampling periods, whereas no trap-dependence effect was found for any group or season (Appendix S9). We did not find any departures from POPAN model assumptions at San Rafael (Appendix S9).

Potential for homing behavior

We recaptured 53 of the 82 marked individuals that were displaced upstream or downstream, at time intervals between four and 353 days after their translocation. A total of 42 of the translocated individuals (51%) eventually returned to their original capture site, suggesting the probability of homing behavior in this species (Table 5). Seven individuals (8%, four adults and three subadults) were found elsewhere along the river, at an average distance of 63 m (range: 10–140) from their original capture site. The remaining 29 individuals (31%) were not recovered by the end of the study. Results of the multiple logistic regression model (Model: Homing ~ Age + Distance + Direction) revealed that the probability of homing was not significantly associated to age class (adults vs. subadults), distance of displacement (10, 30 or 50 m) or the direction of river flow (upstream vs. downstream) (Appendix S10).