Skip to main content
Download PDF

The skin microbiota of the axolotl Ambystoma altamirani is highly influenced by metamorphosis and seasonality but not by pathogen infection

Animal Microbiome volume 4, Article number: 63 (2022) Cite this article

Abstract

Background

Microbiomes have been increasingly recognized as major contributors to host health and survival. In amphibians, bacterial members of the skin microbiota protect their hosts by inhibiting the growth of the fungal pathogen Batrachochytrium dendrobatidis (Bd). Even though several studies describe the influence of biotic and abiotic factors over the skin microbiota, it remains unclear how these symbiotic bacterial communities vary across time and development. This is particularly relevant for species that undergo metamorphosis as it has been shown that host physiology and ecology drastically influence diversity of the skin microbiome.

Results

We found that the skin bacterial communities of the axolotl A. altamirani are largely influenced by the metamorphic status of the host and by seasonal variation of abiotic factors such as temperature, pH, dissolved oxygen and conductivity. Despite high Bd prevalence in these samples, the bacterial diversity of the skin microbiota did not differ between infected and non-infected axolotls, although relative abundance of particular bacteria were correlated with Bd infection intensity.

Conclusions

Our work shows that metamorphosis is a crucial process that shapes skin bacterial communities and that axolotls under different developmental stages respond differently to environmental seasonal variations. Moreover, this study greatly contributes to a better understanding of the factors that shape amphibian skin microbiota, especially in a largely underexplored group like axolotls (Mexican Ambystoma species).

Background

Host associated microbiomes are vital for host health and survival, as they play relevant functions linked to nutrition, reproduction, behavior, defense against pathogens or predators [1,2,3,4,5]. Specifically, some animal associated microbiomes contribute to host health due to their ability to inhibit the growth of pathogens responsible for infectious diseases threatening diverse host species such as bats, snakes, or amphibians [6,7,8]. For instance, it has been shown that some members of the amphibian skin microbiome inhibit the growth of the lethal pathogens Batrachochytrium dendrobatidis (Bd) and B. salamandrivorans [9,10,11,12], which have caused amphibian populations declines and extinctions worldwide [13].

Studies accumulated over the past two decades showed that the amphibian skin microbiome is influenced by host associated factors (host genetics and development) [14,15,16], microhabitat related factors (environmental microorganisms, habitat abiotic conditions and pathogen presence) [17,18,19,20,21], and climatic and geographical factors (seasonality, precipitation, temperature or land use) [14, 22,23,24,25].

In the case of host-associated factors, it has been shown that the skin microbiota of amphibians (specifically frogs) changes across development and particularly before and after metamorphosis [26,27,28]. During metamorphosis amphibians in larval stages transition to adults following a series of physiological rearrangements such as tail reabsorption, limb development and remodeling of muscles, heart, intestine brain, and skin [29]. Metamorphosis also induces immunosuppression in response to thyroid and corticosteroid hormone signaling and eventually the immune system reorganizes and gradually matures in newly metamorphosed adults [30].

Along with physiological rearrangements, many amphibian species go through behavioral and lifestyle changes, while larval stages inhabit aquatic environments, adults become terrestrial and only return to water environments in the reproductive season [31,32,33]. These changes in microhabitat occupancy could influence skin microbial composition since the environmental microbial communities are one of the main sources of microbial diversity [16, 17].

In the case of climatic factors, temporal variation of abiotic factors [34] such as temperature and precipitation have a strong influence over amphibian skin microbial community structure [22, 35]. For example, in tropical regions microbial diversity on the amphibian skin differs between wet and dry seasons [19, 26, 36]. In temperate regions, where the four seasons are well defined through the year, seasonal changes have been linked to the temporal dynamics of the amphibian skin microbiota [22, 37,38,39]. In addition, it has been shown that spatial variation such as elevation gradients [40,41,42] or distinct microhabitats [43] influence the skin microbial diversity of amphibians.

Bd influence over the amphibian skin microbiota has been described in amphibian species with contrasting Bd infection status (infected–non-infected [19] and high Bd prevalence—low Bd prevalence [44]). These studies showed that disruption of skin microbiota following Bd infection can influence host survival and that the final outcome of the infection depends on the interplay between host, microbiome and the environment [21, 23, 45].

Here we analyzed the skin bacterial diversity of the axolotl Ambystoma altamirani, a stream dwelling salamander endemic to conifer and oak-pine forest from the central region of Mexico [46]. A. altamirani is a facultative paedomorphic species in which, metamorphic (without gills) and pre-metamorphic (with gills) individuals inhabit the same streams all year long [47, 48], allowing us to evaluate how metamorphosis and seasonality influence the skin microbiota in a species living in the same aquatic environment across time and development. In addition, we evaluated if skin microbiota differs from environmental bacterial communities and if Bd presence and infection intensity influence the skin microbiota of A. altamirani. We hypothesized that A. altamirani skin microbiota would (a) differ from environmental bacterial communities, (b) vary between metamorphic and pre-metamorphic salamanders, (c) change across seasons and (d) differ according to Bd infection status.

Results

We sampled a total of 279 A. altamirani individuals (85 metamorphic and 194 pre-metamorphic) at four locations across four seasons. Additionally, 159 environmental samples from sediment (80) and water (79) were collected. After quality control and rarefaction at 10,000 reads per sample, 13 samples were discarded and the remaining 438 samples were used to perform all diversity analyses (Table 1). A final table with a total of 72,408 amplicon sequence variants (ASVs) was obtained including all samples.

Table 1 Final list of collected samples that passed bioinformatic filters
Full size table

A. altamirani skin microbiota differs from environmental bacterial communities

When comparing the number of unique and shared ASVs across sample types, we found that each sample type harbored many unique ASVs and only 2408 ASVs (3.32% of the total) were shared among the four sample types (Fig. 1A). Sediment and water samples were the samples with highest numbers of unique ASVs (20,031 and 9902 respectively), while metamorphic and pre-metamorphic samples had 8916, and 6650 unique ASVs respectively. Interestingly only 677 ASVs (1.16% of the total ASVs) were shared between metamorphic and pre-metamorphic salamanders.

Fig. 1
figure 1

Bacterial diversity of A. altamirani skin and environmental samples. A Upset plot illustrating the number of unique and shared ASVs. Numbers aside the color bars indicate how many ASVs were present on each sample type (color bars) and shared between sample types (gray bars). B Alpha Faith’s Phylogenetic diversity (PD) across sample types. C Principal coordinate analysis (PCoA) based on weighted UniFrac distances across sample types. D Beta dispersion using Analysis of multivariate homogeneity of groups dispersions. Letters a–d indicate statistically significant comparisons

Full size image

Taxonomic results showed that, Burkholderiaceae was the most abundant bacterial family in all four sample types accounting for 32.6% and 51.1% of the relative abundance in metamorphic and pre-metamorphic samples respectively, and 14.6% and 40.8% of sediment and water respectively (Additional file 1: Figure S1). For the axolotl samples we found that Chitinophagaceae and Pseudomonadaceae varied in relative abundance according to host metamorphic status, with Chitinophagaceae showing a higher abundance in pre-metamorphic axolotls (metamorphic 2.7%/pre-metamorphic 27%) and Pseudomonadaceae being more abundant in metamorphic samples (metamorphic 18.1%/pre-metamorphic 6.4%).

Bacterial alpha diversity was significantly different between sample types (metamorphic, pre-metamorphic, sediment and water) as measured by ASV richness (Kruskal–Wallis (KW), \(\chi 2\) = 278.46, p-value < 0.001, df = 3), Shannon index (KW, \(\chi 2\) = 276.28, p-value < 0.001, df = 3) and Faith´s phylogenetic diversity (PD) (KW, \(\chi 2\) = 286.91, p-value < 0.001, df = 3) (Fig. 1B). Post hoc pairwise comparisons for each alpha diversity index showed significant differences among all sample types (Additional file 2: Table S1) except for metamorphic salamanders and water in ASV richness (Wilcoxon, p-value = 0.48) and Shannon diversity index (Wilcoxon, p-value = 0.66). Sediment samples showed the highest alpha diversity values while pre-metamorphic salamanders always had the lowest values.

Bacterial community composition based on the weighted UniFrac distance matrix varied significantly among sample types (PERMANOVA, pseudo-F = 64.76, p-value < 0.001, df = 3) (Fig. 1C, Additional file 2: Table S2). Dispersion significantly differed among sample types according to the permutational test (PERMPUTEST, F = 34.5, p-value = 0.001, df = 3) (Fig. 1D, Additional file 2: Table S3).

