Occupancy modelling reveals a highly restricted and fragmented distribution in a threatened montane frog (Philoria kundagungan) in subtropical Australian rainforests

Keywords: climate change, geographic range, habitat preference, montane amphibians, occupancy modelling, Philoria kundagungan, subtropical rainforest, threatened species.

Habitat loss and degradation, overexploitation, invasive species, disease, pollution and, more recently, climate change, have caused a rapid reduction in global biodiversity (Barnosky et al. 2011; Ceballos et al. 2017; Scheele et al. 2019b). Despite efforts to reduce greenhouse gas emissions, the threat of climate change to biodiversity is increasing and will continue for at least the next several hundred years (Pachauri et al. 2014). The documented effects of climate change on species include shifts in geographic range, changes to the timing of breeding activity, increased frequency of catastrophic events and ultimately further population declines and species extinctions (Foden et al. 2013; Pachauri et al. 2014; Harris et al. 2018). Understanding the ecological responses to climate change allows natural resource managers to better anticipate these changes and develop effective conservation strategies (Tilman et al. 2017; Hisano et al. 2018).

Montane habitats contain unique fauna that are likely to be disproportionately threatened under future climate scenarios. Species occupying montane stream habitats are exposed to threats from climate change that are generally not shared with lowland species (Hu and Riveros-Iregui 2016; Palomo 2017; Alexander et al. 2018). Warming is predicted to drive range shifts to higher elevations, as opposed to the latitudinal shifts experienced by species occupying lowland habitats (Chen et al. 2011). The climates of montane habitats are often substantially different from those of adjacent low-elevation environments, which act as dispersal barriers to montane species (Hu and Riveros-Iregui 2016; Palomo 2017; Alexander et al. 2018). Species with poor dispersal abilities are often isolated in these ‘islands’ of habitat and only have the option of up-slope range shifts to escape increasing temperatures (Pounds et al. 1999; Chen et al. 2011; Enriquez-Urzelai et al. 2019). For example, Raxworthy et al. (2008) observed mean increases in the elevational midpoint range of 30 species of reptiles and amphibians in response to warming between 1993 and 2003; the result was a substantial reduction in the available habitat of these species. Montane ecosystems contain a high diversity of species and are considered to be one of the most at-risk ecosystem types to changing climates (Raxworthy et al. 2008; Laurance et al. 2011a; Enriquez-Urzelai et al. 2019). However, due to a lack of effective monitoring, the basic conservation status of these species, such as geographic distribution and population trend, are poorly known, particularly for amphibians (Scheele et al. 2019a).

Because of their higher sensitivity to changes in temperature and moisture than other vertebrates, amphibians present a valuable opportunity to monitor upslope range shifts and range reductions in response to climate change. Amphibians are ectothermic, lay unshelled eggs, and have biphasic life cycles with the timing of breeding events commonly driven by temperature and precipitation cues (Wells 2010). Montane amphibians are likely to be more sensitive to anthropogenic climate change compared with low-elevation species because, on average, they are more narrowly distributed (Laurance et al. 2011a), have smaller population sizes, rendering them more vulnerable to decline when exposed to diseases (Scheele et al. 2019b), and have a greater vulnerability to genetic isolation due to dispersal barriers and abilities (Parmesan 2006). Additionally, amphibians are the most threatened class of vertebrates, with at least 2190 species estimated to be in danger of extinction globally (IUCN 2019). In Australia, 45 (19%) amphibian species are considered to be threatened or extinct according to International Union for Conservation of Nature (IUCN) Red List criteria (Gillespie et al. 2020). Habitat loss, infectious disease (particularly chytridiomycosis), climate change, fire regimes and invasive species are the leading factors driving amphibian population declines in Australia (Gillespie et al. 2020). Monitoring of threatened amphibian species is essential to understand the dynamics of range shifts in montane species in response to major threats such as climate change, invasive species, land-use change and disease, and to develop effective conservation strategies (Stuart et al. 2004; Magurran et al. 2010; Scheele et al. 2019a). However, monitoring of threatened amphibian species in Australia has broad deficiencies and is presently inadequate for informing long-term species conservation (Scheele et al. 2019a).

