We sampled a total of 279 A. altamirani individuals (85 metamorphic and 194 non 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, 13 samples were discarded, and 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
List of samples collected of each sample type by sampling period.
|
Metamorphic
|
Non metamorphic
|
Sediment
|
Water
|
Total Samples (N)
|
Summer
(July 2019)
|
25
|
41
|
19
|
20
|
105
|
Autumn
(October 2019)
|
28
|
29
|
20
|
20
|
97
|
Winter
(January 2020)
|
9
|
66
|
20
|
20
|
115
|
Spring
(April 2020)
|
23
|
58
|
20
|
20
|
121
|
Total Samples (N)
|
85
|
194
|
79
|
80
|
438
|
Number of samples from the skin of A. altamirani individuals (metamorphic and non metamorphic) and environmental samples (sediment and water). Numbers in bold indicate total number of samples collected for each sample type or season. |
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 (Figure 1A). Sediment and water samples were the samples with highest numbers of unique ASVs (20,031 and 9,902 respectively), while metamorphic and non metamorphic samples had 8,916, and 6,650 unique ASVs respectively. Interestingly only 677 ASVs were shared between metamorphic and non metamorphic salamanders.
Taxonomic diversity 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 non metamorphic samples respectively, and 14.6% and 40.8% of sediment and water respectively (Additional file 1, Supplementary Fig. 1). For the axolotl samples we found that Chitinophagaceae and Pseudomonadaceae varied in relative abundance according to host metamorphic status, being Chitinophagaceae (metamorphic 2.7% / non metamorphic 27%) and Pseudomonadaceae more abundant in metamorphic samples (metamorphic 18.1% / non metamorphic 6.4%).
Bacterial alpha diversity was significantly different between sample types (metamorphic, non metamorphic, sediment and water) with Observed ASVs (Kruskal-Wallis (KW), \(\chi 2\) = 278.46, p-value < 0.001), Shannon index (KW, \(\chi 2\) = 276.28, p-value < 0.001) and Faith´s phylogenetic diversity (PD) (KW, \(\chi 2\) = 286.91, p-value < 0.001) (Fig. 1B). Post hoc pairwise comparisons for each alpha diversity index showed significant differences among all sample types (Additional file 2, Supplementary Table 1) except for metamorphic salamanders and water in observed ASVs (Wilcoxon, p-value = 0.48) and Shannon diversity index (Wilcoxon, p-value = 0.66). Sediment samples showed the highest alpha diversity values while non metamorphic salamanders always had the lowest values.
Bacterial beta diversity based on the weighted UniFrac distance matrix varied significantly among sample types (PERMANOVA, pseudo-F = 64.76, p-value < 0.001) (Fig. 1C, Additional file 2, Supplementary Table 2). Sample dispersion also significantly differed among sample types (PERMUTEST, F = 34.5, p-value = 0.001) (Fig. 1D, Additional file 2, Supplementary Table 3) with non metamorphic samples showing the greatest dispersion.
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 non 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 non 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 non metamorphic samples with Burkholderiaceae, Chitinophagaceae being the families with higher W values.
To identify the core microbiota of each A. altamirani metamorphic stage we calculated the bacterial ASVs that were highly prevalent in either metamorphic or non metamorphic individuals but absent from environmental samples. We found that seven bacterial ASVs represent the bacterial core of metamorphic axolotls accounting for a cumulative relative abundance of 17.4%. Meanwhile, four bacterial ASVs represent the bacterial core of non metamorphic axolotls accounting for 52.26% of the relative abundance (Table 2). Interestingly, two of the four core ASVs present in non metamorphic samples accounted for 25.9% and 20.75% of the relative abundance in these samples and belonged to the bacterial families Chitinophagaceae and Burkholderiaceae respectively. In addition, we identified that two ASVs from the Pseudomonadaceae family were part of the core of both metamorphic and non metamorphic samples.
Table 2
Amplicon sequence variants (ASVs) defining the core skin microbiota of metamorphic and non metamorphic A. altamirani.
ASV ID
|
Sample Type
|
Family
|
Relative abundance
|
Persistence
|
9936daae333af6e517a9deb4b9e18ffa
|
Metamorphic
|
Pseudomonadaceae
|
13.83
|
95.12%
|
6d0c9d0395e6a2a7667eb0b07c17a275
|
Metamorphic
|
Burkholderiaceae
|
1.4
|
97.56%
|
17d60505100c3cf44d4f9fad620d1636
|
Metamorphic
|
Pseudomonadaceae
|
0.74
|
93.90%
|
be8eb25874b4202cf98050dbadeeb7ce
|
Metamorphic
|
Burkholderiaceae
|
0.33
|
93.90%
|
8fd9eab61a0f63db6cbb201ce66b484b
|
Metamorphic
|
Burkholderiaceae
|
0.36
|
82.93%
|
6ddbdbae5830daa9d7d08270857b02b8
|
Metamorphic
|
Rhizobiaceae
|
0.27
|
82.93%
|
d21b6c5e4a5c9d5f40bda6b543e4208f
|
Metamorphic
|
Xanthomonadaceae
|
0.21
|
82.93%
|
3c28f0caf9183357de05d1882a943f8e
|
Non metamorphic
|
Chitinophagaceae
|
25.09
|
96.84%
|
ed5a79897d0f82525c3854759d384c26
|
Non metamorphic
|
Burkholderiaceae
|
20.75
|
98.42%
|
9936daae333af6e517a9deb4b9e18ffa
|
Non metamorphic
|
Pseudomonadaceae
|
6.03
|
83.16%
|
17d60505100c3cf44d4f9fad620d1636
|
Non metamorphic
|
Pseudomonadaceae
|
0.39
|
87.37%
|
ASVs were considered part of the skin bacterial core if they were present in ≥ 80% of the skin samples for metamorphic or non metamorphic axolotls.
Seasonality And Location Differentially Influence Skin Bacterial Diversity In Metamorphic And Non 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) and sampling locations (MANOVA, Wilks = 0.0009, p-value < 0.001). While all physicochemical variables varied across seasons, dissolved oxygen was the only variable that did not vary between sampling locations (Additional file 2, Supplementary Table 4).
Alpha PD of metamorphic axolotls varied significantly across seasons (KW, \(\chi 2\) = 13.69, p-value = 0.003) (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, Supplementary Table 5). In contrast, PD of non metamorphic A. altamirani (Fig. 3B) did not differ across consecutive seasons (KW, \(\chi 2\) = 0.21, p-value = 0.97) (Additional file 2, Supplementary Table 5).
Additionally, we found that seasonality significantly influenced skin bacterial beta diversity of metamorphic (PERMANOVA, pseudo-F = 12.37, p-value < 0.001) (Fig. 3D, Additional file 2, Supplementary Table 6) and non metamorphic (PERMANOVA, pseudo-F = 15.69, p-value < 0.001) (Fig. 3E, Additional file 2, Supplementary Table 7) axolotls. Specifically, pairwise PERMANOVAs showed that metamorphic samples differed between winter-spring seasons (PERMANOVA, pseudo-F = 14.92, p-value = 0.001), while non metamorphic skin microbiota differed between autumn-winter (PERMANOVA, pseudo-F = 13.47, p-value < 0.001) and winter-spring seasons (PERMANOVA, pseudo-F = 12.61, p-value < 0.001).
Three bacterial families were identified by ANCOM as differentially abundant in metamorphic samples between winter-spring seasons (Fig. 4). In the case of non metamorphic individuals, ANCOM identified three bacterial families that were differentially abundant between autumn-winter and seven families as differentially abundant between winter-spring (Fig. 4). Psuedomonadaceae, Aquaspirillaceae and Shewanellaceae were shared between metamorphic and non metamorphic axolotls during winter and spring seasons. However, Pseudomonadaceae was differentially enriched according to host metamorphic status, being more abundant in metamorphic axolotls during spring and more abundant in non metamorphic axolotls during winter.
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), however post hoc paired test showed that PD only differed significantly between sites 2 and 3 (Additional file 1, Supplementary Fig. 2A). Bacterial PD of non metamorphic samples also varied significantly between sampling locations (KW, \(\chi 2\) = 40.9, p-value = 6.71e-9). 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, Supplementary Fig. 2C, Additional file 2, Supplementary Table 8). Beta diversity of the skin bacterial communities was also influenced by sampling location in metamorphic (PERMANOVA, pseudo-F = 2.71, p-value = 0.006) and in non metamorphic samples (PERMANOVA, pseudo-F = 31.34, p-value = 0.001) (Additional file 1, Supplementary Fig. 2B, D). Pairwise comparisons showed that bacterial beta diversity only differed between sites 2 and 3 in metamorphic axolotls (Additional file 2, Supplementary Table 9), while bacterial beta diversity differed between all sampling locations for non metamorphic samples (Additional file 2, Supplementary Table 10).
Biotic and abiotic factors influence the skin bacterial community structure of A. altamirani
Our results showed that beta bacterial diversity of A. altamirani skin is influenced by season 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 following biotic and abiotic factors that resulted informative 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 (Table 3, Additional file 1, Supplementary Fig. 3). Permutational analyses (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 4).
Table 3
Variance explained of each canonical axis calculated by dbRDA.
|
F
|
p-value
|
Variance explained
|
Cumulative variance
|
CAP1
|
63.3875
|
0.001
|
0.173
|
0.173
|
CAP2
|
18.4679
|
0.001
|
0.050
|
0.223
|
CAP3
|
9.8065
|
0.001
|
0.026
|
0.250
|
CAP4
|
5.1977
|
0.006
|
0.014
|
0.264
|
CAP5
|
2.5146
|
0.123
|
0.006
|
0.271
|
CAP6
|
1.5675
|
0.387
|
0.004
|
0.275
|
CAP7
|
1.1402
|
0.563
|
0.003
|
0.279
|
CAP8
|
0.7233
|
0.764
|
0.001
|
0.281
|
Columns indicate: F statistic, p-values, variance explaines by each canonical axis, and the cumulative variance calculated by the Permutational like ANOVA. Numbers in bold indicate significant p-values and the cumulative variance for each statistically significant canonical axis.
Table 4
Permutational like ANOVA results of each variable introduced in the dbRDA regression model.
|
F
|
p-value
|
Developmental Stage
|
39.121
|
0.001
|
Delta Temperature
|
19.889
|
0.001
|
pH
|
15.854
|
0.001
|
Mean Temperature
|
12.053
|
0.001
|
Elevation
|
5.334
|
0.001
|
Dissolved Oxygen
|
4.478
|
0.002
|
Conductivity
|
3.470
|
0.005
|
Weight
|
2.604
|
0.018
|
Columns indicate: F statistic, p-values calculated by Permutational like ANOVA for each variable. Numbers in bold indicate significant p-value. |
Skin bacterial diversity of A. altamirani is not influenced by Bd infection status but specific bacterial taxa abundance correlate with infection intensity
Alpha PD did not differ between infected and non-infected samples in both metamorphic (KW, = 0.09, p-value = 0.76) (Figure 5A) and non metamorphic (KW, = 0.51, p-value = 0.47) A. altamirani samples (Figure 5C). 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) (Figure 5B) and non metamorphic salamanders (PERMANOVA, pseudo-F = 2.45, p-value = 0.08) (Figure 5D).
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 non metamorphic samples respectively, correlated with pathogen infection loads (Additional file 1, Supplementary Fig. 4). Specifically, 116 (out of 139) and 52 (out of 128) bacterial ASVs had positive correlations with pathogen infection loads in metamorphic and non metamorphic samples respectively.
All the ASVs that correlated with pathogen load in metamorphic samples had low relative abundances (0.001–0.67%) (Additional file 1, Supplementary Fig. 4A), while in non metamorphic samples the correlated ASVs ranged from 0.001–28.5% (Additional file 1, Supplementary Fig. 4B). Twelve ASVs with significant correlations were shared between metamorphic and non metamorphic samples and six of these showed different types of correlations with Bd loads according to axolotl metamorphic status. Interestingly, six of the twelve shared ASVs, had different correlations (negative vs positive) according to host metamorphic status (Additional file 2, Supplementary Table 11).