Global exposure risk of frogs to increasing environmental dryness

Supplementary Information

This document contains the R workflow for processing and analysing the raw data, and generating figures for the manuscript titled “Global exposure risk of frogs to increasing environmental dryness”.

GitHub Repository: The Rmarkdown file can be downloaded from the Code drop down menu (top right).

Abbreviations

• AI: Aridity Index
• \(e_a\): actual vapour pressure
• \(e_s\): saturation vapour pressure
• EWL: Evaporative water loss
• IUCN: International Union for Conservation of Nature
• PDSI: Palmer Drought Severity Index
• PRISMA: Preferred Reporting Items for Systematic Reviews and Meta-Analyses
• RH: Relative humidity
• \(r_i\): Relative skin resistance or resistance to water loss
• SA: surface-area
• SPP: Shared Socioeconomic Pathway
• VPD: Vapour pressure deficit
• WoS: Web of Science
• WU: water uptake

---

Spatial risk

Descriptors

The following series of code is to estimate the risk of anuran amphibians to increasing aridity (Aridity Index; AI) and extreme drought (Palmer Drought Severity Index; PDSI). The AI was categorised to five categories (Humid, Dry sub-humid, Semi-arid, Arid, and Hyper-arid) and the PDSI was categorised to seven categories (Extremely moist, Very moist, Moderate moist, Normal, Moderate drought, Severe drought, Extreme drought) based on descriptions from Budyko (1961) and Palmer (1965), respectively (Table S1).

Generalised ecological types or ‘ecotype’ of each anuran species were classified based on descriptions from Moen and Wiens (2017) focusing on adult behaviour and microhabitat preferences outside the breeding season, given that many anurans breed in water but are not adapted to live in water all year (Table S2).

Table S1

Table S1 Aridity index (Budyko, 1961) and the Palmer Drought Severity Index (Palmer, 1965) categorised.

Aridity index (AI)	AI value	Palmer Drought Severity Index (PDSI)	PDSI Value
Humid	≥ 0.65	Extremely moist	≥ 4
Dry sub-humid	0.50 – ≤ 0.65	Very moist	3 – ≤ 4
Semi-arid	0.20 – ≤ 0.50	Moderate moist	2 – ≤ 3
Arid	0.05 – ≤ 0.20	Normal	-2 – ≤ 2
Hyper-arid	< 0.05	Moderate drought	-3 – ≤ -2
		Severe drought	-4 – ≤ -3
		Extreme drought	< -4

Table S2

Table S2 Microhabitat preference or ecotype (Moen and Wiens, 2017).

Ecotype	Description
Aquatic	Almost always in water.
Arboreal	Ability to climb vegetation and spends a substantial amount of time above ground (e.g. forests, on top of vegetations in wetlands).
Semi-aquatic	Stay near water bodies, usually permanent water (e.g. ground-level wetlands). Often observed in water and nearby vegetation.
Ground-dwelling	Terrestrial only, not always near water bodies, however always in a moist environment (e.g. forest leaf litter, moist soils, scrublands, grassland). Sometimes above ground and sometimes semi-fossorial. Can be classified as broad.
Fossorial	The non-breeding season is spent underground in burrows. Capable of burrowing and aestivating.
Stream-dwelling	Restricted to near fast-flowing streams usually on rocks or vegetation near streams.

Species richness

Anuran (frogs and toads) distribution shape files were obtained from the International Union for Conservation of Nature’s Red List of Threatened Species (IUCN Red List; extracted on 11/01/2021). Species richness or species assemblages was defined as the sum of species in each grid cell (0.5°), based on the geographic range, and was calculated using the rasterizeIUCN and calcSR function from the from the rasterSp package. Ecotype-specific species richness was also calculated.

filedir <- "rasterSp/"

# Clean raw data
frog_dat <- read.csv(file.path(data_path, "anuran_species_list.csv")) %>%
    dplyr::select(genus:strategy, SVL_cm:binomial_tree_phylo) %>%
    dplyr::mutate(lnSVL = log(SVL_cm), lnMass = log(mass_g), lnArea = log(shape_area),
        IUCN = factor(IUCN, levels = c("Least Concern", "Near Threatened", "Vulnerable",
            "Endangered", "Critically Endangered", "Extinct")), order_name = factor(order_name),
        family_name = factor(family_name), ecotype = factor(ecotype)) %>%
    dplyr::filter(order_name == "ANURA" & ecotype != "" & IUCN != "" & IUCN != "Extinct") %>%
    droplevels()

# Convert shape files into rasters and save to file once (downloaded shape
# files 08/01/2020) only extracted distribution of native range anuran_rast <-
# rasterizeIUCN(dsn = paste0(filedir, 'ANURA/ANURA.shp'), resolution = 0.5,
# seasonal = c(1, 2), origin = 1, presence = c(1,2), save = TRUE, path =
# paste0(rasterSp_path))

# Calculate anuran richness and by ecotype The calcSR function uses a stepwise
# procedure to calculate the sum of species for each grid cell.
anuran_sr <- rasterSp::calcSR(species_names = frog_dat$binomial_IUCN, path = paste0(rasterSp_path))
aquatic_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Aquatic",
    "binomial_IUCN"], path = paste0(rasterSp_path))
arboreal_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Arboreal",
    "binomial_IUCN"], path = paste0(rasterSp_path))
fossorial_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Fossorial",
    "binomial_IUCN"], path = paste0(rasterSp_path))
ground_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Ground-dwelling",
    "binomial_IUCN"], path = paste0(rasterSp_path))
semi_aq_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Semi-aquatic",
    "binomial_IUCN"], path = paste0(rasterSp_path))
stream_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Stream-dwelling",
    "binomial_IUCN"], path = paste0(rasterSp_path))

# Convert raster to matrix then to data frame
anuran_sr_df <- raster::as.data.frame(raster::rasterToPoints(anuran_sr)) %>%
    dplyr::rename(species_n = layer)
aquatic_sr_df <- raster::as.data.frame(raster::rasterToPoints(aquatic_sr)) %>%
    dplyr::rename(aquatic_n = layer)
arboreal_sr_df <- raster::as.data.frame(raster::rasterToPoints(arboreal_sr)) %>%
    dplyr::rename(arboreal_n = layer)
fossorial_sr_df <- raster::as.data.frame(raster::rasterToPoints(fossorial_sr)) %>%
    dplyr::rename(fossorial_n = layer)
ground_sr_df <- raster::as.data.frame(raster::rasterToPoints(ground_sr)) %>%
    dplyr::rename(ground_n = layer)
semi_aq_sr_df <- raster::as.data.frame(raster::rasterToPoints(semi_aq_sr)) %>%
    dplyr::rename(semi_aq_n = layer)
stream_sr_df <- raster::as.data.frame(raster::rasterToPoints(stream_sr)) %>%
    dplyr::rename(stream_n = layer)

There were 5636 anuran species from the IUCN Red List available for the AI and PDSI analysis.

Calculate AI and PDSI

High resolution (~4 km²) global dataset on precipitation (mm/month) and potential evapotranspiration (mm/month) were obtained from Abatzoglou et al. (2018) under 1) the current climate (1970–2000), 2) an intermediate greenhouse gas emission scenario of +2°C (Shared Socioeconomic Pathway 2–4.5; SPP2–4.5), and 3) a high greenhouse gas emission or “business-as-usual” scenario of +4°C (SSP5–8.5) by 2080–2099 (Fig. S1a–f). AI was calculated as the ratio of the total precipitation relative to the potential evapotranspiration under the current climate, SPP2–4.5, and SSP5–8.5 scenario (Fig. S1g–i).

We used a self-calibrated PDSI with Penman–Monteith potential evapotranspiration representing the current climate (1970–2000) and an intermediate and high emission scenario by 2080–2099 (SSP2–4.5 and SSP5-8.5) from Zhao and Dai (2022). The SSP2–4.5 and SSP5–8.5 scenarios were based on the average of 25 CMIP6 models of precipitation, evapotranspiration, soil moisture, and runoff (Eyring et al., 2016), where the mean annual surface temperature is expected to increase by 2.7°C (2.1–3.5°C range) and 4.4°C (3.3–5.7°C range), respectively by 2080–2100 (IPCC, 2021).

# Import the downloaded files
ppt_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate19812010_ppt.nc")),
    anuran_sr)  # Precipitation year 1981-2010 - Pr (mm)
pet_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate19812010_pet.nc")),
    anuran_sr)  # Evapotranspiration 1981-2010 - ET0 (mm)

# +2C and +4C future scenarios
ppt_2C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate2C_ppt.nc")),
    anuran_sr)
pet_2C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate2C_pet.nc")),
    anuran_sr)
ppt_4C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate4C_ppt.nc")),
    anuran_sr)
pet_4C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate4C_pet.nc")),
    anuran_sr)

# Obtain mean values from 1981-2010
ppt_mean_rast <- raster::calc(ppt_rast, fun = mean, na.rm = TRUE)
pet_mean_rast <- raster::calc(pet_rast, fun = mean, na.rm = TRUE)
ppt_2C_mean_rast <- raster::calc(ppt_2C_rast, fun = mean, na.rm = TRUE)
pet_2C_mean_rast <- raster::calc(pet_2C_rast, fun = mean, na.rm = TRUE)
ppt_4C_mean_rast <- raster::calc(ppt_4C_rast, fun = mean, na.rm = TRUE)
pet_4C_mean_rast <- raster::calc(pet_4C_rast, fun = mean, na.rm = TRUE)

# Obtain mean values from 10-year monthly PDSI
PDSI_ssp245_rast <- raster::stack(x = file.path(spatial_path, "pdsisc.monthly.1900-2100.r2.5x2.5.EnsAvg25Models.TP2.ipe-2.ssp245.nc"))
PDSI_ssp585_rast <- raster::stack(x = file.path(spatial_path, "pdsisc.monthly.1900-2100.r2.5x2.5.EnsAvg25Models.TP2.ipe-2.ssp585.nc"))

PDSI_cur_rast <- raster::subset(PDSI_ssp245_rast, grep("X197.*|X198.*|X199.*", names(PDSI_ssp245_rast),
    value = T))
PDSI_2C_rast <- raster::subset(PDSI_ssp245_rast, grep("X207.*|X208.*|X209.*", names(PDSI_ssp245_rast),
    value = T))
PDSI_4C_rast <- raster::subset(PDSI_ssp585_rast, grep("X207.*|X208.*|X209.*", names(PDSI_ssp585_rast),
    value = T))

PDSI_cur_mean_rast <- raster::calc(PDSI_cur_rast, fun = mean, na.rm = TRUE)
PDSI_2C_mean_rast <- raster::calc(PDSI_2C_rast, fun = mean, na.rm = TRUE)
PDSI_4C_mean_rast <- raster::calc(PDSI_4C_rast, fun = mean, na.rm = TRUE)

# Calculate aridity index [Precipitation (ppt) / Evapotranspiration (pet)]
ai_rast <- raster::overlay(x = ppt_mean_rast, y = pet_mean_rast, fun = function(x,
    y) {
    return(x/y)
})
ai_2C_rast <- raster::overlay(x = ppt_2C_mean_rast, y = pet_2C_mean_rast, fun = function(x,
    y) {
    return(x/y)
})
ai_4C_rast <- raster::overlay(x = ppt_4C_mean_rast, y = pet_4C_mean_rast, fun = function(x,
    y) {
    return(x/y)
})

# Calculate change in PDSI Reproject intensity raster to match anuran_sr
PDSI_cur_mean_rast <- raster::projectRaster(PDSI_cur_mean_rast, anuran_sr)
PDSI_2C_mean_rast <- raster::projectRaster(PDSI_2C_mean_rast, anuran_sr)
PDSI_4C_mean_rast <- raster::projectRaster(PDSI_4C_mean_rast, anuran_sr)

PDSI_2C_diff_rast <- raster::overlay(x = PDSI_cur_mean_rast, y = PDSI_2C_mean_rast,
    fun = function(x, y) {
        return(y - x)
    })
PDSI_4C_diff_rast <- raster::overlay(x = PDSI_cur_mean_rast, y = PDSI_4C_mean_rast,
    fun = function(x, y) {
        return(y - x)
    })

# Convert raster to matrix then to data frame
ai_df <- raster::as.data.frame(raster::rasterToPoints(ai_rast))
ai_2C_df <- raster::as.data.frame(raster::rasterToPoints(ai_2C_rast))
ai_4C_df <- raster::as.data.frame(raster::rasterToPoints(ai_4C_rast))

PDSI_2C_diff_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_2C_diff_rast))
PDSI_4C_diff_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_4C_diff_rast))

# Convert aridity index into category
ai_df <- ai_df %>%
    # Recode
dplyr::mutate(category = case_when(is.infinite(layer) ~ "Humid", layer >= 0.65 ~
    "Humid", layer >= 0.5 & layer < 0.65 ~ "Dry sub-humid", layer >= 0.2 & layer <
    0.5 ~ "Semi-arid", layer >= 0.05 & layer < 0.2 ~ "Arid", layer < 0.05 ~ "Hyper-arid")) %>%
    # Convert to ordered factor
dplyr::mutate(category = factor(category, levels = c("Hyper-arid", "Arid", "Semi-arid",
    "Dry sub-humid", "Humid"), ordered = TRUE))

ai_2C_df <- ai_2C_df %>%
    # Recode
dplyr::mutate(category = case_when(is.infinite(layer) ~ "Humid", layer >= 0.65 ~
    "Humid", layer >= 0.5 & layer < 0.65 ~ "Dry sub-humid", layer >= 0.2 & layer <
    0.5 ~ "Semi-arid", layer >= 0.05 & layer < 0.2 ~ "Arid", layer < 0.05 ~ "Hyper-arid")) %>%
    # Convert to ordered factor
dplyr::mutate(category = factor(category, levels = c("Hyper-arid", "Arid", "Semi-arid",
    "Dry sub-humid", "Humid"), ordered = TRUE))

ai_4C_df <- ai_4C_df %>%
    # Recode
dplyr::mutate(category = case_when(is.infinite(layer) ~ "Humid", layer >= 0.65 ~
    "Humid", layer >= 0.5 & layer < 0.65 ~ "Dry sub-humid", layer >= 0.2 & layer <
    0.5 ~ "Semi-arid", layer >= 0.05 & layer < 0.2 ~ "Arid", layer < 0.05 ~ "Hyper-arid")) %>%
    # Convert to ordered factor
dplyr::mutate(category = factor(category, levels = c("Hyper-arid", "Arid", "Semi-arid",
    "Dry sub-humid", "Humid"), ordered = TRUE))

# Convert PDSI to category
PDSI_2C_diff_df <- PDSI_2C_diff_df %>%
    # Recode
dplyr::mutate(change_2C = case_when(layer >= 4 ~ "4", layer >= 3 & layer < 4 ~ "3",
    layer >= 2 & layer < 3 ~ "2", layer >= 1 & layer < 2 ~ "1", layer >= -1 & layer <
        1 ~ "0", layer >= -2 & layer < -1 ~ "-1", layer >= -3 & layer < -2 ~ "-2",
    layer >= -4 & layer < -3 ~ "-3", layer < -4 ~ "-4")) %>%
    # Convert to ordered factor
dplyr::mutate(change_2C = factor(change_2C, levels = c("-4", "-3", "-2", "-1", "0",
    "1", "2", "3", "4"), ordered = TRUE))

PDSI_4C_diff_df <- PDSI_4C_diff_df %>%
    # Recode
dplyr::mutate(change_4C = case_when(layer >= 4 ~ "4", layer >= 3 & layer < 4 ~ "3",
    layer >= 2 & layer < 3 ~ "2", layer >= 1 & layer < 2 ~ "1", layer >= -1 & layer <
        1 ~ "0", layer >= -2 & layer < -1 ~ "-1", layer >= -3 & layer < -2 ~ "-2",
    layer >= -4 & layer < -3 ~ "-3", layer < -4 ~ "-4")) %>%
    # Convert to ordered factor
dplyr::mutate(change_4C = factor(change_4C, levels = c("-4", "-3", "-2", "-1", "0",
    "1", "2", "3", "4"), ordered = TRUE))

# Crop and mask
world      <- sf::read_sf(file.path(spatial_path, "ne_50m_land/ne_50m_land.shp"))
world_sp   <- sf::as_Spatial(world)
world_spdf <- sp::SpatialPolygonsDataFrame(world_sp, world_sp@data) 

# Precipitation
ppt_curr_plot <- raster::as.data.frame(raster::rasterToPoints(ppt_mean_rast)) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = layer)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_continuous_sequential(palette = "YlGnBu", lim = c(0, 650), name = "Precipitation (mm)") +
  ggtitle("Precipitation (1981-2010)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10, barheight = 0.3, 
                                label.position = "bottom")) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10),
        legend.title = element_text(size = 8), 
        legend.position = "bottom") +
  coord_fixed(ratio = 1)  

ppt_2C_plot <- raster::as.data.frame(raster::rasterToPoints(ppt_2C_mean_rast)) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = layer)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_continuous_sequential(palette = "YlGnBu", lim = c(0, 650)) +
  ggtitle("Precipitation (+2°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)

ppt_4C_plot <- raster::as.data.frame(raster::rasterToPoints(ppt_4C_mean_rast)) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = layer)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_continuous_sequential(palette = "YlGnBu", lim = c(0, 650)) +
  ggtitle("Precipitation (+4°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)  

ppt_legend <- cowplot::get_legend(ppt_curr_plot)

ppt_prow <- cowplot::plot_grid(ppt_curr_plot + theme(legend.position = "none"), 
                   ppt_2C_plot +  theme(legend.position = "none"), 
                   ppt_4C_plot + theme(legend.position = "none"), 
                   align = "v", ncol = 3, labels = c('a', 'b', 'c')) 

ppt_plots <- cowplot::plot_grid(ppt_prow, ppt_legend, ncol = 1, rel_heights = c(1, .1))

# Evapotranspiration
pet_curr_plot <- raster::as.data.frame(raster::rasterToPoints(pet_mean_rast)) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = layer)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_continuous_sequential(palette = "Rocket", lim = c(0, 300), name = "Evapotranspiration (mm)") +
  ggtitle("Evapotranspiration (1981-2010)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  guides(fill = guide_colourbar(barwidth = 10, barheight = 0.3, 
                                label.position = "bottom")) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10),
        legend.title = element_text(size = 8), 
        legend.position = "bottom") +
  coord_fixed(ratio = 1)  

pet_2C_plot <- raster::as.data.frame(raster::rasterToPoints(pet_2C_mean_rast)) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = layer)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_continuous_sequential(palette = "Rocket", lim = c(0, 300)) +
  ggtitle("Evapotranspiration (+2°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)

pet_4C_plot <- raster::as.data.frame(raster::rasterToPoints(pet_4C_mean_rast)) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = layer)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_continuous_sequential(palette = "Rocket", lim = c(0, 300)) +
  ggtitle("Evapotranspiration (+4°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)  

pet_legend <- cowplot::get_legend(pet_curr_plot)

pet_prow <- cowplot::plot_grid(pet_curr_plot + theme(legend.position = "none"), 
                   pet_2C_plot +  theme(legend.position = "none"), 
                   pet_4C_plot + theme(legend.position = "none"), 
                   align = "v", ncol = 3, labels = c('d', 'e', 'f')) 

pet_plots <- cowplot::plot_grid(pet_prow, pet_legend, ncol = 1, rel_heights = c(1, .1))

# Aridity Index
arid_col <- c('#8E063B', '#CB6D53', '#E99A2C', '#F5D579', 'white')
ai_curr_plot <- ggplot() +
  geom_raster(data = ai_df, aes(y = y, x = x, fill = category)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = arid_col,
                    guide = guide_legend(reverse = TRUE),
                    name = "Aridity Index") +
  ggtitle("Aridity Index (1981-2010)") +
  ylab(NULL) + xlab(NULL) +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + 
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10),
        legend.title = element_text(size = 8), 
        legend.key.size = unit(0.2, 'cm'),
        legend.position = "bottom") +
  coord_fixed(ratio = 1)

ai_2C_plot <- ggplot() +
  geom_raster(data = ai_2C_df, aes(y = y, x = x, fill = category)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = arid_col,
                    guide = guide_legend(reverse = TRUE),
                    name = "Aridity Index") +
  ggtitle("Aridity Index (+2°C 2080-2100)") +
  ylab(NULL) + xlab(NULL) +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + 
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)

ai_4C_plot <- ggplot() +
  geom_raster(data = ai_4C_df, aes(y = y, x = x, fill = category)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = arid_col,
                    guide = guide_legend(reverse = TRUE),
                    name = "Aridity Index") +
   ggtitle("Aridity Index (+4°C 2080-2100)") +
  ylab(NULL) + xlab(NULL) +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme_void() + 
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)

ai_legend <- cowplot::get_legend(ai_curr_plot)

ai_prow <- cowplot::plot_grid(ai_curr_plot + theme(legend.position = "none"), 
                   ai_2C_plot +  theme(legend.position = "none"), 
                   ai_4C_plot + theme(legend.position = "none"), 
                   align = "v", ncol = 3, labels = c('g', 'h', 'i')) 

ai_plots <- cowplot::plot_grid(ai_prow, ai_legend, ncol = 1, rel_heights = c(1, .1))

cowplot::plot_grid(ppt_plots, pet_plots, ai_plots, ncol = 1)

Fig. S1. Climate data used to calculate the aridity index. Mean precipitation (mm) under (a) the current scenario from 1981-2010, (b) an intermediate emission scenario of +2°C (Shared Socioeconomic Pathways 2 - 4.5; SSP2-4.5), and (c) a high emission scenario of +4°C (SSP5-8.5) by 2080-2100. Mean potential evapotranspiration (mm) under (d) the current scenario from 1981-2010, (e) an intermediate emission scenario of +2°C (SSP2-4.5), and (f) a high emission scenario of +4°C (SSP5-8.5) by 2080-2100. Calculated Aridity Index for (g) the current scenario from 1981-2010, (h) an intermediate emission scenario of +2°C (SSP2-4.5), and (i) a high emission scenario of +4°C (SSP5-8.5) by 2080-2100.

AI risk

To examine the relationship of species richness with aridity, the number of species per gird cell was overlapped with the aridity raster layer, where each grid was assigned an AI category. The change in species richness between the current and projected (either +2 or +4 °C warming) AI category was calculated as the change in AI category grids (resolution 0.5° for decimal degree coordinates) occupied by anurans relative to the future projection. A decrease indicates reduced number of species with the assigned AI category and vice versa.