One group of Australian frogs that are particularly susceptible to climate change are the six species described within the genus Philoria. All are range-restricted species, have montane distributions and are listed as endangered or critically endangered by the IUCN (IUCN 2019). Five species (P. kundagungan, P. loveridgei, P. pughi, P. richmondensis and P. sphagnicola) occur in northern New South Wales and south-eastern Queensland. They are found in association with seepages or the boggy margins of the headwaters of drainage lines, mostly in upland rainforest (Knowles et al. 2004). Philoria frosti occurs in alpine and subalpine habitat on Mount Baw Baw, Victoria (Anstis 2017). All species lay their eggs in small chambers excavated in mud, leaf litter or under rocks (Knowles et al. 2004; Anstis 2017). Adults are secretive and are rarely encountered outside of their subterranean breeding chambers. Detection relies entirely upon the advertisement call of males (from within burrows), which occur diurnally from early spring through to summer (Knowles et al. 2004; Willacy et al. 2015). For most Philoria species, factors determining their distribution, abundance and detectability are very poorly known.

We examine the occupancy probability, detection probability and geographic distribution of the very poorly known Philoria kundagungan, previously documented from only eight sites. Our aims were to identify drivers of site occupancy and detection probability, define the geographic boundaries of this species and thus provide robust baseline data for future species monitoring. We undertook surveys for P. kundagungan using an occupancy framework (MacKenzie et al. 2002) and investigated site occupancy as a function of measured environmental variables. We conducted repeat surveys at each site to account for imperfect detection and investigated the influence of environmental conditions on detection probability. Additionally, we mapped current potentially suitable habitat for P. kundagungan using the modelled relationship between site occupancy probabilities and site covariates. By providing baseline data for this species, future studies will be able to detect changes to the geographic distribution of P. kundagungan in response to both threatening processes and applied management strategies. Although this study is focussed on a single species, future changes to the geographic distribution of P. kundagungan could serve as an indicator for the health of the montane ecosystem it occupies, which is considered one of the most vulnerable ecosystems in Australia (Laurance et al. 2011b).

Philoria kundagungan is a small (up to 28 mm snout–vent length) terrestrial species with dorsal colouration yellow, orange or red against which red and black blotching is present. The ventral surface varies between yellow with red blotching and red with yellow and brown blotching (Ingram and Corben 1975; Knowles et al. 2004; Anstis 2017). It was described in 1975 from four sites (Ingram and Corben 1975), and subsequent published surveys only increased the number of known localities to eight, all within a small montane area between Mistake Mountains at the northern end of Main Range National Park in south-eastern Queensland south to Tooloom National Park in north-eastern New South Wales (Knowles et al. 2004; Anstis 2017). Previous records were mostly of males, calling during spring and summer, from the boggy margins of first- and second-order streams in rainforest and adjoining wet sclerophyll above 700 m in elevation.

This study was conducted in the mountainous regions of south-east Queensland and north-eastern New South Wales, Australia (27.70°–28.70°S, 152.18°–152.70°E). Surveys were conducted within Main Range, Koreelah, Mount Clunie and Tooloom National Parks, and Koreelah and Beaury State Forests (Fig. 1). These national parks and state forests cover ~460 km² and contain all known P. kundagungan localities. Elevation ranges from 50 to 1375 m; however, due to historical land clearing, areas below 600 m are largely cleared or otherwise degraded and were excluded from the study area. Genetic evidence suggests that none of the Philoria spp. occur sympatrically and that Philoria speciation events occurred ~12 million years ago during the break-up of the Gondwanan rainforests (Knowles et al. 2004). Accordingly, we excluded Yabbra, Toonumbar, Mount Nothofagus and Mount Barney National Parks for our study area as they are known to contain other allopatric species of Philoria (Knowles et al. 2004; S. Donnellan, unpubl. data).

Fig. 1.  Map showing the location of the 32 study sites and previously published records of Philora kundagungan in north-eastern New South Wales and south-eastern Queensland. Rainforest data source: DAWE (2019).