The skin bacterial composition of A. altamirani is mainly influenced by metamorphosis

Clear differences in skin bacterial alpha and beta diversity were found between metamorphic and pre-metamorphic salamanders (Fig. 1B, C, D). To look deeper into the bacterial taxa driving these differences we used an analysis of composition of microbiomes (ANCOM) which identified 45 bacterial families (out of 392 families in the axolotl skin samples) that were differentially abundant between metamorphic and pre-metamorphic samples (Fig. 2). Most of these bacterial families (40 out of 45) were enriched in metamorphic samples, being Verrucomicrobiaceae, Caulobacteraceae and Sphingomonadaceae the families with higher W values. In contrast, five bacterial families were enriched in pre-metamorphic samples with Burkholderiaceae, Chitinophagaceae being the families with higher W values.

Fig. 2
figure 2

ANCOM results showing differentially abundant bacterial families between metamorphic and pre-metamorphic axolotls. Left panel shows ANCOM W values, the middle panel shows the relative proportion for each bacterial family, and the right panel shows the best taxonomic assignments according to the SILVA database at order (O), class (C) or family (F) level. Circles and bars are color-coded according to the host metamorphic status

Full size image

We identified the core skin bacteria in metamorphic and pre-metamorphic A. altamirani axolotls, as those ASVs shared in ≥ 90% of the samples on each specific morph. Four ASVs represented the bacterial core of metamorphic axolotls accounting for a cumulative relative abundance of 16.26% of the ASVs. Meanwhile, two ASVs represented the bacterial core of pre-metamorphic axolotls accounting for 45.78% of the relative abundance (Table 2).

Table 2 Amplicon sequence variants (ASVs) defining the bacterial core of the microbiota of metamorphic and pre-metamorphic A. altamirani
Full size table

We also identified the core bacteria of environmental samples (Additional file 2: Table S4) and we found that metamorphic axolotls shared two core ASVs with the water core and another one with the sediment core. Core bacteria of pre-metamorphic axolotls are not present in the core of environmental samples.

Seasonality and location differentially influence skin bacterial diversity in metamorphic and pre-metamorphic axolotls

Physicochemical variables measured at each sampling location (pH, conductivity, dissolved oxygen, maximin, minimum mean and delta seasonal temperatures) varied significantly across seasons (MANOVA, Wilks = 0.002, p-value < 0.001, df = 3) and sampling locations (MANOVA, Wilks = 0.0009, p-value < 0.001, df = 3). While all physicochemical variables varied across seasons, dissolved oxygen was the only variable that did not vary between sampling locations (Additional file 2: Table S5).

Alpha PD of metamorphic axolotls varied significantly across seasons (KW, \(\chi 2\) = 13.69, p-value = 0.003, df = 3) (Fig. 3A) and post-hoc pairwise comparisons showed that only the transition between winter-spring was significant (Wilcoxon, p-value = 0.005) (Additional file 2: Table S6). In contrast, PD of pre-metamorphic A. altamirani (Fig. 3B) did not differ across consecutive seasons (KW, \(\chi 2\) = 0.21, p-value = 0.97, df = 3) (Additional fiie 2: Table S6).

Fig. 3
figure 3

Seasonal influence over metamorphic and pre-metamorphic skin bacterial diversity. A Phylogenetic diversity (PD) across seasons in metamorphic samples. Letters a–d indicate statistically significant comparisons. B PD across seasons in pre-metamorphic samples. C Seasonal variation of pH, delta temperature and mean temperature of the stream water. D Principal coordinate analysis (PCoA) based on weighted UniFrac distances across seasons of metamorphic samples. E PCoA based on weighted UniFrac distances across seasons in pre-metamorphic samples. Circles in D and E panels are color-coded by season

Full size image

Additionally, we found that seasonality significantly influenced skin bacterial community composition (PERMANOVA, pseudo-F = 12.37, p-value < 0.001, df = 3) (Fig. 3D, Additional file 2: Table S7), but not dispersion (PERMUTEST, F = 1.4, p-value = 0.24, df = 3) (Additional file 2: Table S8) of metamorphic axolotls. Seasonality also influenced skin bacterial community composition (PERMANOVA, pseudo-F = 15.69, p-value < 0.001, df = 3) (Fig. 3E, Additional file 2: Table S9) of pre-metamorphic axolotls, in addition we found that dispersion significantly differed between seasons (BETADISPER, F = 2.7, p-value = 0.038, df = 3) (Additional file 2: Table S10). Specifically, pairwise PERMANOVAs showed that metamorphic samples differed between winter-spring seasons (PERMANOVA, pseudo-F = 14.92, p-value = 0.001, df = 1), while pre-metamorphic skin microbiota differed between autumn–winter (PERMANOVA, pseudo-F = 13.47, p-value < 0.001, df = 1) and winter-spring seasons (PERMANOVA, pseudo-F = 12.61, p-value < 0.001, df = 1).

Three bacterial families were identified by ANCOM as differentially abundant in metamorphic samples between winter-spring seasons (Fig. 4). In the case of pre-metamorphic individuals, ANCOM identified three bacterial families that were differentially abundant between autumn–winter and eleven families differentially abundant between winter-spring (Fig. 4). Pseudomonadaceae, Aquaspirillaceae and Shewanellaceae were significantly enriched in both metamorphic and pre-metamorphic axolotls during winter and spring seasons. However, Pseudomonadaceae was more abundant in metamorphic axolotls during spring and more abundant in pre-metamorphic axolotls during winter. Shewanellaceae was more abundant in winter, and Aquaspirillaceae was present in the winter and completely absent in the spring for both metamorphic and pre-metamorphic axolotls.

Fig. 4
figure 4

ANCOM results showing differentially abundant bacterial families in metamorphic and pre-metamorphic axolotls across consecutive seasons: autumn to winter seasons in pre-metamorphic axolotls, and winter to spring seasons for metamorphic and pre-metamorphic axolotls. From left to right: ANCOM comparisons color-coded by season, ANCOM W values, the relative bacterial family proportion and the best taxonomic assignment according to SILVA at order (O), class (C) or family (F) level. Circles and bars are color-coded by season. Shared bacterial families between metamorphic and pre-metamorphic axolotls between winter and spring seasons are shown in bold

Full size image

When analyzing the effect of location in the skin bacterial diversity, we found that PD of metamorphic samples differed significantly between sampling locations (KW, \(\chi 2\) = 9.69, p-value = 0.02, df = 3), however post hoc paired test showed that PD only differed significantly between sites 2 and 3 (Additional file 1: Figure S2A). Bacterial PD of pre-metamorphic samples also varied significantly between sampling locations (KW, \(\chi 2\) = 40.9, p-value = 6.71e-9, df = 3). Post hoc test showed that most pairwise comparisons were significant with the exception of sites 1 and 3 and sites 2 and 3 (Additional file 1: Figure S2C, Additional file 2: Table S11). Skin bacterial community composition was also influenced by sampling location in metamorphic (PERMANOVA, pseudo-F = 2.71, p-value = 0.006, df = 3) and pre-metamorphic samples (PERMANOVA, pseudo-F = 31.34, p-value = 0.001, df = 3) (Additional file 1: Figure S2B, D). Pairwise comparisons showed that bacterial community composition only differed between sites 2 and 3 in metamorphic axolotls (Additional file 2: Table S12), while community composition differed between all sampling locations for pre-metamorphic samples (Additional file 2: Table S13). Interestingly dispersion did not vary across localities for metamorphic axolotls (PERMUTEST, F = 0.29, p-value = 0.8, df = 3) (Additional file 2: Table S14), but we found significant differences in dispersion for pre-metamorphic axolotls (PERMUTEST, F = 6.68, p-value = 0.01, df = 3) (Additional file 2: Table S15).

Biotic and abiotic factors influence the skin bacterial community structure of A. altamirani

Our results showed that bacterial community composition of A. altamirani skin is influenced by seasonality and location. To assess the specific influence of all the biotic and abiotic factors measured in this study we performed a distance-based Redundancy Analysis (dbRDA) on the skin bacterial beta diversity. After forward model selection, that incorporates all the variables measured, only the biotic and abiotic factors that better explained community composition were included in the dbRDA regression model: host metamorphic status, host weight, pH, dissolved oxygen, conductivity, mean temperature, season delta temperature (difference between the maximum and minimum seasonal temperature) and site elevation.

