Field sampling occurred during three field visits: 14 August – 12 September 2013, 4 January – 9 February 2014, and 5 November – 15 December 2014. In total, 1021 microbial samples (989 frog skin swabs, 32 environmental samples) were collected from 10 locations (30 sites) and 96 host species (Figure 1 and Supplementary Tables 1, 2).

Amphibians were captured during day and night surveys with clean nitrile gloves and placed in sterile Whirl-Pak^® bags (Nasco, Fort Atkinson, WI, United States). For skin microbe sampling, individuals were removed from the bag with a clean pair of nitrile gloves and rinsed with 50 ml of sterilized water. After rinsing, each individual was swabbed with a single sterile rayon swab (MW113, Medical Wire Equipment & Co. Ltd., Corsham, United Kingdom), applying 10 strokes on the ventral abdomen, 5 strokes on each ventral thigh, and 5 strokes on each foot. Swabs were stored in microcentrifuge tubes and transported on ice until transfer to a -20°C freezer. Environmental samples were collected from amphibian-associated habitat, including soil, water, and leaves. For soil samples, 1–2 g of soil were collected; for water samples, 60–120 ml of water were hand-pumped through a 0.22 μm filter; for leaf samples, the surface was swabbed with 30 strokes. Environmental samples were stored in 2 ml tubes and transported on ice until transfer to a -20°C freezer. This study was approved by the Institutional Animal Care and Use Committee of James Madison University (protocol #A01-15), and necessary research and access permits were obtained from the Malagasy Direction Générale des Forêts (DGF) and Madagascar National Parks for all sampling.

Bacterial DNA was extracted using the MoBio PowerSoil DNA isolation kit (MoBio Laboratories, Carlsbad, CA, United States) following the manufacturer’s protocol with minor modifications to increase DNA yield. The V4 region of the 16S rRNA gene was PCR-amplified in triplicate with barcoded primers (515f/806r) following Kueneman et al. (2014). Amplicon concentration was quantified with Quant-iT PicoGreen dsDNA Assay kit (Thermofisher, Waltham, MA, United States). Equal concentrations of each sample were pooled, and the pooled amplicons were cleaned using MoBio UltraClean PCR Clean-up kit (MoBio Laboratories, Carlsbad, CA, United States). The pooled barcoded amplicons were sequenced using 2 × 150 pair-end technology on an Illumina MiSeq platform at the BioFrontiers Institute at the University of Colorado.

Sequence reads were filtered and pre-processed in Quantitative Insights Into Microbial Ecology (QIIME) (Caporaso et al., 2010). Forward reads were demultiplexed and filtered with the following criteria to retain only high quality reads: no Ns within the sequence, no barcode errors, and a minimum of three consecutive low-quality base pairs (minimum q = 10) before read truncation. Only forward reads were used because reverse reads typically suffer from lower quality (Kwon et al., 2013). Quality filtered sequences were clustered into operational taxonomic units (OTUs) using the deblur workflow¹. Deblur is a new sub-operational taxonomic unit (sOTU) approach for amplicon sequencing that incorporates known Illumina error profiles and uses Hamming distances along with a greedy algorithm (Amir et al., 2017). Within this workflow, sequences were trimmed to 150 bp, and sOTU clusters (hereafter called OTUs) with less than 25 reads were removed. Taxonomy was assigned with Ribosomal Database Project Classifier (Wang et al., 2007), and a phylogenetic tree was built in QIIME using the fasttree algorithm (Price et al., 2010). Samples were subsequently rarefied at 4,000 reads per sample to normalize read counts across samples. Sequences have been archived in the SRA database (Bioproject accession number: PRJNA394790).

Number of OTUs, Chao1, effective number of species [exp(Shannon index), Jost, 2006], and Faith’s phylogenetic diversity were calculated for all samples as measures of alpha diversity, i.e., species richness and species diversity. General linear models (GLM) were used to test which factors were significant predictors of amphibian skin bacterial richness and diversity in SPSS v24 (IBM Corp, Armonk, NY, United States). The same factors described below for PERMANOVAs were included in the models to represent site parameters, host ecomorphology and host phylogeny.

Beta diversity was calculated as weighted and unweighted Unifrac distances in QIIME. The resulting matrices were used to explore patterns in beta diversity in two main ways: (1) Multiple Regression on distance Matrices (MRM), and (2) Permutational multivariate analysis of variance (PERMANOVA).

MRM was used to simultaneously explore the role of host phylogeny and host ecology on microbial community structure. We made distance matrices to represent (a) evolutionary divergences among host species (i.e., host phylogeny) and (b) ecomorphological differences (i.e., host ecology). Evolutionary divergences among host species (with 2 or more sampled individuals) were calculated as patristic phylogenetic distances from an ultrametric timetree. We first reconstructed a phylogenetic tree of all host species included in our study based on partial sequences of the mitochondrial 16S rRNA gene (Vieites et al., 2009), in MEGA7 (Kumar et al., 2016) under the Maximum Likelihood optimality criterion, with a general time-reversible (GTR + I) substitution model and using SPR branch swapping. We then manually corrected the tree topology for some wrongly reconstructed deep relationships, based on published multigene phylogenies of mantellids (Wollenberg et al., 2011), microhylids (Scherz et al., 2016), and intrafamilial relationships (Roelants et al., 2007). We then entered this topology as the usertree in MEGA7 along with the 16S alignment, and reconstructed an ultrametric tree using the Real-time method, without absolute calibration. Patristic pairwise distances among included amphibian species were calculated from this ultrametric tree in R using the ape and adephylo packages (Paradis et al., 2004; Jombart and Dray, 2008). Ecomorphological distances were calculated on the basis of major species traits that can be hypothesized to influence the cutaneous microbiota: (1) degree of arboreality, (2) degree of water dependence, (3) kind of breeding water body, (4) primary forest dependence, (5) reproductive mode, (6) aquatic or terrestrial egg deposition, and (7) body size (Glaw and Vences, 2007). The traits were coded as either ordered/ordinal (1,2,4,7) or unordered/categorical (3,5,6) characters. For a complete list of all character states see Supplementary Tables 3, 4. Distance matrices among species were calculated using PAUP v. 4b10 (Swofford, 2002). We considered a priori the first three traits as most likely to be important as they directly indicate distinct microhabitats that the frogs are exposed to, while we considered the remaining traits as possibly important, but probably less influential. MRM analyses were completed with the ‘ecodist’ package in R (Goslee and Urban, 2007; R Core Team, 2014), testing the correlation of the phylogenetic and ecomorphological matrices with microbial community structure for the full dataset (89 species) as well as for amphibians occurring at a single hyperdiverse site (38 species), Ranomafana. Microbial communities were represented by weighted Unifrac distance matrices derived from OTU tables averaged by frog species; that is, the rarified OTU table was first averaged by host species, and pair-wise weighted unifrac distances were subsequently calculated.

Topological congruency analysis was used to further explore the possible existence of phylosymbiosis. For this, we quantified congruence between the host phylogenetic tree and microbial dendrograms using the TreeCmp program. Microbial dendrograms were created in QIIME using both Weighted Unifrac and Bray-Curtis distances derived from OTU tables averaged by frog species. Using TreeCmp, we calculated the normalized Robinson–Foulds scores, where values of 0 indicates complete congruence and values of 1 indicate lack of congruence.

Because no topological congruence was found and the MRM analyses suggested a stronger correlation with host ecology than with phylogeny (see Results), we performed PERMANOVAs to further understand how frog ecomorphology as well as the other variables influence the skin bacterial communities. PERMANOVAs were completed in R (R Core Team, 2014) with the ADONIS2 function in the ‘vegan’ package (Oksanen et al., 2017) to test which factors significantly explained the observed variation in microbial community composition and structure. The “margin” option in ADONIS2 was used to assess marginal effects of each term in a model including all other variables. The following factors were included in the models: elevation, latitude, ecomorphological category [hereafter called host ecomorph, with three categories: (a) arboreal, (b) aquatic and semi-aquatic, or (c) terrestrial], as well as an approximate representation of host phylogeny. The phylogeny variable was chosen because the large number of frog species sampled made it infeasible to include host species as a categorical variable in any model, and because a categorical host species variable would not have captured the phylogenetic relationships among the amphibians. We, therefore, performed non-metric multidimensional scaling (nMDS) constrained to 1 dimension on the patristic distance matrix using SPSS v24 (IBM Corp, Armonk, NY, United States). The coordinates of the nMDS axis were subsequently extracted and used as a proxy variable in PERMANOVA models. PERMANOVA was also used to test whether frog skin microbiota differed from that of the environment in PRIMER7 (Clarke and Gorley, 2015).

Unweighted pair group method with arithmetic mean (UPGMA) was used to evaluate clustering patterns across host ecomorphs and Similarity Profile Analysis (SIMPROF) was used to statistically test for significant structure within the created UPGMA dendrogram. Both analyses were performed in PRIMER7 (Clarke and Gorley, 2015).

The Linear discriminant analysis Effect size (LEfSe) method (Segata et al., 2011) was used to identify which bacterial taxa were most likely explaining the observed differences between categories of interest. LEfSe was used to identify differentially abundant taxa between frogs and the environment, and also to identify differentially abundant taxa among host ecomorphs. Default parameters were used with the exception of increasing the LDA score; taxa with LDA scores greater than 3.0 were considered significant.

Bipartite networks were used to visualize the association of bacterial taxa with a given host ecomorph. These networks were calculated in R (Sedlar et al., 2016), and visualized with Gephi (Bastian et al., 2009).

Frog cutaneous microbial communities were less species rich than those of the environment (# of OTUs – Frog: 219.9 3 ± 6.71(SE)/Env: 948.96 ± 72.98; Chao1 – Frog: 265.37 ± 8.57/Env: 1415.79 ± 130.76), and also less diverse than that of the environment (Effective # of species – Frog: 509.8 ± 40.0/Env: 8644.5 ± 1337.8; Faith’s PD – Frog: 37.28 ± 0.81/Env: 106.86 ± 5.88).

Frog skin microbiotas were dominated by Proteobacteria (Gamma – 46.6%, Beta – 15.4%, Alpha – 9.4%, Delta – 1.6%), Bacteriodetes (8.1%), Actinobacteria (7.9%), Firmicutes (3.6%), and Acidobacteria (2.3%). Microbial community structure on frog skin strongly differed from that of the environment (PERMANOVA: Pseudo-F = 33.34, p = 0.001, Figure 2A). Soil environments were comprised predominantly of Acidobacteria (19.7%), Proteobacteria (Alpha – 17.3%, Gamma – 8.7%, Beta – 8.5%, Delta – 6.4%), Bacteriodetes (7.3%), Verrucomicrobia (5.3%), Chloroflexi (5.1%) Actinobacteria (4.8%), and Planctomycetes (3.3%). Water environments were comprised of Proteobacteria (Beta – 23.7%, Alpha – 21.4%, Gamma – 9.5%, Delta – 3.9%), Actinobacteria (11.8%), Bacteriodetes (7.7%), Planctomycetes (3.9%), Firmicutes (3.0%), Acidobacteria (2.7%), and Verrucomicrobia (2.5%). Leaf surfaces were comprised of Proteobacteria (Alpha – 30.2%, Gamma – 12.3%, Beta – 5.2%, Delta – 2.8 %), Bacteriodetes (24.7%), Actinobacteria (10.1%), Acidobacteria (3.7%), Cyanobacteria (2.8%), and Verrucomicrobia (2.7%).

Fifty-eight bacterial taxa were identified as differentially abundant between frog hosts and the environment using the LEfSe method (LDA > 3). More specifically, 47 taxa exhibited greater relative abundance in the environment and 11 taxa exhibited greater relative abundance on frogs (Supplementary Table 5). For example, Synthrophobacteraceae (Deltaproteobacteria) and Chthoniobacteraceae (Spartobacteria) were enriched in the environment, while Pseudomonadaceae (Gammaproteobacteria) and Enterobacteriaceae (Gammaproteobacteria) were enriched on frog skin (Figure 2B).

To simultaneously explore the role of host phylogeny and host ecology we performed MRMs. When considering the entire dataset, both host phylogeny and host ecology were not significant (MRM -Phylo: p = 0.400, Eco: p = 0.515). However, when considering the single hyperdiverse site Ranomafana, host ecomorphology, derived from the 3-character matrix, was significant while host phylogeny was not (MRM – Eco: p = 0.029, Phylo: p = 0.579). Inclusion of additional ecological traits in the calculation of the ecomorphological distance matrix did not improve the correlation with weighted unifrac distances of the microbial communities (Supplementary Table 6), suggesting that arboreality, water dependence and breeding water body are the most important ecological traits influencing the cutaneous bacterial communities of Malagasy amphibians.

Topological congruency analysis also showed a lack of congruence between host phylogenetic trees and microbial dendrograms, suggesting a limited effect of host phylogeny (Full dataset – Weighted Unifrac: normalized RF score = 1, Bray-Curtis: normalized RF score = 0.98; Ranomafana-Weighted Unifrac: normalized RF score = 0.98, Bray-Curtis: normalized RF score = 1).

Using PERMANOVAs, host ecomorph was the strongest predictor of skin bacterial community structure [PERMANOVA – Weighted: Pseudo-F_(2,980) = 30.308, p = 0.001; Unweighted: Pseudo-F_(2,980) = 7.142, p = 0.001, Table 1]; however, host phylogeny (nMDS 1), latitude, and elevation also explained significant portions of the variation (Table 1). Pairwise comparisons of individuals within ecomorphs were lower than the pairwise distances between ecomorphs (Supplementary Table 7), and pairwise comparisons of host ecomorph classes showed that each ecomorph was significantly different from the others (Arb-Ter: t = 6.699, p = 0.001; Arb-Aqu t = 5.625, p = 0.001; Ter-Aqu t = 2.875, p = 0.001; Figure 3). Multivariate dispersion also did not differ between host ecomorphs (PERMDISP- F = 1.699, p = 0.223). In addition, UPGMA clustering showed that host species within the same ecomorph class typically clustered together (Figure 3), and SIMPROF analysis revealed that there was significant structure within the dendrogram (SIMPROF – p = 0.001).

Fifty-eight bacterial taxa were identified to best explain the observed microbial community differences among host ecomorphs (LEfSe, LDA > 3); seven were differentially more abundant on arboreal frogs, 39 were more abundant on terrestrial frogs, and 12 were more abundant on aquatic frogs (Supplementary Table 8). For example, the genus Pigmentiphaga (Alcaligenaceae) was significantly enriched on arboreal frogs, Agrobacterium (Rhizobiaceae) was significantly enriched on terrestrial frogs, and Methylotenera (Methylophilaceae) was significantly enriched on aquatic frogs (Figure 3). In particular, the relationship between Pigmentiphaga and arboreal frogs was seen across host genera (e.g., Boophis and Heterixalus) and locations (Figure 4).

For selected frog hosts, variation through time was explored. While bacterial composition exhibited differences between sampling time points, key bacterial taxa [e.g., Alcaligenaceae in arboreal frogs (Boophis)] were consistently present through time (Figure 5).

Richness and diversity of frog skin microbial communities were also affected by multiple factors, with host ecomorphology exerting the strongest influence in most cases (Table 2). Host phylogeny and latitude also influenced species richness and diversity (Table 2). Pair-wise comparisons between host ecomorphs revealed that the main effect of host ecomorph was driven by arboreal frogs having significantly less diverse bacterial communities than both aquatic and terrestrial frogs, and aquatic frog having slightly less diverse communities than terrestrial frogs (Tables 3, 4).

From a presence-absence perspective, there were 990 (15%) OTUs, 1592 (24.1%) OTUs, and 1776 (26.9%) OTUs (summed across all individuals of a given ecomorph) that were unique to arboreal, aquatic and terrestrial frogs, respectively. There were also OTUs shared between these groups; 15.3% of OTUs were shared across all ecomorphs, 7.8% of OTUs were shared between arboreal and terrestrial frogs, 7.3% were shared between terrestrial and aquatic frogs, and 3.6% were shared between aquatic and arboreal frogs (Figure 6).