Thirty-two sites (Fig. 1) were surveyed three times between September 2016 and January 2017. Sites were randomly selected within the study area using ArcMap 10.3.1 under four constraints: (1) they contained a stream or creek >100 m in length; (2) they were within or <500 m from rainforest, (3) they were spaced >500 m apart to ensure site independence, and (4) they were <2 km from a road or operational fire trail (to facilitate efficient access). At each site, a 100 m long transect commenced at the highest point of the stream containing a soak. All surveys were performed during daylight hours as male calling activity is minimal at night (L. Bolitho, unpubl. data). Each transect was surveyed for a minimum of 15 min and call recordings of a male P. kundagungan were played through a hand-held speaker (Logitech, Ultimate Ears ROLL 2) every 10 m in order to elicit a response if frogs were not spontaneously calling. For each survey, the number of P. kundagungan detected was recorded. To examine the relationship between P. kundagungan detection and meteorological conditions, five covariates were recorded for each survey: (1) air temperature; (2) relative humidity; (3) precipitation level; (4) precipitation in the previous 24 h; and (5) average minimum temperature for the previous five days. Six site covariates were also recorded: (1) elevation; (2) aspect; (3) average annual temperature; (4) average annual precipitation; (5) vegetation type; and (6) presence/absence of rainforest spinach (Elatostema reticulatum), an herbaceous plant that occurs in very moist areas and is considered a possible indicator of suitable P. kundagungan habitat (Knowles et al. 2004).

Elevation and aspect data were extracted from a digital elevation model grid of Australia (Geoscience Australia 2015) using ArcMap 10.3.1. Average annual temperature and average annual precipitation data were extracted from the New South Wales and Australian Capital Territory Regional Climate Modelling dataset (Evans et al. 2014). Precipitation on the day before a survey and the average minimum temperature for five days prior were extracted from a set of daily meteorology models (Jones et al. 2009) using ArcMap 10.3.1. At the time of survey, relative humidity and temperature were measured with a Kestrel 4500 weather meter and precipitation was categorised into one of four levels (heavy, moderate, light and zero precipitation) by the observer.

Methods developed by (MacKenzie et al. 2018) were followed to estimate detection probabilities (p) and the probability of a site being occupied (ψ) as a function of measured site and sample covariates. The most supported model was identified with Akaike’s Information Criterion (AIC) and Akaike weight of evidence (w_i). As our sample size was relatively small, we have used AICc to correct for biases that can arise from a small sample size. All analyses were completed using single-species, single-season models in PRESENCE 11.8 (MacKenzie et al. 2002).

A two-step process was used to address a priori hypotheses relating to P. kundagungan occupancy. First, the six sampling covariates thought to influence detection probability were modelled with site occupancy fixed as constant. Second, the best-supported model for detection probability was combined with candidate models testing hypotheses predicting P. kundagungan occupancy. The best-supported model from this model set represented the final model explaining the relationship between P. kundagungan occupancy, detection probabilities and measured environmental variables. Prior to analysis, all covariates were tested for multicollinearity in SPSS Statistics 24 and all continuous covariates were normalised, by subtracting the mean and dividing by the standard deviation of each covariate. Model averaging parameter estimates were used to plot the relationship between P. kundagungan occupancy probability and site covariates. A goodness of fit test of the most parameterised model indicated overdispersion in the dataset (ĉ = 2.48), which was accounted for by adjusting initial AICc values by ĉ and comparing models via quasilikelihood AICc values (QAICc).

As the occupancy framework allows inference about species status at unsurveyed locations (MacKenzie et al. 2018), we used the modelled relationship between site occupancy probabilities and site covariates to map the geographic distribution of potentially suitable P. kundagungan habitat. The equation explaining the relationship between model averaged parameter estimates of P. kundagungan occupancy probability and site covariates was used to calculate occupancy probability from geographical information system layers of site covariates. The potential habitat map was constructed in ArcMap 10.3.1 using a pixel size of 30 m². As this species is known to occur only in or adjacent to drainage lines, spatial data were clipped to defined watercourses in the study area. Spatial data defining watercourses were obtained from digital 1 : 25 000 topographic map sheets published by the Queensland Department of Natural Resources, Mines and Energy and the New South Wales Government Spatial Services. Streams occurring in vegetation types other than wet sclerophyll and rainforest were excluded. The lower threshold value for areas considered as potential habitat was set as the lowest occupancy probability value of a site we found to contain P. kundagungan. To calculate the area of potential habitat, stream width was set to 20 m, ensuring that both the stream bed and banks were included in the area calculation. Potential habitat was classified into 10 quality classes relating to modelled P. kundagungan occupancy probability and the percentage within national parks was calculated for each class. To assess habitat connectivity, we calculated the number, size and distance between potential habitat fragments using ArcMap 10.3.1. Due to the likely poor dispersal ability of this species (Anstis 2017), areas of potential habitat spaced >1 km apart were considered to be separate habitat fragments.

