Dramatic declines in amphibian populations around the world have been a subject of research in recent years. Since the 1980s, amphibian populations have declined to the point where 42% of the 7854 species described worldwide (AmphibiaWeb, 2018) have experienced some form of population declines (Stuart et al., 2004; Whittaker et al., 2013). Population declines in amphibians have been attributed to several factors, mainly related to human activities such as habitat destruction, pesticide abuse and climate change, among others (Whitfield et al., 2007; Kerby et al., 2011). However, population declines in protected areas, so called “enigmatic declines,” have drawn the attention to the emerging disease chytridiomycosis, caused by the fungus Batrachochytrium dendrobatidis (Bd) (Lips et al., 2006; Whitfield et al., 2007, 2012; Lips, 2016).

Costa Rica has become a focus for amphibian research. In the late 1990s, a study revealed declines of up to 40% of amphibian populations in the Monteverde Cloud Forest Reserve (Pounds et al., 1997). Also, Lips (1998) documented population declines in Costa Rican amphibians at Las Tablas. Another hotspot for amphibian studies in Costa Rica is La Selva Biological Station in Heredia Province. This is a 16 km² private protected area, managed by the Organization of Tropical Studies, where amphibians have been studied for almost 40 years. Whitfield and collaborators analyzed data collected for over 35 years and they found that all species of terrestrial amphibians declined approximately 75% in this site. Interestingly, that study showed a correlation between population declines and a climate driven reduction of the quantity of leaf litter, but no correlation to Bd emergence was detected (Whitfield et al., 2007).

Over the last years, next generation sequencing technologies allowed scientists to explore oceans, soils, humans and a variety of host related microbiota (Jani and Briggs, 2014). These technologies can be implemented to study amphibians, specifically their skin microbiota. Skin is especially important for amphibian physiology; it is responsible for the exchange of respiratory gasses and osmoregulation (Voyles et al., 2009). Furthermore, the skin is in continuous exposure to microorganisms in the environment, hence these potentially symbiotic microbes may be involved in different physiological processes and interacting with pathogens such as Bd (Madison et al., 2017). Knowledge of the ecology of bacterial communities associated with amphibians has only a few years of development (Bletz et al., 2017a). Previous research showed that bacterial communities of amphibian skin were significantly different to their environments and one of the most important factors driving the structure of the microbiota was host species (McKenzie et al., 2012; Kueneman et al., 2013; Sabino-Pinto et al., 2016). More recently, it has been suggested that environmental factors may also play an important role in the variation of the amphibian bacterial microbiota (Muletz et al., 2017; Bird et al., 2018; Varela et al., 2018). It has also been reported that bacterial communities can vary according to geography (Abarca et al., 2018; Rebollar et al., 2018), ecomorphs (Bletz et al., 2017a) and ontogeny (Longo et al., 2015).

The discovery of bacteria with antagonistic activity to Bd (Harris et al., 2006, 2009) opens the possibility of employing new strategies to control the pathogen; however, due to the rapid loss of species and habitat, it is necessary to speed up the search for beneficial bacteria to be employed in captivity or in wildlife (Flechas et al., 2017). To accelerate this search, the use of culture independent techniques can be a very efficient tool to describe local bacterial communities (Rebollar et al., 2016a). Due to the variation in the inhibitory response of bacteria to Bd (Antwis and Harrison, 2018), and because skin bacterial communities differed strongly from those of the surrounding environment (Bletz et al., 2017a), it has been recommended to implement mitigation strategies at the local level, using bacteria to counteract the circulating pathogen in the same region (Flechas et al., 2017). The aim of this study is to explore and compare the diversity of the skin bacterial communities of twelve tropical amphibian species at La Selva Biological Station in Costa Rica. Further, these communities are compared to previous studies with other amphibian species and were examined for potential anti-Bd function. Our goal is to describe this host-associated microbiota and set the foundation to understand microbiota composition and potential interaction with Bd, as well as to describe host specific relationships between amphibians and microbes in Costa Rica.

