Aim: We investigated the relationship between thermal physiology, elevational distribution, and thermal stress among nine closely related dendrobatid frogs during their aquatic stage by employing an integrated approach, combining thermal physiology, environmental temperature modelling, and predictive assessments of current and future exposure to thermal variation.

Location: Ecuador

Taxon: Amphibians; Anura, Dendrobatidae, Epipedobates, Hyloxalus

Methods: We determined the Thermal Performance Curves of larval growth (TPCs) for each species and modelled the thermal variation in contrasting aquatic larval environments for both present and future times. This allowed us to estimate the expected elevational distributions and forecast periods of exposure to stressful temperatures that inhibit larval growth due to elevation and global warming.

Results: We found significant correlations between optimum temperature (Topt), 50 % maximum performance temperature (maxB50), 50% minimum performance temperature (minB50), and cold resistance (survival at 9 ºC) with the current elevational distributions. However, thermal physiology predicted lower than observed distributions for high-elevation dendrobatids and higher than observed maximum elevations for lowland species. Nonetheless, our models predicted that low thermal variability habitats (i.e. streams and deep permanent ponds) can buffer the future increase in temperatures for all taxa, even when considering the most extreme scenario. In contrast, all species within high thermal variation habitats (open forest temporary ponds) are expected to experience stressful temperatures under present conditions.

Main Conclusions: The findings indicate that thermal physiology may not be a limiting factor for dendrobatid frog species' ranges in this equatorial mountain gradient. Highland species may need to adapt to suboptimal performance, while some lowland species could occupy higher elevations. This study emphasizes the importance of habitat buffering to mitigate thermal stress in the face of climate change for amphibians in tropical mountains.

Selection of Species

We sampled nine Ecuadorian frogs belonging to the family Dendrobatidae during their larval stage, covering a wide elevational range in Ecuador ranging from 38 - 1900 m, between latitudes 1ºN - 4ºS). We focused on five out of the six species of Epipedobates found on the coastal side of Ecuador, as well as four species of Hyloxalus from the Amazonian slope out of the twenty-seven that can be found throughout the country (Ron, et al. 2020; for more information see supplementary Table S1 and Fig. S6). One population of each species was collected from their natural aquatic habitats from December 2014 to April 2016 and transported to the experimental facilities in the Pontificia Universidad Católica del Ecuador (PUCE). However, for three species (Epipedobates machalilla, E. tricolor and Hyloxalus nexipus) we sourced specimens from ‘Balsa de los Sapos’ initiative PUCE. It is important to note that the specimens from captive breeding were first- or second-generation breeds, and we assumed that their physiological performance was not affected by captivity (Pintanel, et al. 2020).

Thermal Performance Curves for Growth and Survival (TPC)

To assess how tadpole growth rate and survival were affected by temperature, we established groups comprising 10-15 individuals. These groups were reared under different constant temperature treatments, including 9, 15, 20, 23.5, 27, 29, 31 and 33 ºC, all while maintaining a photoperiod of 12L:12D over a 10-days period, and with ad libitum food access. Throughout the experiment, we completely replaced the water and food three times, approximately, every 2-3 days, and we daily checked survival of all individuals. To regulate water temperature and maintain it at or above 20ºC (20-33ºC), we utilized portable fluid heaters with adjustable controls (patent license ES1135983U; Madrid, Spain). For the colder treatments at 9 ºC and 15 ºC, we employed TECO TK 1000 chillers. We randomly assigned each tadpole to a specific temperature treatment and individually housed them in plastic cups filled with 400 ml of dechlorinated water. We aerated these containers using an air pump system. We measured their weight and determined their developmental stage according to the (Gosner 1960) scale just before and after the experiment. To calculate the relative growth rate (GRr), we used the following formula:

GRr = (weight_final - weight_initial)/(weight_initial x n_days)

Negative values were considered as non-growth and treated as zero (Overgaard, et al. 2014). Given the non-linear growth pattern of anuran larvae, which approaches exponential decay close to metamorphosis (Harris 1999, Richter-Boix, et al. 2011), we selected experimental tadpoles in early-mid-developmental stages (ranging 25-34 Gosner stages) from either field or captive breeding sources that reached as most 38-39 Gosner stages at the end of the growth experiment.

Estimates of Micro- and Macro-Climatic Temperatures along Elevation

We gathered macroclimate data, including maximum, minimum, and mean temperatures (T_MAX, T_MIN, and T_MEAN), from current and projected future climate conditions using WorldClim, with a spatial resolution of 30 arc-seconds (Hijmans, et al. 2005). To assess shifts in thermal suitability for each species under different climate warming scenarios, we considered future projections for two time frames (2050 and 2070) under both low (RCP 4.5) and high (RCP 8.5) emission scenarios as provided by the CCSM4, HadGEM2-ES, and GFDL-CM3 global circulation models. We used the 'extract' function from the raster R-package (Hijmans 2017b) to extract climatic data from WorldClim at each species' georeferenced sampling locations.