Using methods described by Guillera-Arroita and Lahoz-Monfort (2012), we calculated the probability of detecting a significant (at α = 0.05) difference in P. kundagungan occupancy between our study (Season 1) and a reiteration of this study in a future season (Season 2) when occupancy has declined by 30, 50 or 80%. We set the probability of detection in Season 2 as the average probability of detection calculated in Season 1. Under a scenario where occupancy has declined by 30% between seasons, we also calculated the probability of detecting a significant difference in occupancy between seasons if 10, 25 or 50 sites were added to this study.

Philoria kundagungan was detected on 52 occasions at 21 sites, giving a naïve site occupancy of 0.66 across all sites. Detection probability ranged between 0.92 ± 0.06 (s.e.) and 0.712 ± 0.11 over the three survey rounds. The mean number of individuals detected per transect across all survey rounds was 3.24 ± 0.54. The top-ranked detection model (Table 1) indicated that detection decreased with an increase in temperature at the time of the survey (β = –0.867 ± 0.36) (Fig. 2). This model had 5.1 times more support than the next model. The overall top-ranked occupancy model suggested that elevation influenced site occupancy (β = 4.12 ± 1.46) (Table 2) with sites ≤750 m elevation having an estimated occupancy probability of ≤0.05 ± 0.03 and sites ≥930 m having an estimated occupancy probability of ≥0.95 ± 0.04 (Fig. 3). The two top-ranked occupancy models were included in model-averaging parameter estimates.

Table 1.  Top five models describing P. kundagungan detection probability (p)
Reported is Akaike’s Information Criterion with small sample bias correction (AICc), the relative difference in AICc values compared with the best-supported model (ΔAICc), model weight (w_i) and the number of parameters in the model (k). (.) = parameter constant across occasions. Covariates are shown in parentheses. ‘minimum temperature 5’ = the mean minimum temperature for the previous five days

Fig. 2.  The influence of air temperature (°C) at the time of survey on the probability of Philoria kundagungan detection. Shaded area indicates 95% confidence intervals.

Table 2.  Top five models examining variables associated with P. kundagungan occupancy probability (ψ)
Reported is Akaike’s Information Criterion with small sample and over-dispersion bias correction (QAICc), the relative difference in QAICc values compared with the best-supported model (ΔQAICc), model weight (w_i), and the number of parameters in the model (k). (.) = parameter constant across sites. Covariates are shown in parentheses

Fig. 3.  The influence of elevation (metres above Australian Height Datum (m AHD)) on the probability of Philoria kundagungan occupancy. Shaded area indicates 95% confidence intervals.

Across all surveyed sites elevation ranged between 631 and 1042 m, average annual precipitation ranged between 815 and 1182 mm and average annual temperature ranged between 14.7°C and 16.3°C. Vegetation was classified as rainforest for 82.0% of sites, with four sites classified as hoop pine (Araucaria cunninghamii) plantation and two sites classified as eucalypt forest. Models including vegetation type or aspect were much less parsimonious than the constant model. Rainforest spinach (Elatostema reticulatum) was found at 100% of survey sites. Philoria kundagungan were detected in all sampled vegetation types and on all aspects.

The lowest occupancy probability value of a site we found to contain P. kundagungan was 0.05 ± 0.03. The combined length of streams with a >0.05 occupancy probability was estimated to be 551.4 km, giving a total of 11.2 km² of potential habitat (Table 3). Potential habitat comprised 36 fragments spaced ≥1 km apart and 72.8% of potential habitat was located within national parks (Fig. 4).

Table 3.  Summary of modelled potential habitat for P. kundagungan, by stream length, area (using a stream buffer of 20 m) and percentage falling within national parks, within 10 classes of probability of P. kundagungan occupancy

Fig. 4.  The modelled geographic distribution of potential suitable Philoria kundagungan habitat based on calculated occupancy probabilities for streams.