After rarefaction, we obtained a total of 2,621,143 good quality sequences and 6330 OTUs distributed in 38 bacterial phyla (Figure 1A). Eight out of the twelve species under study were dominated by the phylum Proteobacteria (43.95% of the total OTUs). Other abundant phyla were Bacteroidetes (35.59%) and Actinobacteria (8.79%), while Acidobacteria (3.91%), Firmicutes (2.61%), Verrucomicrobia (1.79%) and Planctomycetes (1.04%) were other phyla represented in at least 1% of the OTUs. Among the four species that were not dominated by Proteobacteria, Hypopachus pictiventris, Leptodactylus savagei, and Rhaebo haematiticus presented a majority of Bacteroidetes OTUs, with 89.02, 69.47, and 39.40%, respectively. The bacterial community associated with Rhinella horribilis showed a dominance of Bacteroidetes (43.23%) and Actinobacteria (42.08%), besides this species, only Smilisca manisorum presented a relative abundance of actinomycete OTUs above 10% (23.72%).

We identified 463 bacterial families in our dataset. In terms of relative abundance of OTUs, the Sphingobacteriaceae family presented the highest percentage (13.82%), followed by Comamonadaceae (9.08%), Weeksellaceae (8.14%), Chitinophagaceae (7.95%), and Alcaligenaceae (7.46%) (Supplementary Table 2). The phylogenetic distribution at family and genus level was noteworthy (Figure 1B), for example, showing the dominance of Pedobacter in the skin microbiota of H. pictiventris (84.34%). In C. bransfordii the bacterial community was weakly dominated by OTUs of the genus Chryseobacterium (7.00%) and Agrobacterium (4.51%). The bacterial community of C. fitzingeri was dominated by OTUs of the genus Sphingobacterium (11.34%) and one OTUs in Chitinophagaceae family (6.28%); C. megacephalus had a bacterial community dominated weakly by Agrobacterium (5.10%). The bacterial community of O. pumilio had a dominance of OTUs related to Enterobacteriaceae (16.54%) and the bacterial community of R. haematiticus was dominated by the genus Sphingobacterium (11.98%). The bacterial community of A. callidryas was dominated by Alcaligenaceae OTUs (28.04%). The Cane toad, R. horribilis, showed an abundance of OTUs in the genus Niabella, family Chitinophagaceae (36.75%). P. ridens was dominated by OTUs in the family Alcaligenaceae (55.7%), while the L. savagei community was dominated by Chryseobacterium (37.21%) followed by another OTU belonging to the family Flavobacteriaceae (22.82%).

Figure 2 shows a heatmap of the core microbiota (OTUs that are present in 90% of the samples), as well as the OTUs that have relative abundances significantly different between the bacterial communities. All amphibian species shared OTUs belonging to the genera Chryseobacterium and Pedobacter, as well as several OTUs assigned to the families Alcaligenaceae, Comamonadaceae and Rhizobiaceae. The three most abundant OTUs detected in this study were assigned to the family Enterobacteriacea and to the genera Methylibium and Pedobacter.

The relative abundance of OTUs varied among individuals within each species, as follows: In O. pumilio (n = 25), 12 individuals were dominated by Enterobacteriaceae OTUs, 7 individuals were dominated by Methylibium OTUs, while the other 6 had higher abundance of Chitinophagaceae and Stenotrophomonas. In C. fitzingeri (n = 10), 3 individuals were dominated by Sphingobacterium OTUs, 2 by Chryseobacterium OTUs and the 5 remaining individuals each had a different most abundant OTU. Regarding R. haematiticus (n = 10), 4 individuals were dominated by Sphingobacterium OTUs, 2 by Chryseobacterium OTUs and the 4 remaining individuals each had a different most abundant OTU. In two species, C. bransfordii (n = 5) and A. callidryas (n = 4), all individuals were dominated by a different OTU. Finally, all species with 2 individuals sampled had the same dominant OTU (G. pictiventris, P. ridens and D. diastema), except for C. megacephalus, in which one individual was dominated by Chryseobacterium and the other by Agrobacterium OTUs.