We obtained georeferenced locations for Ecuadorian Epipedobates and Amazonian-side Hyloxalus were obtained from the Museo de Zoología of the Pontificia Universidad Católica del Ecuador database (Ron 2018; see Fig. S6). To address sampling bias, we retained georeferenced points that were separated by more than 1 km, using the ‘distGeo’ function in the geosphere R-package (Hijmans 2017a). We collected in situ microclimatic data from 38 breeding sites, spanning elevations from 23 to 3631 m, representing diverse aquatic habitats with contrasting temperature profiles (see below). We recorded water temperature every 15 minutes using dataloggers (HOBO pendant) deployed at the bottom of each water body. We analysed mean (t_mean), maximum (t_max) and minimum (t_min) daily water temperatures. The number of sampling days varied, ranging from 2 to 546 days depending on the ephemerality and accessibility of the breeding site (see Table S2 in Supporting Information). Finally, we categorized the water bodies into three different thermal variability regimes according to their mean daily thermal range (dr = mean t_max – t_min): low (0.05 – 1 ºC; streams and permanent deep ponds), medium (1 - 2.25 ºC; forested temporary ponds), and high (> 2.39 ºC; temporary opened ponds) (Table S2 in Supporting Information). We selected dr as the discriminating variable for differentiating between thermally contrasting habitats, because it was independent of elevation for all locations and within any thermal range category (R² = 0.02, n = 38, P = 0.891), thus mitigating the confounding effect of decreasing temperatures with elevation due to the adiabatic lapse rate (Dillon, et al. 2006, Sarmiento 1986).

Although most analysed dendrobatid tadpoles are rarely found in open temporary ponds, we included the 'high' category to encompass the full range of environmental temperatures available for frogs and to account for the effects of human-modified tropical forests resulting from anthropic activities such as logging, fragmentation, and habitat conversion (Senior, et al. 2017).

Modelling Microclimatic Temperatures through Species’ Ranges

To gauge the microclimatic temperatures across species distributional ranges, we employed various GLS models to estimate both maximum (t_max) and minimum (t_min) temperatures. These models incorporated T_MAX, T_MIN and T_MEAN, along with their squared values, as response variables, while also considering habitat (refer to Table S6 and S7 in Supporting Information). Following this, we selected the optimal model based on AIC comparison (with a ΔAIC > 3) and extracted coefficients for the intercept and each explanatory variable. These model parameters were then interpolated using current and future temperatures at studied species locations to obtain maximum and minimum mean modeled microenvironmental temperatures (see Supporting Information). Given the low latitude of our sampling locations, which range from 1.165 N to 4.164 S, temperatures exhibit minimal seasonal variation, remaining relatively constant throughout the year. The primary source of temperature variability at our study sites arises from microenvironmental thermal variation (see Fig. S2 in Supporting Information). To estimate realistic daily thermal variation, we interpolated estimated macroclimatic temperatures with field-measured microclimatic temperatures. This interpolation was based on three representative microclimatic datasets with contrasting temperatures, obtained from stream, forested pond, and open pond environments, each collected over 78-100 days (see Fig. S3-5 in Supporting Information).

Relationships of Elevational Distribution with Thermal Sensitivity (TPC) Parameters

Our growth Temperature-Performance Curves (TPCs) construction revealed that the warm extreme of B₅₀ (maxB₅₀), ranging from 24.8ºC and 31.6ºC across species, represents a distinct thermal threshold. Beyond this range, mortality rates increased in TPC experimental treatments between 27ºC and 33ºC (see Fig. 1 and Results). Stressfully cold temperatures were challenging to estimate because at most species, except H. pulchellus, could not survive the coldest treatment at 9ºC over the ten-day experiment. However, the cold extreme of B₅₀ (minB₅₀), ranging between 16.6ºC and 24.3ºC, also represented a critical cold threshold. The number of days tadpoles survived (1 to 10 days) exposed to 9ºC and 15ºC significantly increased at lower minB₅₀ (see Results).

Expected Distributional Ranges and Chronic Thermal Stress across Elevation

Expected elevational ranges per species were determined across habitats, defined as those elevations where potential hours of maximum performance (best growth) exceed 50% of the total hours (see Figures S7-S15 in Supporting Information). Potential hours of best growth per elevation were computed as hours when temperatures ≥ minB₅₀ but ≤ T_opt.

To quantify the impacts of climate change across habitats and throughout the elevational gradient, we projected the present and future percentages of total hours of heat stress (temperatures > maxB₅₀) for each occurrence point. All analyses were conducted using R v3.4.3 (R Core Team 2021).