Bacterial Skin Assemblages of Sympatric Salamanders Are Primarily Shaped by Host Genus

Bacterial assemblages on the skins of amphibians are known to influence pathogen resistance and other important physiological functions in the host. Host-specific factors and the environment play significant roles in structuring skin assemblages. This study used high-throughput 16S rRNA sequencing and multivariate analyses to examine differences in skin-bacterial assemblages from 246 salamanders belonging to three genera in the lungless family Plethodontidae along multiple spatial gradients. Composition and α- and β-diversity of bacterial assemblages were defined, indicator species were identified for each host group, and the relative influences of host- versus environment-specific ecological factors were evaluated. At the broadest spatial scale, host genus, host species, and sampling site were predictive of skin assemblage structure, but host genus and species were more influential after controlling for the marginal effects of site, as well as nestedness of site. Furthermore, assemblage similarity within each host genus did not change with increasing geographic distance. At the smallest spatial scale, site-specific climate analyses revealed different relationships to climatic variables for each of the three genera, and these relationships were determined by host ecomode. Variation in bacterial assemblages of terrestrial hosts correlated with landscape-level climatic variability, and this pattern decayed with increasing water dependence of the host. Results from this study highlight host-specific considerations for researchers studying wildlife diseases in co-occurring, yet ecologically divergent, species.


Introduction
The vertebrate skin-microbial community functions as an extension of the innate immune system, providing physical and chemical barriers to pathogens [1,2]. The amphibian skin microbiome is known to protect its host against pathogens such as Batrachochytrium dendrobatidis (Bd) and B. salamandrivorans (Bsal) through the induction of antimicrobial peptides and the production of antifungal metabolites [3][4][5][6]. Defining the composition and assembly processes affecting the amphibian cutaneous microbiome is important for understanding its role in pathogen defense and other physiological functions in the host [4].
Discrepancies exist among previous studies regarding the relative importance of the metazoan host versus the environment in shaping microbial skin assemblages. Past work using 16S rRNA metabarcoding has revealed that skin assemblage composition differs according to host taxon, life history stage, geographic region, microhabitat type, and even body regions within amphibian individuals [7][8][9][10]. Landscape-level ecological factors known to shape the amphibian skin-bacterial assemblage include humidity, temperature [9], elevation [11], and microbial reservoirs in soil [12,13]. Environmental dispersal influences skin bacterial assemblages of plethodontid salamanders, but hostspecific constituents have also been reported [12][13][14][15]. Significant associations between the skin-bacterial assemblage and amphibian host taxon have been found at the order [16], family [8,9], genus [17,18], species, and even subspecies ranks [19]. Conversely, separate studies comparing terrestrial species of lungless salamanders found that selection of skin bacteria is driven primarily by the host's environment and, secondarily, by its taxonomy [15,20].
Because salamander life history traits inextricably bind them to distinct microhabitats, and thus, potentially unique bacterial reservoirs [21], it is possible that host-specific variation may be more difficult to resolve at lower taxonomic levels due to distinct environmental reservoirs for microbial dispersal. Distance decay is defined as a decrease in assemblage similarity with increasing geographic or environmental distance. Analyses of distance decay, when coupled with an analysis of spatial and climatic data, allow for explicit tests of processes acting on skin assemblage diversity at both the local and regional scales to help disentangle effects of ecological factors and elucidate patterns in community assembly processes [22][23][24]. This study takes advantage of the overlapping geographic distribution, close proximity of salamander microhabitats at a particular site, and evolutionary relatedness of the genera Desmognathus, Eurycea, and Plethodon in order to assess potential factors structuring skin-bacterial assemblages of lungless salamanders of the eastern USA.
The southeastern USA hosts the most biodiverse assemblage of salamander species on planet Earth [25,26]. At least 47 species are found within the state of Tennessee [27]. Tennessee consists of eight physiographic ecoregions (Environmental Protection Agency level III) with different geology, vegetation, and hydrology [28,29]. Three genera of salamanders in the lungless family Plethodontidae (Desmognathus, Eurycea, Plethodon) are found ubiquitously across the three easternmost Tennessee ecoregions (Cumberland Plateau, Ridge and Valley, Blue Ridge Mountains). These genera have distinct life history characteristics and divergent ecomodes. The genus Desmognathus is mostly aquatic and usually found in headwater and secondary order streams [27]. The genus Eurycea is semiaquatic and often found directly adjacent to streams beneath rocks and logs on the riverbank [27]. The genus Plethodon is completely terrestrial and often found just a few meters away from the stream bank beneath moist decaying logs [27]. The sympatric distribution and life histories of these salamander genera make them an ideal study system to understand the relative importance of host characteristics and/or environmental influence in structuring their skin microbiomes.
The objectives of this project were to use high-throughput 16S rRNA metabarcoding and sequencing to (1) characterize the bacterial skin communities of three co-occurring but ecologically divergent genera of plethodontid salamanders across three ecoregions of Tennessee, USA, (2) measure turnover in OTU assemblage structure across multiple gradients (i.e., physiographic ecoregions, climates, salamander genera, and species) with distance decay analyses, and (3) partition variation in assemblage structure according to geography and host genus. We hypothesized that variation in bacterial assemblages would be explained primarily by host genus, but that degree of variation would increase with increasing geographic distance of a single host.