The beta diversity analysis did not show clear clustering among the different bacterial communities; however, some trends are worth noticing. Figure 3 shows H. pictiventris communities separately from the rest of the frogs. Other species did cluster together such as D. diastema, while some samples of O. pumilio and one of A. callidryas grouped close to them. Most O. pumilio samples are widely distributed around the plot, indicating great variation within this species. Beta-diversity analysis showed that host species (ANOSIM R = 0.416, p = 0.01; PERMANOVA PseudoF = 3.388, p = 0.001) and host family (ANOSIM R = 0.392, p = 0.01; PERMANOVA PseudoF = 3.751, p = 0.001) are significant predictors of the variation in bacterial community composition. Contrary, host habitat did not show significant results (ANOSIM R = 0.171, p = 0.11; PERMANOVA PseudoF = 1.866, p = 0.0599). The combined effect of host species and host habitat also influenced significantly the variation in bacterial community composition (ADONIS F = 3.3884, p = 0.001). Regarding alpha-diversity, both Shannon and phylogenetic diversity indexes presented similar results (Figure 4): they did not vary significantly between amphibian species (P > 0.005 in both indexes). Frogs of the genus Craugastor, O. pumilio, R. haematiticus and S. baudinii have the most diverse skin bacterial communities. On the other hand, D. diastema, P. ridens and H. pictiventris showed lower diversity (Supplementary Table 3). Regarding habitat, there was a significant variation between indexes according to habitat (p < 0.005 in both), showing a higher diversity associated to terrestrial species when compared to arboreal species. On the other hand, there were no differences between diversity indexes when we compared terrestrial species versus aquatic species or arboreal species versus aquatic ones (p > 0.005 in all indexes).

Finally, we compared tropical and temperate species, but no significant latitudinal differences were found in alpha diversity; furthermore, two overlapping families, Bufonidae and Hylidae, also presented similar alpha diversity values (Figure 5). The proportion of putatively anti-Bd bacteria showed no significant differences among species (Kruskal–Wallis test, p = 0.149) or between arboreal and terrestrial habitat types (Mann–Whitney U-test, 0.714). However, all 12 species and all samples contained putatively anti-Bd bacterial taxa with an average of 13.3% of the 16S rRNA sequencing reads (Figure 5).

Amphibian skin is an important organ for electrolyte and gas exchange with the environment, and many diseases and contaminants can affect amphibians through it (Campbell et al., 2012). Unfortunately, very little is known about amphibian-associated microbiota in a biodiversity hotspot that has suffered from severe amphibian decline such as Costa Rica. To start bridging this knowledge gap, we studied the skin microbiota of several frog species at La Selva Biological Station. We found that skin bacterial communities differ among amphibian species. Species-level differences in the amphibian skin microbiota have been previously documented in temperate (McKenzie et al., 2012; Kueneman et al., 2013; Walke et al., 2014) and tropical systems (Flechas et al., 2012; Belden et al., 2015; Rebollar et al., 2016b). Furthermore, we found that host family is also a good predictor of bacterial community variation. In tropical regions there are more frog families occupying different microhabitats, this is not the case for temperate zones, where host family was not a strong driver of microbiota when compared to host habitat (Kueneman et al., 2013).

It has been reported that host ecology and habitat are more important factors of variation than host taxonomy (Bletz et al., 2017a; Muletz et al., 2017; Bird et al., 2018), and that most amphibian microbiota comes from the environment (Rebollar et al., 2016b). Despite this, we did not find that habitat was a predictor of variation in the bacterial community of the population under study. This result may be caused by the limited sampling range, or it could reflect amphibian ecology at La Selva, where species transition between different habitats. For example, O. pumilio is principally terrestrial but they climb trees to bring tadpoles to bromeliads, R. haematiticus can be seen in the forest floor as well as in rivers or even in low vegetation and C. fitzingeri can be observed in forest leaf litter, vegetation, near rivers and open areas with pastures (Savage, 2002). However, further research is needed to confirm these results, including analyzing more species from aquatic habitats.