# Merge aridity index, +4C and amphibian species richness
ai_sp_df <- ai_df %>% 
  dplyr::full_join(ai_2C_df, by = c("x" ,"y")) %>%
  dplyr::full_join(ai_4C_df, by = c("x" ,"y")) %>%
  dplyr::full_join(anuran_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(aquatic_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(arboreal_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(fossorial_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(ground_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(semi_aq_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(stream_sr_df, by = c("x" ,"y")) %>% 
  dplyr::rename(aridity = layer.x,
                aridity_2C = layer.y,
                aridity_4C = layer,
                category = category.x,
                category_2C = category.y,
                category_4C = category) %>%
  drop_na(category)

# Calculate grid cells occupied for current climate
ai_sp_sum <- data.frame(ai_sp_df %>% 
                          dplyr::group_by(category) %>% 
                          dplyr::summarise(species_n   = length(species_n[!is.na(species_n)]),
                                           aquatic_n   = length(aquatic_n[!is.na(aquatic_n)]),
                                           arboreal_n  = length(arboreal_n[!is.na(arboreal_n)]),
                                           fossorial_n = length(fossorial_n[!is.na(fossorial_n)]),
                                           ground_n    = length(ground_n[!is.na(ground_n)]),
                                           semi_aq_n   = length(semi_aq_n[!is.na(semi_aq_n)]),
                                           stream_n    = length(stream_n[!is.na(stream_n)]))) %>%
  dplyr::mutate(scenario = "Current",
                all_freq       = species_n / sum(species_n) * 100,
                aquatic_freq   = aquatic_n / sum(aquatic_n) * 100,
                arboreal_freq  = arboreal_n / sum(arboreal_n) * 100,
                fossorial_freq = fossorial_n / sum(fossorial_n) * 100,
                ground_freq    = ground_n / sum(ground_n) * 100,
                semi_aq_freq   = semi_aq_n / sum(semi_aq_n) * 100,
                stream_freq    = stream_n / sum(stream_n) * 100)

# Calculate grid cells occupied for future climate
ai_2C_sp_sum <- data.frame(ai_sp_df %>% 
                             dplyr::group_by(category_2C) %>% 
                             dplyr::summarise(species_n   = length(species_n[!is.na(species_n)]),
                                              aquatic_n   = length(aquatic_n[!is.na(aquatic_n)]),
                                              arboreal_n  = length(arboreal_n[!is.na(arboreal_n)]),
                                              fossorial_n = length(fossorial_n[!is.na(fossorial_n)]),
                                              ground_n    = length(ground_n[!is.na(ground_n)]),
                                              semi_aq_n   = length(semi_aq_n[!is.na(semi_aq_n)]),
                                              stream_n    = length(stream_n[!is.na(stream_n)]))) %>%
  dplyr::mutate(scenario = "SSP245",
                all_freq       = species_n / sum(species_n) * 100,
                aquatic_freq   = aquatic_n / sum(aquatic_n) * 100,
                arboreal_freq  = arboreal_n / sum(arboreal_n) * 100,
                fossorial_freq = fossorial_n / sum(fossorial_n) * 100,
                ground_freq    = ground_n / sum(ground_n) * 100,
                semi_aq_freq   = semi_aq_n / sum(semi_aq_n) * 100,
                stream_freq    = stream_n / sum(stream_n) * 100) %>% 
  dplyr::rename(category = category_2C) 

ai_4C_sp_sum <- data.frame(ai_sp_df %>% 
                             dplyr::group_by(category_4C) %>% 
                             dplyr::summarise(species_n   = length(species_n[!is.na(species_n)]),
                                              aquatic_n   = length(aquatic_n[!is.na(aquatic_n)]),
                                              arboreal_n  = length(arboreal_n[!is.na(arboreal_n)]),
                                              fossorial_n = length(fossorial_n[!is.na(fossorial_n)]),
                                              ground_n    = length(ground_n[!is.na(ground_n)]),
                                              semi_aq_n   = length(semi_aq_n[!is.na(semi_aq_n)]),
                                              stream_n    = length(stream_n[!is.na(stream_n)]))) %>%
  dplyr::mutate(scenario = "SSP585",
                all_freq       = species_n / sum(species_n) * 100,
                aquatic_freq   = aquatic_n / sum(aquatic_n) * 100,
                rboreal_freq   = arboreal_n / sum(arboreal_n) * 100,
                fossorial_freq = fossorial_n / sum(fossorial_n) * 100,
                ground_freq    = ground_n / sum(ground_n) * 100,
                semi_aq_freq   = semi_aq_n / sum(semi_aq_n) * 100,
                stream_freq    = stream_n / sum(stream_n) * 100) %>% 
  dplyr::rename(category = category_4C)

PDSI risk

With a monthly prediction of PDSI from 1950 to 2100 globally available from Zhao and Dai (2022), we categorised future drought risk in three ways: 1) an increase in drought intensity (decrease in PDSI; Fig. S2a–c), 2) an increase in drought frequency (monthly PDSI counts below -2 per year; Fig. S2d–f), and 3) an increase in drought duration (number of consecutive months with PDSI values below -2; Fig. S2g–i). Change in drought intensity (ΔPDSI_[intensity]), frequency (ΔPDSI_[frequency]), and duration (ΔPDSI_[duration]) under a +2 or +4 °C warming scenario (2080–2100) was calculated relative to the 1970–2000 monthly climatology per gird cell (ΔPDSI = PDSI_[future] – PDSI_[current]). Note, the temporal resolution was limited to monthly variation as the PDSI was developed to monitor long-term meteorological drought.

# Change in PDSI intensity 
PDSI_sp_df <- anuran_sr_df %>% 
  dplyr::full_join(PDSI_2C_diff_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(PDSI_4C_diff_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(aquatic_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(arboreal_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(fossorial_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(ground_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(semi_aq_sr_df, by = c("x" ,"y")) %>% 
  dplyr::full_join(stream_sr_df, by = c("x" ,"y")) %>% 
  dplyr::rename(PDSI_2C = layer.x,
                PDSI_4C = layer.y) %>%
  drop_na(change_2C) %>%
  filter(species_n != "NA")

# Calculate grid cells occupied for current climate
PDSI_sp_2C <- data.frame(PDSI_sp_df %>% 
                           dplyr::group_by(change_2C) %>% 
                           dplyr::summarise(species_n   = length(species_n[!is.na(species_n)]),
                                            aquatic_n   = length(aquatic_n[!is.na(aquatic_n)]),
                                            arboreal_n  = length(arboreal_n[!is.na(arboreal_n)]),
                                            fossorial_n = length(fossorial_n[!is.na(fossorial_n)]),
                                            ground_n    = length(ground_n[!is.na(ground_n)]),
                                            semi_aq_n   = length(semi_aq_n[!is.na(semi_aq_n)]),
                                            stream_n    = length(stream_n[!is.na(stream_n)]))) %>%
  dplyr::mutate(all_freq       = species_n / sum(species_n) * 100,
                aquatic_freq   = aquatic_n / sum(aquatic_n) * 100,
                arboreal_freq  = arboreal_n / sum(arboreal_n) * 100,
                fossorial_freq = fossorial_n / sum(fossorial_n) * 100,
                ground_freq    = ground_n / sum(ground_n) * 100,
                semi_aq_freq   = semi_aq_n / sum(semi_aq_n) * 100,
                stream_freq    = stream_n / sum(stream_n) * 100)

PDSI_sp_4C <- data.frame(PDSI_sp_df %>% 
                           dplyr::group_by(change_4C) %>% 
                           dplyr::summarise(species_n   = length(species_n[!is.na(species_n)]),
                                            aquatic_n   = length(aquatic_n[!is.na(aquatic_n)]),
                                            arboreal_n  = length(arboreal_n[!is.na(arboreal_n)]),
                                            fossorial_n = length(fossorial_n[!is.na(fossorial_n)]),
                                            ground_n    = length(ground_n[!is.na(ground_n)]),
                                            semi_aq_n   = length(semi_aq_n[!is.na(semi_aq_n)]),
                                            stream_n    = length(stream_n[!is.na(stream_n)]))) %>%
  dplyr::mutate(all_freq       = species_n / sum(species_n) * 100,
                aquatic_freq   = aquatic_n / sum(aquatic_n) * 100,
                arboreal_freq  = arboreal_n / sum(arboreal_n) * 100,
                fossorial_freq = fossorial_n / sum(fossorial_n) * 100,
                ground_freq    = ground_n / sum(ground_n) * 100,
                semi_aq_freq   = semi_aq_n / sum(semi_aq_n) * 100,
                stream_freq    = stream_n / sum(stream_n) * 100)

# Calculate frequency of moderate to extreme drought (<-2 PDSI) per year.

# Temporal change in PDSI intensity
PDSI_2C_time_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_ssp245_rast))
PDSI_4C_time_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_ssp585_rast))

# Convert wide to long
freq_2C_df <- PDSI_2C_time_df %>%
  tidyr::pivot_longer(!c("x","y"), names_to = "dates", values_to = "PDSI") %>%
  data.frame() %>%
  tidyr::separate(dates, c("year", 'months')) %>%
  dplyr::mutate(year = as.numeric(substring(year, 2))) %>%
  dplyr::filter(year >= 1950 & year != 2100)

freq_4C_df <- PDSI_4C_time_df %>%
  tidyr::pivot_longer(!c("x","y"), names_to = "dates", values_to = "PDSI") %>%
  data.frame() %>%
  tidyr::separate(dates, c("year", 'months')) %>%
  dplyr::mutate(year = as.numeric(substring(year, 2))) %>%
  dplyr::filter(year >= 1950 & year != 2100)

# Summarise average monthly counts <-2 PDSI from 1970 to 1999 or 2080 to 2100
freq_mean_df <- freq_2C_df %>%
  dplyr::filter(year >= 1970 & year <= 1999) %>%
  dplyr::group_by(x, y, year) %>%
  dplyr::summarise(count = sum(PDSI < -2)) %>%
  dplyr::group_by(x, y) %>%
  dplyr::summarise(mean_curr = mean(count))
  
freq_2C_diff_df <- freq_2C_df %>%
  dplyr::filter(year >= 2080 & year != 2100) %>%
  dplyr::group_by(x, y, year) %>%
  dplyr::summarise(count = sum(PDSI < -2)) %>% 
  dplyr::group_by(x, y) %>%
  dplyr::summarise(mean_2C = mean(count)) %>% 
  dplyr::inner_join(freq_mean_df, by = c("x" = "x", "y" = "y")) %>%
  dplyr::mutate(diff_2C = mean_2C - mean_curr)

freq_4C_diff_df <- freq_4C_df %>%
  dplyr::filter(year >= 2080 & year != 2100) %>%
  dplyr::group_by(x, y, year) %>%
  dplyr::summarise(count = sum(PDSI < -2)) %>% 
  dplyr::group_by(x, y) %>%
  dplyr::summarise(mean_4C = mean(count)) %>% 
  dplyr::inner_join(freq_2C_diff_df, by = c("x" = "x", "y" = "y")) %>%
  dplyr::mutate(diff_4C = mean_4C - mean_curr)

# Calculate duration of moderate to extreme drought (<-2 PDSI) from 1970 to 1999
dur_curr_df <- freq_2C_df %>%
  dplyr::mutate(bin = ifelse(PDSI < -2, 1, 0)) %>%
  dplyr::filter(year >= 1970 & year <= 1999) %>%
  dplyr::group_by(x, y) %>% 
  dplyr::summarise(mean_curr = mean(rle(bin)$lengths[rle(bin)$values==1])) %>%
  replace(is.na(.), 0)

# Calculate duration of moderate to extreme drought (<-2 PDSI) from 2080 to 2100
dur_2C_df <- freq_2C_df %>%
  dplyr::mutate(bin = ifelse(PDSI < -2, 1, 0)) %>% 
  dplyr::filter(year >= 2080 & year != 2100) %>%
  dplyr::group_by(x, y) %>% 
  dplyr::summarise(mean_2C = mean(rle(bin)$lengths[rle(bin)$values==1])) %>%
  replace(is.na(.), 0) %>%
  dplyr::inner_join(dur_curr_df, by = c("x" = "x", "y" = "y")) %>%
  dplyr::mutate(diff_2C = mean_2C - mean_curr)

dur_4C_df <- freq_4C_df %>% 
  dplyr::mutate(bin = ifelse(PDSI < -2, 1, 0)) %>%
  dplyr::group_by(x, y) %>% 
  dplyr::filter(year >= 2080 & year != 2100) %>%
  dplyr::summarise(mean_4C = mean(rle(bin)$lengths[rle(bin)$values==1])) %>%
  replace(is.na(.), 0) %>%
  dplyr::inner_join(dur_2C_df, by = c("x" = "x", "y" = "y")) %>%
  dplyr::mutate(diff_4C = mean_4C - mean_curr)

The simultaneous risk of increasing drought intensity, frequency, and duration within a grid cell that are occupied by anurans (species assemblages) was calculated by converting each risk category as binary. Grid cells with a ΔPDSI_[intensity] below -2 (indicating increased drought intensity relative to current scenario) were assigned a ‘1’ binary. Both ΔPDSI_[frequency] and ΔPDSI_[duration] were assigned a binary of ‘1’ if the grid cell has a value of 1 month or higher (indicating increase in frequency or duration relative to current scenario). The number of overlapping binaries were summed up per grid cell. Therefore, a risk factor of 2 indicate species assemblages in the grid cell are at increasing risk of two drought events. We estimated which species assemblages were at risk of experiencing drought events using an arbitrary risk factor scale (species richness × drought risk), where grid cells with high drought risk and high species richness have higher “assemblage-level risk” than grid cells with high drought risk and low species richness (low assemblage-level risk).

# Intensity
names(PDSI_2C_diff_rast) <- "delta_int_2C"
names(PDSI_4C_diff_rast) <- "delta_int_4C"

# Frequency
freq_diff_rast <- rasterFromXYZ(freq_4C_diff_df) # convert to raster
crs(freq_diff_rast) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
freq_diff_rast <- projectRaster(freq_diff_rast, anuran_sr) # match anuran extent
freq_diff_df   <- raster::as.data.frame(raster::rasterToPoints(freq_diff_rast))

PDSI_freq_diff <- subset(freq_diff_rast, 4:5)
names(PDSI_freq_diff) <- c("delta_freq_2C", "delta_freq_4C")

# Duration
dur_diff_rast <- raster::rasterFromXYZ(dur_4C_df) # convert to raster
crs(dur_diff_rast) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
dur_diff_rast <- raster::projectRaster(dur_diff_rast, anuran_sr) # match anuran extent
dur_diff_df <- raster::as.data.frame(raster::rasterToPoints(dur_diff_rast))

PDSI_dur_diff  <- subset(dur_diff_rast, 4:5)
names(PDSI_dur_diff) <- c("delta_dur_2C", "delta_dur_4C")

# Combine relative PDSI metrics
PDSI_risk_comb_rast <- raster::stack(PDSI_2C_diff_rast, PDSI_4C_diff_rast, PDSI_freq_diff, PDSI_dur_diff, resample(anuran_sr, PDSI_4C_diff_rast))
                             
PDSI_risk_comb_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_risk_comb_rast)) %>%
  dplyr::rename("sp_n" = layer) %>%
  dplyr::filter(!is.na(sp_n)) %>%
  dplyr::mutate(delta_int_2C_bin = ifelse(delta_int_2C < -1, 1, 0), # delta_PDSI < -1
         delta_int_4C_bin = ifelse(delta_int_4C < -1, 1, 0), # delta_PDSI < -1
         delta_freq_2C_bin = ifelse(delta_freq_2C >1, 1, 0), # month > 1
         delta_freq_4C_bin = ifelse(delta_freq_4C >1, 1, 0), # month > 1
         delta_dur_2C_bin = ifelse(delta_dur_2C >1, 1, 0), # month > 1
         delta_dur_4C_bin = ifelse(delta_dur_4C >1, 1, 0)) %>% # month > 1 
  dplyr::rowwise() %>%
  dplyr::mutate(count_2C = sum(delta_int_2C_bin, delta_freq_2C_bin, delta_dur_2C_bin, na.rm = T),
         
         count_4C = sum(delta_int_4C_bin, delta_freq_4C_bin, delta_dur_4C_bin, na.rm = T),
         risk_2C  = sp_n * count_2C,
         risk_4C  = sp_n * count_4C,
         count_2C = factor(count_2C, levels = c("3", "2", "1")),
         count_4C = factor(count_4C, levels = c("3", "2", "1"))
         )

# Combine absolute PDSI metrics
PDSI_ab_rast <- raster::stack(PDSI_cur_mean_rast, PDSI_2C_mean_rast, PDSI_4C_mean_rast, freq_diff_rast, dur_diff_rast, resample(anuran_sr, PDSI_4C_mean_rast)) 
PDSI_ab_rast_crop <- raster::mask(crop(PDSI_ab_rast, extent(world)), world) # crop 

PDSI_ab_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_ab_rast_crop)) %>%
  dplyr::rename(PDSI_cur = "layer.1",
                PDSI_2C  = "layer.2",
                PDSI_4C  = "layer.3",
                freq_cur = "mean_curr.1",
                freq_2C  = "mean_2C.1",
                freq_4C  = "mean_4C.1",
                dur_cur  = "mean_curr.2",
                dur_2C   = "mean_2C.2",
                dur_4C   = "mean_4C.2",
                sp_n     = "layer.4") %>%
  filter(sp_n != "NA")

colours_PDSI <- RColorBrewer::brewer.pal(9, "RdBu")

# Intensity
PDSI_cur_plot <- PDSI_ab_df %>%
  dplyr::mutate(PDSI_cur_cat = case_when(
    PDSI_cur >= 2 & PDSI_cur < 3 ~ '2',  
    PDSI_cur >= 1 & PDSI_cur < 2 ~ '1',  
    PDSI_cur >= -1 & PDSI_cur < 1 ~ '0'
  )) %>% 
  dplyr::mutate(PDSI_cur_cat = factor(PDSI_cur_cat,
                                levels = c('0', '1', '2'),
                                ordered = TRUE)) %>%
  filter(PDSI_cur_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = PDSI_cur_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = c("#F7F7F7", "#D1E5F0" ,"#92C5DE")) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean PDSI intensity (1970–2000)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

PDSI_2C_plot <- PDSI_ab_df %>%
  dplyr::mutate(PDSI_2C_cat = case_when(
    PDSI_2C >= 4 ~ "4", 
    PDSI_2C >= 3 & PDSI_2C < 4 ~ '3',
    PDSI_2C >= 2 & PDSI_2C < 3 ~ '2',  
    PDSI_2C >= 1 & PDSI_2C < 2 ~ '1',  
    PDSI_2C >= -1 & PDSI_2C < 1 ~ '0',
    PDSI_2C >= -2 & PDSI_2C < -1 ~ '-1',
    PDSI_2C >= -3 & PDSI_2C < -2 ~ '-2',
    PDSI_2C >= -4 & PDSI_2C < -3 ~ '-3',   
    PDSI_2C < -4 ~ '-4' 
  )) %>% 
  dplyr::mutate(PDSI_2C_cat = factor(PDSI_2C_cat,
                                levels = c('-4', '-3', '-2', '-1', '0', '1', '2', '3', '4'),
                                ordered = TRUE)) %>%
  filter(PDSI_2C_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = PDSI_2C_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = colours_PDSI) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean PDSI intensity (+2°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

PDSI_4C_plot <- PDSI_ab_df %>%
  dplyr::mutate(PDSI_4C_cat = case_when(
    PDSI_4C >= 4 ~ "4", 
    PDSI_4C >= 3 & PDSI_4C < 4 ~ '3',
    PDSI_4C >= 2 & PDSI_4C < 3 ~ '2',  
    PDSI_4C >= 1 & PDSI_4C < 2 ~ '1',  
    PDSI_4C >= -1 & PDSI_4C < 1 ~ '0',
    PDSI_4C >= -2 & PDSI_4C < -1 ~ '-1',
    PDSI_4C >= -3 & PDSI_4C < -2 ~ '-2',
    PDSI_4C >= -4 & PDSI_4C < -3 ~ '-3',   
    PDSI_4C < -4 ~ '-4' 
  )) %>% 
  dplyr::mutate(PDSI_4C_cat = factor(PDSI_4C_cat,
                                levels = c('-4', '-3', '-2', '-1', '0', '1', '2', '3', '4'),
                                ordered = TRUE)) %>%
  filter(PDSI_4C_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = PDSI_4C_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = colours_PDSI, name = "PDSI", guide = guide_legend(reverse = TRUE, nrow = 1)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean PDSI intensity (+4°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10),
        legend.title = element_text(size = 8), 
        legend.key.height= unit(0.2, 'cm'),
        legend.key.width= unit(0.2, 'cm'),
        legend.position = "bottom") +
  coord_fixed(ratio = 1) # fixed ratio

int_legend <- cowplot::get_legend(PDSI_4C_plot)

int_prow <- cowplot::plot_grid(PDSI_cur_plot + theme(legend.position = "none"), 
                   PDSI_2C_plot +  theme(legend.position = "none"), 
                   PDSI_4C_plot + theme(legend.position = "none"), 
                   align = "v", ncol = 3, labels = c('a', 'b', 'c')) 

int_plots <- cowplot::plot_grid(int_prow, int_legend, ncol = 1, rel_heights = c(1, .1))

# Frequency
freq_cur_plot <- PDSI_ab_df %>%
  dplyr::mutate(freq_cur_cat = case_when(
    freq_cur > 0.1 & freq_cur < 1 ~ '<1'
  )) %>%
  filter(freq_cur_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = freq_cur_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = "#ededed") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean drought frequency (1970–2000)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

freq_2C_plot <- PDSI_ab_df %>%
  dplyr::mutate(freq_2C_cat = case_when(
    freq_2C >= 10 ~ "10-12", 
    freq_2C >= 8 & freq_2C < 10 ~ '8-10',
    freq_2C >= 6 & freq_2C < 8 ~ '6-8',  
    freq_2C >= 4 & freq_2C < 6 ~ '4-6',  
    freq_2C >= 2 & freq_2C < 4 ~ '2-4',
    freq_2C >= 1 & freq_2C < 2 ~ '1-2',
    freq_2C > 0.1 & freq_2C < 1 ~ '<1'
  )) %>% 
  dplyr::mutate(freq_2C_cat = factor(freq_2C_cat,
                                     levels = c('10-12', '8-10', '6-8', '4-6', '2-4', '1-2', '<1'),
                                     ordered = TRUE)) %>%
  filter(freq_2C_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = freq_2C_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#FAE9DF", "#F7F7F7")) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean drought frequency (+2°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

freq_4C_plot <- PDSI_ab_df %>%
  dplyr::mutate(freq_4C_cat = case_when(
    freq_4C >= 10 ~ "10-12", 
    freq_4C >= 8 & freq_4C < 10 ~ '8-10',
    freq_4C >= 6 & freq_4C < 8 ~ '6-8',  
    freq_4C >= 4 & freq_4C < 6 ~ '4-6',  
    freq_4C >= 2 & freq_4C < 4 ~ '2-4',
    freq_4C >= 1 & freq_4C < 2 ~ '1-2',
    freq_4C > 0.1 & freq_4C < 1 ~ '<1'
  )) %>% 
  dplyr::mutate(freq_4C_cat = factor(freq_4C_cat,
                                     levels = c('10-12', '8-10', '6-8', '4-6', '2-4', '1-2', '<1'),
                                     ordered = TRUE)) %>%
  filter(freq_4C_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = freq_4C_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#FAE9DF", "#F7F7F7"), name = "Months", guide = guide_legend(reverse = TRUE, nrow = 1)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean drought frequency (+4°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10),
        legend.title = element_text(size = 8), 
        legend.key.height= unit(0.2, 'cm'),
        legend.key.width= unit(0.2, 'cm'),
        legend.position = "bottom") +
  coord_fixed(ratio = 1) # fixed ratio

freq_legend <- cowplot::get_legend(freq_4C_plot)

freq_prow <- cowplot::plot_grid(freq_cur_plot + theme(legend.position = "none"), 
                   freq_2C_plot +  theme(legend.position = "none"), 
                   freq_4C_plot + theme(legend.position = "none"), 
                   align = "v", ncol = 3, labels = c('d', 'e', 'f')) 

freq_plots <- cowplot::plot_grid(freq_prow, freq_legend, ncol = 1, rel_heights = c(1, .1))


# Duration
dur_cur_plot <- PDSI_ab_df %>%
  dplyr::mutate(dur_cur_cat = case_when(
    dur_cur >= 2 & dur_cur < 4 ~ '2-4',
     dur_cur >= 1 & dur_cur < 2 ~ '1-2',
    dur_cur > 0.1 & dur_cur < 1 ~ '<1'
  )) %>% 
  dplyr::mutate(dur_cur_cat = factor(dur_cur_cat,
                                levels = c('2-4', '1-2', '<1'),
                                ordered = TRUE)) %>%
  filter(dur_cur_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = dur_cur_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = c("#FDDBC7", "#FAE9DF", "#F7F7F7")) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean drought duration (1970–2000)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

dur_2C_plot <- PDSI_ab_df %>%
  dplyr::mutate(dur_2C_cat = case_when(
    dur_2C >= 10 ~ ">10", 
    dur_2C >= 8 & dur_2C < 10 ~ '8-10',
    dur_2C >= 6 & dur_2C < 8 ~ '6-8',  
    dur_2C >= 4 & dur_2C < 6 ~ '4-6',  
    dur_2C >= 2 & dur_2C < 4 ~ '2-4',
    dur_2C >= 1 & dur_2C < 2 ~ '1-2',
    dur_2C > 0.1 & dur_2C < 1 ~ '<1'
  )) %>% 
  dplyr::mutate(dur_2C_cat = factor(dur_2C_cat,
                                levels = c('>10', '8-10', '6-8', '4-6', '2-4', '1-2', '<1'),
                                ordered = TRUE))%>%
  filter(dur_2C_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = dur_2C_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#FAE9DF", "#F7F7F7")) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean drought duration (+2°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

dur_4C_plot <- PDSI_ab_df %>%
  dplyr::mutate(dur_4C_cat = case_when(
    dur_4C >= 10 ~ ">10", 
    dur_4C >= 8 & dur_4C < 10 ~ '8-10',
    dur_4C >= 6 & dur_4C < 8 ~ '6-8',  
    dur_4C >= 4 & dur_4C < 6 ~ '4-6',  
    dur_4C >= 2 & dur_4C < 4 ~ '2-4',
    dur_4C >= 1 & dur_4C < 4 ~ '1-2',
    dur_4C > 0.1 & dur_4C < 1 ~ '<1'
  )) %>% 
  dplyr::mutate(dur_4C_cat = factor(dur_4C_cat,
                                levels = c('>10', '8-10', '6-8', '4-6', '2-4', '1-2', '<1'),
                                ordered = TRUE)) %>%
  filter(dur_4C_cat != "NA") %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = dur_4C_cat)) +
  geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#FAE9DF", "#F7F7F7"), name = "Months", guide = guide_legend(reverse = TRUE, nrow = 1)) +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Mean drought duration (+4°C 2080-2100)") +
  scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10),
        legend.title = element_text(size = 8), 
        legend.key.size = unit(0.2, 'cm'),
        legend.position = "bottom") +
  coord_fixed(ratio = 1) # fixed ratio

