Supplementary Information

Authors

Nicholas C Wu

Tomás Villada-Cadavid

Christopher Turbill

Published

26th Jun 2025

GitHub Repository: The Quarto code to produce the HTML 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 titled “Local climate and genetic influence on intraspecific patterns in torpor physiology of a cave-roosting bat”.

Supplementary Methods

Respirometry set-up

Fig. S1 | Schematic diagram of the flow-through respirometry set up to measure oxygen consumption, carbon dioxide production, and evaporative water loss of seven torpored Miniopterus orianae oceanensis and one control blank chamber. The mass flow controllers were placed outside the temperature-controlled incubator to change the flow rate during the experiment manually.

Calculations and conversation

The rate of oxygen consumption ($\dot{V}\textrm{O}_2$; ml O₂ min^-1) was calculated following Lighton (2019) in Equation 1:

\[ \dot{V}\textrm{O}_2 = \frac{\textrm{FR}^\textrm{'}_{\textrm{i}}[(F_{\textrm{i}}\textrm{O}_{2}-F_{\textrm{e}}\textrm{O}_{2})-F_{\textrm{e}}\textrm{O}_{2}(F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}-F^\textrm{'}_{\textrm{i}}\textrm{CO}_{2})]}{(1-F_{\textrm{e}}\textrm{O}_{2})} \tag{1}\]

where $\textrm{FR}^\textrm{'}_{\textrm{i}}$ is the incurrent dry flow rate corrected for water vapour ($\textrm{FR}^\textrm{'}_{\textrm{i}} = \textrm{FR}_{\textrm{i}} \times \textrm{BP}/(\textrm{BP}-\textrm{WVP})$). $F_{\textrm{i}}\textrm{O}_{2}$ and $F_{\textrm{e}}\textrm{O}_{2}$ is the dry incurrent (baseline) and excurrent O₂ gas concentration, respectively. $F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}$ and $F^\textrm{'}_{\textrm{i}}\textrm{CO}_{2}$ is the dry” equivalent of the excurrent and incurrent CO2 gas concentration corrected for water vapour, respectively. The rate of carbon dioxide production ($\dot{V}\textrm{CO}_2$; ml CO₂ min^-1) was calculated as Equation 2:

\[ \dot{V}\textrm{CO}_2 = \frac{\textrm{FR}^\textrm{'}_{\textrm{i}}[(F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}-F^\textrm{'}_{\textrm{i}}\textrm{CO}_{2})-F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}(F_{\textrm{i}}\textrm{O}_{2}-F_{\textrm{e}}\textrm{O}_{2})]}{(1-F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2})} \tag{2}\]

Metabolic rate (MR; kJ h^-1) was calculated from $\dot{V}\textrm{O}_2$ and $\dot{V}\textrm{CO}_2$ as Equation 3:

\[ \textrm{MR} = 0.06 \times (3.941 \times \dot{V}\textrm{O}_2 + 1.106 \times \dot{V}\textrm{CO}_2) \times 4.184 \tag{3}\]

The evaporative water loss (EWL) or $\dot{V}\textrm{H}_2\textrm{O}$ (µg H2O min^-1) was calculated as Equation 4:

\[ \dot{V}\textrm{H}_2\textrm{O} = \textrm{FR}_{\textrm{i}} \left( \frac{\textrm{WVP}}{T \times R_{\textrm{w}}} \right) \tag{4}\]

where WVP is the water vapour pressure (kPa) minus the baseline WVP, $T$ is temperature in Kelvin (°C + 273.15), and $R_{\textrm{w}}$ is the gas constant for water vapour (461.5 J kg⁻¹ K⁻¹). $\dot{V}\textrm{H}_2\textrm{O}$ was converted to mg H₂O h^-1 by multiplying by 0.06. Prior to calculating $\dot{V}\textrm{O}_2$, $\dot{V}\textrm{CO}_2$, and $\dot{V}\textrm{H}_2\textrm{O}$, the baseline drift and time lag of the fractional gas concentrations were corrected for.

Cave microclimate

The local air temperature and relative humidity outside and inside the cave was recorded to monitor differences in microclimate conditions available for cave-roosting bats. An external HOBO logger (MX2302A; Onset, MA, USA) mounted with a radiation shield (RS3-B; Onset) was placed outside ~20 m away at 1.5 m height from the cave entrance to log external microclimate at 30 second intervals. An internal HOBO logger (MX2301A; Onset) was placed inside the cave above ground in the same chamber where bats were observed roosting to log cave microclimate at 30 second intervals. Cave surface temperature and relative humidity was measured with an infrared thermometer (TG56 IR 30:1; FLIR, OR, USA) and HOBO logger on an 8 m extendable pole (MX2302A) every 10 m inside the cave or where bats were found. Cave surface temperature were measured in three replicates at three time points, May, July, and September. Vapour pressure deficit (kPa) was calculated as the difference between the saturation water vapour pressure ($e_s$) and the actual water vapour pressure ($e_a$). $e_s$ was calculated following the Clausius-Clapeyron equation (Stull, 2000).

Fig. S2 | Distribution of measured (a) cave surface temperatures (°C), (b) cave relative humidity (RH; %), and (c) calculated cave vapour pressure deficit (VPD; kPa) from Jenolan and Yessabah during the 2023 Austral winter period (May–September). Surface temperature, RH, and VPD with bats present were coloured black, while surface temperature, RH, and VPD without bats present were coloured grey. Dashed lines represented the experimental exposure temperatures which are within the expected temperature range in the caves.

Process Data

Metadata of the imported respirometry data for analysis.

Name	Description
pop	Unique identifier of population source (Jenolan, Yessabah).
date_exp	The date when the experiment was conducted in day/month/year
ID	Unique identifier of the individual.
ID_gene	Unique identifier of the biopsy sample for genetic analysis
sex	Sex of the individual (female, male)
FA_mm	Forearm length in millimetre (mm).
mass_g	Body mass in grams (g).
temp_exp	Set air temperature exposure in degrees Celcius (°C).
temp_act	Recorded air temperature in degrees Celcius (°C).
RH_treatment	Recorded relative humidity in percentage (%).
chamber_n	Unique identifier of the metabolic chamber
setFR_ml_min	Set flow rate (ml/min).
BP_kPa	Recorded barometric pressure (kPa).
FR_adj_ml_min	Adjusted-flow rate (ml/min).
T_start	Start time of the respirometry measurments as year-month-day hour:min:sec.
T_end	End time of the respirometry measurments as year-month-day hour:min:sec.
t_tor_h	Time between onset of torpor and respirometry measurments.
VO2_ml_min	Calculated rate of oxygen consumption in millilitres per minute (ml/min).
VCO2_ml_h	Calculated rate of carbon dioxide production in millilitres per minute (ml/min).
EWL_mg_h	Calculated rate of evaporative water loss in milligrams per minute (mg/h).
Tsk_C	Calculated skin surface temperature in degrees Celcius (°C).
VO2_ml_h	Calculated rate of oxygen consumption in millilitres per hour (ml/h).
VCO2_ml_h	Calculated rate of carbon dioxide production in millilitres per hour (ml/h).
VO2_ml_g_h	Mass-specific oxygen consumption in ml/g/h.
VCO2_ml_g_h	Mass-specific carbon dioxide production in ml/g/h.
EWL_mg_g_h	Mass-specific evaporative water loss in mg/g/h.
lnMass	The natural logarithm of body mass.
MR_kJ_h	Calculated metabolic rate in kilojoules per hour (kJ/h).
MR_J_h	Calculated metabolic rate in joules per hour (J/h).
MR_W	Calculated metabolic rate in Watts (W).
MR_mW	Calculated metabolic rate in milliWatts (mW).
t_diff	Difference between skin temperature and actual air temperature in degrees Celcius (°C).

Load and clean data

Load and clean the raw data for analysis. Table S1 summarises data on torpor physiology traits by site, sex, and temperature exposure.

Code

# Respirometry data
resp_dat <- read.csv(file.path(data_path, "data/resp_raw_dat.csv")) %>%
  dplyr::select(-c(file_name, FR_ml_min, transmitter_n, freq_hz)) %>%
  dplyr::filter(cond == "dry" & ID != c("MAM-18", "MAM-16")) %>% # remove individuals not in torpor
  dplyr::mutate(pop         = factor(pop),
                VO2_ml_h    = VO2_ml_min * 60,
                VCO2_ml_h   = VCO2_ml_min * 60,
                VO2_ml_g_h  = VO2_ml_h / mass_g,
                VCO2_ml_g_h = VCO2_ml_h / mass_g,
                EWL_mg_g_h  = EWL_mg_h / mass_g,
                lnMass      = log(mass_g),
                T_start     = lubridate::dmy_hm(T_start),
                T_end       = lubridate::dmy_hm(T_end),
                MR_kJ_h     = (0.06 * (3.941 * VO2_ml_min + 1.106 * VCO2_ml_min)) * 4.184, # Convert VO2 and VCO2 to MR kJ h
                MR_J_h      = MR_kJ_h * 1000,
                MR_W        = MR_kJ_h * 0.2777777778, # convert to Watt
                MR_mW       = MR_W * 1000,
                lnVO2       = log(VO2_ml_h),
                lnVCO2      = log(VCO2_ml_h),
                lnMR        = log(MR_J_h),
                lnEWL       = log(EWL_mg_h),
                t_diff      = Tsk_C - temp_act)  

# filter torpor data only
torpor_dat <- resp_dat %>%
  dplyr::filter(Tsk_C < 28)

resp_sum <- torpor_dat %>%
  dplyr::group_by(pop, temp_exp) %>%
  dplyr::summarise(n         = length(VO2_ml_h),
                   EWL_mean  = mean(EWL_mg_h, na.rm = T),
                   EWL_sd    = sd(EWL_mg_h, na.rm = T),
                   Tsk_mean  = mean(Tsk_C, na.rm = T),
                   Tsk_sd    = sd(Tsk_C, na.rm = T),
                   MR_mean   = mean(MR_kJ_h, na.rm = T),
                   MR_sd     = sd(MR_kJ_h, na.rm = T)
                   )

Table S1 | Summary data from the respirometry experiment. Sample size (n), mean body mass (g), mean ± standard deviation rate of oxygen consumption ($\dot{V}\textrm{O}_2$, ml h^-1), rate of carbon dioxide production ($\dot{V}\textrm{CO}_2$, ml h^-1), torpor metabolic rate (TMR, J h^-1), torpor evaporative water loss (TEWL, mg h^-1), the skin surface temperature (T_sk, °C), and Pd detection via qPCR by sex, air temperature exposure (T_a, °C), and sites.

Site	T_a (°C)	n	Sex	Mass (g)	$\dot{V}\textrm{O}_2$ (ml h^-1)	$\dot{V}\textrm{CO}_2$ (ml h^-1)	TMR (J h^-1)	TEWL (mg h^-1)	T_sk (°C)
Jenolan	5	4	female	13.68	2.37±1.43	2.05±1.20	48.61±29.10	4.31±1.45	15.48±1.49
Jenolan	5	3	male	14.27	2.66±1.94	2.26±1.76	54.39±40.20	4.35±1.31	16.67±0.51
Jenolan	8	4	female	13.68	3.27±2.02	2.93±1.98	67.55±42.53	4.80±0.74	16.93±2.60
Jenolan	8	3	male	14.27	4.13±2.50	3.61±2.20	84.86±51.28	4.92±1.08	17.67±0.19
Jenolan	11	6	female	14.03	6.42± 6.36	6.08±6.13	133.98±133.17	5.01±1.04	18.65±2.26
Jenolan	11	5	male	14.40	13.35±10.07	11.89±9.86	275.17±211.41	11.46±5.31	20.23±3.09
Jenolan	15	6	female	13.62	19.95±18.51	16.52±14.75	405.40±373.30	13.75±6.54	21.24±2.58
Jenolan	15	5	male	14.40	24.33±14.73	21.05±12.32	498.53±299.35	21.67±9.04	23.35±2.28
Yessabah	5	7	female	13.26	0.72±0.33	0.69±0.30	15.01± 6.82	3.37±0.67	14.97±0.46
Yessabah	5	8	male	15.19	1.65±0.68	1.54±0.63	34.27±14.15	3.79±0.31	15.19±0.39
Yessabah	8	7	female	13.26	1.19±0.48	1.15±0.49	24.98±10.19	3.83±0.49	16.10±0.44
Yessabah	8	8	male	15.19	2.22±1.17	2.07±1.13	46.23±24.47	6.31±3.68	16.05±0.53
Yessabah	11	7	female	13.26	1.77±0.90	1.67±0.84	36.99±18.73	5.31±1.84	17.13±0.31
Yessabah	11	8	male	15.19	2.61±1.37	2.43±1.35	54.32±28.73	9.42±5.70	17.27±0.46
Yessabah	15	7	female	13.26	3.75±2.05	3.40±1.74	77.56±41.80	21.54±6.12	20.42±1.91
Yessabah	15	8	male	15.19	3.43±1.68	3.18±1.58	71.37±34.82	23.46±6.89	19.65±0.52