The dbRDA calculated eight canonical components for the PCA, however anova.cca (by = axis) showed that only four of these canonical components were statistically significant. These four statistically significant canonical axes explained 26.47% of the variation in the weighted UniFrac distance matrix (Fig. 5, Additional file 2: Table S16). Permutational analysis (anova.cca, by = terms) over each variable in the model showed that the metamorphic status of the host (PERMANOVA, pseudo-F = 39.1, p-value = 0.001) had the greatest effect-size over the variation, followed by seasonal delta temperature (PERMANOVA, pseudo-F = 19.8, p-value = 0.001), pH (PERMANOVA, pseudo-F = 15.85, p-value = 0.001) and seasonal mean temperature (PERMANOVA, pseudo-F = 12.05, p-value = 0.001) (Fig. 3C, Table 3).

Fig. 5
figure 5

Distance based redundancy analysis of A. altamirani skin bacterial communities. Distances in the PCA are based on a weighted UniFrac distance matrix. Vector directions indicate the type of correlation of each predictor variable. Distance of each sample with respect to vectors highlight the weight of the correlation with a given predictor variable. Non quantitative variables are represented as centroids (outlined circles larger). Circles are color-coded by host metamorphic status

Full size image
Table 3 Permutational like ANOVA results of each variable introduced in the dbRDA regression model
Full size table

Skin bacterial diversity of A. altamirani is not influenced by Bd infection status but specific bacterial taxa abundances correlate with infection intensity

Pathogen prevalence and infection intensity were conducted by Basanta et al. [49]. Briefly we found a Bd prevalence of 70.3% across all samples specifically 54 (out of 85) metamorphic and 142 (out of 194) pre-metamorphic axolotls resulted positive for Bd infection.

Alpha PD did not differ between infected and non-infected samples in both metamorphic (KW, \(\chi 2\) = 0.09, p-value = 0.76, df = 1) (Fig. 6A) and pre-metamorphic (KW, \(\chi 2\) = 0.51, p-value = 0.47, df = 1) A. altamirani samples (Fig. 6C). Additionally, beta diversity based on the weighted UniFrac distance matrix did not vary between infected and non-infected samples for metamorphic (PERMANOVA, pseudo-F = 1.37, p-value = 0.19, df = 1) (Fig. 6B) and pre-metamorphic axolotls (PERMANOVA, pseudo-F = 2.45, p-value = 0.08, df = 1) (Fig. 6D).

Fig. 6
figure 6

A. altamirani skin bacterial diversity with respect to Bd infection status. A Alpha phylogenetic diversity (PD) between infected and non-infected in metamorphic axolotls. B Principal coordinate analysis (PCoA) based on weighted UniFrac distances in infected vs non-infected metamorphic axolotls. C PD between infected and non-infected in pre-metamorphic axolotls. D PCoA based on weighted UniFrac distances in infected vs non-infected in pre-metamorphic axolotls. Circles are color-coded by Bd infection status

Full size image

Even though alpha and beta diversity did not vary according to Bd infection status, Kendall’s correlation test showed that the relative abundance of 139 and 129 ASV present in infected metamorphic and pre-metamorphic samples respectively, significantly correlated with pathogen infection loads (Additional file 1: Figure S3). Specifically, 116 (out of 139) and 52 (out of 128) bacterial ASVs had positive correlations with pathogen infection loads in metamorphic and pre-metamorphic samples, respectively.

Almost all the ASVs that correlated with pathogen load in metamorphic samples had low relative abundances ranging from 0.001 to 0.67% (Additional file 1: Figure S3A), while in pre-metamorphic samples, correlated ASVs ranged from 0.001 to 28.5% (Additional file 1: Figure S3B). Among the ASVs that correlated with pathogen infection intensities, twelve of them were shared between metamorphic and pre-metamorphic axolotls and half of these ASVs had a differential type of correlation between morphs (e.g., positive in metamorphic and negative in pre-metamorphic) (Additional file 2: Table S17). Two of these ASVs were classified as members of the Chitinophagaceae family and both were positively correlated with Bd infection intensity in metamorphic axolotls and negatively correlated in pre-metamorphic axolotls.

Discussion

The aim of this study was to evaluate the influence of metamorphosis, seasonality and pathogen presence over the skin microbiota of the axolotl A. altamirani. Since this is the first study exploring the skin microbiota of A. altamirani, we also evaluated if skin bacterial diversity differed from environmental bacterial communities of the streams where this species inhabits.

Consistent with previous studies showing differences between amphibian skin microbiota and their surrounding environmental bacterial communities [20, 50], we found that A. altamirani skin bacterial microbiota significantly differed from environmental samples, and that a great portion of the ASVs were unique to each sample type (ie, sediment, water, metamorphic and pre-metamorphic axolotls, Fig. 1A), supporting the idea that the amphibian skin hosts a distinctive bacterial repertoire compared to the environmental samples [18, 35, 37, 43, 51]. Our results highlight differences on the microbial diversity between the bacterial communities of the skin and the environment. We identified differences in alpha and beta diversity compared to water and sediment. Also, we identified core bacteria that were unique to axolotls or were clearly enriched on their skins compared to the environment.

Several studies have shown that amphibian skin microbiota varies significantly across host development [26, 27, 41]. These studies focused on amphibian species that transition from an aquatic larval stage to a terrestrial adult stage [22, 36, 37, 52], making it difficult to tease apart the effects of host development stage and habitat conditions on skin microbial diversity [17, 18]. For species where adult and larval stages coexist in the same aquatic environment (i.e., newts), host developmental stage has had contrasting results in different species; for example adult and larvae of Lissotriton boscai showed clear differences in skin bacterial community composition, however this pattern was not observed in Triturus marmoratus [52].

In this study, we evaluated the influence of metamorphosis over skin bacterial diversity on a paedomorphic salamander species (axolotl) in which metamorphic and pre-metamorphic stages coexist in permanent streams along their life cycle [47, 53]. Our results showed that A. altamirani skin bacterial communities are strongly shaped by metamorphosis. Specifically, we found that pre-metamorphic individuals harbor less diverse and more variable skin bacterial communities compared to metamorphic individuals.

These differences could be explained by differences in skin mucus composition, immune response, or gene expression before and after metamorphosis as it has been proposed that mucus chemical composition (e.g., production of antimicrobial peptides) play a critical role in shaping the skin microbiota as well as in defense against pathogens [30, 54, 55]. Antimicrobial peptide repertory of the skin changes through development in some frog species [56], and some bacteria can induce the synthesis of specific antimicrobial peptides [51]. In addition, the number and distribution of Leydig cells, which have been associated with the secretion of mucus [57], changes across urodele development [58, 59].

In addition, the core microbiota analysis and ANCOM results shown here highlighted differences in composition between metamorphic and pre-metamorphic axolotls. Specifically, pre-metamorphic skin microbiota was composed by fewer core members and had less differentially abundant bacterial ASVs when compared to metamorphic skin microbiota. It is interesting to highlight that both analyses identified that families Chitinophagaceae and Burkholderiaceae were enriched in pre-metamorphic samples, specially two ASVs from these families conform the core microbiota of pre-metamorphic axolotls which account 45.7% of relative abundance in these samples.

Bacteria from the family Chitinophagaceae and Burkholderiaceae have been isolated from other amphibian hosts and have shown the ability to inhibit Bd [60,61,62]. Moreover, some members of the Chitinophagaceae family such as Chitinophaga pinenis can degrade chitin [63] which is a main component of fungal cell wall. In our study, the high prevalence of these bacterial families on the skin of A. altamirani could suggest that these bacteria play a protective role against chytrid pathogens.

Temporal and spatial dynamics of amphibian skin microbiota have been linked to variation in environmental factors such as temperature, precipitation or elevation [25, 26, 34, 35, 52]. Specifically, temperature fluctuations over short periods of time [22] and seasonal variations (dry–wet) [38] have been linked to differences in bacterial skin diversity on amphibians inhabiting aquatic environments. Our results showed that seasonal variation of temperature (delta temperature and mean temperature), pH, conductivity, and dissolved oxygen influence axolotl skin bacterial diversity.