Our occupancy power analysis results show that the probability of detecting a significant (at α = 0.05) difference in P. kundagungan occupancy between two seasons if occupancy declined by 30.0, 50.0 or 80.0% was 32.4, 73.6 and 99.7% respectively. Under a scenario where occupancy has declined by 30.0% between seasons, the probability of detecting a significant difference in occupancy if 10, 25 or 50 sites were added to this study was 40.6, 51.7 and 67.2% respectively.

Using the approach described, we attained baseline occupancy data for P. kundagungan and a more than 2-fold increase in the number of published localities for this species. Although this species can be challenging to detect, by conducting surveys during daylight hours under optimal environmental conditions and by broadcasting P. kundagungan call recordings, a high detection probability (0.92) can be achieved. Our monitoring approach relied upon the detection of vocalising male P. kundagungan and we did not detect females. This is a common limitation in anuran surveys (Newell et al. 2013) and there is currently no other method available that can effectively detect female P. kundagungan. We found that detection decreased with an increase in temperature at the time of the survey, which is likely to be the result of the influence of air temperature on the seasonal timing of P. kundagungan breeding activity. Similar results were reported for P. richmondensis, for which ambient temperature was found to be the primary variable driving calling activity, with calling activity peaking between 15.0°C and 16.0°C and then decreasing as temperatures warmed into summer months (Willacy et al. 2015). Temperature and recent rainfall affect the calling activity of most frog species (Navas 1996; Van Sluys et al. 2012; Plenderleith et al. 2018), and unfavourable environmental conditions may decrease calling activity, reducing detection probabilities (Pérez-Granados et al. 2020). Further research investigating P. kundagungan calling phenology would ensure surveys are timed to coincide with environmental conditions most favourable to calling activity, thus maximising detection probabilities.

Previous accounts of P. kundagungan state that the species is rare and provide limited information on habitat, for example ‘montane sub-tropical rainforest in very damp situations’ (Ingram and Corben 1975) Our study greatly increased the number of sites from which the species is known and has clarified its habitat requirements. Occupancy modelling shows that, P. kundagungan occurs at high probability at high elevation (>850 m) headwater streams while headwater streams in lowland areas (<700 m) are very rarely occupied. Our modelling found that aspect and average annual precipitation did not influence occupancy. Vegetation type, if directly adjacent to rainforest, also did not influence occupancy, and our surveys detected P. kundagungan in high elevation hoop pine (Araucaria cunninghamii) plantations in Koreelah State Forest, a habitat type previously not documented as supporting this species.

Our mapping of modelled potential habitat showed that the availability of suitable P. kundagungan habitat is extremely small (~11 km²) and highly fragmented. Although the dispersal ability of this frog is not known, the fragmented nature of populations and current distribution suggests that dispersal between sites is unlikely. Known localities of the three most northern species of Philoria reveal that these species are separated by as little as 8 km of lowland areas (Knowles et al. 2004), which suggests that lower elevation areas surrounding P. kundagungan habitats may act as dispersal barriers. Potential habitat was distributed across 36 geographically isolated (>1 km) habitat fragments, with all surveyed populations characterised by a small number of individuals (<4 calling males per 100 m). Of particular concern is a fragment of habitat (currently known to support this species) in the north of Tooloom National Park and the adjacent southern part of Beaury State Forest, which is less than 0.5 km² in size and more than 5 km from the nearest potential P. kundagungan habitat. Due to its isolation, very small size, and the resulting higher risk of local extinction, the Tooloom population should be carefully monitored to ensure the unique genetic variants found in this area (L. Bolitho, unpubl. data) are not lost.

Although most of the modelled potential habitat is located within national parks, there are several areas of habitat that are not afforded this level of protection including: all headwater streams in the high elevation areas of Koreelah State Forest, all tributaries of Crystal Creek in Beaury State Forest, all tributaries of the Condamine River located on Wilsons Peak and several headwater streams adjacent to Main Range National Park. These areas of potential habitat are modelled as very likely to contain P. kundagungan and are potentially at risk from anthropogenic activities, such as logging. Surveys in these areas should be prioritised given the high probability of occurrence and potential impact from known land use activities.

