Habitat Disturbance Linked with Host Microbiome Dispersion and Bd Dynamics in Temperate Amphibians

Anthropogenic habitat disturbances can dramatically alter ecological community interactions, including host–pathogen dynamics. Recent work has highlighted the potential for habitat disturbances to alter host-associated microbial communities, but the associations between anthropogenic disturbance, host microbiomes, and pathogens are unresolved. Amphibian skin microbial communities are particularly responsive to factors like temperature, physiochemistry, pathogen infection, and environmental microbial reservoirs. Through a field survey on wild populations of Acris crepitans (Hylidae) and Lithobates catesbeianus (Ranidae), we assessed the effects of habitat disturbance and connectivity on environmental bacterial reservoirs, Batrachochytrium dendrobatidis (Bd) infection, and skin microbiome composition. We found higher measures of microbiome dispersion (a measure of community variability) in A. crepitans from more disturbed ponds, supporting the hypothesis that disturbance increases stochasticity in biological communities. We also found that habitat disturbance limited microbiome similarity between locations for both species, suggesting greater isolation of bacterial assemblages in more disturbed areas. Higher disturbance was associated with lower Bd prevalence for A. crepitans, which could signify suboptimal microclimates for Bd in disturbed habitats. Combined, our findings show that reduced microbiome stability stemming from habitat disturbance could compromise population health, even in the absence of pathogenic infection.