Previous studies have shown that spatial variation has an influence on skin bacterial diversity of terrestrial salamanders [14, 41, 64]; for example populations of Ensatina eschscholtzii from different geographic locations vary in bacterial community composition [15]. In this study we found that sampling location significantly influences skin bacterial diversity, and this effect is stronger in pre-metamorphic axolotls. Considering that the main source of diversity of the skin microbiota are the environmental microbial communities and that they vary in response to environmental variation [65,66,67] it is likely that the skin microbiota reflect to some extent the environmental variations across localities as seen in the case of pre-metamorphic axolotls. It has been shown that skin bacterial diversity vary in response to precipitation [19, 23] temperature [22] or elevation gradients [24, 41, 42]. However, genetic differences across populations could also explain some of our results, since a previous study showed that A. altamirani populations of sites 2 and 3 are genetically differentiated [68]. Additional work is needed to tease apart the effects of environment and host genetics on the skin microbial diversity of A. altamirani.

Stability as a characteristic of an ecological community could be defined as the response to disturbance, comprising resilience and resistance against external disturbances [69, 70]. It has been shown that the stability of the amphibian skin microbiota can change after the experimental exposure to fungal [71] and viral pathogens [72]. Also, it has been shown that skin microbiomes with higher diversity are less stable to a pathogen induced disturbance [72].

Environmental variation across seasons is a different kind of perturbation for microbial communities, and it has been shown that seasonality influences microbial communities of soil [73,74,75], water [76, 77] and host-associated microbiomes [78, 79]. We found that skin bacterial communities of A. altamirani vary across seasons, particularly in pre-metamorphic axolotls which have a lower bacterial diversity compared to metamorphic axolotls.

Together our results suggest that more diverse bacterial communities (as the ones present in metamorphic axolotls) allow for a more stable microbiota that could be either more resistant or resilient to the environmental variation. Similar patterns of diversity—stability trough time have been described in populations of Rana sierrae [39]. Further studies are needed to evaluate if these patterns of stability across seasons influence the function of the skin bacterial communities of A. altamirani, as it has been shown that less stable bacterial communities show less functional redundancy [80].

Disruption of the skin microbiota following Bd infections has been previously documented in naive amphibian populations before and after Bd infection [20, 21], and in populations with different pathogen intensities where Bd seems to be present in an enzootic stage [23, 44]. Even when Bd was highly prevalent in A. altamirani populations [49], we did not find any significant influence of Bd presence over the skin bacterial diversity.

Of the 279 axolotls sampled only two individuals exhibited clear signs of infection [81] (lethargy, skin ulceration and extreme skin sloughing) and they died soon after we sampled them. Apart from these two cases, the remaining individuals showed no signs of infection. These observations suggest that this population is able to tolerate Bd infection. Further studies testing the survival rates of A. altamirani against Bd are needed to elucidate if this species is resistant or susceptible to chytridiomycosis.

It has been shown that Bd presence have contrasting effects over skin microbiota diversity inducing changes in skin microbiota composition following infection [21, 23, 71] or not influencing diversity of skin microbial communities [42, 44] as we found in this study. However, it also has been shown that relative abundances of some bacterial members of the skin microbiota correlates with chytrid infection intensity [19, 44, 45] and its suggested that according to the type of correlation these groups could act as anti Bd bacteria [19]. We identified several bacteria with positive and negative correlations with Bd infection intensities and most of these ASVs exhibited low relative abundances.

Observations in several amphibian species indicate that certain bacteria with properties such as biofilm formation [82] or putative inhibitory ability [55] are positively or negatively corelated with a decrease of Bd prevalence. Thus, we expect that bacteria with negative correlations to infection intensity could be important in the defense against Bd in A. altamirani. However, these Bd-inhibitory bacteria exhibited reduced abundances over the amphibian skin [83, 84].

Inhibitory potential against Bd has been described for several bacterial isolates mainly form Burkholderiaceae, Yersiniceae, Pseudomonadaceae or Xanthomondaceae families [60, 85,86,87,88], We found that Burkholderiaceae and Chitinophagaceae were highly abundant over A. altamirani skin. In line with our results, high abundance of Burkholderiaceae in Anaxyrus boreas skin microbiota correlated with reduced fungal presence over the skin during early life stages [27]. Additionally, populations of R. sierrae with contrasting Bd loads (high vs low) exhibited differential abundances of Burkholderiaceae [21, 44]. In the case of Chitinophagaceae little is known about their inhibitory ability against Bd with only few isolates considered as Bd-inhibitory strains [25], and further work is needed to elucidate if members of this bacterial family present on A. altamirani skin display inhibitory functions against Bd.

Conclusion

Our results show that host metamorphic status is a major determinant of A. altamirani, influencing diversity and structure of skin symbiotic bacterial communities. To our knowledge this study is the first to address how the effects of environmental variation over the skin bacterial communities are dependent on the amphibian developmental stage; we demonstrate that seasonal environmental variation significantly influences bacterial skin diversity of A. altamirani, and that metamorphic and pre-metamorphic axolotls respond differently to environmental variation. Despite a growing body of literature suggesting that Bd influences skin bacterial diversity we did not find such an effect. Nonetheless, we found that particular bacterial taxa are likely interacting with Bd. Further studies using metagenomics and cultivation techniques could elucidate if changes in skin microbiota across development and across seasons are reflecting functional differences regarding Bd inhibition or other host symbiotic traits [89, 90].

Methods

Sample collection

Skin samples were collected during four sampling periods at three-month intervals (July 2019, October 2019, January 2020, and April 2020) spanning all the seasons of a whole year at four localities at La Sierra de Cruces, Estado de México, México (Table 1, Additional file 2: Table S18). Individuals of A. altamirani were captured at each location using dip nets and held individually in sterile plastic containers filled with stream water until swabbing. Sampling occurred for three consecutive hours across a 150 m transect along each stream. Each captured salamander was manipulated with sterile nitrile gloves, rinsed with 25 ml of sterile deionized water to remove transient microorganisms from the skin and swabbed 30 times (five times in their ventral and dorsal surface each and five times in each limb joint) using sterile rayon swabs (MWE, Corsham UK). Swabs were stored in 1.5 ml microcentrifuge tubes containing 170 μl of DNA/RNA Shield (Zymo Research, Irvine, USA) and kept at 4 °C during field work. Once in the laboratory tubes were stored at − 80 °C until processing. Immediately after swabbing morphometric measurements of weight, tail and body length were registered for each individual. Once all axolotls were swabbed and measured, they were released at the same site of capture. Sampling was approved by Subsecretaría de Gestión para la Protección Ambiental under the permit number: SGPA/DGVS/5673/19.

For the purposes of this work, we classified axolotl samples as metamorphic and pre-metamorphic according to the presence or absence of gills as reported previously [59]. Recognizing that gilled individuals of A. altamirani could be either juvenile or paedomorphic adults, we classified non-gilled axolotls as metamorphic and gilled axolotls as pre-metamorphic respectively in order to evaluate the effect of the metamorphic status of the host.

Additionally, five samples of sediment and water were collected at each location in all sampling periods. Water samples were obtained by submerging a sterile rayon swab at approximately 20 cm deep inside water for 10 s, and sediment samples were obtained by submerging swabs inside the bottom sediment of the stream for 10 s [50].

Environmental characterization

Stream water temperature was recorded at 1 h intervals during one year at each sampling location using Onset HOBO dataloggers (Onset Computer Corporation, Bourne, USA) from June 2019 to April 2020. Additionally, pH, dissolved oxygen and conductivity of the water was registered using a HANNA multiparameter HI98194 (HANNA Instruments, USA) during each sampling. Measurements were taken at each location in triplicate across 10 m transects. To evaluate if these physicochemical variables vary between seasons and sampling location, we applied a two-way MANOVA test in R (v 4.0.2).

DNA extraction and sequencing

Amplicon libraries of the 16S rRNA gene spanning the V4 region were constructed using 515F/806R primers following the Earth Microbiome Project standard protocol (www.earthmicrobiome.org) and previously published studies [50, 91]. In brief, DNA was extracted from skin and environmental swabs using the Qiagen DNeasy Blood and Tissue kit (Qiagen, Valencia, USA) following manufacturer instructions with an initial lysozyme incubation step at 37° for 1 h. Samples were PCR amplified in triplicate plus one negative control per sample, PCR products and negative controls were verified in 1% agarose gels, and PCR products were pooled in one tube per sample. Pools were quantified using a Qubit 4.0 fluorometer (Invitrogen, Thermo Fisher Scientific, Waltham, USA), samples were pooled in two amplicon libraries at a concentration of 240 ng per sample (221 and 217 samples each). Each pool was cleaned using the QIAquick PCR clean up kit (Qiagen, Valencia, USA). 16S amplicon libraries were sequenced in two sequencing runs (250 single end) using v2 Illumina chemistry at Dana-Farber Cancer Institute of Harvard University.