dur_legend <- cowplot::get_legend(dur_4C_plot)

dur_prow <- cowplot::plot_grid(dur_cur_plot + theme(legend.position = "none"), 
                   dur_2C_plot +  theme(legend.position = "none"), 
                   dur_4C_plot + theme(legend.position = "none"), 
                   align = "v", ncol = 3, labels = c('g', 'h', 'i')) 

dur_plots <- cowplot::plot_grid(dur_prow, dur_legend, ncol = 1, rel_heights = c(1, .1))


cowplot::plot_grid(int_plots, freq_plots, dur_plots, ncol = 1)

Fig. S2. Mean monthly Palmer Drought Severity Index (PDSI) self-calibrated with Penman–Monteith potential evapotranspiration from (a) 1970–2000, and from 2080-2100 under (b) under an intermediate emission scenario of +2°C (Shared Socioeconomic Pathways 2 - 4.5; SSP2-4.5), and (c) under a high emission scenario of +4°C (SSP5-8.5). Mean yearly frequency of moderate to extreme drought (PDSI < -2) from (d) 1970–2000, and from 2080-2100 (e) under an intermediate emission scenario of +2°C (SSP2-4.5), and (f) under a high emission scenario of +4°C (SSP5-8.5). Mean duration of moderate to extreme drought (PDSI < -2) from (g) 1970–2000, and from 2080-2100 (h) under an intermediate emission scenario of +2°C (SSP2-4.5), and (i) under a high emission scenario of +4°C (SSP5-8.5).

Summary

Ecotype AI distribution

Mean ± s.d. aridity index occupied by different anuran ecotypes.

Ecotype	AI (mean ± s.d.)
Stream-dwelling	1.46 ± 0.62
Arboreal	1.01 ± 0.62
Semi-aquatic	1.01 ± 0.58
Aquatic	0.98 ± 0.61
Ground-dwelling	0.89 ± 0.64
Fossorial	0.80 ± 0.61

---

2°C AI risk

Percent change in aridity by 2080–2100 (+ 2°C) relative to current climate (1981–2010) for all anurans and also by ecotype.

Ecotype	All (%)	Aquatic 2°C (%)	Arboreal 2°C (%)	Fossorial 2°C(%)	Ground-dwelling 2°C (%)	Semi-aquatic 2°C (%)	Stream-dwelling 2°C (%)
Hyper-arid	-2.3	13.1	-15.2	-37.25	-2.0	4.7	-18.2
Arid	5.0	10.2	6.3	9.90	4.9	10.0	26.8
Semi-arid	1.4	2.1	1.5	0.75	1.6	1.1	13.9
Dry sub-humid	2.5	11.1	9.5	-3.67	8.4	2.8	10.4
Humid	-1.4	-3.2	-2.2	-2.00	-2.4	-1.3	-1.3

---

4°C AI risk

Percent change in aridity by 2080–2100 (+ 4°C) relative to current climate (1981–2010) for all anurans and also by ecotype.

Ecotype	All (%)	Aquatic 2°C (%)	Arboreal 2°C (%)	Fossorial 2°C(%)	Ground-dwelling 2°C (%)	Semi-aquatic 2°C (%)	Stream-dwelling 2°C (%)
Hyper-arid	3.0	28.3	-4.3	-13.7	3.1	11.6	-9.1
Arid	13.2	33.0	20.6	23.1	13.5	27.0	43.7
Semi-arid	8.4	11.6	14.3	4.5	8.8	10.8	53.8
Dry sub-humid	9.0	25.2	16.6	-3.5	16.0	11.6	23.9
Humid	-5.8	-9.8	-7.6	-7.5	-7.6	-5.7	-3.8

---

2°C PDSI risk

Percent change in drought intensity (ΔPDSI), drought frequency (ΔFrequency), and drought duration (ΔDuration) by 2080–2099 (+ 2°C scenario) relative to current climate (1970–2000) in grid cells occupied by anurans and also by ecotype.

ΔPDSI	All (%)	Aquatic (%)	Arboreal (%)	Fossorial (%)	Ground-dwelling (%)	Semi-aquatic (%)	Stream-dwelling (%)
-4	0.01	0.00	0.00	0.00	0.00	0.01	0.00
-3	0.15	0.00	0.00	0.00	0.06	0.19	0.00
-2	2.00	4.22	3.35	4.35	2.18	2.46	5.48
-1	18.55	33.52	28.74	27.06	20.99	21.74	26.88
0	71.05	54.90	55.68	54.07	68.60	66.72	60.07
1	6.40	5.19	8.86	9.94	6.05	6.84	6.19
2	1.59	1.83	2.90	3.99	1.84	1.74	1.23
3	0.25	0.34	0.48	0.59	0.29	0.30	0.14

ΔFrequency (months)	All (%)	Aquatic (%)	Arboreal (%)	Fossorial (%)	Ground-dwelling (%)	Semi-aquatic (%)	Stream-dwelling (%)
10-12	0.14	0.00	0.00	0.00	0.05	0.17	0.00
8-10	0.22	0.03	0.06	0.03	0.18	0.25	0.20
6-8	0.61	0.78	0.70	0.79	0.60	0.70	0.86
4-6	2.15	3.86	3.40	3.01	2.33	2.26	3.71
2-4	12.25	21.14	19.23	16.82	13.42	13.75	26.57
1-2	16.90	24.28	17.56	18.39	18.51	18.26	10.49
0-1	56.26	42.56	46.21	46.64	53.17	54.34	48.46
-0-1	11.45	7.36	12.84	14.32	11.73	10.27	9.72

ΔDuration (months)	All (%)	Aquatic (%)	Arboreal (%)	Fossorial (%)	Ground-dwelling (%)	Semi-aquatic (%)	Stream-dwelling (%)
>10	0.27	0.00	0.00	0.00	0.10	0.34	0.00
8-10	0.06	0.00	0.00	0.00	0.04	0.06	0.01
6-8	0.11	0.00	0.00	0.00	0.09	0.11	0.07
4-6	0.47	0.67	0.61	0.75	0.48	0.55	0.57
2-4	7.79	14.23	11.94	10.29	8.61	9.21	15.78
1-2	38.11	34.96	29.11	23.85	38.14	39.87	27.17
0-1	39.57	41.03	42.90	47.79	38.34	37.31	44.53
-0-1	12.98	8.46	14.49	16.04	13.47	11.90	11.82
-1-2	0.64	0.65	0.95	1.28	0.71	0.64	0.04

---

4°C PDSI risk

Percent change in drought intensity (ΔPDSI), frequency (ΔFrequency), and duration (ΔDuration) by 2080–2099 (+ 4°C scenario) relative to current climate (1970–2000) in grid cells occupied by anurans and by ecotype.

ΔPDSI	All (%)	Aquatic (%)	Arboreal (%)	Fossorial (%)	Ground-dwelling (%)	Semi-aquatic (%)	Stream-dwelling (%)
-4	0.83	1.3	0.96	1.3	0.74	1.0	1.9
-3	4.44	9.9	7.80	9.6	5.04	5.4	15.2
-2	10.88	17.7	16.15	11.0	12.26	12.6	14.4
-1	22.33	26.6	19.04	25.3	24.78	23.6	12.2
0	44.62	31.9	37.41	32.8	42.33	40.2	41.6
1	9.63	4.5	6.79	5.6	6.97	9.5	6.8
2	2.94	2.2	4.17	4.3	2.87	2.9	3.4
3	1.39	1.9	2.30	2.4	1.60	1.4	1.9
4	2.94	3.9	5.38	7.7	3.41	3.4	2.5

ΔFrequency (months)	All (%)	Aquatic (%)	Arboreal (%)	Fossorial (%)	Ground-dwelling (%)	Semi-aquatic (%)	Stream-dwelling (%)
10-12	1.6	2.2	1.9	2.0	1.6	1.9	2.9
8-10	4.9	9.3	7.7	7.3	5.4	5.3	11.6
6-8	7.0	10.9	10.0	8.4	7.8	7.5	11.9
4-6	14.0	17.1	13.2	9.9	15.7	15.4	10.7
2-4	22.8	22.8	17.9	25.1	24.1	22.9	12.4
1-2	12.1	11.2	12.9	12.1	11.4	11.4	16.2
0-1	30.1	20.3	26.5	22.8	25.6	28.4	29.4
-0-1	7.6	6.2	9.8	12.3	8.4	7.2	4.8

ΔDuration (months)	All (%)	Aquatic (%)	Arboreal (%)	Fossorial (%)	Ground-dwelling (%)	Semi-aquatic (%)	Stream-dwelling (%)
>10	2.53	2.8	2.6	2.6	2.46	2.84	3.90
8-10	1.02	2.0	1.6	1.6	1.14	1.14	2.16
6-8	2.33	4.0	3.7	3.1	2.60	2.50	4.47
4-6	11.07	16.6	13.5	9.0	12.48	12.55	14.75
2-4	24.82	23.0	19.5	18.4	26.38	26.19	17.99
1-2	23.42	22.9	20.5	23.1	20.52	22.95	19.89
0-1	25.18	21.7	26.0	26.7	23.70	22.79	29.51
-0-1	8.83	5.8	11.0	13.5	9.79	8.12	7.27
-1-2	0.81	1.2	1.4	1.9	0.93	0.93	0.05

---

Systematic search

Exclusion criteria

• Studies with languages that were not translatable on Google Translate.
• Data from experimental treatments such as pathogen infection, pollutants or pesticide contaminants, and pharmacological manipulation (e.g., antidiuretic hormones).
• Data presented as percentage change in body weight without the initial weight.
• Data present with no time component. E.g. Change in body weight until 20% of water loss was reached (% change in mass).
• Studies where units were not presented.

Preferred Reporting Items for Systematic Reviews and Meta-Analyses

The following Boolean search string for Web of Science (WoS): (frog OR toad OR anur*) AND (water loss OR water uptake OR evaporate* OR skin resist* OR cutaneous evapo*) AND (dehydration OR desiccation OR water OR soil moist*) NOT (tadpole OR larva* OR mice OR rat OR plant, OR seed OR fish) and for Scopus: (TITLE-ABS-KEY(frog OR toad OR anur*) AND TITLE-ABS-KEY(water AND loss OR water AND uptake OR evaporate* OR skin AND resist* OR cutaneous AND evapo*)) were performed on the 21 July 2022, resulting in 449 records from 1938 to 2022 from the WoS, and 51 records from 1951 to 2022 from Scopus. We excluded topics related to chemistry, physical, materials science, medical, geology, fisheries, oceanography, and engineering. Title and abstract screening of the systematic search was conducted in Rayyan (Ouzzani et al., 2016). Forward and backward searches were also conducted on Google Scholar until the 15th August 2022 (Fig. S3) for additional papers, and we also searched for papers from Feder and Burggren (1992), Lillywhite (2006), and Hillman et al. (2009). Lastly, additional studies (n = 3) were provided by the Stellenbosch University librarian (searched 1st February 2023).

Fig. S3. PRISMA flow diagram for the systematic data-collection process. n = number of papers remaining after each stage of selection. k is the number of observations/records, and n in the last row is number of species after processing data.

Definitions

We note, that “water-proof” is a subjective term, and the mechanism is only well understood for the Phyllomedusa genus with the presence of waxy secretion of cutaneous lipids, and the stereotypical wiping behaviour (Blaylock et al., 1976; Gomez et al., 2006). Other arboreal species show dry mucus film on the skin or the presence of cutaneous surface fluid that may contribute to high skin resistance (Kobelt and Linsenmair, 1986). We labelled frogs as “water-proof” based on each study’s definition, i.e., if the authors describe the species as “water-proof”. However, not all “water-proof” frogs excrete water-proof secretions throughout the year, and may only show this adaptation during certain periods such as the dry season (R. P. Bovo, pers. comm.).

Another post-hoc definition for classifying water-proof is whether the skin resistance is above a certain threshold. Hillman et al. (2009) defined anurans as moderately water-proof with skin resistance between 2 to 20 s cm^-1, and highly water-proof with skin resistance >20 s cm^-1. However, even within studies, the same species will show different skin resistance depending on the environmental exposure. e.g. lower skin resistance at higher temperature (Buttemer and Thomas, 2003). We included both definitions in the raw data file, but we used the prior definition for analysis.

EWL measurements labelled with “cocoon” is straight-forward. Cocoon is the presence of accumulated epidermal layers in aestivating anurans. EWL measurements labelled with “hollow” were conducted with frogs in an artificial hollow shelter protecting the frog from the exposed air.

Anurans can also change EWL across the measurement period (i.e. decrease in EWL as they become more dehydrated). However, most studies either took an average value over the experiment duration, or presented the first hour of measurements. When studies provided change in EWL at different time points, we took the first hour EWL measurements. Note, almost all animals were considered fully hydrated during the start of the EWL experiment, where they were submersed or provided water prior to the EWL experiments.

Unpublished data

Additional EWL, \(r_i\), and WU data were obtained from unpublished sources. The Atelopus carrikeri data is available on an online repository (Solano & Albert, 2000), and the Thoropa miliaris data can be obtained from C. Navas (navas@usp.br) upon request.

Briefly, the method for Atelopus carrikeri is described below:

15 Atelopus carrikeri were collected near the basin of Pico Cristòbal Colòn (10°54’05.8”N 73°55’00.7”W) and brought back to the Universidad del Magdalena campus. To estimate evaporative water loss, fully hydrated individuals were briefly dried with absorbent paper and the urine was removed by lightly pressing their abdomen and weighed to obtain the initial weight. Each individual was placed in a wind tunnel at 2.4 s cm^-1, and weighed every five minute until individuals reached 70-65% of their initial weight (Titon Jr et al., 2010). The rate of EWL was calculated from the value of the regression of the weight lost over time and expressed in mg min^-1. EWL corrected for surface area was calculated by dividing by 2/3 of the total surface area which corresponds to the region of the body in contact with the air flow (mg cm^-2 min^-1) (Titon Jr and Gomes, 2015; Titon Jr et al., 2010; Withers et al., 1982).

To estimate rate of WU, Individuals where placed in containers with absorbent towels moistened with enough water to cover the ventral region after the EWL experiment. The individuals were dried and weighed every four minutes for six consecutive times (Anderson et al., 2017; Titon Jr et al., 2010). The WU rate per area was calculated from the value of the regression of the body weight gained as a function of time (mg min^-1), and in turn divided by 1/3 of the total surface area that corresponds to the ventral area in contact with water during absorption, the values were expressed in units of mg cm^-2 min^-1. All experiments were conducted at an average air temperature of 25.55 ± 0.28°C. The relative humidity was not measured but approximated based on the location and time of the experiment (R. Solano pers. comm.)

The method for Thoropa miliaris is described below:

14 Thoropa miliaris were collected around Ubatuba (approximate location 23°26′13″S 45°04′08″W), and 15 individuals were collected in the Serra do Mar State Park, Picinguaba Nucleus (approximate location 23°27’50”S to 23°15’00”S and 45°15’00”W to 44°43’30”W). Animals were brought back to the University of São Paulo campus and kept individually in plastic enclosures (17 x 18 x 107 cm) for 5-7 days prior to estimating EWL. Animals were fed captive cockroaches (Nauphoeta cinérea) twice a week and were provided fresh water ad libitum (tap water) in a plastic cup. To estimate EWL, fully hydrated individuals were briefly dried with absorbent paper and the urine was removed by lightly pressing their abdomen and weighed to obtain the initial weight. Each individual was placed in a wind tunnel at 235 s cm^-1, and weighed every five minute until individuals reached 90% of their initial weight or the animal did not regain posture when turned over. Each individual measurement was replicated twice over the (3 days in between measurements) and averaged. The rate of EWL from the average values was calculated and presented as in mg h^-1. EWL corrected for surface area was calculated by dividing by 2/3 of the total surface area which corresponds to the region of the body in contact with the air flow (mg cm^-2 min^-1). The experiment was conducted at an air temperature between 22-25°C at a relative humidity of 50-60%.

Data provided by R. Parelli Bovo (rpbovo@gmail.com) was originally used for a separate paper which is currently in development. However, he has kindly offered to provide the data for this study pre-submission. The species that were collected by R. Parelli Bovo are provided below:

• Adenomera marmorata
• Ameerega flavopicta
• Aplastodiscus leucopygius
• Boana pardalis
• Boana raniceps
• Bokermannohyla alvarengai
• Bokermannohyla circumdata
• Bokermannohyla saxicola
• Brachycephalus pitanga
• Dendropsophus microps
• Dendropsophus sanborni
• Dermatonotus muelleri
• Elachistocleis cesarii
• Haddadus binotatus
• Ischnocnema guentheri
• Ischnocnema parva
• Leptodactylus fuscus
• Leptodactylus syphax
• Odontophrynus americanus
• Phyllomedusa megacephala
• Physalaemus atlanticus
• Physalaemus olfersii
• Pithecopus azureae
• Proceratophrys appendiculata
• Proceratophrys cururu
• Rhinella diptycha
• Rhinella ornata
• Rhinella rubescens
• Scinax fuscovarius
• Trachycephalus mesophaeus

Adults were collected along elevational gradients from sea level to 1,600 m a.s.l in Brazil’s Atlantic Forest and Cerrado during warm/wet seasons (September to February of 2011–2017). Animals were transported to the laboratory (613 m: -22.397331° S, 47.547799° W). All animals were kept individually in plastic terraria, with shelter and tap water ad libitum, under natural thermal/humidity regimes (daily temperature range of 22–27 °C, air humidity of 45–65%) and photoperiod (light/dark cycles 13.5:10.5).

Measurements were taken once for each individual, in the following sequence: EWL and WU within 3-5 days after field collection. Animals were not fed as the measurements were taken in a fasted state. Before and after any experimental run, animals were checked for motor coordination, skin color, posture, and responsiveness. When individuals failed these checks, they were not include in the dataset and subsequent analyses.

EWL and WU experiments were conducted following Anderson et al. (2017), Gouveia et al. (2019), and Bovo et al. (2023). Briefly, animals’ temperature and hydration state were standardised before the trials by holding each amphibian in individual PVC containers, filled with 0.5 cm (body sizes < 5 g) or 1 cm of water (body sizes > 5 g), and placed inside a climate-controlled incubator (122FC model - Eletrolab) at 25 °C for 1 h. To obtain the initial standardised body mass (100% hydrated), frogs were gently dried with absorbent paper and their abdomen were lightly pressed to remove urine from their urinary bladder. Then, frogs were individually placed in PVC containers (1,000 ml) and EWL rates were measured in a typical open-flow system at 25 °C, containing a mass flow meter (SS-3 Subsampler, Sable Systems) supplying a stable airflow at 21.66 cm³ s^-1 (1,300 ml min^-1), a RH/Dewpoint Controller (DG-4, Sable Systems) standardizing the incurrent relative humidity at 30%, and a water vapor analyzer (RH-300 RH/Dewpoint Analyzer, Sable Systems) quantifying incurrent and excurrent air. All equipment was interfaced to a computer by an analog/digital unit (UI2, Sable Systems) to record changes in airflow and water vapor density (WVD) every 1.0 second. Skin body temperatures were taken at the end of the experiments, typically 30-40 min for frogs < 5 g and 60 min for frogs > 5 g, for quantification of the WVD at saturation at the evaporating (skin) surface. To quantify the EWL rates, the water vapor density deficit was calculated by the difference between an empty container and one containing the animal (Spotila and Berman 1976). Then, total transepithelial EWL was corrected for unit area of exposed skin surface (2/3 of the total surface; McClanahan Jr and Baldwin (1969)) and expressed as mass of water per exposed surface area per time (e.g. μg H₂O cm⁻² s⁻¹). To minimize the animal’s activity during EWL measurements, experiments were performed during the day (opposite to their natural active period) with chambers in the dark. To detect any significant changes in behavior and posture that could affect the records, individuals were often visually inspected during trials. Data were discarded if an animal urinated during the experiment.

