Open habitats increase vulnerability of amphibian tadpoles to climate warming across latitude

Supplementary Information

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

This supplementary file contains the R workflow for processing and analysing the raw data, and creating figures for the manuscript ID: GEB-2022-0268 titled “Open habitats increase vulnerability of amphibian tadpoles to climate warming across latitude”.

---

1 Pond thermodynamic model

Water temperature is computed based on a heat budget considering longwave radiation, solar heating, and standard bulk aerodynamic flux parameterizations for sensible and latent heat (Mitchell et al., 2012). The budget is driven by fundamental meteorological conditions provided in the microclimate model of NicheMapR (Kearney and Porter, 2017), specific to each site, with the Equation (1.1).

\[\begin{equation} (Lp_w\ c_p)\ \ \frac{dT_w}{dt}=p_w\ c_p\ (I\ T_i-O\ T_w+R\ T_{sky}A_L\ )+(H-E-S+\phi_{LW} +\phi_{SW})\ A_L \tag{1.1} \end{equation}\]

where \(T_w\) is the lake temperature, \(T_i\) is the inflow temperature, \(T_{sky}\) is the air temperature (all in °K), \(p_w\) is the density of water (kg m^-3), \(c_p\) is the specific heat capacity of water (J kg^-1 K^-1), \(L\) is the lake volume (m^-3) and \(I\), \(O\), and (\(A_L\ R\)) are the inflow, outflow and rainfall volume fluxes respectively (m^-3 d^-1), where \(A_L\) is the lake surface area (m^-2) used to multiply the rainfall depth (m d^-1). \(H\) (1.2), \(E\) (1.3), \(S\) (1.4) are heat fluxes (W m^-2) for sensible heat, evaporative and soil heat conduction, respectively, calculated according to:

\[\begin{equation} H=p_a\ c_p\ C_H\ u_2 (T_{sky}-T_w) \tag{1.2} \end{equation}\]

\[\begin{equation} E=p_a\ \lambda\ C_E\ u_2\ (e_s-e_a) \tag{1.3} \end{equation}\]

\[\begin{equation} S=\frac{C_S}{z_s}(T_s-T_w) \tag{1.4} \end{equation}\]

where \(C_H\) and \(C_E\) are bulk transfer coefficients for sensible and latent heat respectively (-), \(\lambda\) is the latent heat of vaporization (J kg^-1), \(p_a\) is the density of air (kg m^-3), and \(C_S\) is the diffusivity of heat into the soil below the water (W m^-1 K^-1), and \(z_s\) is the active soil depth over which heat diffusion occurs. \(\phi_{LW}\) (1.5) is the new longwave radiation (W m^-2) and \(\phi_{SW}\) (1.6) is the solar insolation (W m^-2), calculated from:

\[\begin{equation} \phi_{LW}=\varepsilon_a\ \sigma(T_{sky})^4-\varepsilon_w\ \sigma(T_w)^4 \tag{1.5} \end{equation}\]

\[\begin{equation} \phi_{SW}=(1-\alpha)\ {\hat{\phi}}_{SW}\bigg[\frac{A_L}{A_{MAX}}\ k_s \bigg] \tag{1.6} \end{equation}\]

where \(\varepsilon\) is the emissivity of air or water, \(\sigma\) is the Stefan-Boltzmann constant, \(\alpha\) is the shortwave radiation albedo, \({\hat{\phi}}_{SW}\) is the incident solar radiation above vegetation (W m^-2). The last term (1.6) accounts for the effect of wetland vegetation shading on the incident solar intensity, with \(k_s\) defined as a constant to increase the shading effect, and \(A_{MAX}\) is the maximum area of inundation extent (m^-2). In NicheMapR, assumed values for \(C_S\), \(C_H\), \(C_E\), \(\varepsilon_a\), \(k_s\) were 0.5, 0.0013, 0.0013, 0.99 and 0.25, respectively.

---

2 Dataset

2.1 Prepare raw data for analysis

The raw data for analysis is available on GitHub. This section is the import and clean the raw data for analysis.

# Clean raw data
clean_dat <- read.csv("https://raw.githubusercontent.com/nicholaswunz/tadpole-thermal-tolerance/main/data/raw_data.csv", na.strings = c("" , "NA" ), stringsAsFactors = TRUE) %>%
  dplyr::mutate(WT           = CTmax - Tmax,
                region       = factor(region, levels = c('Subtropics', 'Temperate')),
                country      = as.factor(country),
                habitat      = as.factor(collection_habitat),
                species_tree = gsub(" ", "_", species))

# Summarise data by group
summary_dat <- as.data.frame(clean_dat %>%
  dplyr::group_by(region, country, habitat, family, species, sp_abbr, species_tree, breeding_habitat, citation) %>%
  dplyr::summarise(lon        = mean(lon),
                   lat        = mean(lat),
                   CTmax_mean = mean(CTmax),
                   CTmax_SD   = sd(CTmax),
                   Tmax       = mean(Tmax),
                   WT_mean    = mean(WT),
                   WT_SD      = sd(WT),
                   n          = n(),
                   month_start = mean(month_start),
                   month_end   = mean(month_end))) %>%
  dplyr::relocate(c(breeding_habitat, citation), .after = last_col())

2.2 Species data summary

Table S1a

Table S1. List of species included in the study with their sample size (n), mean ± s.d. critical thermal maxima (CT_max), maximum environmental water temperature (T_max), habitat type based on collection site, and spatial coordinates (latitude and longitude as decimal degree).

Species (abbreviation)	n	CTmax ± s.d. (°C)	Tmax (°C)	Collection habitat	Latitude	Longitude
Bufo bankorensis (BB)	15	37.54±0.40	18.79	Forest	23.919	120.883
Microhyla heymonsi (MH)	15	37.44±0.58	23.90	Forest	23.919	120.883
Odorrana swinhoana (OS)	30	35.85±0.51	23.55	Forest	23.988	121.129
Kurixalus berylliniris (KBE)	16	33.54±0.73	24.74	Forest	22.803	121.051
Kurixalus eiffingeri (KE)	15	31.05±1.60	25.01	Forest	23.690	120.791
Polypedates braueri (PB)	12	40.17±0.32	26.77	Forest	23.889	120.891
Zhangixalus moltrechti (ZM)	15	37.36±0.52	24.17	Forest	23.889	120.891
Duttaphrynus melanostictus (DM)	15	41.41±0.22	35.12	Open	24.181	120.609
Fejervarya limnocharis (FL)	20	41.50±0.93	30.12	Open	24.863	121.550
Hyla chinensis (HC)	12	38.23±0.39	29.14	Open	24.594	120.999
Microhyla fissipes (MF)	19	39.71±0.77	29.34	Open	24.890	121.567
Micryletta steinegeri (MS)	12	37.35±0.45	33.03	Open	23.836	120.748
Hylarana latouchii (HL)	12	37.52±0.37	21.34	Open	23.891	120.884
Rana longicrus (RL)	10	35.79±0.51	19.75	Open	25.191	121.590
Rana sauteri (RS)	14	32.44±1.14	24.89	Open	23.902	120.891
Buergeria choui (BC)	20	41.40±1.03	30.12	Open	24.863	121.550
Buergeria robusta (BR)	13	38.48±1.21	25.28	Open	23.879	120.887
Kurixalus idiootocus (KI)	10	37.40±1.35	30.00	Open	23.914	120.896
Rana pirica (RP)	8	33.91±1.48	14.50	Forest	42.595	141.281
Bufo japonicus formosus (BJF)	7	40.10±0.48	32.50	Open	43.143	141.204
Glandirana rugosa (GR)	10	39.15±0.35	32.50	Open	43.002	141.181
Bombina orientalis (BO)	30	39.93±0.67	23.29	Forest	36.813	128.054
Bufo gargarizans (BG)	48	38.04±0.61	17.95	Forest	36.491	127.197
Rana coreana (RC)	42	36.02±1.63	16.99	Forest	36.718	126.618
Rana uenoi (RU)	12	38.23±0.68	21.95	Forest	36.812	128.052
Dryophytes japonicus (DJ)	36	40.91±1.52	32.50	Open	37.599	126.892
Kaloula borealis (KBO)	6	42.65±0.26	36.84	Open	37.521	126.880
Pelophylax nigromaculatus (PN)	8	39.50±3.61	32.50	Open	36.887	126.936
Glandirana emeljanovi (GE)	30	38.10±1.07	29.05	Open	37.832	127.558

• Bufo gargarizans was recently renamed to Bufo schalinensis (Othman et al., 2022).

Table S1b