Bioinformatic pipeline

Sequences were processed using Quantitative Insights Into Microbial Ecology (QIIME v2-2020.2) [92]. A total of 8,434,775 and 8,821,621 demultiplexed raw sequences were obtained from the sequencing facility for each sequencing run respectively. Prior to quality control primers were trimmed from the sequences using the cut-adapt plugin in Qiime2 then sequences were quality filtered and denoised independently for each run using the DADA2 plugin to obtain two single feature table. Feature tables obtained for each sequencing run were merged to generate a final Amplicon Sequence Variant (ASV) table containing 14,415,727 reads with a mean read depth of 32,900 reads per sample.

A phylogenetic tree was generated using the representative sequences of the ASV table using the q2-phylogeny plugin which first uses mafft to perform sequence alignment and then generate a phylogeny using FastTree. Samples were rarefied at 10,000 reads per sample according to observed ASV rarefaction curves in order to preserve the largest number of samples and sequences. After denoising and rarefaction at 10,000 reads per sample, seven axolotl (out of 279, three from metamorphic and four from pre-metamorphic) and eight environmental (out of 159, 3 from sediment and five from water) samples were discarded due to low read counts (< 10,000 reads per sample). The rarefied table containing 10,000 reads per sample was used for all further analyses including to calculate alpha and beta diversity metrics using the q2-diversity plugin. Taxonomy was assigned using a naive Bayesian classifier pre-trained for the V4 region (515F/806R 16 s rRNA) on the SILVA 132 99% database [93].

Microbial diversity and composition analyses

