Skip to main content

Habitat disturbance influences the skin microbiome of a rediscovered neotropical-montane frog

Abstract

Background

The skin microbiome serves as a first line defense against pathogens in vertebrates. In amphibians, it has the potential to protect against the chytrid fungus Batrachochytrium dendrobatis (Bd), a likely agent of amphibian declines. Alteration of the microbiome associated with unfavorable environmental changes produced by anthropogenic activities may make the host more susceptible to pathogens. Some amphibian species that were thought to be “extinct” have been rediscovered years after population declines in the late 1980s probably due to evolved Bd-resistance and are now threatened by anthropogenic land-use changes. Understanding the effects of habitat disturbance on the host skin microbiome is relevant for understanding the health of these species, along with its susceptibility to pathogens such as Bd. Here, we investigate the influence of habitat alteration on the skin bacterial communities as well as specifically the putative Bd-inhibitory bacterial communities of the montane frog Lithobates vibicarius. This species, after years of not being observed, was rediscovered in small populations inhabiting undisturbed and disturbed landscapes, and with continuous presence of Bd.

Results

We found that cutaneous bacterial communities of tadpoles and adults differed between undisturbed and disturbed habitats. The adults from disturbed habitats exhibited greater community dispersion than those from undisturbed habitats. We observed a higher richness of putative Bd-inhibitory bacterial strains in adults from disturbed habitats than in those from undisturbed habitats, as well as a greater number of these potential protective bacteria with a high relative abundance.

Conclusions

Our findings support the microbial “Anna Karenina principle”, in which disturbance is hypothesized to cause greater microbial dispersion in communities, a so-called dysbiosis, which is a response of animal microbiomes to stress factors that decrease the ability of the host or its microbiome to regulate community composition. On the positive side, the high richness and relative abundance of putative Bd-inhibitory bacteria may indicate the development of a defense mechanism that enhances Bd-protection, attributed to a co-occurrence of more than 30-years of host and pathogen in these disturbed habitats. Our results provide important insight into the influence of human-modified landscapes on the skin microbiome and health implications of Bd-survivor species.

Background

All vertebrates are hosts to numerous symbiotic microorganisms that reside on their skin [1, 2]. The host depends on its natural microbial communities for different vital functions and can therefore contribute to maintaining host health (e.g. protection against pathogen invasion and colonization) [3,4,5]. Disruption of microbial communities beyond the natural range (i.e., dysbiosis) associated with environmental stress, such as unfavorable environmental changes caused by anthropogenic activities, may lead to an increased incidence of certain diseases in their hosts [6, 7]. Severe disruption of microbial communities caused by changes in the pattern of microbiota abundance and/or loss of beneficial bacterial members in the community can affect the functionality of the microbiome and make the host more susceptible to pathogens [8, 9]. Thus, dysbiosis caused by adverse factors can disrupt the performance benefits of host-associated microbiota and may contribute to unwanted collateral damage to hosts [7]. This means that it can threaten the survival of individuals and the growth of populations due to a cascade of negative health effects [10,11,12,13]. In fact, several recent studies have documented how skin and gut microbiomes from different wildlife species can be sensitive to habitat fragmentation and contaminants attributed to human activities [14,15,16,17,18].

In amphibians, the environment is a strong determinant of the skin microbiome, but less so than the evolutionary history and development [19, 20]. Different environmental conditions, together with the available pool of potential bacterial symbionts in the environment, have been shown to influence the host microbiome, which may result in distinct host bacterial communities across locations [21, 22]. Large-scale differences attributed to anthropogenic disturbance in the host environment (e.g. conversion of natural land use to agriculture, and habitat fragmentation) may alter abiotic factors such as soil and water properties, which in turn impact on the diversity of environmental microorganisms, including those that may colonize hosts [22]. Furthermore, as human activities also affect the physiology, ecology, and behavior of wildlife [23], this can cause simultaneous changes in the dynamics of the amphibian skin microbiome. These types of adverse conditions could cause a dysbiosis in the host microbiome that could have a negative impact on amphibian health [24]. Analysis of the effects of land use change (e.g., habitat disturbance by human activities) is an approach to understanding how microbial communities associated with hosts respond to adverse environmental conditions (e.g., Amato et al. [11] and Hughey et al. [25]).

Some neotropical amphibians that were thought to be extinct have been rediscovered years after the amphibian populations declined in the late 1980s [26,27,28]. These species have been found in isolated, small populations that persist in the presence of the fungal chytrid pathogen Batrachochytrium dendrobatidis (Bd) [28]. This pathogen is a likely causative agent of the amphibian declines [29, 30]. It is suggested that host recoveries after declines are not caused by attenuation of Bd and may be due to evolutionary changes in host response [31]. These hosts and Bd appear to be in an enzootic state (i.e. a period in which the host and pathogen coexist) that allows amphibians to recover, but they remain vulnerable to this pathogen as Bd remains as pathogenic as it was during catastrophic declines [28, 31]. Previous studies have highlighted the relevance of the host skin microbiome to the defense mechanisms of the amphibian skin against Bd [32, 33]. Therefore, understanding the impact of habitat disturbance on the host skin microbiome is relevant for understanding host resistance and susceptibility to pathogens such as Bd [24]. The latter is of great interest for the conservation of amphibians threatened by Bd that survive after their presumed extinction and are under the constant influence of human activities.

Our study focused on the Green-eyed frog (Lithobates vibicarius), which is a tropical mountain amphibian. This species declined dramatically throughout its historical range in the late 1980s, to the point of being considered possibly extinct, and probably due to an epizootic phase of Bd [28]. Surprisingly, years after its disappearance, some populations of the species have been rediscovered in the Cordillera Central and Cordillera de Tilarán in Costa Rica [28, 34]. The species now appears to persist with Bd in an enzootic phase [28]. In our study area, the species survives in small populations inhabiting undisturbed and disturbed landscapes, and with continuous presence of Bd (Alvarado et al. unpublished observations). Individuals in disturbed areas are exposed to human activities such as livestock and agriculture, which could potentially threaten the health of the species. We previously found that these surviving populations have a low prevalence of Bd-infection and a dynamic skin microbiome with putative Bd-inhibitory bacteria in our study area (see Jiménez et al. [35])

In this study, we examined the potential effects of anthropogenic disturbance on skin bacterial communities and putative Bd-inhibitory bacteria in L. vibicarius. We predicted that the composition of the skin bacterial community would differ between individuals in different habitat types (i.e. undisturbed and disturbed habitats). We also expected that individuals from disturbed habitats would possess a more prominent dispersion of the bacterial community, suggesting a dysbiotic state, as these types of stressors can impact the natural range of the microbiome through different mechanisms [6, 9]. Since L. vibicarius now appears to be able to persist despite the presence of Bd in our study area (see Jiménez et al. [35]) and amphibians harbor diverse bacterial communities that may provide protection against Bd [36, 37], we expected individuals in both habitat types to possess bacterial strains with putative Bd-inhibitory activity. However, as Bd infections in other amphibian species have been detected more frequently in areas with human disturbances [38, 39], we expected individuals from undisturbed and disturbed habitats to possess different distribution patterns of bacterial strains with putative Bd-inhibitory activity. In addition, we predicted a higher prevalence and relative abundance of putative Bd-inhibitory bacteria in individuals from disturbed habitats where Bd is suggested to be more present in amphibians and could be a developed potential mechanism to fight this pathogen [40, 41].

Results

