Changes in physiology and microbial diversity in larval ornate chorus frogs are associated with habitat quality

We examined multiple health measures in ornate chorus frog tadpoles collected from across their range as a function of water quality and land cover characteristics. Water quality, land cover, developed land and location were linked to health indicators, suggesting these measures may aid in identifying populations at increased decline risk.


Introduction
As the human population continues to increase, land-use conversion has altered habitat suitability for many species (Vitousek et al., 1997;McKinney, 2002). Habitat loss, fragmentation and degradation, coupled with climate change and other anthropogenic factors, are among the most significant drivers of population declines and species extinctions (Brook et al., 2008;Mantyka-Pringle et al., 2012;Segan et al., 2016;Ceballos et al., 2017). Landscape-level disturbances such as habitat loss and degradation affect populations by altering habitat connectivity, composition and quality (Fahrig, 2003). For aquatic organisms, local impacts on water quality such as changes in temperature, pH, contaminant, nutrient and sediment levels (Meybeck, 2004) can also pose significant consequences to population health and resilience. Together, the additive or synergistic effects are especially troubling for species, especially those with complex life cycles that require multiple habitats to complete development (Hayes et al., 2010;Blaustein et al., 2011). Examining physiological responses of individuals to such factors can provide insight into the mechanisms by which environmental stressors can lead to population declines (Wikelski and Cooke, 2006). In turn, understanding the consequences of these factors on populations could aid in better management practices for declining species.
Glucocorticoid (GC) hormones can be useful tools to assess the physiological response of organisms to environmental and anthropogenic stressors which, when combined with other metrics (e.g., immune function), can be indicators of individual and population health (Sheriff et al., 2011;Dantzer et al., 2014;Jeffrey et al., 2015). In response to stressors, GC hormones are released by the hypothalamo-pituitary-interrenal (HPI) axis above baseline levels as an adaptive response that assists in energy mobilization, mediates natural changes in physiology and behaviour, kick-starts the immune response and helps return the organism to homeostasis (Dhabhar and McEwen, 1999;Sapolsky et al., 2000;Cyr and Romero, 2009;Romero et al., 2009). However, under prolonged (chronic) stress, vertebrates may have elevated or lowered baseline GCs. In such cases, the HPI axis can become dysregulated (where individuals have a muted response to additional stressors), resulting in reduced physiological health, pathology and ultimately death (Wingfield and Romero, 2001;McEwen and Wingfield, 2003;Cyr and Romero, 2009;Romero et al., 2009). Vertebrates under chronic stress have suppressed growth and development, reduced reproduction and reduced immune function (Sapolsky et al., 2000;McEwen and Wingfield, 2003;Romero, 2004;Rollins-Smith, 2017). Pre-stressor baseline GC levels and the response to acute stressors (HPI axis responsiveness) can be measured to aid in identifying whether populations are under chronic stress, which may be associated with higher risk of declines (Homan et al., 2003;Dantzer et al., 2014;Baugh et al., 2018;. Measuring GC hormones may identify populations to target for management but may not provide a full picture of environmental effects on population health. Using additional metrics such as immune function and microbial diversity is needed for a more holistic view of how environmental stressors impact population health (Breuner et al., 2013).
Amphibians have the highest threat status of all vertebrate classes (Stuart et al., 2004;Foden et al., 2013) and are particularly susceptible to changes in habitat quality. Species that are especially vulnerable use ponds or other water bodies to breed but occupy terrestrial habitats outside the breeding season (Becker et al., 2007). As such, the quality and conditions of both the aquatic and surrounding terrestrial environments play a factor in the health of these amphibians. Environmental changes and conditions associated with lower quality habitat may be perceived as stressors, which elevate or dysregulate GC hormones and affect physiological health (Homan et al., 2003;Janin et al., 2011Janin et al., , 2012Chambers et al., 2013;Gabor et al., 2018). Amphibians release the hormone corticosterone (CORT; the primary amphibian GC) when exposed to various stressors, including increased pond drying, high salinity, extreme pH and temperature and other water quality variables (Denver, 1998;Chambers et al., 2013;Gomez-Mestre et al., 2013;Narayan and Hero, 2014a, b;Burraco and Gomez-Mestre, 2016;Gabor et al., 2018). In addition, landscape characteristics such as the extent of canopy cover and forest fragmentation around breeding ponds, as well as substrate type, are significant predictors of CORT in adult common toads (Bufo bufo; Janin et al., 2011Janin et al., , 2012 and spotted salamanders (Ambystoma maculatum; Homan et al., 2003), and increased CORT in Jollyville Plateau salamanders (Eurycea tonkawae) has been associated with more urbanized streams . Landscapelevel factors may also influence pond quality, thus affecting aquatic larvae indirectly. Acute changes in pond water salinity can be lethal for larval amphibians (Brown and Walls, 2013). Canopy cover is known to influence amphibian assemblages, with higher amphibian diversity and faster growth in opencanopy wetlands (Werner and Glennemeier, 1999;Skelly et al., 2005). However, increased temperatures associated with climate change combined with open canopy wetlands may reduce wetland hydroperiod, which negatively affects juvenile recruitment and larval survival (Semlitsch et al. 1996;Daszak et al. 2005). These findings suggest a link between habitat quality at multiple spatial scales and physiological health in amphibians.
Microbiota play a role in amphibian immune response to disease (Harris et al., 2009;Woodhams et al., 2016) and may vary with abiotic factors (Bletz et al., 2017a;Jani and Briggs, 2018;Varela et al., 2018), across populations (Kueneman et al., 2013;Hernandez-Gomez et al., 2017;Jani and Briggs, 2018) and among species (Kueneman et al., 2013;Belden et al., 2015;Varela et al., 2018). Microbiota living on the skin combined with skin secretions make up the mucosome, or micro-ecosystem of the skin, providing the first line of defence in amphibian disease resistance (Woodhams et al., 2007;Harris et al., 2009;Woodhams et al., 2014). Identifying the bacterial community of amphibians in a population and the ability of their mucosome to fight infection (mucosome function; Woodhams et al., 2014) can be used to predict disease risk across populations (Woodhams et al., 2014). Elevated CORT over prolonged periods is associated with disease (Warne et al., 2011;Gabor et al., 2013b;Peterson et al., 2013;Gabor et al., 2015Gabor et al., , 2018Kindermann et al., 2017) (Rollins-Smith et al., 2011) and affects gut microbial communities in other taxa (Clarke et al., 2014). Further, the microbiome of amphibians shifts with temperature (Kohl and Yahn, 2016), soil pH, precipitation (Varela et al., 2018) and infection intensity (Jani and Briggs, 2018). Therefore, reduced habitat quality and increased stress from environmental stressors may be associated with altered skin bacterial communities and immunity in amphibians, which could lead to increased declines.

2017), inhibits production of amphibian skin peptides
To explore potential underlying causes of population declines in the ornate chorus frog (Pseudacris ornata), we measured environmental variables at multiple spatial scales, including pond water quality and surrounding landscape characteristics, and examined their effect on three health metrics in tadpoles. Because of the known interactions between habitat quality and population health, we measured and examined the relationships of these variables with CORT release rates, mucosome function, and bacterial community diversity from populations of P. ornata across their range. We tested several hypotheses to assess whether environmental quality affects tadpole health: first, we quantified baseline and acute stress-induced CORT release rates for tadpoles collected across the species' range as a measure of the CORT profiles for a pond. We hypothesized that habitat quality at multiple spatial scales affects CORT release rates and therefore the CORT profile in P. ornata larvae, as environmental conditions are associated with altered physiological health in amphibians (Homan et al., 2003;Janin et al., 2011Janin et al., , 2012Chambers et al., 2013;Gabor et al., 2018). We measured pond water quality [temperature, pH, conductivity and total dissolved solids (TDSs)] as well as land cover characteristics (percent canopy, land cover type and percent developed land) at three spatial scales surrounding each pond. We also tested the hypotheses that these environmental conditions and CORT release rates are associated with altered immune function and an altered skin bacterial community, as both environment and stress are associated with altered immune defences (Rollins-Smith, 2017;Bletz et al., 2017a;Jani and Briggs, 2018;Varela et al., 2018). Our examination of multiple metrics of population health allows us to potentially identify links between environmental stress from reduced habitat quality and population declines which could aid in management practices.

Study species
The ornate chorus frog (P. ornata) is endemic to the southeastern Coastal Plain and longleaf pine (Pinus palustris) ecosystem of the southeastern USA (Means, 2006). This frog has a disjunct range that extends from North Carolina south to Florida and west to Louisiana (Enge et al., 2014;Powell et al., 2016). The ornate chorus frog is a longleaf pine specialist common throughout the panhandle and northern regions of Florida in mesic and xeric habitats, but this species has been declining along its southern range edge in the Florida peninsula and western edge in Louisiana (Bartlett and Bartlett, 2011;Enge et al., 2014). Declines are potentially caused by increased temperatures, droughts and reduced habitat quality associated with fire suppression (Means and Simberloff, 1987;Enge et al., 2014). Fire suppression and insufficient management can reduce environmental quality for species endemic to pyrogenic systems, such as P. ornata and the longleaf pine-wiregrass (P. palustris-Aristida spp.) ecosystem of the southeastern USA (Noss, 2013(Noss, , 2018. Hydroperiod is also considered an important factor associated with P. ornata population persistence, as this species requires fishless, seasonally inundated water bodies with a 3-to 4-month hydroperiod for complete tadpole development (Semlitsch et al., 1996;Enge et al., 2014). Though little is known about the adult frog, it is known to be fossorial, using loose, sandy soils to burrow and has been found over 400 m from breeding ponds (Brown and Means, 1984). Adults are winter breeders, with males calling as early as November, with actual breeding occurring December or January through March, when females deposit eggs on submerged vertical vegetation (Dorcas and Gibbons, 2008;Enge et al., 2014). Roadside ditches, flooded fields and marshes adjacent to forested areas, as well as pine and mixed-woodland ponds can serve as aquatic breeding sites (Dorcas and Gibbons, 2008;Bartlett and Bartlett, 2011;Enge et al., 2014;Powell et al., 2016). However, most breeding ponds are found within sandhills, flatwoods, wetlands and pine forest/plantations and are associated with open-canopy sites with herbaceous understory, characteristic of a shortinterval fire regime (Gorman et al., 2013;Noss, 2013;Enge et al., 2014;Noss, 2018).

Field sampling
We collected tadpoles (stages 30-40; Gosner, 1960) across 7 properties and up to 4 ponds (Sites) per property (Property) for a total of 16 ponds throughout the species' range in Florida, Georgia and South Carolina (Table 1 (Table 1). At each pond, we recorded water temperature, conductivity and TDS; pH was also collected from ponds in 2017 (Table 2). Additionally, we conducted a geographic information systems analysis of land cover at each pond using land cover class and percent canopy cover data from the National Land Cover Database (NLCD, 2011). Using ArcMap (ESRI, Redlands, CA), we created a 100, 500 and 1000 m buffer around each pond. Land cover classes were determined for each 30 m pixel within the three spatial scales. The 16 land cover types recognized by the NLCD were combined into six classes for analysis: Agriculture, Developed, Forest, Shrub, Water and Wetlands. The Agriculture class included grasslands, pasture and cultivated crops; Developed included all intensities of developed and barren land; Forest included deciduous, evergreen and mixed forest cover; Shrub consisted of short, woody plants   herbaceous wetlands. We then quantified the percent of each land cover type and the percent canopy cover within the three spatial scales using the zonal statistics as a table tool in ArcMap. We used the primary land cover class, percent developed land and percent canopy cover for each of the three spatial scales in analyses.

Water-borne CORT release rates
In 2016 and 2017, we used a non-invasive water-borne hormone method (Gabor et al., 2013a(Gabor et al., , 2016 to measure CORT release rates from up to 40 individual tadpoles from each pond (Supplementary Table S1). This method measures CORT secreted through the skin, urine and faeces, which provides an integrated measure of the cumulative effects of chronic stress (Sheriff et al., 2011;Dantzer et al., 2014). We captured tadpoles using dipnets and then placed individuals in 250 ml beakers, one individual per beaker, containing 100 ml of bottled spring water and a perforated Nalgene liner. This liner allowed us to remove the tadpole but leave the water sample in the beaker. Half of the individuals were unmanipulated for 1 hour to obtain baseline CORT release rates and the other half were agitated by shaking the beakers for 1 min every 3 min for 1 hour to measure the response to this acute agitation stressor (Gabor et al., 2016;Forsburg et al., 2019). These baseline and agitation CORT release rates represent the treatment categories. After 1 hour, we removed the liner with the individual from each beaker and poured the water containing leached hormones into HDPE plastic cups. Samples were maintained on ice and transported to the laboratory for immediate extraction or frozen at −20 • C to be processed at a later time (Ellis et al., 2004). Frozen samples were extracted within 1 week of collection. We photographed each tadpole from the side with a ruler for scale before releasing it back into the pond. Snout-vent length (SVL, in mm) was measured from photographs using the program ImageJ (Rosband, 1997). In 2017, all tadpoles from which water-borne hormones were collected were also swabbed for skin bacterial community analysis. After the water samples containing leached hormones were collected, each individual was swabbed ventrally from the mouth to the vent using a  single, sterile swab. The cotton tips of each swab were placed in Eppendorf tubes, frozen and transported back to the lab for analysis.

Hormone extraction and validation
Once defrosted, we filtered each water sample through standard coffee filters (equivalent to grade 4 filter paper) to remove large debris and faecal material. All samples were then drawn through C18 solid phase extraction (SPE) columns (SepPak Vac 3 cc/500 mg; Waters, Inc.) primed with 4 ml methanol and 4 ml distilled water (Gabor et al., 2016). After hormone extraction, SPE columns containing the hormone residue were stored at −20 • C for up to 3 months before elution. We eluted columns with 4 ml methanol into borosilicate test tubes. Eluted samples were dried using nitrogen gas flowing through an Evap-O-Rac (Cole-Palmer) while sitting in a hot water bath (37 • C). Each dried sample was then resuspended in a mixture of 5% ethanol and 95% enzyme-immunoassay (EIA) buffer (Cayman Chemical Company, Inc.) to a final volume of 260 μl. We also ran blank controls of bottled spring water resuspended to a final volume of 130 μl. We reconstituted at the minimum volume required to plate each sample in duplicate (with some leftover) because water samples have such low CORT values. The samples were resuspended to 260 μl to dilute the CORT concentrations to within the range of the standard curve. Each sample was run in duplicate on CORT EIA plates (N = 26; # 501320, Cayman Chemical Company, Inc.) and absorbance was read using a spectrophotometer plate reader (BioTek ELX 800) set to 405 nm. Final CORT values (pg/ml) were multiplied by amount resuspended (0.260 ml) and divided by SVL for final unit of pg/mm per hour. Water samples were multiplied by reconstitution volume (0.130 ml) and the relevant spring water values (5.4-15.6 pg/sample) were subtracted from the CORT release rates of each tadpole (adjusted range: 17.8-1957.4 pg/mm per hour).
The water-borne CORT collection method was validated for P. ornata using a pooled sample of hormones from seven non-experimental tadpoles serially diluted (following Gabor et al., 2016). We examined parallelism of the serial dilution curve (1:1-1:32) of the pooled sample to the known standard curve (comparison of slopes, t 11 = 1.304, P = 0.22). We also assessed quantitative recovery by spiking the pooled sample with each of eight standards. The observed recovery ranged from 62.0% to 82.3%, and we found a linear relationship between observed and expected slopes (β = 0.77, F 1,6 = 1817.75, R 2 = 0.997, P < 0.0001). Using a different pooled control sample run in quadruplicate on each plate, we examined the intra-and inter-plate variation of the 26 plates. Intra-plate variation ranged from 0.29-16.33% and mean inter-plate variation was 13.11%. The sensitivity of the CORT EIA plates ranged from 41.47 to 1390.60 pg/ml on average.

Mucosome function
In 2017, 1 ml water was removed from each baseline CORT sample and stored at −20 • C to assess mucosome function. A subset of these samples was used to assess mucosome function of individuals from each pond (Supplementary Table  S1). Using the Cell Titer Glo 2 kit (Promega), a 25-μl water sample-containing host skin peptides, mucosal antibodies and bacterial communities-was combined with a 25-μl solution of Batrachochytrium dendrobatidis (Bd) zoospores with known concentration (approximately 25 000 total zoospores per 25 μl) on 96-well plates (Barnhardt, 2018). Bd zoospores were collected from plate cultures of isolate JEL423, which is within the Global Panzootic Lineage (BdGPL), by flooding 3-to 5-day-old culture plates with autoclaved water to stimulate zoospore release from zoosporangia. This Bd strain was used because it was readily available, easy to grow in the laboratory and detected across North America (Farrer et al., 2011;Rosenblum et al., 2013;Marshall et al., 2019). Each mucosome sample was then plated in triplicate to assess Bd viability in the presence of the mucosome. In addition, six heat-killed Bd standards (0%, 20%, 40%, 60%, 80%, 100%) were plated in triplicate. After incubation for 1 hour at room temperature, 50 μl of Cell Titer Glo reagent was added to each well and placed on an orbital shaker at 200 rpm for 2 min and incubated at room temperature for an additional 15 min. A luminescent plate reader was then used (Biotek Synergy H1) to assess the percent cell viability from a ratio of live:dead Bd to determine the mucosome function against the pathogen. The gain of the plate reader was set at 150.

Skin bacterial diversity
We extracted DNA from skin swabs of approximately 10 tadpoles per pond (Supplementary Table S1) using 50 μl of PrepMan (Applied Biosystems, Inc.) following the manufacturer's protocol. Swabs and extracts were spun down briefly then the swab was inverted with sterile forceps and spun down again at 2 348 x g for 1 min. The swab was removed from the tube with sterile forceps and the remaining extract was centrifuged at 21 130 x g for 5 mins to pellet any precipitates. Without disturbing the pellet, 40 μl of extract was transferred to a new 1.5-ml centrifuge tube (adapted from Becker et al., 2015). The purified extracts were used to generate PCR amplicons of bacterial 16S V4 properties using 515F and GOLAY barcoded 806R primers (Caporaso et al., 2011) following the Earth Microbiome Protocol (Gilbert et al., 2014). Amplicons were quantified using a Qubit dsDNA Broad sensitivity Assay Kit (Invitrogen) and, after normalizing and pooling, amplicons were size selected on a BluePippen (SageScience, Inc.) using a 2% agarose 100-600 Bp cartridge, the resulting fraction was cleaned using AMPure XP magnetic beads (Beckman Coulter, Inc). The pooled library size (∼350 Bp) and concentration were verified on a TapeStation 2200 (Agilent Technologies, Inc.) using a D1000 screen tape and reagents. Sequencing was performed on a MiSeq instrument (Illumina, Inc.) using the 600-cycle MiSeq Reagent Kit v3. The resulting reads from the MiSeq were trimmed of adapters and parsed in BaseSpace (Illumina, Inc.) according to their respective barcode. The dada2 version 1.5.0 pipeline (Callahan et al., 2016) in R (version 3.5.2; R Core Team, 2018) was used to first inspect read quality profiles and then filter and trim reads. Reverse reads were not used in downstream analyses owing to low quality, and thus, paired read merging was skipped after dereplication. Forward sequence reads were trimmed to 220 Bp based on their quality profile. Taxonomy assignment was performed using dada2 formatted FASTA file from the Silva taxonomic training data (version 132). The resulting dada2 tables were merged with the sample metadata and analysed using phyloseq (McMurdie and Holmes, 2013). Using phyloseq, samples were rarified to 25 000 sequences (average number of sequences = 229 663) and alpha diversity indices (Richness, Shannon and Simpson) were calculated. Richness is the number of species present, whereas Shannon diversity and Simpson diversity are based on the proportion of individuals of a particular species within the sample. The Shannon diversity index provides a value of diversity from estimating how likely an unknown individual is from a known species. Simpson diversity places a higher importance on dominant or more abundant species and is the probability two individuals are from different species. In these indices, a higher value (uncertainty/probability) indicates a higher diversity. Measuring all three indices provides a more complete picture of the bacterial diversity and composition on the amphibian skin. The final data set consisted of 8 337 operational taxonomic units (OTUs) across 98 samples.

Statistical analysis
All analyses were run using R (version 3.5.2; R Core Team 2018). First, to determine if there was spatial autocorrelation in CORT profiles across ponds, distance-based Moran's eigenvector mapping (MEM) was used to map longitude and latitude for each pond and build spatial predictors using the following packages: adespatial (Dray et al., 2018), geosphere (Hijmans, 2017) and vegan (Oksanen et al., 2019). Permutation analysis of variance (ANOVA) on the distance-based Moran axes identified two of three spatial axes as significant predictors of CORT variation. These two spatial axes, MEM2 and MEM3, corresponded to scaled distance measures of longitude and latitude, respectively. These axes were extracted and incorporated into model building. Eigenvalues for these two positive spatial axes were −0.78 and 0.24, respectively, indicating they represent autocorrelation between ponds at a fine spatial scale.
We developed linear mixed effect (lme) models in the package nlme (Pinheiro and Bates, 2018) to examine predictors of (i) CORT release rates, (ii) mucosome function and (iii) skin bacterial diversity indices (Richness, Shannon, Simpson) among ponds as separate response variables. All models for all response variables included Site as a ran- dom effect. CORT models examined natural log-transformed CORT release rates standardized by SVL (pg/mm per hour) as the response variable. Prior to building models, we first ran a preliminary analysis of models containing each land cover characteristic calculated from NLCD variables (landcover type, percent urban development and percent canopy cover) at each spatial scale as sole predictors for each of the response variables. To determine the most important scale for the three predictors, the models containing the landcover characteristic at the different scales were ranked according to Akaike's information criterion adjusted for small sample size (AICc) and the scale included in the highest ranked model was retained for subsequent models. Water quality variables (water temperature, conductivity, TDS, pH) from all ponds were combined using a principle components analysis (PCA) to reduce the number of variables in analyses. We ran analyses to assess predictors of CORT release rates for each year (2016 and 2017) separately, as not all ponds were sampled in both years and pH was an additional pond characteristic added in 2017. For 2016, PCA revealed that PC1 accounted for 83.3% of the variation in the data and was driven by conductivity and TDS, whereas PC2 accounted for 16.7% of the variation and was mainly driven by water temperature (Table 3). In 2017, PC1 accounted for 59.6% of the variation and was driven by conductivity and TDS, whereas PC2 accounted for 26.4% of the variation and was driven by water temperature (Table 3). To test the relationship of predictors to each of the three response variables, models incorporated main effects including the dominant land cover type within 100 m of ponds (Landcover100), the percent developed land within 1000 m of ponds (Urban1000), the percent canopy cover within 500 m of ponds (Canopy500), PC1 and PC2 for each year, spatial axes (MEM2 and MEM3), the property each pond was located within (Property), as well as additive models. We created a total of 28 models for analyses (Supplementary Table S2). Additionally, we analysed each treatment (Baseline or Agitation) separately within each year. A similar list of models was used to examine predictors of mucosome function and skin microbial diversity (Supplementary Tables S3 and S4), substituting the correct scale of the landscape predictors as determined by separate preliminary analyses and adding baseline CORT release rates (BCORT) as a sole predictor in an additional model. Because other analyses were run to assess impacts of both Site and Property on skin bacterial diversity indices (PERMANOVA; see below), Property was not included as a predictor in the bacterial diversity models. This resulted in 29 models for mucosome analysis and 28 models for skin bacterial diversity analysis. We used each bacterial diversity index (Richness, Shannon and Simpson) as a separate response variable in separate analyses involving the entire model set. All models were ranked according to Akaike information criterion corrected for small sample size (AIC c ) using the package MuMIn (Barton 2018). We calculated parameters using the maximum-likelihood estimation during the model-selection process. Model-averaged parameter estimates, unconditional standard error (SE) and unconditional 95% CIs were calcu- lated for candidate models ( AICc < 2) using the package AICcmodavg (Mazerolle 2019). In addition, we ran independent t-tests to compare baseline and agitation-induced CORT release rates for each pond and year separately to specifically address which populations might be chronically stressed, indicated by an inability to mount a CORT response above baseline levels. We also ran an ANOVA and Levene's test using the package lawstat (Gastwirth et al., 2017) to examine differences in mean mucosome function and mean variance across Sites.
Means of each of the three alpha diversity indices were compared across Site and Property. We analysed mean Shannon diversity across Sites and Properties using an ANOVA. Because data were non-normal, we used a Kruskal-Wallis test to compare the mean Richness and Simpson diversity across Sites and Properties. The relative contribution of Site and Property to microbiome diversity (beta diversity) were analysed using PERMANOVA (Adonis; 999 permutations) using the package vegan (Oksanen et al., 2019). All pairwise comparisons were assessed for significant factors using the package pairwiseAdonis (Martinez Arbizu, 2019). Beta diversity of skin bacterial communities (OTUs) were plotted using the Bray-Curtis method of non-metric multidimensional scaling (NMDS) and the following packages: phyloseq (McMurdie and Holmes, 2013) and ggplot2 (Wickham, 2016).

Water-borne CORT release rates
AICc model selection indicated six candidate models ( AICc < 2) to predict baseline CORT release rates and five candidate models to predict agitation CORT release rates for sites sampled in 2016 (Supplementary Table S5). Model averaged parameter estimates indicated both PCs and both spatial axes were top predictors of baseline CORT release rates, though no predictor was significant (CI did not overlap 0; Table 4). Model averaged parameter estimates indicated both PCs and MEM3 (scaled latitude) were significant predictors of agitation CORT release rates in 2016 (Table 4)  Only parameters from the top-ranking models ( AICc < 2) are included. Bolded terms indicate parameters used for inference whose 95% CIs do not overlap 0. Fig. 2a) and lower in ponds with lower water temperature (negative loading on PC2; Fig. 2b). Additionally, agitation CORT release rates were higher at lower MEM3 values (Fig. 3). For 2017, two models predicting baseline CORT release rates and two models predicting agitation CORT release rates had a AICc < 2 (Supplementary Table S6). Model averaged parameter estimates indicated the percent developed land within 1000 m of ponds was a significant predictor of both baseline and agitation CORT release rates for ponds sampled in 2017 (Table 5). Baseline and agitation CORT release rates were both higher in ponds with more developed land within 1000 m (Fig. 4). In 2016, the marginal R 2 for the top model explaining baseline and agitation CORT release rates were 0.24 and 0.26, respectively. The conditional R 2 for these same models were 0.34 and 0.32, respectively. The marginal R 2 is calculated using only the fixed effects, whereas the conditional R 2 includes both fixed and random effects. Therefore, the inclusion of random effects and fixed effects explained more variation in CORT release rates than fixed effects alone. The difference between marginal and conditional R 2 values was greater for 2017, where the marginal R 2 for the top model explaining baseline CORT release rates was 0.23, and for agitation was 0.21, with the conditional R 2 values of 0.38 and 0.60, respectively. Average CORT release rates across both years were highest at the site in Lafayette Forest Wildlife Environmental Area (LF1; Fig. 5). Examining CORT release rates within each treatment for each pond indicated 40% (N = 4) of ponds sampled in 2016 and 50% (N = 6) of ponds sampled in 2017 did not show significant differences between baseline and agitation CORT release rates, though only one population (JC3) did not show a difference across both years (Supplementary Table  S7; Fig. 5). However, not all ponds were resampled both years owing to a lack of tadpoles or complete loss of the pond in one case.  Only parameters from the top-ranking models ( AICc < 2) are included. Bolded terms indicate parameters used for inference whose 95% CIs do not overlap 0.

Mucosome function
There was no significant difference in the mean mucosome function across all ponds (ANOVA: F 11,109 = 0.476, P = 0.914; Supplementary Fig. S1). Levene's test further determined no difference in the mean variance across all ponds (W = 1.19; P = 0.302). Mean Bd viability in the presence of the mucosome was 88.1%. Five candidate models ( AICc < 2) were selected from AICc model selection to explain the variation in mucosome function (Supplementary Table S6). These models included PC1, PC2, percent developed land within 100 m (Urban100) and the percent canopy within 500 m (Canopy500) as top predictors, though model averaged parameter estimates did not indicate any significant predictor (Table 5).

Skin bacterial diversity
Alpha diversity of skin bacterial communities differed across Sites and Properties in both Richness and Shannon diversity (Richness, Site: χ 2 = 42.80, df = 9, P < 0.001; Property: χ 2 = 35.60, df = 5, P < 0.001; Shannon, Site: F 9,88 = 2.32, P = 0.022; Property: F 5,92 = 3.18, P = 0.011), but not in Simpson diversity (Simpson, Site: χ 2 = 13.88, df = 9, P = 0.127; Property: χ 2 = 7.61, df = 5, P = 0.179). Measures of diversity tended to be highest at sites within the Jones Center and lowest at sites within Eglin AFB, with the exception of Richness that was lowest within Apalachicola National Forest (Supplementary Table S8). Two candidate models ( AICc < 2) were selected to explain the variation in Richness across sites (Supplementary Table S6). Model averaged parameter estimates indicated dominant land cover within 1000 m of each pond was a significant predictor of bacterial richness (Table 5; Fig. 6), with Richness being highest in forest and lower in large wetlands and in ponds surrounded by shrubby vegetation. The top model explaining Richness had a marginal R 2 of 0.24 and a conditional R 2 of 0.36, indicating some of the variation was attributed to the random effect of Site. Six candidate models were selected to explain Shannon diversity (Supplementary Table S6). Model averaged parameter estimates indicated that the percent developed land within 500 m (Urban500; Supplementary Fig. S2) and the percent canopy cover within 100 m (Canopy100; Fig. 7a) were significant predictors of Shannon diversity (Table 5; .  , with diversity being higher in ponds surrounded by more developed land and lower in ponds with more canopy cover. Additionally, Shannon diversity was higher as MEM2 increased (Table 5; Fig. 7b). The marginal and conditional R 2 for the top model explaining Shannon diversity were 0.07 and 0.10, respectively. Two candidate models were selected to explain Simpson diversity across ponds (Supplementary Table  S6), with MEM2 (scaled longitude) being the only significant predictor (Table 5; Fig. 8). As MEM2 increased, Simpson diversity also increased. The marginal and conditional R 2 for the top model explaining Simpson diversity was 0.14. Prominent bacterial families included Burkholderiaceae, Sphingomonadaceae, Xanthomonadaceae, Pseudomonadaceae, Enterobacteriaceae and Desulfovibrionaceae (Fig. 9). The relative contribution of both Site and Property accounted for significant variation in microbiome beta diversity when analysed with a permutational ANOVA (Site: ADONIS R 2 = 0.47, P = 0.001; Property: ADONIS R 2 = 0.32, P = 0.001). All pairwise comparisons between Sites were significant (P < 0.05) except for the comparison of site AP1 to AP2 (F = 1.46, P = 0.144) and AP1 to SM4 (F = 1.58, P = 0.117; Supplementary Table S9). All pairwise comparisons between Properties were significant (P < 0.001) except for the comparison of the Apalachicola Property to St. Marks (F = 1.80, P = 0.066; Supplementary Table S10). An ordination plot based on the Bray-Curtis method of NMDS showed bacterial community similarities within Sites and across Properties (Fig. 10).

Discussion
We found that both water quality and land cover characteristics were associated with increased CORT release rates and altered skin bacterial communities among populations of P. ornata. We found both baseline and agitation CORT profiles were linked to pond water quality, as well as the amount of developed land within 1000 m of ponds. These results indicate that environmental variation in two different settings-the aquatic and surrounding terrestrial habitatsoperate to influence the physiology (GCs) of larval P. ornata. Similarly, pond water quality, land cover and spatial location of ponds were significant predictors of mucosome function and skin bacterial diversity, indicating environmental variation may affect larval P. ornata immune defences through changes in bacterial communities, though we did not see significant effects of these predictors directly on mucosome function. Therefore, perturbations and changes to environmental conditions at multiple spatial scales may significantly impact population health in this species. Given that some of these populations differed from year to year in their CORT profile, the impact of environmental variation-along with the ability to recover from stressor effects-likely varied over time, showing a need for continued monitoring (Blaustein et al., 2011). Taken together, the effects of environmental variation on GC function and bacterial communities across years suggests that persistent deviations from optimum conditions may affect population health, thus setting the stage for population declines.

Water-borne CORT release rates
Water quality variables and nearby urban development were important predictors of larval ornate chorus frog CORT profiles. In 2016, both baseline and agitation CORT release rates were higher in ponds with higher water temperature, conductivity and TDS, as shown by the relationship with PC1 and PC2, whereas the percent urban development within 1000 m was positively associated with CORT release rates in 2017. Additionally, agitation CORT release rates in 2016 were higher at lower MEM3 values, indicating ponds closest to each other were most similar and ponds at lower latitudes had higher CORT release rates, probably from the warmer water temperatures at the southernmost site (LF1). Similarly, studies involving salamanders also found elevated waterborne CORT release rates associated with higher temperatures (Novarro et al., 2018;Millikin et al., 2019) and higher conductivity (Chambers, 2011); although, conductivity did not significantly affect CORT in Rana sylvatica or Hyla versicolor tadpoles (Chambers, 2011). Both acute and chronic thermal stressors are associated with increased urinary CORT in Rhinella marina, the invasive cane toad (Narayan et al., 2012;Narayan and Hero, 2014a and b). Elevated CORT can assist with energy mobilization and metabolic increases associated with higher temperatures but may also indicate stress (Chambers et al., 2013;Narayan and Hero, 2014a and b;Novarro et al., 2018;Millikin et al., 2019). Several studies show increased developmental rate and body condition of tadpoles in warmer water (Schiesari, 2006;Reading, 2010), though stress is generally associated with lower body condition and reduced growth in amphibians (Glennemeier  Ponds are ordered geographically from west to east. and Denver, 2002;Hu et al., 2008;Denver, 2009;Janin et al., 2011;Crespi and Warne, 2013). Therefore, higher water temperatures may be an environmental stressor for this winter-breeding amphibian, as warmer water temperatures may result in a shorter hydroperiod, which is a known environmental stressor itself (Denver, 1998;Gomez-Mestre et al., 2013). Our results suggest environmental variation in water quality characteristics are important determinants of the CORT profiles of larval P. ornata, potentially affecting their physiological health. Pond breeding frogs may be affected by habitat alterations at multiple scales due to their use of both ponds, as larvae and to breed, and terrestrial habitats as adults. Janin et al. (2011) found that body condition was lower and urinary CORT metabolites were higher in adult Bufo from ponds characterized by low forest availability and high fragmentation within 500 m, suggesting surrounding habitat affects amphibian health. We found that ponds with higher percentages of developed land within 1000 m were associated with higher water-borne CORT release rates in 2017. Gabor et al. (2018) found Jollyville Plateau salamanders (Eurycea tonkawae) had higher water-borne CORT release rates in more urbanized streams as indicated by greater percent impervious cover within the watershed. We also found that CORT release rates were highest in tadpoles from a pond within Lafayette Forest Wildlife Environmental Area (LF1), which had some of the warmest water temperatures and was dominated by shrubby vegetation, a result of infrequent prescribed fires. Of the properties sampled, this land was the most recently acquired and managed by the state of Florida (purchased in 2008, management began in 2010) with the first prescribed fires conducted in 2013 (N. Lambert, Florida Fish and Wildlife Conservation Commission, pers. comm.). Prescribed fire is vital to management of longleaf pine uplands and their embedded wetlands because it controls the growth of invasive, shrubby, mid-story vegetation such as saw palmetto (Serenoa repens) and gallberry (Ilex coriacea; Means, 2006;Noss, 2013;Enge et al., 2014;Noss, 2018). These results suggest habitat quality at a larger scale may affect the CORT profile and subsequent physiological health of P. ornata tadpoles, whether by directly affecting larval health, or indirectly via increased stress on adults.
Elevated CORT in response to an acute stressor is indicative of a healthy endocrine response. Tadpoles that elevate CORT release in response to an acute agitation stressor indicate the HPI axis is still active and not dysregulated due to chronic stress. Therefore, the ability to mount a response to an acute stressor provides one component that can indicate a healthy population. Comparing the average baseline and agitation CORT release rates from each pond over both years showed four ponds in 2016 and six ponds in 2017 without a significant difference in their stress response (Fig. 5). The lack of a significant stress response in these populations suggests that these ponds may be suffering from chronic stress (Romero, 2004;Cyr and Romero, 2009). These populations may be able to recover if environmental suitability increases and/or if individuals show compensatory growth once they leave the pond and possibly find suitable habitat (Metcalfe and Monaghan, 2001). When individuals return as adults to the pond to lay eggs, their offspring may not show health effects if the water and environmental quality have improved. However, stressful conditions during larval development may have negative post-metamorphic carry-over effects on juvenile and adult frogs (Crespi and Warne, 2013). The populations that did not show a healthy CORT response may fail to recover when exposed to consecutive years of low environmental quality, or if the stressful conditions experienced during the larval stage have carry-over effects in the adults. Populations of P. ornata in the peninsular region of Florida are declining and undergoing local extirpations (Enge et al., 2014). We were unable to find P. ornata tadpoles (nor measure CORT) in this region, and it is possible that these declines may be the consequence of an inability to recover from past harsh environmental conditions (or repeated harsh conditions) such as the drought (and associated warmer temperatures) experienced by this region from 2006 to 2008 and from 2010 to 2012 (Enge et al., 2014). Even though our results indicate environmental variation at multiple spatial and environmental scales (i.e. water quality vs. terrestrial habitat) impact the CORT profiles of larval P. ornata, random effects in the models (as indicated by the conditional R 2 ) accounted for some of the variation in CORT release rates. This suggests additional variability between ponds within each property associated with unmeasured factors. The difference between the marginal and conditional R 2 was greater in 2017, indicating more variation in the CORT release rates between these sites was not accounted for by the predictors. This is likely due to natural differences in pond characteristics (i.e. size, depth, hydroperiod) and other habitat factors not examined. This variation contributes to the uncertainty that enshrouds management decisions for P. ornata and ecologically similar species. However, the consequences of development and water quality on the health of P. ornata, indicated in this study, lends additional support for the conservation and appropriate management of longleaf pine forests, including frequent prescribed fire during the growing season, to maintain high-quality pine savannahs for this and other longleaf pine endemic species (Means, 2006;Noss, 2013;Enge et al., 2014;Noss, 2018).

Mucosome function
The amphibian mucosome, including both skin secretions and microbial communities, are important factors in amphibian health and immunity. The ability of the mucosome to fight infection has been used as a holistic measure of amphibian skin defences (Woodhams et al., 2014). We found no differences in the mean mucosome function across ponds, though the Bd cell viability within a pond was quite variable. Additionally, there were no significant predictors of mucosome function in our models, though water quality and landscape variables were included in the top models. These results indicate no difference in the bacterial and mucosal defences across ponds. In this study, mucosome function only reduced Bd cell viability to 88.1%, indicating this species may be susceptible to Bd infection. Indeed, Gervasi et al. (2017) . found that P. ornata had the highest average Bd infection load of 20 species studied and a high rate of mortality following infection. Similarly, Horner et al. (2017) found P. ornata had a high infection intensity and a high Bd prevalence across ponds. However, due to logistics, we were unable to use the strain of Bd to which these frogs are naturally exposed, which may be why we did not see any difference in mucosome function across ponds. Woodhams et al. (2014) found a significant correlation between mucosome function and both field and laboratory Bd infection prevalence for Swiss and European frogs (for a Swiss Bd strain); thus, future studies could repeat exposure experiments using local isolates.

Skin bacterial diversity
Alpha and beta skin bacterial diversity varied across both Sites and Properties. Primary predictors of bacterial diversity at each pond were the dominant land cover type within 1000 m, the percent developed land within 500 m, the percent canopy cover within 100 m and the spatial location of ponds. All alpha diversity indices were highest at sites located within the Jones Center property, and these wetlands were the largest sampled and surrounded by mature longleaf pine forest with short-interval prescribed fires to maintain high-quality pine savannah. Even though a site in Apalachicola National Forest (AP2) had the lowest observed Richness, both Shannon and Simpson diversity indices were lowest within the Eglin Air Force Base (EG) property. The ponds sampled at Eglin were along a powerline right-of-way and thus were composed of more disturbed habitat. Variation in skin bacterial communities are generally associated with pond characteristics including conductivity, pH, temperature, dissolved oxygen and precipitation (Kueneman et al., 2013;Krynak et al., 2015Krynak et al., , 2016Varela et al., 2018). The significant association of bacterial diversity with land cover suggests landscape-scale features also affect bacterial communities. Bacterial richness, the number of species present, was positively associated with ponds surrounded by forest. Interestingly, Shannon diversity, a measure of the abundance and variety of bacteria present, was higher in ponds with more developed land within 500 m and with less canopy cover within 100 m. Enge et al. (2014) found that most wetlands inhabited by P. ornata lacked a canopy or had an open canopy that allowed growth of grasses during dry down periods, a result of short-interval fire management. Our results suggest these open canopy ponds with greater environmental variation may have higher bacterial diversity.
In southeastern USA, land-use conversion, fire suppression and droughts are important factors in the decline of amphibian species, especially those reliant on the savannah-like longleaf pine-wiregrass (P. palustris-Aristida spp.) ecosystem. In addition to its control of invasive shrubs and hardwoods, summer growing-season fire, when ponds are typically dry, is required to maintain grasses and sedges in ephemeral wetlands used by those pond-breeding amphibians that are endemic to this ecosystem (Noss, 2013(Noss, , 2018. Within the pine savannahs of the southeastern Coastal Plain, a short-interval fire regime maintains herbaceous vegetation, reduces canopy cover in breeding ponds and increases amphibian abundance and diversity (Noss, 2013(Noss, , 2018. Our findings that landcover type, the percent developed land and increased canopy cover were associated with altered skin bacterial diversity indicate environmental variation may also contribute to changes in skin bacterial diversity of larval P. ornata. These effects may ultimately impact the immune defences of species like P. ornata that are specialists of the longleaf pine ecosystem. The sole significant predictor of Simpson diversity, and an additional significant predictor of Shannon diversity, was the spatial MEM axis MEM2 (scaled longitude). The eigenvalue for this axis was −0.78, indicating autocorrelation among sites at a fine spatial scale. This autocorrelation indicates both Shannon and Simpson diversities are most similar in nearby ponds, which is expected, and that longitude had an effect on the diversity observed. Our results suggest ponds sampled in the western portion of the range had lower Shannon and Simpson bacterial diversity than ponds sampled in the eastern portion of the range. The more similar diversity among ponds in close proximity to one another is also suggested in the results of our distance-based ordination plot of bacterial diversity (Fig. 10). However, the pairwise comparison of beta diversity between sites showed that all comparisons between sites except two showed significant differences (Table S8) and all comparisons between properties were significantly different except for those that were geographically the closest (Apalachicola National Forest and St. Marks National Wildlife Refuge; Table S9). As we have seen, environmental characteristics can influence bacterial diversity. Additionally, maintaining high-quality ponds and wetlands across a species range can aid in establishing a 14 higher bacterial diversity. Together, these results suggest site location and landscape characteristics around ponds play an important role in the skin bacterial community composition, which can affect immune defences, ultimately influencing the health of these larval amphibians.
Though bacterial communities are known to vary across sites and even among species within the same site (see Belden et al., 2015), some bacterial families have high proportions of Bd inhibitory isolates. Bletz et al. (2017b) examined skin bacteria from Madagascar amphibians and found three families with 75% or more of the isolates having inhibitory effects on Bd infection, suggesting these may be good candidates for broad probiotic treatment. These three families (Enterobacteriaceae, Pseudomonadaceae and Xanthamonadaceae) were also well represented in our site diversity; although, P. ornata have high Bd prevalence (Gervasi et al., 2017) and the mucosome function in our study only reduced Bd cell viability to ∼88% (see Mucosome function under Results above), suggesting other factors may be contributing to the infection of this species. Continued research should consider the skin bacterial community to aid in developing probiotics to fight known pathogens such as B. dendrobatidis and B. salamandrivorans (Woodhams et al., 2016).

Conclusions
Amphibian populations have been declining rapidly worldwide, and this decline is likely tied to their need for multiple habitats and their susceptibility to environmental changes (Stuart et al., 2004;Becker et al., 2007). There is a growing need to monitor population health and factors that influence their health to provide early warning signs for proactive management decisions. Our results suggest the conservation of P. ornata could be enhanced by management actions that address water quality and forest composition to maintain high-quality aquatic and terrestrial habitat. Natural resource managers could consider habitat quality at multiple spatial scales to best manage larval amphibian populations and the bacterial communities present. By examining baseline and stress (agitation) induced CORT profiles and other health metrics across locations and for several years in succession, researchers and managers may be able to use this method proactively to identify populations at increased risk of declines to focus management efforts. Additionally, little is known about the adult life history of this cryptic and fossorial species (Brown and Means, 1984;Bartlett and Bartlett, 2011), warranting additional research and implementing management practices that maintain suitable habitat for both aquatic and terrestrial life stages, as landscape-level factors may have impacts in aquatic larvae directly, as well as indirectly by affecting adults. By examining environmental variables, the CORT profiles and skin bacterial communities of populations of P. ornata across the range of the species, we identified factors at multiple spatial scales that affect population health. Though some populations showed an ability to mount a CORT response in different years, without habitat remediation, other populations may no longer be able to recover from changes in environmental conditions.

Funding
This work was supported by a cooperative agreement between the U.S. Geological Survey and Texas State University (G15AC00457 to C.R.G. and S.C.W.) and start-up funds from Texas State University to D.R.