Statistical analyses for alpha and beta diversity were carried out using the rarefied table at 10,000 sequences per sample; these analyses were computed in R (v 4.0.2) unless otherwise stated. We first perform Kruskal–Wallis (KW) and post hoc Wilcoxon ranks sum test were used to determine differences in alpha diversity (Shannon, Faith’s Phylogenetic Diversity (PD) and ASV richness) between sample types (metamorphic, pre-metamorphic, sediment, and water. In addition we perform KW and post hoc Wilcoxon ranks sum test to evaluate the influence of seasonality, sampling location and Bd infection status over the skin microbiota of metamorphic and pre-metamorphic axolotls individually.

Beta diversity was evaluated using a weighted UniFrac distance matrix generated using the rarefied table at 10,000 sequences per sample to determine differences in bacterial community structure across sample types. In addition, we generated two independent weighted UniFrac distance matrices for metamorphic and pre-metamorphic axolotls to evaluate the influence of seasonality, sampling location and Bd infection status. Statistical comparisons were conducted with permutational multivariate analyses (PERMANOVA) using the q2-diversity plugin in Qiime2 (v 2020.2). Beta diversity dispersion was calculated from the each weighted UniFrac distance matrix using the function betadisper in the vegan package [94], and then we applied PERMUTEST based on 999 permutations to identify significant differences for dispersion, specifically we evaluate dispersion between sample types, as well between seasonality and sampling locations for the metamorphic and pre-metamorphic distance matrices.

ANCOM [95] was used to identify bacterial families that were differentially abundant between metamorphic and pre-metamorphic salamanders and between samples from consecutive seasons (summer-autumn, autumn–winter, winter-spring). Prior to analysis low abundant ASVs (< 50 reads) were filtered out and then we collapsed all ASVs at family level, ANCOM was performed using the q2-composition plugin in Qiime2. Briefly, ANCOM applies a centered log ratio transformation on the relative abundance of each bacterial family and tests the null hypothesis that mean log absolute abundance of each family does not differ between sample types. An internal statistic (W) is calculated each time a taxon rejects this null hypothesis, then ANCOM generates an empirical distribution using W values in order to test which taxon in this case which bacterial families are differentially abundant between samples. ANCOM between consecutive seasons was only applied if PERMANOVA results showed significant differences between consecutive seasons (winter-spring for metamorphic salamanders and autumn–winter and winter-spring for pre-metamorphic salamanders).

Core microbiome was calculated independently for metamorphic and pre-metamorphic axolotls using the feature-table plugin in Qiime2. In brief, we generated four feature tables that contain all the ASVs present in each sample type (metamorphic, pre-metamorphic, sediment and water samples). Then we identify all the ASVs present in ≥ 90% of the samples of each sample type, using the core-features function.

Additionally, correlations between the relative abundance of each ASV of the infected samples and Bd infection intensities were calculated with Kendall rank correlation coefficient correcting for multiple comparisons (Benjamini-Hochberg) using cor.test function of the stats package in R [96]. To generate graphics for all the results Qiime2 artifacts were imported to R using the package qiime2R [97], then figures were generated using packages ggplot2 [98, 99], Fantaxtic [100] and UpSetR [101], color pallet of the figures are colorblind friendly and were selected from the MetBrewer package in R [102].

Biotic and abiotic factors influencing the skin microbial structure

In order to explore the specific influence of biotic (developmental stage, weight, tail length, snout vent length, Bd presence and Bd infection intensity) and abiotic factors (pH, conductivity, dissolved oxygen, mean season temperature, delta season temperature, and elevation) on skin bacterial community composition, we applied a distance-based redundancy analysis (dbRDA) on the weighted UniFrac distance matrix using the capscale function of the vegan package [94]. dbRDA is a canonical ordination method that applies multiple linear regression to a distance matrix and then computes a principal component analysis (PCA) [103]. Prior to analyses non-categorical biotic and abiotic variables were z-scored to control for differences in magnitudes between factors. The ordistep function of the vegan package [94] was used for model selection in both directions with 999 permutations to select the best regression model. Once the dbRDa was obtained anova.cca function was used to perform an ANOVA like permutation test to evaluate the significance of each calculated canonical axis (anova.cca, by = axix) and the specific significance of each factor in the regression model (anova.cca, by = terms).

Availability of data and materials

All 16 s rRNA gene raw data in this study are publicly available at the NCBI SRA under BioProject PRJNA819099. Sample metadata, output data of DADA2, R and Qiime2 scripts for analysis and figures included in this manuscript are available at https://github.com/EmanuelMartinez-Ugalde/A.-altamirani-16S-skin-microbiota.

Abbreviations

ANCOM:

Analysis of composition of microbiomes

ASV:

Amplicon sequence variant

Bd:

Batrachochytrium dendrobatidis

KW:

Kruskal–Wallis

dbRDA:

Distance based redundancy analysis

MANOVA:

Multivariate analysis of variance

PD:

Phylogenetic diversity

PERMANOVA:

Permutational multivariate analysis of variance

PERMUTEST:

Permutation-based test of multivariate homogeneity of group dispersions

References

  1. Bahrndorff S, Alemu T, Alemneh T, Lund NJ. The microbiome of animals: implications for conservation biology. Int J Genomics. 2016;2016:5304028.

    Article  Google Scholar 

  2. Hird SM. Evolutionary biology needs wild microbiomes. Front Microbiol. 2017;8:725.

    Article  Google Scholar 

  3. Vaelli PM, Theis KR, Williams JE, O’Connell LA, Foster JA, Eisthen HL. The skin microbiome facilitates adaptive tetrodotoxin production in poisonous newts. Elife. 2020;9:e53898.

    Article  Google Scholar 

  4. Ross AA, Rodrigues Hoffmann A, Neufeld JD. The skin microbiome of vertebrates. Microbiome. 2019;7:79.

    Article  Google Scholar 

  5. Comizzoli P, Power ML, Bornbusch SL, Muletz-Wolz CR. Interactions between reproductive biology and microbiomes in wild animal species. Anim Microbiome. 2021;3:87.

    Article  CAS  Google Scholar 

  6. Lemieux-Labonté V, Dorville NAS-Y, Willis CKR, Lapointe F-J. Antifungal potential of the skin microbiota of hibernating big brown bats (Eptesicus fuscus) infected with the causal agent of white-nose syndrome. Front Microbiol. 2020;11:1776.

    Article  Google Scholar 

  7. Hill AJ, Leys JE, Bryan D, Erdman FM, Malone KS, Russell GN, et al. Common cutaneous bacteria isolated from snakes inhibit growth of Ophidiomyces ophiodiicola. EcoHealth. 2018;15:109–20.

    Article  Google Scholar 

  8. Harris RN, James TY, Lauer A, Simon MA, Patel A. Amphibian pathogen Batrachochytrium dendrobatidis is inhibited by the cutaneous bacteria of amphibian species. EcoHealth. 2006;3:53.

    Article  Google Scholar 

  9. 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.

    Article  CAS  Google Scholar 

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

    Article  Google Scholar 

  11. Woodhams DC, Bosch J, Briggs CJ, Cashins S, Davis LR, Lauer A, et al. Mitigating amphibian disease: strategies to maintain wild populations and control chytridiomycosis. Front Zool. 2011;8:8.

    Article  Google Scholar 

  12. Bletz MC, Loudon AH, Becker MH, Bell SC, Woodhams DC, Minbiole KPC, et al. Mitigating amphibian chytridiomycosis with bioaugmentation: characteristics of effective probiotics and strategies for their selection and use. Ecol Lett. 2013;16:807–20.

    Article  Google Scholar 

  13. 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.

    Article  CAS  Google Scholar 

  14. 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.

    Article  Google Scholar 

  15. Prado-Irwin SR, Bird AK, Zink AG, Vredenburg VT. Intraspecific variation in the skin-associated microbiome of a terrestrial salamander. Microb Ecol. 2017;74:745–56.

    Article  CAS  Google Scholar 

  16. McKenzie VJ, Bowers RM, Fierer N, Knight R, Lauber CL. Co-habiting amphibian species harbor unique skin bacterial communities in wild populations. ISME J. 2012;6:588–96.

    Article  CAS  Google Scholar 

  17. Loudon AH, Woodhams DC, Parfrey LW, Archer H, Knight R, McKenzie V, et al. Microbial community dynamics and effect of environmental microbial reservoirs on red-backed salamanders (Plethodon cinereus). ISME J. 2014;8:830–40.

    Article  CAS  Google Scholar 

  18. Walke JB, Becker MH, Loftus SC, House LL, Cormier G, Jensen RV, et al. Amphibian skin may select for rare environmental microbes. ISME J. 2014;8:2207–17.

    Article  CAS  Google Scholar 

  19. Familiar López M, Rebollar EA, Harris RN, Vredenburg VT, Hero J-M. Temporal Variation of the skin bacterial community and Batrachochytrium dendrobatidis infection in the terrestrial cryptic frog Philoria loveridgei. Front Microbiol. 2017;8:2535.

    Article  Google Scholar 

  20. Kearns PJ, Fischer S, Fernández-Beaskoetxea S, Gabor CR, Bosch J, Bowen JL, et al. Fight fungi with fungi: antifungal properties of the amphibian mycobiome. Front Microbiol. 2017;8:2494.

    Article  Google Scholar 

  21. Jani AJ, Briggs CJ. The pathogen Batrachochytrium dendrobatidis disturbs the frog skin microbiome during a natural epidemic and experimental infection. Proc Natl Acad Sci. 2014;111:E5049-58.

    Article  Google Scholar 

  22. 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.

    Article  Google Scholar 

  23. Longo AV, Zamudio KR. Environmental fluctuations and host skin bacteria shift survival advantage between frogs and their fungal pathogen. ISME J. 2017;11:349–61.

    Article  Google Scholar 

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

    Article  Google Scholar 

  25. Barnes EM, Kutos S, Naghshineh N, Mesko M, You Q, Lewis JD. Assembly of the amphibian microbiome is influenced by the effects of land-use change on environmental reservoirs. Environ Microbiol. 2021;23:4595–611.

    Article  Google Scholar 

  26. Longo AV, Savage AE, Hewson I, Zamudio KR. Seasonal and ontogenetic variation of skin microbial communities and relationships to natural disease dynamics in declining amphibians. R Soc Open Sci. 2015;2:140377.

    Article  Google Scholar 

  27. Kueneman JG, Woodhams DC, Van Treuren W, Archer HM, Knight R, McKenzie VJ. Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas). ISME J. 2016;10:934–44.

    Article  Google Scholar 

  28. Demircan T, Ovezmyradov G, Yıldırım B, Keskin İ, İlhan AE, Fesçioğlu EC, et al. Experimentally induced metamorphosis in highly regenerative axolotl (Ambystoma mexicanum) under constant diet restructures microbiota. Sci Rep. 2018;8:10974.

    Article  Google Scholar 

  29. Brown DD, Cai L. Amphibian metamorphosis. Dev Biol. 2007;306:20–33.

    Article  CAS  Google Scholar 

  30. Rollins-Smith LA. Metamorphosis and the amphibian immune system. Immunol Rev. 1998;166:221–30.

    Article  CAS  Google Scholar 

  31. Newman RA. Ecological constraints on amphibian metamorphosis: interactions of temperature and larval density with responses to changing food level. Oecologia. 1998;115:9–16.

    Article  Google Scholar 

  32. Wilbur HM, Collins JP. Ecological aspects of amphibian metamorphosis: nonnormal distributions of competitive ability reflect selection for facultative metamorphosis. Science. 1973;182:1305–14.

    Article  CAS  Google Scholar 

  33. Azizi E, Landberg T. Effects of metamorphosis on the aquatic escape response of the two-lined salamander (Eurycea bislineata). J Exp Biol. 2002;205:841–9.

    Article  Google Scholar 

  34. Kueneman J, Bletz M, McKenzie V, Becker CG, Joseph M, Abarca J, et al. Community richness of amphibian skin bacteria correlates with bioclimate at the global scale. Nat Ecol Evol. 2019;3:1.

    Article  Google Scholar 

  35. 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.

    Article  Google Scholar 

  36. Xu L, Xiang M, Zhu W, Zhang M, Chen H, Huang J, et al. The behavior of amphibians shapes their symbiotic microbiomes. mSystems. 2020;5:e00626-20.

    Article  Google Scholar 

  37. Douglas AJ, Hug LA, Katzenback BA. Composition of the _North American wood frog (Rana sylvatica) bacterial skin microbiome and seasonal variation in community structure. Microb Ecol. 2021;81:78–92.

    Article  CAS  Google Scholar 

  38. Nava-González B, Suazo-Ortuño I, López PB, Maldonado-López Y, Lopez-Toledo L, Raggi L, et al. Inhibition of Batrachochytrium dendrobatidis infection by skin bacterial communities in wild amphibian populations. Microb Ecol. 2021;82:666–76.

    Article  Google Scholar 

  39. Ellison S, Knapp R, Vredenburg V. Longitudinal patterns in the skin microbiome of wild, individually marked frogs from the Sierra Nevada. California ISME Commun. 2021;1:45.

    Article  Google Scholar 

  40. Catenazzi A, Flechas SV, Burkart D, Hooven ND, Townsend J, Vredenburg VT. Widespread elevational occurrence of antifungal bacteria in Andean amphibians decimated by disease: a complex role for skin symbionts in defense against chytridiomycosis. Front Microbiol. 2018. https://doi.org/10.3389/fmicb.2018.00465.

    Article  Google Scholar 

  41. Bresciano JC, Salvador CA, Paz-y-Miño C, Parody-Merino AM, Bosch J, Woodhams DC. Variation in the presence of anti-Batrachochytrium dendrobatidis bacteria of amphibians across life stages and elevations in Ecuador. EcoHealth. 2015;12:310–9.

    Article  CAS  Google Scholar 

  42. McKnight DT, Huerlimann R, Bower DS, Schwarzkopf L, Alford RA, Zenger KR. The interplay of fungal and bacterial microbiomes on rainforest frogs following a disease outbreak. Ecosphere. 2022;13:e4037.

    Article  Google Scholar 

  43. Bletz MC, Archer H, Harris RN, McKenzie VJ, Rabemananjara FCE, Rakotoarison A, et al. Host ecology rather than host phylogeny drives amphibian skin microbial community structure in the biodiversity hotspot of Madagascar. Front Microbiol. 2017. https://doi.org/10.3389/fmicb.2017.01530.

    Article  Google Scholar 

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

    Article  Google Scholar 

  45. Muletz-Wolz CR, Fleischer RC, Lips KR. Fungal disease and temperature alter skin microbiome structure in an experimental salamander system. Mol Ecol. 2019. https://doi.org/10.1111/mec.15122.

    Article  Google Scholar 

  46. Woolrich G, Smith G, Lemos-Espinal JA, Zamora ABE, Montoya-Ayala R. Observed localities for three endangered, endemic Mexican ambystomatids (Ambystoma altamirani, A. leorae, and A. rivulare) from central Mexico. Herpetol Bull. 2017;139:12–5.

    Google Scholar 

  47. Lemos-Espinal JA, Smith GR, Ruíz ÁH, Ayala RM. Stream use and population characteristics of the endangered salamander, Ambystoma altamirani, from the Arroyo Los Axolotes, State of Mexico, Mexico. Southwest Nat. 2016;61:28–32. https://doi.org/10.1894/0038-4909-61.1.28.

    Article  Google Scholar 

  48. Camacho ZAV, Smith GR, Ayala RM, Lemos-Espinal JA. Distribution and population structure of Ambystoma altamirani from the Llano de Lobos, State of México, Mexico. West North Am Nat. 2020. https://doi.org/10.3398/064.080.0210.

    Article  Google Scholar 

  49. Basanta MD, Anaya-Morales SL, Martínez-Ugalde E, González Martínez TM, Ávila-Akerberg VD, Vázquez Trejo M, et al. Metamorphosis and seasonality are major determinants of chytrid infection in a paedomorphic salamander. Anim Conserv. 2022. https://doi.org/10.1111/acv.12824.

    Article  Google Scholar 

  50. 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.

    Article  CAS  Google Scholar 

  51. Woodhams DC, Rollins-Smith LA, Reinert LK, Lam BA, Harris RN, Briggs CJ, et al. Probiotics modulate a novel amphibian skin defense peptide that is antifungal and facilitates growth of antifungal bacteria. Microb Ecol. 2020;79:192–202.

    Article  CAS  Google Scholar 

  52. 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.

    Article  CAS  Google Scholar 

  53. Hernández VV, Smith GR, Ayala RM, Lemos-Espinal JA. The relationship between body and substrate color for Ambystoma altamirani (Caudata: Ambystomatidae) from the Arroyo los Axolotes, Mexico. Phyllomedusa J Herpetol. 2020;19:243–51.

    Article  Google Scholar 

  54. Palacios-Martinez J, Caballero-Perez J, Espinal-Centeno A, Marquez-Chavoya G, Lomeli H, Salas-Vidal E, et al. Multi-organ transcriptomic landscape of Ambystoma velasci metamorphosis. Dev Biol. 2020;466:22–35.

    Article  Google Scholar 

  55. Jiménez RR, Carfagno A, Linhoff L, Gratwicke B, Woodhams DC, Chafran LS, et al. Inhibitory bacterial diversity and mucosome function differentiate susceptibility of appalachian salamanders to chytrid fungal infection. Appl Environ Microbiol. 2022;88:e01818-e1821.

    Article  Google Scholar 

  56. Woodhams DC, Bell SC, Bigler L, Caprioli RM, Chaurand P, Lam BA, et al. Life history linked to immune investment in developing amphibians. Conserv Physiol. 2016;4:cow025.

    Article  Google Scholar 

  57. Jarial MS. Fine structure of the epidermal Leydig cells in the axolotl Ambystoma mexicanum in relation to their function. J Anat. 1989;167:95–102.

    CAS  Google Scholar 

  58. Gerling S, D’Haese J, Greven H. Number and distribution of Leydig cells (LC) in the epidermis of the growing axolotl, Ambystoma mexicanum (Amphibia: Urodela). Vertebr Zool. 2012;62(1):97–111.

    Google Scholar 

  59. Ohmura H, Wakahara M. Transformation of skin from larval to adult types in normally metamorphosing and metamorphosis-arrested salamander Hynobius retardatus. Differentiation. 1998;63:237–46.

    Google Scholar 

  60. Bletz MC, Myers J, Woodhams DC, Rabemananjara FCE, Rakotonirina A, Weldon C, et al. Estimating herd immunity to amphibian chytridiomycosis in Madagascar based on the defensive function of amphibian skin bacteria. Front Microbiol. 2017;8:1751.

    Article  Google Scholar 

  61. 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. https://doi.org/10.3389/fmicb.2020.00101.

    Article  Google Scholar 

  62. 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. Ecology. 2015;96:595–595.

    Article  Google Scholar 

  63. Ramakrishna B, Vaikuntapu P, Mallakuntla MK, Bhuvanachandra B, Sivaramakrishna D, Uikey S, et al. Carboxy-terminal glycosyl hydrolase 18 domain of a carbohydrate active protein of Chitinophaga pinensis is a non-processive exochitinase. Int J Biol Macromol. 2018;115:1225–32.

    Article  CAS  Google Scholar 

  64. Belasen AM, Riolo MA, Bletz MC, Lyra ML, Toledo LF, James TY. Geography, host genetics, and cross-domain microbial networks structure the skin microbiota of fragmented Brazilian Atlantic forest frog populations. Ecol Evol. 2021;11:9293–307.

    Article  Google Scholar 

  65. Shen C, Xiong J, Zhang H, Feng Y, Lin X, Li X, et al. Soil pH drives the spatial distribution of bacterial communities along elevation on Changbai mountain. Soil Biol Biochem. 2013;57:204–11.

    Article  CAS  Google Scholar 

  66. Clark JS, Campbell JH, Grizzle H, Acosta-Martìnez V, Zak JC. Soil microbial community response to drought and precipitation variability in the Chihuahuan desert. Microb Ecol. 2009;57:248–60.

    Article  Google Scholar 

  67. Yao Z, Du S, Liang C, Zhao Y, Dini-Andreote F, Wang K, et al. Bacterial community assembly in a typical estuarine marsh with multiple environmental gradients. Appl Environ Microbiol. 2019;85:e02602-e2618.

    Article  CAS  Google Scholar 

  68. Monroy-Vilchis O, Heredia-Bobadilla R-L, Zarco-González MM, Ávila-Akerberg V, Sunny A. Genetic diversity and structure of two endangered mole salamander species of the Trans-Mexican Volcanic Belt. Herpetozoa. 2019;32:237–48.

    Article  Google Scholar 

  69. Shade A, Peter H, Allison S, Baho D, Berga M, Buergmann H, et al. Fundamentals of microbial community resistance and resilience. Front Microbiol. 2012;3:417.

    Article  Google Scholar 

  70. Song H-S, Renslow RS, Fredrickson JK, Lindemann SR. Integrating ecological and engineering concepts of resilience in microbial communities. Front Microbiol. 2015;6:1298.

    Article  Google Scholar 

  71. Jani AJ, Bushell J, Arisdakessian CG, Belcaid M, Boiano DM, Brown C, et al. The amphibian microbiome exhibits poor resilience following pathogen-induced disturbance. ISME J. 2021;15:1628–40.

    Article  CAS  Google Scholar 

  72. 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. https://doi.org/10.3389/fmicb.2019.02883.

    Article  Google Scholar 

  73. Shigyo N, Umeki K, Hirao T. Seasonal dynamics of soil fungal and bacterial communities in cool-temperate montane forests. Front Microbiol. 2019. https://doi.org/10.3389/fmicb.2019.01944.

    Article  Google Scholar 

  74. Luo X, Wang MK, Hu G, Weng B. Seasonal change in microbial diversity and its relationship with soil chemical properties in an Orchard. PLoS ONE. 2019;14:e0215556.

    Article  CAS  Google Scholar 

  75. Madegwa YM, Uchida Y. Land use and season drive changes in soil microbial communities and related functions in agricultural soils. Environ DNA. 2021;3:1214–28.

    Article  CAS  Google Scholar 

  76. Wilhelm SW, LeCleir GR, Bullerjahn GS, McKay RM, Saxton MA, Twiss MR, et al. Seasonal changes in microbial community structure and activity imply winter production is linked to summer hypoxia in a large lake. FEMS Microbiol Ecol. 2014;87:475–85.

    Article  CAS  Google Scholar 

  77. Mestre M, Höfer J, Sala MM, Gasol JM. Seasonal variation of bacterial diversity along the marine particulate matter continuum. Front Microbiol. 2020. https://doi.org/10.3389/fmicb.2020.01590.

    Article  Google Scholar 

  78. Palladino G, Biagi E, Rampelli S, Musella M, D’Amico F, Turroni S, et al. Seasonal changes in microbial communities associated with the jewel anemone Corynactis viridis. Front Mar Sci. 2021. https://doi.org/10.3389/fmars.2021.627585.

    Article  Google Scholar 

  79. Bird S, Prewer E, Kutz S, Leclerc L-M, Vilaça ST, Kyle CJ. Geography, seasonality, and host-associated population structure influence the fecal microbiome of a genetically depauparate Arctic mammal. Ecol Evol. 2019;9:13202–17.

    Article  Google Scholar 

  80. Biggs CR, Yeager LA, Bolser DG, Bonsell C, Dichiera AM, Hou Z, et al. Does functional redundancy affect ecological stability and resilience? a review and meta-analysis. Ecosphere. 2020;11:e03184.

    Article  Google Scholar 

  81. Van Rooij P, Martel A, Haesebrouck F, Pasmans F. Amphibian chytridiomycosis: a review with focus on fungus-host interactions. Vet Res. 2015. https://doi.org/10.1186/s13567-015-0266-0.

    Article  Google Scholar 

  82. Chen MY, Alexiev A, McKenzie VJ. Bacterial biofilm thickness and fungal inhibitory bacterial richness both prevent establishment of the amphibian fungal pathogen Batrachochytrium dendrobatidis. Appl Environ Microbiol. 2022;88:e01604-e1621.

    Article  CAS  Google Scholar 

  83. Keiser CN, Wantman T, Rebollar EA, Harris RN. Tadpole body size and behaviour alter the social acquisition of a defensive bacterial symbiont. R Soc Open Sci. 2019;6:191080.

    Article  Google Scholar 

  84. Muletz Wolz C, Myers J, Domangue R, Herrick J, Harris R. Soil bioaugmentation with amphibian cutaneous bacteria protects amphibian hosts from infection by Batrachochytrium dendrobatidis. Biol Conserv. 2012;152:119–26.

    Article  Google Scholar 

  85. Woodhams DC, LaBumbard BC, Barnhart KL, Becker MH, Bletz MC, Escobar LA, et al. Prodigiosin, violacein, and volatile Organic compounds produced by widespread cutaneous bacteria of amphibians can inhibit two Batrachochytrium fungal pathogens. Microb Ecol. 2018;75:1049–62.

    Article  CAS  Google Scholar 

  86. 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.

    Article  Google Scholar 

  87. Brucker R, Baylor C, Walters R, Lauer A, Harris R, Minbiole K. The Identification of 2,4-diacetylphloroglucinol as an antifungal metabolite produced by cutaneous bacteria of the salamander Plethodon cinereus. J Chem Ecol. 2008;34:39–43.

    Article  CAS  Google Scholar 

  88. Brucker RM, Harris RN, Schwantes CR, Gallaher TN, Flaherty DC, Lam BA, et al. Amphibian chemical defense: antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus. J Chem Ecol. 2008;34:1422–9.

    Article  CAS  Google Scholar 

  89. Rebollar EA, Gutiérrez-Preciado A, Noecker C, Eng A, Hughey MC, Medina D, et al. The skin microbiome of the Neotropical Frog Craugastor fitzingeri: inferring potential bacterial-host-pathogen interactions from metagenomic data. Front Microbiol. 2018;9:466.

    Article  Google Scholar 

  90. SahebKashaf S, Proctor DM, Deming C, Saary P, Hölzer M, et al. Integrating cultivation and metagenomics for a multi-kingdom view of skin microbiome diversity and functions. Nat Microbiol. 2022;7:169–79.

    Article  CAS  Google Scholar 

  91. Belden LK, Hughey MC, Rebollar EA, Umile TP, Loftus SC, Burzynski EA, et al. Panamanian frog species host unique skin bacterial communities. Front Microbiol. 2015. https://doi.org/10.3389/fmicb.2015.01171.

    Article  Google Scholar 

  92. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.

    Article  CAS  Google Scholar 

  93. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–6.

    Article  CAS  Google Scholar 

  94. Jari Oksanen, F. Guillaume Blanchet, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, et al. vegan: community ecology package. R package version 2.5–7. 2020. https://CRAN.R-project.org/package=vegan.

  95. Mandal S, Van Treuren W, White RA, Eggesbø M, Knight R, Peddada SD. Analysis of composition of microbiomes: a novel method for studying microbial composition. Microb Ecol Health Dis. 2015;26:27663.

    Google Scholar 

  96. R Core Team. R: a language and environment for statistical computing. R foundation for statistical computing, Vienna, Austria. 2021. URL https://www.R-project.org/.

  97. Bisanz J. qiime2R: importing QIIME2 artifacts and associated data into R sessions. Jordan E Bisanz. 2018. https://github.com/jbisanz/qiime2R.

  98. Wickham H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, 2016. https://ggplot2.tidyverse.org.

  99. Teunisse HGM. Fantaxtic: fantaxtic plots from Phyloseq data. R package version 0.1.0. 2018. https://github.com/gmteunisse/Fantaxtic.

  100. Gehlenborg N. UpSetR: a more scalable alternative to Venn and Euler diagrams for visualizing intersecting sets. R package version 1.4.0. 2019. https://CRAN.R-project.org/package=UpSetR.

  101. Mills BR. MetBrewer: color palette package in R inspired by works at the Metropolitan Museum of Art in New York. 2022. https://github.com/BlakeRMills/MetBrewer.

  102. Borcard D, Gillet F, Legendre P. Numerical ecology with R. 2011.

