An average richness of 591 OTUs (range 3 – 3,601) was recovered per skin swab. The most abundant bacterial phyla shared by all three salamander genera were Proteobacteria, Verrucomicrobia, Acidobacteria, Actinobacteria, and Bacteroidetes (Fig. 1). Shapiro-Wilks tests indicated that data were not normally distributed within salamander genera (p < 0.05 for all tests). Both measures of α-diversity differed among salamander genera (Shannon index: Kruskal-Wallis χ2 = 22.752, p < 0.001 ; inverse-Simpson index: Kruskal-Wallis χ2 = 32.472, p < 0.001).
At all nine sites and for all three ecoregions, multivariate dispersion of bacterial skin assemblages (betadisper) did not differ among host genera (p > 0.05 for all tests). When the entire dataset was included, the average structure of bacterial assemblages differed only by salamander genus (PERMANOVA p < 0.001, R2 =0.0407, Table 2). There was a significant interaction effect between genus and ecoregion (p < 0.001), indicating that host-specific differences in bacterial assemblage structure vary regionally. Pairwise comparisons indicated that skin assemblages were distinct for each genus (Bonferroni-adjusted p = 0.003 for each pair of host genera, Fig. 2).
Table 2
PERMANOVA results for main effects tests of the variables Genus, Ecoregion, and Season, and Genus-Ecoregion interaction. Terms were added sequentially (first to last). Abbreviations: df = Degrees of freedom; SS = Sum of Squares; MS = Mean Square; Pr(>F) = percent of permutations resulting in values greater than Pseudo-F.
Source
|
df
|
SS
|
MS
|
Pseudo-F
|
R2
|
Pr(>F)
|
Genus
|
2
|
4.179
|
2.08964
|
5.3134
|
0.04072
|
0.001 ***
|
|
Ecoregion
|
2
|
2.704
|
1.35189
|
3.4375
|
0.02635
|
0.645
|
|
Season
|
1
|
0.527
|
0.52677
|
1.3394
|
0.00513
|
0.033
|
|
Genus:Ecoregion
|
4
|
2.400
|
0.60007
|
1.5258
|
0.02339
|
0.001 ***
|
|
Residuals
|
236
|
92.814
|
0.39328
|
|
0.90441
|
|
|
Total
|
245
|
102.624
|
|
|
1.00000
|
|
|
PERMANOVAs, as implemented with function adonis, add factors sequentially during R2 calculation, but the sequence of terms matters when factors are correlated and linear dependency can be assumed, as in the case of host genus and site. Distance-based redundancy analyses (function dbrda) were required to find marginal effects for these two factors [34]. After controlling for the effect of host genus, the marginal effect of site was still significant, and the inverse was also true (dbrda; Genus: F2, 235 = 4.7661, p ≤ 0.001; Site: F8, 235 = 2.4884, p ≤ 0.001). These analyses were repeated for each ecoregion separately. For all three ecoregions, host genus explained a greater proportion of skin assemblage variation than site, whether permutations were constrained within sites (adonis; Cumberland Plateau: Genus: F2, 65 = 2.1565, R2 = 0.05594, p ≤ 0.001; Site: F2, 65 = 1.6793, R2 = 0.04356, p = 0.072; Ridge and Valley: Genus: F2, 73 = 2.6833, R2 = 0.06177, p ≤ 0.001; Site: F2, 73 = 2.4127, R2 = 0.05554, p = 0.981; Blue Ridge: Genus: F2, 82 = 3.9375, R2 = 0.07914, p ≤ 0.001; Site: F2, 82 = 2.2490, R2 = 0.04520, p = 0.035; Table 3) or not constrained within sites (Cumberland Plateau: Genus: F2, 65 = 2.1565, R2 = 0.05594, p ≤ 0.001; Site: F2, 65 = 1.6793, R2 = 0.04356, p ≤ 0.001; Ridge and Valley: Genus: F2, 73 = 2.6833, R2 = 0.06177, p ≤ 0.001; Site: F2, 73 = 2.4127, R2 = 0.05554, p ≤ 0.001; Blue Ridge: Genus: F2, 82 = 3.9375, R2 = 0.07914, p ≤ 0.001; Site: F2, 82 = 2.2490, R2 = 0.04520, p ≤ 0.001).
Table 3
PERMANOVA results for main effects tests of the variables Genus, Site, and Genus-Site interaction within the (a) Cumberland Plateau, (b) Ridge and Valley, and (c) Blue Ridge ecoregions when permutations were constrained within sites. Terms added sequentially (first to last). Abbreviations: Ge = Genus; Res = Residuals; df = Degrees of freedom; SS = Sum of Squares; MS = Mean Square; Pr(>F) = percent of permutations resulting in values greater than Pseudo-F.
(a) Cumberland Plateau
Source
|
df
|
SS
|
MS
|
Pseudo-F
|
R2
|
Pr(>F)
|
Ge
|
2
|
1.6841
|
0.84205
|
2.1565
|
0.05594
|
0.001 ***
|
|
Site
|
2
|
1.3114
|
0.65572
|
1.6793
|
0.04356
|
0.072
|
|
Ge x Site
|
4
|
1.7283
|
0.43207
|
1.1065
|
0.05741
|
0.070
|
|
Res
|
65
|
25.3805
|
0.39047
|
|
0.84309
|
|
|
Total
|
73
|
30.1043
|
|
|
1.00000
|
|
|
(b) Ridge and Valley
Source
|
df
|
SS
|
MS
|
Pseudo-F
|
R2
|
Pr(>F)
|
Ge
|
2
|
2.0670
|
1.03336
|
2.6833
|
0.06177
|
0.001 ***
|
|
Site
|
2
|
1.8580
|
0.92918
|
2.4127
|
0.05554
|
0.981
|
|
Ge x Site
|
3
|
1.4210
|
0.47374
|
1.2301
|
0.04248
|
0.015 *
|
|
Res
|
73
|
28.1113
|
0.38511
|
|
0.84022
|
|
|
Total
|
80
|
33.4600
|
|
|
1.00000
|
|
|
(c) Blue Ridge
Source
|
df
|
SS
|
MS
|
Pseudo-F
|
R2
|
Pr(>F)
|
Ge
|
2
|
2.7830
|
1.39147
|
3.9375
|
0.07914
|
0.001 ***
|
|
Site
|
2
|
1.5900
|
0.79478
|
2.2490
|
0.04520
|
0.035 *
|
|
Ge x Site
|
4
|
1.8130
|
0.45323
|
1.2825
|
0.05156
|
0.011 *
|
|
Res
|
82
|
28.978
|
0.35339
|
|
0.82410
|
|
|
Total
|
90
|
35.163
|
|
|
1.00000
|
|
|
No significant distance-decay relationship was found for any measure of β-diversity (GLM; SOR: z = 0.002, p > 0.05; SIM: z = 0.05, p > 0.05, SNE: z = -0.082, p > 0.05; Fig. 3), and there was no significant difference in the rate of decay between the three genera for all three measures of β-diversity (ANOVA; SOR: p > 0.05, SIM: p > 0.05, SNE: p >0.05).
|
For analyses of climate data, total β-diversity (SOR), measured as multivariate dispersion, was significantly different across precipitation seasonality values for all three genera (betadisper; Plethodon: F3, 70 = 7.6037, p ≤ 0.01, Eurycea: F3, 61 = 13.899, p ≤ 0.01, Desmognathus: F3, 100 = 8.4202, p ≤ 0.01). Additionally, β-diversity measured as multivariate dispersion (betadispr) varied across values for mean diurnal range for both Plethodon and Eurycea (Plethodon: F6, 67 = 5.4202, p ≤ 0.01, Eurycea: F7, 57 = 2.2618, p > 0.01), but not for Desmognathus microbial assemblages (F7, 96 =2.1205, p > 0.01). Interestingly, all other climatic variables were only significantly predictive of Plethodon salamander microbial assemblages (isothermality: F4, 69 = 4.4332, p ≤ 0.01, precipitation of the driest quarter: F7, 66 = 3.7282, p ≤ 0.01, maximum temperature of warmest month: F6, 67 = 3.3851, p ≤ 0.01), but not for Eurycea (isothermality: F5, 59 = 2.3785, p > 0.01, precipitation of the driest quarter: F8, 56 = 1.8114, p > 0.01, maximum temperature of warmest month: F6, 58 = 4.0724, p ≤ 0.01) or Desmognathus (isothermality: F5, 98 = 2.4632, p > 0.01, precipitation of the driest quarter: F8, 95 = 1.3033, p > 0.01, maximum temperature of warmest month: F6, 97 = 1.9748, p > 0.01).