SNPs and genetic relatedness

The following code was written with Ian Brennan to process, filter and trim the raw SNP data (containing 46,224 loci, 4,075 nucleotide sites) using the dartR package (Mijangos et al., 2022). Note: the raw SNP data are being used as part of a larger study examining the genetic connectivity and population dynamics of Australian bent-winged bats. The SNP data will be available upon it’s publication (Planckaert et al., In preparation).

Code

# Load DArT file
# Visually check SNP data
dartR.base::gl.smearplot(bat_gl_wu)

nLoc(bat_gl_wu) # 46224 loci 
nInd(bat_gl_wu) # 76 individuals

## FILTERING ## ------------------------------------------
# Check raw data
dartR.base::gl.report.callrate(bat_gl_wu) # loci callrate 
dartR.base::gl.report.callrate(bat_gl_wu, method = 'ind') # individual callrate
dartR.base::gl.report.reproducibility(bat_gl_wu) # reproducibility

# strict threshold on repeatability
bat_gl2 <- dartR.base::gl.filter.reproducibility(bat_gl_wu, threshold = 0.999)

# remove invariant (monomorphic) sites, we will have to do this
# again later, once we've filtered down the dataset and potentially
# caused new monomorphs
bat_gl3 <- dartR.base::gl.filter.monomorphs(bat_gl2)

# adjust threshold on callrate as necessary
dartR.base::gl.report.callrate(bat_gl3) # loci callrate 
bat_gl4 <- dartR.base::gl.filter.callrate(bat_gl3, method = "loc", threshold = 0.9)

# filter out multiple SNPs from the same (linked) locus
bat_gl5 <- dartR.base::gl.filter.secondaries(bat_gl4, method = "best", v = 2)

# filter by the minor allele frequency, I wouldn't exceed 0.05
# in small datasets you worry about losing loci, so can be ~0.04
# in larger datasets should probably be closer to ~0.01
bat_gl6 <- dartR.base::gl.filter.maf(bat_gl5, 
                                     threshold = (2 / nInd(bat_gl5)), 
                                     v = 5)

# one more pass to clean up callrate and remove monomorphs introduced
dartR.base::gl.report.callrate(bat_gl6, method = 'ind') # individual callrate
bat_gl7 <- dartR.base::gl.filter.callrate(bat_gl6, 
                                          method = "ind", 
                                          threshold = 0.95, 
                                          mono.rm = T, 
                                          recalc = T, 
                                          recursive = T, 
                                          v = 5)

# Export a Concatenated SNP Alignment
# Export with amiguities (A/G = R; C/T=Y) for heterozygous sites
dartR.base::gl2fasta(bat_gl7, 
                     method = 4,
                     outfile = "bat_gl7_mth4.fas",
                     outpath = "local_file",
                     verbose = 5)
# Methods 2 and 4 randomly resolve heterozygous sites

# read in the alignment as a DNAbin object via ape
seq.tmp <- seqinr::read.fasta("local_file",
                      seqtype = "DNA",
                      forceDNAtolower = F)

for (i in 1:length(seq.tmp)){
  # remove any empty entries in the alignment
  seq.tmp[[i]] <- seq.tmp[[i]][-which(seq.tmp[[i]] == " ")]
}
# make a dictionary for ambiguity codes
amb.codes <- list(M = c("A","C"),R = c("A","G"), W = c("A","T"),S = c("C","G"),Y = c("C","T"),K = c("G","T"))
# loop through and replace any ambiguous bases by randomly choosing one or the other appropriate base
for (i in 1:length(seq.tmp)){
  seq.trash <- seq.tmp[[i]]
  rep.trash <- which(seq.trash %in% names(amb.codes))
  for (j in rep.trash){
    seq.trash[j] <- sample(amb.codes[[which(names(amb.codes) == seq.trash[j])]], 1)
  }
  seq.tmp[[i]] <- seq.trash
}
# By randomly resolving the heterozygous positions, some sites in the alignment
# are rendered invariable. Let's get rid of them. We'll do this manually (not using gl.filter.monomorphs)
# because our sequence data are in DNAbin format, not genlight (GL).
del.pos <- c()
for (i in 1:length(seq.tmp[[1]])){
  pos.trash <- unlist(purrr::map(seq.tmp,i))
  length.trash <- length(unique(pos.trash))
  if(length.trash == 1){
    del.pos <- c(del.pos,i)
  } else {
    if(length.trash ==2 & "N" %in% pos.trash){
      del.pos <- c(del.pos,i)
    }
  }
}
for(i in 1:length(seq.tmp)){
  seq.tmp[[i]] <- seq.tmp[[i]][-del.pos]
}

seqinr::write.fasta(seq.tmp,
                    names = names(seq.tmp),
                    nbchar = length(seq.tmp[[1]]),
                    file.out = "local_file")

The resulting fasta file was processed through the IQ-TREE software to construct the genetic relatedness tree with ascertainment bias correction (Nguyen et al., 2015).

Prepare tree for analysis

Code

phylo_tree <- ape::read.tree("http://raw.githubusercontent.com/nicholaswunz/bat-resp/refs/heads/main/files/bat_gl7_mth3.treefile")

tree_tip_label <- phylo_tree$tip.label # extract tree tip names
ID_list   <- unique(torpor_dat$ID_gene) # extract individual name from respirometry dataset
#as.data.frame(setdiff(tree_tip_label, ID_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, ID_list)) # prune phylo_tree to keep species from the main dataset

# check if lengths match for both data and tree
#length(unique(pruned_tree$tip.label)) # 27
#length(unique(torpor_dat$ID_gene)) # 27

# 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 covariance matrix for analysis
phylo_cor <- ape::vcv(pruned_tree_cal, cor = F)

Analysis

Run phylogenetic multi-level model via brm function. We constructed four chains with 5,000 steps per chain, including 2,500-step warm-up periods, hence a total of 10,000 steps were retained to estimate posterior distributions [i.e., (5,000 − 2,500) × 4 = 10,000]. An alpha delta was .99 to decrease the number of divergent transitions. Fixed effects were assigned weakly informative priors following a Gaussian distribution to speed up model convergence, and Student’s t prior with three degrees of freedom was used for group-level, hierarchical effects. The degree of convergence of model was deemed to be achieved when the Gelman–Rubin statistic (Gelman and Rubin, 1992), was 1.

Code

Tsk_mod <- brms::brm(Tsk_C ~ temp_exp * pop + sex + lnMass + tor_start + (1 | ID) +
    (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor),
    chains = 4, cores = 4, iter = 10000, warmup = 5000, control = list(adapt_delta = 0.99))

MR_mod <- brms::brm(lnMR ~ temp_exp * pop + sex + lnMass + tor_start + (1 | ID) +
    (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor),
    chains = 4, cores = 4, iter = 10000, warmup = 5000, control = list(adapt_delta = 0.99))

MR_mod2 <- brms::brm(lnMR ~ temp_exp * pop + sex + lnMass + Tsk_C + tor_start + (1 |
    ID) + (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor),
    chains = 4, cores = 4, iter = 10000, warmup = 5000, control = list(adapt_delta = 0.99))

EWL_mod <- brms::brm(lnEWL ~ temp_exp * pop + sex + lnMass + tor_start + (1 | ID) +
    (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor),
    chains = 4, cores = 4, iter = 10000, warmup = 5000, control = list(adapt_delta = 0.99))

Fig. S3 | Scatterplots of the observed data (y) vs the average simulated data (y_rep) from the posterior predictive distribution for the (a) T_sk model, (b) the TMR model, and (c) the TEWL model. Dashed line represents a slope of 1.

Model output

Table S2 | Summary statistics testing the effect of air temperature on the bat skin surface temperature (T_sk). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%), which includes the intercept ($\beta_0$), exposure temperature (T_a), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass, the time since onset of torpor, and the interaction between T_a and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is 68.32 ± 2.88%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is 87.01 ± 1.56%.

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
$\beta_0$	9.27	15.52	-21.18	40.32
T_a	0.52	0.05	0.43	0.63
Population - Yessabah	-1.70	1.02	-3.71	0.28
Sex - Male	0.20	0.91	-1.62	1.99
lnMass	1.99	5.92	-9.81	13.68
Torpor onset	-0.13	0.15	-0.44	0.17
T_a: Population - Yessabah	-0.04	0.06	-0.16	0.08
Group-level effects
$\sigma_{chamber}^2$	0.44	0.40	0.02	1.46
$\sigma_{ID}^2$	1.57	0.34	0.99	2.32
$\sigma_{phylogeny}^2$	0.61	0.55	0.02	2.06
Heritability
$h^2$	0.25	0.24	0.00	0.81

Table S3 | Summary statistics testing the effect of air temperature on torpor metabolic rate (TMR). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%), which includes the intercept ($\beta_0$), exposure temperature (T_a), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass, the time since onset of torpor, and the interaction between T_a and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is 66.31 ± 4.18%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is 86.77 ± 1.56%.

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
$\beta_0$	-8.32	6.07	-20.08	3.71
T_a	0.14	0.02	0.10	0.18
Population - Yessabah	-0.93	0.41	-1.74	-0.11
Sex - Male	-0.04	0.35	-0.74	0.62
lnMass	4.56	2.32	-0.02	9.08
Torpor onset	-0.07	0.06	-0.19	0.05
T_a: Population - Yessabah	-0.03	0.02	-0.08	0.02
Group-level effects
$\sigma_{chamber}^2$	0.23	0.21	0.01	0.80
$\sigma_{ID}^2$	0.56	0.17	0.18	0.89
$\sigma_{phylogeny}^2$	0.42	0.35	0.01	1.32
Heritability
$h^2$	0.41	0.30	0.00	0.92

Table S4 | Summary statistics testing the effect of air temperature on torpor evaporative water loss (TEWL). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%), which includes the intercept ($\beta_0$), exposure temperature (T_a), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass,the time since onset of torpor, and the interaction between T_a and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is 73.62 ± 2.62%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is 77.36 ± 2.25%.

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
$\beta_0$	-1.77	2.17	-6.03	2.49
T_a	0.13	0.02	0.09	0.16
Population - Yessabah	-0.47	0.25	-0.96	0.02
Sex - Male	0.20	0.13	-0.05	0.45
lnMass	0.96	0.83	-0.67	2.60
Torpor onset	-0.04	0.02	-0.08	0.01
T_a: Population - Yessabah	0.05	0.02	0.01	0.09
Group-level effects
$\sigma_{chamber}^2$	0.10	0.08	0.00	0.31
$\sigma_{ID}^2$	0.08	0.06	0.00	0.21
$\sigma_{phylogeny}^2$	0.15	0.09	0.01	0.35
Heritability
$h^2$	0.15	0.13	0.00	0.49