Download references

Acknowledgements

Emanuel. Martinez-Ugalde is a doctoral student from the Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional Autónoma de México (UNAM), and has received a Consejo Nacional de Ciencia y Tecnología (CONACyT) fellowship (781100). We thank Paola Palacios, Juan Pablo Pérez, Honoria Rueda-González, Inés Rosas-Chávez, Luis Zárate, and Bienes Comunales of Santiago Tlazala for the field assistance and Dr. Rafael Peña Miller for the assistance in the preparation of field data loggers housing. We thank Molly Bletz and Doug Woodhams for providing barcoded primers used during library construction.

Funding

This work was supported by PAPIIT/UNAM Grant No. IA201419 granted to ER.

Author information

Authors and Affiliations

  1. Centro de Ciencias Genómicas, Universidad Nacional Autónoma de México, Cuernavaca, Mexico

    Emanuel Martínez-Ugalde, Sara Lucia Anaya-Morales & Eria A. Rebollar

  2. Instituto de Ciencias Agropecuarias y Rurales, Universidad Autónoma del Estado de México, Toluca, Mexico

    Víctor Ávila-Akerberg

  3. Facultad de Ciencias, Universidad Nacional Autónoma de México, Mexico City, Mexico

    Tanya M. González Martínez, Montserrat Vázquez Trejo & Dalia Zavala Hernández

  4. Department of Biology, University of Mississippi, Oxford, MS, USA

    Sara Lucia Anaya-Morales