Proteobacteria was the most commonly detected phylum and was more abundant in leaf litter species. This is the main phylum recovered from soils (Kim et al., 2014), tree leaves (Kembel et al., 2014) and amphibians (Belden et al., 2015; Bletz et al., 2017a; Varela et al., 2018) in the tropics. The other three dominant phyla agree with the main four phyla of bacteria associated with amphibians in temperate or tropical areas (Lauer et al., 2008; McKenzie et al., 2012; Walke et al., 2014; Becker et al., 2014; Belden et al., 2015). However, we found that Comamonadaceae and Alcaligenaceae were the most common Proteobacteria families in the bacterial communities of the amphibians sampled, in contrast with other studies in both tropical and temperate zones that identified Pseudomonadaceae as the most abundant family (Lam et al., 2010; Rebollar et al., 2016b, 2018).

Bacteroidetes was the second most abundant phylum in skin of the frogs at La Selva. This phylum is also widely distributed in the environment, including soils, sediments and sea water (Thomas et al., 2011). Bacteroidetes is common in four aquatic and terrestrial frog species. Important genera within this phylum present in the microbiota of amphibians at La Selva were Pedobacter, Chryseobacterium and Niabella. High dominance of Pedobacter OTUs was found in H. pictiventris. The abundance of this genus is important because it has been shown to have inhibitory activity against Bd (Harris et al., 2006; Lauer et al., 2008). Pedobacter has been isolated from multiple habitats and it can dominate aerobic communities in different aquatic and terrestrial environments, due to its high metabolic diversity (Gordon et al., 2009). H. pictiventris spends most of its life underground but in the breeding season they congregate in puddles in the forest (Savage, 2002). This sudden change of microclimate could affect the skin bacterial community allowing Pedobacter to thrive over other genera. However, further research is needed to establish whether these two organisms maintain a symbiotic relationship.

Craugastor spp., R. haematiticus, R. horribilis, and O. pumilio had higher overall bacterial diversity compared with the rest of species; P. ridens and H. pictiventris have the least bacterial diversity. At La Selva, alpha diversity metrics as Shannon and Phylogenetic diversity indexes show that typical leaf-litter dwelling families (Craugastoridae, Bufonidae and Dendrobatidae) have greater bacterial diversity than arboreal or aquatic families such as Hylidae, Phyllomedusidae, Leptodactylidae, Microhylidae or Eleutherodactylidae. It is possible that more soil bacteria can colonize amphibians in the leaf litter than in arboreal or aquatic habitats, this pattern agrees with previous studies (Belden et al., 2015; Rebollar et al., 2016b; Bletz et al., 2017a).

Our results show differences with similar studies performed in nearby Panama. In those studies, different frog species were analyzed (except for C. fitzingeri and A. callidryas), suggesting that differences in bacterial composition between countries may be due to different host species tested. One important difference is that Actinobacteria was the second most abundant phylum in frogs of Panama (Belden et al., 2015; Rebollar et al., 2016b), while in Costa Rica this group was in third place. For shared species, however, the alpha diversity indexes appear to be similar, with higher values for C. fitzingeri. The microbiota of C. fitzingeri from La Selva is enriched with OTUs of the Sphingobacteriaceae family, while those of Panama are enriched with other families such as Comamonadaceae and Moraxellaceae (Rebollar et al., 2018). The Shannon index of C. fitzingeri in Costa Rica is similar to that of C. fitzingeri in “Bd-naive” sites in Panamá (Rebollar et al., 2016b), although La Selva is a “Bd-endemic” site where it has been reported that C. fitzingeri along with other species of amphibians carry Bd (Whitfield et al., 2012). Additionally, in our study, in every species with more than 4 amphibians sampled, the most abundant OTU was not the same in all of them. Both observations suggest that there are intraspecific variations in amphibian bacterial communities, these variations are possibly associated with environmental factors, presence of Bd or internal population variations of amphibian species, as has been reported by other studies (Muletz et al., 2017; Abarca et al., 2018; Rebollar et al., 2018).