Table S5 | Summary statistics testing the effect of air temperature and skin surface temperature on torpor metabolic rate (TMR). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%). The model includes the intercept ($\beta_0$), exposure temperature (T_a), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass, T_sk, the time since onset of torpor, and the interaction between T_a and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is 66.31 ± 4.18%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is 86.77 ± 1.56%.

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
$\beta_0$	-10.10	4.44	-18.87	-1.32
T_a	0.03	0.03	-0.02	0.08
Population - Yessabah	-0.53	0.33	-1.18	0.11
Sex - Male	-0.08	0.25	-0.58	0.41
lnMass	4.05	1.70	0.70	7.43
T_sk	0.21	0.04	0.14	0.29
Torpor onset	-0.04	0.05	-0.14	0.05
T_a: Population - Yessabah	-0.02	0.02	-0.06	0.02
Group-level effects
$\sigma_{chamber}^2$	0.25	0.20	0.01	0.76
$\sigma_{ID}^2$	0.38	0.12	0.14	0.62
$\sigma_{phylogeny}^2$	0.30	0.23	0.01	0.88
Heritability
$h^2$	0.36	0.27	0.00	0.86

Figures

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

Figure 1 - Study sites

Download elevation map from rnaturalearth and convert to dataframe.

Code

# Australia map
aus_map <- rnaturalearth::ne_states(country = "Australia", returnclass = "sf")

# Elevation map
aus_elev_sp <- elevatr::get_elev_raster(locations = aus_map, z = 7, clip = "locations")  # 1km resolution or 30 arc sec
# aus_elev_sp <- readRDS(file.path(data_path, '/aus_elev_sp.rds'))

aus_elev_df <- raster::as.data.frame(aus_elev_sp, xy = TRUE)
colnames(aus_elev_df)[3] <- "elevation"
aus_elev_df <- aus_elev_df[complete.cases(aus_elev_df), ]  # Remove NAs

aus_elev_df_2 <- aus_elev_df %>%
    dplyr::filter(x >= 143 & y <= -24)

Produce eastcoast map with study sites, and temperature profiles.

Code

# Miniopterus distribution
range_map <- raster::shapefile(file.path(data_path, "data/Miniopterus_orianae_oceanensis/Miniopterus_orianae_oceanensis.shp"))
range_map <- sp::spTransform(range_map, raster::crs("+proj=longlat +datum=WGS84 +no_defs"))

# Sites
resp_sites_df <- data.frame(x = c(150.018564, 152.688690), 
                            y = c(-33.803127, -31.093943), 
                            site = c("Jenolan","Yessabah"))

map_plot <- ggplot() +
  geom_polygon(data = range_map, aes(x = long, y = lat, group = group), colour = NA, fill = "#f595a6") +
  geom_raster(data = aus_elev_df_2, aes(x = x, y = y, fill = elevation), alpha = 0.5, show.legend = F) +
  scale_fill_gradient2(low = "white", mid = "white", high = "black", midpoint = 0) +
  geom_sf(data = aus_map, color = "#8C8C8C", fill = NA) +
  geom_point(data = resp_sites_df, aes(x = x, y = y), shape = 21, colour = "black", fill = "white", size = 3) +
  xlim(143, 155) + ylim(-39.5, -24) +
  coord_sf() +
  theme_void()

# Load all files 
HOBO_list  <- list.files(path = file.path(data_path, "data/"), pattern = "\\HOBO.csv$", full.names = T)
HOBO_list2 <- data.frame(file_name = HOBO_list, id = as.character(1:length(HOBO_list))) # id for joining

HOBO_dat <- HOBO_list %>% 
  lapply(read_csv) %>% # read all files at once
  dplyr::bind_rows(.id = "id") %>% # bind all tables into one object, and give id for each
  data.frame() %>%
  dplyr::mutate(date_time = as.POSIXct(date_time, format = "%m/%d/%Y %H:%M"),
                date      = as.POSIXct(date_time, format = "%m/%d/%Y"),
                time      = strftime(date_time, format = "%H:%M:%S"),
                jul_day   = as.POSIXlt(date_time)$yday,
                month     = format(as.Date(date_time, format = "%m/%d/%Y"),"%m"),
                month_cat = format(as.Date(date_time, format = "%m/%d/%Y"),"%B"),
                location = factor(location, levels = c("EXT", "deep"))) 

temp_plot <- ggplot() +
  geom_line(data = HOBO_dat %>% dplyr::filter(location == "EXT"), 
            aes(x = jul_day, y = air_temp), size = 0.5, colour = "#f595a6") +
  geom_line(data = HOBO_dat %>% dplyr::filter(location == "deep"), 
            aes(x = jul_day, y = air_temp), size = 1.2, colour = "black") +
  ylab("Air temperature (°C)") + xlab(NULL) +
  mytheme() +
  theme(legend.position = "bottom") +
  facet_wrap(~ site, ncol = 1)

cowplot::plot_grid(map_plot, temp_plot, ncol = 2, rel_widths = c(1, 0.5))

Figure 2 - T_sk, TMR, and TEWL

Code

Tsk_me <- as.data.frame(ggeffects::ggpredict(Tsk_mod, terms = c("temp_exp", "pop", "tor_start [5]", "lnMass [2.7]")))

Tsk_plot <- ggplot() +
  # raw data
  geom_line(data = torpor_dat, aes(x = temp_exp, y = Tsk_C, group = ID), 
            colour = "lightgrey",
            alpha = 0.5,
            position = position_dodge(0.2)) +
  geom_point(data = torpor_dat, 
             aes(x = temp_exp, y = Tsk_C, shape = pop), 
             size = 1, position = position_dodge(0.2), alpha = 0.5) +
  # summary data
  geom_errorbar(data = Tsk_me, aes(x = x, ymin = conf.low, ymax = conf.high, group = group), 
                width = 0.1, position = position_dodge(0.2)) +
  geom_line(data = Tsk_me, aes(x = x, y = predicted, group = group, linetype = group), position = position_dodge(0.2)) +
  geom_point(data = Tsk_me, aes(x = x, y = predicted, shape = group, fill = group), size = 3, position = position_dodge(0.2)) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("white", "black")) +
  xlab(expression(italic("T")["a"]~"(°C)")) +
  ylab(expression(italic("T")["sk"]~"(°C)")) +
  mytheme()

MR_me <- as.data.frame(ggeffects::ggpredict(MR_mod, terms = c("temp_exp", "pop", "tor_start [5]", "lnMass [2.7]"))) %>%
  dplyr::mutate(est = exp(predicted),
                low = exp(conf.low),
                high = exp(conf.high))

MR_plot <- ggplot() +
  # raw data
  geom_line(data = torpor_dat, aes(x = temp_exp, y = MR_J_h, group = ID), 
            colour = "lightgrey",
            alpha = 0.5,
            position = position_dodge(0.2)) +
  geom_point(data = torpor_dat, 
             aes(x = temp_exp, y = MR_J_h, shape = pop), 
             size = 1, position = position_dodge(0.2), alpha = 0.5) +
  # summary data
  geom_errorbar(data = MR_me, aes(x = x, ymin = low, ymax = high, group = group), 
                width = 0.1, position = position_dodge(0.2)) +
  geom_line(data = MR_me, aes(x = x, y = est, group = group, linetype = group), position = position_dodge(0.2)) +
  geom_point(data = MR_me, aes(x = x, y = est, shape = group, fill = group), size = 3, position = position_dodge(0.2)) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("white", "black")) +
  xlab(expression(italic("T")["a"]~"(°C)")) +
  ylab(expression("TMR (J h"^"-1"*")")) +
  scale_y_continuous(breaks = c(1, 10, 25, 50, 100, 200, 400, 800), trans = "sqrt") +
  mytheme()

EWL_me <- as.data.frame(ggeffects::ggpredict(EWL_mod, terms = c("temp_exp", "pop", "tor_start [5]", "lnMass [2.7]"))) %>%
  dplyr::mutate(est = exp(predicted),
                low = exp(conf.low),
                high = exp(conf.high))

EWL_plot <- ggplot() +
  # raw data
  geom_line(data = torpor_dat, aes(x = temp_exp, y = EWL_mg_h, group = ID), 
            colour = "lightgrey",
            alpha = 0.5,
            position = position_dodge(0.2)) +
  geom_point(data = torpor_dat, 
             aes(x = temp_exp, y = EWL_mg_h, shape = pop), 
             size = 1, position = position_dodge(0.2), alpha = 0.5) +
  # summary data
  geom_errorbar(data = EWL_me, aes(x = x, ymin = low, ymax = high, group = group), 
                width = 0.1, position = position_dodge(0.2)) +
  geom_line(data = EWL_me, aes(x = x, y = est, group = group, linetype = group), position = position_dodge(0.2)) +
  geom_point(data = EWL_me, aes(x = x, y = est, shape = group, fill = group), size = 3, position = position_dodge(0.2)) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("white", "black")) +
  xlab(expression(italic("T")["a"]~"(°C)")) +
  ylab(expression("TEWL (mg H"["2"]*"O"~"h"^"-1"*")")) +
  scale_y_continuous(breaks = c(1, 2.5, 5, 10, 20, 30), trans = "sqrt") +
  mytheme()

MR_me2 <- as.data.frame(ggeffects::ggpredict(MR_mod2, terms = c("Tsk_C", "pop", "tor_start [5]", "lnMass [2.7]"))) %>%
  dplyr::mutate(est = exp(predicted),
                low = exp(conf.low),
                high = exp(conf.high))

MR_plot2 <- ggplot() +
  # raw data
  geom_line(data = torpor_dat, aes(x = Tsk_C, y = MR_J_h, group = ID), 
            colour = "lightgrey",
            alpha = 0.5,
            position = position_dodge(0.2)) +
  geom_point(data = torpor_dat, 
             aes(x = Tsk_C, y = MR_J_h, shape = pop), 
             size = 1, position = position_dodge(0.2), alpha = 0.5) +
  # summary data
  geom_ribbon(data = MR_me2, aes(x = x, ymin = low, ymax = high,  group = group), alpha = 0.1) +
  geom_line(data = MR_me2, aes(x = x, y = est, group = group, linetype = group)) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("white", "black")) +
  xlab(expression(italic("T")["sk"]~"(°C)")) +
  ylab(expression("MR (J h"^"-1"*")")) +
  mytheme()

cowplot::plot_grid(MR_plot + theme(legend.position = "none"), 
                   EWL_plot + theme(legend.position = "none"), 
                   Tsk_plot + theme(legend.position = "none"), 
                   MR_plot2 + theme(legend.position = "none"), 
                   nrow = 2, labels = c('a', 'b', 'c', 'd'))

Figure 3 - Phylogeny

Code

MR_mod_post  <- posterior_samples(MR_mod) # extract posterior
EWL_mod_post <- posterior_samples(EWL_mod) # extract posterior

# Loop to extract within ID-level posterior of intercept
ID_unique  <- sort(unique(torpor_dat$ID)) # pop ID names
out_list   <- vector(mode = 'list', length = length(ID_unique))
for (j in seq_along(ID_unique)) {
  out_list[[j]]                 <- data.frame(matrix(NA, nrow(MR_mod_post), 1))
  dat_names                     <-  paste0(ID_unique[j])
  names(out_list[[j]])          <-  dat_names
  out_list[[j]][, dat_names[1]] <-  MR_mod_post[, paste0('r_ID[', ID_unique[j], ',Intercept]')]
}

MR_Int_out  <- do.call('cbind.data.frame', out_list)
MR_Int_post <- sample_n(MR_Int_out, 4000) # Extract 4000 random rows

out_list <- vector(mode = 'list', length = length(ID_unique))
for (j in seq_along(ID_unique)) {
  out_list[[j]]                 <- data.frame(matrix(NA, nrow(EWL_mod_post), 1))
  dat_names                     <-  paste0(ID_unique[j])
  names(out_list[[j]])          <-  dat_names
  out_list[[j]][, dat_names[1]] <-  EWL_mod_post[, paste0('r_ID[', ID_unique[j], ',Intercept]')]
}