Sample Collection
Adult salamanders belonging to the family Plethodontidae were targeted for skin swab sampling, including the genera Plethodon (n = 74), Eurycea (n = 63), and Desmognathus (n = 109). Species were selected based on similar ecological traits and syntopy. The three desmognathine species D. fuscus, D. monticola, and D. santeetlah are primarily aquatic, stream-dwelling species which can occupy a common microhabitat at sites where they co-occur [27,30]. The salamanders Eurycea cirrigera and E. wilderae are semiaquatic, riparian zone species belonging to the two-lined salamander complex [27]. The congeneric long-tailed salamander E. longicauda is more troglophilic, but its life history traits overlap with that of the two-lined species [31]. Lastly, the genus Plethodon typifies the fully terrestrial ecomode. The species used in this study, P. glutinosus and P. jordani, also occupy comparable ecological niches within their respective geographic ranges [27]. See Appendix A for detailed species identification and genotyping protocols.
Nine protected sites across three physiographic ecoregions of Tennessee were sampled-the Cumberland Plateau, Ridge and Valley, and Blue Ridge Mountains. Five of the salamander species included in this study are ubiquitous across all three physiographic ecoregions, but three species (D. santeetlah, E. wilderae, and P. jordani) are restricted to the Blue Ridge Mountains [32]. Members of all three salamander genera were sampled at each site, except for one site within the Ridge and Valley ecoregion where the genus Plethodon does not occur. To rule out potential impacts of seasonal variation on microbial assemblages, samples were collected during concurrent weeks in both 2016 and 2017, specifically, May 13-21 and July 21-30.
Each salamander was captured by hand, placed into a clean plastic bag with one corner removed for water drainage, and rinsed for 60 s with 2-h-autoclaved Millipore water to remove transient microbes. A sterile rayon-tipped swab was rolled across all body surfaces for 15 strokes and stored in an empty, sterile 2-mL tube to be used for characterization of cutaneous bacterial assemblages via 16S rRNA metabarcoding and high-throughput sequencing. A full-process negative control sample was collected at field sites by rinsing and swabbing the inside of an empty plastic bag (as above), so that contaminating bacteria could be sequenced and removed during bioinformatics processing of sequence data, described below. Salamanders were released the same day at the same location where they were captured. Collection permits were obtained, and institutional guidelines for ethical use of animals in research were observed. All samples were promptly stored in a − 12° C portable freezer, transported to the laboratory, and stored at − 20° C until processing.

