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 km2) 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

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 (yrep) 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 H2O 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 (Tskin, °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