Xavier A. Harrison1,2*†
Stephen J. Price1,3
Kevin Hopkins1William T. M. Leung1
Chris Sergeant1
Trenton W. J. Garner1- 1Institute of Zoology, Zoological Society of London, London, United Kingdom
- 2Centre for Ecology and Conservation, University of Exeter, Exeter, United Kingdom
- 3UCL Genetics Institute, University College London, London, United Kingdom
Variation among animals in their host-associated microbial communities is increasingly recognized as a key determinant of important life history traits including growth, metabolism, and resistance to disease. Quantitative estimates of the factors shaping the stability of host microbiomes over time at the individual level in non-model organisms are scarce. Addressing this gap in our knowledge is important, as variation among individuals in microbiome stability may represent temporal gain or loss of key microbial species and functions linked to host health and/or fitness. Here we use controlled experiments to investigate how both heterogeneity in microbial species richness of the environment and exposure to the emerging pathogen Ranavirus influence the structure and temporal dynamics of the skin microbiome in a vertebrate host, the European common frog (Rana temporaria). Our evidence suggests that altering the bacterial species richness of the environment drives divergent temporal microbiome dynamics of the amphibian skin. Exposure to ranavirus effects changes in skin microbiome structure irrespective of total microbial diversity, but individuals with higher pre-exposure skin microbiome diversity appeared to exhibit higher survival. Higher diversity skin microbiomes also appear less stable over time compared to lower diversity microbiomes, but stability of the 100 most abundant (“core”) community members was similar irrespective of microbiome richness. Our study highlights the importance of extrinsic factors in determining the stability of host microbiomes over time, which may in turn have important consequences for the stability of host-microbe interactions and microbiome-fitness correlations.
Introduction
Animals are host to diverse communities of microbes, collectively referred to as the microbiome. Variation among individuals in their microbiomes has been linked to variation in host resistance to pathogens (Ford and King, 2016; King et al., 2016; Villarino et al., 2016; Antwis and Harrison, 2018; Warne et al., 2019), and disruption of the microbiome by external stressors (e.g., antibiotics) can have long term negative effects on host health (Theriot et al., 2014; Knutie et al., 2017; Warne et al., 2019). Though there is growing evidence that perturbation of the microbiome can have deleterious effects on host physiology, an understanding of the drivers of individual microbiome dynamics over time, and resistance to perturbation, remain relatively scarce in non-model organisms (Loudon et al., 2014b, 2016; Videvall et al., 2019). Addressing this shortfall in our knowledge is of fundamental importance to understanding the adaptive value of microbiomes for host health and fitness, as microbiome-health correlations may not be stable over time if microbiome flux represents loss of key microbial species and/or genes critical for optimal host physiology. Variation among individuals in their resistance to microbiome perturbation, and resilience following perturbation, could be a critical determinant of the distribution and stability of traits such as resistance to pathogens in natural populations.
The amphibian skin microbiome is rapidly becoming established as a model system for understanding the tripartite relationships between host, microbiome, and pathogens (e.g., Harris et al., 2009; Longo et al., 2015; Kueneman et al., 2016; Bates et al., 2018; Campbell et al., 2018a, b, 2019; Ross et al., 2019). Production of metabolites by skin-associated bacteria is a crucial component of immune defense against lethal fungal pathogens such as Batrachochytrium dendrobatidis (Bd) (e.g., Brucker et al., 2008; Becker et al., 2009) and Batrachochytrium salimandrivorans (Muletz-Wolz et al., 2017). Anti-fungal metabolite production by bacteria increases dramatically when the bacteria are co-cultured (Loudon et al., 2014a), suggesting that microbiome-mediated host protection is likely a function of synergistic interactions among community members. Greater microbiome diversity may therefore offer increased protection from pathogens (e.g., Piovia-Scott et al., 2017; Antwis and Harrison, 2018; Greenspan et al., 2019; but see Becker et al., 2019), but the ecological processes structuring and maintaining microbial diversity on amphibian skin remain relatively understudied, especially at the level of the individual (Loudon et al., 2014b, 2016; Longo and Zamudio, 2017; Hughey et al., 2019). For example, the diversity-stability hypothesis predicts that more diverse communities should be more resistant to disturbance, and several empirical studies support this hypothesis in plant community assemblages (McCann, 2000; Costello et al., 2012), but it is unclear whether this theory is also relevant at the scale of host-associated microbial communities (Costello et al., 2012, but see Koskella et al., 2017). Though several studies have sought to measure the influence of pathogenic infection on host microbiome structure (Jani and Briggs, 2014; Longo et al., 2015; Longo and Zamudio, 2017), investigations of whether the magnitude of microbiome disruption for infected hosts is modulated by initial microbiome state remains relatively scarce (see Jani et al., 2017; Jani and Briggs, 2018).
Here, we use experiments to examine how both the diversity of the environmental microbial reservoir and exposure to the lethal pathogen ranavirus influence skin microbial community dynamics in a native United Kingdom amphibian species, the European Common frog (Rana temporaria). The emerging infectious disease (EID) ranavirosis represent a significant threat to ectothermic vertebrate health, and infection with ranaviruses is associated with mass mortality, population extirpations and declines in biodiversity at a global scale (Jancovich et al., 2005; Fox et al., 2006; Bigarre et al., 2008; Ariel et al., 2009; Une et al., 2009; Whittington et al., 2010; Jensen et al., 2011; Allender et al., 2013; Earl and Gray, 2014; Price et al., 2014, 2019; Stark et al., 2014; Brunner et al., 2015; George et al., 2015; Miaud et al., 2016; Rijks et al., 2016; Rosa et al., 2017). Ranavirus was responsible for multi-species amphibian declines in continental Europe (Price et al., 2014), and of the common frog in the United Kingdom (Teacher et al., 2010), but also alters the age structure of remnant United Kingdom common frog populations (Campbell et al., 2018a). The frequency and severity of disease outbreaks are predicted to worsen alongside human-mediated range expansion of ranaviruses (Jancovich et al., 2005; Schloegel et al., 2009; Price et al., 2016, 2019). To manipulate environmental microbiome diversity, we assembled experimental units that either contained a complex natural bacterial reservoir (complex habitats, containing a soil substrate and leaf litter) or simplified ones (simple habitats, containing stony terrestrial substrates and no leaf litter). We performed two sequential experiments. In the first experiment, we group-housed 96 R. temporaria metamorphs in blocks of six individuals (n = 48 individuals per habitat treatment). For the second experiment, we individually housed 48 individuals in habitat treatments (n = 24 per habitat). Detailed experimental protocols are listed in the section “Materials and Methods” below. Both experiments allowed us to measure the influence of environmental microbiome on host microbiome structure and disruption of the host microbiome by pathogen exposure. Experiment 1 was designed to allow us to measure habitat-dependent mortality following exposure to ranavirus. Conversely, individual housing of frogs in Experiment 2 allowed us to track individual habitat- and pathogen-dependent microbiome trajectories over time, as well as within-individual changes in microbiome stability. Specifically, we sought to test whether (i) more diverse environmental bacterial reservoirs elicited more diverse frog skin microbiomes, (ii) more diverse skin microbiomes were more stable over time; and (iii) whether microbiome diversity predicted differences in resistance to ranavirus, manifesting as lower infection burdens and/or higher survival following exposure.
Materials and Methods
Ethics Statement
All experimental procedures and husbandry methods were approved by the ZSL Ethics Committee before any work was undertaken and was done under licensing by the United Kingdom Home Office (PPL 70/7830, P8897246A). Animal health and welfare was monitored daily during both the rearing and experimental periods and all animals were fed ad libitum (Tetra Tabimin for tadpoles, small crickets dusted with calcium and the Vetark Nutrobal vitamin supplement for metamorphosed frogs) throughout.
Experimental Protocols
Animal Rearing
Rana temporaria metamorphs were reared from tadpoles hatched from clutches sourced from United Kingdom garden ponds where ranavirosis had not been reported to the Garden Wildlife Health project1. Animals that completed metamorphosis were cohoused in large groups (no more than 30 per enclosure) in 460 × 300 × 170 mm Exo Terra Faunaria containing cleaned pea gravel, a large, cork bark cover object and sloped to accommodate a small aquatic area. Experimental animals were haphazardly selected from four group enclosures.
Preparation of Habitat Treatment Enclosures
The general layout of both habitat types was shared in that they both contained a filled, plastic PCR tip box (terrestrial platform) with a cover object, elevated above an aquatic area filled with aged tapwater and autoclaved pea gravel formed into a slope leading from the aquatic area to the platform. The two key differences were that; (i) the terrestrial platforms in complex habitats contained garden compost as a substrate, whilst the terrestrial platforms in simple habitats contained standard and autoclaved pea shingle and; (ii) leaf litter collected from Regents Park, London, was added to the aquatic area in the complex habitats. Complex habitat enclosures were left uncovered and outdoors for 2 weeks prior to the start of experiments, while simple habitat enclosures were prepared the day before frogs were transferred into replicates. During the experiment, uneaten cricket corpses were removed from simple habitat enclosures, but left in complex habitat enclosures. Experiment 1 comprised 16 replicate blocks, each housing six recently metamorphosed frogs (8 blocks/48 frogs per habitat treatment). Experiment 2 comprised 48 smaller units each housing an individual frog (24 frogs per habitat treatment). Following rearing in an outdoor facility, animals were moved to a procedure room and housed individually for 7 days in Perspex boxes with a cover object and damp paper towel as substrate to acclimatize to experimental conditions prior to any manipulations. Individuals were randomly assigned to experimental replicates and treatments (complex or simple habitats) using a script written in R (R Core Team, 2019). We note that our habitat manipulation altered both the bacterial richness in the environment and the structural composition of habitats (e.g., pea shingle vs. soil as a terrestrial substrate). Though this could have influenced the results, for example by changing the dynamics of host contact with the environment, these differences may be more representative of natural variation in microbiome-habitat relationships, where we would expect habitat heterogeneity to covary with microbiome structure.
Swabbing Protocols
For both experiments, we rinsed individuals in sterilized aged tap water to remove transient environmental microbes, and then swabbed the skin of the body and limbs of frogs with MW100 DrySwabs (Medical Wire Equipment, United Kingdom). In experiment 1, all animals were swabbed on Day 1 immediately preceding transfer to experimental units, then again on day 14, the latter referred to as the “pre-exposure” swab. Following the day 14 swab, we exposed individuals to either ranavirus or the control (see protocol below), and then swabbed all individuals again on Day 17 to measure the effect of ranavirus exposure (“post-exposure swab”). We swabbed all animals alive at the end of the experiment on Day 30, but do not include these data here as sample size per habitat-treatment group combination was low and unbalanced. For experiment 2, individuals were swabbed more frequently at Day 1, Day 7, Day 14 (pre-exposure), and Day 16 (post-exposure). For experiment 1 we present the pre- and post-exposure swabs as a 2-level time variable, whereas for experiment 2 we present all four time points as a time series.
Environmental swab samples (two per experimental unit, one terrestrial and one aquatic) were also collected on day 14 preceding pathogen exposure procedures. Terrestrial swabs were taken by running the swab over the terrestrial substrate and inside the cover objects twice. Aquatic swabs were taken by submerging the swab in the aquatic portion of the tank. These swabs allow us to assess how environmental microbiome diversity influences host skin microbiome diversity.
Ranavirus Exposure
Experimental units were randomly assigned to pathogen treatment group (ranavirus or sham) for both experiments using a script written in R. Prior to this, Ranavirus (FV3-like isolate RUK13, Cunningham et al., 2007) was cultured in EPC cells at 27°C, harvested after the cell layer had completely cleared, subjected to three rounds of freeze-thaw and then cleared of cells and cellular debris by centrifugation at 800g for 10 min and discarding the cell pellet. Virus titer was estimated using a 50% Tissue culture Infective Dose assay (TCID50) and calculated following the method of Reed and Muench (1938). Sham exposure media was produced by harvesting the supernatant of a pure culture of EPCs after the same 800g, 10 min spin. For exposures, animals were transferred either as co-housed groups (Experiment 1) or individually (Experiment 2) to 90 mm petri dishes containing 19 mL of aged tap water. Depending on treatment, either 1 mL of stock virus culture at 2 × 106 TCID50/mL (giving a final exposure concentration of 1 × 105 TCID50/mL) or 1 mL of sham media was added to the petri dish. Animals were exposed in petri dishes for 6 h before being returned to their habitat treatment enclosures. We used daily health and welfare checks throughout the experiment to monitor survival rates. We also used daily checks to monitor for signs of disease commonly associated with ranavirosis (see below: Price et al., 2016). We ended Experiment 1 on day 30 when all surviving frogs appeared physically healthy and when mortality had subsided, and Experiment 2 on day 16 following the post-exposure swab.
16S Sequencing and Bioinformatics
16S metagenetic library preparation was carried out using a modified version of the protocol detailed in Kozich et al. (2013) that amplifies the v4 section of the 16S rRNA gene. Sequencing was performed using 250 bp paired-end reads on an Illumina Miseq using a v2 chemistry 500 cycle cartridge (detailed information in Supplementary File “Detailed Amplicon Sequencing Methods”). Experiment 1 and 2 were processed on separate MiSeq runs, but all comparisons and statistical tests are made among samples within runs, so negating batch effects and inter-run variability. We processed raw 16S reads in the DADA2 pipeline (Callahan et al., 2016), using standard parameters as per the online tutorial. We used phyloseq (McMurdie and Holmes, 2013) for downstream sequence processing. In both experiments we removed amplicon sequence variants (ASVs) present in the no-template controls (Experiment 1: 54 ASVs of 11,640; Expt. 2: 466 of 14,963). To focus on differences in high abundance ASVs and to remove any potential bias introduced by small differences in low-abundance reads, we removed all ASVs from the dataset with fewer than 100 reads (e.g., Longo and Zamudio, 2017), leaving 5,796,063 reads of 1446 ASVs for Experiment 1, and 7,068,790 reads of 1969 ASVs for Experiment 2 used in downstream analysis. Reads per sample ranged from 6237–66993 (Experiment 1) to 16406–59607 (Experiment 2). We rarefied data to the minimum per-experiment sequencing depth prior to analysis.
Viral Load Quantification
Liver samples were extracted with DNeasy Blood & Tissue kits (Qiagen) following the manufacturer’s protocol. We quantified viral loads in all individuals using the qPCR method of Leung et al. (2017), which normalizes viral DNA quantities relative to host DNA in the sample.
Statistical Analysis
We conducted all statistical analyses in R. Due to differences in experimental design, the sets of analyses employed vary by experiment. For example, we did not track individual ID in group-housing in Experiment 1 and so do not examine drivers of within-individual changes in microbiome stability, but do so in Experiment 2. Likewise we did not assay survival in Experiment 2, but do present survival analyses for Experiment 1.
We fitted mixed effects models in the R package lme4 (Bates et al., 2015) and ranked competing models by AICc using the R package MuMIn (Bartoñ, 2019). We considered all models within six AICc units of the best supported AICc model to have relatively equal support in the data. To remove overly complex models from consideration we also applied the nesting rule (see Richards, 2008; Harrison et al., 2018) to remove models that were more complex versions of models with better AIC support. Where we refer to the “top model set,” we refer to the delta-6-AICc model set after the nesting rule has been applied. Where appropriate, we refitted models in a Bayesian framework using the Stan computational framework2 accessed with the brms package (Bürkner, 2017, 2018). The advantage of the Bayesian framework is that it allows quantification of uncertainty in parameters such as slopes and r2 values. Where appropriate, we specified mildly informative priors for parameters such as the correlation between random effects and slopes to speed up sampling and optimize convergence. We assessed convergence of chains using the Gelman-Rubin statistic, and inspected plots of posterior draws to verify adequate mixing of chains and sampling. Detailed descriptions of all statistical analyses and code are provided as an R Markdown document.
Diversity Indices
We calculated two metrics of alpha diversity: (i) richness as the exponent of the Shannon diversity index, also referred to as the effective number of species; and (ii) evenness, measured as the Shannon index divided by the log of the number of observed sequences in a sample. To derive measures of beta diversity, we performed Non-Metric Multidimensional Scaling (NMDS) ordinations on Bray–Curtis distance among bacterial community ASV abundances distances using the R package vegan (Oksanen et al., 2017). We also extracted NMDS1 values from these ordinations for analysis in statistical models (see below).
Experiment 1
We fitted a model containing the three-way interaction among time (pre- vs. post-exposure), habitat (Complex vs. Simple) and exposure (ranavirus vs. control) as predictors of alpha diversity, with separate models for richness and evenness. All models included a random intercept term for block ID (experimental tank) and used a Negative Binomial error structure. We performed PERMANOVA analysis in the R package vegan to test for differences among samples in beta diversity, also containing the time:habitat:pathogen three way interaction, and marginalizing the effect of block ID.
Experiment 2
We fitted a model containing day, day2 (to permit non-linear effects of time), habitat (Complex vs. Simple) and exposure (ranavirus vs. control) as well as an interaction between day2 and habitat as predictors of alpha diversity. All models included a random intercept for individual, and a random slope for day given individual. More complex models could not be fitted given the data available and produced convergence warnings. We used a Negative Binomial error structure for alpha diversity models to control for overdispersion (see Harrison, 2014) and a Gaussian error structure for beta diversity models. We performed PERMANOVA on Bray-Curtis distances among samples using the R package vegan to test for differences in beta diversity. We fitted a model containing a 3-way interaction between habitat, day and ranavirus exposure, permutated 999 times to derive p values for effects. We also fitted a linear mixed effects model to examine factors predicting NMDS1 variation among individuals, and included habitat, day, day2 and pathogen exposure as main effects, as well as habitat:pathogen exposure, habitat:day and habitat:day2 as predictors. All models included a random intercept for individuals. We could not include a random slope for day given individual as this produced convergence warnings.
Survival Analysis
We used the R package coxme (Therneau, 2015) to examine differences in survival dependent on habitat and pathogen exposure whilst controlling for block ID in Experiment 1. Sample size for this analysis was 85 individuals (42 in Simple Habitats and 43 in Complex habitats) across eight habitat blocks per habitat type. We censored eight individuals because they died prior to exposure. We ranked survival models by AICc to derive a top model set.
Predicted Functional Analysis
We used the ASV abundance matrices from Day 7 to predict functional profiles of microbial communities using PIPHILLIN (Iwai et al., 2016) and tested for differences in functional profiles dependent on habitat using Constrained Correspondence Analysis (CCA) in the R package vegan. We used the May 2017 release of the KEGG database and 97% identity cutoff. We visualized differences in predicted functional repertoire by plotting the axes of a CCA model fitted in vegan where we specified the two-level habitat predictor as the constrained variable. We performed predicted functional analysis only on Experiment 2 data as controlling for block effects in DESeq2 (Love et al., 2014) is difficult, and individual hosing of Experiment 2 obviates the need for this and so should more tightly control the false positive rate.
Stability Over Time
We calculated microbiome stability as the correlation between ASV abundances across two time points, following Lahti et al. (2014). That is, microbiome stability over time is estimated as the correlation between the two vectors of microbial community abundances from an individual for two time points, where stronger correlations indicate greater stability. We used Day 7 and 13 in Experiment 2 to quantify baseline stability, and tested variation in stability dependent on habitat using a t-test. We also calculated stability following exposure to a pathogen using the Day 13 and Day 16 ASV abundances and tested for a correlation between pre- and post-infection stability using Spearman’s correlation tests. We also calculated change in stability across the two time points by subtracting pre-infection stability from post infection stability.
We used ANOVA to test whether change in stability was explained by habitat or pathogen treatment. We repeated the above analyses restricting the dataset to the top 100 most abundant ASVs in each habitat to represent the “core” microbiome.
Results
Environment and Pathogen Exposure Modify Skin Microbiome Structure (Experiment 1)
Alpha Diversity
Bacterial richness and evenness of common frog skin was directly influenced by the complexity of the bacterial species reservoir in the environment. Individuals in habitats with higher environmental bacterial species richness possessed greater mean skin bacterial diversity (r = 0.82, p = 0.001; Figure 1A and Supplementary Figure S1). There was some evidence that overall effective number of species increased over time, but only weak evidence that this effect was dependent on habitat treatment (Figure 1B and Supplementary Table S1). There was no evidence of an effect of the interaction between time, habitat and pathogen exposure, or a main effect of ranavirus exposure on overall species richness of the microbiome (Supplementary Table S1). The top model investigating drivers of differences in community evenness contained interactions between ranavirus exposure and time, as well as habitat and time (Supplementary Table S2a). As for species richness, these results indicated that community evenness was lower in Simple Habitats prior to ranavirus exposure (Supplementary Table S2b). There was also some evidence that evenness increased over time for complex habitats, but differences due to ranavirus exposure were not clear (Supplementary Figure S2 and Supplementary Table S2b).
Figure 1. Environment modifies host skin bacterial community structure. (A) Effect of environmental bacterial richness on host skin microbiome diversity (Experiment 1, Time Point 2 [pre-exposure]). Individuals in habitats with more species rich environmental bacterial reservoirs also had higher skin bacterial richness (Complex) compared to individuals in habitats with lower diversity bacterial reservoirs (Simple). (B) There was no effect of exposure to ranavirus on mean levels of bacterial richness in either habitat type (Experiment 1, Time Point 3 [post-exposure]). (C) Skin bacterial community structure (beta diversity) differed significantly based on habitat and time (pre- vs. post-exposure). We detected a significant three-way interaction between time, habitat and pathogen exposure, suggesting that ranavirus exposure causes shifts in community structure dependent on habitat complexity. Note that no individuals were exposed to ranavirus in the Pre-Exposure panel, but individual points are colored by pathogen treatment (ranavirus vs. control) in both panels to allow comparison of groups across time.
Beta Diversity
PERMANOVA analysis controlling for block identified the dominant source of variation in microbial communities to be habitat (simple vs. complex, r2 = 15.9%, p = 0.001; Figure 1C). Exposure to ranavirus also influenced microbial community structure (infection main effect, r2 = 2.8%, p = 0.001), but critically operated via habitat:infection and infection:time interactions. The habitat:infection:time interaction was not significant (p = 0.08; Supplementary Table S3).
Survival Following Exposure to Ranavirus (Experiment 1)
Individuals in simple habitats exposed to ranavirus exhibited higher rates of mortality (68.4%) than individuals in complex habitats exposed to ranavirus (52.2%). The best-supported model contained effects of both habitat complexity and disease treatment on survival (Figure 2A and Supplementary Table S4). A model containing only disease treatment received marginally less support (ΔAICc = 0.22). Though the model containing the interaction between habitat and treatment was in the Δ6 AIC model set, it was a more complex version of a simpler model with better AIC support and so was removed under the nesting rule (Richards, 2008). Model averaged coefficients [and 95% confidence intervals] from the survival model were: Habitat 0.57 [−0.167,1.3] and ranavirus exposure 2.26 [1.2,3.33].
Figure 2. Survival following exposure to ranavirus. (A) Survival data for frogs exposed to Ranavirus or Control, in both Complex (red shading) and Simple (blue shading) habitats. The top model explaining variation in survival contained effects of both habitat and pathogen exposure. The model containing the habitat:pathogen interaction was not retained under the nesting rule. (B) Ranaviral infection loads following exposure to ranavirus, split by whether individuals died following exposure or were still alive on Day 30 at the close of the experiment. There was no difference in ranaviral infection burdens based on habitat treatment, but individuals that died had significantly higher infection loads. Three individuals in the control group exhibited weak infection.
There was no difference between habitats in likelihood of exhibiting gross signs of disease (Binomial GLMM, mean probability of exhibiting signs of disease [95% credible intervals]: complex 0.48 [0.11,0.82]; simple 0.5 [0.1,0.85]; pMCMC = 0.92) or in severity of visible signs of disease (Ordinal GLMM, mean probability of being scored category 0 [95% credible intervals]; complex 0.51 [0.12,0.9]; simple 0.46 [0.06,0.93]; pMCMC = 0.88). Individuals that died following exposure to ranavirus had higher viral loads than those that were still alive at the end of the experiment (Figure 2B). The best supported model examining variation in viral loads contained only the main effect of mortality (r2 38.7% [95% CI 15.1–56.1%]), as all other models with weaker support were removed under the nesting rule (Supplementary Table S5). Three individuals in the Control treatments died following exposure and exhibited weak ranavirus infections; these were inconsistent with the higher infection loads observed in other individuals that died after exposure (Figure 2B).
Environment Alters Host Microbiome Dynamics Over Time (Experiment 2)
Alpha Diversity
All individuals had similar bacterial species richness on Day 0 when they entered the experimental habitats, but the dynamics of host microbiome bacterial species difference over time differed markedly depending on habitat treatment (Figure 3A). The best supported model explaining differences in richness contained an interaction between day2 and habitat. When marginalizing the effects of time (sampling day) and variation among individuals in their change in diversity over time, individuals in complex habitats had greater skin bacterial diversity than those in simple habitats (Supplementary Figure S3). The top model explained 18.91% of variation in alpha diversity (95% credible interval 7.38–32.17%). There was no evidence that ranavirus exposure altered the dynamics of richness over time (Supplementary Tables S6, S7). When considering only the pre-exposure (Day 13) and post-exposure (Day 16) time points, the modal response was an increase in richness across the two time points. The top model examining factors predicting microbial community evenness contained only effects of day and day2, but no habitat main effect or interactions. The null model was also retained in the top model set (Supplementary Table S8). These data corroborate those from experiment 1 indicating a change in evenness over time. However, it is the environment (habitat) that appears to drive changes predominantly in the dynamics of microbial richness of amphibian skin over time.
Figure 3. Dynamics of frog skin bacterial communities over time. (A) Trends over time in bacterial alpha diversity (effective number of species) depending on Habitat treatment (rows) and pathogen exposure (columns). Dynamics of alpha diversity over time were significantly different in Complex habitats. When marginalizing the effects of time (day), individuals in Complex habitats possessed higher bacterial species richness compared to individuals in Simple habitats (see Supplementary Figure S2). (B) Trends over time in bacterial community structure (beta diversity). Individuals were exposed to ranavirus or control between days 13 and 16, but points are colored by disease treatment at all time points to allow tracking of beta diversity over time for different groups. (C) Temporal trends in primary axis of NMDS ordination (beta diversity) depending on Habitat treatment (rows) and pathogen exposure (columns). As with (A), there was strong support in the data for an interaction between day and habitat on beta diversity trajectories over time.
Beta Diversity
PERMANOVA performed on Bray–Curtis distances revealed significant effects of habitat, day and a habitat:day interaction (all p = 0.001) on bacterial beta diversity (Figure 3B and Supplementary Table S9). Collectively these terms explained roughly 15% of the variation in variation among individuals in bacterial community structure. There was no evidence that exposure to ranavirus modified the structure of bacterial communities (all interaction terms containing an effect of ranavirus exposure, p > 0.05), nor evidence of a ranavirus main effect (p = 0.28, Supplementary Table S9). Linear modeling of factors predicting NMDS1 revealed clear evidence of habitat-dependent variation in beta diversity trajectory over time for all individuals (Figure 3C). The only model in the top model set explaining predictors of NMDS1 contained an interaction between habitat type and day2 (Figure 3C and Supplementary Figure S4), corroborating the results of the PERMANOVA above.
Functional Traits
Predicted functional analysis using PIPHILLIN revealed distinct separation in the functional repertoires of the amphibian skin bacterial microbiome based on habitat after 7 days (CCA analysis, effect of habitat F(1,45) = 3.15, p = 0.01, Supplementary Figure S5). Analysis using DESeq2 revealed 12 pathways that were significantly more abundant in simple Habitats, and 8 pathways more abundant in complex habitats (Supplementary Table S10).
Viral Load Data
Viral loads of frogs following exposure to ranavirus in Experiment 2 were weak; mean viral load was 0.0013 viral copies per host cell [range 0.0001–0.01]. There was no significant difference in the mean viral load between animals in Complex and Simple habitats (Wilcoxon rank sum test, W = 88, p = 0.19, Supplementary Figure S6). No control animals in either habitat treatment had detectable levels of virus (Supplementary Figure S5).
Patterns of Microbiome Stability Varied by Habitat (Experiment 2)
When considering all ASVs, complex habitats exhibited decreased stability over time prior to infection when compared to simple habitats (t = 5.8, df = 43.3, p < 0.001; Figure 4A). However, an individual’s microbiome stability appeared consistent over time when comparing stability prior to pathogen exposure and stability following pathogen exposure (Figure 4B and Supplementary Figure S7A).
Figure 4. Microbiome stability over time. (A) Microbiome stability, measured as the correlation between ASV abundances across two times points, prior to pathogen exposure. Frogs in Simple habitats appear to have more stable microbial communities than those in Complex Habitats. (B) Scatterplot of microbiome stability over two sampling points prior to pathogen exposure (x axis) and two sampling points either side of pathogen exposure (y axis). Individual microbiome stability appears relatively consistent over time, irrespective of habitat or pathogen exposure. Dashed line represents 1:1 line of perfect correlation. Plots (C,D) are identical to plots (A,B), but use only the top 100 most abundant ASVs for each habitat type, representing a “core microbiome.” There is no difference between habitat types in stability for the core microbiome (C), and the correlation between stability values over time remains, though the relationship is weaker (D).
When considering only the top 100 ASVs, the difference in stability over time pre-infection was no longer apparent (t = 1.2, df = 42, p = 0.21; Figure 4C). Forty-seven ASVs from five Phyla were common to both sets of 100 most abundant ASVs by habitat (Supplementary Table S11). The most common Phylum of shared ASVs was Proteobacteria, comprising 32 of the 47 ASVs (68%). Actinobacteria and Bacteroidetes accounted for 13% each of the shared ASV taxonomy. At the genus level, notable shared ASVs were classified as Citrobacter, Acinetobacter, Chryseobacterium, and Stenotrophomonas, all of which have been associated with production of metabolites that inhibit other amphibian pathogens like B. dendrobatidis (e.g., Antwis and Harrison, 2018).
Both habitats still exhibit consistent levels of stability either side of exposure to the pathogen (Figure 4D and Supplementary Figure S7B), though the correlation is weaker. There was no evidence that habitat treatment, pathogen exposure or their interaction affected the magnitude of change in stability over time, for either all ASVs or the analysis restricted to the top 100 most abundant ASVs (ANOVA, all p > 0.27). Individual trajectories of microbiome richness and beta diversitty are provided in Supplementary Figures S8, S9.
Discussion
Though our knowledge of the factors shaping the structure of the host-associated microbiota is increasing, studies directed at understanding the predictors of longitudinal variation of the microbiome in non-model organisms are relatively scarce (e.g., Smith et al., 2015; Videvall et al., 2019). Our results from two experiments suggest that the structure and temporal dynamics of the amphibian skin microbiome are influenced by both the environment and exposure to a lethal pathogen. Overall alpha diversity of the microbiome appeared to influence temporal stability, where more “species-rich” microbial communities were less stable over time compared to less diverse communities. Crucially, this effect disappeared when considering only the top 100 most abundant bacterial taxa, suggesting “core microbiome” stability may be relatively uniform irrespective of total diversity. Finally, our survival data suggest that higher skin microbiome diversity may correlate with greater survival following exposure to the lethal pathogen ranavirus. Our results have important implications for our understanding of factors driving variation among individuals in the stability of both their microbiomes and the strength of host-microbe interactions over time, and in turn how both traits may be compromised by external stressors such as exposure to pathogens.
Environmental Microbial Diversity Influence Temporal Microbiome Dynamics
By manipulating the microbial reservoir in the environment, we elicited differential patterns of microbiome diversity on the skin of common frogs. Microbial communities differing in diversity also exhibited distinctive signatures of change over time. Higher diversity skin microbiomes appeared less stable over time, an effect driven primarily by weak correlations over time in the abundances of rarer bacterial taxa. The “core microbiome” of the most abundant ASVs in each habitat type appeared stable irrespective of overall diversity. Most strikingly, microbiome stability itself appeared conserved over time: stability between the first two time points correlated strongly with the two time points bracketing exposure to the pathogen. Our data therefore support the idea of consistent variation among individuals in microbiome stability over time, where structure itself is a function of the environment that an individual inhabits. These data support previous work indicating that environmental context and complexity is a key determinant of the assembly and stochasticity dynamics of host microbiomes, with important consequences for host resistance ot disease (Becker et al., 2017). It is notable that 47 of the 100 most abundant per-habitat ASVs were common to both habitats, alluding to a constrained core microbiome structure irrespective of habitat microbial diversity and structure.
A major outstanding question is what are the consequences of consistent within and among-individual variation in microbiome stability over time? Predicted functional analysis highlighted that differences in overall microbial community diversity also reflected differences in functional repertoire (see also Bletz et al., 2016) which could reflect variation in the presence or strength of key interactions between microbes and hosts such as production of antimicrobial metabolites that defend the host from pathogens (e.g., Brucker et al., 2008; Becker et al., 2009; Antwis and Harrison, 2018). Measures of temporal stability or stochasticity of microbiomes, and the processes that drive these traits, are critical for understanding the consistency of microbe-mediated functions over time. Most of the data we have on within-individual microbiome dynamics and stability come from human or model organism studies (e.g., Fink et al., 2013; Kelly et al., 2016; Gilbert et al., 2018; Schirmer et al., 2018), with relatively few on non-model organisms (e.g., Antwis et al., 2019; Videvall et al., 2019), and fewer still from controlled experiments (Kueneman et al., 2016; Grottoli et al., 2018). Though several studies have measured seasonal dynamics of microbiome in species such as aphids (Smith et al., 2015) and mosquitoes (Novakova et al., 2017); they rely on population-based metrics of microbiome structure that may mask substantial among-individual variation in microbiome dynamics. Amphibians and their skin microbiomes provide a model for understanding the processes shaping the forces of colonization, competition and coexistence of microbial species on a vertebrate host, and quantifying the emergent functional properties of these microbial communities and their consequences for the host. The properties of this system make it well suited to testing the applicability of established ecological theory derived from eukaryotic communities to prokaryotic assemblages associated with animals, plants and soils, including the relationship between biodiversity and ecosystem function (e.g., see Koskella et al., 2017; Greenspan et al., 2019).
Exposure to Ranavirus Disrupts the Host Skin Microbiome
Our data from experiment 1 revealed that exposure to ranavirus elicited subtle but significant changes to the structure of the amphibian skin microbiome after 48 h. We predicted that more diverse microbial communities should be more resistant to perturbation by the ranavirus, but our data suggest that the skin microbiomes of individuals in both habitat treatments were affected by the pathogen. This supports previous work showing that pathogens like Bd can destabilize host microbiomes (e.g., Jani and Briggs, 2014; Walke et al., 2015; Longo and Zamudio, 2017). Though we didn’t detect a similar time:habitat:pathogen interaction in Experiment 2, this can be explained by the relatively low infection burdens in this experiment. Results from our experimental work here are well supported by counterpart investigations into the structure of the microbiota of wild common frogs. These studies have illustrated distinct differences in bacterial community structure at sites suffering mass mortality events due to ranavirus compared to sites where no such outbreaks have been detected (Campbell et al., 2018b, 2019), even after accounting for differences among populations (Campbell et al., 2019). These correlative data from wild frogs could represent bacterial communities in some populations associated with protection of the host from viral infection, or marked shifts in microbiome structure in populations suffering ranaviral infection. Use of pesticides has been associated with increased prevalence of ranaviruses (North et al., 2015) and could theoretically be mediated by disruption of the both environmental and host-associated bacterial communities. Considering these data with results from our experiments hints that both processes may be responsible for the patterns observed in nature. Disruption of the host microbiome by pathogens of wild vertebrates is likely to be far more common than the existing literature suggests. The scarcity of studies directed at quantifying microbiome disruption by pathogens means we currently lack the ability to compare the magnitude of the perturbation effect among host species and both host and pathogen taxonomic groups.
Links Between Microbiome and Survival Following Ranavirus Exposure
Our controlled infection experiment revealed that individuals with less diverse microbiomes exhibited higher mortality following exposure to ranavirus compared to individuals with higher diversity microbial communities, consistent with our predictions. In our models, greater resistance to pathogenic infection as an emergent property of microbiome diversity would be evidenced by a diversity (habitat) by pathogen exposure interaction term. We note that though there was reasonable support for a model containing this interaction in our top model set, it was not retained under the nesting rule. As such, there exists some model section uncertainty regarding the effect of microbiome diversity on resistance to ranavirus infection. Several studies have provided evidence consistent with a correlation between overall microbiome diversity and susceptibility to infectious disease and costs associated with host responses to pathogen exposure (e.g., Cariveau et al., 2014; Kueneman et al., 2016), though these effects are not always consistent (e.g., Becker et al., 2019; Ma et al., 2019). Disruption of the normal microbiome by administration of antibiotics to laboratory mice can permit successful infection of Clostridium difficile (Theriot et al., 2014), loss of microbiome diversity in amphibians can increase susceptibility to the fungal pathogen Bd (Kueneman et al., 2016), and disruption of the microbiome in early life can increase downstream susceptibility to parasites (Knutie et al., 2017). Notably, augmentation of low diversity skin microbiomes with key taxa from the more diverse wild-type microbiome can reverse the observed increase in susceptibility to a lethal pathogen like Bd (Kueneman et al., 2016). Our habitat treatments differed in overall physical structure as well as microbial diversity, as complex habitats contained different terrestrial substrate as well as leaf litter in the water. Traits such as host microbiome richness, temporal dynamics, and resistance to disease will be governed by both extrinsic processes (microbial diversity present in the environment capable of colonizing the host) and intrinsic factors, such as host immunogenetic variation (e.g., Bolnick et al., 2014) and physiological stress (e.g., Lokmer and Wegner, 2015). Though variation in host stress due to structural heterogeneity between habitat treatments could have influenced microbiome dynamics, we believe this effect would be minimal in our data as all individuals were reared under captive conditions as tadpoles, and so were acclimated to conditions found in the simple habitat treatments. Nevertheless, future work will standardize the environmental structure of the habitats and manipulate only the microbial reservoir to remove the potential for such differences between treatments.
The mechanisms underpinning diversity-disease relationships in amphibians warrant further investigation. Microbiome diversity alone cannot be considered a beneficial trait for hosts; rather diversity itself is an emergent property of ecological processes playing out within the host (Shade, 2017) that underpin the true mechanism. More diverse microbiomes could be more likely to contain species producing antiviral compounds such as bacteriocins (see Drider et al., 2016), or to prime the host immune system to produce anti-microbial peptides (Woodhams et al., 2019) that can inactivate ranavirus virions (Chinchar et al., 2004). As expected, ranavirus-exposed individuals that died during the experiment had higher viral loads than those that survived. Higher survival in individuals with more diverse microbiomes could represent microbe-mediated defense preventing infection burdens from reaching lethal thresholds. Indeed, our predicted functional analysis of skin bacterial microbiomes revealed distinct differences dependent on diversity. Though complex skin microbiomes were predicted to differ in relative abundance of pathways linked to human viral infections, the relevance of such differences to amphibian defense against ranavirus remains to be determined. An important priority for future work is to quantify the true functional genetic repertoire of amphibian skin microbiomes to permit identification of potential metabolic pathways linked to disease, and examine how their relative abundance changes in concert with overall microbiome diversity and microbial species composition. Addressing this knowledge gap requires integration of further ‘omic tools such as shotgun metagenomics and metabolomics with more common amplicon sequencing metagenetics (Rebollar et al., 2016). Finally, Warne et al. (2019) recently showed that disruption of the gut microbiome in early life can influence host metabolism and susceptibility to ranavirus in later life. Given that the oral cavity and alimentary canal are major routes of infection for ranaviruses (e.g., Robert et al., 2011; Saucedo et al., 2019), one possibility is that measurements of skin microbiome diversity in amphibians are also reflective of gut microbiome diversity. Future work should quantify this covariation between multi-site microbiome dynamics and seek to understand the functional consequences of increased skin and gut microbiome diversity in hosts vulnerable to ranavirus.
Author’s Note
This manuscript has been released as a preprint at bioRxiv (Harrison et al., 2017).
Data Availability Statement
R markdown scripts and data sets permitting full reproduction of all analyses are hosted on FigShare at https://doi.org/10.6084/m9.figshare.c.4607198. Sequences have been uploaded to the NCBI Sequence Read Archive under BioProject accession numbers PRJNA559513 (Experiment 1) and PRJNA559522 (Experiment 2).
Ethics Statement
The animal study was reviewed and approved by the Institute of Zoology Ethics Committee.
Author Contributions
XH and TG designed the experiment. XH, SP, WL, CS, and TG conducted the experiment. XH and KH sequenced the microbial data. SP and WL produced ranaviral infection load data. XH analyzed the data. XH, SP, and TG wrote the manuscript with input from all authors.
Funding
XH was supported by an Institute of Zoology Independent Research Fellowship and a Royal Society Research grant (RG130550). TG acknowledges funding from NERC (NE/M000338/1) that supported this research.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
XH thanks Matt Gray for helpful discussion about the manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.02883/full#supplementary-material
Footnotes
References
Allender, M. C., Mitchell, M. A., McRuer, D., Christian, S., and Byrd, J. (2013). Prevalence, clinical signs, and natural history characteristics of frog virus 3-like infections in eastern box turtles (Terrapene carolina carolina). Herpetol. Conserv. Biol. 8, 308–320.
Antwis, R. E., Edwards, K. L., Unwin, B., Walker, S. L., and Shultz, S. (2019). Rare gut microbiota associated with breeding success, hormone metabolites and ovarian cycle phase in the critically endangered eastern black rhino. Microbiome 7:27. doi: 10.1186/s40168-019-0639-0
Antwis, R. E., and Harrison, X. A. (2018). Probiotic consortia are not uniformly effective against different amphibian chytrid pathogen isolates. Mol. Ecol. 27, 577–589. doi: 10.1111/mec.14456
Ariel, E., Kielgast, J., Svart, H. E., Larsen, K., Tapiovaara, H., Jensen, B. B., et al. (2009). Ranavirus in wild edible frogs Pelophylax kl. esculentus in Denmark. Dis. Aquat. Organ. 85, 7–14. doi: 10.3354/dao02060
Bates, D., Maechler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48. doi: 10.18637/jss.v067.i01
Bates, K. A., Clare, F. C., O’Hanlon, S., Bosch, J., Brookes, L., Hopkins, K., et al. (2018). Amphibian chytridiomycosis outbreak dynamics are linked with host skin bacterial community structure. Nat. Commun. 9:693. doi: 10.1038/s41467-018-02967-w
Becker, C. G., Bletz, M. C., Greenspan, S. E., Rodriguez, D., Lambertini, C., Jenkinson, T. S., et al. (2019). Low-load pathogen spillover predicts shifts in skin microbiome and survival of a terrestrial-breeding amphibian. Proc. R. Soc. B 286:20191114. doi: 10.1098/rspb.2019.1114
Becker, C. G., Longo, A. V., Haddad, C. F., and Zamudio, K. R. (2017). Land cover and forest connectivity alter the interactions among host, pathogen and skin microbiome. Proc. R. Soc. B Biol. Sci. 284:20170582. doi: 10.1098/rspb.2017.0582
Becker, M. H., Brucker, R. M., Schwantes, C. R., Harris, R. N., and Minbiole, K. P. (2009). The bacterially produced metabolite violacein is associated with survival of amphibians infected with a lethal fungus. Appl. Environ. Microbiol. 75, 6635–6638. doi: 10.1128/AEM.01294-09
Bigarre, L., Cabon, J., Baud, M., and Castric, J. (2008). Ranaviruses associated with high mortalities in catfish in France. Bull. Eur. Assoc. Fish Pathol. 28, 163–168.
Bletz, M. C., Goedbloed, D. J., Sanchez, E., Reinhardt, T., Tebbe, C. C., Bhuju, S., et al. (2016). Amphibian gut microbiota shifts differentially in community structure but converges on habitat-specific predicted functions. Nat. Commun. 7:13699. doi: 10.1038/ncomms13699
Bolnick, D. I., Snowberg, L. K., Caporaso, J. G., Lauber, C., Knight, R., and Stutz, W. E. (2014). Major histocompatibility C omplex class II b polymorphism influences gut microbiota composition and diversity. Mol. Ecol. 23, 4831–4845. doi: 10.1111/mec.12846
Brucker, R. M., Harris, R. N., Schwantes, C. R., Gallaher, T. N., Flaherty, D. C., Lam, B. A., et al. (2008). Amphibian chemical defense: antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus. J. Chem. Ecol. 34, 1422–1429. doi: 10.1007/s10886-008-9555-7
Brunner, J. L., Storfer, A., Gray, M. J., and Hoverman, J. T. (2015). “Ranavirus ecology and evolution: from epidemiology to extinction,” in Ranaviruses, eds M. J. Gray, and V. G. Chinchar, (Berlin: Springer International Publishing), 71–104. doi: 10.1007/978-3-319-13755-1_4
Bürkner, P. (2017). brms: an R Package for bayesian multilevel models using stan. J. Statist. Softw. 80, 1–28. doi: 10.18637/jss.v080.i01
Bürkner, P. (2018). Advanced bayesian multilevel modeling with the R Package brms. R J. 10, 395–411. doi: 10.32614/RJ-2018-017
Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., and Holmes, S. P. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13:581. doi: 10.1038/nmeth.3869
Campbell, L. J., Garner, T. W., Tessa, G., Scheele, B. C., Griffiths, A. G., Wilfert, L., et al. (2018a). An emerging viral pathogen truncates population age structure in a European amphibian and may reduce population viability. PeerJ 6:e5949. doi: 10.7717/peerj.5949
Campbell, L. J., Hammond, S. A., Price, S. J., Sharma, M. D., Garner, T. W., Birol, I., et al. (2018b). A novel approach to wildlife transcriptomics provides evidence of disease-mediated differential expression and changes to the microbiome of amphibian populations. Mol. Ecol. 27, 1413–1427. doi: 10.1111/mec.14528
Campbell, L. J., Garner, T. W. J., Hopkins, K., Griffiths, A. G. F., and Harrison, X. A. (2019). Outbreaks of an emerging viral disease covary with differences in the composition of the skin microbiome of a wild UK amphibian. Front. Microbiol. 10:1245. doi: 10.3389/fmicb.2019.01245
Cariveau, D. P., Powell, J. E., Koch, H., Winfree, R., and Moran, N. A. (2014). Variation in gut microbial communities and its association with pathogen infection in wild bumble bees (Bombus). ISME J. 8, 2369–2379. doi: 10.1038/ismej.2014.68
Chinchar, V. G., Bryan, L., Silphadaung, U., Noga, E., Wade, D., and Rollins-Smith, L. (2004). Inactivation of viruses infecting ectothermic animals by amphibian and piscine antimicrobial peptides. Virology 323, 268–275. doi: 10.1016/j.virol.2004.02.029
Costello, E. K., Stagaman, K., Dethlefsen, L., Bohannan, B. J., and Relman, D. A. (2012). The application of ecological theory toward an understanding of the human microbiome. Science 336, 1255–1262. doi: 10.1126/science.1224203
Cunningham, A. A., Hyatt, A. D., Russell, P., and Bennett, P. M. (2007). Emerging epidemic diseases of frogs in Britain are dependent on the source of Ranavirus agent and the route of exposure. Epidemiol. Infect. 135, 1200–1212. doi: 10.1017/s0950268806007679
Drider, D., Bendali, F., Naghmouchi, K., and Chikindas, M. L. (2016). Bacteriocins: not only antibacterial agents. Probiotics Antimicrob. Proteins 8, 177–182. doi: 10.1007/s12602-016-9223-0
Earl, J. E., and Gray, M. J. (2014). Introduction of Ranavirus to isolated wood frog populations could cause local extinction. Ecohealth 11, 581–592. doi: 10.1007/s10393-014-0950-y
Fink, C., Staubach, F., Kuenzel, S., Baines, J. F., and Roeder, T. (2013). Noninvasive analysis of microbiome dynamics in the fruit fly Drosophila melanogaster. Appl. Environ. Microbiol. 79, 6984–6988. doi: 10.1128/AEM.01903-13
Ford, S. A., and King, K. C. (2016). Harnessing the power of defensive microbes: evolutionary implications in nature and disease control. PLoS Pathog. 12:e1005465. doi: 10.1371/journal.ppat.1005465
Fox, S. F., Greer, A. L., Torres-Cervantes, R., and Collins, J. P. (2006). First case of Ranavirus-associated morbidity and mortality in natural populations of the South American frog Atelognathus patagonicus. Dis. Aquat. Organ. 72, 87–92. doi: 10.3354/dao072087
George, M. R., John, K. R., Mansoor, M. M., Saravanakumar, R., Sundar, P., and Pradeep, V. (2015). Isolation and characterization of a Ranavirus from koi, Cyprinus carpio L., experiencing mass mortalities in India. J. Fish Dis. 38, 389–403. doi: 10.1111/jfd.12246
Gilbert, J. A., Blaser, M. J., Caporaso, J. G., Jansson, J. K., Lynch, S. V., and Knight, R. (2018). Current understanding of the human microbiome. Nat. Med. 24:392. doi: 10.1038/nm.4517
Greenspan, S. E., Lyra, M. L., Migliorini, G. H., Kersch-Becker, M. F., Bletz, M. C., Lisboa, C. S., et al. (2019). Arthropod–bacteria interactions influence assembly of aquatic host microbiome and pathogen defense. Proc. R. Soc. B. 286:20190924. doi: 10.1098/rspb.2019.0924
Grottoli, A. G., Dalcin Martins, P., Wilkins, M. J., Johnston, M. D., Warner, M. E., Cai, W.-J., et al. (2018). Coral physiology and microbiome dynamics under combined warming and ocean acidification. PLoS One 13:e0191156. doi: 10.1371/journal.pone.0191156
Harris, R. N., Brucker, R. M., Walke, J. B., Becker, M. H., Schwantes, C. R., Flaherty, D. C., et al. (2009). Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824. doi: 10.1038/ismej.2009.27
Harrison, X. A. (2014). Using observation-level random effects to model overdispersion in count data in ecology and evolution. PeerJ 2:e616. doi: 10.7717/peerj.616
Harrison, X. A., Donaldson, L., Correa-Cano, M. E., Evans, J., Fisher, D. N., Goodwin, C. E., et al. (2018). A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6:e4794. doi: 10.7717/peerj.4794
Harrison, X. A., Price, S. J., Hopkins, K., Leung, W., Sergeant, C., and Garner, T. W. J. G. (2017). Host microbiome richness predicts resistance to disturbance by pathogenic infection in a vertebrate host. bioRxiv [Preprint].
Hughey, M. C., Sokol, E. R., Walke, J. B., Becker, M. H., and Belden, L. K. (2019). Ecological correlates of large-scale turnover in the dominant members of Pseudacris crucifer skin bacterial communities. Microb. Ecol. 4, 1–11. doi: 10.1007/s00248-019-01372-0
Iwai, S., Weinmaier, T., Schmidt, B. L., Albertson, D. G., Poloso, N. J., Dabbagh, K., et al. (2016). Piphillin: improved prediction of metagenomic content by direct inference from human microbiomes. PLoS One 11:e0166104. doi: 10.1371/journal.pone.0166104
Jancovich, J. K., Davidson, E. W., Parameswaran, N., Mao, J., Chinchar, V. G., Collins, J. P., et al. (2005). Evidence for emergence of an amphibian iridoviral disease because of human-enhanced spread. Mol. Ecol. 14, 213–224. doi: 10.1111/j.1365-294x.2004.02387.x
Jani, A. J., and Briggs, C. J. (2014). The pathogen Batrachochytrium dendrobatidis disturbs the frog skin microbiome during a natural epidemic and experimental infection. Proc. Natl. Acad. Sci. U.S.A. 111, E5049–E5058. doi: 10.1073/pnas.1412752111
Jani, A. J., and Briggs, C. J. (2018). Host and aquatic environment shape the amphibian skin microbiome but effects on downstream resistance to the pathogen Batrachochytrium dendrobatidis are variable. Front. Microbiol. 9:487. doi: 10.3389/fmicb.2018.00487
Jani, A. J., Knapp, R. A., and Briggs, C. J. (2017). Epidemic and endemic pathogen dynamics correspond to distinct host population microbiomes at a landscape scale. Proc. R. Soc. B 284:20170944. doi: 10.1098/rspb.2017.0944
Jensen, B. B., Holopainen, R., Tapiovaara, H., and Ariel, E. (2011). Susceptibility of pike-perch Sander lucioperca to a panel of Ranavirus isolates. Aquaculture 313, 24–30. doi: 10.1016/j.aquaculture.2011.01.036
Kelly, B. J., Imai, I., Bittinger, K., Laughlin, A., Fuchs, B. D., Bushman, F. D., et al. (2016). Composition and dynamics of the respiratory tract microbiome in intubated patients. Microbiome 4:7. doi: 10.1186/s40168-016-0151-8
King, K. C., Brockhurst, M. A., Vasieva, O., Paterson, S., Betts, A., Ford, S. A., et al. (2016). Rapid evolution of microbe-mediated protection against pathogens in a worm host. ISME J. 10:1915. doi: 10.1038/ismej.2015.259
Knutie, S. A., Wilkinson, C. L., Kohl, K. D., and Rohr, J. R. (2017). Early-life disruption of amphibian microbiota decreases later-life resistance to parasites. Nat. Commun. 8:86. doi: 10.1038/s41467-017-00119-0
Koskella, B., Hall, L. J., and Metcalf, C. J. E. (2017). The microbiome beyond the horizon of ecological and evolutionary theory. Nat. Ecol. Evol. 1:1606. doi: 10.1038/s41559-017-0340-2
Kozich, J. J., Westcott, S. L., Baxter, N. T., Highlander, S. K., and Schloss, P. D. (2013). Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl. Env. Microbiol. 79, 5112–5120. doi: 10.1128/AEM.01043-13
Kueneman, J. G., Woodhams, D. C., Harris, R., Archer, H. M., Knight, R., and McKenzie, V. J. (2016). Probiotic treatment restores protection against lethal fungal infection lost during amphibian captivity. Proc. R. Soc. B 283:20161553. doi: 10.1098/rspb.2016.1553
Lahti, L., Salojärvi, J., Salonen, A., Scheffer, M., and De Vos, W. M. (2014). Tipping elements in the human intestinal ecosystem. Nat. Commun. 5:4344. doi: 10.1038/ncomms5344
Leung, W. T., Thomas-Walters, L., Garner, T. W., Balloux, F., Durrant, C., and Price, S. J. (2017). A quantitative-PCR based method to estimate Ranavirus viral load following normalisation by reference to an ultraconserved vertebrate target. J. Virol. Methods 249, 147–155. doi: 10.1016/j.jviromet.2017.08.016
Lokmer, A., and Wegner, K. M. (2015). Hemolymph microbiome of Pacific oysters in response to temperature, temperature stress and infection. ISME J. 9, 670–682. doi: 10.1038/ismej.2014.160
Longo, A. V., Savage, A. E., Hewson, I., and Zamudio, K. R. (2015). Seasonal and ontogenetic variation of skin microbial communities and relationships to natural disease dynamics in declining amphibians. R. Soc. Open Sci. 2:140377. doi: 10.1098/rsos.140377
Longo, A. V., and Zamudio, K. R. (2017). Environmental fluctuations and host skin bacteria shift survival advantage between frogs and their fungal pathogen. ISME J. 11, 349–361. doi: 10.1038/ismej.2016.138
Loudon, A. H., Venkataraman, A., Van Treuren, W., Woodhams, D. C., Parfrey, L. W., and McKenzie, V. J. (2016). Vertebrate hosts as islands: dynamics of selection, immigration, loss, persistence, and potential function of bacteria on salamander skin. Front. Microbiol. 7:333. doi: 10.3389/fmicb.2016.00333
Loudon, A. H., Holland, J. A., Umile, T. P., Burzynski, E. A., Minbiole, K. P., and Harris, R. N. (2014a). Interactions between amphibians’ symbiotic bacteria cause the production of emergent anti-fungal metabolites. Front. Microbiol. 5:441. doi: 10.3389/fmicb.2014.00441
Loudon, A. H., Woodhams, D. C., Parfrey, L. W., Archer, H., Knight, R., McKenzie, V., et al. (2014b). Microbial community dynamics and effect of environmental microbial reservoirs on red-backed salamanders (Plethodon cinereus). ISME J. 8, 830–840. doi: 10.1038/ismej.2013.200
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 15:550.
Ma, Z. S., Li, L., and Gotelli, N. J. (2019). Diversity-disease relationships and shared species analyses for human microbiome-associated diseases. ISME J. 13, 1911–1919. doi: 10.1038/s41396-019-0395-y
McMurdie, P. J., and Holmes, S. (2013). phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217. doi: 10.1371/journal.pone.0061217
Miaud, C., Pozet, F., Gaudin, N. C. G., Martel, A., Pasmans, F., and Labrut, S. (2016). Ranavirus causes mass die-offs of alpine amphibians in the southwestern alps (France). J. Wildl. Dis. 52, 242–252. doi: 10.7589/2015-05-113
Muletz-Wolz, C. R., Almario, J. G., Barnett, S. E., DiRenzo, G. V., Martel, A., Pasmans, F., et al. (2017). Inhibition of fungal pathogens across genotypes and temperatures by amphibian skin bacteria. Front. Microbiol. 8:1551. doi: 10.3389/fmicb.2017.01551
North, A. C., Hodgson, D. J., Price, S. J., and Griffiths, A. G. F. (2015). Anthropogenic and ecological drivers of amphibian disease (Ranavirosis). PLoS One 10:e0127037. doi: 10.1371/journal.pone.0127037
Novakova, E., Woodhams, D. C., Rodríguez-Ruano, S. M., Brucker, R. M., Leff, J. W., Maharaj, A., et al. (2017). Mosquito microbiome dynamics, a background for prevalence and seasonality of West Nile virus. Front. Microbiol. 8:526. doi: 10.3389/fmicb.2017.00526
Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O’Hara, R. B., et al. (2017). vegan: Community Ecology Package. R package version 2.4-3.
Piovia-Scott, J., Rejmanek, D., Woodhams, D. C., Worth, S. J., Kenny, H., McKenzie, V., et al. (2017). Greater species richness of bacterial skin symbionts better suppresses the amphibian fungal pathogen Batrachochytrium dendrobatidis. Microb. Ecol. 74, 1–10. doi: 10.1007/s00248-016-0916-4
Price, S. J., Garner, T. W., Cunningham, A. A., Langton, T. E., and Nichols, R. A. (2016). Reconstructing the emergence of a lethal infectious disease of wildlife supports a key role for spread through translocations by humans. Proc. R. Soc. B 283:20160952. doi: 10.1098/rspb.2016.0952
Price, S. J., Garner, T. W., Nichols, R. A., Balloux, F., Ayres, C., de Alba, A. M., et al. (2014). Collapse of amphibian communities due to an introduced Ranavirus. Curr. Biol. 24, 2586–2591. doi: 10.1016/j.cub.2014.09.028
Price, S. J., Leung, W. T., Owen, C. J., Puschendorf, R., Sergeant, C., Cunningham, A. A., et al. (2019). Effects of historic and projected climate change on the range and impacts of an emerging wildlife disease. Glob. Chang. Biol. 25, 2648–2660. doi: 10.1111/gcb.14651
R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Rebollar, E. A., Antwis, R. E., Becker, M. H., Belden, L. K., Bletz, M. C., Brucker, R. M., et al. (2016). Using “Omics” and integrated multi-omics approaches to guide probiotic selection to mitigate chytridiomycosis and other emerging infectious diseases. Front. Microbiol. 7:68. doi: 10.3389/fmicb.2016.00068
Reed, L. J., and Muench, H. (1938). A simple method of estimating fifty per cent endpoints. Am. J. Epidemiol. 27, 493–497. doi: 10.1093/oxfordjournals.aje.a118408
Richards, S. A. (2008). Dealing with overdispersed count data in applied ecology. J. Appl. Ecol. 45, 218–227. doi: 10.1111/j.1365-2664.2007.01377.x
Rijks, J. M., Saucedo, B., Sluijs, A. S., der, Wilkie, G. S., Asten, A. J. A. M., et al. (2016). Investigation of amphibian mortality events in wildlife reveals an on-going Ranavirus epidemic in the north of the Netherlands. PLoS One 11:e0157473. doi: 10.1371/journal.pone.0157473
Robert, J., George, E., Andino, F. D. J., and Chen, G. (2011). Waterborne infectivity of the Ranavirus frog virus 3 in Xenopus laevis. Virology 417, 410–417. doi: 10.1016/j.virol.2011.06.026
Rosa, G. M., Sabino-Pinto, J., Laurentino, T. G., Martel, A., Pasmans, F., Rebelo, R., et al. (2017). Impact of asynchronous emergence of two lethal pathogens on amphibian assemblages. Sci. Rep. 7:43260. doi: 10.1038/srep43260
Ross, A. A., Hoffmann, A. R., and Neufeld, J. D. (2019). The skin microbiome of vertebrates. Microbiome 7:79. doi: 10.1186/s40168-019-0694-6
Saucedo, B., Garner, T. W., Kruithof, N., Allain, S. J., Goodman, M. J., Cranfield, R. J., et al. (2019). Common midwife toad RANAVIRUSes replicate first in the oral cavity of smooth newts (Lissotriton vulgaris) and show distinct strain-associated pathogenicity. Sci. Rep. 9:4453. doi: 10.1038/s41598-019-41214-0
Schirmer, M., Franzosa, E. A., Lloyd-Price, J., McIver, L. J., Schwager, R., Poon, T. W., et al. (2018). Dynamics of metatranscription in the inflammatory bowel disease gut microbiome. Nat. Microbiol. 3:337. doi: 10.1038/s41564-017-0089-z
Schloegel, L. M., Picco, A. M., Kilpatrick, A. M., Davies, A. J., Hyatt, A. D., and Daszak, P. (2009). Magnitude of the US trade in amphibians and presence of Batrachochytrium dendrobatidis and Ranavirus infection in imported North American bullfrogs (Rana catesbeiana). Biol. Conserv. 142, 1420–1426. doi: 10.1016/j.biocon.2009.02.007
Shade, A. (2017). Diversity is the question, not the answer. ISME J. 11, 1–6. doi: 10.1038/ismej.2016.118
Smith, A. H., Łukasik, P., O’Connor, M. P., Lee, A., Mayo, G., Drott, M. T., et al. (2015). Patterns, causes and consequences of defensive microbiome dynamics across multiple scales. Mol. Ecol. 24, 1135–1149. doi: 10.1111/mec.13095
Stark, T., Laurijssens, C., Weterings, M., Sluijs, A. S., der Martel, A., and Pasmans, F. (2014). Death in the clouds: Ranavirus associated mortality in assemblage of cloud forest amphibians in Nicaragua. Acta Herpetol. 9, 125–127.
Teacher, A. G. F., Cunningham, A. A., and Garner, T. W. J. (2010). Assessing the long-term impact of Ranavirus infection in wild common frog populations. Anim. Conserv. 13, 514–522. doi: 10.1111/j.1469-1795.2010.00373.x
Theriot, C. M., Koenigsknecht, M. J., Carlson, P. E. Jr., Hatton, G. E., Nelson, A. M., Li, B., et al. (2014). Antibiotic-induced shifts in the mouse gut microbiome and metabolome increase susceptibility to Clostridium difficile infection. Nat. Commun. 5:3114. doi: 10.1038/ncomms4114
Une, Y., Sakuma, A., Matsueda, H., Nakai, K., and Murakami, M. (2009). Ranavirus outbreak in North American bullfrogs (Rana catesbeiana), Japan, 2008. Emerg. Infect. Dis. 15, 1146–1147. doi: 10.3201/eid1507.081636
Videvall, E., Song, S. J., Bensch, H. M., Strandh, M., Engelbrecht, A., Serfontein, N., et al. (2019). Major shifts in gut microbiota during development and its relationship to growth in ostriches. Mol. Ecol. 28, 2653–2667. doi: 10.1111/mec.15087
Villarino, N. F., LeCleir, G. R., Denny, J. E., Dearth, S. P., Harding, C. L., Sloan, S. S., et al. (2016). Composition of the gut microbiota modulates the severity of malaria. Proc. Natl. Acad. Sci. 113, 2235–2240. doi: 10.1073/pnas.1504887113
Walke, J. B., Becker, M. H., Loftus, S. C., House, L. L., Teotonio, T. L., Minbiole, K. P., et al. (2015). Community structure and function of amphibian skin microbes: an experiment with bullfrogs exposed to a chytrid fungus. PLoS One 10:e0139848. doi: 10.1371/journal.pone.0139848
Warne, R. W., Kirschman, L., and Zeglin, L. (2019). Manipulation of gut microbiota during critical developmental windows affect host physiological performance and disease susceptibility across ontogeny. J. Anim. Ecol. 88, 845–856.
Whittington, R. J., Becker, J. A., and Dennis, M. M. (2010). Iridovirus infections in finfish – critical review with emphasis on Ranaviruses. J. Fish Dis. 33, 95–122. doi: 10.1111/j.1365-2761.2009.01110.x
Keywords: amphibian conservation, microbiome stability, host-microbe interactions, amphibian disease, ranavirus, FV3-like ranavirus
Citation: Harrison XA, Price SJ, Hopkins K, Leung WTM, Sergeant C and Garner TWJ (2019) Diversity-Stability Dynamics of the Amphibian Skin Microbiome and Susceptibility to a Lethal Viral Pathogen. Front. Microbiol. 10:2883. doi: 10.3389/fmicb.2019.02883
Received: 19 August 2019; Accepted: 29 November 2019;
Published: 20 December 2019.
Edited by:
Eria Alaide Rebollar, National Autonomous University of Mexico (UNAM), MexicoReviewed by:
Andrea J. Jani, University of Hawai‘i at Mânoa, United StatesC. Guilherme Becker, The University of Alabama, United States
Copyright © 2019 Harrison, Price, Hopkins, Leung, Sergeant and Garner. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xavier A. Harrison, x.harrison@exeter.ac.uk; xav.harrison@gmail.com
†ORCID: Xavier A. Harrison orcid.org/0000-0002-2004-3601