Authors
  1. Emanuel Martínez-Ugalde
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Víctor Ávila-Akerberg
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Tanya M. González Martínez
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Montserrat Vázquez Trejo
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Dalia Zavala Hernández
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Sara Lucia Anaya-Morales
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Eria A. Rebollar
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

ER conceptualized and designed the study. ER and VAA secured the funding. All authors contributed to sample collection. EM-U performed the molecular work and bioinformatic analysis. EM-U and ER wrote the manuscript. All authors participated in the improvement of the manuscript.

Corresponding author

Correspondence to Eria A. Rebollar.

Ethics declarations

Ethics approval and consent to participate

Our research was approved by the ethical standards of Universidad Nacional Autónoma de México, additionally capture and sampling of A. altamirani was approved by Subsecretaría de Gestión para la Protección Ambiental under the permit number: SGPA/DGVS/5673/19.

Consent for publication

Not applicable.

Competing interests

The authors declare 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.

Supplementary Figures.

Additional file 2.

Supplementary Tables.

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/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Martínez-Ugalde, E., Ávila-Akerberg, V., González Martínez, T.M. et al. The skin microbiota of the axolotl Ambystoma altamirani is highly influenced by metamorphosis and seasonality but not by pathogen infection. anim microbiome 4, 63 (2022). https://doi.org/10.1186/s42523-022-00215-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42523-022-00215-7

Keywords

Animal Microbiome

ISSN: 2524-4671