Skin microbiota differs drastically between co-occurring frogs and newts

Diverse microbial assemblages inhabit amphibian skin and are known to differ among species; however, few studies have analysed these differences in systems that minimize confounding factors, such as season, location or host ecology. We used high-throughput amplicon sequencing to compare cutaneous microbiotas among two ranid frogs (Rana dalmatina, R. temporaria) and four salamandrid newts (Ichthyosaura alpestris, Lissotriton helveticus, L. vulgaris, Triturus cristatus) breeding simultaneously in two ponds near Braunschweig, Germany. We found that bacterial communities differed strongly and consistently between these two distinct amphibian clades. While frogs and newts had similar cutaneous bacterial richness, their bacterial composition strongly differed. Average Jaccard distances between frogs and newts were over 0.5, while between species within these groups distances were only 0.387 and 0.407 for frogs and newts, respectively. At the operational taxonomic unit (OTU) level, 31 taxa exhibited significantly different relative abundances between frogs and newts. This finding suggests that chemical or physical characteristics of these amphibians' mucosal environments provide highly selective conditions for bacterial colonizers. Multi-omics analyses of hosts and their microbiota as well as directed efforts to understand chemical differences in the mucosal environments (e.g. pH), and the specificities of host-produced compounds against potential colonizers will help to better understand this intriguing pattern.


Introduction
The mucosal environment of amphibian skin provides a habitat in which microbes can live. These microbial assemblages form communities that interact and can provide functional benefits for the host, including protection against skin pathogens, such as Batrachochytrium dendrobatidis (Bd) [1]. Numerous studies have 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

DNA extraction and sequencing
Whole-community DNA was extracted from swabs with the MoBio PowerSoil-htp 96-well DNA isolation kit (MoBio, Carlsbad, CA, USA) following the manufacturer's protocol with minor adjustments, i.e. a 10-minute incubation at 65°C after C1 addition, increased centrifuge time to account for slower rotor speed available and 10-minute incubation at room temperature after C6 addition to the spin column. DNA extracts were stored at −20°C until further processing. The dual-index approach was used to PCRamplify the V4 region of the bacterial 16S rRNA gene using the 515F and 806R primers [22]. Both the forward and reverse primers included the Illumina adapter and a unique barcode tag. Polymerase chain reaction (PCR) methods followed those used in a previous study [23]. Samples were pooled together in approximately equal concentration (as determined by gel band strength), and then cleaned using a Qiagen MiniElute Gel Extraction Kit. A Qubit fluorometer with a broad-range dsDNA kit was used to quantify the DNA after clean-up. The pooled PCR amplicons of all samples were sequenced using paired-end 2 × 250 v2 chemistry on an Illumina MiSeq sequencer at the Helmholtz Center for Infection Biology in Braunschweig, Germany.