Introduction
Anthropogenic habitat alterations such as deforestation can lead to instability and decline of animal populations by disrupting ecological interactions [1][2][3][4][5][6][7]. For example, habitat disturbance can influence interactions between predators and prey [8], plants and pollinators [9], and hosts and pathogens [10,11]. The effects of habitat disturbance on disease dynamics are often nuanced and the driving mechanisms may vary among landscapes and pathogen transmission types [12,13]. Some of these factors include increased population density which could concentrate susceptible hosts, reduced ecosystem complexity which could increase encounter rate between susceptible hosts and pathogens, and increased host stress which could reduce immune function [14][15][16]. Additionally, habitat alterations can lead to changes in terrestrial and aquatic microbial communities, with downstream effects on host-associated microbiomes, defined as microbial communities and their metabolic products [17][18][19]. Temperature variation, habitat connectivity, and aquatic soil influx, which are key variables affected by anthropogenic habitat disturbance, have all been correlated with host microbiome composition [20][21][22]. Consequential changes to host microbiome composition following habitat disturbances can have sub-lethal negative impacts on host immunity, potentially increasing disease susceptibility and threatening host health [23].
Host microbiomes play vital roles in organism functions like metabolism and immunity [24][25][26]. While many studies have found that higher microbiome diversity confers beneficial effects to the host through increases in host immune function and direct defense against invading pathogens [27][28][29][30], more diverse microbiomes are also more likely to include pathogenic taxa, or may indicate a disrupted state [31][32][33]. Microbiomes can enter a state of dysbiosis when beneficial microbes are depleted, pathogenic microbes are enriched, or overall diversity is reduced [27,31,34,35], but classifying microbiome dysbiosis can be difficult as it varies between hosts and systems. In some systems, higher variance in the microbiome of individuals can be used to identify microbiomes in a state of dysbiosis [31]. Variance can be represented by greater dispersion, which is the average distance from data points to the group centroid in the consolidated microbial community data represented by principal coordinates. This has been dubbed the "Anna Karenina principle," where microbial communities become more variable between hosts under stress, leading to greater spread along principal component axes [31]. This principle is supported by recent studies on corals [18], anemones [36], frogs [21,37], and humans [38,39].
Relatively few studies have examined interactions between habitat disturbance, host microbiome composition, and invading pathogens. Amphibians often require multiple habitat types to complete their life cycle. They are thus particularly responsive to anthropogenic habitat disturbances that may alter biotic and abiotic environmental conditions and influence their movement. A recent study found that proximity to anthropogenic development was associated with higher dispersion and altered function of the microbiome of a neotropical frog species [37]. Recent work has also revealed associations between the composition of the amphibian skin microbiome and the outcome of infection with Batrachochytrium dendrobatidis (Bd), a waterborne fungal pathogen that is transmitted between hosts sharing the same habitat [26,30,[40][41][42][43][44]. In line with work on the microbiome of other organisms, classifying protective aspects of the amphibian microbiome is complex and often conflicting between host species and study systems [45][46][47]. The outcome of pathogen infection may, therefore, be determined by the function of key members in the skin microbiome, and the composition of the microbiome as a whole. Disturbed habitats likely impose multiple correlated and contrasting challenges on amphibians, their skin microbiomes, and disease resistance, so studies combining multiscale measures of disturbance are needed.
In this study, we use a Bd-enzootic system in southeastern USA to determine the influences of habitat disturbance on microbiome-Bd dynamics. In line with recent studies [21,31,[36][37][38][39], we hypothesized that habitat disturbance, connectivity, and associated pond water physiochemistry would drive bacterial assembly in the environment, which together would shape the composition of the amphibian skin microbiome. We proposed two contrasting predictions for how habitat disturbance would impact Bd: (1) disturbance could directly reduce Bd prevalence and infection loads by shifting optimal microclimates for Bd [48,49] or (2) disturbance could increase stochasticity in the microbiome, which could reduce host-protective effects of the microbial community and therefore indirectly increase Bd prevalence and loads. Through this study, we shed light on the links connecting landscape change to microbiome disruption and disease in temperate amphibian populations.

Sample Sites and Study Design
From 22 April to 8 May 2019, we sampled 30 permanent ponds in Tuscaloosa County, AL (Fig. S1), collecting data on frog skin microbial diversity, Bd infection, and environmental parameters. We captured sufficient sample sizes of frog species at 19 of the sites, but we used all 30 ponds for environmental analyses (Table S1). For each sampling site, we calculated the percent of anthropogenic and natural landcover within a 500-m radius of each pond using Arc-GIS (version 10.5.1). We used GAP Land Cover Data [50] to classify land cover as either natural (aquatic vegetation, forest and woodland, introduced and semi natural vegetation, nonvascular and sparse vascular rock vegetation, or shrubland and grassland), disturbed (agricultural vegetation, developed and other human use, or recently disturbed or modified), or open water at 30-m resolution. For statistical analysis, we used continuous (percentage) and binary (high > 45% or low < 45% as in [20]) measures of habitat disturbance. Of the 19 sites, 10 had relatively low levels of disturbance and 9 had high levels of disturbance.
Based on initial surveys, we focused on 2 common frog species: the northern cricket frog (Acris crepitans) and the American bullfrog (Lithobates catesbeianus). These two species are both native to Alabama and while occupying similar habitats are from distinct families (Hylidae vs. Ranidae, respectively). This provides the opportunity to study the responses of diverse taxa to disturbance. At each site, we sampled up to 20 frogs of the most abundant species, resulting in 156 A. crepitans from 10 sites (6 high disturbance, 4 low disturbance) and 100 L. catesbeianus from 9 sites (3 high disturbance, 6 low disturbance), with no overlapping sites between frog species. We captured an additional 7 A. crepitans at a local lake that we use in our habitat connectivity analysis but excluded from other analyses due to the differences in lake water parameters.
We conducted nocturnal surveys by searching along pond edges and captured frogs using dipnets or by hand. Sampling was random with regard to sex and size. All A. crepitans were similar in size, while L. catesbeianus varied (see Supplemental Data). Dipnets were rinsed between individuals and new plastic bags were placed over the searcher's hand.
Field equipment (boots, waders, and nets) were sterilized with bleach between sample sites to prevent the spread of microorganisms. We temporarily held frogs in plastic bags until processing. Each individual was weighed, measured (snout-vent length), and swabbed, then released at the location of capture. We rinsed amphibians with sterile distilled water prior to swabbing to reduce sampling of unassociated environmental microbes. Individuals were swabbed using a sterile swab (Medical Wire) 5 times on the underside of each foot and 5 times on each side of the venter for a total of 30 passes [51]. We dry-stored swabs in sterile tubes on ice while in the field and moved them to a − 20 °C freezer upon return to the lab, where they remained until DNA extraction.
During each sampling trip, we also collected a range of water quality parameters from surface water at 3 locations in each pond (Table S1). Using a YSI Pro Plus multiparameter meter (Yellow Springs Inc., Yellow Springs, OH), we measured temperature, pH, dissolved oxygen (DO), and conductivity. We measured turbidity using the Aquafluor handheld fluorometer/turbidimeter (Turner Designs), which was rinsed between each measurement. From the same 3 locations where we measured water parameters, we also collected 500 mL of surface water to assess nutrient levels and characterize environmental bacteria communities. We collected water in sterile, acid-washed, and autoclaved polyethylene bottles (Nalgene) and flushed the bottles 3 times in pond water prior to sample collection. Water samples were stored on ice until returning to the lab where they were refrigerated and filtered within 24 h. In the lab, we homogenized each water sample separately through gentle inversion, then vacuum filtered 200 mL through a 0.45-µm mixed cellulose ester filter using acid-washed Erlenmeyer flasks and magnetic filter funnels (Pall). Filters were cut into four wedges using sterile scissors and stored in sterile microcentrifuge tubes at − 20 °C until DNA extraction. A 0.45-µm filter size was chosen due to water turbidity levels, which led to clogging of 0.22-µm filters after only 50 mL of water filtration. We compared 3 water samples filtered through 0.22 µm to the same 3 samples filtered through 0.45 µm and found more OTUs recovered by filtering more water through the larger filter size (140.3 ± 22 vs 150 ± 24, respectively). Additionally, 69% of OTUs were recovered by both filters while 18% were only recovered by the 0.45-µm filter and 13% were only recovered by the 0.22-µm filter. Following filtration, a sample of the filtered water was stored in a 15-mL falcon tube at − 20 °C for later nutrient analysis. Nitrate, nitrite, orthophosphate, and ammonia levels were measured using a SEAL digestion block (SEAL Analytical Inc., Mequon, WI).

Molecular Methods
We extracted DNA from skin swabs and water filters using the Qiagen DNeasy Blood & Tissue Kit following a modification of the standard protocol. After the addition of proteinase K, we incubated at 56 °C for 12 h instead of 4 h to increase DNA yield [52], and we doubled all vortex times for water filters to increase the removal of bacteria from filter fibers. We included negative controls of sterile water to monitor potential contamination of extraction reagents.
For qPCR analysis of Bd from skin swab DNA, we diluted DNA 1:10 to reduce reaction inhibition and quantified Bd loads using Taqman qPCR assays on Bd-specific ITS and 5.8S genes [53] and used gBlock synthetic Bd standards (Integrated DNA Technologies) diluted from 10 6 to 10 2 gene copies (g.c.). We ran plates in duplicate, with mismatching samples run in triplicate. Only samples that were positive on 2 plates were recorded as positive. Bd loads were averaged across the duplicate plates and log10-transformed to correct for non-normal residual distributions. Because our two host species vary substantially in average body mass, we divided Bd loads by individual body mass to reduce the influence of swab area on Bd loads.
For metabarcoding bacterial communities from skin swabs and water samples, we followed the Earth Microbiome Project 16S Illumina Amplicon Protocol [54,55], which targets the V4 region of the bacterial 16S rRNA gene using a dual-index approach with 515F and 806R barcoded primers. We PCR-amplified extracted bacterial DNA in duplicate using the following recipe per sample: 12.2 µL of UltraPure water, 4 µL of 5X Phire reaction buffer (Thermo Scientific), 0.4 µL of 10-mM dNTPs (Invitrogen), 0.4 µL of Phire Hot Start II DNA Polymerase (Thermo Scientific), 0.5 µL each of 10-µM barcoded forward and reverse primers (Integrated DNA Technologies), and 2 µL of sample DNA. We ran duplicate PCR plates on SimpliAmp thermal cyclers (Thermo Scientific) according to the following protocol: 98 °C for 3 min, 38 cycles of 98 °C for 5 s, 50 °C for 5 s, and 72 °C for 15 s, then 72 °C for 3 min before holding at 12 °C. After optimization trials, we diluted DNA from water filters 1:100 prior to PCR amplification to reduce inhibition. PCR inhibition was not observed with skin swabs, negating the need for dilution. We included negative controls (water without template DNA) to monitor any potential contamination of PCR reagents. We combined duplicate plates and visualized amplicons in 1% agarose gel to confirm DNA amplification, then pooled 2 µL of each sample into two amplicon libraries, with samples randomly allocated to each library. We purified the libraries using the QIAquick Gel Extraction Kit (Qiagen), then measured DNA concentration using the Qubit 2.0 fluorometer with a dsDNA Broad-Range Assay Kit (Invitrogen). Concentrations of the purified libraries were 42.5 and 21.8 ng/µL. The 16S libraries were sequenced using Illumina MiSeq at Tufts University Core Facility (TUCF Genomics), Boston, MA, USA. All bacterial sequences were deposited in the NCBI Sequence Read Archive (BioProject PRJNA761342).
We received demultiplexed bacterial sequences, then imported forward reads for each sample into Quantitative Insights into Microbial Ecology II (QIIME2 version 2019.10). We used QIIME2 to generate OTU tables and extract metrics of alpha and beta diversity for frog and water bacterial microbiomes. Prior to analyzing sequence data, we trimmed sequences to 150 bp, filtered by quality score [56], and then used the deblur pipeline to cluster sequences into operational taxonomic units (OTUs) based on 97% similarity. Taxonomy was assigned to OTUs using the Greengenes 13.8 reference sequence database. We removed chloroplast and mitochondrial sequences and removed reads from any OTUs contributing less than 0.005% of total reads [56], then rarefied the OTU table to 4000 reads based on rarefaction curves (Fig. S2). This resulted in 99 of 414 samples being excluded, including all extraction and PCR negative controls. For frog and water samples, we used the Shannon diversity index and the number of OTUs as metrics of alpha diversity and Bray-Curtis dissimilarity matrices as a metric of beta diversity. We used principal coordinates (PCo) to visualize patterns in beta diversity. The number of OTUs was highly correlated with Faith's phylogenetic diversity index (R 2 = 0.95, P < 0.0001) and moderately correlated with unweighted UniFrac PCo1 (R 2 = 0.39, P < 0.0001), so patterns in the number of OTUs represent patterns in all three metrics.

Statistical Analyses
We calculated microbiome dispersion (distance from centroid) for each sample using the betadisper function in Program R (vegan package; [57]). We calculated dispersion based on Bray-Curtis dissimilarity matrices grouped by host species. Higher values indicate samples that are further from the average community composition. Sampling was random by sex and body mass. Thus, this metric provides an indication of disruption to the frog skin microbiome, with higher values indicating higher stochasticity, a higher likelihood of containing transient pathogens, and perhaps, lower community stability [21,31,34,58].
We reduced our measures of water quality to two variables using principal component analysis (PCA) in JMP [59]. Water quality PC1 comprised 38% of the total variability and was positively associated with pH, orthophosphate, DO, and conductivity and negatively associated with ammonia. Water quality PC2 comprised 25.6% of the total variability and was positively associated with nitrate, ammonia, conductivity, and orthophosphate and negatively associated with DO. Thus, PC1 generally indicates changes in basic water physiochemistry, while PC2 indicates a gradient in nitrogen levels (Fig. S3).
To compare metrics of diversity (Shannon index, number of OTUs, Bray-Curtis dispersion) between host species, we ran independent general linear mixed models (GLMM) with sample site as a random effect.
We used model selection interfaces in JMP to parse out key variables explaining Bd infection and microbiome dispersion. We analyzed data separately for host species at the "site" level, averaging all pathogen, microbiome, and environmental variables by sample site (10 sites for A. crepitans, 9 for L. catesbeianus). We used Bd prevalence, loads, and microbiome dispersion (Bray-Curtis) separately as responses and included the following predictors in each model: number of OTUs on frogs, environmental bacteria PCo1 (Bray-Curtis), water quality PC1, water quality PC2, habitat disturbance (continuous %), and temperature. We looked at all possible models explaining Bd infection or microbiome dispersion with a limit of up to 5 predictors. We then selected the model with the lowest AICc score, and if more than one model was within 2 ΔAICc's of the top model, we selected the model with the highest R 2 coefficient (Table S2; [60,61]). We used a generalized linear model with a binomial distribution for the model on Bd prevalence and general linear models with a Gaussian distribution for Bd loads and microbiome dispersion. Pearson correlations were used to compare Bd infection and microbiome dispersion.
Using the same variables analyzed using model selection, we conducted structural equation models in R (piece-wiseSEM package; [62]) to test for significant associations while taking into account indirect effects among variables. We looked at links between Bd prevalence, microbiome dispersion, environmental bacteria PCo1, water quality PC1 and PC2, habitat disturbance (continuous %), and temperature. For each analysis, we started with the full model containing all ecologically relevant links and consecutively pruned the model to remove unsupported links and improve AICc model fit.
We used PermANOVA in R (adonis function in vegan package; [57]) to analyze differences in microbiome composition between sites with high (> 45% disturbed) and low (< 45% disturbed) habitat disturbance for A. crepitans, L. catesbeianus, and water samples separately. We also analyzed differences between frog and water bacterial communities. For each model, we controlled for sample site using the "strata" argument.
We calculated spatial connectivity between sites using the same landcover data acquired to measure habitat disturbance [50], and the ArcMap Circuitscape add-on (version 4.0.5; [63]). We coded disturbed habitat pixels with resistance of 5 (moderate resistance) and natural habitat pixels as 1 (low resistance). We then calculated pairwise resistance to movement between sites for both frog species and aquatic bacteria individually. We used partial mantel tests (mantel.partial function in the vegan package, [57]) to analyze the correlation between microbial similarity matrices (Bray-Curtis) and landscape resistance matrices for both species, controlling for Euclidean distances between sites.

Discussion
We found strong links between habitat disturbance, microbiome composition, and Bd infection in A. crepitans. Habitat disturbance was associated with higher dispersion of the skin microbiome, but also lower Bd prevalence. We originally predicted that habitat disturbance would either directly reduce Bd prevalence and loads as expected based on shifts in microclimate, or indirectly increase Bd prevalence through disruption to the skin microbiome. Our results indicate that despite the strong influence of habitat disturbance on microbiome dispersion, Bd prevalence was primarily driven by direct effects of disturbance on abiotic conditions, consistent with previous studies [48,49,64,65]. We also found lower Bd prevalence at sites with warmer water temperatures for L. catesbeianus. As water temperatures at our sites ranged from 20 to 28 °C, which exceeds the known thermal maximum of Bd in culture (25 °C; [66,67]), this is further support of a strong role of abiotic habitat in determining Bd infection dynamics.
Generally, Bd exhibits low pathogenicity in amphibians of the southeastern USA [68,69], indicating that frogs in this temperate climate are tolerant to infection or Bd in this system may have low virulence. This is further supported by high pathogen loads recorded in Acris sp. in past studies [70,71] and high Bd prevalence in our two aquatic-associated species without signs of disease. Even in the absence of Bd-associated disease, however, microbiome disruption through proxy of dispersion can be seen as an indicator of negative health in situations of environmental stress like habitat degradation. This is in line with the Anna Karenina principle put forth by Zaneveld et al. [31], proposing that microbiomes under stress will exhibit higher variability than microbiomes of healthy individuals. This principle has been demonstrated in many different taxa, including the coral species Aiptasia where heat stress was found to increase microbiome dispersion [36]. Similarly, another study found that temperature stress increased dispersion in the gut microbiome of tadpoles, which negatively impacted growth [21]. Furthermore, a recent study on endangered frogs from Costa Rica found skin microbiomes of frogs from sites with greater human disturbance had higher dispersion than those of frogs from natural sites [37]. Jiménez et al. [37] also found higher numbers of Bd-inhibitory OTUs on frogs at more disturbed sites, possibly indicating stress-associated selection on the microbiome to increase protective effects. Our findings reveal the potential for habitat disturbance to increase microbiome variability, providing support for the Anna Karenina principle when host health is subsequently affected. Our results were inconsistent between host species, likely due to morphological and life history differences between A. crepitans and L. catesbeianus. These two host species use the same habitat and are active during the same time of day, although Acris spp. are known to avoid forest cover [72]. Additionally, L. catesbeianus are largebodied and known to be resistant and tolerant to Bd infection [73][74][75]. Although not significant, we found a positive trend between water quality PC2, associated with nutrient While our metric of habitat disturbance was not directly correlated with nutrient levels, elevated nutrient levels in aquatic systems are often linked with disturbance-associated runoff [19]. Therefore, elevated nutrient levels in our study could indicate indirect habitat disturbance beyond our direct measurement of habitat composition. Past findings show that elevated nutrient runoff can lead to pathogenic bacterial blooms and subsequently elevated disease in frogs [33]. As dispersion is often used as a metric of dysbiosis [31], elevated microbiome dispersion in our study could point to similar health concerns.
Disturbance was also linked with distinct bacterial community assembly, indicating that in this instance higher dispersion signifies a shift in community composition in addition to increased stochasticity. While we found no relationship between microbes in frog versus water samples, we did find significant differences in aquatic microbes between sites with high and low disturbance. This indicates that disturbance could be influencing the microbial pools available for recruitment to frog skin microbiomes. As we only sampled the aquatic bacteria at three aquatic points within the 500-m radius, it is possible that other sources of bacteria, such as soil reservoirs, follow the expected trend with disturbance. Amphibians in the temperate southeast breed in late April through May, months when temperatures are still relatively mild. Based on temperature logging a week before sampling, average water temperatures at all sites (20-24 °C) were within the optimal range for Bd (15 and 25 °C; [66,67]). Through cyclic seasonal patterns in infection [76], with higher infections in the fall and spring months, amphibians in this area likely have immune responses primed for Bd infection [77].
We found strong support for microbiome differentiation based on habitat connectivity for both A. crepitans and L. catesbeianus, even after controlling for the relationship between microbiome dissimilarity and Euclidean distance between sites. These species are closely tied to bodies of water during the breeding season [78][79][80], which indicates that aquatic reservoirs of bacteria are the primary determinants of skin microbial composition. The link between microbiome similarity and habitat connectivity could therefore indicate that natural habitats are allowing greater connectivity between environmental reservoirs of bacteria [81], leading to greater similarity on frogs. Our analysis comparing aquatic bacterial similarity to habitat connectivity did not support this hypothesis, which indicates that both reservoir connectivity and host movement between habitats may be influential in microbiome assembly.
With exponentially increasing development of natural habitats, understanding the impacts of habitat disturbance on wildlife populations is vital. Past studies have shown that Bd infection risk is generally lower for generalist frogs that use disturbed habitats [48,49], but amphibians may still be negatively affected despite pathogen reduction in warmer open environments [82]. Our results suggest that even in the absence of disease, microbiome disruption can occur in conjunction with habitat disturbance. Other factors not measured in this study like immune function and stress response are likely influenced by habitat integrity [70,83] and may be mediating links between habitat disturbance and the skin microbiome. Additional experimental work tracking recruitment of environmental bacteria to the microbiome would be valuable in disentangling which habitat parameters are most likely connected to changes in the skin microbiome. Studies focusing on immune modulation by habitat change would also be valuable in determining how much of the microbiome response is mediated by the immune system and vice versa. Our work adds support to a growing body of literature showing that disturbance in many forms can influence the assembly of animal microbiomes in a stochastic manner [21,31,[36][37][38][39]. With so many contemporary threats facing amphibian populations, understanding the influences of habitat disturbance on all aspects of amphibian health is vital in detecting cryptic threats to seemingly stable populations.

Data and Code Availability
Raw data and the R code used for analyses are available online in supplemental files.