Compensatory growth partially mediates the effects of stressors associated with climate change on amphibian physiology

Human-induced climate change and urbanization are predicted to dramatically impact landscape hydrology, which can have devastating impacts on aquatic organisms. For amphibians that rely on aquatic environments to breed and develop, it is essential to understand how the larval environment impacts development, condition, and performance later in life. Two important predicted impacts of climate change and urbanization are reduced hydroperiod and increased larval density. Here we explored how larval density and hydroperiod affect development, morphology, physiology, and immunity at metamorphosis and 35 days post metamorphosis in the frog Rana pipiens. We found that high density had a large negative impact on development and morphology, which resulted in longer larval periods, reduced likelihood of metamorphosis, smaller size at metamorphosis, shorter femur-to-body length ratio, and reduced microbiome species evenness compared to animals that developed in low-density conditions. However, animals from the high-density larval conditions experienced compensatory growth post-metamorphosis, demonstrating accelerated growth in body size and relative femur length compared to animals from the low-density treatments. We also saw an increase in relative gut length and relative liver size in animals that had developed in the high-density treatment, as well as increased immune function, and greater jump distances relative to their leg length across different temperatures, which further supports post-metamorphic growth compensation following a stressful developmental period. Finally, metabolic rate was higher overall but especially at higher temperatures for animals that developed under high-density conditions. A higher metabolic rate at warmer temperatures might indicate a reduced ability to acclimate quickly to increasing environmental temperatures.

Methods
Mesocosm design

The experiment took place at the Pymatuning Laboratory of Ecology in Linesville, Pennsylvania, United States. Sixteen large cattle tanks (770L, Rubbermaid, grey plastic) were cleaned with bleach and allowed to dry for 24h, then filled to a depth of 41cm (600L) with well water 9 weeks prior to the start of the experiment. Lids for the mesocosms were constructed using 50% shade cloth and weighted tubing to ensure no insect predators or other amphibians could inhabit the mesocosms. Leaf litter was gathered from a nearby wooded area and dried, and 250g was added to each mesocosm 6 weeks prior to the start of the experiment. Pond water (1L per mesocosm) with concentrated zooplankton caught with a zooplankton net was added to each mesocosm 7 weeks prior to the start of the experiment. Each mesocosm was fitted with two temperature loggers (HOBO temperature pendants, Onset etc.) set to record the water temperature every 30 minutes, with one 5cm above the bottom on the tank, and one 5cm below the water surface of each tank. One day before the start of the experiment we added 15g of rabbit chow, to provide supplementary food for the animals.

Five egg masses were laid on June 10, 2017 from commercially purchased animals (Connecticut Supply) following hormonal injections of leuprolide acetate (Brannelly, Ohmer, & Richards-Zawacki, 2019). Eggs were allowed to develop and hatch. On June 15, 2017, day 0 of the experiment, tadpoles (stage 25, Gosner, 1960) were evenly sourced from these five egg masses and placed inside the mesocosms. Egg masses were produced by 12 captive adult R. pipiens housed communally for breeding purposes. These 12 adults were swabbed for skin microbiome of bacteria and fungi (rinsed with sterile artificial pond water, swabbed with a rayon-tipped swab, MS113, Medical and Wire Co.) on the first day of hormonal induction and one day after arrival to our facility. The swab samples were stored at -20°C until processing.

The starting density for the high-density treatment was 0.12 tadpoles/L, and the starting density of the low-density treatment was 0.06 tadpoles/L. The larval developmental portion of the experiment ended on day 105. On this day the fast-drying regime was due to reach a water depth of 0cm, and the slow-drying treatment reached a depth of 17cm (approximately 249L). Before removing water on this final day we removed and euthanized (in 0.1% buffered tricaine methanesulfonate, MS-222, Sigma Aldrich) the remaining tadpoles (n=32 from the high-density, fast-drying treatment at an average tadpole density of 0.18 tadpoles/L and n=44 from the high-density slow-drying treatment at an average tadpole density of 0.04 tadpoles/L). All tadpoles in the low-density treatments metamorphosed and escaped the mesocosm before the end of the experiment (low-density fast-drying, day 63 and approximately 249L; low-density slow-drying, day 88 and approximately 352L). All animals that metamorphosed before the end of the experiment were included in the metamorphosis analyses (size at metamorphosis, rate of metamorphosis etc). Because the entire experiment ended 105 days after tadpoles entered the mesocosms, only animals that metamorphosed by day 70 were included in the survival analyses and juvenile condition analyses, so that we could ensure that these animals were at least 35 days old at the time of the assessments.