Sequence processing
Quantitative Insights Into Microbial Ecology (MacQIIME v. 1.9.1) was used to process all sequence data, unless otherwise stated [24]. Forward and reverse reads from each sample were joined using Fastq-join with default settings [25]. Joined reads were quality-filtered to remove low-quality sequences with the QIIME defaults. Sequences were also filtered by read length to include only the reads with a length between 250 and 253 bp (usegalaxy.org). Chimeras were identified on a per-sample basis using usearch61 denovo-based detection within QIIME (http://drive5.com/usearch/usearch_ docs.html) [26], and subsequently removed. After filtering, 644 334 sequences remained for analysis. Sequences were clustered into operational taxonomic units (OTUs) at 97% similarity using an open reference OTU-picking strategy ( [27], http://qiime.org/tutorials/open_reference_illumina_processing. html). Sequences were first matched to the SILVA 119 (24 July 2014) release (https://www.arb-silva.de) and non-matching reads were subsequently clustered de novo using UCLUST [28]. The most abundant sequence from each OTU was selected as a representative sequence. These representative sequences were aligned using PyNAST [29] and taxonomy was assigned using the RDP classifier [30] with the SILVA

Results
Species richness and phylogenetic diversity differed across amphibian species and not between locations (figure 1). This species effect was driven by R. dalmatina exhibiting lower richness and diversity than the other species (table 1). There was no main effect of amphibian order for Chao1, OTU richness or phylogenetic diversity; however, there was an effect for Shannon diversity (ANOVA, F 1,138 = 9.479, p = 0.002). Newts harboured a greater diversity with respect to the Shannon Index.
The cutaneous bacterial community structure differed between frog and newt hosts at both Elm and Kleiwiesen (figure 2, PERMANOVA Pseudo-F = 31.072 (Elm), 8.4498 (Kleiwiesen), p (MC) = 0.001), while no significant differences were observed between species within their respective amphibian orders (figure 2, PERMANOVA, Pseudo-F = 1.327 (Elm), 0.8123 (Kleiwiesen), p (MC) > 0.05). The average Jaccard distance between samples of amphibian groups was 0.546 ± 0.067 (s.d.) and 0.528 ± 0.069 (s.d.) for Elm and Kleiwiesen, respectively. Within each amphibian group, the average Jaccard distance across both localities was 0.407 ± 0.095 (s.d.) and 0.387 ± 0.07 (s.d.) among newt samples and frog samples, respectively. More closely related newt species (i.e. those from the same genus) were not more similar to each other than to other species from different genera. Pair-wise distances between species are provided in table 2.
Frog and newt bacterial communities were composed of many of the same families, but the relative abundances of these taxa differed (figure 2b). Flavobacteriaceae (Frog/Newt: 33.4/15.8%), Enterobacteriaceae (12.3/1.7%) and Moraxellaceae (12.7/1.2%) exhibited greater relative abundance on frogs, while Pseudomonadaceae (8.4/20.7%), Comamonadaceae (4.2/10.8%), Xanthomonadaceae (2.0/7.9%) and Oxalobacteraceae (1.2/3.6%) were more abundant on newts. The pond water was composed of the following dominant groups: Flavobacteriaceae (Elm/Kleiwiesen: 17.1/20.2%), Figure 1. Richness and diversity of skin microbial communities across species at the two sampling sites, Elm and Kleiwiesen (both close to Braunschweig, Germany), using four metrics: OTU richness, Chao1, Shannon diversity and Faith's phylogenetic diversity. Amphibian species are presented with the species sampled at both sites first, followed by those only sampled at one location within each amphibian order and site. Main effect ANOVA results for species and location are provided below each plot. Of the 46 bacterial OTUs found on 100% of individuals of at least one amphibian species, LEfSe analysis revealed that 31 OTUs were differentially abundant between frog and newt skin communities across both sites. Nineteen OTUs were more abundant on newts and 12 OTUs were more abundant on frogs, and these differences were consistent across the two sampled locations (figure 3). In multiple cases, such as an Empedobacter sp. (frogs), Acinetobacter sp. (frogs), Elizabethkingia sp. (frogs) and a Phyllobacteriaceae sp. (newts), the LEfSe-detected OTUs were nearly absent from the communities of the other amphibian group; that is, for example, the Empedobacter OTU was significantly more abundant on frogs and absent from the communities of newts on all but four individuals that showed very low relative abundances (figure 3; electronic supplementary material, figures S1 and S2).  Relative abundances of bacterial taxa at the family level for hosts at the two sampling locations, Elm and Kleiwiesen (close to Braunschweig, Germany). Taxonomic information for dominant groups is listed at the side. For the bar plots, amphibian species are presented with the species sampled at both sites first, followed by those only sampled at one location within each amphibian order and site.
We defined core OTUs as those present on 90% of individuals of the respective species at a given location. The 90%-core microbiomes were as follows, for Elm and Kleiwiesen respectively: R. dalmatina-  (table 3). Many of these core members were also found to be differentially abundant between frogs and newts with the LEfSe analysis (table 3). However, three OTUs were members of the core for all amphibian species across both sites, including a Pseudomonas fluorescens, an Albidiferax sp. and a Chryseobacterium sp. (Table 3).

Discussion
Cohabiting amphibian species, especially during their aquatic phase, offer the opportunity to test the effects of host factors in structuring skin microbial communities by removing extrinsic environmental factors. We hypothesized that amphibian species would harbour distinct bacterial communities, which was, in part, upheld by our results. When sampled at the same time-point, distinct communities were found to inhabit frog and newt hosts, but within these amphibian orders, the species did not have significantly different community structures. This pattern was consistent across both Elm and Kleiwiesen. We did not observe any phylogenetic patterns with respect to the microbiota of newt species (i.e. more closely related newt species did not exhibit greater similarity in their microbial communities).  Relative abundances within the full community of four selected differential OTUs across samples from both locations. Two OTUs exhibit greater relative abundance in newts and two OTUs exhibit greater relative abundance in frogs. Each bar corresponds to an individual amphibian and the colours denote the sampling locations, Elm (black) and Kleiwiesen (grey) (both near Braunschweig, Germany). Electronic supplementary material, figures S1 and S2 provide relative abundance plots for the remaining 29 LEfSe-detected differentially abundant OTUs.