Immediately after the EWL trials, where animals usually lost 98-70 % of their initial body masses (see dataset), each individual was placed in a container (a Petri dish for body sizes smaller than 5 g, or a circular PVC container for body sizes larger than 5 g) filled with tap water at a depth sufficient to cover their ventral region. Animals were taken from the container and carefully blotted with paper tissue and weighed (± 0.0001 g or 0.01 g) every 2 minutes for six consecutive times in a room at 25 °C. WU rates were calculated from the linear regression between body mass increments against time. Then, using the estimated surface area in contact with water (1/3 of the total surface, McClanahan Jr and Baldwin (1969)), the WU rate was calculated and expressed per unit area (e.g. μg H₂O cm⁻² s⁻¹).

Data for Breviceps montanus provided by S. Clusella-Trullas (sct333@sun.ac.za) will be used for a manuscript that is currently in submission. However, she has kindly offered to provide the data for this study pre-submission.

The raw data, including collection sites, are provided on GitHub.

---

Calculations and conversions

Evaporative water loss

The relative skin resistance, \(r_i\), is calculated as:

\[\begin{equation} r_i = r_t - r_b, \tag{1} \end{equation}\]

where \(r_t\) is the total resistance (s cm^-1) and \(r_b\) is the boundary layer resistance based on an equivalent agar model (s cm^-1) or hypothetical biophysical model (4). \(r_t\) is calculated following Spotila and Berman (1976) and Hillman et al. (2009):

\[\begin{equation} r_t = \frac{\rho}{\textrm{CWL}} = (\rho_{v[skin]} - \rho_v \times \textrm{RH}_{[excurrent]}) \times \textrm{CWL}, \tag{2} \end{equation}\]

where \(\rho\) is the vapour density gradient at the surface of the animal (g cm^-3), \(\textrm{CWL}\) is the cutaneous rates of evaporative water loss (\(\textrm{EWL}\)) by surface area (g s^-1 cm^-2), \(\rho_v\) is the water vapour density of water-saturated air (g cm^-3), \(\rho_{v[skin]}\) is the \(\rho_v\) at the skin surface, and \(\textrm{RH}\) is the fractional saturation (0-1). \(\rho\) is calculated as the difference in the water vapour density at the temperature of the evaporating surface by the water vapour density in the air at \(\textrm{RH}\) as fractional saturation (Feder and Burggren, 1992).

\(\textrm{EWL}\) (g s^-1) can be measured gravimetrically (change in mass over time) or changes in water vapour pressure, \(e\) (kPa), or \(\textrm{RH}\) between the incurrent and excurrent air for a flow-through system. In a flow-through system, \(\textrm{EWL}\) can be calculated from Hillman et al. (2009) and Riddell et al. (2017):

\[\begin{equation} \textrm{EWL} = \frac{e}{T \times R_v} \times \textrm{FR} = [(\rho_v \times \textrm{RH}_{[incurrent]}) - (\rho_v \times \textrm{RH}_{[excurrent]})] \times \textrm{FR}, \tag{3} \end{equation}\]

where \(T\) is temperature in Kelvin (K), \(R_v\) is the gas constant for water vapor (461.5 J K^-1 kg^-1), and \(\textrm{FR}\) is the flow rate (ml s^-1)

Boundary layer resistance

The boundary layer resistance is the layer of air in which an object (or organism) exchanges heat and mass with its surrounding environment. The boundary layer adds to the resistance of water loss by changing the physical conditions of the microclimate around the organism (Senzano et al., 2022). If \(r_b\) was not reported or not estimated by a representative agar model, we used the theoretical \(r_b\) with known forced convection (\(V\)) from Riddell et al. (2017):

\[\begin{equation} r_b = \frac{0.93 \times \rho \times C_{\rho}}{h_c}, \tag{4} \end{equation}\]

where, \(\rho\) is the density of the air (kg m^-3), \(C_{\rho}\) is the specific heat of air (kg^-1 K^-1) calculated as:

\[\begin{equation} C_{\rho} = \frac{1004.84 + (1846.40 \times r_w)}{1 + r_w}, \tag{5} \end{equation}\]

where \(r_w\) is mixing ratio of water vapour (Tracy et al., 2019) estimated using:

\[\begin{equation} r_w = \frac{0.62197 \times 1.0053 \times e}{p - 1.0053 \times e}, \tag{6} \end{equation}\]

where \(e\) is the vapour pressure (kPa) and \(p\) is the atmospheric pressure (kPa). \(h_c\) (4) is the coefficient of convective-heat transfer under forced convection approximated by:

\[\begin{equation} h_c = 0.0923(V^{0.0333} \times D^{-0.666}), \tag{7} \end{equation}\]

where \(V\) is the wind speed (m s^-1) and \(D\) is the characteristic dimension (m). \(D\) is the cross-sectional radius, assuming the dimensions of a sphere for a frog, based on the mass and the total surface area.

If the study was conducted under free convection, we used a modified (7) from Riddell et al. (2017):

\[\begin{equation} h_c = \frac{0.48 \times \textrm{Gr}^{0.25} \times k}{D}, \tag{8} \end{equation}\]

where \(\textrm{Gr}\) is the Grashof number, \(k\) is the thermal conductivity of the fluid (W m^-1 K^-1) \(D\) is the characteristic dimension (m). \(\textrm{Gr}\) describes the relative strength of the buoyant forces to viscous force of the air surrounding the object calculated as:

\[\begin{equation} \textrm{Gr} = \frac{g \times \beta \times \Delta{T} \times D^{3}}{v^{2}}, \tag{9} \end{equation}\]

where \(g\) is the acceleration by gravity (9.80 m s^-2), \(\beta\) is the coefficient of volumetric expansion (3.67 \(\times\) 10^-3 °C^-1), \(D\) is the characteristic dimension, \(v\) is the kinematic viscosity from Riddell et al. (2017), and \(\Delta{T}\) is estimated by:

\[\begin{equation} \Delta{T} = T_0(1+0.38 \times e_0 / p) - T(1+0.38 \times e/p), \tag{10} \end{equation}\]

where \(T_0\) is the temperature (K) of the surface, \(e_0\) is the water vapour pressure directly above the surface, and \(p\) is the atmospheric pressure (Monteith and Unsworth, 2013).

Vapour pressure deficit

The vapour pressure deficit (VPD) was calculated as:

\[\begin{equation} VPD = e_s - e_a, \tag{11} \end{equation}\]

where \(e_s\) is the saturation vapour pressure (kPa) at a given temperature and \(e_a\) is the actual vapour pressure in the air (kPa). \(e_s\) was calculated using the Clausius-Clapeyron equation Stull and Ahrens (2000):

\[\begin{equation} e_s = e_0 \times \exp \left[\frac{L}{R_v} \left( \frac{1}{T_0} - \frac{1}{T} \right) \right], \tag{12} \end{equation}\]

where \(e_0\) = 0.611 kPa and \(T_0\) – 273 K are constant parameters, and \(R_v\) – 461.5 J K^-1 Kg^-1 and \(L\) = 2.5 x 106 J kg^-1 are gas constant for water vapour and the latent heat of vapourization, respectively. \(T\) (as K) is the actual temperature measured. \(e_a\) was approximated as:

\[\begin{equation} e_a = \frac{\textrm{RH} \times e_s}{100}, \tag{13} \end{equation}\]

where \(\textrm{RH}\) is the measured relative humidity (%).

Surface area

For studies that presented EWL and WU without correcting for surface area, we corrected to surface area by dividing the absolute rates by 2/3 of the total surface area for EWL, which corresponds to the area exposed to air when anurans keep the water conservation posture (Withers et al., 1982; Withers et al., 1984). For WU, we corrected to surface area by dividing the absolute rates by 1/3 of the total surface area because WU primarily occurs in the pelvic ‘seat patch’ region of the ventral skin (Baldwin, 1974; McClanahan Jr and Baldwin, 1969; Willumsen et al., 2007). Total surface area (cm²) was estimated via family-specific mass-surface area scaling relationship from Klein et al. (2016):

\[\begin{equation} SA = \beta_0 M^{\beta_1}, \tag{14} \end{equation}\]

where \(\beta_0\) is the intercept, \(M\) is the animal body mass (g), and \(\beta_1\) is the slope. If the family-specific mass-surface area scaling relationship was not presented, we used the power equation for all Anura (9.8537\(M\)^0.6745) as a conservative estimate.

Percentage mass change to actual mass change

For studies that presented % change in mass (\(\Delta M\)) with raw data available on initial mass, we converted percentage change in mass (% h^-1) back to absolute mass change (g h^-1) following:

\[\begin{equation} \textrm{EWL} = M \times \left(\frac{\Delta M}{100}\right), \tag{15} \end{equation}\]

where \(M\) is the initial body mass (g), and \(\Delta M\) is the percentage change in mass (% h^-1). Note, studies with % change in body mass without the initial body mass presented were not included in the data extraction.

---

Dataset

Load and clean the raw dataset raw_data.csv.

# Load and clean raw data
raw_dat <- read.csv(file.path(data_path, "raw_data.csv")) %>%
  dplyr::select(study_ID:unit, r_s_est, calculated) %>%
  dplyr::mutate(ecotype  = factor(ecotype),
         family   = factor(family),
         origin   = factor(origin),
         strategy = factor(case_when(strategy == "" ~ "none", TRUE ~ as.character(strategy))),
         strategy = fct_relevel(strategy, "none", "water-proof", "cocoon", "hollow"),
         trait    = factor(trait),
         response = factor(response),
         lnMass   = log(mean_mass_g),
         lnFlow   = log(airflow_cm_s + 1),
         es_kPa   = ifelse(trait == "water loss", 0.611 * exp(2500000 / 461.5 * (1 / 273 - 1 / (trt_temp + 273.15))), NA), # saturation vapor pressure (kPa) at a given temperature
         ea_kPa   = ifelse(trait == "water loss", RH_perc * es_kPa / 100, NA), # actual vapor pressure (kPa)
         VPD_kPa  = es_kPa - ea_kPa,
         lnVPD    = log(VPD_kPa)) %>%
  dplyr::filter(species_phylo != "") # remove rows with no species

ewl_dat <- raw_dat %>%
  dplyr::filter(response == "evaporative water loss" & !is.na(unit_corrected_mean)) %>%
  dplyr::mutate(mg_h_mean = unit_corrected_mean * dors_SA_cm2,
                mg_h_sd   = unit_corrected_sd * dors_SA_cm2,    
                lnMean    = log(mg_h_mean),
                v         = mg_h_sd^2 / sample_size, # sampling variance (v)
                sei       = sqrt(v), # standard error (SE)
                inv       = 1 / sei, # precision (inverse of SE) )
                w         = 1 / v, # weight (inverse of variance) 
                h_70      = (mean_mass_g - (mean_mass_g * 0.7)) / (mg_h_mean * 0.001)) 

resist_dat <- raw_dat %>%
  dplyr::filter(!is.na(r_s_est)) %>%
  dplyr::mutate(lnMean = log(r_s_est + 1))

wu_dat <- raw_dat %>%
  dplyr::filter(trait == "water gain" & !is.na(unit_corrected_mean)) %>%
  dplyr::mutate(mg_h_mean = unit_corrected_mean * vent_SA_cm2,
                mg_h_sd   = unit_corrected_sd * vent_SA_cm2,
                lnMean    = log(mg_h_mean),
                v         = mg_h_sd^2 / sample_size, # sampling variance (v)
                sei       = sqrt(v), # standard error (SE)
                inv       = 1 / sei) # precision (inverse of SE)

#ewl_dat %>%
  #group_by(strategy) %>%
  #summarise(mean = mean(h_70, na.rm = TRUE))

There are 238 species with data on evaporative water loss (102 studies), 224 species with skin resistance data (81 studies), and 121 species with water uptake data (51 studies). We used the evaporative water loss to analyse differences between ecotype as there are more species represented compared to the skin resistance data. There were three species that were excluded from the analysis because they were not listed in the IUCN Red List and were not present in the phylogenetic tree from Jetz and Pyron (2018): Brachycephalus pitanga, Elachistocleis cesarii, and Leptodactylus luctator.

Geographical bias in hydroregulation studies

Fig. S4. Spatial distribution of hydroregulation studies using wild caught anurans (study n = 113). 19 studies used captive raised anurans. There are large regions around central Africa and Eurasia with high amphibian diversity, but no hydrological studies conducted (Fig. 1b).

---

Prepare phylogeny for analysis

The phylogeny was obtained from Jetz and Pyron (2018) comprising of 7,238 species. The tree was pruned to match the 256 species extracted for the subsequent analysis.

# Load tree
phylo_tree <- ape::read.tree(file.path(data_path, "amph_shl_new_Consensus_7238.tre"))

# Pruning data and phylogeny
tree_tip_label <- phylo_tree$tip.label # extract tree tip names
sp_list        <- raw_dat$species_phylo # extract species name from mean data
pruned_tree    <- ape::drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, sp_list)) # prune phylo_tree to keep species from the raw data
pruned_tree    <- phytools::force.ultrametric(pruned_tree, method = "extend") # ultrametricize the tree
phylo_cor      <- vcv(pruned_tree, cor = T)

setdiff(raw_dat$species_phylo, rownames(phylo_cor))
setdiff(rownames(phylo_cor), raw_dat$species_phylo)

#ape::is.ultrametric(pruned_tree)

Analysis

The original model incorporated origin (lab-raised or wild-caught) and whether the bladder was void of urine prior to the experiment. However, the model’s bulk effective samples size (ESS) was too low, indicating posterior means and medians may be unreliable. Since origin and whether the bladder was emptied did not influence \(r_i\), they were excluded in the final model. For transparency, on average, wild-caught anurans had, on average, higher \(r_i\) relative to lab-raised anurans (0.08 [-0.47:0.31]), and anurans with their bladder voided of urine had higher \(r_i\) (0.17 [-0.59: 0.25]). However there is substantial variability between the wild-caught and lab-raised groups, and whether anurans with their bladder voided or not.

# options(brms.backend = 'cmdstanr') # Error from using Rstans 'error in
# unserialize(socklist[[n]]) : error reading from connection'. Used cmdstanr
# around it.

priors <- c(prior(normal(0, 3), "b"), prior(normal(0, 3), "Intercept"), prior(student_t(3,
    0, 10), "sd"), prior(student_t(3, 0, 10), "sigma"))

ri_model <- brms::brm(lnMean ~ ecotype + strategy + lnMass + lnVPD + lnFlow + (1 |
    study_ID) + (1 + lnMass | species_iucn) + (1 | gr(species_phylo, cov = phylo)),
    data = resist_dat, family = gaussian(), data2 = list(phylo = phylo_cor), prior = priors,
    chains = 4, cores = 4, iter = 5000, warmup = 2500, control = list(adapt_delta = 0.99,
        max_treedepth = 15))

wu_model <- brms::brm(lnMean ~ ecotype + lnMass + trt_temp + hydration + origin +
    (1 | study_ID) + (1 + lnMass | species_iucn) + (1 | gr(species_phylo, cov = phylo)),
    data = wu_dat %>%
        dplyr::group_by(ecotype) %>%
        dplyr::filter(n() >= 5), family = gaussian(), data2 = list(phylo = phylo_cor),
    prior = priors, chains = 4, cores = 4, iter = 5000, warmup = 2500, control = list(adapt_delta = 0.99,
        max_treedepth = 15))

Fig. S5. Scatterplots of the observed data (y) vs the average simulated data (y_rep) from the posterior predictive distribution for the (a) resistance to water loss model, and (b) the water uptake model. Dashed line represents a slope of 1.

Model output

Table S3 - \(r_i\) model

Table S3. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the \(r_i\) model, which includes the intercept (\(\beta_0\)), ecotype, strategy, the natural logarithm of body mass (lnMass), the natural logarithm of vapour pressure deficit (lnVPD), and the natural logarithm of the experimental flow rate (lnFlow). Group-level effects include the standard deviations (\(\sigma\)) for study-level observations (\(\sigma_{study}^2\)), phylogenetic relatedness (\(\sigma_{phylogeny}^2\)), and the correlation among species (\(\sigma_{species}^2\)) and body mass (\(\sigma_{lnMass}^2\)). \(R_{marginal}^2\) represents the variance explained by fixed effects, while \(R_{conditional}^2\) represents the variance explained by both fixed effects and group-level effects.

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
ln\(\beta_0\)	1.11	0.55	0.01	2.18
Arboreal	0.73	0.54	-0.36	1.79
Fossorial	0.05	0.55	-1.01	1.13
Ground-dwelling	0.15	0.54	-0.90	1.22
Semi-aquatic	-0.02	0.54	-1.08	1.05
Stream-dwelling	0.01	0.58	-1.11	1.16
Water-proof	1.88	0.17	1.54	2.23
Cocoon	3.03	0.14	2.76	3.31
Hollow	1.65	0.44	0.79	2.51
lnMass	0.11	0.03	0.04	0.17
lnVPD	-0.01	0.05	-0.11	0.10
lnFlow	-0.22	0.06	-0.34	-0.10
Group-level effects
\(\sigma_{species}^2\)	0.35	0.12	0.06	0.57
\(\sigma_{lnMass}^2\)	0.11	0.05	0.01	0.19
cor(\(\sigma_{species}^2\), \(\sigma_{lnMass}^2\))	-0.58	0.42	-0.99	0.65
\(\sigma_{phylogeny}^2\)	0.35	0.13	0.07	0.61
\(\sigma_{study}^2\)	0.82	0.09	0.67	1.01
Phylogenetic signal
\(\lambda\)	0.37	0.17	0.03	0.65
Variance
\(R_{marginal}^2\)	64.18		59.67	67.79
\(R_{conditional}^2\)	88.99		86.95	90.83

---

Table S4 - WU model

Table S4. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the water uptake model, which includes the intercept (\(\beta_0\)), ecotype, the natural logarithm of body mass (lnMass), treatment temperature, initial hydration level, and origin (lab-raised or wild caught). Group-level effects include the standard deviations (\(\sigma\)) for study-level observations (\(\sigma_{study}^2\)), phylogenetic relatedness (\(\sigma_{phylogeny}^2\)), and the correlation among species (\(\sigma_{species}^2\)) and body mass (\(\sigma_{lnMass}^2\)). \(R_{marginal}^2\) represents the variance explained by fixed effects, while \(R_{conditional}^2\) represents the variance explained by both fixed effects and group-level effects.

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
ln\(\beta_0\)	8.90	0.85	7.22	10.56
Fossorial	-0.42	0.29	-0.98	0.16
Ground-dwelling	-0.19	0.22	-0.62	0.23
Semi-aquatic	-0.17	0.29	-0.72	0.41
ecotypeStream-dwelling	-0.62	0.30	-1.22	-0.02
lnMass	0.80	0.05	0.70	0.91
Treatment temperature	0.01	0.01	-0.01	0.04
Initial hydration	-0.05	0.01	-0.06	-0.04
Origin - Wild	0.02	0.42	-0.83	0.84
Group-level effects
\(\sigma_{species}^2\)	0.24	0.14	0.01	0.52
\(\sigma_{lnMass}^2\)	0.09	0.05	0.01	0.20
cor(\(\sigma_{species}^2\), \(\sigma_{lnMass}^2\))	0.11	0.52	-0.86	0.95
\(\sigma_{phylogeny}^2\)	0.58	0.19	0.18	0.94
\(\sigma_{study}^2\)	1.03	0.13	0.80	1.31
Phylogenetic signal
\(\lambda\)	0.77	0.16	0.30	0.93
Variance
\(R_{marginal}^2\)	69.65		62.90	75.05
\(R_{conditional}^2\)	97.68		96.51	98.52

---

Fig. S6 - \(r_i\) model

Fig. S6. Differences in resistance to water loss, \(r_i\) (s cm^-1) by (a) ecotype and (b) water-conserving strategies. Note, when plotting by ecotype, the \(r_i\) excludes behavioural water-conserving strategies such as during cocoon-forming and inside hollows. Mean estimates ± 95% CI presented in black points and error bars, while raw values were presented as grey points.

---

Fig. S7 - WU model

Fig. S7. Differences in cutaneous water uptake (mg H₂O h^-1) by ecotype. Mean estimates ± 95% CI presented in black points and error bars, while raw values were presented as grey points. The size of the grey points indicates study precision (inverse of standard error).

---

NicheMapR

To estimate the influence of warming and drought on activity of a hypothetical frog, we simulated a water and heat energy exchange model (Fig. S8a) and its interaction with a simulated local microclimate using the NicheMapR package (Kearney and Porter, 2017; Kearney and Porter, 2020).

Fig. S8. Summary water and energy exchange model from Tracy (1976) integrated into NicheMapR and water conserving strategies. (a) Schematic summary of the exchanges of energy and water between a frog and its environment used to develop the transient-state model of water exchange. In respect to water exchange, the net water loss represents water loss from respiratory, cutaneous, ocular, and cloaca evaporation as well as urinary and faecal water (Pirtle et al., 2019). (b) Behavioural, morphological, and physiological strategies employed by frogs on land to reduce water loss (Hillman et al., 2009).

Simulate water-conserving behaviours

We used the modularized version of NicheMapR’s v3.2.1 ectotherm model (i.e. ectoR_devel) and coded behavioral functions to account for the influence of hydration on activity (behav_functions function) which can be found in the behav_functions.R file on GitHub.

When the animal is not active, it is simulated to go to an underground retreat. It selects the shallowest depth with temperatures between Tmax and Tmin. If water = T, the frog selects the shallowest node with temperatures between Tmax and Tmin and a water potential >= -72.5, which was reported as a soil water potential from which Rana pipiens could absorb water. When the animal is belowground, it re-hydrates if the soil water potential is >= -72.5, at a rate specified in hyd.rate and proportional to the level of dehydration, following:

\[\begin{equation} \frac{Hyd_{max} - Hyd_{i}}{Hyd_{max}} \times hyd.rate, \tag{16} \end{equation}\]

It is important to note that currently frogs are not re-hydrating when active aboveground.

The “skinwet” term (\(p_{wet}\)) determines the proportion of the total surface area used in the calculation of mass transfer of water from the surface. Here, \(p_{wet}\) was calculated from empirically derived skin resistance following Pirtle et al. (2019):

\[\begin{equation} p_{wet} = \frac{1}{h_D \times r_i + (\frac{P_r}{S_c})^\frac{2}{3}}, \tag{17} \end{equation}\]