EWL_Int_out  <- do.call('cbind.data.frame', out_list)
EWL_Int_post <- sample_n(EWL_Int_out, 4000) # Extract 4000 random rows

# Tidy data
MR_Int_post_long <- MR_Int_post %>%
  tidyr::pivot_longer(cols = everything(), names_to = "ID", values_to = "Intercept")

EWL_Int_post_long <- EWL_Int_post %>%
  tidyr::pivot_longer(cols = everything(), names_to = "ID", values_to = "Intercept")

MR_Int_post_sum <- MR_Int_post_long %>%
  dplyr::group_by(ID) %>%
  dplyr::summarise(mean_MR = mean(Intercept))

EWL_Int_post_sum <- EWL_Int_post_long %>%
  dplyr::group_by(ID) %>%
  dplyr::summarise(mean_EWL = mean(Intercept))

phylo_dat <- torpor_dat %>%
  dplyr::group_by(ID_gene, ID, sex, pop) %>%
  summarise(mean_mass = mean(mass_g))

x <- tidytree::as_tibble(pruned_tree_cal)
y <- dplyr::full_join(phylo_dat, MR_Int_post_sum, by = 'ID')
y <- dplyr::full_join(y, EWL_Int_post_sum, by = 'ID')

tree_dat  <- tidytree::as.treedata(dplyr::full_join(x, y %>% rename(label = ID_gene), by = 'label'))

tree_plot <- ggtree::ggtree(tree_dat) + 
  ggtree::geom_tippoint(aes(shape = pop, fill = pop), size = 2, show.legend = F) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("white", "black")) 

tree_plot2 <- tree_plot + ggtreeExtra::geom_fruit(geom = geom_bar,
                                    mapping = aes(x = mean_MR, fill = pop),
                                    colour = "black",
                                    size = 0.1,
                                    stat = "identity", 
                                    orientation = "y", 
                                    axis.params = list(axis = "x", text.size = 3), 
                                    offset = 0.3, 
                                    pwidth = 0.4,
                                    show.legend = F)

tree_plot2 + ggtreeExtra::geom_fruit(geom = geom_bar,
                                     mapping = aes(x = mean_EWL, fill = pop),
                                     colour = "black",
                                     size = 0.1,
                                     stat = "identity", 
                                     orientation = "y", 
                                     axis.params = list(axis = "x", text.size = 3), 
                                     offset = 0.3, 
                                     pwidth = 0.4,
                                     show.legend = F)

References

Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science 7, 457–472.

Lighton, J. R. B. (2019). Measuring metabolic rates: A manual for scientists. ed. 2. Oxford University Press.

Mijangos, J. L., Berry, O. F., Pacioni, C. and Georges, A. (2022). dartR v2: An accessible genetic analysis platform for conservation, ecology and agriculture. Methods in Ecology and Evolution 18, 691–699.

Nguyen, L.-T., Schmidt, H. A., Von Haeseler, A. and Minh, B. Q. (2015). IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution 32, 268–274.

Stull, R. B. (2000). Meteorology for scientists and engineers: A technical companion book with ahrens’ meteorology today.

Session Information

R version 4.5.0 (2025-04-11)

Platform: aarch64-apple-darwin20

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

other attached packages: MCMCglmm(v.2.36), coda(v.0.19-4.1), Matrix(v.1.7-3), rstan(v.2.32.7), StanHeaders(v.2.32.10), brms(v.2.22.0), Rcpp(v.1.0.14), raster(v.3.6-32), elevatr(v.0.99.0), sp(v.2.2-0), sf(v.1.0-21), rnaturalearth(v.1.0.1), ggtreeExtra(v.1.18.0), ggtree(v.3.16.0), ape(v.5.8-1), dartR.data(v.1.0.8), adegenet(v.2.1.11), ade4(v.1.7-23), cowplot(v.1.1.3), lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.4), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.3.0), ggplot2(v.3.5.2), tidyverse(v.2.0.0) and knitr(v.1.50)

loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), tensorA(v.0.36.2.1), rstudioapi(v.0.17.1), jsonlite(v.2.0.0), datawizard(v.1.1.0), magrittr(v.2.0.3), nloptr(v.2.2.1), farver(v.2.1.2), rmarkdown(v.2.29), fs(v.1.6.6), vctrs(v.0.6.5), minqa(v.1.2.8), terra(v.1.8-54), htmltools(v.0.5.8.1), haven(v.2.5.5), distributional(v.0.5.0), gridGraphics(v.0.5-1), KernSmooth(v.2.23-26), htmlwidgets(v.1.6.4), plyr(v.1.8.9), TMB(v.1.9.17), igraph(v.2.1.4), mime(v.0.13), lifecycle(v.1.0.4), iterators(v.1.0.14), pkgconfig(v.2.0.3), R6(v.2.6.1), fastmap(v.1.2.0), rbibutils(v.2.3), shiny(v.1.10.0), numDeriv(v.2016.8-1.1), digest(v.0.6.37), aplot(v.0.2.7), colorspace(v.2.1-1), ggnewscale(v.0.5.2), patchwork(v.1.3.1), vegan(v.2.7-1), labeling(v.0.4.3), progressr(v.0.15.1), timechange(v.0.3.0), httr(v.1.4.7), abind(v.1.4-8), mgcv(v.1.9-3), compiler(v.4.5.0), proxy(v.0.4-27), pander(v.0.6.6), bit64(v.4.6.0-1), withr(v.3.0.2), backports(v.1.5.0), inline(v.0.3.21), DBI(v.1.2.3), performance(v.0.14.0), QuickJSR(v.1.8.0), pkgbuild(v.1.4.8), MASS(v.7.3-65), classInt(v.0.4-11), corpcor(v.1.6.10), loo(v.2.8.0), permute(v.0.9-7), tools(v.4.5.0), units(v.0.8-7), httpuv(v.1.6.16), glue(v.1.8.0), nlme(v.3.1-168), promises(v.1.3.3), grid(v.4.5.0), checkmate(v.2.3.2), cluster(v.2.1.8.1), reshape2(v.1.4.4), generics(v.0.1.4), glmmTMB(v.1.1.11), seqinr(v.4.2-36), gtable(v.0.3.6), tzdb(v.0.5.0), class(v.7.3-23), data.table(v.1.17.6), hms(v.1.1.3), foreach(v.1.5.2), pillar(v.1.10.2), vroom(v.1.6.5), yulab.utils(v.0.2.0), posterior(v.1.6.1), later(v.1.4.2), splines(v.4.5.0), treeio(v.1.32.0), lattice(v.0.22-7), bit(v.4.6.0), SparseM(v.1.84-2), tidyselect(v.1.2.1), reformulas(v.0.4.1), gridExtra(v.2.3), stats4(v.4.5.0), xfun(v.0.52), bridgesampling(v.1.1-2), matrixStats(v.1.5.0), stringi(v.1.8.7), boot(v.1.3-31), lazyeval(v.0.2.2), ggfun(v.0.1.9), yaml(v.2.3.10), pacman(v.0.5.1), evaluate(v.1.0.4), codetools(v.0.2-20), ggplotify(v.0.1.2), cli(v.3.6.5), RcppParallel(v.5.1.10), Rdpack(v.2.6.4), xtable(v.1.8-4), ggeffects(v.2.3.0), parallel(v.4.5.0), rstantools(v.2.4.0), bayestestR(v.0.16.0), cubature(v.2.1.4), bayesplot(v.1.13.0), Brobdingnag(v.1.2-9), lme4(v.1.1-37), mvtnorm(v.1.3-3), tidytree(v.0.4.6), scales(v.1.4.0), e1071(v.1.7-16), insight(v.1.3.0), crayon(v.1.5.3) and rlang(v.1.1.6)