High-Throughput Sequencing and Bioinformatics
Total genomic DNA was extracted from swab tips using the QIAamp PowerFecal Kit (QIAGEN, Inc., Hilden, Germany). For each DNA extraction, a single sterile swab was processed as above to serve as a negative control (DNA extraction blank). The V4 region of the 16S rRNA gene was PCR-amplified with MC Lab I-5™ 2X Hi-Fi Master Mix (Molecular Cloning Laboratories, South San Francisco, CA) and primers 515F and 806R containing i5 and i7 Illumina flow cell adapter sequences [as in 33]. During amplicon PCR, autoclaved Millipore water was used in place of DNA in a single well of the PCR plate, to serve as an additional negative control sample. Amplicons were dual-indexed with Nextera XT Index Kit Set B, normalized, pooled, and sequenced on the Illumina MiSeq platform (Illumina, Inc., San Diego, CA) using 2 × 250 bp paired-end reads.
To characterize bacterial assemblage composition, mothur v. 1.40.3 was used to screen sequences according to standard filtering parameters, cluster sequences at 97% identity into operational taxonomic units (OTUs), assign taxonomy to OTUs, and remove non-target and contaminant sequences (n = 814 OTUs, roughly 4.7% of the total) [33]. OTUs were used in downstream analyses rather than ASVs (amplicon sequence variants), as negligible differences have been noted between approaches in ecological studies [34], and ASVs have been shown to artificially inflate diversity estimates due to splitting single genomes into separate clusters [35]. Alternatively, these data could be reassessed using ASVs to reflect a more refined taxonomic level not subject to a distance-based threshold. The taxonomy of each OTU was determined by comparison to the SILVA v.132 reference database [36]. Final sequence counts for samples were compared after quality control and found to be significantly different across libraries (Kruskal-Wallis χ 2 = 48.462, df = 3, p < 0.001). Therefore, the final dataset was subsampled at 2009 sequences per sample to account for library size differences [37]. For all downstream statistical analyses described below, OTUs occurring fewer than twice in the dataset were excluded, to avoid undue weighting of these OTUs in the presence/absence matrices of community structure, such as the one used in this study (discussed in the "Statistical Analyses" section). The final OTU × sample matrix consisted of 16,381 OTUs and 246 skin swab samples from salamanders belonging to the genera Desmognathus (n = 109), Eurycea (n = 63), and Plethodon (n = 74) ( Table 1). Note that, according to Sanger sequencing results for a single mitochondrial locus (see Appendix A), we found seven D. santeetlah individuals at a site within the Ridge and Valley where the species has not been documented. However, sequencing of additional loci could lead to a new species assignment for these individuals.