where \(r_i\) is the skin’s resistance to water vapor transfer (s m^-1), \(P_r\) is the Prandtl number (dimensionless), \(S_c\) is the Schmidt number (dimensionless), \(h_D\) is the mass transfer coefficient (m s^-1). To calculate \(p_{wet}\), the mass transfer coefficient (\(h_D\)) must be known. Mass transfer refers to the movement of a substance though a fluid interface, driven by changes in the concentration gradient. The mass transfer coefficient controls the rate of diffusion and is dependent on velocity, temperature, and physical properties of the interphase. Similar to the heat transfer coefficient, the mass transfer coefficient also has two components derived from free and forced convection. We consider forced convection only (as free convention should be negligible within a measurement chamber), and this can be calculated from the heat transfer coefficient:

\[\begin{equation} h_D = (\frac{h_C}{C_p \times \rho}) \times (\frac{P_r}{S_c}) ^ \frac{2}{3}, \tag{18} \end{equation}\]

where \(P_r\) is the Prandtl number (dimensionless), \(S_c\) is the Schmidt number (dimensionless), \(h_C\) is the heat transfer coefficient (J s^-1 cm^-2 K^-1), \(\rho\) is the density of dry air (g cm^-3), and \(C_p\) is the specific heat of air (1.01 J g^-1 K^-1).

Simulate drought and warming

Microclimates represent the physical environments experienced by an organism. They are a necessity for mechanistic niche modelling because it is the environment experienced at the scale of the individual that needs to be provided to the equations of energy and mass balance. The microclimate model was implemented as described by Kearney et al. (2014). Specifically, it was driven by historical 0.05° grid (~5 km) daily weather input layers (air temperature, vapor pressure, wind speed, and cloud cover).

We simulated meteorological drought (defined as less than average rainfall) from the micro_era5 function from NicheMapR by decreasing the rainmult function by 0.5 (Fig. S9). A factor of 0.5 means 50% of the original rainfall. To simulate frog activity under drought and warming scenarios, we constructed four climate conditions using the micro_era5 function:

1. current normal scenario representing the mean annual air temperature and rainfall from 1981–2010,
2. current drought scenario representing the mean annual air temperature from 1981-2010 and annual rainfall based on 2017-2020 drought in Australia (Fig. S9),
3. warming normal scenario representing the “business-as-usual” scenario, where the global surface temperature is estimated to increase by 4°C by 2080–2100 with no effect of drought (+4°C only), and
4. warming drought scenario with a +4°C increasing in air temperature and annual rainfall based on 2017-2020 drought in Australia.

First, download the mcera5 microclimate data to a local directory to run the micro_era5 function faster.

# get ERA5 data with package mcera5 
# assign your credentials (register here: https://cds.climate.copernicus.eu/user/register)
uid         <- "$$$$$$"
cds_api_key <- "$$$$$$$$-$$$$-$$$$-$$$$-$$$$$$$$$$$$"
ecmwfr::wf_set_key(user = uid, key = cds_api_key, service = "cds")

# bounding coordinates (in WGS84 / EPSG:4326)
c(-41.190, -20.186) 
xmn <- -42
xmx <- -41
ymn <- -21
ymx <- -20

xmn <- 152
xmx <- 154
ymn <- -28
ymx <- -26

# temporal extent
st_time <- lubridate::ymd("2016:01:01") # earliest sampling date
en_time <- lubridate::ymd("2018:12:31") # latest sampling date

# filename and location for downloaded .nc files
file_prefix <- "era5"
op <- "YOUR DIRECTORY"

# build a request (covering multiple years)
req <- mcera5::build_era5_request(xmin = xmn, xmax = xmx,
                                  ymin = ymn, ymax = ymx,
                                  start_time   = st_time,
                                  end_time     = en_time,
                                  outfile_name = file_prefix)

mcera5::request_era5(request = req, uid = uid, out_path = op)

Next, run the micro_era5 function for the four scenarios in 2017 (representing typical rainfall year), and one for 2019 to validate the observed rainfall data. Note, we ran the microclimate models from 2016 to 2017 to ‘spin up’ the soil moisture for 2016, and only used the microclimate data for 2017. This also applied to the 2019 microclimate, where we ran the microclimate models from 2018 to 2019, and used the microclimate data for 2019.

# Load microclimates
micro_curr_wet <- readRDS("micro_output/Australia_Current normal.rds")
micro_curr_dry <- readRDS("micro_output/Australia_Current dry.rds")
micro_warm_wet <- readRDS("micro_output/Australia_Warming normal.rds")
micro_warm_dry <- readRDS("micro_output/Australia_Warming dry.rds")

# Simulate rainfall for 2019 drought
act_dry <- micro_era5(loc = c(153.09249, -27.6235), # Karawatha, QLD.
                      runshade = 1, minshade = 0, # shade parameters
                      warm = 0, # current climate
                      dstart = "01/01/2018", 
                      dfinish = "31/12/2019",
                      spatial = '/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/Drought project/Spatial data/karawatha/era5', # change to your location
                      save = 0)

To verify our microclimate models, we plotted the simulated the daily rainfall for 2017 (representing typical rainfall) and 2019 (recent historical drought) with the observed daily rainfall in Australia for the following location: Karawatha, Southeast Queensland, Australia (153.09249, -27.6235). Karawatha provides a good case study because there are many ecotypes found in this location: ground-dwelling (e.g. Rhinella marina), arboreal (e.g. Litoria caerulea), fossorial (e.g. Cyclorana alboguttata), semi-aquatic (e.g. Litoria nasuta), and this area has experienced drought recently (2019). We extracted rainfall data from a weather station next to Karawatha which experienced “very much below average” rainfall from the Australian Bureau of Meteorology.

# Load observed rainfall for 2019
obs_rainfall <- read.csv(file.path(data_path, "obs_rainfall.csv")) %>%
  tibble::rowid_to_column("DOY") %>%
  dplyr::mutate(date = lubridate::make_date(year = 2019, month = month_num, day = day))

# Merge observed and predicted rainfall together
model_rain_wet <- as.data.frame(micro_curr_wet$RAINFALL) %>%
  tail(-366) %>% # remove 2016 simulation
  tibble::rowid_to_column("DOY") %>%
  rename(rainfall_mm = "micro_curr_wet$RAINFALL") %>%
  merge(obs_rainfall, by = "DOY", all.x = TRUE)

model_2019 <- as.data.frame(act_dry$RAINFALL) %>%
  tail(-366) %>% # remove 2018 simulation
  tibble::rowid_to_column("DOY") %>%
  rename(rainfall_mm = "act_dry$RAINFALL") %>%
  merge(obs_rainfall, by = "DOY", all.x = TRUE)

# 2017 normal validate
rain_wet_plot <- model_rain_wet %>%
  ggplot() +
  geom_line(aes(x = DOY, y = karawatha_17), colour = "grey") +
  geom_line(aes(x = DOY, y = rainfall_mm * 1.2), colour = "red", alpha = 0.5) +
  labs(x = "Day of the year", y = expression("Daily rainfall (mm day"^{"-1"}*")")) +
  ggtitle("2017 Normal") +
  scale_y_continuous(lim = c(0, 100), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(0, 365, 50), expand = c(0, 0)) + 
  mytheme() + theme(legend.position = "bottom")

# 2019 drought validate
rain_dry_plot <- model_2019 %>%
  ggplot() +
  geom_line(aes(x = DOY, y = karawatha_19), colour = "grey") +
  geom_line(aes(x = DOY, y = rainfall_mm), colour = "red", alpha = 0.5) +
  labs(x = "Day of the year", y = expression("Daily rainfall (mm day"^{"-1"}*")")) +
  ggtitle("2019 Drought") +
  scale_y_continuous(lim = c(0, 100), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(0, 365, 50), expand = c(0, 0)) + 
  mytheme() + theme(legend.position = "bottom")

cowplot::plot_grid(rain_wet_plot, rain_dry_plot, ncol = 2)

Fig. S9. Simulated rainfall from the micro_era5 function (red) and the observed rainfall for Karawatha, QLD from the Australian Bureau of Meteorology (grey). The total yearly observed rainfall for 2017 was as 1100 mm, and the total simulated yearly rainfall from NicheMapR was 1097.2 mm. The total yearly observed rainfall for the 2019 drought was 596 mm, and the total simulated yearly rainfall from NicheMapR was 556.22 mm. The average yearly rainfall across 1981-2010 is around 1100 mm per year.

Simulate activity - ecotypes

We simulated the potential number of hours for activity in a year (\(t_{act}\)) which represents the suitable thermal and hydric conditions for the animal to move beyond their retreat to either catch prey or finding mates (Kearney and Porter, 2020). The thermoregulatory and hydroregulatory sequence in the model assumes that frogs will only be ‘active’ (i.e. moving beyond their retreat, catching prey, finding mates) between the minimum and maximum body temperature thresholds for foraging, Tmin and Tmax (°C), and the minimum tolerated hydration, min.hyd (%). In addition, if water.act = T, the frog will only be active if temperatures are within the suitable temp range, and the frog is not expected to go below allowed hydration levels (e.g. if %hydration is not below min.hyd).

To estimate the potential \(t_{act}\) for one year, we simulated the hypothetical frog to be active during the day and night (24 hours). While most frogs are nocturnal, there are some frog species found moving during the day. The water hydric parameters (e.g. skin resistance, water uptake rate, dehydration tolerance) were based on the average values of the dataset collected from the PRISMA search. The thermal parameters (minimum and maximum foraging and critical temperature) were based on Rhinalla marina which has been verified in Kearney et al. (2008).

Three water-saving strategies were constructed to broadly reflect each ecotype that use either behavioural strategies such as microhabitat selection, or physiological strategies such as increased skin resistance or skin thickness (Table S5). For example, many amphibians (especially arid specialist) seek or burrow underground to regulate body temperature and water balance without exhibiting thicker skin or higher skin resistance. This is because underground burrows are generally cooler, less varied temperature fluctuation, and more moist relative to the surface climate.

Table S5. Modified parameters for the sim.ecto function for each hypothetical frog model. The skin resistance values were based on the average empirical measurements for a typical frog and a waterproof frog in the meta-analysis. Seek shade represents microclimate simulated under shaded conditions (0-90% shade). Retreat underground represents the ability for the frog to seek suitable microclimates underground. Can climb represents the ability for the frog to seek suitable microclimates above ground level (e.g. trees, cliffs). Note, frogs are not active when underground, and water uptake only occurs during inactivity.

Water-saving strategy	Skin resistance	Seek shade?	Retreat underground?	Can climb?	Ecotype
Shade only	Low	Yes	No	No	Ground-dwelling
Waterproof	High	Yes	No	Yes	Arboreal
Fossorial	Low	Yes	Yes	No	Fossorial

It is important to note that we focus only on the direct effects of environmental change on water budgets by predicting changes in EWL rates. The conclusions we draw do not account for the indirect effects of environmental change on water budgets, such as potential changes to thermal tolerance as a result of dehydration, or the thermoregulatory difficulties that may arise as a result of habitat modification.

source("/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/Drought project/code/behav_functions.R")
# Construct frog model
# compute the heat exchange by convection (extract mass transfer coefficient, Prandtl number and Schmidt number)
CONV_out <- NicheMapR::CONV_ENDO(TS     = 19, # skin temperature (°C)
                                 TENV   = 20, # fluid temperature (°C)
                                 SHAPE  = 4, # 4 is ellipsoid
                                 SURFAR = mean(resist_dat$dors_SA_cm2, na.rm = TRUE) / 10000,  # surface area for convection, m2
                                 FLTYPE = 0, # fluid type: 0 = air
                                 FURTST = 0, # test of presence of fur (length x diameter x density x depth) (-)
                                 D      = mean(resist_dat$D, na.rm = TRUE), 
                                 TFA    = 20, # initial fur/air interface temperature
                                 VEL    = mean(resist_dat$airflow_cm_s, na.rm = TRUE) / 100, # wind speed (m/s)
                                 ZFUR   = 0, # fur depth, mean (m)
                                 BP     = 101325, # barometric pressure at sea level
                                 ELEV   = 0) # elevation (m)

# basic parameters for 30 g frog
Ww_g         <- exp(mean(log(raw_dat$mean_mass_g), na.rm = TRUE)) # geometric mean wet weight of animal (g), account for uneven distribution
r_s_low      <- resist_dat %>% dplyr::filter(strategy == "none") %>% dplyr::select(r_s_est) # skin resistance
pct_wet_high <- 1 / (CONV_out[5] * mean(r_s_low$r_s_est, na.rm = TRUE) + (CONV_out[11] / CONV_out[13]) ^ 0.6666666) * 100  # % of surface area acting as a free-water exchanger (Pirtle et al 2017)

# parameters for a water-proof frog
r_s_high    <- resist_dat %>% dplyr::filter(strategy == "water-proof") %>% dplyr::select(r_s_est) # skin resistance
pct_wet_low <- 1 / (CONV_out[5] * max(r_s_high$r_s_est, na.rm = TRUE) + (CONV_out[11] / CONV_out[13]) ^ 0.6666666) * 100 # % of surface area acting as a free-water exchanger (Pirtle et al 2017)

# Thermal traits based on Rhinella marina
Tmin   <- 13.7 # minimum Tb at which activity occurs (Kearney et al 2008)
Tmax   <- 35.4 # maximum Tb at which activity occurs (Kearney et al 2008, Tingley & Shine 2011)
#T_pref <- 24 # preferred Tb (Kearney et al 2008)
CTmax  <- 37 # critical thermal minimum (affects choice of retreat) Tracy et al 2012
CTmin  <- 6 # critical thermal maximum (affects choice of retreat) Kolbe et al 2010, McCann et al 2014

# Water balance traits
min_hyd <- 80 # minimum tolerated hydration before activity declines (% of fully hydrated animals)
hyd.death <- 50 # minimum tolerated hydration before death (% of fully hydrated animals)
wu_rate <- wu_dat %>% dplyr::filter(strategy == "none") %>% dplyr::select(mg_h_mean)
hyd_rate <- exp(mean(log(wu_rate$mg_h_mean), na.rm = TRUE)) / 1000 # geometric mean rehydration rate (g/h), account for uneven distribution
# depends on current and max hydration like this: hyd.rate * ((hyd - hyd.current) / hyd)

# behav = 'diurnal', 'nocturnal' or 'both'
# water; does the frog select depth according to water potential? (TRUE or FALSE)
# water.act; does the activity depend on water loss? (TRUE or FALSE)

# SHADE MODEL
shad_curr_wet_mod <- sim.ecto(micro_curr_wet, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = FALSE, water.act = TRUE)

shad_curr_dry_mod <- sim.ecto(micro_curr_dry, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = FALSE, water.act = TRUE)

shad_warm_wet_mod <- sim.ecto(micro_warm_wet, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = FALSE, water.act = TRUE)

shad_warm_dry_mod <- sim.ecto(micro_warm_dry, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = FALSE, water.act = TRUE)

# WATER-PROOF MODEL
tree_curr_wet_mod <- sim.ecto(micro_curr_wet, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = TRUE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_low, 
                              water = TRUE, water.act = TRUE)

tree_curr_dry_mod <- sim.ecto(micro_curr_dry, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = TRUE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_low, 
                              water = TRUE, water.act = TRUE)

tree_warm_wet_mod <- sim.ecto(micro_warm_wet, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = TRUE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_low, 
                              water = TRUE, water.act = TRUE)

tree_warm_dry_mod <- sim.ecto(micro_warm_dry, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = FALSE, climb = TRUE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_low, 
                              water = TRUE, water.act = TRUE)

# BURROWING MODEL
burr_curr_wet_mod <- sim.ecto(micro_curr_wet, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = TRUE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = TRUE, water.act = TRUE)

burr_curr_dry_mod <- sim.ecto(micro_curr_dry, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = TRUE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = TRUE, water.act = TRUE)

burr_warm_wet_mod <- sim.ecto(micro_warm_wet, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = TRUE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = TRUE, water.act = TRUE)

burr_warm_dry_mod <- sim.ecto(micro_warm_dry, Ww_g = Ww_g, shape = 4, 
                              Tmax = Tmax, Tmin = Tmin, CTmin = CTmin, CTmax = CTmax,
                              behav = 'both', in.shade = TRUE, burrow = TRUE, climb = FALSE,
                              min.hyd = min_hyd, hyd.death = hyd.death,
                              hyd.rate = hyd_rate, pct_wet = pct_wet_high, 
                              water = TRUE, water.act = TRUE)

# SHADE MODEL
shad_curr_wet_df <- data.frame(shad_curr_wet_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = shad_curr_wet_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_wet = length(active[active == TRUE]))

shad_curr_dry_df <- data.frame(shad_curr_dry_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = shad_curr_dry_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_dry = length(active[active == TRUE]))

shad_warm_wet_df <- data.frame(shad_warm_wet_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = shad_warm_wet_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_wet = length(active[active == TRUE]))

shad_warm_dry_df <- data.frame(shad_warm_dry_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = shad_warm_dry_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_dry = length(active[active == TRUE]))

shad_model <- shad_curr_wet_df %>%
  merge(shad_curr_dry_df, by = "day") %>%
  merge(shad_warm_wet_df, by = "day") %>%
  merge(shad_warm_dry_df, by = "day") %>%
  pivot_longer(!day, names_to = "condition", values_to = "hours") %>%
  dplyr::mutate(condition = factor(condition, levels = c("warm_dry", "warm_wet", "curr_dry", "curr_wet")),
                season = case_when(
    day >= 335 ~ "summer",
    day >= 1 & day <= 61 ~ "summer",
    day >= 62 & day <= 153 ~ "autumn",
    day >= 154 & day <= 244 ~ "winter",
    day >= 245 & day <= 334 ~ "spring"
    ))

# TREE MODEL
tree_curr_wet_df <- data.frame(tree_curr_wet_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = tree_curr_wet_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_wet = length(active[active == TRUE]))

tree_curr_dry_df <- data.frame(tree_curr_dry_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = tree_curr_dry_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_dry = length(active[active == TRUE]))

tree_warm_wet_df <- data.frame(tree_warm_wet_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = tree_warm_wet_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_wet = length(active[active == TRUE]))

tree_warm_dry_df <- data.frame(tree_warm_dry_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = tree_warm_dry_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_dry = length(active[active == TRUE]))

tree_model <- tree_curr_wet_df %>%
  merge(tree_curr_dry_df, by = "day") %>%
  merge(tree_warm_wet_df, by = "day") %>%
  merge(tree_warm_dry_df, by = "day") %>%
  pivot_longer(!day, names_to = "condition", values_to = "hours") %>%
  dplyr::mutate(condition = factor(condition, levels = c("warm_dry", "warm_wet", "curr_dry", "curr_wet")),
                season = case_when(
    day >= 335 ~ "summer",
    day >= 1 & day <= 61 ~ "summer",
    day >= 62 & day <= 153 ~ "autumn",
    day >= 154 & day <= 244 ~ "winter",
    day >= 245 & day <= 334 ~ "spring"
    ))

# FOSSORIAL MODEL
burr_curr_wet_df <- data.frame(burr_curr_wet_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = burr_curr_wet_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_wet = length(active[active == TRUE]))

burr_curr_dry_df <- data.frame(burr_curr_dry_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = burr_curr_dry_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_dry = length(active[active == TRUE]))

burr_warm_wet_df <- data.frame(burr_warm_wet_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = burr_warm_wet_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_wet = length(active[active == TRUE]))

burr_warm_dry_df <- data.frame(burr_warm_dry_mod$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = burr_warm_dry_mod.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_dry = length(active[active == TRUE]))

burr_model <- burr_curr_wet_df %>%
  merge(burr_curr_dry_df, by = "day") %>%
  merge(burr_warm_wet_df, by = "day") %>%
  merge(burr_warm_dry_df, by = "day") %>%
  pivot_longer(!day, names_to = "condition", values_to = "hours") %>%
  dplyr::mutate(condition = factor(condition, levels = c("warm_dry", "warm_wet", "curr_dry", "curr_wet")),
                season = case_when(
    day >= 335 ~ "summer",
    day >= 1 & day <= 61 ~ "summer",
    day >= 62 & day <= 153 ~ "autumn",
    day >= 154 & day <= 244 ~ "winter",
    day >= 245 & day <= 334 ~ "spring"
    ))

Model output - ecotypes

Table S6a - Summary

Table S6a. Total potential activity hours per year (\(t_{act}\)) for a 8.81 g frog under different warming (current or warming) and drought (normal or drought) conditions with different water-saving strategies.

data.frame(strategy = c("Shade only", "Shade only", "Shade only", "Shade only",
                        "Waterproof", "Waterproof", "Waterproof", "Waterproof",
                        "Burrowing", "Burrowing", "Burrowing", "Burrowing"),
           temp = c("current", "current", "warming", "warming"),
           rain = c("normal", "drought"),
           t_act_h = c(sum(shad_curr_wet_df$curr_wet), sum(shad_curr_dry_df$curr_dry), sum(shad_warm_wet_df$warm_wet), sum(shad_warm_dry_df$warm_dry),
                       sum(tree_curr_wet_df$curr_wet), sum(tree_curr_dry_df$curr_dry), sum(tree_warm_wet_df$warm_wet), sum(tree_warm_dry_df$warm_dry),
                       sum(burr_curr_wet_df$curr_wet), sum(burr_curr_dry_df$curr_dry), sum(burr_warm_wet_df$warm_wet), sum(burr_warm_dry_df$warm_dry)),
           t_act_per = c(sum(shad_curr_wet_df$curr_wet) / 8760 * 100, sum(shad_curr_dry_df$curr_dry) / 8760 * 100, sum(shad_warm_wet_df$warm_wet) / 8760 * 100, sum(shad_warm_dry_df$warm_dry) / 8760 * 100,
                       sum(tree_curr_wet_df$curr_wet) / 8760 * 100, sum(tree_curr_dry_df$curr_dry) / 8760 * 100, sum(tree_warm_wet_df$warm_wet) / 8760 * 100, sum(tree_warm_dry_df$warm_dry) / 8760 * 100,
                       sum(burr_curr_wet_df$curr_wet) / 8760 * 100, sum(burr_curr_dry_df$curr_dry) / 8760 * 100, sum(burr_warm_wet_df$warm_wet) / 8760 * 100, sum(burr_warm_dry_df$warm_dry) / 8760 * 100)) %>%
  knitr::kable(col.names = c("Water-saving strategy", "Warming simulation", "Drought simulation", "$t_{act}$ (h)", "$t_{act}$ (%)"))

Water-saving strategy	Warming simulation	Drought simulation	\(t_{act}\) (h)	\(t_{act}\) (%)
Shade only	current	normal	3934	45
Shade only	current	drought	3729	43
Shade only	warming	normal	4268	49
Shade only	warming	drought	4022	46
Waterproof	current	normal	4724	54
Waterproof	current	drought	4482	51
Waterproof	warming	normal	4887	56
Waterproof	warming	drought	4596	52
Burrowing	current	normal	4271	49
Burrowing	current	drought	4152	47
Burrowing	warming	normal	4652	53
Burrowing	warming	drought	4535	52

---

Table S6b - Relative change year

Table S6b. Change in \(t_{act}\) (%) for each water-saving strategy relative to the current normal scenario under warming only, drought only, and warming and drought combined.