The skin microbiome of frog samples from undisturbed and disturbed habitats were used for investigating the influence of habitat disturbance on the skin bacterial communities in tadpoles and adults of L. vibicarius. We prepared two microbial 16S rRNA datasets according to the analysis, i.e., 1) tadpole dataset and 2) adult dataset (Additional file 1: Table S1). The number of sequences per sample in the “tadpole dataset” ranged from 10,845 to 56,159 and in the “adult dataset” from 9378 to 54,704. The tadpole dataset contained samples collected from five sites (undisturbed sites = 2, disturbed sites = 3) in 2016 and 2017, while adult dataset contained samples collected at six study sites (undisturbed sites = 3, disturbed sites = 3) in 2016 and 2017 (Additional file 1: Table S1).

Skin bacterial alpha and beta diversity in tadpoles and adults of Lithobates vibicarius between undisturbed and disturbed habitats

Habitat type explained significant variation in the number of observed ASVs (χ2 = 6.54, p = 0.01) and the Shannon diversity (χ2 = 6.83, p = 0.008) in tadpoles (Additional file 2: Table S2). Both measures of alpha diversity were higher in undisturbed habitats than in disturbed habitats (number of observed ASVs: odds ratio (OR) = 1.36, 95% confidence intervals (95% CI) not overlapping 1, Shannon diversity: OR = 1.16, 95% CI not overlapping 1, Additional file 2: Fig. S1A-B). The phylogenetic diversity of the tadpoles was similar between undisturbed and disturbed habitats (χ2 = 1.08, p = 0.29, Additional file 2: Table S2, Additional file 2: Fig. S1C). We observed that year of survey influenced the number of observed ASVs and phylogenetic diversity in tadpoles (Additional file 2: Table S2).

We found no differences in alpha diversity of adults between undisturbed and disturbed habitats (number of observed ASVs: χ2 = 0.07, p = 0.79, Shannon diversity: χ2 = 1.61, p = 0.20, phylogenetic diversity, χ2 = 0.51, p = 0.47, Additional file 2: Table S2, Additional file 2: Fig. S1D-F). The three alpha diversity metrics in adults were not influenced by year of the survey (Additional file 2: Table S2).

Regarding beta diversity, we observed differences in the ASV presence/absence and ASV abundance pattern when accounting phylogenetic lineages (unweighted UniFrac: R2 = 0.07, p = 0.001, weighted UniFrac: R2 = 0.14, p = 0.001, Additional file 2: Table S3, Fig. 1a-b) in tadpoles between undisturbed and disturbed habitats. In addition, non-phylogenetic abundance-based pattern was distinct in tadpoles between the two habitat types (Bray-Curtis dissimilarity: R2 = 0.14, p = 0.001, Additional file 2: Table S3, Fig. 1c). The year of the survey explained some variation for the three beta diversity metrics (Additional file 2: Table S3). Additionally, we observed that study sites explained part of the differentiation on ASVs abundance pattern of tadpoles (Bray-Curtis dissimilarity: R2 = 0.36, p = 0.001, Additional file 2: Table S4, Additional file 2: Fig. S2A). In adults, the skin bacteria presence/absence and abundance pattern when accounting for phylogenetic lineages differed between undisturbed and disturbed habitats (unweighted UniFrac: R2 = 0.08, p = 0.004, weighted UniFrac: R2 = 0.12, p = 0.049, Additional file 2: Table S3, Fig. 1d-e). Further, we observed a difference based on abundance pattern between habitat types when phylogenetic lineages were not considered (Bray-Curtis dissimilarity: R2 = 0.07, p = 0.001, Additional file 2: Table S3, Fig. 1f). The year of the survey explained part of the variation for the three beta diversity metrics (Additional file 2: Table S3). The NMDS plot of ASVs abundance pattern showed that adults clustered by study sites, and adults in sites without habitat disturbance clustered more tightly than those in sites with habitat disturbance (Additional file 2: Fig. S2B).

Fig. 1
figure1

Beta diversity of L. vibicarius skin microbiota between undisturbed/disturbed habitats. The Non-Metric Multidimensional Scaling plots (NMDS) of the beta diversity of the microbiota of tadpoles (a-c) and adults (d-f) were based on unweighted UniFrac, weighted UniFrac and Bray-Curtis distance matrices. Each light point represents the bacterial community of an individual; point color indicates habitat types (green - undisturbed habitats and orange - disturbed habitats) and shape indicates year of survey (square - year 2016 and triangle - year 2017)

Tadpoles across both habitat types did not show significantly different dispersion of the skin bacterial community, according to the unweighted UniFrac distance measure (χ2 = 0.04, p = 0.84, Fig. 2a) and weighted UniFrac distance measure (χ2 = 0.31, p = 0.57, Fig. 2b), as opposed to the Bray-Curtis dissimilarity, which showed a marginally different dispersion (χ2 = 3.41, p = 0.05, Fig. 2c). The skin microbiome of tadpoles from the disturbed habitats showed a more dispersed community among replicates than those from the undisturbed habitats (Bray-Curtis dissimilarity: OR = 1.12, 95% CI = 1.00–1.27). Adults from undisturbed and disturbed habitats showed a clear difference in the community dispersion of bacterial skin consortia, based on the unweighted UniFrac distance measure (χ2 = 4.08, p = 0.04), the weighted UniFrac distance measure (χ2 = 19.94, p < 0.0001) and the Bray-Curtis dissimilarity (χ2 = 12.49, p = 0.001). Specifically, adults from disturbed habitats showed a greater community dispersion than those from undisturbed habitats (Fig. 2d-f). We found a stronger difference of community dispersion between habitat types with the weighted UniFrac distance measure (OR = 1.67, 95% CI = 1.30–2.14) and Bray-Curtis dissimilarity (OR = 1.29, 95% CI = 1.11–1.49) than with the unweighted UniFrac distance measure (OR = 1.10, 95% CI = 1.00–1.21). The latter suggests a stronger effect of habitat disturbance on community dispersion for the ASV abundance pattern than for the ASV presence/absence.

Fig. 2
figure2

Beta diversity dispersion of L. vibicarius skin microbiota between undisturbed/disturbed habitats. Inter-individual distances from tadpoles (ac) and adults (df) was calculated based on unweighted UniFrac, weighted UniFrac and Bray-Curtis distance matrices. Each light point represents the bacterial community of an individual; point color indicates habitat types (green - undisturbed habitats and orange - disturbed habitats)

Differential abundance of bacterial ASVs in the skin of tadpoles and adults of Lithobates vibicarius across habitat types

EdgeR analysis detected 75 and 106 ASVs in the skin of tadpoles and adults, respectively, which differed significantly in abundance between undisturbed and disturbed habitats (Additional file 3: Table S5 and Table S6). In tadpoles, we observed an increase in abundance of 39 ASVs in the undisturbed habitats and 36 ASVs in the disturbed habitats. Some ASVs with higher abundance in tadpoles in the undisturbed habitats belonged to the order Burkholderiales, the families Planococcaceae and Victivallaceae, and to genera PW3 and Petrovella (Fig. 3, Additional file 3: Table S5). On the contrary, the ASVs that show a significantly higher abundance in tadpoles in the disturbed habitats belong to several families such as Comamonadaceae and Ruminococcaceae, and to the genera Limnohabitans, Geobacter and Rhodospirillum (Fig. 3, Additional file 3: Table S5).

Fig. 3
figure3

Fold-changes of bacterial ASVs detected in tadpoles that differ between undisturbed/disturbed habitats. Bacterial ASV with a log2 fold change less than 0 were more abundant in the undisturbed habitat indicated on the left and colored on green, whereas those with a log2 fold change higher than 0 were more abundant in disturbed habitat indicated on the right and colored on orange. Bacterial ASVs (circles) are sized by mean relative abundance across samples. Only the ASVs identified at genus and species level are shown. See Additional file 3: Table S5 for complete results of edgeR analysis and detailed taxonomic affiliation