The densities tested in this study are within the range of what *R. pipiens *experiences in the wild. Tadpole densities vary dramatically in the wild depending on pond type (low-density in permanent or semi-permanent ponds, to extremely high-density in drying ephemeral ponds). Of note, we have seen tadpoles at extremely high-density (hundreds of tadpoles in small ephemeral forest pools (estimated <600 L; pers obs), and this is supported by the literature in nearby locations (average density of 0.05 tadpoles/L, with a maximum of 0.5 tadpoles/L; Smith et al., 2003). Therefore, the densities observed in this study are within the observed range for this species in this part of the world, and do not represent extreme densities. Additionally, the densities used in this experiment are slightly less dense but within the range other experiments that have used in mesocosm designes assessing the impacts of rearing density on this species (e.g., low- to high-density of 0.07 - 0.2 tadpoles/L Dare et al., 2006; Distel & Boone, 2011).

Microbiome extraction

Skin swabs were extracted using the DNeasy 96 blood and tissue extraction (69581, Qiagen) with the following modifications during the lysis and incubation steps: 180μL of warmed lysozyme lysis buffer (20mg/mL, Lysozyme, 89833, ThermoFisher Scientific) was added to the samples and incubated at 37°C for 1hr. Then 25μL of proteinase K plus 200μL of buffer AL was added and incubated at 52°C for 30min. then 200μL of 100% ethanol was added, and the rest of the kit protocol was followed. We eluted the samples in 100μL of elution buffer after a 5 min incubation at room temperature for a final extract volume of 100uL. The extracts were frozen and stored until sequencing.

Samples were PCR processed to amplify the v4 region (515f – 806r) of the 16S rRNA gene of bacteria and the ITS-1 and ITS-2 region of the fungal ITS gene following the Earth Microbiome protocol developed by Caporaso et al. (2012). Individual barcoded primers were used for PCR amplification (Parada et al., 2016) and the reaction mixture consisted of 11μL of PCR-grade water, 10μL of PCR master mix (5PRIME HotMasterMix) 1μL of forward and reverse primers and 2μL of DNA sample were added to the 16S PCR reactions. For ITS PCR, 10μL of PCR-grade water, 10μL of PCR master mix (5PRIME HotMasterMix), 1μL of forward and reverse primers, and 3μL of DNA template were added to each PCR reaction. 16S thermocycler conditions consisted of 94°C for 3 minutes and then 35-cycles of 94°C for 45 seconds, 50°C for 60 seconds, and 72°C for 90 seconds, followed by 72°C for 10 minutes and then a 4°C hold. Fungal ITS conditions consisted of 94°C for 3 minutes, and then 35-cycles of 94°C for 30 seconds, 52°C for 30 seconds, and 68°C for 30 seconds, followed by 68°C for 10 minutes and then a 4°C hold. All samples were run in duplicate. Duplicate reactions were then pooled and purified using the Thermo Scientific GeneJET PCR Purification Kit (K0702, Thomas Scientific). Finally, samples were analysed for DNA content using a Nanodrop and standardized to 50ng/μL. Sequencing was performed individually for the 16S and ITS libraries using a MiSeq Reagent Kit v2 (300-cycle) from Illumina.

Bacterial killing ability of the whole blood

The BKA assay was based on Liebl and Martin (2009) and modified from the protocol provided at www.ecoimmunology.org, which is funded by NSF-0947177), Savage et al. (2016), and Brannelly et al (2019). The test used 6µL of whole blood from each individual, and we used *Escherichia coli *in the killing assay. We activated one BactiDisk pellet (8739; ATCC®, Manassas, VA, USA) in 10mL phosphate-buffered saline (PBS) at 30˚C for 30 min then diluted the sample to 5 x 104 bacteria per mL and immediately refrigerated it. We vortexed whole blood for 10 sec and added 3µL to 69µL of a solution sterile CO2-independent media (ThermoFisher, Waltham, MA, USA) containing 4mM L-glutamine (ThermoFisher, Waltham, MA, USA) broth. A second tube of blood and CO2-independent media was prepared for each sample as a sample blank. We added 25µL of *E. coli *to each sample, and 25µL of sterile PBS to each sample blank, incubated the samples for 1h at 21˚C, and then added 500µL of sterile tryptic soy broth (TSB; SigmaAldrich, St. Louis, MO, USA) to each sample. We prepared three positive controls and control blanks for each plate, with positive controls containing 500µL of TSB broth, 25µL *E. coli, *and 72µL PBS and blanks containing 500µL of TSB broth and 92µL of sterile PBS. All samples were then vortexed and incubated at 30˚C for 12h (Savage *et al., *2016). We ran each sample, sample blank, control and control blank in triplicate at a 100µL volume and measured absorbance at 405 nm (using an Epoch BioTek plate reader, Winooski, VT, USA) in a 96 well sterile cell culture plate (Costar® Corning Inc., Corning, NY, USA).