We found similarities in Shannon diversity index and phylogenetic diversity between temperate zones and tropical zones. At bacterial composition level, temperate zones frogs appear to have more Firmicutes than tropical ones (Kueneman et al., 2013; Loudon et al., 2014; Belden et al., 2015). It is important to point out that the temperate studies included mostly aquatic species, and only one terrestrial family: Bufonidae, while at La Selva we studied mostly terrestrial and arboreal species; therefore, we can only compare two families between Costa Rica and temperate zones: Bufonidae and Hylidae. Those families showed similar richness and phylogenetic diversity. In Bufonidae, Kueneman et al. (2013) found similar abundance of three phyla, the same found in R. horribilis at La Selva. Both in Kueneman et al. (2013) and in Belden et al. (2015) the aquatic temperate species was Lithobates catesbeianus, and they found different microbiota composition between tropical and temperate populations. To establish a better comparison, it is necessary to characterize the bacterial community of a fully aquatic tropical frog, ideally in the Family Ranidae. Even though the tropics present very different conditions than temperate zones, it is possible that the skin texture and habits of certain species may represent a more homogeneous microhabitat for some bacterial genera regardless of geographical location.

We recovered several sequences related to putatively anti-Bd bacterial strains in the Antifungal Isolates Database (Woodhams et al., 2015). This database has been used to determine proportions of OTUs of anti-Bd bacteria among other species of amphibians (Kueneman et al., 2016), as well as to study associations between the degrees of infection and the presence of OTUs related to inhibitory bacteria (Jani et al., 2017). Chryseobacterium and Pseudomonas comprised the majority of putatively anti-Bd bacteria from all species examined. These two genera have been reported as the most abundant among inhibitory Bd bacteria within their phylum (Park et al., 2014; Becker et al., 2015). In this study over 13% of sequences were related to putatively anti-Bd strains. This percentage of putatively anti-Bd bacteria is similar to that found using culture-based techniques in Costa Rican amphibians (11.1% in cell-free supernatant assays) (Madison et al., 2017). The presence of OTUs associated with putatively anti-Bd bacteria suggests the existence of bacteria that can help fight the disease, but this is still a hypothetical value. For this reason, inhibition tests with cultivable bacteria from each amphibian species (Antwis et al., 2015; Bletz et al., 2017b; Burkart et al., 2017) are necessary to confirm if some amphibians indeed carry higher proportions of anti Bd bacteria.

Understanding the microbial community in each region is one of the first steps to implement mechanisms to mitigate chytridiomycosis, including employing probiotic bacteria (Woodhams et al., 2016). Here we find that skin bacterial community of amphibians at La Selva varies according to species and family. Within the bacterial community there are potential symbiotic bacteria that should be the focus of future studies. It is important to know and describe local bacterial communities associated with amphibian skin to understand the mechanisms of interaction between these microorganisms and thus implement efficient bioaugmentation strategies (Flechas et al., 2012). This work represents an effort to characterize the bacterial communities of amphibians in Costa Rica. La Selva Biological Station is a tropical forest research center of worldwide recognition, and this research can serve as a foundation for further studies aimed at designing mitigation strategies to control the threat of Bd and human-associated factors as drivers of amphibian extinction.

This study was carried out in accordance with the recommendations of the Institutional Commission on Biodiversity at the University of Costa Rica. Our work with amphibians has been reviewed and approved by the Animal Care and Ethics Committee from the University of Costa Rica (Resolution CICUA-029-15).

SW, JK, VM and AP-T conceived and designed the experiments. JA, GV, IZ, and CM-C performed the experiments. JA, GV, IZ, SW, DW, JK, and VM analyzed the data. JK, VM, and CM-C contributed reagents, materials, and analysis tools. JA, GV, IZ, DW, VM, and AP-T wrote the paper. SW, DW, JK, and CM-C reviewed draft on the paper.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer MV declared a past co-authorship with several of the authors VM and DW to the handling Editor.