In adults, 39 ASVs were more abundant in undisturbed habitats and 67 ASVs in disturbed habitats. Some of the ASVs with significantly higher abundance in adults in the undisturbed habitats belonged to the families Ruminococcaceae and Nostocaceae, and to the genera Streptococcus and Actinomyces (Fig. 4, Additional file 3: Table S6). On the other hand, the ASVs showing significantly higher abundance in the disturbed habitats belonged to the families Lachnospiraceae, Bacteroidaceae and Comamonadaceae, and to the genera Bacteroides, Acinetobacter, Odoribacter and Pseudomonas (Fig. 4, Additional file 3: Table S6).

Fig. 4
figure4

Fold-changes of bacterial ASVs detected in adults that differ between undisturbed/disturbed habitats. Bacterial ASV with a log2 fold change less than 0 were more abundant in the undisturbed habitats indicated on the left and colored on green, whereas those with a log2 fold change higher than 0 were more abundant in disturbed habitats indicated on the right and colored on orange. Bacterial ASVs (circles) are sized by mean relative abundance across samples. Only the ASVs identified at genus and species level are shown. See Additional file 3: Table S6 for complete results of edgeR analysis and detailed taxonomic affiliation

Putative Bd-inhibitory bacterial community of adults of Lithobates vibicarius between undisturbed and disturbed habitats

The putative Bd-inhibitory bacterial community followed a similar pattern in alpha diversity as the entire bacterial community across habitat types, except for the number of observed ASVs with potential Bd inhibitory function, which showed a difference between undisturbed and disturbed habitats (Fig. 5a-c, Additional file 4: Table S7). Specifically, adults from the disturbed habitats showed a higher number of observed ASVs with potential Bd inhibitory function than those from undisturbed habitats (OR = 1.33, 95% CI = 1.02–1.75, Fig. 5a). The year of survey explained part of the variation for the observed ASVs and the Shannon diversity (Additional file 4: Table S7). Regarding beta diversity, the putative Bd-inhibitory bacterial community between habitat types differed when unweighted UniFrac (R2 = 0.02, p = 0.001) and Bray-Curtis dissimilarity (R2 = 0.06, p = 0.001) were used, while it remained similar when weighted UniFrac was used (R2 = 0.04, p = 0.09) (Fig. 5d-f, Additional file 4: Table S8). The year of survey explained part of the variation for the unweighted UniFrac and Bray-Curtis dissimilarity (Additional file 4: Table S8). The NMDS plot showed that adults in sites without habitat disturbance clustered more closely than those in sites with habitat disturbance (Additional file 4: Fig. S4).

Fig. 5
figure5

Diversity of the putative Bd-inhibitory bacterial community of L. vibicarius between undisturbed/disturbed habitats. Putative Bd-inhibitory bacterial alpha diversity of adults between habitat types (a-c). Lines indicate 95% confidence intervals (95% CIs). The Non-Metric Multidimensional Scaling plots (NMDS) of putative Bd-inhibitory bacterial beta diversity of adults regarding habitat types (d-f). Each light point represents the Bd-inhibitory bacterial ASV community of an individual; point color indicates habitat types (green - undisturbed habitats and orange - disturbed habitats) and shape indicates year of survey (square - year 2016 and triangle – year 2017). Ellipses show 95% CIs for each habitat type

Noteworthy, the edgeR analysis revealed important differences in the abundance pattern of putative Bd-inhibitory ASVs between adults in undisturbed and disturbed habitats (Fig. 6). We detected 18 putative Bd-inhibitory ASVs with significantly higher abundance in adults in disturbed habitats, whereas only five showed this difference in adults in undisturbed habitats (Fig. 6). Part of these ASVs of adults in the disturbed habitats belong to the genera Acinetobacter, Pseudomonas and Arthrobacter, and to the species Acinetobacter guillouiae, Pseudomonas umsongensis, P. viridiflava.

Fig. 6
figure6

Fold-changes of putative Bd-inhibitory bacterial ASVs detected in adults that differ between undisturbed/disturbed habitats. Bacterial ASV with a log2 fold change less than 0 were more abundant in the undisturbed habitats indicated on the left and colored on green, whereas those with a log2 fold change higher than 0 were more abundant in disturbed habitats indicated on the right and colored on orange. Bacterial ASVs are sized by mean relative abundance across samples. The highest possible taxonomic assignment (maximal to the species level) is shown for each ASV. G, Genus; F, Family; O, Order

Discussion

Using high-throughput amplicon sequencing, we provide evidence that habitat disturbance, which involves cattle grazing lands, agriculture activities, and loss of natural forest coverage can lead to changes in the natural range of variation of the amphibian skin microbiome, in an endangered species with suspected evolution of Bd resistance. Due to the involvement of the microbiome in host health, understanding how anthropogenic disturbance modifies host microbiomes and the extent to which this can impact the host health, provide us with relevant information to create healthier environments for amphibians [35].

We hypothesized that habitat disturbance would influence the skin microbiome of L. vibicarius. Our results show that the alpha and beta diversity of the tadpole and adult skin bacterial communities differed between undisturbed and disturbed habitats. In addition, the relative abundance of certain bacterial taxa on the skin of both life stages differed between the two habitat types. These microbiome changes are congruent with other studies with different taxa assessing the effects of environmental changes (e.g., the gut microbiota of Darwin’s finches in the Galapagos islands [42]) [11, 18, 42]. It is known that human activities have an impact on the abiotic conditions of the environment affecting the diversity of host bacterial communities. For example, water and soil characteristics have shown to influence the amphibian skin bacterial communities across distinct habitat types [37, 43]. Further, environmental conditions cause changes to the environmental microorganisms (e.g., soil and water bacteria) [22], including those that can colonize hosts. Therefore, the environmental conditions in our undisturbed and disturbed habitats may be directly and/or indirectly driving the diversity in the skin microbiome of L. vibicarius. Human activities can also increase host stress, affecting corticosterone levels that can induce changes in the composition of the microbial communities [23]. These patterns could be responsible for generating the observed differences in L. vibicarius skin microbiomes between undisturbed and disturbed habitats. Our results support the idea that habitat disturbance influences the skin microbiome structure of an endangered neotropical-montane frog, but due to low sample size (e.g., number of study sites), particularly for tadpoles, and the correlative nature of our data, our results need further elaboration. The role of abiotic factors and environmental bacterial symbionts in the amphibian skin microbiota variation are relevant for consideration in future work towards a better understanding of environmental changes in the skin microbiome variation of L. vibicarius.

We found that the number of observed ASVs and Shannon diversity in individuals of tadpoles from the disturbed habitats was lower than in those from the undisturbed habitats. Further, our data revealed significant differences of the beta diversity of the skin microbial communities in tadpoles and adults between habitat types. The low number of bacterial strains and reduced diversity of the microbiome may contribute to an increased susceptibility to opportunistic pathogens and poor health [41, 44, 45]. According to Bletz et al. [46], differences in the overall diversity of the microbial community also reflected differences in functional performance that may reflect variation in both presence and strength of host-microbe interactions (e.g., production of antimicrobial metabolites to defend the host from pathogens). If this is the case with L. vibicarius in this study, further detriment to habitat quality and the introduction of novel microbiota by livestock (e.g. potential pathogens) could cause health detriment to the host.