Table S1. Continuation of species summary with additional information on Family, and notes on the species typical breeding habitat based on the associated references. Species grouped as forest primarily breed in forest habitats. Species grouped as forest edge typically breed at the edges of forests (more exposure to sunlight). Species grouped as open primarily breed in open habitats. Species grouped as generalist typically breed in both forest and open habitats.

Family	Species	Breeding habitat	Citation
Bufonidae	Bufo bankorensis	Generalist	Huang et al. 1996
Microhylidae	Microhyla heymonsi	Generalist	Chuang pers. comm.
Ranidae	Odorrana swinhoana	Forest	Lai et al. 2007
Rhacophoridae	Kurixalus berylliniris	Forest	Wu et al. 2016
Rhacophoridae	Kurixalus eiffingeri	Forest	Lin et al. 2008
Rhacophoridae	Polypedates braueri	Forest	Hsu et al. 2012
Rhacophoridae	Zhangixalus moltrechti	Forest edge	Chang et al. 2014
Bufonidae	Duttaphrynus melanostictus	Open	Huang et al. 1997
Dicroglossidae	Fejervarya limnocharis	Open	Kuan et al. 2011
Hylidae	Hyla chinensis	Generalist	Chuang pers. comm.
Microhylidae	Microhyla fissipes	Open	Chuang pers. comm.
Microhylidae	Micryletta steinegeri	Open	Chuang pers. comm.
Ranidae	Hylarana latouchii	Generalist	Huang et al. 2004
Ranidae	Rana longicrus	Open	Kam et al. 1995
Ranidae	Rana sauteri	Generalist	Jang-Liaw et al. 2009
Rhacophoridae	Buergeria choui	Generalist	Tominaga et al. 2015
Rhacophoridae	Buergeria robusta	Generalist	Lin et al. 2012
Rhacophoridae	Kurixalus idiootocus	Generalist	Hou et al. 2008
Ranidae	Rana pirica	Forest	Matsui et al. 2018; Haramura pers. comm.
Bufonidae	Bufo japonicus formosus	Open	Matsui et al. 2018; Haramura pers. comm.
Ranidae	Glandirana rugosa	Open	Matsui et al. 2018; Haramura pers. comm.
Bombinatoridae	Bombina orientalis	Forest edge	Yang et al. 2001
Bufonidae	Bufo gargarizans	Forest edge	Yang et al. 2001; Ambu et al. 2022
Ranidae	Rana coreana	Generalist	Yang et al. 2001; Ambu et al. 2022
Ranidae	Rana uenoi	Generalist	Yang et al. 2001; Ambu et al. 2022
Hylidae	Dryophytes japonicus	Generalist	Yang et al. 2001; Roh et al. 2014
Microhylidae	Kaloula borealis	Open	Yang et al. 2001; Jung et al. 2013
Microhylidae	Pelophylax nigromaculatus	Open	Yang et al. 2001; Yoo et al. 2019
Ranidae	Glandirana emeljanovi	Generalist	Yang et al. 2001

• Bufo gargarizans was recently renamed to Bufo schalinensis (Othman et al., 2022).
  
• Reference list: Ambu et al. (2022), Chang et al. (2014), Hou et al. (2008), Hsu et al. (2012), Huang et al. (1996), Huang et al. (1997), Huang et al. (2004), Jang-Liaw and Lee (2009), Jung et al. (2013), Kam et al. (1995), Kuan and Lin (2011), Lai et al. (2007), Lin et al. (2012), Lin and Kam (2008), Matsui and Maeda (2018), Roh et al. (2014), Tominaga et al. (2015), Wu et al. (2016), Yang et al. (2001), Yoo et al. (2019).

---

3 Phylogeny

3.1 Species ascension gene number

Table S2. Accession number of genes of all species (ingroups and outgroups) in the analysis from the National Center for Biotechnology Information (NCBI).

Species	12S	16S	COI	cytb
Ingroup
Bufo gargarizans	DQ275350	DQ275350	DQ275350	DQ275350
Bufo japonicus formosus	LC061223	LC061223		AB713498
Bombina orientalis	AY585338	AY585338	AY585338	AY585338
Buergeria choui	DQ283055	DQ283055	MH034165	KC151119
Buergeria robusta	AB530075	AF026370	GU244379	JF802879
Bufo bankorensis	AF160768	AF160786	HQ650558	AB159260
Dryophytes japonicus	AB303949	AB303949	AB303949	AB303949
Duttaphrynus melanostictus	AY458592	AY458592	AY458592	AY458592
Fejervarya limnocharis	AY158705	AY158705	AY158705	AY158705
Glandirana emeljanovi	NC_030211	NC_030211	NC_030211	NC_030211
Glandirana rugosa	KF771341	KF771341	JQ844517	KF771341
Hyla chinensis	AY458593	AY458593	AY458593	AY458593
Hylarana latouchii	KF771284	AB058880	JN700815	EU034751
Kaloula borealis	NC_020044	NC_020044	NC_020044	NC_020044
Kurixalus berylliniris		DQ468669	DQ468677
Kurixalus eiffingeri	AB933305	AF026363	DQ468680
Kurixalus idiootocus	AB933306	DQ468674	DQ468682	GQ204503
Microhyla fissipes	AB201175	AB201185	KR087802	AB201219
Microhyla heymonsi	AY458596	AY458596	AY458596	AY458596
Micryletta steinegeri	AB634638	AB634696
Odorrana swinhoana	AB200929	KF185045	HQ650557
Pelophylax nigromaculatus	KT878718	KT878718	KT878718	KT878718
Polypedates braueri	AB728016	AB728003	KR087860
Rana coreana	NC_024548	NC_024548	NC_024548	NC_024548
Rana longicrus		AB058881		JF939130
Rana pirica	KX269184	KX269184	KX024946	KX269331
Rana sauteri	AB685767	AB211495		EU034956
Rana uenoi	KX269177	KX024885	KX024945	KX024969
Zhangixalus moltrechti	AF458145	DQ468676	HQ650556	U00710
Outgroup
Ascaphus truei	X86225	DQ283116	AJ871087	AF277330
Batrachuperus pinchonii	AY916007	DQ283340	KM201393	AY593142
Gyrinophilus porphyriticus	EU336432	DQ283255	NC_006341	NC_006341

3.2 Prepare phylogeny for analysis

This section is the import the phylogenetic tree (see phylogenetic reconstruction in main text) comprising of 29 species, match names with the data set, and create a correlation matrix for subsequent analysis.

# Load phylogeny data
phylo_tree <- ape::read.tree("https://raw.githubusercontent.com/nicholaswunz/tadpole-thermal-tolerance/main/data/tree2021_12.nwk")

tree_tip_label <- phylo_tree$tip.label # extract tree tip names
species_list   <- unique(clean_dat$species_tree) # extract species name from main dataset
missing_sp     <- as.data.frame(setdiff(tree_tip_label, species_list)) # check what names from the dataset not found in phylo_tree
# Three outgroups in tree

pruned_tree    <- ape::drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_list)) # prune phylo_tree to keep species from the main dataset

# check if lengths match for both data and tree
sp_n <- length(unique(pruned_tree$tip.label)) # 29
#length(unique(clean_dat$species_tree)) # 29

# Check for polytomies
#ape::is.binary(pruned_tree) # TRUE

# Compute branch lengths
pruned_tree_cal <- ape::compute.brlen(pruned_tree, method = "Grafen", power = 1)
#ape::is.ultrametric(pruned_tree_cal) # TRUE

# Create correlation matrix for analysis
phylo_cor <- vcv(pruned_tree_cal, cor = T)

Fig. S1. Phylogenetic reconstruction of all 29 amphibian species used in this study. Species that were collected from closed-forest habitats shown with filled triangles, while species that were collected from open habitats shown with open circles.

---

4 Microclimate models and future scenario

4.1 NicheMapR parameters

Table S3. Parameter values used to estimate pond temperature using NicheMapR. Parameters for which settings differ from default values are shown in bold. Default values for other parameters in the model remained.