We adjusted for plate inconstancies (delta=405nm absorbance reading - 560nm absorbance reading) and calculated BKA index (Allen et al., 2009):

This index assigns large positive numbers to animals that had high bacterial killing ability (i.e., where many bacteria were killed) and low or negative values when there was little bacterial killing ability within the blood. This index also accounts for bacterial die-off or plate inconsistencies within the assay that is unrelated to the blood’s bactericidal ability (Allen et al., 2009). Individual reaction wells with extreme outliers were removed from the analysis (n=2 reaction wells) if the BKA index value exceeded three standard deviations from the mean. Although two reaction wells were excluded from analyses because they were extreme outliers, all animals (n=137) were included in the analysis.

Jumping performance

Each animal was tested at all six temperatures over 3–4 days (n=156). We strategically chose five different temperature orders to sequentially test the animals. The temperature orders (in ºC) were A) 20, 15, 10, 30, 25, 5, 20, B) 15, 20, 25, 30, 5, 10, 15 C) 5, 15, 20, 25, 30, 10, 5, D) 10, 15, 20, 5, 25, 30, 10, E) 20, 25, 15, 10, 30, 5, 20. We randomly placed the animals in those these temperature groups, and these temperature orders were included as a covariate in the analysis. Each frog was tested seven times total, where the first temperature was repeated at the end of the temperature regime to ensure animals were not fatigued. This jump trial number (1–7) was included as a covariate in the analysis. Animals were included in the analysis if the mean jump distance was at least 70% of the mean jump distance of the first jump temperature (excluded n=60 animals). For each jumping trial, we removed the animals from the chambers and encouraged them to jump by lightly touching the urostyle up to five times. The first three successful jumps per temperature recorded. Trials in which a frog did not jump successfully at the first or second encouragement were not included in the analysis for that jump temperature, nor were those that successfully jumped only once for that jump temperature (n=96 used in final analysis). Body temperature was recorded using an infrared temperature gun just prior to each jump attempt, and this temperature was used in the analysis. The start and end position of the first jump an animal made from rest was marked and measured to the nearest 1mm.

Metabolic rate

Prior to the metabolic rate trials we fasted the animals for 5 days before testing. We calculated rate of oxygen consumption using the equation (Vleck, 1987):

where VO2=rate of oxygen consumption, and V=Volume of air in the chamber (calculated as V=Vchamber – Vpaper towel – Vanimal). We measured mass of the animal and paper towel subtracted and calculated the volume of each (assuming 0.001 L/g density). Fi is the initial dry O2 fractional reading (i.e., the machine output) and Fe is end dry O2 fractional reading. VH2O is the change in the fractional water vapor concentration and assumed to be equal to 0 because the champer it is constant and saturated (damp paper towel was used as substrate).

Microbiome analysis

Raw 16S bacterial sequence data were processed through the QIIME2 pipeline version 2018.4 (Bolyen et al., 2018). Sequences were demultiplexed following the DADA2 pipeline within QIIME2, the sequence (McDonald et al., 2012) reads were filtered, processed and assigned to operational taxonomic features (OTFs). Singleton observable taxonomic features were removed, and a phylogenetic tree using FASTTREE (Prince et al., 2010) was built and taxonomy was classified using the q2-feature-classifier (Bokulich et al., 2018) classify-sklearn naïve Bayes taxonomy classifier against Greengenes reference sequences (McDonald et al., 2012), and sequences identified as chloroplasts or mitochondria were removed from downstream analysis. Observable taxonomic features tables were rarefied to 10,000 reads. 71 individual samples were processed and sequenced, and 8 samples were not included after rarefication due to low observable taxonomic features reads; therefore 63 samples were analysed in downstream analyses. A further two metamorph samples were removed due to missing metadata, for a total of 11 adults and 50 metamorph samples. Alpha diversity metrics were used to measure the bacterial community diversity within each rarefied sample. Alpha diversity measures were calculated within QIIME2: the number of observable taxonomic features (qualitative measure of richness), Shannon diversity (quantitative measure of community richness), and Pielou’s evenness (a measure of community evenness). Beta diversity metric outputs from QIIME2 to compare bacterial community composition were weighted UniFrac distances (community structure, which compares samples on the basis of presence, absence and relative abundance of bacterial taxonomic features), and unweighted UniFrac distances (community membership, which compares samples on the basis of presence and absence of bacterial observable taxonomic features).

The raw ITS fungal sequence was processed through a similar QIIME2 pipeline, where the taxonomy classifier was built against UNITE reference sequences (Abarenkov et al., 2020). Observable taxonomic features tables were rarefied to 600 reads. Ten samples were not included after rarefication because of low observable taxonomic features reads, and one sample was not included due to missing metadata; therefore, 52 metamorphs and 9 adults were analysed.