We observed that individuals in the disturbed habitats, particularly adult frogs, showed a bacterial community with higher dispersion than those found in the undisturbed habitats. This pattern led us to explore the “Anna Karenina Principle” (AKP) for disturbed microbiomes [9], which aims at a better understanding of host resilience, as the conformation of the microbiome can influence hosts fitness and health. The pattern of increased dispersal of L. vibicarius in disturbed habitats is suggestive of the effects of microbial AKP. These are characterized by a greater variation in individual profiles that may be to an increased hosts stress due to the presence of environmental stressors in the disturbed habitats that alters a stable state of the microbiome in unpredictable ways that affect host health [7, 9]. Previous studies have also shown AKP effects on microbiomes of wildlife species (e.g., the gut microbiome of macaques in contact with humans [47]). We cannot definitively determine whether AKP effects occur without longitudinal sampling, however, our initial findings suggest these effects. Since the skin microbiota of tadpoles and adults in disturbed habitats appear to be in a dysbiotic state relative to those in undisturbed habitats, they could have an increased susceptibility to pathogenic infections that threaten the health of the host [41, 48]. This is particularly important because there is increasing evidence that skin and non-skin pathogens affect the skin microbiota of amphibians [45, 47, 49], and amphibians in pastures are infested by a greater richness and abundance of parasites [50]. Although these explanations are intriguing and very plausible, a semi-controlled field experiment would be necessary to corroborate this hypothesis.

