Managing flows for frogs: wetland inundation extent and duration promote wetland-dependent amphibian breeding success

Study area

The Gwydir Wetlands and Macquarie Marshes are located in the northern portion of the MDB in south-eastern Australia (Fig. 1), a semiarid region characterised by low average rainfall and high evaporation rates (Quentin Grafton et al. 2022). The Gwydir Wetlands are formed on the floodplains of the lower Gwydir River catchment and the Macquarie Marshes comprise a large floodplain wetland system formed by the breakdown of the Macquarie River in the lower catchment (Fig. 1). Both wetland systems are examples of mosaic landscapes consisting of diverse vegetation types that provide habitat for dependent fauna at critical life stages, sometimes in high abundance, when they are variably inundated from river flows (Fig. 1). Recognised for their high biodiversity values, both wetlands include protected areas under the NSW reserve estate, with parts recognised as internationally important under the Ramsar Convention. Owing to water resource development, both wetlands have been substantially affected by alteration to the natural flow regime and changes in land use that has caused declines in the health and extent of wetland habitat (Kingsford 2000). To mitigate these impacts, both wetlands have a long history of environmental water management targeting the inundation of wetland habitats (Foerster 2008; Berney and Hosking 2016).

Study species

Fifteen frog taxa have been recorded across different wetland habitats in the study area. All have wide distributions and are not considered threatened. The species display a range of morphological and behavioural traits, with prior work identifying three broad groupings based on those traits that may inform their flow–ecology relationship (see Ocock et al. 2014, 2016; Sarker et al. 2022). Non-burrowing ground frogs, such as the barking marsh frog (Limnodynastes fletcheri), have limited ability to withstand drying and are reliant on floodplain habitats, including wetlands, waterholes and creeks, for refuge. They also exhibit breeding behaviour in response to inundation of wetlands by river flows in the absence of rainfall. Arboreal species and burrowing species have various adaptations that reduce cutaneous water loss, including forming a cocoon after burrowing underground (e.g. water-holding frog, Cyclorana platycephala) and adopting water-conserving postures (e.g. green tree frog, Litoria caerulea). General activity and breeding behaviour in these species is primarily in response to local rainfall and temperatures, and unrelated to flow characteristics, such as timing, extent and duration of wetland inundation, though inundated floodplain wetlands may provide refuge and breeding habitat.

Frog and tadpole surveys

Frog surveys were completed at 15 wetland sites across the Gwydir Wetlands system and at 14 sites in the Macquarie Marshes (Fig. 1) between austral spring 2015 to summer 2020. The survey sites covered a range of habitat types including vegetated mixed marsh lagoons, open waterholes and riverine habitats. All sites were surveyed each September (early spring) and November (late spring) from 2015 to 2020 with some exceptions due to site access constraints. In most years, the early spring surveys coincided with the initial arrival of flows filling the wetlands (Fig. 2 and 3). The late spring survey provided opportunity for detection of breeding success such as tadpoles or recently metamorphosed frogs and coincided with the main activity window for most frog species considered flow-dependent. Additional late summer surveys were conducted in the Gwydir Wetlands (only) to capture outcomes of managed events that inundated wetlands from late spring into late summer in February 2017 (all sites) and March 2019 (six sites only), and to capture responses to the first significant flows after a period of extended drought in March 2020 (seven sites only) (Fig. 2).

Fig. 2.

Daily rainfall (a) at Moree Airport and daily flow rates (b) from an upstream gauging station (gauge number 418004) indicative of wetland inundation extent in the Gwydir Wetlands over the study period showing timing of surveys (closed circles).

Fig. 3.

Daily rainfall (a) at Quambone Station and daily flow rates (b) from upstream gauging stations (gauge numbers 421088 and 421090 combined) indicative of wetland inundation extent in the Macquarie Marshes over the study period, showing timing of surveys (closed circles).

All field surveys were approved under NSW National Parks and Wildlife Service scientific licence SL101491 and NSW Department of Planning and Environment, Environment & Heritage Group animal research authority 141209-03.