As the climate in north-eastern New South Wales and south-east Queensland continues to warm over the next half-century (Laurance et al. 2011b; Evans et al. 2014), P. kundagungan populations will likely become further isolated and exposed to novel conditions. These areas are forecast to experience 0.6–1.2° of warming (relative to the 1986–2005 average) by 2030 (Laurance et al. 2011b; Evans et al. 2014). Due to dispersal barriers, P. kundagungan populations will only be able to track suitable thermal conditions upslope, leading to a reduction in geographic distribution and further fragmentation of existing habitat. As P. kundagungan calling activity (a surrogate for breeding) appears to be linked with temperature, forecast temperature increases (Evans et al. 2014) are also likely to disrupt seasonal breeding patterns and thus affect recruitment rates.

Climate change is likely to impact P. kundagungan through progressively increasing the severity and length of drought events and fire seasons. Climate models indicate a trend towards lower annual rainfall in eastern Australia (Laurance et al. 2011b; Evans et al. 2014; Pachauri et al. 2014) and moisture inputs from cloud-stripping which amount to ~40% of annual rainfall in the region (Hutley et al. 1997) are expected to be reduced through a forecast rise in the average cloud base elevation (Oliveira et al. 2014). As average annual moisture inputs decline, P. kundagungan habitat may also be threatened by desiccation, leading to a reduction in available P. kundagungan habitat and recruitment rates. Reduced annual moisture inputs coupled with increased temperatures, will increase the risk of fire in P. kundagungan habitat, particularly in rainforest that was previously thought to have a low fire risk.

In November–December 2019, severe wildfires burnt through the study area, potentially impacting >50% of P. kundagungan habitat (DAWE 2020). Under the combined pressures of increased temperatures, reduced moisture inputs and increased fire risk, P. kundagungan populations are likely to be impacted in the coming decades. By regularly monitoring P. kundagungan occupancy rates and distribution using the approach we have identified here, management strategies to mitigate the effects of climate change on this species such as captive husbandry and assisted translocation (Newell 2018) could be implemented if required.

Our baseline dataset offers a new opportunity for future studies to detect changes through time in the occupancy and distribution of P. kundagungan. This approach has been successfully used for long-term monitoring of numerous amphibian species and has revealed essential information regarding their basic conservation parameters (MacKenzie et al. 2002; Mazerolle et al. 2007; MacKenzie et al. 2018). For example, Adams et al. (2013) found that occupancy modelling was capable of assessing population trends on a continental scale for a taxonomically varied group of amphibians. Ongoing monitoring will be particularly important for P. kundagungan as it is apparent that this species is under multiple threats and conservation efforts for this species will require careful planning underpinned by up-to-date knowledge of this species’ occupancy and distribution status across the landscape. Our occupancy power analysis results show that adding at least 10 survey sites to reiterations of this study would substantially increase the probability of detecting a significant difference in P. kundagungan occupancy between two seasons. Due to the interannual variability observed in the occupancy of amphibian species (Collins and Halliday 2005), we suggest that our surveys should be repeated annually initially to understand short-term dynamics of occupancy and then at 3-year intervals to guide long-term conservation strategies.

This study provides important baseline occupancy and distribution data for a poorly known endemic montane frog, which is a critical first step in developing an effective species monitoring program. This study revealed that, within its range, P. kundagungan more frequently occupies high-elevation headwater streams while headwater streams in areas <700 m are very rarely occupied. We find that the availability of potentially suitable habitat is extremely small and highly fragmented. Our occurrence probability map provides a practical tool to design P. kundagungan conservation actions; however, without continued monitoring, the future impact of climate change and other potential threats to P. kundagungan could go unnoticed. Reiterations of this study would yield the crucial data required to guide the improvement of conservation efforts for this highly threatened montane species.

The authors declare no conflicts of interest.

This research was undertaken in accordance with relevant state requirements. In New South Wales the study was approved by the Southern Cross University Animal Care and Ethics Committee (authority ARA-18-077) with fieldwork authorised by New South Wales Scientific Licence SL100492. In Queensland, Animal Care and Ethics approval was via Department of Agriculture and Fisheries (authority SA 2016/08/565) and fieldwork was undertaken in collaboration with Queensland Parks and Wildlife Service (supervised by HBH).