It is known that pathogens can cause both exacerbated microbial variability/dispersion and increased susceptibility to pathogenic infections (e.g., Salmonella typhimurium in the mouse [51] and Bd in the bullfrog [52]. Infection by other pathogens beyond Bd (e.g., Aeromonas hydrophila, Saprolegnia spp.), and even co-infections (e.g., Bd and Ranaviruses [53]) may increase the level of cutaneous microbial dispersion in L. vibicarius in disturbed habitats, causing unpredictable host health impacts. Further research with multiple surveys across the different reproductive periods of L. vibicarius will be essential to study these interactions, which will increase considerable our understanding of disease dynamics under a new approach, the disease pyramid (host – pathogen – microbiome – environment [54]).

According to our predictions, we observed different patterns in the diversity of the putative Bd-inhibitory bacterial community of adults between the two habitat types. The adults in disturbed habitats showed a greater richness of putative Bd-inhibitory bacteria than those in undisturbed habitats. Also, adults in disturbed habitats had a higher number of putative Bd-inhibitory bacteria with higher relative abundance than those in undisturbed habitats. We cannot conclude that these patterns of the bacterial community are part of the defense mechanism against Bd. However, a previous study found a higher proportion of pathogen-resistant bacteria in frogs inhabiting farms suggesting a higher occurrence rate of harmful pathogens in the farmland habitat [16]. Therefore, it is also possible that L. vibicarius remains at risk from Bd in the disturbed habitats. Besides that, previous studies indicate that richer and more diverse communities show greater functionality in protecting the host against Bd [40, 41, 55]. Furthermore, it is known that the relative abundance patterns of different protective bacterial taxa influence the protective function against invading pathogens, and the bacterial diversity patterns can determine in part the production of some metabolites against Bd by these symbionts [32, 56]. Based on this evidence, we suggest that both the high richness and the increased number of Bd-inhibitory bacterial strains with high relative abundance can be enhancing protection against Bd in adults of L. vibicarius in disturbed habitats. We should also consider the possibility of a functional redundancy of these Bd-inhibitory bacterial communities between adults in the two habitat types. A previous study found that the composition of the Bd-inhibitory bacterial community, using only bacterial isolates, of the salamander Plethodon cinereus changed along an anthropogenic land-use gradient [57]. The findings of Barnes et al. [57] indicate that this type of pattern does not correspond to a change in Bd-inhibitory function, on the contrary it is suggestive of a functionally redundant characteristic along a gradient of habitat change. Our results highlight the need of future research to investigate the anti-Bd functionality of the microbiome of L. vibicarius using a transcriptomic approach.

The presence of adult individuals in the disturbed habitats with the highest richness of putative Bd-inhibitory bacteria may help to understand why anthropogenic habitat loss is negatively associated with Bd in amphibian populations, along with the hypothesis that suboptimal microclimate conditions for Bd in disturbed habitats is a potential mechanism for reducing Bd infections, as suggested by Becker and Zamudio [58]. However, this should be interpreted with caution, because this effect could only be applied within species that can so far tolerate anthropogenic activities, and the level of disturbance in our study sites may not have seriously affected the functionality of the microbiome to protect against Bd. If the level of anthropogenic environmental stressors increases and the environmental conditions of the habitat change dramatically, it could be that the functionality of the Bd-protective bacterial community is disrupted, along with the proper functioning of the host immune system. This could therefore cause the spread of hosts that develop chytridiomycosis and increase the risk of extinction.

Besides the influence of habitat disturbance on the skin bacterial community of tadpoles and adults of L. vibicarius, we found temporal variation in these symbionts (i.e., between years of survey). The effect of time playing a role shaping the microbial communities is consistent with other studies showing a variation in the bacterial communities across short and long time scales [59, 60]. The observed community variation over time is likely associated with fluctuations in external environmental conditions, such as rainfall, humidity and temperature [60, 61]. Here, we also observed temporal variation in the putative Bd-inhibitory community that could suggest a disruption of the protective function, but also that functional redundancy could be occurring. Extrinsic factors from the environment are known to determine the stability of the host microbiome over time, which may therefore have relevant consequences for the stability of host-microbiome relationships and microbiome-fitness correlations [45]. More research is needed to elucidate the relevance of temporal variation in the skin bacterial communities to better understand factors influencing host-associated skin bacterial communities in the context of habitat disturbance and climate change. For example, an increment of livestock grazing across time can increase the amount of cattle feces, which can enhance the environment for novel pathogenic microorganisms that amphibians can acquire into their microbial community. Due to potential increment of anthropogenic pressures in the disturbed habitats of our study sites across years, a long-term microbiome surveillance will be relevant for exploring the dynamics (year-to-year variation) of putative Bd-inhibitory ASVs, in terms of survival and proliferation, especially for amphibian conservation efforts.

Conclusions

Our study suggests that habitat alteration influences the skin microbiome of L. vibicarius. Our findings are in line with the Anna Karenina Principle, which is an important response of animal microbiomes to stress factors that decrease the ability of the host or its microbiome to regulate community composition. The observed pattern of microbiome dispersion may be due to the presence of environmental stressors apart from pathogens, which can perturb a stable state leaving animals more susceptible to pathogen infections even beyond Bd. Our results further suggest that despite signs of a dysbiotic state in the disturbed habitats, the species inhabiting these habitats seems to harbor a Bd-inhibitory community with characteristics that enhance protection against Bd. The latter could be explained by the development of a defense mechanism to improve Bd-protection, attributed to the co-occurrence of host and pathogen for more than 30 years in these habitats [62]. However, we suspect that if the level of anthropogenic stressors increases, the Bd-protective bacterial community could be altered along with its functionality, putting the health of the host at an increased risk again. Despite limitations of the study due to low sample sizes, our results emphasize the importance of careful managements of landscapes outside protected areas to ensure the health of the endangered amphibian populations, particularly because of the potential introduction of new Bd bacterial strains and novel pathogens. Finally, our study contributes to efforts to develop microbial indices and baselines that can be used to manage or manipulate the amphibian skin microbiome to better protect host health across localities with different human use.

Methods

Study sites and sampling

The study was conducted at six sites in the high-altitude cloud forest of Juan Castro Blanco National Park and adjacent private lands within the province of Alajuela, Costa Rica (10° 16′ 44.4“ N, 84° 19’ 58.9” W). The elevation of the study sites ranges from 1850 to 2100 m.a.s.l. We have not provided the exact coordinates of the study sites to discourage people from visiting them, because we hope to prevent the possible introduction of pathogens and illegal collecting.

The study sites were monitored for the presence of L. vibicarius for two consecutive years (2016–2017) during the rainy season (September–November). Three study sites (Lagunillas, Monjes and Congo) represent areas in disturbed landscapes adjacent to the National Park. Each of these sites contains a pond and small streams within the pastures; the pastures are managed habitats for livestock grazing and agricultural activities and are surrounded by secondary and fragmented mountain rainforest. In addition, the loss of natural forest cover, the presence of cow feces, antibiotic treatments to cows and the use of pesticides (e.g. herbicides) on these sites could adversely affect the quality of habitat for amphibians. The other three study sites (Pozo verde, Pozo seco and Tamara) are in the protected and largely undisturbed National Park. Each of these sites comprises a large pond surrounded by well-preserved late successional montane forests.

We captured tadpoles and adults of L. vibicarius and collected two skin swabs to analyze the skin bacterial community and the Bd-infection status, using previously published procedures [35]. We released the animals at the site of captured after swabbing. The swabs for bacterial analysis were stored in sterile 1.5 ml Eppendorf tubes with 300 μl of DNA/RNA Shield (Zymo Research). The swabs for Bd analysis were placed directly in dry 1.5 ml sterile Eppendorf tubes. All swabs were stored at − 20 °C until further processing. We obtained unequal sample sizes between the study sites for each life stage due to unfavorable conditions during fieldwork. Details on sample sizes per study site are summarized in Additional file 1: Table S1. We used only the samples that were negative for Bd. The presence of Bd DNA from swabs was assessed according to previously published procedures [53]. In a previous study, a subsample of the individuals was used to understand the variation of life history in the skin microbiome [35], whereas the focus of the present study was to investigate the influence of habitat disturbance on the skin microbiome.

Molecular methods and sequencing

The bacterial genomic DNA from swabs was extracted using the NucleoSpin Soil (MACHEREY-NAGEL) following the manufacturer’s instructions and 16S rRNA gene amplicon sequencing methods were performed as described previously [35]. Briefly, we amplified the hypervariable V4 region of the 16S rRNA gene (291 bp) with the primer pair 515 F (5-GTGCCAGCMGCCGCGGTAA-3) and 806 R (5-GGACTACHVGGGTWTCTAAT-3). After purification (NucleoMag NGS Beads, Macherey-Nagel) and quantification (QuantiT™ PicoGreen® kit, Invitrogen-Life Technologies) of barcoded samples, the pooled sample library was sequenced as paired-end run on Illumina MiSeq platform at the Institute of Evolutionary Ecology and Conservation Genomics, Ulm University.

Sequence data processing

The initial processing of the sequence reads was done using the open source software QIIME2 (version 2019.1) according to previously published procedures [35]. The sequences were assigned to amplicon sequence variants (ASVs) and assigned to the taxonomy using the Greengenes database (version 13_8) as reference (http://greengenes.lbl.gov). Any sequence classified as chloroplast, mitochondria, archaea, eukaryota and unclassified phylum was excluded. We imported our data into R environment (version 3.4.4) [63] for further processing of the sequences using the R package “phyloseq” [64]. Samples with less than 8000 sequences were excluded. We also eliminated any ASV with less than 10 sequences. For alpha and beta diversity analyses, we normalized the read counts between samples by rarefying the datasets according to the sample with the lowest number of reads among samples.

Statistical analysis

To assess alpha diversity, we calculated three different metrics for each sample using the R packages “phyloseq” and “picante”: skin microbial community richness (i.e., the number of observed ASVs), ASVs diversity (i.e., the Shannon diversity), and phylogenetic diversity (i.e., Faiths phylogenetic diversity). We first performed Generalized Linear Models (GLMs) with a Gaussian distribution to test whether the alpha diversity of the tadpole skin bacterial communities differs between undisturbed and disturbed habitats. We incorporated the year of survey as a covariate in the GLMs to account for annual variation. Then, to assess whether the alpha diversity of the adults differs between undisturbed and disturbed habitats, we performed Generalized Linear Mixed Models (GLMMs) with a Gaussian error structure. Again, we incorporated the year of the study as a covariate and implemented study site identity as a random factor to control for non-independence of samples from the same site. Alpha diversity measures were log-transformed prior to analysis to improve normality. We determined the significance of covariates using the Anova function with χ2 as the test statistic. We used the odds ratio (OR) and 95% confidence intervals (95% CI) [43] to quantify the deviation (i.e., the effect size) of the alpha diversity metrics across habitats. Odds ratios and 95% CI were estimated using the R package “emmeans” [65].

We calculated three separate beta diversity metrics: unweighted UniFrac (phylogenetic presence/absence-based differences), weighted UniFrac (phylogenetic abundance-based differences) and Bray-Curtis dissimilarity (non-phylogenetic abundance-based differences) to investigate the influence of habitat disturbance on beta diversity of tadpole and adult skin microbiomes. We excluded ASVs with 10% prevalence or less, as beta diversity analyses are sensitive to rare ASVs variation. These measurements were calculated within the R package “phyloseq”. We fitted Permutational Multivariate Analysis of Variance (PERMANOVA) using the adonis function of the R package “vegan” [66] to test whether the measure of beta diversity differs significantly between habitat types. In the tadpole and adult analyses, we incorporated the year of survey as a covariate to account its variation. In the adult analyses, we included study site identity in the strata argument of the adonis function to consider the nesting effect. Due to the lack of tadpole samples at some sites during the two years of study, we could not include study site identity in the strata argument for the PERMANOVA analyses with tadpole samples. However, we performed a PERMANOVA using the Bray-Curtis dissimilarity and included study sites and year of study as covariates. Moreover, we used Non-Metric Multidimensional Scaling (NMDS) plots to visualize beta-diversity patterns using the R package “phyloseq”.

We performed multivariate of homogeneity of group dispersion analysis to evaluate how habitat types influences dispersal in bacterial communities among individuals in the same group (inter-individual distances) using the three beta diversity metrics. For each metric, we used a GLM with a Gaussian error structure to test the differences in the point-to-center distances between the microbial communities on the skin of individuals from undisturbed and disturbed habitats. We used the ORs to estimate the effect size between the habitats as described above.

Empirical Analysis of Digital Gene Expression Data (hereafter edgeR analysis) within the R package “edgeR” [67] was used to identify bacterial ASVs in both tadpoles and adults that significantly differed in abundance between undisturbed and disturbed habitat (FDR-corrected p-values at p < 0.001). For these analyses, we used the unrarefied dataset and the method Trimmed Mean of M-values (TMM) for normalization of samples.

We examined the community of putative Bd-inhibitory bacterial ASVs among habitat types as we did previously for the entire adult skin bacterial community of adults. For alpha and beta diversity, we followed the same statistical framework as above. To explore the variation of these bacterial strains between undisturbed and disturbed habitats, we performed an edgeR analysis to identify ASVs with differential abundance between habitat types. To identify putative Bd-inhibitory bacterial ASVs, we filtered the ASVs tables to retain only those reads that matched 100% sequence identity to those from the database of cultivable anti-Bd bacteria “Antifungal Isolates Database” [68]. Note that this does not necessarily mean that these bacterial strains possess anti-Bd properties, but that they might be candidates [69].

Availability of data and materials

The 16S rRNA sequences used in this study are available in the Dryad Digital Repository https://doi.org/10.5061/dryad.n34035p

Abbreviations

ASV:

Amplicon sequence variant

Bd :

Batrachochytrium dendrobatidis

AKP:

Anna Karenina principle

GLM:

Generalized Linear Models

GLMM:

Generalized Linear Mixed Models

PERMANOVA:

Permutational Multivariate Analysis of Variance

TMM:

Trimmed Mean of M-values

NMDS:

Non-Metric Multidimensional Scaling

OR:

odds ratio

References

  1. 1.

    McFall-Ngai M, Hadfield MG, Bosch TCG, Carey HV, Domazet-Lošo T, Douglas AE, et al. Animals in a bacterial world, a new imperative for the life sciences. PNAS. 2013;110:3229–36.

    CAS  PubMed  Google Scholar 

  2. 2.

    Pita L, Rix L, Slaby BM, Franke A, Hentschel U. The sponge holobiont in a changing ocean: from microbes to ecosystems. Microbiome. 2018;6:46.

  3. 3.

    Ubeda C, Djukovic A, Isaac S. Roles of the intestinal microbiota in pathogen protection. Clin Transl Immunol. 2017;6:e128.

    Google Scholar 

  4. 4.

    Grogan LF, Cashins SD, Skerratt LF, Berger L, McFadden MS, Harlow P, et al. Evolution of resistance to chytridiomycosis is associated with a robust early immune response. Mol Ecol. 2018;27:919–34.

    PubMed  Google Scholar 

  5. 5.

    McKenney EA, Koelle K, Dunn RR, Yoder AD. The ecosystem services of animal microbiomes. Mol Ecol. 2018;27:2164–72.

    CAS  PubMed  Google Scholar 

  6. 6.

    Trevelline BK, Fontaine SS, Hartup BK, Kohl KD. Conservation biology needs a microbial renaissance: a call for the consideration of host-associated microbiota in wildlife management practices. P Roy Soc B-Biol Sci. 2019;286:20182448.

    Google Scholar 

  7. 7.

    West AG, Waite DW, Deines P, Bourne DG, Digby A, McKenzie VJ, et al. The microbiome in threatened species conservation. Biol Conserv. 2019;229:85–98.

    Google Scholar 

  8. 8.

    Petersen C, Round JL. Defining dysbiosis and its influence on host immunity and disease. Cell Microbiol. 2014;16:1024–33.

    CAS  PubMed  PubMed Central  Google Scholar 

  9. 9.

    Zaneveld JR, McMinds R, Vega TR. Stress and stability: applying the Anna Karenina principle to animal microbiomes. Nat Microbiol. 2017;2:17121.

    CAS  PubMed  Google Scholar 

  10. 10.

    Kinross JM, Darzi AW, Nicholson JK. Gut microbiome-host interactions in health and disease. Genome Med. 2011;3:14.

    PubMed  PubMed Central  Google Scholar 

  11. 11.

    Amato KR, Yeoman CJ, Kent A, Righini N, Carbonero F, Estrada A, et al. Habitat degradation impacts black howler monkey (Alouatta pigra) gastrointestinal microbiomes. ISME J. 2013;7:1344–53.

    CAS  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Althani AA, Marei HE, Hamdi WS, Nasrallah GK, El Zowalaty ME, Al Khodor S, et al. Human microbiome and its association with health and diseases. J Cell Physiol. 2016;231:1688–94.

    CAS  PubMed  Google Scholar 

  13. 13.

    Fackelmann G, Sommer S. Microplastics and the gut microbiome: how chronically exposed species may suffer from gut dysbiosis. Mar Pollut Bull. 2019;143:193–203.

    CAS  PubMed  Google Scholar 

  14. 14.

    Barelli C, Albanese D, Donati C, Pindo M, Dallago C, Rovero F, et al. Habitat fragmentation is associated to gut microbiota diversity of an endangered primate: implications for conservation. Sci Rep. 2015;5:14862.

    CAS  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Bennett G, Malone M, Sauther ML, Cuozzo FP, White B, Nelson KE, et al. Host age, social group, and habitat type influence the gut microbiota of wild ring-tailed lemurs (Lemur catta). Am J Primatol. 2016;78:883–92.

    CAS  PubMed  Google Scholar 

  16. 16.

    Chang C-W, Huang B-H, Lin S-M, Huang C-L, Liao P-C. Changes of diet and dominant intestinal microbes in farmland frogs. BMC Microbiol. 2016;16:33.

    PubMed  PubMed Central  Google Scholar 

  17. 17.

    Huang B-H, Chang C-W, Huang C-W, Gao J, Liao P-C. Composition and functional specialists of the gut microbiota of frogs reflect habitat differences and agricultural activity. Front Microbiol. 2018;8:2670.

  18. 18.

    Assis AB, Bevier CR, Barreto CC, Navas CA. Environmental influences on and antimicrobial activity of the skin microbiota of Proceratophrys boiei (Amphibia, Anura) across forest fragments. Ecol Evol. 2020;10:901–13.

    Google Scholar 

  19. 19.

    Kueneman JG, Parfrey LW, Woodhams DC, Archer HM, Knight R, McKenzie VJ. The amphibian skin-associated microbiome across species, space and life history stages. Mol Ecol. 2014;23:1238–50.

    PubMed  Google Scholar 

  20. 20.

    Rebollar EA, Hughey MC, Medina D, Harris RN, Ibáñez R, Belden LK. Skin bacterial diversity of Panamanian frogs is associated with host susceptibility and presence of Batrachochytrium dendrobatidis. ISME J. 2016;10:1682–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Muletz-Wolz CR, Myers JM, Domangue RJ, Herrick JB, Harris RN. Soil bioaugmentation with amphibian cutaneous bacteria protects amphibian hosts from infection by Batrachochytrium dendrobatidis. Biol Conserv. 2012;152:119–26.

    Google Scholar 

  22. 22.

    Wang H, Marshall CW, Cheng M, Xu H, Li H, Yang X, et al. Changes in land use driven by urbanization impact nitrogen cycling and the microbial community composition in soils. Sci Rep. 2017;7:1–12.

    PubMed  PubMed Central  Google Scholar 

  23. 23.

    Noguera JC, Aira M, Pérez-Losada M, Domínguez J, Velando A. Glucocorticoids modulate gastrointestinal microbiome in a wild bird. R Soc Open Sci. 2018;5:171743.

    PubMed  PubMed Central  Google Scholar 

  24. 24.

    Jiménez RR, Sommer S. The amphibian microbiome: natural range of variation, pathogenic dysbiosis, and role in conservation. Biodivers Conserv. 2017;26:763–86.

    Google Scholar 

  25. 25.

    Hughey MC, Pena JA, Reyes R, Medina D, Belden LK, Burrowes PA. Skin bacterial microbiome of a generalist Puerto Rican frog varies along elevation and land use gradients. PeerJ. 2017;5:e3688.

    PubMed  PubMed Central  Google Scholar 

  26. 26.

    González-Maya JF, Belant JL, Wyatt SA, Schipper J, Cardenal J, Corrales D, et al. Renewing hope: the rediscovery of Atelopus varius in Costa Rica. Amphib-Reptile. 2013;34:573–8.

    Google Scholar 

  27. 27.

    Jiménez RR, Alvarado G. Craugastor escoces (Anura: Craugastoridae) reappears after 30 years: rediscovery of an “extinct” neotropical frog. Amphib-Reptile. 2017;38:257–9.

    Google Scholar 

  28. 28.

    Whitfield SM, Alvarado G, Abarca J, Zumbado H, Zuñiga I, Wainwright M, et al. Differential patterns of Batrachochytrium dendrobatidis infection in relict amphibian populations following severe disease-associated declines. Dis Aquat Org. 2017;126:33–41.

    PubMed  Google Scholar 

  29. 29.

    Stuart SN, Chanson JS, Cox NA, Young BE, Rodrigues ASL, Fischman DL, et al. Status and trends of amphibian declines and extinctions worldwide. Science. 2004;306:1783–6.

    CAS  PubMed  Google Scholar 

  30. 30.

    Scheele BC, Pasmans F, Skerratt LF, Berger L, Martel A, Beukema W, et al. Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science. 2019;363:1459–63.

    CAS  PubMed  Google Scholar 

  31. 31.

    Voyles J, Woodhams DC, Saenz V, Byrne AQ, Perez R, Rios-Sotelo G, et al. Shifts in disease dynamics in a tropical amphibian assemblage are not due to pathogen attenuation. Science. 2018;359:1517–9.

    CAS  PubMed  Google Scholar 

  32. 32.

    Harris RN, Brucker RM, Walke JB, Becker MH, Schwantes CR, Flaherty DC, et al. Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 2009;3:818–24.

    CAS  PubMed  Google Scholar 

  33. 33.

    Muletz-Wolz CR, Almario JG, Barnett SE, DiRenzo GV, Martel A, Pasmans F, et al. Inhibition of fungal pathogens across genotypes and temperatures by amphibian skin bacteria. Front Microbiol. 2017;8:1551.

    PubMed  PubMed Central  Google Scholar 

  34. 34.

    Castro-Cruz A, García-Fernández F. Report on Lithobates vibicarius (Cope 1894) (Anura: Ranidae) in Parque Nacional de Agua Juan Castro Blanco, Alajuela. Costa Rica Froglog. 2012;100:69–70.

    Google Scholar 

  35. 35.

    Jiménez RR, Alvarado G, Estrella J, Sommer S. Moving beyond the host: unraveling the skin microbiome of endangered Costa Rican amphibians. Front Microbiol. 2019;10:2060.

    PubMed  PubMed Central  Google Scholar 

  36. 36.

    Muletz-Wolz CR, DiRenzo GV, Yarwood SA, Campbell Grant EH, Fleischer RC, Lips KR. Antifungal bacteria on woodland salamander skin exhibit high taxonomic diversity and geographic variability. Appl Environ Microbiol. 2017;83:e00186.

  37. 37.

    Muletz-Wolz CR, Yarwood SA, Grant EHC, Fleischer RC, Lips KR. Effects of host species and environment on the skin microbiome of plethodontid salamanders. J Anim Ecol. 2018;87:341–53.

    PubMed  Google Scholar 

  38. 38.

    Siddons Spencer R, Bray Marin C, Searle CL. Higher infection prevalence in amphibians inhabiting human-made compared to natural wetlands. J Wildl Dis. 2020;56.

  39. 39.

    Murray KA, Retallick RWR, Puschendorf R, Skerratt LF, Rosauer D, McCallum HI, et al. Assessing spatial patterns of disease risk to biodiversity: implications for the management of the amphibian pathogen, Batrachochytrium dendrobatidis. J App Ecol. 2011;48:163–73.

    Google Scholar 

  40. 40.

    Piovia-Scott J, Rejmanek D, Woodhams DC, Worth SJ, Kenny H, McKenzie V, et al. Greater species richness of bacterial skin symbionts better suppresses the amphibian fungal pathogen Batrachochytrium dendrobatidis. Microb Ecol. 2017;74:217–26.

    PubMed  Google Scholar 

  41. 41.

    Ellison S, Knapp RA, Sparagon W, Swei A, Vredenburg VT. Reduced skin bacterial diversity correlates with increased pathogen infection intensity in an endangered amphibian host. Mol Ecol. 2018;28:127–40.

    PubMed  Google Scholar 

  42. 42.

    Knutie SA, Chaves JA, Gotanda KM. Human activity can influence the gut microbiota of Darwin’s finches in the Galapagos Islands. Mol Ecol. 2019;28:2441–50.

    PubMed  Google Scholar 

  43. 43.

    Krynak KL, Burke DJ, Benard MF. Landscape and water characteristics correlate with immune defense traits across Blanchard’s cricket frog (Acris blanchardi) populations. Biol Conserv. 2016;193:153–67.

    Google Scholar 

  44. 44.

    Carroll IM, Ringel-Kulka T, Siddle JP, Ringel Y. Alterations in composition and diversity of the intestinal microbiota in patients with diarrhea-predominant irritable bowel syndrome. Neurogastroenterol Motil. 2012;24:521–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Harrison XA, Price SJ, Hopkins K, Leung WTM, Sergeant C, Garner TWJ. Diversity-stability dynamics of the amphibian skin microbiome and susceptibility to a lethal viral pathogen. Front Microbiol. 2019;10:2883.

  46. 46.

    Bletz MC, Goedbloed DJ, Sanchez E, Reinhardt T, Tebbe CC, Bhuju S, et al. Amphibian gut microbiota shifts differentially in community structure but converges on habitat-specific predicted functions. Nat Commun. 2016;7:1–12.

    Google Scholar 

  47. 47.

    Grant ET, Kyes RC, Kyes P, Trinh P, Ramirez V, Tanee T, et al. Fecal microbiota dysbiosis in macaques and humans within a shared environment. PLoS One. 2019;14:e0210679.

  48. 48.

    Lozupone CA, Stombaugh JI, Gordon JI, Jansson JK, Knight R. Diversity, stability and resilience of the human gut microbiota. Nature. 2012;489:220–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Federici E, Rossi R, Fidati L, Paracucchi R, Scargetta S, Montalbani E, et al. Characterization of the skin microbiota in Italian stream frogs (Rana italica) infected and uninfected by a cutaneous parasitic disease. Microbes Environ. 2015;30:262–9.

    PubMed  PubMed Central  Google Scholar 

  50. 50.

    McKenzie VJ. Human land use and patterns of parasitism in tropical amphibian hosts. Biol Conserv. 2007;137:102–16.

    Google Scholar 

  51. 51.

    Rausch P, Steck N, Suwandi A, Seidel JA, Künzel S, Bhullar K, et al. Expression of the blood-group-related gene B4galnt2 alters susceptibility to Salmonella infection. PLoS Pathog. 2015;11:e1005008.

    PubMed  PubMed Central  Google Scholar 

  52. 52.

    Walke JB, Becker MH, Loftus SC, House LL, Teotonio TL, Minbiole KPC, et al. Community structure and function of amphibian skin microbes: an experiment with bullfrogs exposed to a chytrid fungus. PLoS One. 2015;10:e0139848.

    PubMed  PubMed Central  Google Scholar 

  53. 53.

    Whitfield SM, Geerdes E, Chacon I, Ballestero Rodriguez E, Jimenez RR, Donnelly MA, et al. Infection and co-infection by the amphibian chytrid fungus and ranavirus in wild Costa Rican frogs. Dis Aquat Org. 2013;104:173–8.

    PubMed  Google Scholar 

  54. 54.

    Bernardo-Cravo AP, Schmeller DS, Chatzinotas A, Vredenburg VT, Loyau A. Environmental factors and host microbiomes shape host–pathogen dynamics. Trends Parasitol. 2020;36:616–33.

    CAS  PubMed  Google Scholar 

  55. 55.

    Bell SC, Garland S, Alford RA. Increased numbers of culturable inhibitory bacterial taxa may mitigate the effects of Batrachochytrium dendrobatidis in Australian wet tropics frogs. Front Microbiol. 2018;9:1604.

  56. 56.

    Harris RN, Lauer A, Simon MA, Banning JL, Alford RA. Addition of antifungal skin bacteria to salamanders ameliorates the effects of chytridiomycosis. Dis Aquat Org. 2009;83:11–6.

    PubMed  Google Scholar 

  57. 57.

    Barnes EM, Carter EL, Lewis JD. Predicting microbiome function across space is confounded by strain-level differences and functional redundancy across taxa. Front Microbiol. 2020;11:101.

  58. 58.

    Becker CG, Zamudio KR. Tropical amphibian populations experience higher disease risk in natural habitats. Proc Natl Acad Sci U S A. 2011;108:9893–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Sabino-Pinto J, Galán P, Rodríguez S, Bletz MC, Bhuju S, Geffers R, et al. Temporal changes in cutaneous bacterial communities of terrestrial- and aquatic-phase newts (Amphibia). Environ Microbiol. 2017;19:3025–38.

    CAS  PubMed  Google Scholar 

  60. 60.

    Estrada A, Hughey MC, Medina D, Rebollar EA, Walke JB, Harris RN, et al. Skin bacterial communities of neotropical treefrogs vary with local environmental conditions at the time of sampling. PeerJ. 2019;7:e7044.

  61. 61.

    Bletz MC, Perl RGB, Bobowski BT, Japke LM, Tebbe CC, Dohrmann AB, et al. Amphibian skin microbiota exhibits temporal variation in community structure but stability of predicted Bd-inhibitory function. ISME J. 2017;11:1521–34.

    PubMed  PubMed Central  Google Scholar 

  62. 62.

    De León ME, Zumbado-Ulate H, García-Rodríguez A, Alvarado G, Sulaeman H, Bolaños F, et al. Batrachochytrium dendrobatidis infection in amphibians predates first known epizootic in Costa Rica. PLoS One. 2019;14:e0208969.

    PubMed  PubMed Central  Google Scholar 

  63. 63.

    R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2019.

    Google Scholar 

  64. 64.

    McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. Plos One. 2013;8:e61217.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Lenth R. emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.4.7; 2019.

    Google Scholar 

  66. 66.

    Oksanen J, Guillaume Blanchet F, Kindt R, Legendre P, Minchin PR, O’Hara, R. B., et al. vegan: community ecology package. R package version 2.5-4. 2019.

    Google Scholar 

  67. 67.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    CAS  PubMed  Google Scholar 

  68. 68.

    Woodhams DC, Alford RA, Antwis RE, Archer H, Becker MH, Belden LK, et al. Antifungal isolates database of amphibian skin-associated bacteria and function against emerging fungal pathogens. Ecol. 2015;96:595.

    Google Scholar 

  69. 69.

    Becker MH, Walke JB, Murrill L, Woodhams DC, Reinert LK, Rollins-Smith LA, et al. Phylogenetic distribution of symbiotic bacteria from Panamanian amphibians that inhibit growth of the lethal fungal pathogen Batrachochytrium dendrobatidis. Mol Ecol. 2015;24:1628–41.

    PubMed  Google Scholar 

Download references

Acknowledgments

We are grateful to the APANAJUCA, the SINAC, Douglas Vargas, Erick Ballestero, Michael Méndez, Jeremy Quirós, Alexa Morales, Francisco Quesada, and Rolando Campos for supporting field logistics. We thank Stefan Brändel, Luis Víquez, Alexander Heni, and especially Kerstin Wilhelm for technical and laboratory assistance. We also extend our thanks to Mark Gillingham and Karina Montero for the bioinformatic support. We thank Steven Whitfield and Brian Folt for scientific input and language corrections during manuscript preparation. We thank the reviewers for their thoughtful comments as these have helped to improve this paper.

Funding

This study was supported by the grants from the Auckland Zoo, Rufford Foundation, Cleveland Metroparks Zoo, Cleveland Zoological Society, and start-up funding to SS by the University of Ulm. RRJ was supported by the National Council for Scientific and Technological Research of Costa Rica and the Ministry of Science and Technology of Costa Rica. Open Access funding enabled and organized by Projekt DEAL.

Author information

Affiliations

Authors

Contributions

RRJ, SS, and GA contributed to the research design. RRJ, GA, and JS performed the sample collection. RRJ carried out laboratory work on microbiome samples. RRJ performed bioinformatic analysis and statistical analysis. RRJ wrote the manuscript with input from SS, GA, and JS. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Randall R. Jiménez.

Ethics declarations

Ethics approval and consent to participate

Sampling protocols and the animal care protocol were reviewed and approved by National Commission for the Biodiversity Management of Costa Rica (CONAGEBio, R-055-2015-OT-CONAGEBio) and Costa Rica MINAE-SINAC (SINAC-ACAHN-PI-R-004-2015).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1 : Table S1

. Samples of Lithobates vibicarius used for microbiome analysis. Samples with > 8000 reads.

Additional file 2 : Table S2.

Summary of GLMs and GLMMs predicting alpha diversity metrics in tadpoles and adults, respectively. Habitat types and year of survey were use as covariates. Study site was used as a random factor in the GLMMs. Significant p values (< 0.05) are shown in bold. Fig. S1. Alpha diversity metrics among tadpoles (A-C) and adults (D-F) regarding habitat types. Each light point represents the bacterial skin community of an individual sample; point color indicates habitat type (green - undisturbed habitat and orange - disturbed habitat). Lines indicate 95% confidence intervals (95% CIs). Table S3. Summary of PERMANOVAs of beta diversity metrics in tadpoles and adults. Habitat types and year of survey were use as covariates and study site were used as covariates. Study site was used as a random factor in the analysis of adults. Significant p values (< 0.05) are shown in bold. Table S4. Summary of PERMANOVAs of beta diversity in tadpoles. Study types and year of survey were use as covariates. Significant p values (< 0.05) are shown in bold. Fig. S2. Beta diversity of L. vibicarius skin microbiota among study sites and year of survey. The Non-Metric Multidimensional Scaling plots (NMDS) of the beta diversity of the microbiota of tadpoles (A) and adults (B) were based on Bray-Curtis dissimilarity. Each point represents the bacterial community of an individual; point color indicates study sites and shape indicates year of survey.

Additional file 3 : Table S5.

Fold-changes of bacterial ASVs detected on the skin of tadpoles that differ between undisturbed/disturbed habitats. Table S6. Fold-changes of bacterial ASVs detected on the skin of adults that differ between undisturbed/disturbed habitats.

Additional file 4 : Table S7.

Summary of the GLMMs predicting putative Bd-inhibitory bacterial alpha diversity metrics in adults. Habitat types and year of survey were use as covariates and study site was used as a random factor in the models. Significant p values (< 0.05) are shown in bold. Table S8. Summary of PERMANOVAs of the putative Bd-inhibitory bacterial beta diversity in adults. Habitat types and year of survey were use as covariates and study site was used as a random factor in the models. Significant p values (< 0.05) are shown in bold. Fig. S4. Beta diversity of putative Bd-inhibitory bacterial community on adults across sites and year of survey. The Non-Metric Multidimensional Scaling plots (NMDS) of the beta diversity of the microbiota of adults were based on Bray-Curtis dissimilarity. Each point represents the bacterial community of an individual; point color indicates study sites and shape indicates year of survey.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Jiménez, R.R., Alvarado, G., Sandoval, J. et al. Habitat disturbance influences the skin microbiome of a rediscovered neotropical-montane frog. BMC Microbiol 20, 292 (2020). https://doi.org/10.1186/s12866-020-01979-1

Download citation

Keywords

  • Skin bacterial diversity
  • 16S rRNA amplicon sequencing
  • Land-use change
  • Anthropogenic disturbance
  • Batrachochytrium dendrobatidis (Bd)
  • Amphibian health implications and conservation