data.frame(strategy = c("Shade only","Waterproof", "Burrowing"),
           curr_wet = c(sum(shad_curr_wet_df$curr_wet), sum(tree_curr_wet_df$curr_wet), sum(burr_curr_wet_df$curr_wet)),
           curr_dry = c(sum(shad_curr_dry_df$curr_dry), sum(tree_curr_dry_df$curr_dry), sum(burr_curr_dry_df$curr_dry)),
           warm_wet = c(sum(shad_warm_wet_df$warm_wet), sum(tree_warm_wet_df$warm_wet), sum(burr_warm_wet_df$warm_wet)),
           warm_dry = c(sum(shad_warm_dry_df$warm_dry), sum(tree_warm_dry_df$warm_dry), sum(burr_warm_dry_df$warm_dry))) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  knitr::kable(col.names = c("Water-saving strategy", "$\\Delta$ warming only (%)", "$\\Delta$ drought only (%)", "$\\Delta$ warming and drought (%)"))

Water-saving strategy	\(\Delta\) warming only (%)	\(\Delta\) drought only (%)	\(\Delta\) warming and drought (%)
Shade only	8.5	-5.2	2.2
Waterproof	3.5	-5.1	-2.7
Burrowing	8.9	-2.8	6.2

data.frame(strategy = c("Shade only","Waterproof", "Burrowing"),
           curr_wet = c(sum(shad_curr_wet_df$curr_wet), sum(tree_curr_wet_df$curr_wet), sum(burr_curr_wet_df$curr_wet)),
           curr_dry = c(sum(shad_curr_dry_df$curr_dry), sum(tree_curr_dry_df$curr_dry), sum(burr_curr_dry_df$curr_dry)),
           warm_wet = c(sum(shad_warm_wet_df$warm_wet), sum(tree_warm_wet_df$warm_wet), sum(burr_warm_wet_df$warm_wet)),
           warm_dry = c(sum(shad_warm_dry_df$warm_dry), sum(tree_warm_dry_df$warm_dry), sum(burr_warm_dry_df$warm_dry))) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet),
                delta_dry      = (curr_dry - curr_wet),
                delta_warm_dry = (warm_dry - curr_wet)) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  knitr::kable(col.names = c("Water-saving strategy", "$\\Delta$ warming only (%)", "$\\Delta$ drought only (%)", "$\\Delta$ warming and drought (%)"))

Water-saving strategy	\(\Delta\) warming only (%)	\(\Delta\) drought only (%)	\(\Delta\) warming and drought (%)
Shade only	334	-205	88
Waterproof	163	-242	-128
Burrowing	381	-119	264

---

Table S6c - Relative change season

Table S6c. Change in \(t_{act}\) (%) for each water-saving strategy relative to the current normal scenario under warming only, drought only, and warming and drought combined.

# summer
shad_model_summer <- shad_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tree_model_summer <- tree_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