We completed a daytime habitat assessment and tadpole survey, followed by a night-time frog survey at each site during each survey period with all surveys completed within 4 days and nights. Nocturnal transect surveys for frogs began at least 30 min after sunset. At each site, we used two approaches. We first completed a 5-min audio survey to record each frog species heard calling. Each species calling was assigned a category based on an index of calling abundance as we could not reliably discern the precise number of individuals calling particularly as the chorus size increased, with ‘rare’, 1–5 individuals calling; ‘common’, 6–10; ‘abundant’, 11–20; and ‘very abundant’, greater than 20 individuals calling (Royle 2004). We then did a visual and auditory encounter survey (VAES) equating to a total of 1-person h where two experienced field staff searched across the wetland area for a duration of 30 min. Spotlights (with a maximum of 1000 lumens) were used to search for frogs in inundated areas, along the wetland edge and into the surrounding terrestrial habitat. For each individual observed, we recorded species, method of detection (seen only, heard calling in close proximity to the surveyor, or seen and heard), and life stage (adult or recently metamorphosed). When needed, individuals were caught by hand with dampened plastic gloves and measured with a soft ruler or callipers. Hygiene protocols were followed, and all gear disinfected between survey areas.

We used two methods to distinguish metamorphs from adult frogs. First, the individual was classified as a metamorph if it had a visible residual tail stub. Second, for flow-dependent frog species only, if the individual did not have a tail stub but was equal to or smaller than the size recorded for metamorphs, namely less than ~26 mm for barking marsh frogs and salmon-striped frogs (Limnodynastes salmini), ~23 mm or less for striped marsh frogs (Limnodynastes tasmaniensis) and 12 mm or less for eastern sign-bearing froglets (Crinia parinsignifera) (following Anstis 2013 and J. F. Ocock, unpubl. data). Species taxonomy follow Anstis (2013).

Habitat assessment

A rapid assessment of habitat structure, inundation extent and water depth at each site, was completed during each daytime site survey. Three vegetation quadrats were located at each site. For each quadrat, the percentage open water, percentage terrestrial vegetation, percentage bare ground and percentage cover of each of the aquatic vegetation categories (submerged, floating, short and tall emergent) was estimated. Estimates of the vegetation categories could total more than 100% as layers of different aquatic vegetation types can overlap each other. Water depth was estimated (cm) in the centre of each quadrat. Data were then averaged for the eight habitat covariates (Table 1).

Environmental predictors