ID	Value	Variable name	Variable unit
latitude	Provided in table S1	latitude coordinate	\(^{\circ}\)
longitude	Provided in table S1	latitude coordinate	\(^{\circ}\)
nyears	2	number of years to run	#
soiltype	clay (modified*)	soil type	-
REFL	0.15	soil solar reflectance	decimal %
slope	0	slope	\(^{\circ}\)
aspect	0	aspect	\(^{\circ}\)
DEP	0, 2, 5, 10, 15, 20, 30, 45, 90, 200	soil depth	cm
minshade	0 (unshaded) or 90 (shaded)	minimum shade level	%
maxshade	10 (unshaded) 100 (shaded)	maximum shade level	%
runshade	TRUE	run model with shade?	-
clearsky	FALSE	run with clear sky?	-
SLE	0.95	substrate longwave IR emissivity	decimal %
Thcond	2.5	soil mineral thermal conductivity	W m\(\\^{-1}\) K\(\\^{-1}\)
Density	2.56	soil mineral density	kg m\(\\^{-3}\)
SpecHeat	870	soil mineral heat capacity	J kg\(\\^{-1}\) K\(\\^{-1}\)
BulkDensity	1.4	soil bulk density	kg m\(\\^{-3}\)
hori	0	horizon angles	\(^{\circ}\)
cap	FALSE	organic surface layer?	-
PE	-2.731491	air entry potential	J kg\(\\^{-1}\)
KS	6.41E-05	saturated conductivity	kg s m\(\\^{-3}\)
BB	16.1914	Campbell\('\)s soil \('\)b\('\) parameter	-
BD	1.4	soil bulk density	mg m\(\\^{-3}\)
DD	2.56	soil density	mg m\(\\^{-3}\)
maxpool	500	maximum pool depth	mm
rainmult	10	rainfall multiplier	-
LAI	0	leaf area index	-
snowmodel	FALSE	run snow model?	-
soilgrids	FALSE	query soilgrids?	-
windfac	TRUE	wind multiplier	-
warm	0 (current) or 4 (warming)	warming offset vector	\(^{\circ}\)C

* textural soil properties for clay: depth (2.5, 7.5, 22.5, 45, 80, 150), bulk density 1.4 mg m\(\\^3\), 90% clay, 0% silt, 10% sand.

4.2 Define soil parameters

Define soil parameters to enable construction of hypothetical water body in NicheMapR. Parameters modified presented in Table S3.

# Define textural soil properties for clay
soilpro <- matrix(data = 0, nrow = 6, ncol = 5)
colnames(soilpro) <- c('depth', 'blkdens', 'clay', 'silt', 'sand')
soilpro <- as.data.frame(soilpro) %>%
  dplyr::mutate(depth   = c(2.5, 7.5, 22.5, 45, 80, 150),
                blkdens = 1.4, # bulk density (Mg/m3)
                clay    = 90, # % clay
                silt    = 0.0, # % silt
                sand    = 10) # % sand
DEP     <- c(0, 2, 5, 10, 15, 20, 30, 40, 50, 100) # Soil nodes (cm)

# Get hydraulic properties for given soil
soil.hydro  <- NicheMapR::pedotransfer(soilpro = soilpro, DEP = DEP)
PE <- soil.hydro$PE
BB <- soil.hydro$BB
BD <- soil.hydro$BD
KS <- soil.hydro$KS
BulkDensity <- BD[seq(1, 19, 2)]

rainmult  <- 10 # Rainfall multiplier to reflect catchment (1 will result in unadjusted rainfall)
maxpool   <- 500 # maximium pooling depth (mm), i.e. wetland depth
LAI       <- 0 # leaf area index (zero so no transpiration)
cap       <- 0 # do not simulate the organic cap layer
soilgrids <- 0 # don't use SoilGrids soil database for thermal properties

4.3 Simulate pond temperature

Run NichemapR to simulate pond temperatures under four hypothetical scenarios: (a) current climate with open habitat, (b) current climate with shaded habitat, (c) warming climate with open habitat, and (d) warming climate with shaded habitat.

# Simulate dates for 2 years
dates <- NicheMapR::micro_global(nyear = 2)$dates

# Create loop for each coordinates
est_Tmax_open <- NULL
for( i in 1:nrow(summary_dat)){ 
  result <- data.frame(
    NicheMapR::micro_global(loc      = c(summary_dat$lon[i], summary_dat$lat[i]), 
                            runshade = 0, minshade = 0, maxshade = 10, # shade parameters
                            PCTWET   = 100, # model surface soil wetness (%), 
                            PE = PE, BB = BB, BD = BD, KS = KS, rainmult = rainmult, 
                            maxpool = maxpool, LAI = LAI, cap = cap, soilgrids = soilgrids, # soil parameters to simulate pool
                            warm = 0, # current climate
                            nyear = 2,)$soil) %>% # extract soil parameters
    dplyr::mutate(dates = dates) %>%
    dplyr::filter(dates >= summary_dat$month_start[i] & dates <= summary_dat$month_end[i]) %>% 
    dplyr::summarise(est_Tmax_open = max(D0cm))
  print(result)
  results_2 <- cbind(summary_dat$species[i], result)
  est_Tmax_open <- rbind(est_Tmax_open, results_2)
}

est_Tmax_shade <- NULL
for( i in 1:nrow(summary_dat)){ 
  result <- data.frame(
    NicheMapR::micro_global(loc      = c(summary_dat$lon[i], summary_dat$lat[i]), 
                            runshade = 0, minshade = 95, maxshade = 100, # shade parameters
                            PCTWET   = 100, # model surface soil wetness (%), 
                            PE = PE, BB = BB, BD = BD, KS = KS, rainmult = rainmult, maxpool = maxpool, LAI = LAI, cap = cap, soilgrids = soilgrids, # soil parameters to simulate pool
                            warm = 0, # current climate
                            nyear = 2)$soil) %>% # Extract soil parameters
    dplyr::mutate(dates = dates) %>%
    dplyr::filter(dates >= summary_dat$month_start[i] & dates <= summary_dat$month_end[i]) %>% 
    dplyr::summarise(est_Tmax_shade = max(D0cm))
  print(result)
  est_Tmax_shade <- rbind(est_Tmax_shade, result)
}

est_Tmax_open_4C <- NULL
for( i in 1:nrow(summary_dat)){ 
  result <- data.frame(
    NicheMapR::micro_global(loc      = c(summary_dat$lon[i], summary_dat$lat[i]), 
                            runshade = 0, minshade = 0, maxshade = 10, # shade parameters
                            PCTWET   = 100, # model surface soil wetness (%), 
                            PE = PE, BB = BB, BD = BD, KS = KS, rainmult = rainmult, maxpool = maxpool, LAI = LAI, cap = cap, soilgrids = soilgrids, # soil parameters to simulate pool
                            warm = 4, # future climate
                            nyear = 2,)$soil) %>% # extract soil parameters
    dplyr::mutate(dates = dates) %>%
    dplyr::filter(dates >= summary_dat$month_start[i] & dates <= summary_dat$month_end[i]) %>% 
    dplyr::summarise(est_Tmax_open_4C = max(D0cm))
  print(result)
  est_Tmax_open_4C <- rbind(est_Tmax_open_4C, result)
}

est_Tmax_shade_4C <- NULL
for( i in 1:nrow(summary_dat)){ 
  result <- data.frame(
    NicheMapR::micro_global(loc      = c(summary_dat$lon[i], summary_dat$lat[i]), 
                            runshade = 0, minshade = 95, maxshade = 100, # shade parameters
                            PCTWET   = 100, # model surface soil wetness (%), 
                            PE = PE, BB = BB, BD = BD, KS = KS, rainmult = rainmult, maxpool = maxpool, LAI = LAI, cap = cap, soilgrids = soilgrids, # soil parameters to simulate pool
                            warm = 4, # future climate
                            nyear = 2)$soil) %>% # extract soil parameters
    dplyr::mutate(dates = dates) %>%
    dplyr::filter(dates >= summary_dat$month_start[i] & dates <= summary_dat$month_end[i]) %>% 
    dplyr::summarise(est_Tmax_shade_4C = max(D0cm))
  print(result)
  est_Tmax_shade_4C <- rbind(est_Tmax_shade_4C, result)
}

proj_data <- cbind(summary_dat, est_Tmax_shade, est_Tmax_open[,2], est_Tmax_shade_4C, est_Tmax_open_4C) %>% 
  dplyr::rename(est_Tmax_open = "est_Tmax_open[, 2]")

4.4 Verify predicted values

Predicted T_max from the heat budget models were validated with a linear model (below) and visually in Fig. S2e.

predict_obs_model <- lm(Tmax ~ est_Tmax_shade + habitat, data = proj_data)

verify    <- summary(predict_obs_model)$r.squared
model_mae <- Metrics::mae(proj_data$Tmax, predict(predict_obs_model))

Fig. S2 Maximum water temperature (T_max as degree Celsius) output from the microclimate model. Example pond temperature output over a one year period for (a) Bombina orientalis, a temperate species in Korea (36.813 °N, 128.054 °E), (b) Pelophylax nigromaculatus, a temperate species in Korea (36.887 °N, 126.936 °E), (c) Kurixalus eiffingeri, a subtropics species in Taiwan (23.914 °N, 120.896 °E), and (d) Duttaphrynus melanostictus, a subtropics species in Taiwan (24.181 °N, 120.608 °E). The solid and dashed line represents the estimated water temperature under shaded conditions (95–100% vegetation shade) and under open conditions (1–10% vegetation shade), respectively. The filled area was the environmental data collection period, and the single dashed horizontal grey line is the observed T_max from the study. (e) The observed and estimated T_max for each species, matching the study location and time. Observed temperatures were recorded by HOBO Water Temperature loggers while the estimated T_max was estimated from NicheMapR under open (filled triangles) and shaded habitats (open circles). Dashed line represents 1:1 scale, and R² = 0.7137031, with a mean absolute error of 2.3863238.

---

5 Relationship between CT_max, T_max, and warming tolerance

Run analysis looking at the relationship between (a) CT_max and T_max (CT_max–T_max model), (b) warming tolerance (WT) and T_max (WT–T_max model), and (c) WT and CT_max (WT–CT_max model).

set.seed(10)

# Set priors
prior <- c(
  brms::prior(normal(0, 5), "b"), # mean of 0 and SD of 10 (wide distribution)
  brms::prior(student_t(3, 0, 10), "sd"), # class of random effect deviation to calculate - has to be positive (no negative SD)
  brms::prior(student_t(3, 0, 10), "sigma")) # residual SD parameter

# Run CTmax_Tmax_model
CTmax_Tmax_model <- brms::brm(CTmax ~ Tmax + habitat + lat + (1 | species) + (1 | gr(species_tree, cov = phylo)),
                              data    = clean_dat,
                              family  = gaussian,
                              data2   = list(phylo = phylo_cor),
                              prior   = prior,
                              iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                              control = list(adapt_delta = 0.99, max_treedepth = 18))

# Run WT_Tmax_model
WT_Tmax_model <- brms::brm(WT ~ Tmax + habitat + lat + (1 | species) + (1 | gr(species_tree, cov = phylo)),
                           data    = clean_dat,
                           family  = gaussian,
                           data2   = list(phylo = phylo_cor),
                           prior   = prior,
                           iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                           control = list(adapt_delta = 0.99, max_treedepth = 18))

# Run WT_CTmax_model (grouped by species because of co-correlation of Ctmax and WT with Tmax)
WT_CTmax_model <- brms::brm(WT_mean ~ CTmax_mean + habitat + lat + (1 | gr(species_tree, cov = phylo)),
                            data    = summary_dat,
                            family  = gaussian,
                            prior   = prior,
                            data2   = list(phylo = phylo_cor),
                            iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                            control = list(adapt_delta = 0.99, max_treedepth = 18))

Fig. S3. Scatterplots of the observed data (y) vs the average simulated data (y_rep) from the posterior predictive distribution for the (a) CT_max–T_max model, (b) WT–T_max model, and (c) WT–CT_max model. Dashed line represents a slope of 1.

---

5.1 Model output

Table S4

Table S4. Point estimates and 95% credible intervals (as determined using Bayesian methods) for fitted parameters estimated for the correlation between CT_max and T_max, which includes fixed-effect parameters for \(\beta_0\), the average across-species intercept; \(\beta_1\), the average across-species slope of T_max; habitat (closed and open habitat); and latitude. Group-level effects include the standard deviations (sd) for species-level variation in the intercept (\(\Delta_{\beta_0[\zeta]}\)), and phylogenetic variation in the intercept (\(\Delta_{\beta_0[\phi]}\)). Phylogenetic signal was calculated as \(\lambda = \sigma_\phi^2 / (\sigma_\phi^2 + \sigma^2)\), where \(\sigma_\phi^2\) is the estimated phylogenetic variance (i.e. sd\((\Delta_{\beta_0[\phi]})^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	29.388	3.241	23.026	35.762
Tmax, \(\beta_1\)	0.241	0.117	0.009	0.471
Habitat - Open	0.268	1.220	-2.146	2.644
Latitude	0.081	0.065	-0.048	0.210
Group-level effects
sd\((\Delta_{\beta_0[\zeta]})\)	2.086	0.504	1.011	3.064
sd\((\Delta_{\beta_0[\phi]})\)	1.812	1.409	0.066	5.277
Phylogenetic signal
\(\lambda\)	0.580	0.311	0.004	0.961

Table S5

Table S5. Point estimates and 95% credible intervals (as determined using Bayesian methods) for fitted parameters estimated for the correlation between warming tolerance and T_max, which includes fixed-effect parameters for \(\beta_0\), the average across-species intercept; \(\beta_1\), the average across-species slope of T_max; habitat (closed and open habitat); and latitude. Group-level effects include the standard deviations (sd) for species-level variation in the intercept (\(\Delta_{\beta_0[\zeta]}\)), and phylogenetic variation in the intercept (\(\Delta_{\beta_0[\phi]}\)). Phylogenetic signal was calculated as \(\lambda = \sigma_\phi^2 / (\sigma_\phi^2 + \sigma^2)\), where \(\sigma_\phi^2\) is the estimated phylogenetic variance (i.e. sd\((\Delta_{\beta_0[\phi]})^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	29.445	3.334	22.853	36.016
Tmax, \(\beta_1\)	-0.761	0.117	-0.992	-0.531
Habitat - Open	0.303	1.236	-2.140	2.752
Latitude	0.081	0.065	-0.050	0.209
Group-level effects
sd\((\Delta_{\beta_0[\zeta]})\)	2.088	0.500	1.029	3.071
sd\((\Delta_{\beta_0[\phi]})\)	1.863	1.429	0.081	5.334
Phylogenetic signal
\(\lambda\)	0.589	0.310	0.006	0.963

Table S6

Table S6. Point estimates and 95% credible intervals (as determined using Bayesian methods) for fitted parameters estimated for the correlation between warming tolerance and CT_max at the species level, which includes fixed-effect parameters for \(\beta_0\), the average across-species intercept; \(\beta_1\), the average across-species slope of CT_max; habitat (closed and open habitat); and latitude. Group-level effects include the standard deviations (sd) for phylogenetic variation in the intercept (\(\Delta_{\beta_0[\phi]}\)). Phylogenetic signal was calculated as \(\lambda = \sigma_\phi^2 / (\sigma_\phi^2 + \sigma^2)\), where \(\sigma_\phi^2\) is the estimated phylogenetic variance (i.e. sd\((\Delta_{\beta_0[\phi]})^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	4.957	12.022	-18.649	28.594
CTmax, \(\beta_1\)	0.278	0.336	-0.384	0.933
Habitat - Open	-5.589	1.611	-8.700	-2.343
Latitude	-0.028	0.120	-0.262	0.213
Group-level effect
sd\((\Delta_{\beta_0[\phi]})\)	3.240	2.162	0.175	8.526
Phylogenetic signal
\(\lambda\)	0.416	0.289	0.002	0.947

---

6 Habitat and region analysis

Run analysis looking at the relationship between (a) CT_max and (b) warming tolerance (WT) between habitat (closed-forest and open) and region (subtropic and temperate).

CTmax_model <- brms::brm(CTmax ~ habitat + region + (1 | species) + (1 | gr(species_tree, cov = phylo)),
                         data    = clean_dat,
                         family  = gaussian,
                         data2   = list(phylo = phylo_cor),
                         prior   = prior,
                         iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                         control = list(adapt_delta = 0.99, max_treedepth = 18))

WT_model <- brms::brm(WT ~ habitat * region + (1 | species) + (1 | gr(species_tree, cov = phylo)),
                      data    = clean_dat,
                      family  = gaussian,
                      data2   = list(phylo = phylo_cor),
                      prior   = prior,
                      iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                      control = list(adapt_delta = 0.99, max_treedepth = 18))

WT_lat_model <- brms::brm(WT ~ lat * habitat + (1 | species) + (1 | gr(species_tree, cov = phylo)),
                      data    = clean_dat,
                      family  = gaussian,
                      data2   = list(phylo = phylo_cor),
                      prior   = prior,
                      iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                      control = list(adapt_delta = 0.99, max_treedepth = 18))

Fig. S4. Scatterplots of the observed data (y) vs the average simulated data (y_rep) from the posterior predictive distribution for the (a) CT_max model, (b) WT model, and (c) WT–latitude model. Dashed line represents a slope of 1.

---

Model output

Table S7

Table S7. Point estimates and 95% credible intervals (as determined using Bayesian methods) for fitted parameters estimated for the relationship between CT_max with habitat and region, which includes fixed-effect parameters for \(\beta_0\), the average across-species intercept; habitat (closed and open habitat); and region (subtropics, temperate). Group-level effects include the standard deviations (sd) for species-level variation in the intercept (\(\Delta_{\beta_0[\zeta]}\)), and phylogenetic variation in the intercept (\(\Delta_{\beta_0[\phi]}\)). Phylogenetic signal was calculated as \(\lambda = \sigma_\phi^2 / (\sigma_\phi^2 + \sigma^2)\), where \(\sigma_\phi^2\) is the estimated phylogenetic variance (i.e. sd\((\Delta_{\beta_0[\phi]})^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	36.548	1.378	33.888	39.434
Habitat - Open	2.133	0.928	0.321	3.963
Region - Temperate	1.688	0.987	-0.272	3.611
Group-level effects
sd\((\Delta_{\beta_0[\zeta]})\)	2.105	0.519	1.075	3.141
sd\((\Delta_{\beta_0[\phi]})\)	2.264	1.394	0.203	5.527
Phylogenetic signal
\(\lambda\)	0.697	0.258	0.036	0.965

Table S8

Table S8. Point estimates and 95% credible intervals (as determined using Bayesian methods) for fitted parameters estimated for the correlation between warming tolerance and T_max, which includes fixed-effect parameters for \(\beta_0\), the average across-species intercept; \(\beta_1\), the average across-species slope of T_max; habitat (closed and open habitat); and latitude. Group-level effects include the standard deviations (sd) for species-level variation in the intercept (\(\Delta_{\beta_0[\zeta]}\)), and phylogenetic variation in the intercept (\(\Delta_{\beta_0[\phi]}\)). Phylogenetic signal was calculated as \(\lambda = \sigma_\phi^2 / (\sigma_\phi^2 + \sigma^2)\), where \(\sigma_\phi^2\) is the estimated phylogenetic variance (i.e. sd\((\Delta_{\beta_0[\phi]})^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	13.087	1.773	9.685	16.783
Habitat - Open	-3.005	1.524	-6.064	-0.011
Region - Temperate	3.756	1.920	-0.142	7.410
Habitat - Open: Region - Temperate	-6.241	2.349	-10.699	-1.529
Group-level effects
sd\((\Delta_{\beta_0[\zeta]})\)	3.184	0.603	2.072	4.453
sd\((\Delta_{\beta_0[\phi]})\)	2.102	1.644	0.083	6.297
Phylogenetic signal
\(\lambda\)	0.625	0.307	0.006	0.973

7 Warming and deforesation analysis

Run analysis looking at the relationship between climate (current and warming), with the interaction between category (open and shaded simulated habitat), habitat (closed and open habitat - breeding) and region (subtropics, temperate) on warming tolerance.

# Calculate WT
proj_data <- proj_data %>%
  dplyr::mutate(WT_shade_cur = CTmax_mean - est_Tmax_shade,
                WT_open_cur  = CTmax_mean - est_Tmax_open,
                WT_shade_4C  = CTmax_mean - est_Tmax_shade_4C,
                WT_open_4C   = CTmax_mean - est_Tmax_open_4C)

# Convert wide to long format
proj_long <- proj_data %>%
  dplyr::select(region, species, CTmax_mean, habitat, lat, WT_shade_cur, WT_open_cur, WT_shade_4C, WT_open_4C) %>%
  tidyr::pivot_longer(!c(region, species, CTmax_mean, habitat, lat), names_to = "category", values_to = "values") %>%
  dplyr::mutate(climate = case_when(endsWith(category, "cur") ~ "Current",
                                    endsWith(category, "4C")~ "Warming"),
                category = dplyr::recode(category, 
                                         WT_open_cur = "Open habitat",
                                         WT_shade_cur = "Shaded habitat",
                                         WT_open_4C = "Open habitat",
                                         WT_shade_4C = "Shaded habitat"))

# Add phylogeny to dataset
proj_long$species_tree <- summary_dat$species_tree[match(proj_long$species, summary_dat$species)]

# Run brms model
proj_model <- brms::brm(values ~ climate + category * habitat * region + (1 | species) + (1 | gr(species_tree, cov = phylo)),
                        data    = proj_long,
                        family  = gaussian,
                        data2   = list(phylo = phylo_cor),
                        iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                        control = list(adapt_delta = 0.99, max_treedepth = 18))

Fig. S5. Scatterplot of the observed data (y) vs the average simulated data (yrep) from the posterior predictive distribution for the climate and deforestation projection model. Dashed line represents a slope of 1.

7.1 Model output

Table S9. Point estimates and 95% credible intervals (as determined using Bayesian methods) for fitted parameters estimated for the relationship between warming tolerance with climate change and deforestation, which includes fixed-effect parameters for \(\beta_0\), the average across-species intercept; climate (current, warming); and the three-way interaction with shaded (open and shaded simulated habitat), habitat (closed and open habitat - breeding) and region (subtropics, temperate). Group-level effects include the standard deviations (sd) for species-level variation in the intercept (\(\Delta_{\beta_0[\zeta]}\)), and phylogenetic variation in the intercept (\(\Delta_{\beta_0[\phi]}\)). Phylogenetic signal was calculated as \(\lambda = \sigma_\phi^2 / (\sigma_\phi^2 + \sigma^2)\), where \(\sigma_\phi^2\) is the estimated phylogenetic variance (i.e. sd\((\Delta_{\beta_0[\phi]})^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	9.097	2.239	4.688	13.798
Climate - Warming	-3.266	0.145	-3.554	-2.983
Shaded	6.845	0.292	6.266	7.412
Habitat - Open	-0.072	1.863	-3.681	3.610
Region - Temperate	-3.538	2.552	-8.637	1.472
Shaded: Habitat - Open	0.068	0.375	-0.656	0.803
Shaded: Region - Temperate	4.977	0.452	4.098	5.856
Habitat - Open: Region - Temperate	2.713	3.232	-3.559	9.076
Shaded: Habitat - Open: Region - Temperate	-2.206	0.597	-3.376	-1.036
Group-level effects
sd\((\Delta_{\beta_0[\zeta]})\)	3.333	0.855	1.286	4.901
sd\((\Delta_{\beta_0[\phi]})\)	2.793	2.390	0.091	8.739
Phylogenetic signal
\(\lambda\)	0.735	0.302	0.014	0.992

---

8 Figures for main text

Raw figures produced were modified in Adobe Illustrator for publication.

8.1 Figure 2 - CT_max, T_max, warming tolerance relationship

# Extract marginal effects
CTmax_Tmax_me <- as.data.frame(ggeffects::ggpredict(CTmax_Tmax_model, terms = c("Tmax[sample = 25]"))) 
WT_Tmax_me    <- as.data.frame(ggeffects::ggpredict(WT_Tmax_model, terms = c("Tmax[sample = 25]"))) 

# Fig 2
ctmax_tmax_plot <- summary_dat %>%
  ggplot() +
  geom_ribbon(data = CTmax_Tmax_me, aes(x = x, ymin = conf.low, ymax = conf.high), fill = "#FE945C",alpha = 0.1) +
  geom_line(data = CTmax_Tmax_me, aes(x = x, y = conf.low), colour = "#FE945C", linetype = "dashed") +
  geom_line(data = CTmax_Tmax_me, aes(x = x, y = conf.high), colour = "#FE945C", linetype = "dashed") +
  geom_line(data = CTmax_Tmax_me, aes(x = x, y = predicted), colour = "#DD3B24", size = 1) +
  geom_point(aes(x = Tmax, y = CTmax_mean, colour = region, shape = habitat), size = 2, fill = "white") +
  scale_shape_manual(values = c(17, 21)) +
  scale_colour_manual(values = c("#DD3B24", "#FE945C")) +
  xlab(expression(T["max"]~"(°C)")) +
  ylab(expression("CT"["max"]~"(°C)")) +
  mytheme()

wt_tmax_plot <- summary_dat %>%
  ggplot() +
  geom_ribbon(data = WT_Tmax_me, aes(x = x, ymin = conf.low, ymax = conf.high), fill = "#FE945C",alpha = 0.1) +
  geom_line(data = WT_Tmax_me, aes(x = x, y = conf.low), colour = "#FE945C", linetype = "dashed") +
  geom_line(data = WT_Tmax_me, aes(x = x, y = conf.high), colour = "#FE945C", linetype = "dashed") +
  geom_line(data = WT_Tmax_me, aes(x = x, y = predicted), colour = "#DD3B24", size = 1) +
  geom_point(aes(x = Tmax, y = WT_mean, colour = region, shape = habitat), size = 2, fill = "white") +
  scale_shape_manual(values = c(17, 21)) +
  scale_colour_manual(values = c("#DD3B24", "#FE945C")) +
  xlab(expression(T['max']~"(°C)")) +
  ylab("Warming tolerance") +
  mytheme()

wt_ctmax_plot <- summary_dat %>%
  ggplot() +
  geom_point(aes(x = CTmax_mean, y = WT_mean, colour = region, shape = habitat), size = 2, fill = "white") +
  scale_shape_manual(values = c(17, 21)) +
  scale_colour_manual(values = c("#DD3B24", "#FE945C")) +
  xlab(expression(CT['max']~"(°C)")) +
  ylab("Warming tolerance") +
  mytheme()

cowplot::plot_grid(ctmax_tmax_plot + theme(legend.position = "none"), 
                   wt_tmax_plot + theme(legend.position = "none"), 
                   wt_ctmax_plot + theme(legend.position = "none"),
                   ncol = 3, align = "hv", axis = "tblr", labels = c("a", "b", "c"))

8.2 Figure 3 - CT_max and warming tolerance across habitat and region

Latitude as category (Fig. 3a and 3b)

# Fig 3 production
CTmax_plot <- summary_dat %>%
  # regrouping
  dplyr::group_by(region) %>%
  dplyr::arrange(habitat, CTmax_mean) %>%
  mutate(x = row_number()) %>% # change to number to order by habitat and CTmax
  ungroup() %>%
  # plotting
  ggplot(aes(x = x, y = CTmax_mean, colour = region, shape = habitat, label = sp_abbr)) +
  stat_mean_line(aes(linetype = habitat),  alpha = 0.5) +
  geom_pointrange(aes(ymin = CTmax_mean - CTmax_SD, ymax = CTmax_mean + CTmax_SD, linetype = habitat)) +
  geom_point(size = 3, fill = "white") +
  scale_shape_manual(values = c(17, 21)) +
  scale_colour_manual(values = c("#DD3B24", "#FE945C")) +
  geom_text(nudge_y = 1.5, size = 3, show.legend = FALSE) +
  xlab(NULL) +
  ylab(expression("CT"["max"]~"(°C)")) +
  facet_grid(. ~ region, scales = 'free_x') +
  mytheme() +
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_blank(),
        axis.ticks.x = element_blank())

WT_plot <- summary_dat %>%
  # regrouping
  dplyr::group_by(region) %>%
  dplyr::arrange(habitat, WT_mean) %>%
  mutate(x = row_number()) %>% # change to number to order by habitat and WT
  ungroup() %>%
  # plotting
  ggplot(aes(x = x, y = WT_mean, colour = region, shape = habitat, label = sp_abbr)) +
  stat_mean_line(aes(linetype = habitat), alpha = 0.5) +
  geom_pointrange(aes(ymin = WT_mean - WT_SD, ymax = WT_mean + WT_SD, linetype = habitat)) +
  geom_point(size = 3, fill = "white") +
  scale_shape_manual(values = c(17, 21)) +
  scale_colour_manual(values = c("#DD3B24", "#FE945C")) +
  geom_text(nudge_y = 2, size = 3, show.legend = FALSE) +
  ylim(0, 25) +
  xlab("Species") +
  ylab("Warming tolerance") +
  facet_grid(. ~ region, scales = 'free_x') +
  mytheme() +
  theme(axis.title.x = element_blank(),
        axis.text.x  = element_blank(),
        axis.ticks.x = element_blank())

cowplot::plot_grid(CTmax_plot, WT_plot, nrow = 2, align = "hv", axis = "tblr", labels = c("a", "b"))

Latitude as continuous (Fig. 3c)

WT_lat_me <- as.data.frame(ggeffects::ggpredict(WT_lat_model, terms = c("lat [sample = 26]", "habitat")))

summary_dat %>%
  ggplot() +
  geom_ribbon(data = WT_lat_me, aes(x = x, ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.1) +
  geom_line(data = WT_lat_me, aes(x = x, y = conf.low, group = group), colour = "#DD3B24", linetype = "dashed") +
  geom_line(data = WT_lat_me, aes(x = x, y = conf.high, group = group), colour = "#DD3B24", linetype = "dashed") +
  geom_line(data = WT_lat_me, aes(x = x, y = predicted, linetype = group), colour = "#DD3B24", size = 1) +
  geom_point(aes(x = lat, y = WT_mean, shape = habitat), colour = "#DD3B24", size = 2, fill = "white") +
  scale_shape_manual(values = c(17, 21)) +
  scale_fill_manual(values = c("#FE945C", "white")) +
  xlab("Latitude (°)") +
  ylab("Warming tolerance") +
  mytheme()

8.3 Figure 4 - Warming and deforestation projection

# Extract estimate +/- 95% CI
conditions <- make_conditions(proj_model, c("climate", "region"))
predict_me <- brms::conditional_effects(proj_model, "category:habitat", conditions = conditions)
predict_dat <- as.data.frame(predict_me[[1]]) %>% dplyr::rename(estimate = estimate__, ci.lb = lower__, ci.ub = upper__)

proj_long %>%
  dplyr::mutate(category = forcats::fct_relevel(as.factor(category),"Shaded habitat", "Open habitat")) %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_point(aes(x = category, y = values, colour = region, shape = habitat), size = 2, position = position_dodge(0.6), fill = "grey", alpha = 0.5) +
  geom_pointrange(data = predict_dat, aes(x = category, y = estimate, ymin = ci.lb, ymax = ci.ub, colour = region, shape = habitat, linetype = habitat), 
                  position = position_dodge(0.6), size = 0.7, fill = "white") +
  ylab("Warming tolerance") + xlab(NULL) +
  scale_shape_manual(values = c(17, 21)) +
  scale_colour_manual(values = c("#DD3B24", "#FE945C")) +
  scale_fill_manual(values = c("#DD3B24", "#FE945C")) +
  ylim(-10, max(proj_long$values)) +
  facet_grid(. ~ climate) +
  mytheme() +
  theme(legend.position = "bottom")

---

9 [20/01/2023 CORRECTION] Microclimate models and future scenario

The micro_global function did not have the soil moisture function on by default. After consultation with the NicheMapR developer Michael Kearney, we re-estimated the pond temperatures using the micro_era5 function which provides fine-scale resolution (~1 km weather model forecast.

The model does not confirm that pond drying will ‘not occur’ as stated in the main text. However, the higher rain multiplier (rainmult) reduces the likelihood of pond drying.

The equations from the Pond thermodynamic model do not apply to the NicheMapR models. We followed the methods from Enriquez-Urzelai et al. (2019) to estimate pond temperature.

9.1 Download microclimate spatial map

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)
xmn <- min(summary_dat$lon) - 1 
xmx <- max(summary_dat$lon) + 1 
ymn <- min(summary_dat$lat) - 1
ymx <- max(summary_dat$lat) + 1

# temporal extent
st_time <- lubridate::ymd("2013:01:01") # earliest sampling date
en_time <- lubridate::ymd("2020: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)

9.2 Simulate pond temperature

Run NichemapR to simulate pond temperatures under four hypothetical scenarios: (a) current climate with open habitat, (b) current climate with shaded habitat, (c) warming climate with open habitat, and (d) warming climate with shaded habitat. To ensure the pond does not dry up during the sampling period, we increased rainmult to 20.

sum_dat <- as.data.frame(clean_dat %>%
  dplyr::group_by(region, country, habitat, family, species) %>%
  dplyr::summarise(lon        = mean(lon),
                   lat        = mean(lat),
                   Tmax       = mean(Tmax),
                   year_start = min(year_start),
                   year_end   = min(year_finish),
                   month_start = mean(month_start),
                   month_end   = mean(month_end)))

# Create loop for each coordinates
est_Tmax <- NULL
for( i in 1:nrow(sum_dat)){ 
  micro_out <- 
    NicheMapR::micro_era5(loc = c(sum_dat$lon[i], sum_dat$lat[i]), 
                          runshade = 1, minshade = 0, # shade parameters
                          PE = PE, BB = BB, BD = BD, KS = KS, rainmult = 20, 
                          maxpool = maxpool, LAI = LAI, cap = cap, soilgrids = soilgrids, # soil parameters to simulate pool
                          warm = 0, # current climate
                          DEP = DEP,
                          dstart = paste("01/01/", sum_dat$year_start[i], sep =""), 
                          dfinish = paste("31/12/", sum_dat$year_end[i], sep =""),
                          runmoist = 1, 
                          spatial = '/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/NicheMapR/era5', # change to your location
                          save = 0)
  
  metout <- data.frame(micro_out$metout) %>% # extract soil parameters
  dplyr::mutate(dates = row_number() / 730) %>% # hours to month
  dplyr::filter(dates >= sum_dat$month_start[i] & dates <= sum_dat$month_end[i]) 
  
  # extract water temp in open habitat
    open <- data.frame(micro_out$soil) %>% # extract soil parameters
      dplyr::mutate(dates = row_number() / 730) %>% # hours to month
      merge(metout, by = "dates") %>%
      dplyr::filter(dates >= sum_dat$month_start[i] & dates <= sum_dat$month_end[i] & POOLDEP > 400) %>% 
      dplyr::summarise(est_Tmax_open = max(D50cm))
  print(open)

  # extract water temp in shaded habitat
  closed <- data.frame(micro_out$shadsoil) %>% # extract soil parameters
      dplyr::mutate(dates = row_number() / 730) %>% # hours to month
      merge(metout, by = "dates") %>%
    dplyr::filter(dates >= sum_dat$month_start[i] & dates <= sum_dat$month_end[i] & POOLDEP > 400) %>% 
    dplyr::summarise(est_Tmax_shade = max(D50cm))
  print(closed)
  results <- cbind(sum_dat$species[i], open, closed)
  est_Tmax <- rbind(est_Tmax, results)
}

est_Tmax_4C <- NULL
for( i in 1:nrow(sum_dat)){ 
  micro_out <- 
    NicheMapR::micro_era5(loc = c(sum_dat$lon[i], sum_dat$lat[i]), 
                          runshade = 1, minshade = 0, # shade parameters
                          PE = PE, BB = BB, BD = BD, KS = KS, rainmult = 20, 
                          maxpool = maxpool, LAI = LAI, cap = cap, soilgrids = soilgrids, # soil parameters to simulate pool
                          warm = 4, # future climate
                          DEP = DEP,
                          dstart = paste("01/01/", sum_dat$year_start[i], sep =""), 
                          dfinish = paste("31/12/", sum_dat$year_end[i], sep =""),
                          runmoist = 1, 
                          spatial = '/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/NicheMapR/era5', # change to your location
                          save = 0)
  metout <- data.frame(micro_out$metout) %>% # extract soil parameters
  dplyr::mutate(dates = row_number() / 730) %>% # hours to month
  dplyr::filter(dates >= sum_dat$month_start[i] & dates <= sum_dat$month_end[i]) 
  
  # extract water temp in open habitat
    open <- data.frame(micro_out$soil) %>% # extract soil parameters
      dplyr::mutate(dates = row_number() / 730) %>% # hours to month
      merge(metout, by = "dates") %>%
      dplyr::filter(dates >= sum_dat$month_start[i] & dates <= sum_dat$month_end[i] & POOLDEP > 400) %>% 
      dplyr::summarise(est_Tmax_open_4C = max(D50cm))
  print(open)
  # extract water temp in shaded habitat
  closed <- data.frame(micro_out$shadsoil) %>% # extract soil parameters
      dplyr::mutate(dates = row_number() / 730) %>% # hours to month
      merge(metout, by = "dates") %>%
    dplyr::filter(dates >= sum_dat$month_start[i] & dates <= sum_dat$month_end[i] & POOLDEP > 400) %>% 
    dplyr::summarise(est_Tmax_shade_4C = max(D50cm))
  print(closed)
  results <- cbind(sum_dat$species[i], open, closed)
  est_Tmax_4C <- rbind(est_Tmax_4C, results)
}

proj_data <- cbind(summary_dat, est_Tmax[,2], est_Tmax[,3], est_Tmax_4C[,2], est_Tmax_4C[,3]) %>% 
  dplyr::rename(est_Tmax_open = "est_Tmax[, 2]", 
                est_Tmax_shade = "est_Tmax[, 3]",
                est_Tmax_open_4C = "est_Tmax_4C[, 2]", 
                est_Tmax_shade_4C = "est_Tmax_4C[, 3]")

9.3 Verify predicted values

Predicted T_max from the heat budget models were validated with a linear model (below) and visually in Updated Fig. S2.

proj_test <- proj_data %>%
  mutate(est_Tmax_test = ifelse(habitat == "Forest", est_Tmax_shade, est_Tmax_open)) %>%
  filter_all(all_vars(!is.infinite(.)))

predict_obs_model <- lm(Tmax ~ est_Tmax_test, data = proj_test)

verify2    <- summary(predict_obs_model)$r.squared
model_mae2 <- Metrics::mae(proj_test$Tmax, predict(predict_obs_model))

## downloading DEM via package elevatr 
## extracting weather data locally from /Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/NicheMapR/era5 
## computing radiation and elevation effects with package microclima 
## Downscaling radiation and wind speed 
## Calculating meso-scale terrain effects 
## running microclimate model for  365  days from  2014-01-01  to  2014-12-31 23:00:00  at site  long 121.1295 lat 23.98788 
## Note: the output column `SOLR` in metout and shadmet is for unshaded solar radiation adjusted for slope, aspect and horizon angle 
##    user  system elapsed 
##   5.642   0.075   5.777

Updated Fig. S2 The observed and estimated T_max for each species, matching the study location and time. Observed temperatures were recorded by HOBO Water Temperature loggers while the estimated T_max was estimated from NicheMapR under open (filled triangles) and shaded habitats (open circles). Dashed line represents 1:1 scale, and R² = 0.5621393, with a mean absolute error of 2.9461639.

The updated estimated pond temperatures did not differ much relative to the original simulations.

---

10 References

Ambu, J., Andersen, D. and Borzée, A. (2022). Spawning site selection and segregation at the landscape, habitat and microhabitat scales for five syntopic asian amphibians. Amphibia-Reptilia 43, 169–180.

Chang, Y., Tseng, W., Chen, C., Huang, C., Chen, Y. and Hatch, K. (2014). Winter breeding and high tadpole densities may benefit the growth and development of tadpoles in a subtropical lowland treefrog. Journal of Zoology 294, 154–160.

Hou, W.-S., Chang, Y.-H. and Wang, H.-W. (2008). Climatic effects and impacts of lakeshore bank designs on the activity of Chirixalus idiootocus in Yilan, Taiwan. Ecological Engineering 32, 52–59.

Hsu, J.-L., Kam, Y.-C. and Fellers, G. M. (2012). Overwintering tadpoles and loss of fitness correlates in Polypedates braueri tadpoles that use artificial pools in a lowland agroecosystem. Herpetologica 68, 184–194.

Huang, W.-S., Lin, J.-Y. and Yu, J. Y.-L. (1996). The male reproductive cycle of the toad, Bufo bankorensis, in Taiwan. Zoological Studies 35, 128–137.

Huang, W.-S., Lin, J.-Y. and Yu, J. Y.-L. (1997). Male reproductive cycle of the toad Bufo melanostictus in Taiwan. Zoological Science 14, 497–503.

Huang, W.-S., Cheng, Y.-S. and Tu, H.-Y. (2004). Reproductive patterns of two sympatric ranid frogs, Rana latouchii and R. sauteri, with comments on anuran breeding seasons in Taiwan. Collection and Research 17, 1–10.

Jang-Liaw, N.-H. and Lee, T.-H. (2009). Intraspecific relationships of populations of the brown frog Rana sauteri (Ranidae) on Taiwan, inferred from mitochondrial cytochrome b sequences. Zoological Science 26, 608–616.

Jung, Y. S., Park, M. O. and Koo, B. H. (2013). A substitute habitat planning for Kaloula borealis based on wetlands. Journal of the Korean Society of Environmental Restoration Technology 16, 1–15.

Kam, Y.-C., Wang, C.-S. and Lin, Y.-S. (1995). Reproduction and diet of the brown frog Rana longicrus in Taiwan. Zoological Studies 34, 193–201.

Kearney, M. R. and Porter, W. P. (2017). NicheMapR – an r package for biophysical modelling: The microclimate model. Ecography 40, 664–674.

Kuan, S. and Lin, Y. (2011). Bigger or faster? Spring and summer tadpole cohorts use different life‐history strategies. Journal of Zoology 285, 165–171.

Lai, Y., Shieh, B. and Kam, Y. (2007). Population patterns of a riparian frog (Rana swinhoana) before and after an earthquake in subtropical Taiwan. Biotropica 39, 731–736.

Lin, Y.-S. and Kam, Y.-C. (2008). Nest choice and breeding phenology of an arboreal-breeding frog, Kurixalus eiffingeri (Rhacophoridae), in a bamboo forest. Zoological Studies 47, 129.

Lin, H.-D., Chen, Y.-R. and Lin, S.-M. (2012). Strict consistency between genetic and topographic landscapes of the brown tree frog (Buergeria robusta) in Taiwan. Molecular Phylogenetics and Evolution 62, 251–262.

Matsui, M. and Maeda, N. (2018). Encyclopedia of japanese frogs Japan book. Japan: Bun-ichiSogo Shuppan.

Mitchell, N., Hipsey, M. R., Arnall, S., McGrath, G., Tareque, H. B., Kuchling, G., Vogwill, R., Sivapalan, M., Porter, W. P. and Kearney, M. R. (2012). Linking eco-energetics and eco-hydrology to select sites for the assisted colonization of australia’s rarest reptile. Biology 2, 1–25.

Othman, S. N., Litvinchuk, S. N., Maslova, I., Dahn, H., Messenger, K. R., Andersen, D., Jowers, M. J., Kojima, Y., Skorinov, D. V. and Yasumiba, K. (2022). From Gondwana to the Yellow sea, evolutionary diversifications of true toads Bufo sp. In the Eastern Palearctic and a revisit of species boundaries for Asian lineages. Elife 11, e70494.

Roh, G., Borzée, A. and Jang, Y. (2014). Spatiotemporal distributions and habitat characteristics of the endangered treefrog, Hyla suweonensis, in relation to sympatric H. japonica. Ecological Informatics 24, 78–84.

Tominaga, A., Matsui, M., Eto, K. and Ota, H. (2015). Phylogeny and differentiation of wide-ranging Ryukyu Kajika frog Buergeria japonica (Amphibia: Rhacophoridae): Geographic genetic pattern not simply explained by vicariance through strait formation. Zoological Science 32, 240–247.

Wu, S.-P., Huang, C.-C., Tsai, C.-L., Lin, T.-E., Jhang, J.-J. and Wu, S.-H. (2016). Systematic revision of the Taiwanese genus Kurixalus members with a description of two new endemic species (Anura, Rhacophoridae). ZooKeys 121.

Yang, S. Y., Kim, J. B., Min, M. S., Suh, J. H. and Kang, Y. J. (2001). Monograph of Korean Amphibia. Seoul, Republic of Korea: Academy book Press.

Yoo, N., Do, M. S., Nam, H.-K., Choi, G., Son, S.-J. and Yoo, J.-C. (2019). Habitat characteristics of anuran species inhabiting rice fields of western mid-south Korea - In the case of Daeho reclamation agricultural land by farming practices. Korean Journal of Ecology and Environment 52, 366–377.

---

10.1 Session Information

R version 4.2.2 (2022-10-31)

Platform: aarch64-apple-darwin20 (64-bit)

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: microclima(v.0.1.0), RNetCDF(v.2.6-1), mcera5(v.0.3.0), raster(v.3.6-13), sp(v.1.5-1), pander(v.0.6.5), NicheMapR(v.3.2.1), diversitree(v.0.9-16), phytools(v.1.2-0), maps(v.3.4.1), ape(v.5.6-2), rstan(v.2.21.7), StanHeaders(v.2.21.0-7), brms(v.2.18.0), Rcpp(v.1.0.9), colorspace(v.2.0-3), forcats(v.0.5.2), stringr(v.1.5.0), dplyr(v.1.0.10), purrr(v.1.0.1), readr(v.2.1.3), tidyr(v.1.2.1), tibble(v.3.1.8), tidyverse(v.1.3.2), Matrix(v.1.5-1), cowplot(v.1.1.1), ggplot2(v.3.4.0) and knitr(v.1.41)

loaded via a namespace (and not attached): utf8(v.1.2.2), tidyselect(v.1.2.0), htmlwidgets(v.1.6.1), hoardr(v.0.5.2), grid(v.4.2.2), combinat(v.0.0-8), munsell(v.0.5.0), units(v.0.8-1), codetools(v.0.2-18), DT(v.0.26), miniUI(v.0.1.1.1), withr(v.2.5.0), Brobdingnag(v.1.2-9), progressr(v.0.13.0), highr(v.0.10), rstudioapi(v.0.14), stats4(v.4.2.2), Metrics(v.0.1.4), bayesplot(v.1.10.0), labeling(v.0.4.2), emmeans(v.1.8.3), mnormt(v.2.1.1), optimParallel(v.1.0-2), farver(v.2.1.1), bridgesampling(v.1.1-2), coda(v.0.19-4), vctrs(v.0.5.1), generics(v.0.1.3), clusterGeneration(v.1.3.7), xfun(v.0.36), timechange(v.0.2.0), slippymath(v.0.3.1), R6(v.2.5.1), markdown(v.1.4), fields(v.14.1), cachem(v.1.0.6), assertthat(v.0.2.1), promises(v.1.2.0.1), scales(v.1.2.1), RNCEP(v.1.0.10), googlesheets4(v.1.0.1), gtable(v.0.3.1), ncmeta(v.0.3.5), processx(v.3.8.0), phangorn(v.2.10.0), spam(v.2.9-1), rlang(v.1.0.6), scatterplot3d(v.0.3-42), rgdal(v.1.6-4), gargle(v.1.2.1), broom(v.1.0.2), checkmate(v.2.1.0), inline(v.0.3.19), yaml(v.2.3.6), reshape2(v.1.4.4), abind(v.1.4-5), modelr(v.0.1.10), threejs(v.0.3.3), crosstalk(v.1.2.0), backports(v.1.4.1), httpuv(v.1.6.8), tensorA(v.0.36.2), tools(v.4.2.2), tcltk(v.4.2.2), bookdown(v.0.31), ellipsis(v.0.3.2), jquerylib(v.0.1.4), posterior(v.1.3.1), RColorBrewer(v.1.1-3), proxy(v.0.4-27), plyr(v.1.8.8), progress(v.1.2.2), base64enc(v.0.1-3), classInt(v.0.4-8), elevatr(v.0.4.2), ps(v.1.7.2), prettyunits(v.1.1.1), rpart(v.4.1.19), tgp(v.2.4-21), viridis(v.0.6.2), deSolve(v.1.34), zoo(v.1.8-11), cluster(v.2.1.4), haven(v.2.5.1), fs(v.1.5.2), crul(v.1.3), magrittr(v.2.0.3), colourpicker(v.1.2.0), reprex(v.2.0.2), googledrive(v.2.0.0), mvtnorm(v.1.1-3), matrixStats(v.0.63.0), hms(v.1.1.2), shinyjs(v.2.1.0), mime(v.0.12), evaluate(v.0.19), xtable(v.1.8-4), XML(v.3.99-0.13), shinystan(v.2.6.0), readxl(v.1.4.1), gridExtra(v.2.3), ggeffects(v.1.1.4), rstantools(v.2.2.0), compiler(v.4.2.2), tidync(v.0.3.0), ncdf4(v.1.21), KernSmooth(v.2.23-20), crayon(v.1.5.2), htmltools(v.0.5.4), later(v.1.3.0), tzdb(v.0.3.0), expm(v.0.999-6), RcppParallel(v.5.1.5), lubridate(v.1.9.0), DBI(v.1.1.3), sjlabelled(v.1.2.0), dbplyr(v.2.2.1), subplex(v.1.8), rappdirs(v.0.3.3), MASS(v.7.3-58.1), sf(v.1.0-9), cli(v.3.6.0), quadprog(v.1.5-8), parallel(v.4.2.2), insight(v.0.18.8), dotCall64(v.1.0-2), igraph(v.1.3.5), pkgconfig(v.2.0.3), numDeriv(v.2016.8-1.1), terra(v.1.6-47), xml2(v.1.3.3), dygraphs(v.1.1.1.6), bslib(v.0.4.2), estimability(v.1.4.1), rvest(v.1.0.3), distributional(v.0.3.1), callr(v.3.7.3), digest(v.0.6.31), httpcode(v.0.3.0), rmarkdown(v.2.19), cellranger(v.1.1.0), fastmatch(v.1.1-3), curl(v.5.0.0), shiny(v.1.7.4), gtools(v.3.9.4), lifecycle(v.1.0.3), nlme(v.3.1-160), jsonlite(v.1.8.4), viridisLite(v.0.4.1), fansi(v.1.0.3), pillar(v.1.8.1), lattice(v.0.20-45), loo(v.2.5.1), fastmap(v.1.1.0), httr(v.1.4.4), plotrix(v.3.8-2), pkgbuild(v.1.4.0), glue(v.1.6.2), xts(v.0.12.2), shinythemes(v.1.2.0), class(v.7.3-20), stringi(v.1.7.12), sass(v.0.4.4), rnoaa(v.1.3.8), e1071(v.1.7-12) and maptree(v.1.4-8)

---

1. Hawkesbury Institute for the Environment, Western Sydney University, NSW 2753, Australia, nicholas.wu.nz@gmail.com↩︎

2. I would like to thank J. Kong (Trinity College Dublin) for showing me the way of Rmarkdown.↩︎