--- title: "Supplementary Information" date: today date-format: "Do MMM YYYY" author: - name: Nicholas C Wu email: nicholas.wu.nz@gmail.com orcid: 0000-0002-7130-1279 - name: Tomás Villada-Cadavid email: t.cadavid@westernsydney.edu.au orcid: 0000-0002-1743-0694 - name: Christopher Turbill email: c.turbill@westernsydney.edu.au orcid: 0000-0001-9810-7102 csl: ./bib/the-journal-of-experimental-biology.csl bibliography: ./bib/batresp_ref.bib link-citations: true #- Nicholas C. Wu^[Hawkesbury Institute for the Environment, Western Sydney University, NSW 2753, Australia, nicholas.wu.nz@gmail.com] #- Tomás Villada Cadavid^[Hawkesbury Institute for the Environment, Western Sydney University, NSW 2753, Australia, t.cadavid@westernsydney.edu.au] #- Christopher Turbill^[Hawkesbury Institute for the Environment, Western Sydney University, NSW 2753, Australia, c.turbill@westernsydney.edu.au], ^[School of Science, Western Sydney University, NSW 2753, Australia] format: html: embed-resources: true code_download: yes reference-location: margin toc: true code-fold: true code-tools: true revealjs: theme: - journal - journal_mod.scss linkcolor: "#e12744" fontsize: 95%; --- [GitHub Repository](https://github.com/nicholaswunz/bat-resp): The Quarto code to produce the HTML 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 titled "**Local climate and genetic influence on intraspecific patterns in torpor physiology of a cave-roosting bat**". ```{r setup, include=FALSE} library(knitr) knitr::opts_chunk$set(echo = TRUE, cache = FALSE, tidy = TRUE, message = FALSE, warning = FALSE) options(dplyr.width = Inf, knitr.kable.NA = "", digits = 2) # Load library pacman::p_load(tidyverse, cowplot, dartR.base, ape, ggtree, ggtreeExtra, rnaturalearth, sf, sp, elevatr, raster, brms, rstan, MCMCglmm) # Functions mytheme <- function() { theme_bw() + theme(panel.border = element_rect(fill = NA, colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_blank(), axis.ticks = element_line(colour = "black"), axis.text = element_text(size = 10, colour = "black"), axis.title = element_text(size = 10), axis.title.y = element_text(vjust = 3), axis.title.x = element_text(vjust = -1), panel.background = element_rect(fill = NA), plot.background = element_rect(fill = NA, color = NA), plot.margin = unit(c(0.2, 0.2, 0.2, 0.2), units = , "cm"), legend.background = element_rect(fill = NA, color = NA), legend.box.background = element_rect(fill = NA, color = NA), strip.text.x = element_text(size = 10, color = "black", face = "bold"), strip.background = element_rect(fill = NA, color = NA), plot.title = element_text(size = 12) ) } # set up plot theme addbracket <- function(x){paste("[", x, "]")} # add brackets between words # Set options in Rstan set.seed(10) rstan_options(auto_write = TRUE) # translate to STAN platform for running Bayesian model options(mc.cores = parallel::detectCores()) # detects how many cores available to use # Set directory data_path <- file.path("/Users/nicholaswu/Library/CloudStorage/OneDrive-MurdochUniversity/Projects - Old/LP200100331 - WNS - WSU/M - Bat respirometry/") ``` *** # Supplementary Methods ## Respirometry set-up ![](https://raw.githubusercontent.com/nicholaswunz/bat-resp/refs/heads/main/files/Fig%20S1%20-%20Respirometry%20setup.png) **Fig. S1** | Schematic diagram of the flow-through respirometry set up to measure oxygen consumption, carbon dioxide production, and evaporative water loss of seven torpored *Miniopterus orianae oceanensis* and one control blank chamber. The mass flow controllers were placed outside the temperature-controlled incubator to change the flow rate during the experiment manually. ## Calculations and conversation The rate of oxygen consumption ($\dot{V}\textrm{O}_2$; ml O~2~ min^-1^) was calculated following @Lighton2019 in @eq-vo2: $$ \dot{V}\textrm{O}_2 = \frac{\textrm{FR}^\textrm{'}_{\textrm{i}}[(F_{\textrm{i}}\textrm{O}_{2}-F_{\textrm{e}}\textrm{O}_{2})-F_{\textrm{e}}\textrm{O}_{2}(F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}-F^\textrm{'}_{\textrm{i}}\textrm{CO}_{2})]}{(1-F_{\textrm{e}}\textrm{O}_{2})} $$ {#eq-vo2} where $\textrm{FR}^\textrm{'}_{\textrm{i}}$ is the incurrent dry flow rate corrected for water vapour ($\textrm{FR}^\textrm{'}_{\textrm{i}} = \textrm{FR}_{\textrm{i}} \times \textrm{BP}/(\textrm{BP}-\textrm{WVP})$). $F_{\textrm{i}}\textrm{O}_{2}$ and $F_{\textrm{e}}\textrm{O}_{2}$ is the dry incurrent (baseline) and excurrent O~2~ gas concentration, respectively. $F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}$ and $F^\textrm{'}_{\textrm{i}}\textrm{CO}_{2}$ is the dry” equivalent of the excurrent and incurrent CO2 gas concentration corrected for water vapour, respectively. The rate of carbon dioxide production ($\dot{V}\textrm{CO}_2$; ml CO~2~ min^-1^) was calculated as @eq-vco2: $$ \dot{V}\textrm{CO}_2 = \frac{\textrm{FR}^\textrm{'}_{\textrm{i}}[(F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}-F^\textrm{'}_{\textrm{i}}\textrm{CO}_{2})-F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2}(F_{\textrm{i}}\textrm{O}_{2}-F_{\textrm{e}}\textrm{O}_{2})]}{(1-F^\textrm{'}_{\textrm{e}}\textrm{CO}_{2})} $$ {#eq-vco2} Metabolic rate (MR; kJ h^-1^) was calculated from $\dot{V}\textrm{O}_2$ and $\dot{V}\textrm{CO}_2$ as @eq-mr: $$ \textrm{MR} = 0.06 \times (3.941 \times \dot{V}\textrm{O}_2 + 1.106 \times \dot{V}\textrm{CO}_2) \times 4.184 $$ {#eq-mr} The evaporative water loss (EWL) or $\dot{V}\textrm{H}_2\textrm{O}$ (µg H2O min^-1^) was calculated as @eq-ewl: $$ \dot{V}\textrm{H}_2\textrm{O} = \textrm{FR}_{\textrm{i}} \left( \frac{\textrm{WVP}}{T \times R_{\textrm{w}}} \right) $$ {#eq-ewl} where WVP is the water vapour pressure (kPa) minus the baseline WVP, $T$ is temperature in Kelvin (°C + 273.15), and $R_{\textrm{w}}$ is the gas constant for water vapour (461.5 J kg^−1^ K^−1^). $\dot{V}\textrm{H}_2\textrm{O}$ was converted to mg H~2~O h^-1^ by multiplying by 0.06. Prior to calculating $\dot{V}\textrm{O}_2$, $\dot{V}\textrm{CO}_2$, and $\dot{V}\textrm{H}_2\textrm{O}$, the baseline drift and time lag of the fractional gas concentrations were corrected for. ## Cave microclimate The local air temperature and relative humidity outside and inside the cave was recorded to monitor differences in microclimate conditions available for cave-roosting bats. An external HOBO logger (MX2302A; Onset, MA, USA) mounted with a radiation shield (RS3-B; Onset) was placed outside ~20 m away at 1.5 m height from the cave entrance to log external microclimate at 30 second intervals. An internal HOBO logger (MX2301A; Onset) was placed inside the cave above ground in the same chamber where bats were observed roosting to log cave microclimate at 30 second intervals. Cave surface temperature and relative humidity was measured with an infrared thermometer (TG56 IR 30:1; FLIR, OR, USA) and HOBO logger on an 8 m extendable pole (MX2302A) every 10 m inside the cave or where bats were found. Cave surface temperature were measured in three replicates at three time points, May, July, and September. Vapour pressure deficit (kPa) was calculated as the difference between the saturation water vapour pressure ($e_s$) and the actual water vapour pressure ($e_a$). $e_s$ was calculated following the Clausius-Clapeyron equation [@Stull2000]. ```{r Fig S2, fig.align='center', fig.height=5, fig.width=9, echo = FALSE} temp_plot <- read.csv("https://raw.githubusercontent.com/nicholaswunz/bat-resp/refs/heads/main/data/cave_survey.csv") %>% ggplot(aes(x = suf_temp, fill = bats)) + geom_histogram(colour = NA, position = "identity", binwidth = 0.3) + geom_vline(xintercept = c(5, 8, 11, 15), linetype = "dashed") + xlab("Cave surface temperature (°C)") + ylab("Count") + scale_fill_manual(values= c("lightgrey", "black")) + scale_x_continuous(expand = c(0, 0), limits = c(0, 20)) + scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) + guides(fill = guide_legend(title = "Bats present?")) + mytheme() + facet_wrap(~ region, nrow = 2) RH_plot <- read.csv("https://raw.githubusercontent.com/nicholaswunz/bat-resp/refs/heads/main/data/cave_survey.csv") %>% ggplot(aes(x = RH, fill = bats)) + geom_histogram(colour = NA, position = "identity", binwidth = 1.4) + xlab("Cave relative humidity (%)") + ylab("Count") + scale_fill_manual(values= c("lightgrey", "black")) + scale_x_continuous(expand = c(0, 0), limits = c(0, 100)) + scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) + guides(fill = guide_legend(title = "Bats present?")) + mytheme() + facet_wrap(~ region, nrow = 2) VPD_plot <- read.csv("https://raw.githubusercontent.com/nicholaswunz/bat-resp/refs/heads/main/data/cave_survey.csv") %>% dplyr::mutate(es_kPa = 0.611 * exp(2500000 / 461.5 * (1 / 273 - 1 / (temp_C + 273.15))), ea_kPa = RH * es_kPa / 100, VPD_kPa = es_kPa - ea_kPa) %>% ggplot(aes(x = VPD_kPa, fill = bats)) + geom_histogram(colour = NA, position = "identity", binwidth = 0.03) + xlab(" Cave VPD (kPa)") + ylab("Count") + scale_fill_manual(values= c("lightgrey", "black")) + scale_x_continuous(expand = c(0, 0), limits = c(0, 1.8)) + scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) + guides(fill = guide_legend(title = "Bats present?")) + mytheme() + facet_wrap(~ region, nrow = 2) cowplot::plot_grid(temp_plot + theme(legend.position = "none"), RH_plot + theme(legend.position = "none"), VPD_plot + theme(legend.position = "none"), ncol = 3, labels = c('a', 'b', 'c')) ``` **Fig. S2** | Distribution of measured (**a**) cave surface temperatures (°C), (**b**) cave relative humidity (RH; %), and (**c**) calculated cave vapour pressure deficit (VPD; kPa) from Jenolan and Yessabah during the 2023 Austral winter period (May–September). Surface temperature, RH, and VPD with bats present were coloured <span style="color:#000000;">**black**</span>, while surface temperature, RH, and VPD without bats present were coloured <span style="color:#999999;">**grey**</span>. Dashed lines represented the experimental exposure temperatures which are within the expected temperature range in the caves. # Process Data Metadata of the imported respirometry data for analysis. ```{r metadata, echo=FALSE, message=FALSE} data.frame(Name = c("pop", "date_exp", "ID ", "ID_gene","sex", "FA_mm", "mass_g", "temp_exp", "temp_act", "RH_treatment", "chamber_n", "setFR_ml_min", "BP_kPa", "FR_adj_ml_min", "T_start", "T_end", "t_tor_h", "VO2_ml_min", "VCO2_ml_h", "EWL_mg_h", "Tsk_C","VO2_ml_h", "VCO2_ml_h", "VO2_ml_g_h", "VCO2_ml_g_h", "EWL_mg_g_h", "lnMass", "MR_kJ_h", "MR_J_h", "MR_W", "MR_mW", "t_diff"), Description = c("Unique identifier of population source (Jenolan, Yessabah).", "The date when the experiment was conducted in day/month/year", "Unique identifier of the individual.", "Unique identifier of the biopsy sample for genetic analysis", "Sex of the individual (female, male)", "Forearm length in millimetre (mm).", "Body mass in grams (g).", "Set air temperature exposure in degrees Celcius (°C).", "Recorded air temperature in degrees Celcius (°C).", "Recorded relative humidity in percentage (%).", "Unique identifier of the metabolic chamber", "Set flow rate (ml/min).", "Recorded barometric pressure (kPa).", "Adjusted-flow rate (ml/min).", "Start time of the respirometry measurments as year-month-day hour:min:sec.", "End time of the respirometry measurments as year-month-day hour:min:sec.", "Time between onset of torpor and respirometry measurments.", "Calculated rate of oxygen consumption in millilitres per minute (ml/min).", "Calculated rate of carbon dioxide production in millilitres per minute (ml/min).", "Calculated rate of evaporative water loss in milligrams per minute (mg/h).", "Calculated skin surface temperature in degrees Celcius (°C).", "Calculated rate of oxygen consumption in millilitres per hour (ml/h).", "Calculated rate of carbon dioxide production in millilitres per hour (ml/h).", "Mass-specific oxygen consumption in ml/g/h.", "Mass-specific carbon dioxide production in ml/g/h.", "Mass-specific evaporative water loss in mg/g/h.", "The natural logarithm of body mass.", "Calculated metabolic rate in kilojoules per hour (kJ/h).", "Calculated metabolic rate in joules per hour (J/h).", "Calculated metabolic rate in Watts (W).", "Calculated metabolic rate in milliWatts (mW).", "Difference between skin temperature and actual air temperature in degrees Celcius (°C).")) %>% knitr::kable() ``` ## Load and clean data Load and clean the raw data for analysis. **Table S1** summarises data on torpor physiology traits by site, sex, and temperature exposure. ```{r clean, message=FALSE} # Respirometry data resp_dat <- read.csv(file.path(data_path, "data/resp_raw_dat.csv")) %>% dplyr::select(-c(file_name, FR_ml_min, transmitter_n, freq_hz)) %>% dplyr::filter(cond == "dry" & ID != c("MAM-18", "MAM-16")) %>% # remove individuals not in torpor dplyr::mutate(pop = factor(pop), VO2_ml_h = VO2_ml_min * 60, VCO2_ml_h = VCO2_ml_min * 60, VO2_ml_g_h = VO2_ml_h / mass_g, VCO2_ml_g_h = VCO2_ml_h / mass_g, EWL_mg_g_h = EWL_mg_h / mass_g, lnMass = log(mass_g), T_start = lubridate::dmy_hm(T_start), T_end = lubridate::dmy_hm(T_end), MR_kJ_h = (0.06 * (3.941 * VO2_ml_min + 1.106 * VCO2_ml_min)) * 4.184, # Convert VO2 and VCO2 to MR kJ h MR_J_h = MR_kJ_h * 1000, MR_W = MR_kJ_h * 0.2777777778, # convert to Watt MR_mW = MR_W * 1000, lnVO2 = log(VO2_ml_h), lnVCO2 = log(VCO2_ml_h), lnMR = log(MR_J_h), lnEWL = log(EWL_mg_h), t_diff = Tsk_C - temp_act) # filter torpor data only torpor_dat <- resp_dat %>% dplyr::filter(Tsk_C < 28) resp_sum <- torpor_dat %>% dplyr::group_by(pop, temp_exp) %>% dplyr::summarise(n = length(VO2_ml_h), EWL_mean = mean(EWL_mg_h, na.rm = T), EWL_sd = sd(EWL_mg_h, na.rm = T), Tsk_mean = mean(Tsk_C, na.rm = T), Tsk_sd = sd(Tsk_C, na.rm = T), MR_mean = mean(MR_kJ_h, na.rm = T), MR_sd = sd(MR_kJ_h, na.rm = T) ) ``` **Table S1** | Summary data from the respirometry experiment. Sample size (*n*), mean body mass (g), mean ± standard deviation rate of oxygen consumption ($\dot{V}\textrm{O}_2$, ml h^-1^), rate of carbon dioxide production ($\dot{V}\textrm{CO}_2$, ml h^-1^), torpor metabolic rate (TMR, J h^-1^), torpor evaporative water loss (TEWL, mg h^-1^), the skin surface temperature (*T*~sk~, °C), and *Pd* detection via qPCR by sex, air temperature exposure (*T*~a~, °C), and sites. ```{r tableS1, echo = FALSE} tableS1 <- torpor_dat %>% group_by(pop, temp_exp, sex) %>% dplyr::summarise(n = length(VO2_ml_h), mass = mean(mass_g, na.rm = T), O2_mean = mean(VO2_ml_h, na.rm = T), O2_sd = sd(VO2_ml_h, na.rm = T), CO2_mean = mean(VCO2_ml_h, na.rm = T), CO2_sd = sd(VCO2_ml_h, na.rm = T), MR_mean = mean(MR_J_h, na.rm = T), MR_sd = sd(MR_J_h, na.rm = T), EWL_mean = mean(EWL_mg_h, na.rm = T), EWL_sd = sd(EWL_mg_h, na.rm = T), Tsk_mean = mean(Tsk_C, na.rm = T), Tsk_sd = sd(Tsk_C, na.rm = T) ) %>% dplyr::mutate(mass = format(round(mass, 2), nsmall = 2), O2_mean = format(round(O2_mean, 2), nsmall = 2), O2_sd = format(round(O2_sd, 2), nsmall = 2), CO2_mean = format(round(CO2_mean, 2), nsmall = 2), CO2_sd = format(round(CO2_sd, 2), nsmall = 2), MR_mean = format(round(MR_mean, 2), nsmall = 2), MR_sd = format(round(MR_sd, 2), nsmall = 2), EWL_mean = format(round(EWL_mean, 2), nsmall = 2), EWL_sd = format(round(EWL_sd, 2), nsmall = 2), Tsk_mean = format(round(Tsk_mean, 2), nsmall = 2), Tsk_sd = format(round(Tsk_sd, 2), nsmall = 2), O2 = paste(O2_mean, O2_sd, sep = "±"), CO2 = paste(CO2_mean, CO2_sd, sep = "±"), MR = paste(MR_mean, MR_sd, sep = "±"), EWL = paste(EWL_mean, EWL_sd, sep = "±"), Tsk = paste(Tsk_mean, Tsk_sd, sep = "±")) %>% dplyr::select(pop, temp_exp, n, sex, mass, O2, CO2, MR, EWL, Tsk) tableS1 %>% dplyr::rename( "Site" = pop, "*T*~a~ (°C)" = temp_exp, "*n*" = n, "Sex" = sex, "Mass (g)" = mass, "$\\dot{V}\\textrm{O}_2$ (ml h^-1^)" = O2, "$\\dot{V}\\textrm{CO}_2$ (ml h^-1^)" = CO2, "TMR (J h^-1^)" = MR, "TEWL (mg h^-1^)" = EWL, "*T*~sk~ (°C)" = Tsk) %>% knitr::kable(digits = 2) ``` ## SNPs and genetic relatedness The following code was written with Ian Brennan to process, filter and trim the raw SNP data (containing 46,224 loci, 4,075 nucleotide sites) using the [`dartR`](https://green-striped-gecko.github.io/dartR/) package [@dartR2022]. Note: the raw SNP data are being used as part of a larger study examining the genetic connectivity and population dynamics of Australian bent-winged bats. The SNP data will be available upon it's publication (Planckaert *et al.*, In preparation). ```{r SNP, message=FALSE, eval=FALSE} # Load DArT file # Visually check SNP data dartR.base::gl.smearplot(bat_gl_wu) nLoc(bat_gl_wu) # 46224 loci nInd(bat_gl_wu) # 76 individuals ## FILTERING ## ------------------------------------------ # Check raw data dartR.base::gl.report.callrate(bat_gl_wu) # loci callrate dartR.base::gl.report.callrate(bat_gl_wu, method = 'ind') # individual callrate dartR.base::gl.report.reproducibility(bat_gl_wu) # reproducibility # strict threshold on repeatability bat_gl2 <- dartR.base::gl.filter.reproducibility(bat_gl_wu, threshold = 0.999) # remove invariant (monomorphic) sites, we will have to do this # again later, once we've filtered down the dataset and potentially # caused new monomorphs bat_gl3 <- dartR.base::gl.filter.monomorphs(bat_gl2) # adjust threshold on callrate as necessary dartR.base::gl.report.callrate(bat_gl3) # loci callrate bat_gl4 <- dartR.base::gl.filter.callrate(bat_gl3, method = "loc", threshold = 0.9) # filter out multiple SNPs from the same (linked) locus bat_gl5 <- dartR.base::gl.filter.secondaries(bat_gl4, method = "best", v = 2) # filter by the minor allele frequency, I wouldn't exceed 0.05 # in small datasets you worry about losing loci, so can be ~0.04 # in larger datasets should probably be closer to ~0.01 bat_gl6 <- dartR.base::gl.filter.maf(bat_gl5, threshold = (2 / nInd(bat_gl5)), v = 5) # one more pass to clean up callrate and remove monomorphs introduced dartR.base::gl.report.callrate(bat_gl6, method = 'ind') # individual callrate bat_gl7 <- dartR.base::gl.filter.callrate(bat_gl6, method = "ind", threshold = 0.95, mono.rm = T, recalc = T, recursive = T, v = 5) # Export a Concatenated SNP Alignment # Export with amiguities (A/G = R; C/T=Y) for heterozygous sites dartR.base::gl2fasta(bat_gl7, method = 4, outfile = "bat_gl7_mth4.fas", outpath = "local_file", verbose = 5) # Methods 2 and 4 randomly resolve heterozygous sites # read in the alignment as a DNAbin object via ape seq.tmp <- seqinr::read.fasta("local_file", seqtype = "DNA", forceDNAtolower = F) for (i in 1:length(seq.tmp)){ # remove any empty entries in the alignment seq.tmp[[i]] <- seq.tmp[[i]][-which(seq.tmp[[i]] == " ")] } # make a dictionary for ambiguity codes amb.codes <- list(M = c("A","C"),R = c("A","G"), W = c("A","T"),S = c("C","G"),Y = c("C","T"),K = c("G","T")) # loop through and replace any ambiguous bases by randomly choosing one or the other appropriate base for (i in 1:length(seq.tmp)){ seq.trash <- seq.tmp[[i]] rep.trash <- which(seq.trash %in% names(amb.codes)) for (j in rep.trash){ seq.trash[j] <- sample(amb.codes[[which(names(amb.codes) == seq.trash[j])]], 1) } seq.tmp[[i]] <- seq.trash } # By randomly resolving the heterozygous positions, some sites in the alignment # are rendered invariable. Let's get rid of them. We'll do this manually (not using gl.filter.monomorphs) # because our sequence data are in DNAbin format, not genlight (GL). del.pos <- c() for (i in 1:length(seq.tmp[[1]])){ pos.trash <- unlist(purrr::map(seq.tmp,i)) length.trash <- length(unique(pos.trash)) if(length.trash == 1){ del.pos <- c(del.pos,i) } else { if(length.trash ==2 & "N" %in% pos.trash){ del.pos <- c(del.pos,i) } } } for(i in 1:length(seq.tmp)){ seq.tmp[[i]] <- seq.tmp[[i]][-del.pos] } seqinr::write.fasta(seq.tmp, names = names(seq.tmp), nbchar = length(seq.tmp[[1]]), file.out = "local_file") ``` The resulting fasta file was processed through the [IQ-TREE](http://www.iqtree.org/) software to construct the genetic relatedness tree with ascertainment bias correction [@Nguyen2015]. ## Prepare tree for analysis ```{r tree} phylo_tree <- ape::read.tree("http://raw.githubusercontent.com/nicholaswunz/bat-resp/refs/heads/main/files/bat_gl7_mth3.treefile") tree_tip_label <- phylo_tree$tip.label # extract tree tip names ID_list <- unique(torpor_dat$ID_gene) # extract individual name from respirometry dataset #as.data.frame(setdiff(tree_tip_label, ID_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, ID_list)) # prune phylo_tree to keep species from the main dataset # check if lengths match for both data and tree #length(unique(pruned_tree$tip.label)) # 27 #length(unique(torpor_dat$ID_gene)) # 27 # 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 covariance matrix for analysis phylo_cor <- ape::vcv(pruned_tree_cal, cor = F) ``` # Analysis Run phylogenetic multi-level model via [`brm`](https://www.rdocumentation.org/packages/brms/versions/0.6.0/topics/brm) function. We constructed four chains with 5,000 steps per chain, including 2,500-step warm-up periods, hence a total of 10,000 steps were retained to estimate posterior distributions [i.e., (5,000 − 2,500) × 4 = 10,000]. An alpha delta was .99 to decrease the number of divergent transitions. Fixed effects were assigned weakly informative priors following a Gaussian distribution to speed up model convergence, and Student’s t prior with three degrees of freedom was used for group-level, hierarchical effects. The degree of convergence of model was deemed to be achieved when the Gelman–Rubin statistic [@Gelman1992], was 1. ```{r mod, message=FALSE, cache=TRUE, results="hide"} Tsk_mod <- brms::brm(Tsk_C ~ temp_exp * pop + sex + lnMass + tor_start + (1 | ID) + (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor), chains = 4, cores = 4, iter = 1e4, warmup = 5e3, control = list(adapt_delta = 0.99)) MR_mod <- brms::brm(lnMR ~ temp_exp * pop + sex + lnMass + tor_start + (1 | ID) + (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor), chains = 4, cores = 4, iter = 1e4, warmup = 5e3, control = list(adapt_delta = 0.99)) MR_mod2 <- brms::brm(lnMR ~ temp_exp * pop + sex + lnMass + Tsk_C + tor_start + (1 | ID) + (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor), chains = 4, cores = 4, iter = 1e4, warmup = 5e3, control = list(adapt_delta = 0.99)) EWL_mod <- brms::brm(lnEWL ~ temp_exp * pop + sex + lnMass + tor_start + (1 | ID) + (1 | chamber_n) + (1 | gr(ID_gene, cov = phylo)), data = torpor_dat, data2 = list(phylo = phylo_cor), chains = 4, cores = 4, iter = 1e4, warmup = 5e3, control = list(adapt_delta = 0.99)) ``` ```{r Fig S3, echo = FALSE, fig.width = 9, fig.height = 3, fig.align = 'center'} bayesplot::color_scheme_set("darkgray") Tsk_mod_pp <- brms::pp_check(Tsk_mod, type = "scatter_avg") + mytheme() MR_mod_pp <- brms::pp_check(MR_mod, type = "scatter_avg") + mytheme() EWL_mod_pp <- brms::pp_check(EWL_mod, type = "scatter_avg") + mytheme() cowplot::plot_grid(Tsk_mod_pp, MR_mod_pp, EWL_mod_pp, ncol = 3, labels = c('a', 'b', 'c', 'd')) ``` **Fig. S3** | Scatterplots of the observed data (y) vs the average simulated data (y~rep~) from the posterior predictive distribution for the (**a**) *T*~sk~ model, (**b**) the TMR model, and (**c**) the TEWL model. Dashed line represents a slope of 1. ## Model output ::: {.panel-tabset .nav-pills} ### Table S2 - *T*~sk~ model **Table S2** | Summary statistics testing the effect of air temperature on the bat skin surface temperature (*T*~sk~). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%), which includes the intercept ($\beta_0$), exposure temperature (*T*~a~), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass, the time since onset of torpor, and the interaction between *T*~a~ and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is `r as.character(round(data.frame(performance::r2_bayes(Tsk_mod))[2,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(Tsk_mod))[2,2], 4) * 100)`%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is `r as.character(round(data.frame(performance::r2_bayes(Tsk_mod))[1,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(Tsk_mod))[1,2], 4) * 100)`%. ```{r tableS2, echo = FALSE} # Fixed effect fef <- brms::fixef(Tsk_mod) %>% as.data.frame(.) %>% tibble::rownames_to_column("Parameter") %>% dplyr::mutate(Parameter = str_replace(Parameter, "M", "-")) %>% dplyr::mutate(Parameter = case_when( Parameter == "Intercept" ~ "$\\beta_0$", Parameter == "temp_exp" ~ "*T*~a~", Parameter == "popYessabah" ~ "Population - Yessabah", Parameter == "sexmale" ~ "Sex - Male", Parameter == "ln-ass" ~ "lnMass", Parameter == "tor_start" ~ "Torpor onset", Parameter == "temp_exp:popYessabah" ~ "*T*~a~: Population - Yessabah", )) %>% tibble::add_row(Parameter = "**Fixed effects**", .before = 1) # Random effect ref <- summary(Tsk_mod)$random %>% bind_rows() %>% # unlist tibble::rownames_to_column("Parameter") %>% dplyr::select(1:5) %>% dplyr::mutate(Parameter = c("$\\sigma_{chamber}^2$", "$\\sigma_{ID}^2$", "$\\sigma_{phylogeny}^2$")) %>% dplyr::rename(Q2.5 = 4, Q97.5 = 5) %>% tibble::add_row(Parameter = "**Group-level effects**", .before = 1) # Heritability Tsk_h2 <- hypothesis(Tsk_mod, "sd_ID_gene__Intercept^2 / (sd_ID_gene__Intercept^2 + sigma^2) = 0", class = NULL)$hypothesis %>% dplyr::select(-c(1,6:8)) %>% dplyr::mutate(Parameter = "$h^2$") %>% dplyr::rename(Q2.5 = CI.Lower, Q97.5 = CI.Upper) %>% tibble::add_row(Parameter = "**Heritability**", .before = 1) # Render table bind_rows(fef, ref, Tsk_h2) %>% remove_rownames() %>% knitr::kable(digits = 2) ``` ### Table S3 - TMR model **Table S3** | Summary statistics testing the effect of air temperature on torpor metabolic rate (TMR). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%), which includes the intercept ($\beta_0$), exposure temperature (*T*~a~), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass, the time since onset of torpor, and the interaction between *T*~a~ and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[2,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[2,2], 4) * 100)`%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[1,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[1,2], 4) * 100)`%. ```{r tableS3, echo = FALSE} # Fixed effect fef <- brms::fixef(MR_mod) %>% as.data.frame(.) %>% tibble::rownames_to_column("Parameter") %>% dplyr::mutate(Parameter = str_replace(Parameter, "M", "-")) %>% dplyr::mutate(Parameter = case_when( Parameter == "Intercept" ~ "$\\beta_0$", Parameter == "temp_exp" ~ "*T*~a~", Parameter == "popYessabah" ~ "Population - Yessabah", Parameter == "sexmale" ~ "Sex - Male", Parameter == "ln-ass" ~ "lnMass", Parameter == "tor_start" ~ "Torpor onset", Parameter == "temp_exp:popYessabah" ~ "*T*~a~: Population - Yessabah", )) %>% tibble::add_row(Parameter = "**Fixed effects**", .before = 1) # Random effect ref <- summary(MR_mod)$random %>% bind_rows() %>% # unlist tibble::rownames_to_column("Parameter") %>% dplyr::select(1:5) %>% dplyr::mutate(Parameter = c("$\\sigma_{chamber}^2$", "$\\sigma_{ID}^2$", "$\\sigma_{phylogeny}^2$")) %>% dplyr::rename(Q2.5 = 4, Q97.5 = 5) %>% tibble::add_row(Parameter = "**Group-level effects**", .before = 1) # Heritability MR_h2 <- hypothesis(MR_mod, "sd_ID_gene__Intercept^2 / (sd_ID_gene__Intercept^2 + sigma^2) = 0", class = NULL)$hypothesis %>% dplyr::select(-c(1,6:8)) %>% dplyr::mutate(Parameter = "$h^2$") %>% dplyr::rename(Q2.5 = CI.Lower, Q97.5 = CI.Upper) %>% tibble::add_row(Parameter = "**Heritability**", .before = 1) # Render table bind_rows(fef, ref, MR_h2) %>% remove_rownames() %>% knitr::kable(digits = 2) ``` ### Table S4 - TEWL model **Table S4** | Summary statistics testing the effect of air temperature on torpor evaporative water loss (TEWL). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%), which includes the intercept ($\beta_0$), exposure temperature (*T*~a~), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass,the time since onset of torpor, and the interaction between *T*~a~ and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is `r as.character(round(data.frame(performance::r2_bayes(EWL_mod))[2,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(EWL_mod))[2,2], 4) * 100)`%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is `r as.character(round(data.frame(performance::r2_bayes(EWL_mod))[1,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(EWL_mod))[1,2], 4) * 100)`%. ```{r tableS4, echo = FALSE} # Fixed effect fef <- brms::fixef(EWL_mod) %>% as.data.frame(.) %>% tibble::rownames_to_column("Parameter") %>% dplyr::mutate(Parameter = str_replace(Parameter, "M", "-")) %>% dplyr::mutate(Parameter = case_when( Parameter == "Intercept" ~ "$\\beta_0$", Parameter == "temp_exp" ~ "*T*~a~", Parameter == "popYessabah" ~ "Population - Yessabah", Parameter == "sexmale" ~ "Sex - Male", Parameter == "ln-ass" ~ "lnMass", Parameter == "tor_start" ~ "Torpor onset", Parameter == "temp_exp:popYessabah" ~ "*T*~a~: Population - Yessabah", )) %>% tibble::add_row(Parameter = "**Fixed effects**", .before = 1) # Random effect ref <- summary(EWL_mod)$random %>% bind_rows() %>% # unlist tibble::rownames_to_column("Parameter") %>% dplyr::select(1:5) %>% dplyr::mutate(Parameter = c("$\\sigma_{chamber}^2$", "$\\sigma_{ID}^2$", "$\\sigma_{phylogeny}^2$")) %>% dplyr::rename(Q2.5 = 4, Q97.5 = 5) %>% tibble::add_row(Parameter = "**Group-level effects**", .before = 1) # Heritability EWL_h2 <- hypothesis(EWL_mod, "sd_ID_gene__Intercept^2 / (sd_ID_gene__Intercept^2 + sigma^2) = 0", class = NULL)$hypothesis %>% dplyr::select(-c(1,6:8)) %>% dplyr::mutate(Parameter = "$h^2$") %>% dplyr::rename(Q2.5 = CI.Lower, Q97.5 = CI.Upper) %>% tibble::add_row(Parameter = "**Heritability**", .before = 1) # Render table bind_rows(fef, ref, EWL_h2) %>% remove_rownames() %>% knitr::kable(digits = 2) ``` ### Table S5 - TMR + *T*~sk~ model **Table S5** | Summary statistics testing the effect of air temperature and skin surface temperature on torpor metabolic rate (TMR). Mean parameter estimates, estimate error, and 95% Bayesian credible intervals (lower 2.5% and upper 97.5%). The model includes the intercept ($\beta_0$), exposure temperature (*T*~a~), population (Jenolan, Yessabah), sex (female, male), the logarithm-transformed body mass, *T*~sk~, the time since onset of torpor, and the interaction between *T*~a~ and population. Group-level effects include the standard deviations ($\sigma$) for chamber-level effects ($\sigma_{chamber}^2$), individual-level effects ($\sigma_{ID}^2$), and phylogenetic relatedness ($\sigma_{phylogeny}^2$). Heritability was calculated as $h^2 = \sigma_{phylogeny}^2 / (\sigma_{phylogeny}^2 + \sigma^2)$. The $R_{marginal}^2$ representing variance explained by fixed effects only is `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[2,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[2,2], 4) * 100)`%, while the $R_{conditional}^2$ representing variance explained by both fixed and group-level effects is `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[1,1], 4) * 100)` ± `r as.character(round(data.frame(performance::r2_bayes(MR_mod))[1,2], 4) * 100)`%. ```{r tableS5, echo = FALSE} # Fixed effect fef <- brms::fixef(MR_mod2) %>% as.data.frame(.) %>% tibble::rownames_to_column("Parameter") %>% dplyr::mutate(Parameter = str_replace(Parameter, "M", "-")) %>% dplyr::mutate(Parameter = case_when( Parameter == "Intercept" ~ "$\\beta_0$", Parameter == "temp_exp" ~ "*T*~a~", Parameter == "popYessabah" ~ "Population - Yessabah", Parameter == "sexmale" ~ "Sex - Male", Parameter == "ln-ass" ~ "lnMass", Parameter == "Tsk_C" ~ "*T*~sk~", Parameter == "tor_start" ~ "Torpor onset", Parameter == "temp_exp:popYessabah" ~ "*T*~a~: Population - Yessabah", )) %>% tibble::add_row(Parameter = "**Fixed effects**", .before = 1) # Random effect ref <- summary(MR_mod2)$random %>% bind_rows() %>% # unlist tibble::rownames_to_column("Parameter") %>% dplyr::select(1:5) %>% dplyr::mutate(Parameter = c("$\\sigma_{chamber}^2$", "$\\sigma_{ID}^2$", "$\\sigma_{phylogeny}^2$")) %>% dplyr::rename(Q2.5 = 4, Q97.5 = 5) %>% tibble::add_row(Parameter = "**Group-level effects**", .before = 1) #Heritability MR2_h2 <- hypothesis(MR_mod2, "sd_ID_gene__Intercept^2 / (sd_ID_gene__Intercept^2 + sigma^2) = 0", class = NULL)$hypothesis %>% dplyr::select(-c(1,6:8)) %>% dplyr::mutate(Parameter = "$h^2$") %>% dplyr::rename(Q2.5 = CI.Lower, Q97.5 = CI.Upper) %>% tibble::add_row(Parameter = "**Heritability**", .before = 1) # Render table bind_rows(fef, ref, MR2_h2) %>% remove_rownames() %>% knitr::kable(digits = 2) ``` ::: *** # Figures Code to produce the main document figures is provided below. Figures produced were further modified in [Adobe Illustrator](https://www.adobe.com/au/products/illustrator.html) for publication. ## Figure 1 - Study sites Download elevation map from `rnaturalearth` and convert to dataframe. ```{r map, cache=TRUE, results="hide"} # Australia map aus_map <- rnaturalearth::ne_states(country = 'Australia', returnclass = 'sf') # Elevation map aus_elev_sp <- elevatr::get_elev_raster(locations = aus_map, z = 7, clip = "locations") # 1km resolution or 30 arc sec #aus_elev_sp <- readRDS(file.path(data_path, "/aus_elev_sp.rds")) aus_elev_df <- raster::as.data.frame(aus_elev_sp, xy = TRUE) colnames(aus_elev_df)[3] <- "elevation" aus_elev_df <- aus_elev_df[complete.cases(aus_elev_df),] # Remove NAs aus_elev_df_2 <- aus_elev_df %>% dplyr::filter(x >= 143 & y <= -24) ``` Produce eastcoast map with study sites, and temperature profiles. ```{r Fig 1, fig.align='center', fig.height=4, fig.width=8, fig.show='hide'} # Miniopterus distribution range_map <- raster::shapefile(file.path(data_path, "data/Miniopterus_orianae_oceanensis/Miniopterus_orianae_oceanensis.shp")) range_map <- sp::spTransform(range_map, raster::crs("+proj=longlat +datum=WGS84 +no_defs")) # Sites resp_sites_df <- data.frame(x = c(150.018564, 152.688690), y = c(-33.803127, -31.093943), site = c("Jenolan","Yessabah")) map_plot <- ggplot() + geom_polygon(data = range_map, aes(x = long, y = lat, group = group), colour = NA, fill = "#f595a6") + geom_raster(data = aus_elev_df_2, aes(x = x, y = y, fill = elevation), alpha = 0.5, show.legend = F) + scale_fill_gradient2(low = "white", mid = "white", high = "black", midpoint = 0) + geom_sf(data = aus_map, color = "#8C8C8C", fill = NA) + geom_point(data = resp_sites_df, aes(x = x, y = y), shape = 21, colour = "black", fill = "white", size = 3) + xlim(143, 155) + ylim(-39.5, -24) + coord_sf() + theme_void() # Load all files HOBO_list <- list.files(path = file.path(data_path, "data/"), pattern = "\\HOBO.csv$", full.names = T) HOBO_list2 <- data.frame(file_name = HOBO_list, id = as.character(1:length(HOBO_list))) # id for joining HOBO_dat <- HOBO_list %>% lapply(read_csv) %>% # read all files at once dplyr::bind_rows(.id = "id") %>% # bind all tables into one object, and give id for each data.frame() %>% dplyr::mutate(date_time = as.POSIXct(date_time, format = "%m/%d/%Y %H:%M"), date = as.POSIXct(date_time, format = "%m/%d/%Y"), time = strftime(date_time, format = "%H:%M:%S"), jul_day = as.POSIXlt(date_time)$yday, month = format(as.Date(date_time, format = "%m/%d/%Y"),"%m"), month_cat = format(as.Date(date_time, format = "%m/%d/%Y"),"%B"), location = factor(location, levels = c("EXT", "deep"))) temp_plot <- ggplot() + geom_line(data = HOBO_dat %>% dplyr::filter(location == "EXT"), aes(x = jul_day, y = air_temp), size = 0.5, colour = "#f595a6") + geom_line(data = HOBO_dat %>% dplyr::filter(location == "deep"), aes(x = jul_day, y = air_temp), size = 1.2, colour = "black") + ylab("Air temperature (°C)") + xlab(NULL) + mytheme() + theme(legend.position = "bottom") + facet_wrap(~ site, ncol = 1) cowplot::plot_grid(map_plot, temp_plot, ncol = 2, rel_widths = c(1, 0.5)) ``` ## Figure 2 - *T*~sk~, TMR, and TEWL ```{r Fig 2, fig.align='center', fig.height=6, fig.width=8, fig.show='hide'} Tsk_me <- as.data.frame(ggeffects::ggpredict(Tsk_mod, terms = c("temp_exp", "pop", "tor_start [5]", "lnMass [2.7]"))) Tsk_plot <- ggplot() + # raw data geom_line(data = torpor_dat, aes(x = temp_exp, y = Tsk_C, group = ID), colour = "lightgrey", alpha = 0.5, position = position_dodge(0.2)) + geom_point(data = torpor_dat, aes(x = temp_exp, y = Tsk_C, shape = pop), size = 1, position = position_dodge(0.2), alpha = 0.5) + # summary data geom_errorbar(data = Tsk_me, aes(x = x, ymin = conf.low, ymax = conf.high, group = group), width = 0.1, position = position_dodge(0.2)) + geom_line(data = Tsk_me, aes(x = x, y = predicted, group = group, linetype = group), position = position_dodge(0.2)) + geom_point(data = Tsk_me, aes(x = x, y = predicted, shape = group, fill = group), size = 3, position = position_dodge(0.2)) + scale_shape_manual(values = c(21, 24)) + scale_fill_manual(values = c("white", "black")) + xlab(expression(italic("T")["a"]~"(°C)")) + ylab(expression(italic("T")["sk"]~"(°C)")) + mytheme() MR_me <- as.data.frame(ggeffects::ggpredict(MR_mod, terms = c("temp_exp", "pop", "tor_start [5]", "lnMass [2.7]"))) %>% dplyr::mutate(est = exp(predicted), low = exp(conf.low), high = exp(conf.high)) MR_plot <- ggplot() + # raw data geom_line(data = torpor_dat, aes(x = temp_exp, y = MR_J_h, group = ID), colour = "lightgrey", alpha = 0.5, position = position_dodge(0.2)) + geom_point(data = torpor_dat, aes(x = temp_exp, y = MR_J_h, shape = pop), size = 1, position = position_dodge(0.2), alpha = 0.5) + # summary data geom_errorbar(data = MR_me, aes(x = x, ymin = low, ymax = high, group = group), width = 0.1, position = position_dodge(0.2)) + geom_line(data = MR_me, aes(x = x, y = est, group = group, linetype = group), position = position_dodge(0.2)) + geom_point(data = MR_me, aes(x = x, y = est, shape = group, fill = group), size = 3, position = position_dodge(0.2)) + scale_shape_manual(values = c(21, 24)) + scale_fill_manual(values = c("white", "black")) + xlab(expression(italic("T")["a"]~"(°C)")) + ylab(expression("TMR (J h"^"-1"*")")) + scale_y_continuous(breaks = c(1, 10, 25, 50, 100, 200, 400, 800), trans = "sqrt") + mytheme() EWL_me <- as.data.frame(ggeffects::ggpredict(EWL_mod, terms = c("temp_exp", "pop", "tor_start [5]", "lnMass [2.7]"))) %>% dplyr::mutate(est = exp(predicted), low = exp(conf.low), high = exp(conf.high)) EWL_plot <- ggplot() + # raw data geom_line(data = torpor_dat, aes(x = temp_exp, y = EWL_mg_h, group = ID), colour = "lightgrey", alpha = 0.5, position = position_dodge(0.2)) + geom_point(data = torpor_dat, aes(x = temp_exp, y = EWL_mg_h, shape = pop), size = 1, position = position_dodge(0.2), alpha = 0.5) + # summary data geom_errorbar(data = EWL_me, aes(x = x, ymin = low, ymax = high, group = group), width = 0.1, position = position_dodge(0.2)) + geom_line(data = EWL_me, aes(x = x, y = est, group = group, linetype = group), position = position_dodge(0.2)) + geom_point(data = EWL_me, aes(x = x, y = est, shape = group, fill = group), size = 3, position = position_dodge(0.2)) + scale_shape_manual(values = c(21, 24)) + scale_fill_manual(values = c("white", "black")) + xlab(expression(italic("T")["a"]~"(°C)")) + ylab(expression("TEWL (mg H"["2"]*"O"~"h"^"-1"*")")) + scale_y_continuous(breaks = c(1, 2.5, 5, 10, 20, 30), trans = "sqrt") + mytheme() MR_me2 <- as.data.frame(ggeffects::ggpredict(MR_mod2, terms = c("Tsk_C", "pop", "tor_start [5]", "lnMass [2.7]"))) %>% dplyr::mutate(est = exp(predicted), low = exp(conf.low), high = exp(conf.high)) MR_plot2 <- ggplot() + # raw data geom_line(data = torpor_dat, aes(x = Tsk_C, y = MR_J_h, group = ID), colour = "lightgrey", alpha = 0.5, position = position_dodge(0.2)) + geom_point(data = torpor_dat, aes(x = Tsk_C, y = MR_J_h, shape = pop), size = 1, position = position_dodge(0.2), alpha = 0.5) + # summary data geom_ribbon(data = MR_me2, aes(x = x, ymin = low, ymax = high, group = group), alpha = 0.1) + geom_line(data = MR_me2, aes(x = x, y = est, group = group, linetype = group)) + scale_shape_manual(values = c(21, 24)) + scale_fill_manual(values = c("white", "black")) + xlab(expression(italic("T")["sk"]~"(°C)")) + ylab(expression("MR (J h"^"-1"*")")) + mytheme() cowplot::plot_grid(MR_plot + theme(legend.position = "none"), EWL_plot + theme(legend.position = "none"), Tsk_plot + theme(legend.position = "none"), MR_plot2 + theme(legend.position = "none"), nrow = 2, labels = c('a', 'b', 'c', 'd')) ``` ## Figure 3 - Phylogeny ```{r Fig 3, fig.align='center', fig.height=5, fig.width=8, fig.show='hide'} MR_mod_post <- posterior_samples(MR_mod) # extract posterior EWL_mod_post <- posterior_samples(EWL_mod) # extract posterior # Loop to extract within ID-level posterior of intercept ID_unique <- sort(unique(torpor_dat$ID)) # pop ID names out_list <- vector(mode = 'list', length = length(ID_unique)) for (j in seq_along(ID_unique)) { out_list[[j]] <- data.frame(matrix(NA, nrow(MR_mod_post), 1)) dat_names <- paste0(ID_unique[j]) names(out_list[[j]]) <- dat_names out_list[[j]][, dat_names[1]] <- MR_mod_post[, paste0('r_ID[', ID_unique[j], ',Intercept]')] } MR_Int_out <- do.call('cbind.data.frame', out_list) MR_Int_post <- sample_n(MR_Int_out, 4000) # Extract 4000 random rows out_list <- vector(mode = 'list', length = length(ID_unique)) for (j in seq_along(ID_unique)) { out_list[[j]] <- data.frame(matrix(NA, nrow(EWL_mod_post), 1)) dat_names <- paste0(ID_unique[j]) names(out_list[[j]]) <- dat_names out_list[[j]][, dat_names[1]] <- EWL_mod_post[, paste0('r_ID[', ID_unique[j], ',Intercept]')] } EWL_Int_out <- do.call('cbind.data.frame', out_list) EWL_Int_post <- sample_n(EWL_Int_out, 4000) # Extract 4000 random rows # Tidy data MR_Int_post_long <- MR_Int_post %>% tidyr::pivot_longer(cols = everything(), names_to = "ID", values_to = "Intercept") EWL_Int_post_long <- EWL_Int_post %>% tidyr::pivot_longer(cols = everything(), names_to = "ID", values_to = "Intercept") MR_Int_post_sum <- MR_Int_post_long %>% dplyr::group_by(ID) %>% dplyr::summarise(mean_MR = mean(Intercept)) EWL_Int_post_sum <- EWL_Int_post_long %>% dplyr::group_by(ID) %>% dplyr::summarise(mean_EWL = mean(Intercept)) phylo_dat <- torpor_dat %>% dplyr::group_by(ID_gene, ID, sex, pop) %>% summarise(mean_mass = mean(mass_g)) x <- tidytree::as_tibble(pruned_tree_cal) y <- dplyr::full_join(phylo_dat, MR_Int_post_sum, by = 'ID') y <- dplyr::full_join(y, EWL_Int_post_sum, by = 'ID') tree_dat <- tidytree::as.treedata(dplyr::full_join(x, y %>% rename(label = ID_gene), by = 'label')) tree_plot <- ggtree::ggtree(tree_dat) + ggtree::geom_tippoint(aes(shape = pop, fill = pop), size = 2, show.legend = F) + scale_shape_manual(values = c(21, 24)) + scale_fill_manual(values = c("white", "black")) tree_plot2 <- tree_plot + ggtreeExtra::geom_fruit(geom = geom_bar, mapping = aes(x = mean_MR, fill = pop), colour = "black", size = 0.1, stat = "identity", orientation = "y", axis.params = list(axis = "x", text.size = 3), offset = 0.3, pwidth = 0.4, show.legend = F) tree_plot2 + ggtreeExtra::geom_fruit(geom = geom_bar, mapping = aes(x = mean_EWL, fill = pop), colour = "black", size = 0.1, stat = "identity", orientation = "y", axis.params = list(axis = "x", text.size = 3), offset = 0.3, pwidth = 0.4, show.legend = F) ``` *** # References {-} <div id="refs"></div> <br> *** ## Session Information ```{r sessioninfo, echo = FALSE} pander::pander(sessionInfo(), locale = FALSE) ```

Supplementary Methods

Respirometry set-up

Calculations and conversation

Cave microclimate

Process Data

Load and clean data

SNPs and genetic relatedness

Prepare tree for analysis

Analysis

Model output

Figures

Figure 1 - Study sites

Figure 2 - Tsk, TMR, and TEWL

Figure 3 - Phylogeny

References

Session Information

Figure 2 - T_sk, TMR, and TEWL