Predicting the distribution of Australian frogs and their overlap with Batrachochytrium dendrobatidis under climate change

Amphibians, with over 40% of assessed species listed as threatened, are disproportionately at risk in the global extinction crisis. Among the many factors implicated in the ongoing loss of amphibian biodiversity are climate change and the disease chytridiomycosis, caused by the fungus Batrachochytrium dendrobatidis (Bd). These threats are of particular concern in Australia, where Bd has been implicated in the declines of at least 43 frog species, and climate change is emerging as an additional threat. We explore how climate change is likely to affect the distributions of Australian frog species and Bd to the year 2100, and how the spatial and climatic niche overlap between chytridiomycosis‐declined frogs and Bd could shift.


| INTRODUC TI ON
The Earth currently faces an extinction crisis of such magnitude that it has been suggested a mass extinction event is underway (Barnosky et al., 2011;Ceballos et al., 2015Ceballos et al., , 2017. Among vertebrates, amphibians are at the highest risk of extinction, with 41% of the over 7000 currently assessed species listed as threatened (IUCN, 2021), with many more species also likely to be threatened, given the ~8500 currently described species (AmphibiaWeb, 2022).
Amphibians face threats that are common to all terrestrial vertebrates, including habitat loss and modification (Becker et al., 2007) and predation by invasive species (Kats & Ferrer, 2003). Unlike other vertebrates, however, amphibians face the additional threat of a global disease pandemic (chytridiomycosis) caused by the chytrid fungi Batrachochytrium dendrobatidis (Bd) (Berger et al., 1998) and less commonly, B. salamandrivorans (Bsal) (Martel et al., 2013). These pathogens have had a substantial impact on amphibian diversity and are associated with the decline and extinction of hundreds of species worldwide (Scheele et al., 2019).
As with other taxa, climate change is also recognized as an emerging threat to amphibians (Duan et al., 2016;Ryan et al., 2014).
Specifically, amphibian population persistence could be threatened by higher variability in seasonal precipitation (Sodhi et al., 2008;Walls et al., 2013), increased aridity driven by higher temperatures , and prolonged periods of drought (Scheele et al., 2012), which can disrupt phenological cycles and cause recruitment failure (Cohen et al., 2018;Todd et al., 2011). For many species, populations may become unviable in parts of their distribution, leading to contraction of geographic ranges and elevated risk of extinction. Predicting the threat that climate change may pose to species' future persistence is essential for effective conservation planning (Groves et al., 2012). One tool to examine species responses to climate change is species distribution modelling (SDM) (Araújo et al., 2005;Rodríguez et al., 2007;Wiens et al., 2009). The basic premise of SDM is that through comparison of the environmental variables between places a species is known to occur and presumed to not occur, one can define the current realized niche of that species in environmental space (Elith & Leathwick, 2009). These associations can then be combined with predicted future climatic conditions to predict the way a species' distribution may shift under alternative projected climate change scenarios (Araújo et al., 2005).
This information can then be used to help guide effective conservation planning through the diversion of resources to locations where a species has the greatest chance of persistence (Reside et al., 2018).
In the case of amphibians, predicting species' future distributions is complicated by the potential for climate change to also shift the distribution of Bd (Rödder et al., 2010;Xie et al., 2016). The distribution of Bd, the more common of the two pathogens that cause chytridiomycosis, is strongly climate-limited: Bd favours relatively cool and humid environments, becoming less common in drier and hotter climates (Berger et al., 2005). The spatial extent of Bd is predicted to be reduced under most climate change scenarios (Rödder et al., 2010), which could lead to favourable outcomes for many amphibian species if this reduces their geographic and environmental niche overlap with Bd (Xie et al., 2016). However, it is also possible that a common directional pressure of climate change may cause consistent shifts in the distribution of both amphibians and Bd, such that proportional geographic and niche overlap is maintained. Furthermore, if amphibian distributions are predicted to contract away from regions that are currently unfavourable for Bd, the relative degree of overlap could increase, elevating the risk for these species. Risk from increased exposure to Bd would be compounded for species that also suffer from reduced range size. Projections of changes to interactions between amphibians and Bd remain relatively unexplored, and the future extinction risk of many species could be higher or lower than that predicted from the direct effects of climate change alone.
In this paper, we examined the predicted effects of climate change, both direct and indirect via altered niche overlap with Bd, on Australia's amphibian fauna. Australia supports ~240 amphibian species, all of which belong to the order Anura (except for one introduced caudate, which we did not consider in our study) (Anstis, 2017).
The majority of these species are endemic, and 45 are currently considered threatened (Gillespie et al., 2020). Batrachochytrium dendrobatidis was first detected in Australia in 1978; its distribution rapidly expanded and now appears to be approaching its climatic limits (Murray et al., 2011). In eastern Australia, chytridiomycosis has taken a large toll, with the disease implicated in population declines of at least 43 species, including seven apparent extinctions (Scheele et al., 2017). Finally, climate change has been identified as a major potential threat to Australia's frog fauna (Gillespie et al., 2020), with increasing aridity, prolonged drought, and lower annual precipitation in certain regions of Australia likely to be drivers of future declines (Blaustein et al., 2010;Chiew et al., 2009;Nicholls, 2004).
Conservation planning and prioritization usually emphasise the protection of currently threatened species. However, a forwardlooking, proactive approach to prioritization requires consideration not only of threatened species, but also currently secure species at risk of becoming threatened (Cardillo et al., 2004(Cardillo et al., , 2006. To inform current and future conservation priorities, our study aimed to answer the following questions: (1) How may the effects of climate change upon the distributions of both Australian frog species and Bd influence the makeup of Australia's threatened frog species by the year 2100? (2) How many species are predicted to experience substantially reduced range sizes, increased geographic and niche overlap with Bd, or both? (3) Which biomes and geographic regions of Australia are predicted to experience more severe loss of amphibian species? 2 | ME THODS

| Data processing
Most frog species occurrence data for this project were acquired from the Atlas of Living Australia (ALA, 2021), the Global Biodiversity Information Facility (GBIF.org, 2021) and FrogID (Rowley & Callaghan, 2020). In addition to publicly available records from FrogID, access to non-public records for threatened species were granted by the FrogID curators under a licence. All records were checked for misnomers, incorrect spelling and changed species names, then filtered for quality such that records were only retained if they were located in Australia, were correct to at least two decimal places, had a coordinate uncertainty of less than 1000 m, and had a year associated with the record, following Gueta and Carmel (2016).
For each species, records were then cropped if they were located outside of the species' currently recognized distribution, using spatial polygons as defined by the Frogs of Australia application (Hoskin et al., 2015). To account for potential distributional inaccuracies, polygons were first buffered by either 0.5 degrees (~50 km) or the mean distance from the polygon's centroid to boundary, whichever was smaller, such that records lying immediately outside the polygon were retained, under the assumption they were still valid records.
Although the majority of Australian frogs occur only on the Australian continent, there are a number of species with ranges outside Australia (both natural and introduced). We included records from outside Australia in our models for these species, and used the methods described above for processing these records, although we used range boundary data from the IUCN (2021). Given Bd occurs globally, we also obtained a global data set for these records, and applied the same filtering procedures, although without the use of range boundaries. Occurrence data for Bd was extracted from Murray et al. (2010) and the 'Bd Maps' data set (Olson et al., 2013). All initial filtering was conducted in R (R Core Team, 2021), using the raster (Hijmans, 2021), sp (Pebesma & Bivand, 2005) and rgdal (Bivand et al., 2021) packages.
We accessed climatic data from the WorldClim 2.0 data set at a 2.5 arcminute resolution, as this was the finest resolution available for both current and future time periods (Fick & Hijmans, 2017). We chose to use annual mean temperature (BIO1), mean diurnal range (BIO2), annual precipitation (BIO12), precipitation in the wet quarter (BIO16), precipitation in the dry quarter (BIO17) and elevation as our predictor variables, due to their likely roles in determining ecological suitability for frogs at a general level, across a large geographic region including a vast range of environments (Anstis, 2017). We also retrieved raster layers for future climate predictions under two shared socio-economic pathways (SSPs), namely SSP126 (a lower emissions, "optimistic" scenario) and SSP370 (a moderate to high emissions, more likely, scenario) for the above-mentioned variables, for the periods 2021-2040, 2041-2060, 2061-2080 and 2081-2100. We downloaded these layers for all available climate models (BCC-CSM2-MR, CanESM5, CNRM-CM6-1, CNRM-ESM2-1, GFDL-ESM4, IPSL-CM6A-LR, MIROC-ES2L, MIROC6, and MRI-ESM2-0), so as to reduce the influence of climate model choice on our projections. We ensured these variables were not collinear (variance inflation factor <10) using the usdm package (Naimi et al., 2014).

| Species distribution modelling
We performed species distribution modelling using maximum entropy (Maxent) version 3.4.4 (Phillips et al., 2004(Phillips et al., , 2017, implemented in the dismo package (Hijmans et al., 2020). Maxent describes correlations between species occurrences and selected climatic and/or environmental covariates to define a species' fundamental niche, which can then be projected into geographic space (Elith et al., 2011). Prior to modelling frog distributions, we created a raster layer to address potential sampling biases, adapting the background manipulation method proposed by Fourcade et al. (2014).
Using the MASS package (Ripley et al., 2013), we created a Gaussian two-dimensional kernel density raster using the occurrences of all frogs in Australia and used this weighted raster to select background "pseudo-absence" points, whereby background points were more likely to be selected in areas with higher numbers of recorded occurrences (typically urban areas). This raster was generated at a resolution of 2.5 arcminutes. Using this bias layer, each species was then independently modelled. For species with ranges extending outside Australia, we created a separate bias layer for the records which occurred outside Australia following the same procedure, although setting the extent to the IUCN-described distribution. For these species, we then scaled the two density layers such that they were comparable (i.e. the cells with the highest density in Australia had an equal density to those with the highest density outside Australia).
Initially, for each species, occurrence records were thinned such that one record per grid cell of the bias layer was retained.
We defined the background extent for each species using the buffered polygon described above (i.e. either 0.5 degrees or the mean distance from the species' distributional centroid to its boundary, whichever was smaller), and cropped the bias layer and all environmental layers to this extent. Given we were modelling many species with vastly differing distribution extents, we chose against using a uniform number of background points, such as the 10,000 often employed in SDM studies. Rather, for each species, we selected 20% of cells in the background extent as background points, with selection weighted according to the values of the previously described bias layer. We then used the ENMeval package (Muscarella et al., 2014) to select the optimal model tuning parameters for each species.
We ran the "ENMevaluate" function using the occurrences for each species, the biased background points, and the environmental predictors, trialling regularisation multiplier (RM) values of 0.5, 1, 1.5 and 2, and feature classes (FC) of "linear," "hinge," "quadratic," and all combinations of these, under the "block" method. We chose the block method, whereby presence and background points are equally partitioned into four spatially derived blocks of equal numbers (Kass et al., 2021), to account for any potential spatial autocorrelation presented in the data, as suggested by Roberts et al. (2017). Our final model selection followed guidelines proposed by Kass et al. (2021), whereby the model with the highest area under the receiver operating characteristic curve (AUC) and lowest omission rate were selected.
We also created null models, using the "ENMnulls" function in the ENMeval package. We created a distribution of 100 models, using the optimal model parameters as suggested during the previous step. We used this distribution to determine whether our optimised model performed better than a model built with random data, and only continued analyses for models with a significantly (p < .05) higher AUC value and lower omission rate. We discarded any models with AUC <0.7, or an omission rate >0.25. These values were chosen such that poor performing models were excluded, though also so a suitable number of species for continent-wide analyses were retained. For species whose models met these thresholds, we projected the model to Australia, using the climatic predictors as previously stated (we did not project to areas outside of Australia). We also projected the models to future time periods, for each combination of climate model, time period and emissions scenario stated previously, amounting to 73 total projections for each species.
We used a similar approach to model the spatial distribution of Bd, using all occurrence records from our global data set to build the species distribution model, before projecting it to Australia. We created a spatial bias layer as before, using all Bd occurrence records to build the density raster, and completed the "ENMeval" step to define the optimal model parameters. Given the vastly larger spatial extent and number of presence records, we selected 100,000 background points to build the model. The same environmental variables were used as for frogs, at the same resolution. We projected the Bd model to potential future environments as we did for frogs.

| Projecting distributions and niche overlap with Bd
As climatic conditions change, species may disperse to new areas (Travis et al., 2013). However, given the taxonomic coverage of our study, it was unfeasible to incorporate explicit dispersal capabilities in a meaningful way into our analyses, given the differing dispersal capacities amongst frog species (Berry, 2001;Parris, 2004), landscape barriers to dispersal, and lack of knowledge of the dispersal capacity of many species. However, projections allowing for no dispersal would predispose species to patterns of decline.
To incorporate dispersal, we first converted the spatial projection of the current time period for each species to a binary presence/absence raster, such that all cells with a suitability value equal to or higher than the maxSSS threshold computed for that species were deemed "suitable," and all other cells deemed "unsuitable," following methods proposed by Liu et al. (2013). We chose max-SSS, a threshold that maximizes the sum of the true-positive rate and the true-negative rate of predicted points of occurrence, as this has been shown to be one of the most reliable and accurate measures for such conversions (Liu et al., 2013). We calculated the area of this raster to determine the present extent of occurrence (EOO) for each species. We then created a polygon cropped to encompass those cells deemed suitable-this essentially represented our modelled distribution for each species. We buffered this polygon by 0.1 degrees (approximately 10 km) and cropped the projection for the 2021-2040 period to this extent. Additionally, we only allowed frogs to "expand" into regions with the same or higher elevation (as both downwards dispersal by frogs is unlikely, and refugia of mountain-top species are likely to persist only at or above current elevation; see Reside et al., 2019), and removed any cells in the expansion zone that were at a lower elevation than any cells within the original predicted distribution. We then deemed all cells above the maxSSS threshold for the projection within this boundary to be suitable, and those outside unsuitable, before calculating the new EOO for each species at this time. We then used this projection to allow limited dispersal to the next time period, and so on, until the final temporal distribution (2081-2100). We conducted this process for each climate model and emissions scenario independently. Although this method likely underestimates the dispersal capacities for some particularly labile species, it identified regions in which species had capacity to expand; given the scope of our study, we deemed identification of which species were likely to remain stable (or expand) to be sufficient. We did not crop any of the predictions for Bd: given how Bd has expanded in Australia during the past four decades, if regions become suitable, it is likely to expand into these regions. We calculated estimated EOO for Bd using the maxSSS threshold as described previously.
Following this, we then performed a multivariate environmental similarity surface (MESS) analysis upon our predictions in order to test the analogy between current and future environments, for each species . We used the "mess" function in the dismo package to calculate the proportion of our predicted distributions for each frog species at the year 2100 that occurred in extrapolated environments. We used these results as a further selection criterion in addition to those mentioned earlier for models to be included in the final analysis and removed species with predicted highly extrapolated environments (a mean between climate models of over 50%) in either emissions scenario.
We modelled the potential environmental niche overlap between frog species and Bd using four metrics. First, we modelled geographic spatial overlap by converting each projection of Bd's distribution to a binary presence/absence raster, and overlaid each frog species' binary prediction for each time period and emissions scenario onto this raster, calculating the proportion of a frog's distribution that was shared with Bd. Because infection dynamics may differ within a species' distribution (Savage et al., 2011), we also calculated niche overlap, that is the degree of similarity between the niches of two species, between each frog species and Bd using two indices implemented in the ENMTools package (Warren et al., 2010): Schoener's D (Schoener, 1968), and "I," a variation of Hellinger's Distance (Van der Vaart, 2000), proposed by Warren et al. (2008), hereafter termed "Warren's I" for clarity. We chose to use both of these indices for comparative purposes, given their slightly differing methods of calculation. For these functions, we first cropped each frog's initial spatial projections to a polygon produced from their binary prediction, and also cropped the spatial Bd projection to this polygon. We then used the "raster.overlap" function to calculate niche overlap metrics for all species for all periods and scenarios. Finally, we calculated a niche margin index (NMI) between each applicable frog species and Bd, adapting methods proposed by Broennimann et al. (2021).
Broennimann's method was created in the context of predicting success of biological invasions by calculating an index of similarity between a species' native range and invaded range. We adapted this such that a species' "native range" was that predicted for each relevant frog species, and the "invaded range" was the predicted distribution of Bd, for each relevant scenario. We followed the script provided by Broennimann et al., available on GitHub at https://github.com/ecosp at/NMI, making use of their custom "NMI_function" R function.

| Preparation for statistical analyses
To quantify where the impacts of climate change on Australian frogs are most acute, we compared predicted changes across bioregions, as defined by the Interim Biogeographic Regionalisation for Australia (IBRA7; Australian Government Department of Agriculture, Water, and the Environment, 2012). Specifically, we considered how environmental space that is climatically suitable for frog species within each bioregion was affected, using grid cells as a quantifiable unit. To do this, we created a "majority" prediction raster for 2100 for each emissions scenario, by first overlaying each of the predictions for each climate model, and then retaining cells for which at least five of the nine models predicted presence. We then calculated the sum of grid cells occupied by each species in each bioregion at the starting (i.e. present) time period and under future scenarios. We assigned each species to a particular terrestrial biome (Olson et al., 2001) based on which biome was predicted to have the greatest number of suitable cells for that species.
We converted each species' current and predicted binary distributions into polygons and calculated the centroid of each polygon using the rgeos package (Bivand & Rundel, 2020) to investigate any spatial patterns in distributional shift under climate change. We conducted a number of linear mixed models during analysis (described as reported in the Results), employing the lmerTest package (Kuznetsova et al., 2017), and used the report package (Makowski et al., 2020) to aid the reporting of statistical results.

| Model performance
Of the ~240 frog species in Australia, we were able to produce adequate models (i.e. AUC > 0.7, OR < 0.25, p-value vs. null model <.05, and extrapolation <50%) for 141 species (mean AUC = 0.855 ± 0.004, mean omission rate = 0.120 ± 0.004, mean extrapolation = 0.095 ± 0.002, with no significant difference between climate models or emissions scenarios, as determined by linear mixed models), with frogs from each terrestrial biome of Australia represented (see Appendix S1 for all model details and species level results, and https://doi.org/10.5061/dryad.h44j0 zpn8 for plots of all included species' predicted distributions). The distributions of the species modelled broadly reflects patterns of species richness in Australian frogs, although species from the south-west and the Wet Tropics were proportionally underrepresented (Figure 1).

| Predictions for Australian frog at 2100
An ANOVA performed on a linear mixed model showed climate models to result in statistically different predictions of area, though with a small effect size (F(8) = 5.95, p < .001; Eta 2 (partial) = .02, 90% CI [8.90e −03 , 0.03]). This result highlights the importance of interrogating multiple climate models, such that results are not arbitrarily influenced F I G U R E 1 Species richness of Australian frogs represented in our study. Species richness is represented as the number of species present in each bioregion, as determined in the 7th edition of the Interim Biogeographic Regionalisation for Australia. Biomes, as defined by the Terrestrial Ecoregions of the World (Olson et al., 2001), are indicated by internal boundaries in grey. Bioregions hosting the greatest number of frog species are consistently found on the eastern seaboard, as well as Northern Australia by choosing a single climate model. Accordingly, we averaged predictions across the nine climate models for all species to negate this variation. Under both the low and the high emissions scenarios, we found substantial reductions in available climatically suitable cells for Australian frogs by the year 2100, with the higher emissions scenario resulting in more extensive declines in all regions of Australia ( Figure 2).
These predicted losses were not uniformly spread across the continent, with the tropical savannas of northern Australia retaining most available climatically suitable cells under both scenarios. The south-eastern and south-western regions of Australia experienced the most severe declines, with large proportions of climatically suitable cells predicted to be lost under the more severe SSP370 scenario in some of the most species rich bioregions of Australia (the overall net effect is somewhat masked in southern bioregions, as cells which become unsuitable for some species are offset by a southward shift of species presently residing to the north). Where we predicted the centroid of a species' distribution to shift to a large extent, this was generally in a southward direction, suggesting that for many species, cells in the northern region of their distributions became unsuitable, as cells in the south either remained or became suitable (Figure 2). Species on the eastern seaboard, particularly the south-east, were predicted to experience the greatest shifts in distribution centroid. Based on a linear mixed model, predicted centroid shifts were more extreme (β = 12.68 km, 95% CI although the same number of species are predicted to lose at least some of their climatically suitable space by 2100 under the SSP126 and SSP370 scenarios (72 species, the same in each scenario), the number of species predicted to lose more than 30% of their total extent of occurrence under the lower emissions scenario (23 species) is predicted to more than double (47 species) under the higher emissions scenario (Table 1).
Finally, we performed a chi-squared test to determine whether a species' current binary threat status (threatened/not threatened) was correlated with whether a species was predicted to expand or decline in range. For both emissions scenarios, currently threatened species were more likely to experience a range contraction, χ 2 (1, N = 141) = 7.81, p = .005. When only considering species predicted to decline in distribution, however, whether a species was listed as threatened was not a predictor of the extent to which the decline occurred, as shown by a Welch two sample t-test, for either sce-

| Batrachochytrium dendrobatidis threat at 2100
Our distribution model for Bd proved to be adequate, with an AUC of 0.70, omission rate of 0.13, and p-value (compared with a null model) <.001. We predicted areas along the entire eastern seaboard to be highly suitable for Bd, particularly in mountainous regions and the Wet Tropics (consistent with reported Bd-related declines in these regions) (Figure 4). The south-west region of Australia was also predicted to be suitable (although no declines have been reported from this region).
Under both emissions scenarios, the distribution of Bd was predicted to greatly decline by the year 2100, with reductions to the extent of occurrence of 14.3 ± 0.01% (SSP126) and 25.7 ± 0.02% (SSP370). These losses were mostly observed in the northern regions of the pathogen's distribution (excluding the Wet Tropics), and the south-west (Figure 4). We found no evidence that the distribution of Bd would shift by the year 2100 under any emissions scenario in such a way that novel contact between Bd and currently non-exposed species was likely. Given this finding, as well as the fact that many species currently in contact with Bd have not suffered declines attributable to chytridiomycosis, we restricted niche overlap analyses to species that have been reported to have experienced chytridiomycosis-related declines in the past, as defined by Scheele et al. (2019). This amounted to 24 species in our data set (Appendix S1).
Although the geographic extent of Bd is predicted to decline, the degree of niche overlap between the pathogen and species that have experienced chytridiomycosis-related declines is expected to remain relatively unchanged; that is, the predicted declines in frog distributions track with the declining distribution of Bd. This was true for all metrics of niche overlap ( Figure 5; Table S2.2), determined by performing linear mixed models, predicting each metric independently with year (current or 2100), and fitting the climate model used and the emissions scenario as random effects. All scenarios and measures showed no statistical difference by the year 2100, suggesting that the spatial interactions between Bd and chytridiomycosis-affected frog species is unlikely to be significantly altered. Although the mean for each metric centres around 0, it should be noted that some individual species, particularly when considering geographic overlap and NMI, experience relatively large changes, ranging from ~50% reductions to ~30% increases; these species, particularly those experiencing increases to these metrics, should be monitored regardless of the overall tendency for little change to niche overlap to occur.

| DISCUSS ION
Our results show that even under the most optimistic of our modelled emissions scenarios, many Australian frog species will be at risk of declines in suitable climatic space by the year 2100, with these declines being substantial for a number of species. Our results provide empirical support for the suggestion that climate change will be a significant threatening factor for amphibians over the coming decades, both in Australia (Gillespie et al., 2020) and globally (Bishop et al., 2012). While the geographic extent of Bd is predicted to contract, we found no evidence for widespread reductions in the proportional geographic or niche overlap with Bd for species that have experienced chytridiomycosis-related declines. Given there is no evidence to suggest a reduction in niche overlap, and the relative distribution size for many species afflicted by chytridiomycosis is predicted to contract, there is potential for the relative impact of Bd on some species to increase by the year 2100.
Of the 141 frog species for which we were able to fit models, 23 and 47 species are predicted to lose at least 30% of their current area of climatic suitability by 2100, under low and high emissions scenarios, respectively. However, this may not equate to an equivalent reduction in the geographic range size of these species, as some species may have a broader tolerance to climatic conditions than can be estimated from correlative SDMs, and we have not been able to account for potential phenotypic plasticity or adaptation potential in our models (Bush et al., 2016). Additionally, some species may have the capacity to disperse rapidly enough to track shifting climates (Travis et al., 2013); indeed, although we did not attempt to explicitly model potential range expansions, for a number of species across the continent, as the climates to the north of their distributions become unsuitable, climates to the south become suitable (Figure 2, Figure S2.1). We show a number of species, particularly in the Tropical Savanna in the north of Australia to have potential to increase their distributions substantially as the climate becomes more favourable to the south of their distributions ( Table 1).
On the other hand, the regions of Australia in which we predict climate-space displacement and reductions to be most severe (the east coast, south-east and south-west) are also regions where TA B L E 1 Summary of predictions for each biome relative to current modelled frog biodiversity by the year 2100. For each biome and emissions scenario, the number of frogs predicted to gain any newly available suitable climatic space (resulting in a net increase to their distribution) is shown, as well as the number of species predicted to experience a significant (>30%) net loss of suitable climate space (numbers in parentheses show the proportion of individuals this number represents, relative to the modelled species richness of that biome) Biome Species richness SSP126 SSP370 Montane 0 n/a n/a n/a n/a

F I G U R E 3 (a)
The proportion of total species for each threat status category for which adequate models were produced (and are thus included in this study). (b) The proportion of species within each threat status group predicted to experience any loss of distribution. (c) The mean proportion of remaining distribution for species predicted to experience a decline to their current predicted distributions, arranged by current threat status climate change is predicted to be more extreme, making adaptation more difficult and meaning that physiological limits are likely to be reached more quickly. Furthermore, in these regions, habitat has been heavily modified for agriculture, with remaining habitat often fragmented and limited in extent (Lindenmayer et al., 2012), likely restricting both the potential area of occupancy of species within their envelope of suitable climate and their ability to disperse to newly available environments. As such, our results suggest that we should expect a large proportion of Australia's frog species to decline in their geographic extent, and even where we predict climatic space to become suitable as other regions become unsuitable, a species' ability to take advantage of this may be hampered by both habitat fragmentation and the rate at which the climate is altered (Loarie et al., 2009). Accordingly, it is likely that for some of these species their decline will be sufficient to meet the IUCN Red List criteria for being listed as threatened.
Once a species is listed as threatened, and conservation management plans are enacted, the species' chance of recovery (or stabilization) is greatly enhanced (Hoffmann et al., 2010). However, of particular concern is the fact that in our models, under the SSP126 and SSP370 scenarios, respectively, of the species with a predicted reduction in climate space greater than 30%, 18 species (78% of total) and 36 species (79%) are currently not listed as threatened.
This means that few of these species have conservation action plans, monitoring programmes, or other targeted measures currently in place to secure their populations. There is a risk that the strong emphasis on the protection of species already threatened could leave other species vulnerable to decline. This discrepancy between current and predicted extinction risk has been termed "latent extinction risk," and it has been urged that this should become a focus for more proactive conservation planning (Cardillo et al., 2006). In the context of Australian frogs, this proactive planning could be realized by considering the impact of climate change upon all species, not only those that are currently threatened, so that conservation measures may be implemented in regions where frogs have the greatest chance of persisting under climate change. This being said, given species currently listed as threatened are more likely to be threatened by climate change in the future, and are underrepresented in our models, it is essential that the potential impacts of climate change to these species be given emphasis in conservation management plans.

F I G U R E 4
The predicted distribution of Bd and its environmental overlap with susceptible frog species by the year 2100 under two emissions scenarios (top panel), and the predicted area of Bd under each climate model to the year 2100, with black showing the mean prediction (bottom panel). In the top panel, purple shows regions were Bd was predicted to occur currently, though is predicted to become climatically unsuitable. Orange regions show where the climate is likely to remain suitable for the pathogen. The overall extent of Bd in Australia is predicted to decline by the year 2100, and to a larger extent in the higher emissions scenario. However, significant proportions of Australia, particularly around the eastern seaboard, where chytridiomycosis affects the greatest proportion of species, are likely to remain suitable

| Geographic variation in frog responses to climate change
Our results show the potential loss of future frog biodiversity to be highly uneven geographically, with the effects being felt more strongly in the biomes of southern Australia, particularly the Temperate, Montane and Mediterranean regions. In these biomes, climate change is expected to be mainly realized in the form of lower annual precipitation, reduced environmental water availability, higher temperatures and extended periods of drought (Chiew et al., 2009). These changes are likely to disrupt reproduction and drive phenological shifts, threatening population persistence (Mac Nally et al., 2017;Sheridan et al., 2018). In northern Australia, however, we predict the effects of climate change on frog distributions to generally be far less pronounced; this region is projected to see slight increases in rainfall and improved wetland connectivity (Karim et al., 2015;Whetton, 2011), which would likely favour the majority of frog  (Groves et al., 2012;Reside et al., 2018). Given the relatively lower predicted impact of climate change to frogs in Australia's north, however, different conservation priorities may be pursued.
Compared with other areas of the continent (and indeed globally), ecosystems across Australia's north remain largely intact (Woinarski et al., 2007). Additionally, the full extent of frog diversity is not yet understood in northern Australia, with studies often revealing cryptic diversity within currently described species (Catullo & Keogh, 2014). Although we predict climate change to (generally) not be a significant threat to species in this region, conservation measures remain important here; the opportunity to fully describe and maintain frog biodiversity in a way that may not be possible in south-eastern Australia should not be missed.
An important caveat of our results on the geographic distribution of predicted decreases in suitability is that many small-ranged species from the Wet Tropics were not modelled due to a lack of records and the spatial limitations of SDM (Araújo & Guisan, 2006

| Geographic and niche overlap between frogs and Bd under climate change
Although we predict the spatial extent of Bd to reduce under climate change, corroborating the results of previous studies from F I G U R E 5 Changes to various metrics of niche overlap between Bd and chytridiomycosis-affected species by the year 2100. Green displays the milder SSP126 scenario, whereas red shows the higher emissions SSP370 scenario. Black shows the mean and 95% confidence intervals for each metric in each scenario. In all cases, the mean never significantly deviates from 0, suggesting that niche overlap between Bd and afflicted frog species is unlikely to be significantly altered by climate change other regions (Rödder et al., 2010;Xie et al., 2016), our findings show that this reduction is not likely to reduce the potential impact of chytridiomycosis on susceptible species, as has been hypothesised (Rollins-Smith, 2017). Our finding that niche equivalence for chytridiomycosis-declined species and Bd generally does not decrease suggests that the relative threat posed by Bd to these species is unlikely to decline. In fact, given that the extent of occurrence of many susceptible species is predicted to contract, the proportional risk of chytridiomycosis to these frogs may increase (Becker & Zamudio, 2011). Much of the conservation management for chytridiomycosis-declined species is related to mitigating other threats to these species, particularly in refugia where Bd does not occur (Scheele et al., 2014). However, the general predicted reduction in the extent of occurrence of chytridiomycosis-declined species is likely to result in fewer chytrid-free refugia within frog distributions. Fewer refugia and reduced distributions means species may be less able to recover following mortality events associated with chytridiomycosis (Heard et al., 2015;Scheele et al., 2017;Skerratt et al., 2016).
Our analyses were unable to account for some ecological variables that will affect frog distributions over the coming decades.
Although we allowed a small amount of dispersal, there is of course potential for species to expand more than we predicted; the signatures of range expansion that we captured could result in quite substantial increases into the future as previously unsuitable environments become suitable (Travis et al., 2013). Another caveat is the assumption that the niche a species currently occupies represents its extreme physiological limits, and that adaptation to changing environments is not possible (Bush et al., 2016). Whilst the velocity of climate change is expected to outstrip the ability of many species to adapt (Huey et al., 2012), the possibility of phenotypic plasticity and climate adaptation cannot be ruled out, and in some species, selection may occur for advantageous phenotypes that reduce the overall effect of climate change upon that species (Bestion et al., 2015). However, given the magnitude of range reductions we predict for many species, it appears unlikely these responses will eliminate the threat of climate change for the majority of those likely to be significantly affected.
Although not uniformly spread across the continent, climate change is likely to drive significant declines in Australian frogs by the year 2100, particularly in the continent's south and south-eastern regions. In addition to reductions in distributions attributable to climatic change, the significant threat that Bd poses to species is unlikely to be diminished. As such, conservation management must consider the adverse effects that climate change is likely to bear upon Australian frogs, and action should be taken to secure populations in regions with the greatest chance of persistence under climate change to protect Australia's frog biodiversity over the coming decades.

ACK N OWLED G EM ENTS
B.C.S. was supported by an Australian Research Council grant (DE200100121).

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interests to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
All occurrence data and R scripts are available at https://doi. org/10.5061/dryad.h44j0 zpn8. Records for threatened species, procured from FrogID, are unable to be shared, thus have been provided at 10 km within the provided occurrence data.