burr_model_summer <- burr_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_summer <- data.frame(strategy = c("Shade only","Waterproof", "Burrowing"),
           curr_wet = c(shad_model_summer$hours[4], tree_model_summer$hours[4], burr_model_summer$hours[4]),
           curr_dry = c(shad_model_summer$hours[3], tree_model_summer$hours[3], burr_model_summer$hours[3]),
           warm_wet = c(shad_model_summer$hours[2], tree_model_summer$hours[2], burr_model_summer$hours[2]),
           warm_dry = c(shad_model_summer$hours[1], tree_model_summer$hours[1], burr_model_summer$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(strategy = "**Summer**", .before = 1)

# autumn
shad_model_autumn <- shad_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tree_model_autumn <- tree_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

burr_model_autumn <- burr_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_autumn <- data.frame(strategy = c("Shade only","Waterproof", "Burrowing"),
           curr_wet = c(shad_model_autumn$hours[4], tree_model_autumn$hours[4], burr_model_autumn$hours[4]),
           curr_dry = c(shad_model_autumn$hours[3], tree_model_autumn$hours[3], burr_model_autumn$hours[3]),
           warm_wet = c(shad_model_autumn$hours[2], tree_model_autumn$hours[2], burr_model_autumn$hours[2]),
           warm_dry = c(shad_model_autumn$hours[1], tree_model_autumn$hours[1], burr_model_autumn$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(strategy = "**Autumn**", .before = 1)

# winter
shad_model_winter <- shad_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tree_model_winter <- tree_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

burr_model_winter <- burr_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_winter <- data.frame(strategy = c("Shade only","Waterproof", "Burrowing"),
           curr_wet = c(shad_model_winter$hours[4], tree_model_winter$hours[4], burr_model_winter$hours[4]),
           curr_dry = c(shad_model_winter$hours[3], tree_model_winter$hours[3], burr_model_winter$hours[3]),
           warm_wet = c(shad_model_winter$hours[2], tree_model_winter$hours[2], burr_model_winter$hours[2]),
           warm_dry = c(shad_model_winter$hours[1], tree_model_winter$hours[1], burr_model_winter$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(strategy = "**Winter**", .before = 1)

# spring
shad_model_spring <- shad_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tree_model_spring <- tree_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

burr_model_spring <- burr_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_spring <- data.frame(strategy = c("Shade only","Waterproof", "Burrowing"),
           curr_wet = c(shad_model_spring$hours[4], tree_model_spring$hours[4], burr_model_spring$hours[4]),
           curr_dry = c(shad_model_spring$hours[3], tree_model_spring$hours[3], burr_model_spring$hours[3]),
           warm_wet = c(shad_model_spring$hours[2], tree_model_spring$hours[2], burr_model_spring$hours[2]),
           warm_dry = c(shad_model_spring$hours[1], tree_model_spring$hours[1], burr_model_spring$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(strategy = "**Spring**", .before = 1)

# Render table
bind_rows(tab_summer, tab_autumn, tab_winter, tab_spring) %>% 
  remove_rownames() %>% 
  knitr::kable(col.names = c("Water-saving strategy", "$\\Delta$ warming only (%)", "$\\Delta$ drought only (%)", "$\\Delta$ warming and drought (%)"))

Water-saving strategy	\(\Delta\) warming only (%)	\(\Delta\) drought only (%)	\(\Delta\) warming and drought (%)
Summer
Shade only	-8.39	-13.92	-21.82
Waterproof	-8.61	-13.39	-22.78
Burrowing	-7.16	-6.58	-12.76
Autumn
Shade only	4.75	-1.84	3.26
Waterproof	1.14	-1.99	-0.54
Burrowing	5.87	-1.43	4.99
Winter
Shade only	72.88	-1.91	65.47
Waterproof	38.05	-1.56	33.52
Burrowing	70.55	-0.99	66.60
Spring
Shade only	0.77	-2.80	-4.44
Waterproof	-2.15	-3.65	-7.79
Burrowing	2.30	-1.20	0.83

---

Simulate activity - biomes

To test if changes in activity under different climate scenarios are generalisable, we also ran the shade-only model under three additional representative biome types where anurans typically occupy and are predicted to increase in drought: tropic (Salvador, Brazil), semi-arid (Tankwa Karoo, South Africa), and Mediterranean (Seville, Spain) (Fig. S10). The biome locations were classified based on the Köppen climate classification. All microclimate simulations ran for 2016-2017 and for the same hypothetical frog.

We first downloaded the microclimates locally with the function below (edited by M. Kearney).

source('../micro_era5_local2.R') # a version of micro_era5 that reads from locally downloaded tiles of the whole planet

# Define locations and scenarios
locations <- list(
  list(name = "Brazil", longlat = c(-39.16553, -13.05029)),
  list(name = "South Africa", longlat = c(19.8566, -32.4297)),
  list(name = "Spain", longlat = c(-5.5538, 37.04143)),
  list(name = "Australia", longlat = c(153.09249, -27.6235))
)

scenarios <- list(
  list(name = "Current normal", rainmult = 1.2, warm = 0),
  list(name = "Current dry", rainmult = 0.5, warm = 0),
  list(name = "Warming normal", rainmult = 1.2, warm = 4),
  list(name = "Warming dry", rainmult = 0.5, warm = 4)
)

# Function to run micro_era5 and save output
run_micro_era5 <- function(location=location, scenario=scenario) {
  # Extract parameters
  longlat <- location$longlat
  name <- paste(location$name, scenario$name, sep = "_")
  rainmult <- scenario$rainfact
  warm <- scenario$warm
  cat(name, rainmult, warm, '\n')
  # Run micro_era5
  result <- micro_era5_local(loc = longlat,
                       rainmult = rainmult,
                       warm = warm,
                       dstart = "01/01/2016", 
                       dfinish = "31/12/2017",
                       spatial = "/srv/6300-predecol/data/ERA5/time"
                       )

  # Save output
  output_file <- paste0(name, ".rds")
  saveRDS(result, file = output_file)
}

# Loop over locations and scenarios
for (i in 1:length(locations)) {
  for (j in 1:length(scenarios)) {
    run_micro_era5(locations[i][[1]], scenarios[j][[1]])
  }
}

Due to the long processing time for downloading microclimate for each site, we saved the microclimate locally, and loaded them here.

micro_curr_wet_br <- readRDS("micro_output/Brazil_Current normal.rds")
micro_curr_dry_br <- readRDS("micro_output/Brazil_Current dry.rds")
micro_warm_wet_br <- readRDS("micro_output/Brazil_Warming normal.rds")
micro_warm_dry_br <- readRDS("micro_output/Brazil_Warming dry.rds")

micro_curr_wet_sa <- readRDS("micro_output/South Africa_Current normal.rds")
micro_curr_dry_sa <- readRDS("micro_output/South Africa_Current dry.rds")
micro_warm_wet_sa <- readRDS("micro_output/South Africa_Warming normal.rds")
micro_warm_dry_sa <- readRDS("micro_output/South Africa_Warming dry.rds")

micro_curr_wet_sp <- readRDS("micro_output/Spain_Current normal.rds")
micro_curr_dry_sp <- readRDS("micro_output/Spain_Current dry.rds")
micro_warm_wet_sp <- readRDS("micro_output/Spain_Warming normal.rds")
micro_warm_dry_sp <- readRDS("micro_output/Spain_Warming dry.rds")

model_rain_wet_br <- as.data.frame(micro_curr_wet_br$RAINFALL) %>%
  tail(-366) %>% # remove 2016 simulation
  tibble::rowid_to_column("DOY") %>%
  rename(rainfall_mm = "micro_curr_wet_br$RAINFALL")

model_rain_wet_sa <- as.data.frame(micro_curr_wet_sa$RAINFALL) %>%
  tail(-366) %>% # remove 2016 simulation
  tibble::rowid_to_column("DOY") %>%
  rename(rainfall_mm = "micro_curr_wet_sa$RAINFALL")

model_rain_wet_sp <- as.data.frame(micro_curr_wet_sp$RAINFALL) %>%
  tail(-366) %>% # remove 2016 simulation
  tibble::rowid_to_column("DOY") %>%
  rename(rainfall_mm = "micro_curr_wet_sp$RAINFALL")

p1 <- ggplot() + 
  geom_line(data = model_rain_wet, aes(x = DOY, y = rainfall_mm * 1.2), colour = "darkgrey") +
  ylab("Daily rainfall (mm)") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(0, 100), expand = c(0,0), position = "right") + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

p2 <- ggplot() + 
  geom_line(data = model_rain_wet_br, aes(x = DOY, y = rainfall_mm * 1.2), colour = "darkgrey") +
  ylab("Daily rainfall (mm)") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(0, 100), expand = c(0,0), position = "right") + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

p3 <- ggplot() + 
  geom_line(data = model_rain_wet_sa, aes(x = DOY, y = rainfall_mm * 1.2),  colour = "darkgrey") +
  ylab("Daily rainfall (mm)") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(0, 100), expand = c(0,0), position = "right") + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

p4 <- ggplot() + 
  geom_line(data = model_rain_wet_sp, aes(x = DOY, y = rainfall_mm * 1.2),  colour = "darkgrey") +
  ylab("Daily rainfall (mm)") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(0, 100), expand = c(0,0), position = "right") + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

t1 <- ggplot() + 
  geom_line(data = data.frame(micro_curr_wet$shadsoil), aes(x = DOY, y = D0cm)) +
  ylab(expression("Shaded"~italic("T")["soil"]*" (°C)")) +
  ggtitle("Subtropic", subtitle = "Brisbane, Australia") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(lim = c(-5, 50), expand = c(0,0)) + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

t2 <- ggplot() + 
  geom_line(data = data.frame(micro_curr_wet_br$shadsoil), aes(x = DOY, y = D0cm)) +
  ylab(expression("Shaded"~italic("T")["soil"]*" (°C)")) +
  ggtitle("Tropic", subtitle = "Salvador, Brazil") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(-5, 50), expand = c(0,0)) + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

t3 <- ggplot() + 
  geom_line(data = data.frame(micro_curr_wet_sa$shadsoil), aes(x = DOY, y = D0cm)) +
  ylab(expression("Shaded"~italic("T")["soil"]*" (°C)")) +
  ggtitle("Semi-arid", subtitle = "Tankwa Karoo, South Africa") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(-5, 50), expand = c(0,0)) + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

t4 <- ggplot() + 
  geom_line(data = data.frame(micro_curr_wet_sp$shadsoil), aes(x = DOY, y = D0cm)) +
  ylab(expression("Shaded"~italic("T")["soil"]*" (°C)")) +
  ggtitle("Mediterranean", subtitle = "Seville, Spain") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(lim = c(-5, 50), expand = c(0,0)) + 
  coord_cartesian(xlim = c(0, 365)) +
  mytheme() +
  theme(legend.position = "bottom")

aligned_plots <- align_plots(p1, t1, align = "hv", axis = "tblr")
au_plot <- ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

aligned_plots <- align_plots(p2, t2, align = "hv", axis = "tblr")
br_plot <- ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

aligned_plots <- align_plots(p3, t3, align = "hv", axis = "tblr")
sa_plot <- ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

aligned_plots <- align_plots(p4, t4, align = "hv", axis = "tblr")
sp_plot <- ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

cowplot::plot_grid(au_plot, br_plot, sa_plot, sp_plot, ncol = 2, labels = c('a', 'b', 'c', 'd'))

Fig. S10. Differences in shaded soil temperature (T_skin, °C; black lines) and daily rainfall (mm; grey lines) that anurans are expected to experience under a represenative (a) subtropic biome (total rainfall = 1097.2 mm), (b) tropic biome (total rainfall = 1218.21 mm), (c) semi-arid biome (total rainfall = 83.37 mm), and (d) Mediterranean biome (total rainfall = 385.83 mm).

# BRAZIL MODEL
curr_wet_mod_br <- sim.ecto(micro_curr_wet_br, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

curr_dry_mod_br <- sim.ecto(micro_curr_dry_br, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

warm_wet_mod_br <- sim.ecto(micro_warm_wet_br, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

warm_dry_mod_br <- sim.ecto(micro_warm_dry_br, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

# SOUTH AFRICA MODEL
curr_wet_mod_sa <- sim.ecto(micro_curr_wet_sa, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

curr_dry_mod_sa <- sim.ecto(micro_curr_dry_sa, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

warm_wet_mod_sa <- sim.ecto(micro_warm_wet_sa, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

warm_dry_mod_sa <- sim.ecto(micro_warm_dry_sa, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

# SPAIN MODEL
curr_wet_mod_sp <- sim.ecto(micro_curr_wet_sp, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

curr_dry_mod_sp <- sim.ecto(micro_curr_dry_sp, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

warm_wet_mod_sp <- sim.ecto(micro_warm_wet_sp, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

warm_dry_mod_sp <- sim.ecto(micro_warm_dry_sp, Ww_g = Ww_g, shape = 4, Tmax = Tmax,
    Tmin = Tmin, CTmin = CTmin, CTmax = CTmax, behav = "both", in.shade = TRUE, burrow = FALSE,
    climb = FALSE, min.hyd = min_hyd, hyd.death = hyd.death, hyd.rate = hyd_rate,
    pct_wet = pct_wet_high, water = FALSE, water.act = TRUE)

# BRAZIL MODEL
curr_wet_br_df <- data.frame(curr_wet_mod_br$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = curr_wet_mod_br.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_wet = length(active[active == TRUE]))

curr_dry_br_df <- data.frame(curr_dry_mod_br$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = curr_dry_mod_br.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_dry = length(active[active == TRUE]))

warm_wet_br_df <- data.frame(warm_wet_mod_br$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = warm_wet_mod_br.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_wet = length(active[active == TRUE]))

warm_dry_br_df <- data.frame(warm_dry_mod_br$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = warm_dry_mod_br.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_dry = length(active[active == TRUE]))

br_model <- curr_wet_br_df %>%
  merge(curr_dry_br_df, by = "day") %>%
  merge(warm_wet_br_df, by = "day") %>%
  merge(warm_dry_br_df, by = "day") %>%
  pivot_longer(!day, names_to = "condition", values_to = "hours") %>%
  dplyr::mutate(condition = factor(condition, levels = c("warm_dry", "warm_wet", "curr_dry", "curr_wet")),
                season = case_when(
    day >= 335 ~ "summer",
    day >= 1 & day <= 61 ~ "summer",
    day >= 62 & day <= 153 ~ "autumn",
    day >= 154 & day <= 244 ~ "winter",
    day >= 245 & day <= 334 ~ "spring"
    ))

# SOUTH AFRICA MODEL
curr_wet_sa_df <- data.frame(curr_wet_mod_sa$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = curr_wet_mod_sa.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_wet = length(active[active == TRUE]))

curr_dry_sa_df <- data.frame(curr_dry_mod_sa$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = curr_dry_mod_sa.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_dry = length(active[active == TRUE]))

warm_wet_sa_df <- data.frame(warm_wet_mod_sa$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = warm_wet_mod_sa.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_wet = length(active[active == TRUE]))

warm_dry_sa_df <- data.frame(warm_dry_mod_sa$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = warm_dry_mod_sa.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_dry = length(active[active == TRUE]))

sa_model <- curr_wet_sa_df %>%
  merge(curr_dry_sa_df, by = "day") %>%
  merge(warm_wet_sa_df, by = "day") %>%
  merge(warm_dry_sa_df, by = "day") %>%
  pivot_longer(!day, names_to = "condition", values_to = "hours") %>%
  dplyr::mutate(condition = factor(condition, levels = c("warm_dry", "warm_wet", "curr_dry", "curr_wet")),
                season = case_when(
    day >= 335 ~ "summer",
    day >= 1 & day <= 61 ~ "summer",
    day >= 62 & day <= 153 ~ "autumn",
    day >= 154 & day <= 244 ~ "winter",
    day >= 245 & day <= 334 ~ "spring"
    ))

# SPAIN MODEL
curr_wet_sp_df <- data.frame(curr_wet_mod_sp$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = curr_wet_mod_sp.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_wet = length(active[active == TRUE]))

curr_dry_sp_df <- data.frame(curr_dry_mod_sp$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = curr_dry_mod_sp.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(curr_dry = length(active[active == TRUE]))

warm_wet_sp_df <- data.frame(warm_wet_mod_sp$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = warm_wet_mod_sp.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_wet = length(active[active == TRUE]))

warm_dry_sp_df <- data.frame(warm_dry_mod_sp$act) %>%
  tail(-8784) %>% # remove previous year 
  tibble::rowid_to_column("hour") %>%
  dplyr::rename(active = warm_dry_mod_sp.act) %>%
  dplyr::mutate(day = ceiling(1:8760/24)) %>%
  dplyr::group_by(day) %>%
  dplyr::summarise(warm_dry = length(active[active == TRUE]))

sp_model <- curr_wet_sp_df %>%
  merge(curr_dry_sp_df, by = "day") %>%
  merge(warm_wet_sp_df, by = "day") %>%
  merge(warm_dry_sp_df, by = "day") %>%
  pivot_longer(!day, names_to = "condition", values_to = "hours") %>%
  dplyr::mutate(condition = factor(condition, levels = c("warm_dry", "warm_wet", "curr_dry", "curr_wet")),
                season = case_when(
    day >= 335 ~ "winter",
    day >= 1 & day <= 61 ~ "winter",
    day >= 62 & day <= 153 ~ "spring",
    day >= 154 & day <= 244 ~ "summer",
    day >= 245 & day <= 334 ~ "autumn"
    ))

Model output - biomes

Table S7a - Summary

Table S7a. Total potential activity hours per year (\(t_{act}\)) for a 8.81 g frog under different warming (current or warming) and drought (normal or drought) conditions across represenative biomes.

data.frame(biome = c("Subtropic", "Subtropic", "Subtropic", "Subtropic",
                     "Tropic", "Tropic", "Tropic", "Tropic",
                     "Semi-arid", "Semi-arid", "Semi-arid", "Semi-arid",
                     "Mediterranean", "Mediterranean", "Mediterranean", "Mediterranean"),
           temp = c("current", "current", "warming", "warming"),
           rain = c("normal", "drought"),
           t_act_h = c(sum(shad_curr_wet_df$curr_wet), sum(shad_curr_dry_df$curr_dry),
                       sum(shad_warm_wet_df$warm_wet), sum(shad_warm_dry_df$warm_dry), 
                       sum(curr_wet_br_df$curr_wet), sum(curr_dry_br_df$curr_dry), 
                       sum(warm_wet_br_df$warm_wet), sum(warm_dry_br_df$warm_dry), 
                       sum(curr_wet_sa_df$curr_wet), sum(curr_dry_sa_df$curr_dry), 
                       sum(warm_wet_sa_df$warm_wet), sum(warm_dry_sa_df$warm_dry), 
                       sum(curr_wet_sp_df$curr_wet), sum(curr_dry_sp_df$curr_dry),
                       sum(warm_wet_sp_df$warm_wet), sum(warm_dry_sp_df$warm_dry)),
           t_act_per = c(sum(shad_curr_wet_df$curr_wet) / 8760 * 100, 
                         sum(shad_curr_dry_df$curr_dry) / 8760 * 100, 
                         sum(shad_warm_wet_df$warm_wet) / 8760 * 100, 
                         sum(shad_warm_dry_df$warm_dry) / 8760 * 100,
                       sum(curr_wet_br_df$curr_wet) / 8760 * 100, 
                       sum(curr_dry_br_df$curr_dry) / 8760 * 100, 
                       sum(warm_wet_br_df$warm_wet) / 8760 * 100, 
                       sum(warm_dry_br_df$warm_dry) / 8760 * 100,
                       sum(curr_wet_sa_df$curr_wet) / 8760 * 100, 
                       sum(curr_dry_sa_df$curr_dry) / 8760 * 100, 
                       sum(warm_wet_sa_df$warm_wet) / 8760 * 100, 
                       sum(warm_dry_sa_df$warm_dry) / 8760 * 100,
                       sum(curr_wet_sp_df$curr_wet) / 8760 * 100, 
                       sum(curr_dry_sp_df$curr_dry) / 8760 * 100, 
                       sum(warm_wet_sp_df$warm_wet) / 8760 * 100, 
                       sum(warm_dry_sp_df$warm_dry) / 8760 * 100)) %>%
  knitr::kable(col.names = c("Biome", "Warming simulation", "Drought simulation", "$t_{act}$ (h)", "$t_{act}$ (%)"))

Biome	Warming simulation	Drought simulation	\(t_{act}\) (h)	\(t_{act}\) (%)
Subtropic	current	normal	3934	44.91
Subtropic	current	drought	3729	42.57
Subtropic	warming	normal	4268	48.72
Subtropic	warming	drought	4022	45.91
Tropic	current	normal	5561	63.48
Tropic	current	drought	4355	49.71
Tropic	warming	normal	5373	61.34
Tropic	warming	drought	4113	46.95
Semi-arid	current	normal	42	0.48
Semi-arid	current	drought	43	0.49
Semi-arid	warming	normal	79	0.90
Semi-arid	warming	drought	77	0.88
Mediterranean	current	normal	372	4.25
Mediterranean	current	drought	277	3.16
Mediterranean	warming	normal	907	10.35
Mediterranean	warming	drought	585	6.68

---

Table S7b - Relative change year

Table S7b. Change in \(t_{act}\) (%) f to the current normal scenario under warming only, drought only, and warming and drought combined for a hypothetical frog in different biomes.

data.frame(biome = c("Subtropic","Tropic", "Semi-arid", "Mediterranean"),
           curr_wet = c(sum(shad_curr_wet_df$curr_wet), sum(curr_wet_br_df$curr_wet), sum(curr_wet_sa_df$curr_wet), sum(curr_wet_sp_df$curr_wet)),
           curr_dry = c(sum(shad_curr_dry_df$curr_dry), sum(curr_dry_br_df$curr_dry), sum(curr_dry_sa_df$curr_dry), sum(curr_dry_sp_df$curr_dry)),
           warm_wet = c(sum(shad_warm_wet_df$warm_wet), sum(warm_wet_br_df$warm_wet), sum(warm_wet_sa_df$warm_wet), sum(warm_wet_sp_df$warm_wet)),
           warm_dry = c(sum(shad_warm_dry_df$warm_dry), sum(warm_dry_br_df$warm_dry), sum(warm_dry_sa_df$warm_dry), sum(warm_dry_sp_df$warm_dry))) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  knitr::kable(col.names = c("Biome", "$\\Delta$ warming only (%)", "$\\Delta$ drought only (%)", "$\\Delta$ warming and drought (%)"))

Biome	\(\Delta\) warming only (%)	\(\Delta\) drought only (%)	\(\Delta\) warming and drought (%)
Subtropic	8.5	-5.2	2.2
Tropic	-3.4	-21.7	-26.0
Semi-arid	88.1	2.4	83.3
Mediterranean	143.8	-25.5	57.3

---

Table S7c - Relative change season

Table S7c. Change in \(t_{act}\) (%) relative to the current normal scenario under warming only, drought only, and warming and drought combined for a hypothetical frog in different biomes.

# summer
shad_model_summer <- shad_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

br_model_summer <- br_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sa_model_summer <- sa_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sp_model_summer <- sp_model %>%
  dplyr::filter(season == "summer") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_summer <- data.frame(biome = c("Subtropic","Tropic", "Semi-arid", "Mediterranean"),
           curr_wet = c(shad_model_summer$hours[4], br_model_summer$hours[4], sa_model_summer$hours[4], sp_model_summer$hours[4]), 
           curr_dry = c(shad_model_summer$hours[3], br_model_summer$hours[3], sa_model_summer$hours[3], sp_model_summer$hours[3]),
           warm_wet = c(shad_model_summer$hours[2], br_model_summer$hours[2], sa_model_summer$hours[2], sp_model_summer$hours[2]),
           warm_dry = c(shad_model_summer$hours[1], br_model_summer$hours[1], sa_model_summer$hours[1], sp_model_summer$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(biome = "**Summer**", .before = 1)

# autumn
shad_model_autumn <- shad_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

br_model_autumn <- br_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sa_model_autumn <- sa_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sp_model_autumn <- sp_model %>%
  dplyr::filter(season == "autumn") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_autumn <- data.frame(biome = c("Subtropic","Tropic", "Semi-arid", "Mediterranean"),
           curr_wet = c(shad_model_autumn$hours[4], br_model_autumn$hours[4], sa_model_autumn$hours[4], sp_model_autumn$hours[4]), 
           curr_dry = c(shad_model_autumn$hours[3], br_model_autumn$hours[3], sa_model_autumn$hours[3], sp_model_autumn$hours[3]),
           warm_wet = c(shad_model_autumn$hours[2], br_model_autumn$hours[2], sa_model_autumn$hours[2], sp_model_autumn$hours[2]),
           warm_dry = c(shad_model_autumn$hours[1], br_model_autumn$hours[1], sa_model_autumn$hours[1], sp_model_autumn$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(biome = "**Autumn**", .before = 1)

# winter
shad_model_winter <- shad_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

br_model_winter <- br_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sa_model_winter <- sa_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sp_model_winter <- sp_model %>%
  dplyr::filter(season == "winter") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_winter <- data.frame(biome = c("Subtropic","Tropic", "Semi-arid", "Mediterranean"),
           curr_wet = c(shad_model_winter$hours[4], br_model_winter$hours[4], sa_model_winter$hours[4], sp_model_winter$hours[4]), 
           curr_dry = c(shad_model_winter$hours[3], br_model_winter$hours[3], sa_model_winter$hours[3], sp_model_winter$hours[3]),
           warm_wet = c(shad_model_winter$hours[2], br_model_winter$hours[2], sa_model_winter$hours[2], sp_model_winter$hours[2]),
           warm_dry = c(shad_model_winter$hours[1], br_model_winter$hours[1], sa_model_winter$hours[1], sp_model_winter$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(biome = "**Winter**", .before = 1)

# spring
shad_model_spring <- shad_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

br_model_spring <- br_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sa_model_spring <- sa_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

sp_model_spring <- sp_model %>%
  dplyr::filter(season == "spring") %>%
  dplyr::group_by(condition) %>%
  dplyr::summarise(hours = sum(hours))

tab_spring <- data.frame(biome = c("Subtropic","Tropic", "Semi-arid", "Mediterranean"),
           curr_wet = c(shad_model_spring$hours[4], br_model_spring$hours[4], sa_model_spring$hours[4], sp_model_spring$hours[4]), 
           curr_dry = c(shad_model_spring$hours[3], br_model_spring$hours[3], sa_model_spring$hours[3], sp_model_spring$hours[3]),
           warm_wet = c(shad_model_spring$hours[2], br_model_spring$hours[2], sa_model_spring$hours[2], sp_model_spring$hours[2]),
           warm_dry = c(shad_model_spring$hours[1], br_model_spring$hours[1], sa_model_spring$hours[1], sp_model_spring$hours[1])) %>%
  dplyr::mutate(delta_warm     = (warm_wet - curr_wet) / curr_wet * 100,
                delta_dry      = (curr_dry - curr_wet) / curr_wet * 100,
                delta_warm_dry = (warm_dry - curr_wet) / curr_wet * 100) %>%
  dplyr::select(-c(curr_wet:warm_dry)) %>%
  tibble::add_row(biome = "**Spring**", .before = 1)

# Render table
bind_rows(tab_summer, tab_autumn, tab_winter, tab_spring) %>% 
  remove_rownames() %>% 
  knitr::kable(col.names = c("Biome", "$\\Delta$ warming only (%)", "$\\Delta$ drought only (%)", "$\\Delta$ warming and drought (%)"))

Biome	\(\Delta\) warming only (%)	\(\Delta\) drought only (%)	\(\Delta\) warming and drought (%)
Summer
Subtropic	-8.39	-13.9	-21.8
Tropic	-2.06	-13.7	-12.9
Semi-arid	-9.09	9.1	0.0
Mediterranean	-7.41	-3.7	-11.1
Autumn
Subtropic	4.75	-1.8	3.3
Tropic	-3.93	-20.6	-25.9
Semi-arid	46.67	0.0	46.7
Mediterranean	-0.81	-35.5	-28.2
Winter
Subtropic	72.88	-1.9	65.5
Tropic	-2.37	-21.6	-30.9
Semi-arid	1300.00	0.0	1250.0
Mediterranean	601.47	-16.2	301.5
Spring
Subtropic	0.77	-2.8	-4.4
Tropic	-5.02	-28.6	-29.7
Semi-arid	35.71	0.0	21.4
Mediterranean	84.31	-25.5	30.1

---

Figures

Code to produce the main document figures are detailed below. Figures produced were further modified in Adobe Illustrator for publication.

Figure 1 - AI risk

Create plots for main text figure 1.

# Reproject maps
rob_proj   <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" 
world_rob  <- spTransform(world_spdf, CRSobj = rob_proj)

# land
data("land", package = "tmap")
land_raster <- as(land, "Raster") # convert from star to raster
land_df <- raster::as.data.frame(raster::rasterToPoints(land_raster))

elevation_raster <- projectRaster(land_raster$elevation, crs = rob_proj)
slope_raster     <- raster::terrain(elevation_raster, opt = 'slope')
aspect_raster    <- raster::terrain(elevation_raster, opt = 'aspect')
hill_raster      <- hillShade(slope_raster, aspect_raster, 40, 270) #setting the elevation angle (of the sun) to 40 and the direction angle of the light to 270
hill_m           <- rasterToPoints(hill_raster)
hill_df          <-  data.frame(hill_m)
colnames(hill_df) <- c("lon", "lat", "hill")

# Fig 1a - AI current
AI_rob    <- raster::projectRaster(ai_rast, crs = rob_proj)
AI_df_rob <- raster::as.data.frame(raster::rasterToPoints(AI_rob)) 
AI_df_rob <- AI_df_rob %>% 
  # Recode
  dplyr::mutate(category = case_when(
    is.infinite(layer) ~ 'Humid',
    layer >= 0.65 ~ 'Humid',                       
    layer >= 0.5 & layer < 0.65 ~ 'Dry sub-humid', 
    layer >= 0.2 & layer < 0.5 ~ 'Semi-arid',      
    layer >= 0.05 & layer < 0.2 ~ 'Arid',         
    layer < 0.05 ~ 'Hyper-arid'                    
  )) %>% 
  # Convert to ordered factor
  dplyr::mutate(category = factor(category,
                                  levels = c('Hyper-arid', 'Arid', 'Semi-arid',
                                             'Dry sub-humid', 'Humid'),
                                  ordered = TRUE))

AI_2C_rob    <- projectRaster(ai_2C_rast, crs = rob_proj)
AI_2C_df_rob <- raster::as.data.frame(raster::rasterToPoints(AI_2C_rob)) 
AI_2C_df_rob <- AI_2C_df_rob %>% 
  # Recode
  dplyr::mutate(category = case_when(
    is.infinite(layer) ~ 'Humid',
    layer >= 0.65 ~ 'Humid',                       
    layer >= 0.5 & layer < 0.65 ~ 'Dry sub-humid', 
    layer >= 0.2 & layer < 0.5 ~ 'Semi-arid',     
    layer >= 0.05 & layer < 0.2 ~ 'Arid',         
    layer < 0.05 ~ 'Hyper-arid'                
  )) %>% 
  # Convert to ordered factor
  dplyr::mutate(category = factor(category,
                                  levels = c('Hyper-arid', 'Arid', 'Semi-arid',
                                             'Dry sub-humid', 'Humid'),
                                  ordered = TRUE))

AI_4C_rob    <- projectRaster(ai_4C_rast, crs = rob_proj)
AI_4C_df_rob <- raster::as.data.frame(raster::rasterToPoints(AI_4C_rob)) 
AI_4C_df_rob <- AI_4C_df_rob %>% 
  # Recode
  dplyr::mutate(category = case_when(
    is.infinite(layer) ~ 'Humid',
    layer >= 0.65 ~ 'Humid',                       
    layer >= 0.5 & layer < 0.65 ~ 'Dry sub-humid', 
    layer >= 0.2 & layer < 0.5 ~ 'Semi-arid',      
    layer >= 0.05 & layer < 0.2 ~ 'Arid',         
    layer < 0.05 ~ 'Hyper-arid'                    
  )) %>% 
  # Convert to ordered factor
  dplyr::mutate(category = factor(category,
                                  levels = c('Hyper-arid', 'Arid', 'Semi-arid',
                                             'Dry sub-humid', 'Humid'),
                                  ordered = TRUE))

arid_col <- c('#8E063B', '#CB6D53', '#E99A2C', '#F5D579', 'white')
aridity_plot <- ggplot() +
  geom_raster(data = AI_df_rob, aes(y = y, x = x, fill = category)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = arid_col,
                    guide = guide_legend(reverse = TRUE)) +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Aridity Index (1981-2010)") +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1)  

# Fig 1b - Species richness
anuran_rob    <- projectRaster(anuran_sr, crs = rob_proj)
anuran_df_rob <- raster::as.data.frame(raster::rasterToPoints(anuran_rob)) %>% dplyr::rename(species_n = layer)

sp_breaks = c(1, 5, 25, 100)
anuran_plot <- ggplot() +
  geom_raster(data = anuran_df_rob, aes(y = y, x = x, fill = species_n)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  colorspace::scale_fill_continuous_sequential(palette = "BluGrn", trans = "log",
                                               breaks = sp_breaks, labels = sp_breaks) +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle("Anuran species richness") +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed()

AI_merged <- merge(AI_df_rob, AI_2C_df_rob, by = c("x", "y"), all.x = T) %>%
  merge(AI_4C_df_rob, by = c("x", "y"), all.x = T) %>%
  dplyr::rename(layer_current = layer.x,
                category_current = category.x,
                layer_2C  = layer.y,
                category_2C = category.y,
                layer_4C  = layer,
                category_4C = category) %>%
  dplyr::mutate(change_2C = (layer_2C - layer_current / layer_current) * 100,
                change_2C_AI = case_when(
                  change_2C >= -2 & change_2C <2 ~ 'No change to wetter',
                  change_2C >= -10 & change_2C < -2 ~ '>-2 to -10%',
                  change_2C >= -20 & change_2C < -10 ~ '>-10 to -20%',
                  change_2C >= -40 & change_2C < -20 ~ '>-20 to -40%',
                  change_2C >= -80 & change_2C < -40 ~ '>-40 to -80%',
                  change_2C < -80 ~ '<-80%'),
                  change_4C = (layer_4C - layer_current / layer_current) * 100,
                change_4C_AI = case_when(
                  change_4C >= -2 & change_4C <2 ~ 'No change to wetter',
                  change_4C >= -10 & change_4C < -2 ~ '>-2 to -10%',
                  change_4C >= -20 & change_4C < -10 ~ '>-10 to -20%',
                  change_4C >= -40 & change_4C < -20 ~ '>-20 to -40%',
                  change_4C >= -80 & change_4C < -40 ~ '>-40 to -80%',
                  change_4C < -80 ~ '<-80%'
                )) %>% 
  filter(category_current != "Hyper-arid") %>%
  dplyr::mutate(change_2C_AI = factor(change_2C_AI,
                                   levels = c('<-80%', '>-40 to -80%','>-20 to -40%',
                                              '>-10 to -20%', '>-2 to -10%', 'No change to wetter'),
                                   ordered = TRUE),
                change_4C_AI = factor(change_4C_AI,
                                   levels = c('<-80%', '>-40 to -80%','>-20 to -40%',
                                              '>-10 to -20%', '>-2 to -10%', 'No change to wetter'),
                                   ordered = TRUE))

change_col <- c('#9e5b13', '#bf832e', '#dab367', '#e8ca84', '#f0e0ad', 'white')

# Fig 1C - 2C Difference
change_2C_plot <- ggplot() +
  geom_raster(data = AI_merged, aes(y = y, x = x, fill = change_2C_AI)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = change_col, na.value = "white") +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void()+ ylab(NULL) + xlab(NULL) +
  ggtitle(expression("Change in dryness ("*Delta*"AI +2°C)")) +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) 

# Fig 1d - 4C Difference
change_4C_plot <- ggplot() +
  geom_raster(data = AI_merged, aes(y = y, x = x, fill = change_4C_AI)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_manual(values = change_col, na.value = "white") +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void()+ ylab(NULL) + xlab(NULL) +
  ggtitle(expression("Change in dryness ("*Delta*"AI +4°C)")) +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) 

# Fig 1e - Dryness and species richness
arid_col_2 <- c("#8E063B", "#CB6D53", "#E99A2C", "#F5D579", "light grey")
species_ai_plot <- ai_sp_df %>% 
  ggplot(aes(x = aridity, y = species_n, colour = category)) +
  geom_point(size = 2) +
  coord_trans(x = "sqrt") +
  scale_colour_manual(values = arid_col_2,
                    guide = guide_legend(reverse = TRUE)) +
  scale_x_continuous(limits = c(0,3),
                     breaks = c(0, 0.05, 0.2, 0.5, 0.65)) +
  xlab("Precipitation / Evapotranspiration") +
  ylab("Species richness") +
  mytheme() + theme(legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank())

# Fig 1f - Dryness and species richness by ecotype
ecotype_ai_plot <- ai_sp_df %>%
  dplyr::filter(aridity <=3) %>%
  dplyr::select(aridity, species_n:stream_n) %>%
  tidyr::pivot_longer(!aridity, names_to = "ecotype", values_to = "mean") %>%
  dplyr::filter(ecotype != "species_n") %>%
  drop_na(mean) %>%
  dplyr::mutate(ecotype = factor(ecotype, 
                                 levels = c("fossorial_n", "ground_n", "aquatic_n", 
                                            "arboreal_n", "semi_aq_n", "stream_n"))) %>%
  ggplot(aes(x = aridity, y = ecotype, group = ecotype, fill = ecotype)) +
  ggridges::stat_density_ridges(scale = 1, rel_min_height = 0.01, alpha = 0.5,
                                show.legend = FALSE) +
  scale_fill_viridis_d() +
  scale_x_continuous(trans = "sqrt") +
  labs(x = NULL, y = NULL) +
  mytheme()

# Fig 1g - Change in species AI
# Merge summary data by row
arid_sp_comb          <- data.frame(bind_rows(ai_sp_sum, ai_2C_sp_sum, ai_4C_sp_sum))
arid_sp_comb$scenario <- forcats::fct_relevel(arid_sp_comb$scenario, "Current")

species_occ_plot <- arid_sp_comb %>% 
  ggplot(aes(x = category, y = log(all_freq), group = scenario, colour = scenario, shape = scenario)) +
  geom_point(position = position_dodge(width = .5), size = 1.5) +
  scale_colour_manual(values = c("#F5D579", "#CB6D53", "#8E063B")) +
  scale_y_continuous(breaks = pretty(log(arid_sp_comb$all_freq)), labels = pretty(arid_sp_comb$all_freq)) +
  ylab("% species") + xlab(NULL) +
  mytheme() + theme(legend.position = "none")

# Fig 1
plot_grid(species_ai_plot, ecotype_ai_plot, species_occ_plot, align = "v", ncol = 3) # plot in between.

cowplot::plot_grid(aridity_plot + theme(legend.position = "bottom"), 
                   anuran_plot + theme(legend.position = "bottom"), 
                   change_2C_plot + theme(legend.position = "bottom"),
                   change_4C_plot + theme(legend.position = "bottom"),
                   align = "h", axis = "bt",
                   ncol = 2)

Figure 2 - PDSI risk

Create plots for main text figure 2.

PDSI_risk_rob <- raster::projectRaster(raster::stack(PDSI_2C_diff_rast, PDSI_4C_diff_rast, PDSI_freq_diff, PDSI_dur_diff,
                                                     resample(anuran_sr, PDSI_4C_diff_rast)), 
                                       crs = rob_proj)
                             
PDSI_risk_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_risk_rob)) %>%
  dplyr::rename("sp_n" = layer) %>%
  dplyr::filter(!is.na(sp_n)) %>%
  dplyr::mutate(delta_int_2C_bin = ifelse(delta_int_2C < -1, 1, 0), # delta_PDSI < -1
         delta_int_4C_bin = ifelse(delta_int_4C < -1, 1, 0), # delta_PDSI < -1
         delta_freq_2C_bin = ifelse(delta_freq_2C >1, 1, 0), # month > 1
         delta_freq_4C_bin = ifelse(delta_freq_4C >1, 1, 0), # month > 1
         delta_dur_2C_bin = ifelse(delta_dur_2C >1, 1, 0), # month > 1
         delta_dur_4C_bin = ifelse(delta_dur_4C >1, 1, 0)) %>% # month > 1 
  dplyr::rowwise() %>%
  dplyr::mutate(count_2C = sum(delta_int_2C_bin, delta_freq_2C_bin, delta_dur_2C_bin, na.rm = T),
         
         count_4C = sum(delta_int_4C_bin, delta_freq_4C_bin, delta_dur_4C_bin, na.rm = T),
         risk_2C  = sp_n * count_2C,
         risk_4C  = sp_n * count_4C,
         count_2C = factor(count_2C, levels = c("3", "2", "1")),
         count_4C = factor(count_4C, levels = c("3", "2", "1"))
         ) %>%
  data.frame()

# Plot
PDSI_risk_2C_plot <- ggplot() +
  geom_raster(data = PDSI_risk_df, aes(y = y, x = x, fill = count_2C)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_discrete_sequential(palette = "RedOr", rev = F) +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle(expression("+2°C by 2080-2100")) +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

PDSI_risk_4C_plot <- ggplot() +
  geom_raster(data = PDSI_risk_df, aes(y = y, x = x, fill = count_4C)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_discrete_sequential(palette = "RedOr", rev = F) +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  ggtitle(expression("+4°C by 2080-2100")) +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

sp_2C_plot <- PDSI_risk_df %>%
  dplyr::mutate(risk_2C_cat = case_when(
    risk_2C == 0 ~ '0',
    risk_2C > 0 & risk_2C <=1 ~ '<1',
    risk_2C > 1 & risk_2C < 5 ~ '>1 to 5',
    risk_2C >= 5 & risk_2C < 10 ~ '>5 to 10',
    risk_2C >= 10 & risk_2C < 20 ~ '>10 to 20',
    risk_2C >= 20 & risk_2C < 100 ~ '>20 to 100',
    risk_2C >= 100 & risk_2C < 200 ~ '>100 to 200',
    risk_2C > 200 ~ '>200'),
    risk_2C = na_if(risk_2C, 0),
    risk_2C_cat = factor(risk_2C_cat,
                                   levels = c('>200','>100 to 200', '>20 to 100','>10 to 20', '>5 to 10', '>1 to 5', '<1'))) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = risk_2C_cat)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_discrete_sequential(palette = "Heat", rev = F) +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

sp_4C_plot <- PDSI_risk_df %>%
  dplyr::mutate(risk_4C_cat = case_when(
    risk_4C == 0 ~ '0',
    risk_4C > 0 & risk_4C <=1 ~ '<1',
    risk_4C > 1 & risk_4C < 5 ~ '>1 to 5',
    risk_4C >= 5 & risk_4C < 10 ~ '>5 to 10',
    risk_4C >= 10 & risk_4C < 20 ~ '>10 to 20',
    risk_4C >= 20 & risk_4C < 100 ~ '>20 to 100',
    risk_4C >= 100 & risk_4C < 200 ~ '>100 to 200',
    risk_4C > 200 ~ '>200'),
    risk_4C = na_if(risk_4C, 0),
    risk_4C_cat = factor(risk_4C_cat,
                                   levels = c('>200','>100 to 200', '>20 to 100','>10 to 20', '>5 to 10', '>1 to 5', '<1'))) %>%
  ggplot() +
  geom_raster(aes(y = y, x = x, fill = risk_4C_cat)) +
  geom_polygon(data = world_rob, aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
  scale_fill_discrete_sequential(palette = "Heat", rev = F) +
  ggnewscale::new_scale("fill") +
  geom_raster(data = hill_df %>% filter(hill <= 0.645), aes(lon, lat, fill = hill, alpha = hill), show.legend = FALSE) +
  scale_fill_gradient(low = "black", high = "grey") +
  scale_alpha_continuous(trans = "reverse") +
  theme_void() + ylab(NULL) + xlab(NULL) +
  scale_y_continuous(limits = c(-61e5, 85e5), expand = c(0, 0)) +
  scale_x_continuous(limits = c(-15e6, 16e6), expand = c(0, 0)) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 10)) +
  coord_fixed(ratio = 1) # fixed ratio

cowplot::plot_grid(PDSI_risk_2C_plot + theme(legend.position = "none"),
                   PDSI_risk_4C_plot + theme(legend.position = "none"),
                   sp_2C_plot + theme(legend.position = "none"),
                   sp_4C_plot + theme(legend.position = "none"),
                   ncol = 2,
                   align = "h", axis = "bt", labels = c('a', 'b', 'c', 'd'))

#PDSI_risk_df %>%
  #group_by(count_2C) %>%
  #summarise(length = length(count_2C)) %>%
  #mutate(freq = length / length(PDSI_risk_df$count_2C) * 100)

#PDSI_risk_df %>%
  #group_by(count_4C) %>%
  #summarise(length = length(count_4C)) %>%
  #mutate(freq = length / length(PDSI_risk_df$count_4C) * 100)

Figure 3 - Water loss risk

Run the ectotherm function from NicheMapR for all 12 months from the Terra Climate raster files containing the mean wind speed (m s^-1), maximum air temperature (°C), and vapour pressure deficit (VPD; kPa) for the current (1981 to 2010), +2°C and +4°C scenario. From the ecotherm model, we extracted the EWL (g h^-1) for each raster pixel to determine the spatial variation in EWL.

## Current ## ---------------------------------------
elev      <- getData("worldclim", var = "alt", res = 2.5)
ws_full   <-  raster::brick(file.path(spatial_path, 'TerraClimate19812010_ws.nc')) # 10m, m/s
tmax_full <- raster::brick(file.path(spatial_path, 'TerraClimate19812010_tmax.nc')) # 2m, deg C
vpd_full  <- raster::brick(file.path(spatial_path, 'TerraClimate19812010_vpd.nc')) # 2m, kpa
Tbs       <- crop(tmax_full, extent(elev))
EWLs      <- Tbs 

# Run ectotherm model for each month and save on local drive (current scenario)
for(month in 1:12){
  # Temperature
  Ta.w   <- crop(tmax_full[[month]], extent(elev))
  esat   <- VAPPRS(Ta.w)
  # RH
  e      <- esat - crop(vpd_full[[month]], extent(elev)) * 1000
  RH.w   <- e / esat * 100
  RH.w[RH.w < 0] <- 0
  # Wind speed
  Wind.w <- crop(ws_full[[month]], extent(elev))
  Vstar  <- 0.4 * Wind.w / log((100 / 0.15) + 1) # friction velocity
  Wind.w <- 2.5 * Vstar * log((1 / 0.15) + 1) # correct for reference height (10 m) to ground level (1 cm) assuming level ground
  # compute atmospheric pressure (kPa)
  P.w <- 101325 * ((1 - (0.0065 * elev / 288)) ^ (1 / 0.190284)) / 1000
  
  # Get raster values
  Ta  <- getValues(Ta.w)
  RH  <- getValues(RH.w)
  P   <- getValues(P.w)
  Pa  <- getValues(e)
  ws  <- getValues(Wind.w)
  P[!is.na(Ta) & is.na(P)] <- 101.325
  Ta2 <- Ta[!is.na(Ta)]
  RH2 <- RH[!is.na(Ta)]
  P2  <- P[!is.na(Ta)]
  ws2 <- ws[!is.na(Ta)]
  Pa2 <- Pa[!is.na(Ta)]
  
  div <- 10
  Twb <- Ta2[1:floor(length(Ta2) / div)]
  micro_base <- micro_global()
  
  for(k in 1:div){
    micro <- micro_base
    if(k == 1){
      extra <- round(24 - (length(Twb) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb, rep(Ta2[1], extra))
      start <- 1
      end   <- length(Twb)
    }
    if(k > 1 & k < div){
      fact  <- (length(Twb) * (k - 1) + k - 1)
      Twb4  <- Ta2[fact:(fact + length(Twb))]
      extra <- round(24 - (length(Twb4) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb4, rep(Ta2[1], extra))
      start <- fact
      end   <- fact + length(Twb)
    }        
    if(k == div){
      fact  <- (length(Twb) * (k - 1) + k - 1)
      Twb4  <- Ta2[fact:length(Ta2)]
      extra <- round(24 - (length(Twb4) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb4, rep(Ta2[1], extra))
      start <- fact
      end   <- length(Ta2)
    } 
    
    soilnames     <- dimnames(micro$soil)
    metoutnames   <- dimnames(micro$metout)
    tcondnames    <- dimnames(micro$tcond)
    specheatnames <- dimnames(micro$specheat)
    densitnames   <- dimnames(micro$densit)
    plantnames    <- dimnames(micro$plant)
    
    micro$metout   <- t(array(micro$metout[1,], dim = c(ncol(micro$metout), length(Twb4))))
    dimnames(micro$metout) <- metoutnames
    micro$soil     <- t(array(micro$soil[1,], dim = c(ncol(micro$soil), length(Twb4))))
    dimnames(micro$soil) <- soilnames
    micro$tcond    <- t(array(micro$tcond[1,], dim = c(ncol(micro$tcond), length(Twb4))))
    dimnames(micro$tcond) <- tcondnames
    micro$specheat <- t(array(micro$specheat[1,], dim = c(ncol(micro$specheat), length(Twb4))))
    dimnames(micro$specheat) <- specheatnames  
    micro$densit   <- t(array(micro$densit[1,], dim = c(ncol(micro$densit), length(Twb4))))
    dimnames(micro$densit) <- densitnames      
    micro$plant    <- t(array(micro$plant[1,], dim = c(ncol(micro$plant), length(Twb4))))
    dimnames(micro$plant) <- plantnames 
    
    micro$shadmet <- micro$metout
    micro$shadsoil <- micro$soil
    micro$humid <- micro$soil
    micro$shadhumid <- micro$soil
    micro$soilpot <- micro$soil
    micro$shadpot <- micro$soil
    micro$soilmoist <- micro$soil
    micro$shadmoist <- micro$soil
    micro$shadtcond <- micro$tcond
    micro$shaddensit <- micro$densit
    micro$shadspecheat <- micro$specheat
    micro$shadplant <- micro$plant
    micro$RAINFALL <- rep(0, length(Twb4)/24)
    micro$minshade <- rep(0, length(Twb4)/24)
    micro$maxshade <- rep(90, length(Twb4)/24)
    micro$nyears <- ceiling(length(Twb4)/24/365)
    
    # replace first line of microclimate input with air temperature and humidity at which to compute Tw
    micro$metout[1:length(Twb4),c(3,4,14)] <- c(Ta2[start:end], rep(Ta2[1], extra)) # set air and sky temperature
    micro$metout[1:length(Twb4),c(5,6)] <- c(RH2[start:end], rep(RH2[1], extra)) # set relative humidity
    # micro$metout[1:length(Twb3),c(7,8)] <- 0.05 # set wind speed to low (sensitive to this)
    micro$soil[1:length(Twb4),3] <- c(Ta2[start:end], rep(Ta2[1], extra)) # set surface soil temperature to air temperature
    micro$metout[1:length(Twb4), 2] <- rep(seq(0, 23 * 60, 60), length(Ta2)/24 + 1)[1:length(Twb4)]
    #micro$metout[1:length(Twb3), 1] <- seq(0, 23 * 60, 60)[1:length(Twb3)]
    preshr <- c(P2[start:end], rep(P2[1], extra))*1000          
    
    micro$metout[1:length(Twb4),c(7,8)] <- c(ws2[start:end], rep(ws2[1], extra)) # set wind speed    
    
    # run the ectotherm model 
    message(paste('running ectotherm simulation for month ',month,' for ', length(Twb4),' sites \n'))
    ptm <- proc.time() # Start timing
    gc()
    ecto <- ectotherm(Ww_g = 8.7,
                      live = 0,
                      pct_wet = 87,
                      shape = 4, # frog
                      preshr = preshr)
    message(paste0('runtime ', (proc.time() - ptm)[3], ' seconds')) # Stop the clock
    gc()
    environ <- as.data.frame(ecto$environ)
    masbal  <- as.data.frame(ecto$masbal)
    if(k == 1){
      Tb.ectotherm  <- environ$TC[1:(nrow(environ) - extra)]
      EWL.ectotherm <- masbal$H2OCut_g[1:(nrow(masbal) - extra)]
    }else{
      Tb.ectotherm  <- c(Tb.ectotherm, environ$TC[1:(nrow(environ) - extra)])
      EWL.ectotherm <- c(EWL.ectotherm, masbal$H2OCut_g[1:(nrow(masbal) - extra)])
    }
  }
  Tb2     <- Ta
  Tb2[!is.na(Ta)] <- Tb.ectotherm
  Tb.grid <- Ta.w
  Tb.grid <- setValues(Tb.grid, Tb2)
  #plot(Tb.grid)
  #map('world', add = TRUE)
  
  EWL2 <- Ta
  EWL2[!is.na(Ta)] <- EWL.ectotherm
  EWL.grid <- Ta.w
  EWL.grid <- setValues(EWL.grid, EWL2)
  #plot(EWL.grid)
  EWL.grid[Tb.grid <= 0] <- NA
  #map('world', add = TRUE)

  writeRaster(EWL.grid, file = paste0('EWL_',month,'.nc'), overwrite = TRUE)
} 

## +2C ## ---------------------------------------
tmax_full_2C <- raster::brick(file.path(spatial_path, 'TerraClimate2C_tmax.nc')) # 2m, deg C
es_kPa_2C    <- tmax_full_2C
values(es_kPa_2C) <- 0.611 * exp(2500000 / 461.5 * ( 1 / 273 - 1 / (values(tmax_full_2C) + 273.15)))

# Actual vapour pressure
ea_kPa <- raster::brick(file.path(spatial_path, 'TerraClimate19812010_vap.nc'))
ea_kPa_2C <- ea_kPa
values(ea_kPa_2C) <- values(ea_kPa_2C) * 0.9

# Calculate VPD (kPa)
vpd_full_2C <- ea_kPa_2C
values(vpd_full_2C) <- values(es_kPa_2C) - values(ea_kPa_2C)

Tbs <- crop(tmax_full_2C, extent(elev))
EWLs <- Tbs 

# Run ectotherm model for each month and save on local drive (+2C scenario)
for(month in 1:12){
  Ta.w   <- crop(tmax_full_2C[[month]], extent(elev))
  esat   <- VAPPRS(Ta.w)
  e      <- esat - crop(vpd_full_2C[[month]], extent(elev)) * 1000
  RH.w   <- e / esat * 100
  RH.w[RH.w < 0] <- 0
  Wind.w <- crop(ws_full[[month]], extent(elev))
  Vstar  <- 0.4 * Wind.w / log((100 / 0.15) + 1) # friction velocity
  Wind.w <- 2.5 * Vstar * log((1 / 0.15) + 1) # correct for reference height (10 m) to ground level (1 cm) assuming level ground
  # compute atmospheric pressure (kPa)
  P.w <- 101325 * ((1 - (0.0065 * elev / 288)) ^ (1 / 0.190284)) / 1000
  
  Ta  <- getValues(Ta.w)
  RH  <- getValues(RH.w)
  P   <- getValues(P.w)
  Pa  <- getValues(e)
  ws  <- getValues(Wind.w)
  P[!is.na(Ta) & is.na(P)] <- 101.325
  Ta2 <- Ta[!is.na(Ta)]
  RH2 <- RH[!is.na(Ta)]
  P2  <- P[!is.na(Ta)]
  ws2 <- ws[!is.na(Ta)]
  Pa2 <- Pa[!is.na(Ta)]
  
  div <- 10
  Twb <- Ta2[1:floor(length(Ta2) / div)]
  micro_base <- micro_global()
  
  for(k in 1:div){
    micro <- micro_base
    if(k == 1){
      extra <- round(24 - (length(Twb) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb, rep(Ta2[1], extra))
      start <- 1
      end   <- length(Twb)
    }
    if(k > 1 & k < div){
      fact  <- (length(Twb) * (k - 1) + k - 1)
      Twb4  <- Ta2[fact:(fact + length(Twb))]
      extra <- round(24 - (length(Twb4) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb4, rep(Ta2[1], extra))
      start <- fact
      end   <- fact + length(Twb)
    }        
    if(k == div){
      fact  <- (length(Twb) * (k - 1) + k - 1)
      Twb4  <- Ta2[fact:length(Ta2)]
      extra <- round(24 - (length(Twb4) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb4, rep(Ta2[1], extra))
      start <- fact
      end   <- length(Ta2)
    } 
    
    soilnames     <- dimnames(micro$soil)
    metoutnames   <- dimnames(micro$metout)
    tcondnames    <- dimnames(micro$tcond)
    specheatnames <- dimnames(micro$specheat)
    densitnames   <- dimnames(micro$densit)
    plantnames    <- dimnames(micro$plant)
    
    micro$metout   <- t(array(micro$metout[1,], dim = c(ncol(micro$metout), length(Twb4))))
    dimnames(micro$metout) <- metoutnames
    micro$soil     <- t(array(micro$soil[1,], dim = c(ncol(micro$soil), length(Twb4))))
    dimnames(micro$soil) <- soilnames
    micro$tcond    <- t(array(micro$tcond[1,], dim = c(ncol(micro$tcond), length(Twb4))))
    dimnames(micro$tcond) <- tcondnames
    micro$specheat <- t(array(micro$specheat[1,], dim = c(ncol(micro$specheat), length(Twb4))))
    dimnames(micro$specheat) <- specheatnames  
    micro$densit   <- t(array(micro$densit[1,], dim = c(ncol(micro$densit), length(Twb4))))
    dimnames(micro$densit) <- densitnames      
    micro$plant    <- t(array(micro$plant[1,], dim = c(ncol(micro$plant), length(Twb4))))
    dimnames(micro$plant) <- plantnames 
    
    micro$shadmet      <- micro$metout
    micro$shadsoil     <- micro$soil
    micro$humid        <- micro$soil
    micro$shadhumid    <- micro$soil
    micro$soilpot      <- micro$soil
    micro$shadpot      <- micro$soil
    micro$soilmoist    <- micro$soil
    micro$shadmoist    <- micro$soil
    micro$shadtcond    <- micro$tcond
    micro$shaddensit   <- micro$densit
    micro$shadspecheat <- micro$specheat
    micro$shadplant    <- micro$plant
    micro$RAINFALL     <- rep(0, length(Twb4)/24)
    micro$minshade     <- rep(0, length(Twb4)/24)
    micro$maxshade     <- rep(90, length(Twb4)/24)
    micro$nyears       <- ceiling(length(Twb4)/24/365)
    
    # replace first line of microclimate input with air temperature and humidity at which to compute Tw
    micro$metout[1:length(Twb4),c(3,4,14)] <- c(Ta2[start:end], rep(Ta2[1], extra)) # set air and sky temperature
    micro$metout[1:length(Twb4),c(5,6)]    <- c(RH2[start:end], rep(RH2[1], extra)) # set relative humidity
    # micro$metout[1:length(Twb3),c(7,8)] <- 0.05 # set wind speed to low (sensitive to this)
    micro$soil[1:length(Twb4),3] <- c(Ta2[start:end], rep(Ta2[1], extra)) # set surface soil temperature to air temperature
    micro$metout[1:length(Twb4), 2] <- rep(seq(0, 23 * 60, 60), length(Ta2)/24 + 1)[1:length(Twb4)]
    #micro$metout[1:length(Twb3), 1] <- seq(0, 23 * 60, 60)[1:length(Twb3)]
    preshr <- c(P2[start:end], rep(P2[1], extra))*1000          
    
    micro$metout[1:length(Twb4),c(7,8)] <- c(ws2[start:end], rep(ws2[1], extra)) # set wind speed    
    
    # run the ectotherm model 
    message(paste('running ectotherm simulation for month ',month,' for ', length(Twb4),' sites \n'))
    ptm <- proc.time() # Start timing
    gc()
    ecto <- ectotherm(Ww_g = 8.7,
                      live = 0,
                      pct_wet = 87,
                      shape = 4, # frog
                      preshr = preshr)
    message(paste0('runtime ', (proc.time() - ptm)[3], ' seconds')) # Stop the clock
    gc()
    environ <- as.data.frame(ecto$environ)
    masbal <- as.data.frame(ecto$masbal)
    if(k == 1){
      Tb.ectotherm <- environ$TC[1:(nrow(environ) - extra)]
      EWL.ectotherm <- masbal$H2OCut_g[1:(nrow(masbal) - extra)]
    }else{
      Tb.ectotherm <- c(Tb.ectotherm, environ$TC[1:(nrow(environ) - extra)])
      EWL.ectotherm <- c(EWL.ectotherm, masbal$H2OCut_g[1:(nrow(masbal) - extra)])
    }
  }
  Tb2 <- Ta
  Tb2[!is.na(Ta)] <- Tb.ectotherm
  Tb.grid <- Ta.w
  Tb.grid <- setValues(Tb.grid, Tb2)
  #plot(Tb.grid)
  #map('world', add = TRUE)
  
  EWL2 <- Ta
  EWL2[!is.na(Ta)] <- EWL.ectotherm
  EWL.grid <- Ta.w
  EWL.grid <- setValues(EWL.grid, EWL2)
  #plot(EWL.grid)
  EWL.grid[Tb.grid <= 0] <- NA
  #map('world', add = TRUE)
  
  writeRaster(EWL.grid, file = paste0('EWL_',month,'_2C.nc'), overwrite=TRUE)
}  

## +4C ## ---------------------------------------
# Saturated vapour pressure
tmax_full_4C <- raster::brick(file.path(spatial_path, 'TerraClimate4C_tmax.nc')) # 2m, deg C
es_kPa_4C    <- tmax_full_4C
values(es_kPa_4C) <- 0.611 * exp(2500000 / 461.5 * ( 1 / 273 - 1 / (values(tmax_full_4C) + 273.15)))

# Actual vapour pressure
ea_kPa      <- raster::brick(file.path(spatial_path, 'TerraClimate19812010_vap.nc'))
ea_kPa_4C   <- ea_kPa
values(ea_kPa_4C) <- values(ea_kPa_4C) * 0.8

# Calculate VPD (kPa)
vpd_full_4C <- ea_kPa_4C
values(vpd_full_4C) <- values(es_kPa_4C) - values(ea_kPa_4C)

Tbs  <- crop(tmax_full_4C, extent(elev))
EWLs <- Tbs 

# Run ectotherm model for each month and save on local drive (+4C scenario)
for(month in 1:12){
  Ta.w   <- crop(tmax_full_4C[[month]], extent(elev))
  esat   <- VAPPRS(Ta.w)
  e      <- esat - crop(vpd_full_4C[[month]], extent(elev)) * 1000
  RH.w   <- e / esat * 100
  RH.w[RH.w < 0] <- 0
  Wind.w <- crop(ws_full[[month]], extent(elev))
  Vstar  <- 0.4 * Wind.w / log((100 / 0.15) + 1) # friction velocity
  Wind.w <- 2.5 * Vstar * log((1 / 0.15) + 1) # correct for reference height (10 m) to ground level (1 cm) assuming level ground
  # compute atmospheric pressure (kPa)
  P.w <- 101325 * ((1 - (0.0065 * elev / 288)) ^ (1 / 0.190284)) / 1000
  
  Ta  <- getValues(Ta.w)
  RH  <- getValues(RH.w)
  P   <- getValues(P.w)
  Pa  <- getValues(e)
  ws  <- getValues(Wind.w)
  P[!is.na(Ta) & is.na(P)] <- 101.325
  Ta2 <- Ta[!is.na(Ta)]
  RH2 <- RH[!is.na(Ta)]
  P2  <- P[!is.na(Ta)]
  ws2 <- ws[!is.na(Ta)]
  Pa2 <- Pa[!is.na(Ta)]
  
  div <- 10
  Twb <- Ta2[1:floor(length(Ta2) / div)]
  micro_base <- micro_global()
  
  for(k in 1:div){
    micro <- micro_base
    if(k == 1){
      extra <- round(24 - (length(Twb) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb, rep(Ta2[1], extra))
      start <- 1
      end   <- length(Twb)
    }
    if(k > 1 & k < div){
      fact  <- (length(Twb) * (k - 1) + k - 1)
      Twb4  <- Ta2[fact:(fact + length(Twb))]
      extra <- round(24 - (length(Twb4) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb4, rep(Ta2[1], extra))
      start <- fact
      end   <- fact + length(Twb)
    }        
    if(k == div){
      fact  <- (length(Twb) * (k - 1) + k - 1)
      Twb4  <- Ta2[fact:length(Ta2)]
      extra <- round(24 - (length(Twb4) / 24)%%1 * 24, 0)
      Twb4  <- c(Twb4, rep(Ta2[1], extra))
      start <- fact
      end   <- length(Ta2)
    } 
    
    soilnames     <- dimnames(micro$soil)
    metoutnames   <- dimnames(micro$metout)
    tcondnames    <- dimnames(micro$tcond)
    specheatnames <- dimnames(micro$specheat)
    densitnames   <- dimnames(micro$densit)
    plantnames    <- dimnames(micro$plant)
    
    micro$metout   <- t(array(micro$metout[1,], dim = c(ncol(micro$metout), length(Twb4))))
    dimnames(micro$metout) <- metoutnames
    micro$soil     <- t(array(micro$soil[1,], dim = c(ncol(micro$soil), length(Twb4))))
    dimnames(micro$soil) <- soilnames
    micro$tcond    <- t(array(micro$tcond[1,], dim = c(ncol(micro$tcond), length(Twb4))))
    dimnames(micro$tcond) <- tcondnames
    micro$specheat <- t(array(micro$specheat[1,], dim = c(ncol(micro$specheat), length(Twb4))))
    dimnames(micro$specheat) <- specheatnames  
    micro$densit   <- t(array(micro$densit[1,], dim = c(ncol(micro$densit), length(Twb4))))
    dimnames(micro$densit) <- densitnames      
    micro$plant    <- t(array(micro$plant[1,], dim = c(ncol(micro$plant), length(Twb4))))
    dimnames(micro$plant) <- plantnames 
    
    micro$shadmet