Statistical Analyses
We followed the statistical framework defined by Anderson et al. [38] and used their dual definitions of β-diversity to (1) measure turnover in OTU assemblage structure across three spatial gradients (i.e., physiographic ecoregions, salamander genera, and salamander species) and (2) partition variation in assemblage structure according to geography and host genus. All analyses were performed in R ver. 4.0.3 [39]. Shannon and inverse Simpson indices of α-diversity were calculated for salamander genera and species, and Shapiro-Wilk tests for homogeneity of variance were conducted on these α-diversity metrics. Kruskal-Wallis tests were used to discern differences in α-diversity among groups. Sørensen dissimilarity values of total β-diversity were calculated from a matrix of OTU presence/absence data [as in 40], and betadisper tests of multivariate dispersion in package vegan were conducted on Sørensen indices. Bray-Curtis dissimilarity values were generated from OTU relative abundance data. To determine if β-diversity was similar among co-occurring salamander genera, differences in dispersion were assessed at each of the nine field sites and both within and among the three Tennessee ecoregions. Differences in dispersion were compared using principal coordinates analysis (PCoA). To determine the effect of host genus, ecoregion, and site on skin microbiome structure, the adonis function in package vegan was used to conduct permutated multivariate analyses of variance (PERMANO-VAs) of both Sørensen indices and Bray-Curtis dissimilarities, with Bonferroni-corrected α = 0.0125. We assessed broad-scale patterns in the skin microbiome by running a PERMANOVA for all samples, and ecoregion-specific patterns using a second set of PERMANOVAs for each ecoregion individually. Permutations were constrained within sites to account for nestedness of sites within ecoregions (strata = site). Distance-based redundancy analyses were conducted post hoc to disentangle marginal effects of putatively correlated factors. To visualize differences in bacterial assemblage structure, Bray-Curtis dissimilarity values generated from OTU presence/absence data were plotted by sample type in two-dimensional nonmetric multidimensional scaling (NMDS) ordinations. Indicator species were identified for each host genus and species from OTU abundance data, using the R function indicators in the package indicspecies, as well as the mothur command indicator, at α = 0.05. The indicator value, reported as IndVal in mothur and square root transformed in R to sqrtIV, is the product of the positive predictive value (called "Component A") and the sensitivity value (called "Component B") [41,42]. For the purposes of this study, we interpreted the positive predictive value (Component A) as the estimated conditional probability that the sampled salamander belongs to a specific host group, given that the indicator was found there. We interpreted sensitivity (Component B, also called fidelity) as the estimated conditional probability of finding the indicator on a salamander host belonging to a specific host group. Within the indicators function, we set a threshold of 0.6 for Component A, and a threshold of 0.3 for Component B, sensitivity (fidelity). Indicators not meeting these parameters were excluded.
To assess the influence of geographic distance, we compared total β-diversity (SOR), turnover (SIM), and nestedness (SNE) components of Sørensen diversity across geographic distances, following the methods of Grisnik et al. [43]. Specifically, we calculated the geographic distances between sample points (calculated as Euclidian distances) using the dist function in the vegan package. Due to issues with singular fit in the mixed models, we removed nestedness of our data by averaging beta diversity metrics (SOR, SIM, SNE) by the dummy variable "site contrast." The site contrast variable described the pairwise site-level comparisons. We then used a generalized linear model (GLM; function glm), to compare average assemblage structure and geographic distance. The GLM assumed a binomial distribution, with average geographic distance and host genera set as fixed effects and a Bonferroni-adjusted p-value of 0.016. To determine if there was a difference in the rate of distance decay across the three genera, we used function anova (package car) with a type II sum-of-squares to account for unequal sample sizes across groups.
Prior to investigating climate-based effects on β-diversity for each salamander genus, we verified the appropriateness of using host genus as the study unit (rather than species) by re-running all of the above α-and β-diversity analyses without samples from the three range-restricted species P. jordani, E. wilderae, and D. santeetlah. Because results were unaffected by the exclusion of these host species, we proceeded with the following analysis with all samples included: We acquired climatic data for each site from the WorldClim database (http:// www. world clim. org) at 30-s resolution (1 km 2 ). These data represent a series of 19 bioclimatic variables derived from global temperature and precipitation grids [44]. We ran a correlation analysis and removed highly correlated variables (Pearson's correlation coefficient > 0.75), keeping the variables of hypothesized importance, resulting in five bioclimatic variables [45]. The five variables included in this analysis were precipitation seasonality (Bio15), mean diurnal temperature range (Bio2), isothermality (day-to-night temperature variability across seasons; Bio3), precipitation of driest quarter (Bio17), and maximum temperature of the warmest month (Bio5). We then extracted the values for each climatic variable using the extract function in the raster package [46]. To elucidate how climatic variables influenced β-diversity across sites, we compared variation in multivariate dispersion (betadisper) across the five selected climatic variables for each host genus. We used a Bonferroni-adjusted p-value of 0.01 to account for multiple comparisons.

Results
An average richness of 591 OTUs (range 3-3601) was recovered per skin swab. The most abundant bacterial orders shared by all three salamander genera were Betaproteobacteriales, Verrucomicrobiales, Micrococcales, Flavobacteriales, and Rhizobiales (Fig. 1). Shapiro-Wilk tests indicated that α-diversity data were not normally distributed within any of the three salamander genera or within five of the eight salamander species (p < 0.05 for all tests). However, Shannon indices of α-diversity were normally distributed within the 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 by salamander genus but not by ecoregion (PERMANOVA; Table S2, Appendix C; Fig. 2). These results held true whether OTU abundances were included in the β-diversity metric or not (i.e., Bray-Curtis dissimilarity and Sørensen index, respectively). However, there was a significant interaction effect between genus and ecoregion (p < 0.001; Table S2, Appendix C), indicating that host genus-specific differences in bacterial assemblage structure vary regionally. Pairwise comparisons indicated that skin assemblages were distinct for each pair of host genera (Bonferroni-adjusted p = 0.003).
PERMANOVA results at the host species level were consistent with those found at the genus level. The average structure of bacterial assemblages differed by salamander species but not by ecoregion (PERMANOVA; Table S3, Appendix C). The significant interaction effect between species and ecoregion (p < 0.002; Table S3, Appendix C) indicated that host species-specific differences in bacterial assemblage structure vary regionally. Pairwise comparisons indicated that skin assemblages did not differ between any pair of host species at Bonferroni-adjusted p = 0.003.
PERMANOVAs, as implemented with function adonis, add factors sequentially during R 2 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 [47]. After controlling for the effect of host genus, the marginal effect of site was still significant, and the inverse was also true (dbrda; genus, F 2, 235 = 4.7661, p ≤ 0.001; site, F 8, 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; Table S4, Appendix C) or not constrained within sites (Table S5, Appendix C). Similar distance-based redundancy analysis and PERMANOVA results were found at the host species level (dbrda; species,  Tables S6 and S7, Appendix C).
Indicator species analyses selected multiple unclassified strains from family Burkholderiaceae (phylum Proteobacteria) as highly significant indicator species for genus Desmognathus, as well as Eurycea longicauda (sq. rt. indicator values > 0.84; p-values = 0.005). For genus Eurycea, indicators belonged to the families Polyangiaceae (phylum Proteobacteria; sq. rt. indicator value = 0.57; p-value = 0.010) and Micromonosporaceae (phylum Actinobacteria; sq. rt. indicator value = 0.054; p-value = 0.005). For genus Plethodon, five of the seven OTUs with the highest sq. rt. indicator values (> 0.6) were soil-associated bacteria belonging to the class Acidobacteriia. See Appendix D for detailed results of indicator species analyses.

Discussion
With the impending North American invasion of the chytrid pathogen Bsal, which threatens rare and endemic salamander populations [48], researchers continue to uncover ecological assembly processes of the cutaneous microbiome. This study provides evidence that, within and among sites, sympatric species harbor distinct skin-bacterial assemblages, and this distinctiveness is due, in part, to their differing modes of interaction with the environment. Not surprisingly, both host genus and sampling site were identified as factors that are influential in structuring skin assemblages. At the broadest geospatial scale, when all samples from all regions were included, it was not possible to disentangle the relative influences of host genus and site on bacterial assemblage structure. Looking within each ecoregion, however, a clearer pattern emerged. At the regional level, there was a stronger effect of host genus on skin β-diversity compared to site. These results, consistent with the findings of Buttimer et al. [17], support the conclusion that the influence of the environment varies by host genus. In other words, host genus determines the degree to which the bacterial assemblage is shaped by the environment, and this effect seems to be scale-dependent, with the strength of the environmental signal increasing at smaller spatial scales.
We did not find a distance decay in similarity of salamander cutaneous microbial assemblages, suggesting a lack of dispersal limitation or species sorting driving assembly processes [22,23]. However, we found a significant interaction between host genus and ecoregion which, when coupled with the lack of distance decay, supports the hypothesis that the local environment shapes the cutaneous microbial assemblages of salamanders, likely through species sorting mechanisms [23]. Indicator species analyses offered some additional support for this conclusion, in that many Plethodon indicators have also been isolated from soil habitats (See Appendix D).
Bacterial symbionts that are commonly found in skin assemblages of a particular host taxon, known as indicator species, may play important ecological roles, such as pathogen defense or structural stabilization of microbial communities (e.g., as keystone species). Some similarities were found between the indicators of Desmognathus and Eurycea salamander skin assemblages, especially with respect to antifungal strains. Notably, strains from the family Burkholderiaceae were identified as indicators for salamanders of the genus Desmognathus and for the species Eurycea longicauda. The Burkholderiaceae are core members of skin assemblages on distantly related amphibians [8,9], and they may contribute to Bsal resistance. Several members of this family are potent chytrid antagonists, including strains of Janthinobacterium, Ralstonia, and Burkholderia [49]. In addition, OTUs from the family Chitinophagaceae were identified as indicators for D. monticola and E. wilderae. We presume these chitin-consuming bacteria may protect hosts from fungal pathogens, as seen in the rhizosphere microbiome [50]. Putatively antifungal Flavobacterium strains [see 50] were also found to be significant indicators of the Bsalresistant genus Desmognathus and the species D. fuscus, but not for other host groups. The Verrucomicrobia are underrepresented in culture [51], and, consequently, antifungal properties are unknown; however, this phylum appears to be ecologically important. Strains of Verrucomicrobia have been documented as indicator species and core components of skin assemblages across multiple salamander host groups, including Desmognathus and Eurycea species (Fig. 1, Appendix D), as well as newts (Lissotriton vulgaris and Triturus cristatus) [52]. Symbiotic associations between other indicators (i.e., Xanthomonadales, Cytophagales, Polyangiaceae, and Micromonosporaceae spp.) and their hosts were difficult to assess, as these bacteria have been isolated from a wide range of habitats [53,54].
As expected for terrestrial hosts, many soil-associated bacteria were pinpointed as important members of Plethodon skin assemblages. Two strains from the phylum WPS-2 (Candidatus Eremiobacterota, previously found in a variety of terrestrial environments) [55] were found to be indicator species of Plethodon jordani, as were a number of Acidobacteria for both P. jordani and P. glutinosus. Thus, it appears that soil serves as an important contributor of symbiotic bacteria to Plethodon assemblages. Comparisons among these assemblages and soil samples collected from capture locations could yield evidence of bacterial exchanges between soil and skin [12,13].
Site-specific climate analyses provided a more fine-scale view of the role of the environment in shaping skin assemblages and further reinforced the conclusion that host-related factors and the environment are interacting to shape the cutaneous microbial assemblage. We observed a gradient of climate-driven variation within sites that decayed with increasing water dependence of the host, and this gradient may explain the genus-ecoregion interaction. For the terrestrial genus Plethodon and, to a smaller extent, the semiaquatic genus Eurycea, skin-bacterial assemblages varied with landscape-level climatic factors. Desmognathus assemblages, on the other hand, were resistant against changing climate conditions, perhaps due to buffering by the aquatic habitat. Seasonal temperature fluctuations, for example, are milder in aquatic environments and may have a lesser effect on bacterial assemblages of stream-dwelling salamanders. Conversely, terrestrial salamanders of the genus Plethodon are less capable of seeking out cooler microclimates, so their skin assemblages likely covary with temperature.
Host-specific assemblages have been linked to variation in susceptibility to Bsal infection observed across plethodontid species [18,56] and are therefore important to consider during a response to a Bsal outbreak. Chytridiomycosis is known to be lethal in Eurycea wilderae, although susceptibility varies by population [57]. Plethodon glutinosus are tolerant to infection [56], but the sympatric species Desmognathus fuscus and D. monticola are resistant [56,58]. We found that Desmognathus skin-bacterial assemblages were relatively resistant to landscape-level climatic variation and rich in putatively antifungal bacteria from family Burkholderiaceae and genus Flavobacterium. Future work should investigate whether either of these features of Desmognathus assemblages (resistance against changes in climate and/or richness of antifungal strains) facilitates pathogen defense. Although a high prevalence of antifungal skin bacteria is undoubtedly beneficial for the metazoan host [18], we surmise that it is the structural stability of the cutaneous microbiome, rather than species composition, which likely governs its immune function, and consequently, manifestation of disease in the host [3,10].