We derived a set of environmental predictors based on known or suspected influence on breeding attempts and success of flow-dependent frogs (see Table 1 for details on how each was calculated). Air temperature during the survey, rainfall (immediate and in the preceding 1, 7 and 15 days), seasonality, vegetation structure and water depth are known to influence frog calling and success (Ocock et al. 2016; Littlefair et al. 2021). We represented wetland inundation extent and inundation duration with metrics based on antecedent conditions using satellite-derived inundation maps and river flow gauge data (Thomas et al. 2015). The river gauge data were standardised for each gauge across the whole period of the study by converting all values to a relative rank (i.e. percentile, range 0–1) to allow comparison between the catchments in the global model. The data pre-processing was performed using R package dplyr (ver. 1.1.4, H. Wickham, R. François, L. Henry, K. Müller, D. Vaughan and Posit Software, see https://CRAN.R-project.org/package=dplyr). To test if breeding responses varied geographically, we included catchment (Gwydir and Macquarie) as predictors. Overall, 22 explanatory covariates were derived (Table 1).

Analyses

We first collated the VAES data to calculate the total abundance for the three guilds: non-burrowing ground frogs (‘flow-dependent’), predominantly arboreal frogs and burrowing frogs. We also collated the number of species heard calling and number of surveys where metamorphs were recorded.

We undertook detailed analyses to investigate the importance of the various hydrological and habitat covariates as drivers on breeding attempts (calling abundance index) and breeding success (number of metamorphs) in flow-dependent frog species using a random forest machine-learning algorithm (Fig. 4). The Random Forest is an ensemble learning method and details of the algorithm is provided in Breiman (2001). This analysis approach can handle both linear and non-linear relationships and the presence of multi-collinearity among the input covariates. We used the ‘ranger’ package (ver. 0.16.0, see https://cran.r-project.org/package=ranger; Wright et al. 2017) in the software program R, which allows fast implementation of the random forest.

Fig. 4.

Data analysis approach for identifying key drivers of flow-dependent frog breeding attempts (male calling) and breeding success (number of recently metamorphosed frogs).

We ran separate models for breeding attempts and breeding success as we had hypothesised suitable conditions for male frogs to call and attempt breeding may be different to factors that increased breeding success (e.g. number of individuals successfully completing larval development to the terrestrial stage). We modelled breeding attempts using a calling abundance index of the flow-dependent species as the response variable (see Table 2 for the list of species). To do this, we converted the index categories for each flow-dependent species to a number where we used the minimum number of the range the index represented (i.e. ‘rare’ was converted to 1, ‘common’ to 6, ‘abundant’ to 11 and ‘very abundant’ to 20) and then summed the converted number for all flow-dependent species to obtain one abundance total for each site survey (potential max. value = 120). To evaluate the important drivers of breeding success, we modelled total abundance of metamorphs from all flow-dependent species as the response. We combined counts from the two catchments to create a global dataset (calling: Global and metamorphs: Global). Nineteen covariates were common to both models. The breeding attempts model (20 covariates in total) included survey month (i.e. September, November or March) as we were interested if season was an important influence on calling. The breeding success models (21 covariates total) included 90- and 120-day cumulative river flow metrics as we hypothesised these longer time periods, namely lagged dataset, would be relevant to tadpoles completing metamorphosis but not to conditions for calling males.

The dataset was split into calibration (training data 80%) and validation (test data 20%) using a conditional Latin hypercube sampling (cLHC) strategy. Latin hypercube sampling is a procedure that provides an efficient way of sampling covariates from their multivariate distributions using a stratified random approach. It enables a full coverage of the range of each feature to be selected by maximally stratifying the marginal distribution (Minasny and McBratney 2006). Prior to the model fitting process, covariate selection was performed using the Boruta algorithm, a wrapper function built with the random forest algorithm (Wimalathunge and Bishop 2019). Boruta trains a random forest classifier on the extended dataset and applies an importance measure to evaluate the importance of each covariate where a higher value means the covariate (predictor) is more important (see Fig. 4).

Once the optimum covariates were identified, each model was optimised for key random forest model hyperparameters (mtry, n-nodes, variance type) through 10-fold cross-valuation techniques. Once the optimum model hyperparameters were obtained, those models were validated using leave-one-out cross-validation (LOOCV). Model quality was assessed using: (1) root mean squared error (RMSE) as a measure of model accuracy; (2) normalised RMSE; and (3) Lin’s concordance correlation (LCC), which measures the agreement between observation and predictions of the same variable (the 1:1 line indicates perfect agreement in a validation plot), where low RMSE and high LCC indicate better quality model. This is also indicated by a low value in the variation in the magnitude of errors. Additionally, we calculated mean error (ME) as a measure of model bias, Pearson’s correlation (r) and R² as most published work reports on model evaluation based on coefficient of determination (R²). The variable importance plots (VIP) were used to determine the explanatory variables with the greatest influence on each response variable. The relationship between the response variable and input features is not directly explainable in complex machine-learning models. We used Shapley Additive Explanations (SHAP), an interpretive machine-learning (IML) approach based on game theory, to identify the relationship between the final features included in the random forest model and the considered response variable. The use of SHAP values allowed us to unlock the ‘black-box’ nature machine-learning algorithm when explaining how the considered covariates influenced the two response variables. We visualised the direction of the relationships between the five most important factors and the two responses variables using the ‘fastshap’ package (ver. 0.1.0, B. M. Greenwell, see https://cran.r-project.org/package=fastshap) in the software program R.

Climatic conditions

River flows and rainfall varied considerably over the study period in both catchments (Fig. 2 and 3). The highest flows into the wetlands were in late 2016 leading to widespread inundation until early 2017, which was followed by severe drought conditions from spring 2017 until early 2020, with 2019 as the driest and warmest year on record (Bureau of Meterology 2020). Water managed for the environment was delivered to the wetlands in both catchments during the study period, though the volumes and flow release timing differed among the years depending on condition of the wetland and water availability (NSW Office of Environment and Heritage 2018, 2019, 2020; NSW Department of Planning, Industry and Environment – Environment, Energy and Science 2021, 2022).

Species abundances

Over the 6 years, 13 frog species were observed from 343 site surveys across 95 survey nights (148 and 195 site surveys in the Marshes and Gwydir, respectively) and over 11 000 individual frogs were recorded during the VAES (sum of heard-only, seen-only and heard-and-seen) (Table 2). Six species comprised 99% of all observations and were all non-burrowing species (n = 11 304). Each of these species were recorded widely across the sites and survey nights and were heard calling at all 29 sites (5-min audio survey). Recently metamorphosed frogs from five of these species were recorded at 28 sites across 81 site surveys (Table 2).

The remaining seven species were rarely recorded, comprising less than 1% of all counted frogs (Table 2). This group consisted of two arboreal and five burrowing species. Overall, these species were recorded (seen or heard) at 18 sites, on 30 survey nights and 41 site surveys, and heard only at 11 sites on 15 different survey nights, nearly all when at least 5 mm of rain had fell in the preceding 7 days. Recently metamorphosed frogs from four of these species were seen in five individual surveys.

Random forests: model section and validation

Boruta selected 12 and 11 variables as the most important for breeding attempts and breeding success respectively. Model validation suggested the breeding attempts model was slightly better than breeding success though both models performed reasonably well, indicating the predictors were suitable for predicting the observed values (Table 3, Fig. 5). The LCC was closer to 1 than 0 and the RMSE was low in both models, though both showed some bias (ME and MAE, Table 3). The RMSE, which shows how far predictions are from the true values, was relatively low for both models indicating a small deviation between residuals and observed data. The higher RMSE in the breeding success model is due to the outliers (Fig. 5). The breeding attempts model also had lower variation in the magnitude of errors (difference between RMSE and ME).

Fig. 5.

Validation plot of the measured v. predicted counts for (a) breeding attempts (calling abundance index) and (b) breeding success measured by number of recently metamorphosed frogs.

Influence of flow and habitat on breeding activity and breeding success

The covariate representing inundation extent was the most important influence on the calling abundance index response variable, followed by inundation timing, total volume of flows in the preceding 30 days and water depth (Fig. 6a). The seven other variables included in the important factors had comparatively low influence on breeding responses (Fig. 6a). Inundated area at the time of the survey positively corresponded with calling activity, with the SHAP calling response value increasing linearly 10-fold as wetland inundation increased from 0.5 to 1.0 (Fig. 7). For the covariates relating to timing, calling was higher in September compared to November, and at the start of a flood event compared to other stages (Fig. 7). Increasing water depths up to ~15 cm corresponded with increasing calling, but for further increases in water depth, calling activity decreased. Flow volumes exceeding 0.7 of peak flows (standardised for the catchment over the survey period) had a strongly linear positive influence on calling response (Fig. 7).

Fig. 6.

The variable importance plots generated from the best performing random forest model developed from the data generated on (a) breeding activity and (b) breeding success showing the top selected variables (see Table 2 for full explanation of the variable abbreviations).

Fig. 7.

SHAP dependence plots demonstrating the influence of the top five important variables on (a) breeding activity and (b) breeding success.

All four covariates representing total flows in the preceding months and timing of the flows were an important influence on the total number of metamorphs (Fig. 6b). The relationships for three of the four flow totals (60, 90 and 120 days preceding the surveys) were positive, where the SHAP value for the number of metamorphs increased sharply, from 1 to 20, when flows increased above 0.75 of the peak standardised flow (Fig. 7). In contrast, flows above 0.8 of standardised flows in the preceding 30 days corresponded with decreasing numbers of metamorphs. The timing of flow events (FloodTiming) was the third most important variable, with the highest numbers of recently metamorphosed frogs associated with the middle stage of a flow event. Although catchment was included in the final set of important factors, there was no clear difference in frog responses between the Macquarie and Gwydir systems.