Moving Beyond the Host: Unraveling the Skin Microbiome of Endangered Costa Rican Amphibians

Some neotropical amphibians, including a few species in Costa Rica, were presumed to be “extinct” after dramatic population declines in the late 1980s but have been rediscovered in isolated populations. Such populations seem to have evolved a resistance/tolerance to Batrachochytrium dendrobatidis (Bd), a fungal pathogen that causes a deadly skin disease and is considered one of the main drivers of worldwide amphibian declines. The skin microbiome is an important component of the host’s innate immune system and is associated with Bd-resistance. However, the way that the bacterial diversity of the skin microbiome confers protection against Bd in surviving species remains unclear. We studied variation in the skin microbiome and the prevalence of putatively anti-Bd bacterial taxa in four co-habiting species in the highlands of the Juan Castro Blanco National Park in Costa Rica using 16S rRNA amplicon sequencing. Lithobates vibicarius, Craugastor escoces, and Isthmohyla rivularis have recently been rediscovered, whereas Isthmohyla pseudopuma has suffered population fluctuations but has never disappeared. To investigate the life stage at which the protective skin microbiome is shaped and when shifts occur in the diversity of putatively anti-Bd bacteria, we studied the skin microbiome of tadpoles, juveniles and adults of L. vibicarius. We show that the skin bacterial composition of sympatric species and hosts with distinct Bd-infection statuses differs at the phyla, family, and genus level. We detected 94 amplicon sequence variants (ASVs) with putative anti-Bd activity pertaining to distinct bacterial taxa, e.g., Pseudomonas spp., Acinetobacter johnsonii, and Stenotrophomonas maltophilia. Bd-uninfected L. vibicarius harbored 79% more putatively anti-Bd ASVs than Bd-infected individuals. Although microbiome composition and structure differed across life stages, the diversity of putative anti-Bd bacteria was similar between pre- and post-metamorphic stages of L. vibicarius. Despite low sample size, our results support the idea that the skin microbiome is dynamic and protects against ongoing Bd presence in endangered species persisting after their presumed extinction. Our study serves as a baseline to understand the microbial patterns in species of high conservation value. Identification of microbial signatures linked to variation in disease susceptibility might, therefore, inform mitigation strategies for combating the global decline of amphibians.


INTRODUCTION
All vertebrates harbor a diverse and dynamic microbial community on their skin -the so-called skin microbiomewith which they have co-evolved (McFall-Ngai et al., 2013;Pita et al., 2018). The skin microbiome is considered an important component of the host's innate immune system, being one of the first lines of defense against pathogenic infections and the mediator of disease susceptibility (Harris et al., 2006;King et al., 2016). In vertebrates, adequate protection against pathogens has been linked with distinct characteristics of the microbiome, such as bacterial species richness, microbial community assemblage, and the presence of bacteria that produce certain antimicrobial compounds that inhibit pathogen growth (Harris et al., 2006;Lemieux-Labonté et al., 2017;Ellison et al., 2018).
In amphibians, the skin microbiome can vary across species and populations, being shaped, for example, by ontogeny, geography, and seasonality (Kueneman et al., 2014;Longo et al., 2015;Bletz et al., 2017;Prest et al., 2018). Notably, it can protect against the infection of the fungal pathogen Batrachochytrium dendrobatidis (hereafter Bd), which causes an infectious skin-disease called chytridiomycosis that can kill amphibians through the production of noxious compounds and disruption of cutaneous osmoregulation (Berger et al., 1998;Voyles et al., 2009). In particular, skin bacteria can inhibit Bd colonization and pathogenicity in amphibians (Harris et al., 2009;Muletz-Wolz et al., 2017a). Previous research has shown that the presence of Bd-inhibitory bacteria (e.g., Acinetobacter, Lysobacter, Janthinobacterium) provides a protective function against Bd in certain hosts (Brucker et al., 2008;Woodhams et al., 2018;Rebollar et al., 2019). Characterization of these microbiome patterns across hosts should guide us to a better understanding of co-evolutionary processes and may provide means to protect amphibians against pathogens such as Bd.
Amphibians, especially those in the Neotropics, have experienced severe population declines and extinctions since the 1980s (Whitfield et al., 2007(Whitfield et al., , 2016. In Costa Rica, one of the most amphibian-diverse countries in the world, numerous amphibian populations from the highlands drastically declined during the 1980s and 1990s, to the extent of the possible extirpation of entire species (Bolaños, 2009). A well-known example is the enigmatic extinction of the endemic golden toad (Incilius periglenes) from the protected Monteverde cloud forest (Crump et al., 1992;Pounds et al., 1997). This horrifying situation has been mainly attributed to outbreaks of Bd-disease together with factors such as habitat alteration, environmental contaminants, and climate change (Collins and Storfer, 2003;Skerratt et al., 2007;Whitfield et al., 2016). Molecular evidence suggests that Bd was present in Costa Rica prior to the observed declines, and distribution models have shown that Bd has spread around the country to the detriment of many amphibian populations (Puschendorf et al., 2006(Puschendorf et al., , 2009De León et al., 2018).
Some species that suffered sharp reductions in the abundance of their populations have since shown signs of recovery. Others that were presumed extinct after the extirpation of their known populations, surprisingly, have been rediscovered years later in small and isolated populations (so-called "relict populations"), despite Bd still being present (Pounds et al., 1997;González-Maya et al., 2013;Chaves et al., 2014;Whitfield et al., 2017). Recovering and relict populations seem to have evolved resistance/tolerance to Bd. A protective skin microbiome has been suggested to be one of the most important factors providing this resistance to Bd (Woodhams et al., 2007b(Woodhams et al., , 2016Grogan et al., 2018). Among these resistant populations are four sympatric amphibian species from the highlands of Costa Rica above 1000 m.a.s.l. Three of them, namely Lithobates vibicarius, Craugastor escoces, and Isthmohyla rivularis, persist despite Bd presence after their presumed disappearance/extinction. L. vibicarius declined drastically and apparently disappeared across its range between the late 1980s and mid-1990s, but surprisingly, since 2002, relict populations have been re-encountered in Costa Rica at a few sites coexisting with Bd (Castro-Cruz and García-Fernández, 2012;Whitfield et al., 2016) with relatively stable numbers of adults. C. escoces, which was once abundant across its distribution range, was declared "Extinct" by the IUCN, but a single individual was recently found by our research group after 30 years without sightings and 12 years of being declared extinct (Jiménez and Alvarado, 2017). It is noteworthy that its last detection had preceded the last sighting of the famous golden toad, I. periglenes. Two years of continuous surveys following the rediscovery of the single female have revealed only one new adult (Douglas Vargas, personal communication). The frog I. rivularis also declined in the late 1980s, but since 2007 has occasionally been reported in distinct locations, such as the Monteverde cloud forest (Solís et al., 2010) and our study area, the Juan Castro Blanco National Park (Jiménez et al., 2019). The other species, Isthmohyla pseudopuma, exists with apparently stable populations, although fluctuations in abundance were observed before and after the declines in the populations from Monteverde (Pounds et al., 1997). These Bdsurvivor species, particularly those that were presumed extinct, provide a unique opportunity to improve our understanding of the selection and diversity of Bd-protective microbiomes in neotropical amphibians.
In the present study, we examined the four above-mentioned sympatric frog species to describe the variation in their skin microbiomes and investigate the prevalence of putatively anti-Bd bacterial members according to Bd-infection status. Previous research has shown that the amphibian skin microbiome is species-specific because of the traits with which it associates, such as skin compounds (e.g., antimicrobial peptides), ecology, and behavior (McKenzie et al., 2011;Kueneman et al., 2014;Belden et al., 2015). Therefore, we hypothesized that skin bacterial communities are distinct among adults of different study species. Furthermore, recent evidence has shown distinct microbiome diversity between Bd infected/uninfected amphibians (Ellison et al., 2018), as well as direct interactions between Bd and skin microbes, whereby some bacterial taxa can inhibit Bd growth, colonization, and establishment on the skin of certain hosts (Harris et al., 2009;Park et al., 2014;Becker et al., 2015;Muletz-Wolz et al., 2017a,b). Consequently, if the skin microbiome plays a significant role in Bd-resistance, we would expect a different microbiome between Bd-infected and Bd-uninfected adults, and Bd-uninfected individuals should harbor a higher diversity of putatively anti-Bd bacterial members.
In order to determine the life stage at which the protective skin microbiome is shaped, as well as when shifts occur in the abundance of putatively anti-Bd bacteria, we focused on L. vibicarius. It is well known that dramatic changes occur during amphibian metamorphosis, such as the keratinization of skin tissue and the development of the immune system (Wells, 2007;Rollins-Smith, 2009). Keratinized tissue in amphibians is restricted to the mouthparts of tadpoles and the skin of juveniles and adults (Alibardi, 2001;Wells, 2007). The mucus glands present before and after metamorphosis produce distinct secretions that act as a physical protective barrier and substrate for skin microbial communities, which limit and/or promote the growth of certain bacterial taxa (Rollins-Smith, 2009;Woodhams et al., 2014). Bd infects both larval and post-metamorphic stages but it usually does not produce lethal infections until after metamorphosis, when skin becomes keratinized (Blaustein et al., 2005;Langhammer et al., 2014). Thus, we hypothesized that the skin microbiome differs between the life stages of L. vibicarius. Given that juveniles and adults are usually more susceptible to Bd than tadpoles (Blaustein et al., 2005;Langhammer et al., 2014) presumably due to increased keratin distribution and abundance after the larval stage -we expected an increase in the diversity of putatively anti-Bd bacteria after metamorphosis in L. vibicarius.

Study Sites
We surveyed four study sites in the Juan Castro Blanco National Park and adjacent private lands within Alajuela province, Costa Rica (10 • 16 44.4 N, 84 • 19 58.9 W) for the presence of amphibians over three consecutive years [2015 (year 1), 2016 (year 2), and 2017 (year 3)] during the same time period within the rainy season (September-November). We have not provided the exact coordinates of the study sites to discourage people from visiting them, because we hope to prevent the potential introduction of pathogens and illegal collecting.
Three sites (Lagunillas, Monjes, and Congo) are disturbed landscapes adjacent to the national park. Each of these sites contains a pond and small streams within pastures for cattle and are surrounded by secondary and fragmented montane humid forest. The presence of cow feces, the antibiotic treatment of cows, and the use of pesticides (e.g., herbicides) at these sites might negatively affect habitat quality for local amphibians. The national park site (hereafter NP) includes three closely located large ponds surrounded by late successional and continuous montane humid forest. The altitudes of the study sites range from 1830 to 1900 m.a.s.l.
All four study sites have been monitored since 2015 for the presence of Bd by our research group (see below). L. vibicarius ceased to occur at these sites in the late 1980s (Douglas Vargas, personal communication), coinciding with the amphibian population declines and extirpations thought to be associated with Bd during the late 1990s.

Focal Species
We detected four amphibian species in study sites ( Figure 1A). The most abundant species, L, vibicarius (IUCN status: Vulnerable), is a semi-aquatic montane ranid frog formerly common throughout the mountain ranges of Tilarán, Cordillera Central, and Talamanca in Costa Rica and western Panama (Savage, 2002). This species forages in dense montane forests and breeds in ponds and shallow bodies of water located in closed and open areas (Savage, 2002). The IUCN status of the species was changed to "Critically Endangered" and "Possibly Extinct" after its population declines and disappearance across the distribution range in the late 1990s (IUCN SSC Amphibian Specialist Group and NatureServe, 2013). The causes of its apparent disappearance are still not clear; however, the declines coincided with the arrival of Bd, suggesting Bd as one of the main drivers of the declines. Six years after disappearing, the species was re-encountered at a few sites in Costa Rica and reproduced in high numbers despite its co-existence with Bd (Whitfield et al., 2017). Since 2013, its IUCN status has been changed to "Vulnerable." Craugastor escoces [IUCN status: "Extinct"; rediscovered in 2016 (Jiménez and Alvarado, 2017)] is a stream-breeding endemic frog that used to be extremely common throughout the Cordillera Central mountain range in Costa Rica (Savage, 2002). The cause of its presumed extinction is not fully understood; however, chytridiomycosis is suspected to be one of the main causes (Bolaños and Chaves, 2004). This species belongs to the Central American Craugastor punctariolus group (Hedges et al., 2008) that exemplifies the severe population declines of montane riparian amphibians throughout Central America. Interestingly, the decline and extirpation of another species from this group, C. punctariolus, has been associated with a Bd-outbreak (Ryan et al., 2008).
Isthmohyla rivularis (IUCN status: "Critically Endangered") is an arboreal frog associated with fast-flowing highland streams (Savage, 2002). The species was formerly common throughout the mountain ranges of Tilarán, Cordillera Central, and Talamanca in Costa Rica and western Panama, but populations were extirpated. I. rivularis was presumed to be extinct until 2007 when it was rediscovered in the Monteverde cloud forest (Wheelwright, 2000). It was recently found in our study site as well (Jiménez et al., 2019). The cause of the decline of this species is not known; however, chytridiomycosis has been suggested as a possible factor involved in population extirpations (Solís et al., 2010).
Lastly, I. pseudopuma (IUCN status: "Least Concern") is an arboreal frog that can be found in both disturbed and undisturbed habitats, and its distribution is similar to that of L. vibicarius. It is an explosive breeder, congregating at permanent and temporal ponds in the highlands (Savage, 2002). Individuals infected with Bd have been identified (Puschendorf et al., 2009), but the species seems to have evaded population declines (Pounds et al., 1997).

Sampling and Molecular Detection of Bd
In order to describe the variation in the skin microbiome of sympatric species and specifically to investigate the prevalence of bacterial genera with putatively anti-Bd members according to Bd-infection status, we sampled adults of L vibicarius, I. pseudopuma, I. rivularis, and C. escoces ("species data set"). In order to understand the life stage at which (1) the protective skin microbiome is shaped and (2) shifts occur in the abundance of putatively Bd-inhibitory bacteria, we sampled tadpoles, juveniles, and adults of L. vibicarius ("lifestage data set"). We obtained unequal sample sizes among study sites for each life stage due to adverse conditions during fieldwork, with either few rainy days to stimulate the frogs' activity or difficulty accessing study sites after heavy rains. Sample sizes for microbiome analyses are given in Supplementary Table S1.
Adult and juvenile frogs were captured and handled with clean gloves, and each frog was kept in an individual clean plastic bag or small container. Tadpoles of similar size and developmental stage (Gosner stage 26-29) were captured with a dip net. Swabbing consisted of moving a sterile rayon-tipped swab (Peel Pouch Dryswab TM Fine Tip) across the animal's skin following a standard protocol. The swabbing protocol for adults and juveniles consisted of 10 strokes along the ventral side, five strokes along each thigh, and five strokes across each foot. For tadpoles, the protocol consisted of 10 strokes on each side (along body and tail), 10 strokes on the dorsal surface of the body, and 10 strokes on the mouth. Before swabbing, we rinsed every individual with distilled autoclaved water to remove sediment, dirt, and transient bacteria. We swabbed individuals twice using new swabs each time; the first swab was used for identifying skin bacteria and the second to test for Bd infection (Familiar López et al., 2017). Swabs for bacterial analysis were stored in 1.5 ml sterile Eppendorf tubes with 300 µl DNA/RNA Shield (Zymo Research), and swabs for Bd analysis were placed directly in dry 1.5 ml sterile Eppendorf tubes. All swabs were stored at −20 • C until further processing.
To test for the presence of Bd in individual samples, we followed previously described qPCR protocols (Boyle et al., 2004;Whitfield et al., 2013). We extracted DNA from swabs using 50 mL of PrepMan Ultra (Applied Biosystems, Foster City, CA, United States) and detected the presence of Bd DNA using a TaqMan qPCR assay in a QuantStudio R 3 Real-Time PCR System. Samples were run with one negative control and two positive controls. Positive samples were analyzed in triplicate to verify the result (Whitfield et al., 2013).

Microbiome Sequencing
We extracted bacterial genomic DNA from swabs using standard protocols for the NucleoSpin Soil Kit (Macherey-Nagel). To yield more DNA, we preheated SE buffer to 37 • C and conducted two consecutive elution steps (2 × 45 µl). We followed the 4-primer amplicon tagging scheme of Fluidigm (Access Array System for Illumina Sequencing Systems, Fluidigm Corporation) in which amplicon PCR and barcoding PCR occur simultaneously. For amplification, we used the primers 515F (5 -GTGCCAGCMGCCGCGGTAA-3 ) and 806R (5 -GGACTACHVGGGTWTCTAAT-3 ), which target a 291-bp fragment of the hypervariable V4 region of the 16S rRNA gene. However, these primers had to be adapted according to the Fluidigm protocol, and so tagged target-specific primers (CS1-TS-515F and CS2-TS-806R) were combined with samplespecific primer pairs that contained a barcoding sequence and the adaptor sequences used by the Illumina sequencing systems. Final 15-µl volumes for the target specific 4-primer amplicon tagging reaction contained 10 ng/µl DNA, 1X FastStart PCR grade nucleotide mix buffer without MgCl 2 (Roche), 4.5 mM MgCl 2 (Roche), 200 µM of each PCR grade nucleotide (Roche), 0.05 U/µl FastStart high fidelity enzyme blend (Roche), 1X access array loading reagent (Fluidigm), 400 nM access array barcode primers for Illumina (Fluidigm), 5% DMSO (Roche), 2.4% PCR-certified water, and 50 nM target-specific primers (TS-515F and TS-806R). The barcoded samples were purified by using NucleoMag NGS Beads (Macherey-Nagel) with a 1:1 ratio of amplicons to beads and quantified with PicoGreen on Tecan F200. We then pooled all samples to an equal amount of 12 ng DNA and diluted the pool down to 6 nM. Finally, the library was loaded at 7.5 ppm and sequenced on an Illumina MiSeq instrument by using a 250-bp paired-end strategy with 10% PhiX added to account for low base diversity. Raw sequencing data is available in the Dryad Digital Repository: https://doi.org/10. 5061/dryad.n34035p.

Bioinformatic Analyses
The initial processing of sequencing reads was performed as previously described (Šrut et al., 2019). Sequencing reads were pre-processed with the open-source QIIME2 (version 2019.1) (Bolyen et al., 2018). For dada2 analysis, we trimmed the first bases of each read both to remove primers (-p-trim-leftf 23, -p-trim-left-r 20) and to truncate forward and reverse reads to 200 bp because of decreasing average quality scores of the sequences at the end. Finally, we collapsed reads into representative sequences, known as amplicon sequence variants (ASVs). The ASV approach has been shown to outperform Operational Taxonomic Unit (OTU) clustering; unlike OTU clustering, ASVs are not compared with a reference database and are retained when they meet an arbitrary similarity cutoff (Callahan et al., 2017). We generated a mid-point-rooted bacterial phylogenetic tree for further diversity analyses by using MAFFT (Katoh and Standley, 2013) and constructed a phylogenetic tree using Fast Tree 2 (Price et al., 2010). We assigned taxonomy to the resulting ASVs with the Greengenes database (version 13_8) as a reference 1 . Then, we removed sequences classified as chloroplast, mitochondria, archaea, Eukaryota, and unclassified phylum.
In order to be able to undertake further analyses (i.e., subsequent data filtering and statistical analyses) in the R environment (R Core Team, 2018), the QIIME2 artifacts "feature-table.qza, " "phylogenetic-tree.qza, " and "taxonomy.qza" were converted to .biom, .tsv, and .nwk files, respectively. Finally, we created a new feature-table.biom file (feature-tabletaxonomy.biom) by adding the taxonomy information from the "taxonomy.tsv" file to the "feature-table.biom" file.
We used the R package "phyloseq" (McMurdie and Holmes, 2013) to import and merge the "feature-table-taxonomy.biom" file, the bacterial phylogenetic tree, and a text file containing the metadata. We prepared separate datasets according to the analysis, i.e., a "species dataset" and a "life-stage dataset." In downstream analyses, we excluded ASVs that contained fewer than 50 reads. We also excluded samples with less than 9000 sequences, except for four samples in the "species dataset" [one C. escoces sample (4026 reads) and three L. vibicarius samples (7081, 8416 and 8999 reads)] because of their relevance for the study. The number of reads per sample ranged from 4026 to 52,154 in the "species dataset" and 9195 and 53,957 in the "life-stage dataset." For alpha and beta diversity analyses, we normalized read counts across samples by rarefying datasets according to the sample with the lowest read number using the R package "phyloseq."

Data Visualization and Statistical Analysis in the R Environment
In order to visualize differences in the relative abundance of skin bacterial taxa, we generated stacked bar charts at the phylum, family, and genus levels. These charts were produced using the R packages "phyloseq" and "ggplot2" (Wickham, 2016).
We calculated two alpha diversity indices, "number of observed ASVs" and "Shannon Index" (Shannon, 1948;Spellerberg and Fedor, 2003), for each skin microbiome sample using the R package "phyloseq." The Shannon Index is a metric that weights the number of observed ASVs by their relative evenness across the community (Reese and Dunn, 2018).
To test the effect of life stage on the number of observed ASVs and Shannon Index, we performed Generalized Linear Models (GLMs) with Gamma and Gaussian distributions, respectively. We also included site and year as covariates in the models. Both alpha diversity measures were log-transformed prior to analysis to improve normality. We used the odds ratio (OR) and 95% confidence intervals (95% CI) (Nakagawa and Cuthill, 2007), estimated with the R package "emmeans" (Lenth, 2019), to quantify the deviation (i.e., effect size) of the alpha diversity measures between life stages. In pairwise comparisons, an overlap between the 95% CI of the OR estimates and the value "1" indicates that there is no difference.
We assessed the beta diversity of skin bacterial communities across life stages using the unweighted UniFrac distance (based on presence-absence and phylogenetic distance) and the weighted UniFrac distance (based on abundance and phylogenetic distance). These distance matrices were calculated using the R package "phyloseq." Because these diversity metrics are sensitive to ASVs with very low abundance, we removed low abundance ASVs from the rarefied data by filtering out those with a prevalence lower than 10 (defined as the number of samples in which a taxon appears), and then calculated the distance matrices. We fitted permutational multivariate analyses of variance (PERMANOVAs) using the adonis function of the R package "vegan" (Oksanen et al., 2014) to evaluate the variation of the skin bacterial community composition and structure across life stages. We incorporated site and year in PERMANOVAs to account for their variation. We noted the amount of variation explained by the covariates using the R 2 from the output of the adonis function. In addition, pairwise PERMANOVAs with Bonferroni corrections were performed to further investigate the effect of life stage on skin bacterial beta diversity. Furthermore, we quantified the extent of the difference in the distances between group centroids (life stages) as a measure of effect size by calculating Cohen's d and 95% CI with the R package "compute.es" (Del Re, 2013). If the 95% CI of Cohen's d estimates overlap with "0, " differences are not considered significant in pairwise comparisons. We performed principal coordinate analysis (PCoA) to visualize the beta diversity distances between life stages.
DESeq2 within the R package "phyloseq" was used to identify bacterial taxa (at family and genus levels) that significantly differed in abundance between L. vibicarius life stages (adjusted p-value < 0.01). We calculated the geometric means of bacterial abundance counts before analysis using the unrarefied dataset. DESeq2 fits negative binomial GLMs for group comparisons and accounts for variation in library size between samples. We used the tax_glom function of the R package "phyloseq" to group ASVs at family and genus levels.
Among the four species as well as among L. vibicarius life stages, we explored the presence and abundance of putative anti-Bd bacterial ASVs. To do so, we filtered the ASVs table to retain only the reads that matched with members that are known to have anti-Bd properties using a database of culturable anti-Bd bacteria (Antifungal Isolates Database; Woodhams et al., 2015) as a reference. This database contains 16S rRNA gene sequences isolated from amphibian skin. We retained ASVs with 100% sequence identity match to those in the mentioned database following the methods outlined by Muletz-Wolz et al. (2019). We caution that this does not necessarily mean that these bacterial strains possess anti-Bd properties, rather that they are strong candidates. The use of this database to detect putatively anti-Bd bacterial strains has been used in previous amphibian microbiome studies (e.g., Kueneman et al., 2015;Abarca et al., 2018;Muletz-Wolz et al., 2019). We visualized the presence and relative abundance patterns with a Venn diagram and heatmaps using the R packages "Vennerable" and "ggplot2, " respectively.

Variation in the Skin Microbiome and Prevalence of Putatively Anti-Bd ASVs in Adults of Four Sympatric Species
Amplicon sequence variants of the phylum Proteobacteria dominated the skin microbiome of both Bd-and Bd+ L. vibicarius adult frogs (50 and 40%, respectively) and of Bd+ I. rivularis (82%) and Bd-C. escoces (48%) (Figure 1B). In contrast, ASVs of the phylum Bacteroidetes were most abundant on the skin of Bd-I. pseudopuma (61%) (Figure 1B).
Among the four sympatric frog species, we detected a total of 94 putative anti-Bd ASVs (Supplementary Table S2), from which 77.7 and 31.9% were identified at genus and species levels, respectively. The Bd-L. vibicarius showed the highest number of putatively anti-Bd ASVs among all samples (Supplementary Table S2). Putatively anti-Bd ASVs assigned to the families Enterobacteriaceae and Comamonadaceae, as well as to the genera Pseudomonas and Stenotrophomonas, were the most frequent and abundant on the skin of both Bd+ and Bd-frogs (Figure 2). Specifically, in L. vibicarius, we observed 91 (96.80%) of these putatively anti-Bd ASVs in the skin microbiome of Bd-individuals and 19 (20.21%) in Bd+ individuals (Supplementary Table S2). Interestingly, putative anti-Bd ASVs assigned to Janthinobacterium lividum, Flavobacterium succinicans, A. johnsonii, and Serratia marcescens were present only on Bd-L. vibicarius (Figure 2). On the skin of both Bd-I. pseudopuma and Bd+ I. rivularis, we observed 17 (18.08%) of these associated potential Bd-protective bacterial taxa (Supplementary Table S2). In Bd-I. pseudopuma, putatively anti-Bd ASVs assigned to Chryseobacterium spp. showed the highest abundance, whereas on Bd+ I. rivularis, potential Bdprotective ASVs assigned to Pseudomonas spp. and P. veronii and A. johnsonii were the most abundant (Figure 2). Strikingly, the skin of Bd-C. escoces harbored a relatively low number of putatively anti-Bd ASVs, with only eight (8.50%) taxa detected (Supplementary Table S2), and those assigned to the family Aeromonadaceae and the genus Pseudomonas were detected with the highest abundance (Figure 2).

Shifts in the Skin Microbiome and Putatively Anti-Bd ASVs Across Life Stages of L. vibicarius
The skin microbiomes of L. vibicarius tadpoles, juveniles, and adults were mainly dominated by ASVs from the phylum Proteobacteria (42, 62, and 50%, respectively). Moreover, ASVs from the phylum Firmicutes and Bacteroidetes were highly abundant (>10%) in both tadpoles and adults but were less frequent in juveniles (Figure 3A). At the family level, Comamonadaceae was the most prominent in the microbiome of tadpole skin (17%), whereas Pseudomonadaceae and Alcaligenaceae were most abundant in juveniles (9 and 7.2%, respectively). In the microbiome of adult skin, the most abundant identified family was Alcaligenaceae (17%) followed by Bacteroidaceae (5.5%) ( Figure 3B).
With regard to the alpha diversity of the skin microbiome, we found a strong influence of life stage on the number of observed ASVs (p < 0.001) (Supplementary Figure S2A), but not on the Shannon index (p = 0.95) (Supplementary Figure  S2B). The number of observed ASVs was higher in tadpoles than in juveniles (OR = 1.54, 95% CIs not overlapping 1) and adults (OR = 1.43, 95% CIs not overlapping 1), whereas juveniles and adults showed a similar diversity measure (OR = 0.93, 95% CIs overlapping 1). Both alpha diversity measurements were influenced by year, but not by the study site (Table 1).
Concerning beta diversity, the PERMANOVA models revealed an effect of life stage on the bacterial community composition (Unweighted UniFrac: R 2 = 0.20, p = 0.001) and abundance pattern (Weighted UniFrac: R 2 = 0.17, p = 0.001) ( Table 2). Moreover, the study site and year explained part of the   variation in the skin bacterial communities ( Table 2). Pairwise PERMANOVA tests and Cohen's d effect sizes showed that the R 2 and effect sizes between juveniles and adults were much lower, indicating smaller differences between these stages ( Table 3). The effects of life stage on the beta diversity indices are visualized by PCoA ordination plots (Figure 4). DESeq2 analyses revealed important differences in the prevalence and abundance pattern of bacterial families and genera. Families detected in different life stages were in most cases significantly more abundant in tadpoles (Figure 5). In particular, we detected 34 families with significantly different abundances in the skin microbiome of tadpoles and juveniles ( Figure 5A). Of these, 31 had a higher abundance in tadpoles, whereas only three exhibited this difference in juveniles. In tadpoles, this aspect especially concerned the families Coriobacteriaceae (log2 fold change: −24.83), Deferribacteraceae (−12.79), Dethiosulfovibrionaceae (−12.29), and S24-7 (−11.97) and, in juveniles, the families Sphingobacteriaceae (4.52), Alcaligenaceae (2.75), and Methylocystaceae (1.91) (Figure 5A). Comparing the skin microbiome of tadpoles and adults, we identified 17 families with higher abundances in tadpoles (e.g., Dethiosulfovibrionaceae, Deferribacteraceae, Fusobacteriaceae, and RFP12), whereas adults showed 12 families with higher abundances, such as [Odoribacteraceae] (log2 fold  Figure 5B). Lastly, we observed only the bacterial family Oxalobacteraceae (log2 fold change: −2.62) as being more common in juveniles than in adults, whereas 15 families were more abundant in adults, such as Bacteroidaceae (29.02) and Odoribacteraceae (26.38) (Figure 5C). When looking at the genus level, we identified 27 genera with significantly different abundance between tadpoles and juveniles ( Figure 6A). Of these, 26 genera revealed a higher abundance in tadpoles and only one in juveniles. For instance, the genera Bacteroides (log2 fold change: −27.0), Bilophila (−26.32), Parabacteroides (−25.45), and Bifidobacterium (−25.12) showed a very high abundance in tadpoles relative to juveniles ( Figure 6A). Furthermore, we identified 18 genera that were more abundant in tadpoles (e.g., Clostridium and Mucispirillum) than in adults. However, adults showed a higher abundance of 15 genera, such as Odoribacter (log2 fold change: 27.0), Butyricimonas (24.82), and Eubacterium (27.46) ( Figure 6B). Moreover, in tadpoles, we found one genus (Stenotrophomonas) whose members have been previously associated with possessing anti-Bd activity that differed greatly in abundance relative to adults, whereas four genera with anti-Bd activity members were detected in adults (Chryseobacterium, Hymenobacter, Erwinia, and Citrobacter). The juvenile skin microbiome exhibited a significantly higher abundance of only two genera, namely Arthrobacter (log2 fold change: −25.19) and Chromobacterium (−9.87), relative to adults, whereas adults harbored higher abundances of 11 genera, such as Akkermansia (26.80) and Parabacteroides (26.78) (Figure 6C).
Putatively anti-Bd bacterial ASVs were abundant in both pre-and post-metamorphic stages of L. vibicarius. We detected 104 ASVs associated with inhibiting Bd-growth across all life stages (Figure 7 and Supplementary Table S3). We found 74 of these in the skin microbiome of tadpoles, whereas 58 and 94 were observed in the skin microbiomes of juveniles and adults, respectively (Figure 7 and Supplementary Table S3). Interestingly, 46 ASVs (48%) with putative anti-Bd activity were present in all three life stages, and tadpoles and adults had the highest overlap (66 ASVs, 69%) (Figure 7). The putatively anti-Bd ASVs assigned to the family Enterobacteriaceae, the genera Pseudomonas, and the species P. veronii and A. johnsonii were present in all life stages with the highest abundance (Figure 8).

DISCUSSION
Despite growing evidence that cutaneous bacteria play an important role in amphibian health (reviewed by Jiménez and Sommer, 2017), our current understanding of the microbiome in amphibians from tropical systems remains limited. This is the first study using high-throughput sequencing to investigate the skin microbiome of amphibians that persist in co-existence with Bd after their presumed disappearance/extinction from the highlands of Costa Rica. We hypothesized that (1) skin bacterial communities would differ between co-occurring species and would be distinct in Bd-infected and Bd-uninfected hosts, (2) Bd-uninfected individuals would harbor a higher diversity of putatively Bd-inhibitory bacterial ASVs, (3) the skin microbiome would differ between life stages of L. vibicarius, and (4) the diversity of putatively Bd-inhibitory bacteria would be higher after the metamorphosis of L. vibicarius. With this study, we have begun to fill the gap in our knowledge of the skin microbiome of species with relict populations and suspected resistance to Bd.
As hypothesized, we found distinctly different skin microbiomes between adults of co-habiting species. We observed differences in the composition of ASVs at the phyla, family, and genus levels. For instance, with regard to the microbiome across our Bd-amphibians, the most dominant family on the skin of L. vibicarius is Alcaligenaceae, whereas on I. pseudopuma and C. escoces, the most prominent families are Weeksallaceae and Planococcaceae, respectively. Moreover, the genus Chrysobacterium is present in higher proportions in adult. Bacterial families with a log2 fold change less than 0 were more abundant in the life stage indicated on the left, whereas those with a log2 fold change higher than 0 were more abundant in the indicated life stage on the right. Bacterial families are colored by phylum and sized by mean relative abundance across samples. adult. Bacterial genera with a log2 fold change less than 0 were more abundant in the life stage indicated on the left, whereas those with a log2 fold change higher than 0 were more abundant in the indicated life stage on the right. Bacterial genera are colored by phylum and sized by mean relative abundance across samples. I. pseudopuma than in L. vibicarius and C. escoces. Interestingly, we observed a higher proportion of ASVs of the genus Bacteroides in L. vibicarius compared to the other species. This genus plays an important role in maintaining gut homeostasis (Wexler and Goodman, 2017) and therefore could be relevant for supporting the health of L. vibicarius. This species-specific microbiome pattern is congruent with other studies in the tropics, such as in low-and mid-elevations of Costa Rica and Panama Rebollar et al., 2016b;Abarca et al., 2018) and in the temperate region (McKenzie et al., 2011;Kueneman et al., 2014;Bletz et al., 2017). The family Alcaligenaceae and genus Chryseobacterium were also the most common taxa in the bacterial communities of Agalychnis callidryas, Leptodactylus savagei, and Pristimantis ridens from the lowlands of Costa Rica (Abarca et al., 2018). In this study, we have not evaluated the mechanisms driving variation in the skin microbiome; however, we suspect that a distinct antimicrobial peptide composition among our species, as has been observed in other amphibians, plays a role as a selective force that determines which microbial genotypes can grow on each host's skin (Woodhams et al., 2007a;Franzenburg et al., 2013). This could be driven by differences in ecological and behavioral traits across species (e.g., forest forager vs. arboreal and stream forager). Such differences lead to exposure to distinct environmental conditions (e.g., pH, water temperature), which may result in the stimulation or inhibition of the growth of certain bacterial strains on the skin. Furthermore, the presence of different environmental microbes across the microhabitats used by the species could drive distinct microbial colonization of the skin (Loudon et al., 2014;Bletz et al., 2017).
Considering the Bd-infection status, we observed, as expected, a different bacterial composition at the phylum and family levels between Bd infected/uninfected adults of the same species, for example in Bd-versus Bd+ of L. vibicarius. We also observed a higher diversity of ASVs with putatively anti-Bd activity on the skin of L. vibicarius Bd-individuals than on the skin of Bd+ individuals. More than half of these ASVs are absent in Bd+ individuals, and most of the shared bacterial taxa, such as P. veronii, are more prevalent on the skin of Bd-individuals than on Bd+ individuals. The bacterial composition and higher diversity of putatively beneficial bacteria detected on the skin of Bd-L. vibicarius suggests that the bacterial community can promote the health of Bd-individuals. This observation might be linked to a reduction of the bacterial richness on the skin of highly Bd-infected amphibians, as observed in a previous study (Ellison et al., 2018). Further research is needed to confirm this relationship. We should also consider that the observed differences in the skin microbiome might be a direct effect of Bd, as found in previous studies on the skin microbiome of Rana sierrae and Lithobates catesbeianus (Jani and Briggs, 2014;Walke et al., 2015;Jani et al., 2017).
Furthermore, we found bacterial ASVs with putatively anti-Bd activity in all four species, regardless of their Bd-infection status. Putatively anti-Bd ASVs assigned to the genera Pseudomonas and Chryseobacterium are among the most abundant Bd-protective bacteria on Bd-I. pseudopuma and Bd-C. escoces, respectively. These two genera have been reported with a high proportion of isolates inhibiting Bd-growth . Notably, we observed less than half of these beneficial bacteria (e.g., putatively anti-Bd ASVs assigned to Pseudomonas spp. and A. johnsonii) from the overall Bd-protective bacteria detected across our samples on the skin of Bd+ I. rivularis. These ASVs might limit the impacts of Bd-infection, given that no signs of disease have (so far) been observed in this species. We also detected potentially beneficial ASVs assigned to the genus Janthinobacterium on Bd-L. vibicarius, Bd-I. pseudopuma, Bd-C. escoces, and Bd+ I. rivularis. One species in this genus, J. lividum, detected in our study species, is known to inhibit Bd efficiently on the salamander Plethodon cinereus (Brucker et al., 2008), although this inhibition appears not to be effective in the frog Atelopus zeteki (Becker et al., 2011). Given the devastating impact of Bd on our endangered focal species [and also in other populations in Costa Rica (Wheelwright, 2000;Whitfield et al., 2017)], these findings support the hypothesis that Bd-persistent species possess a skin microbiome with potential for protecting against Bd. Because Bd-C. escoces has a low percentage of ASVs with putatively anti-Bd activity, one can speculate that this species has a skin microbiome with a lower ability to fight Bd. Notwithstanding, this microbiome conformation might be sufficient to protect the species in conjunction with other mechanisms, such as antimicrobial peptides (Woodhams et al., 2007a). Moreover, the protective role of skin bacteria through the production of antifungal metabolites can be a consequence of cooperation and competition among bacteria (Foster et al., 2017). Nevertheless, because C. escoces is a critically endangered species and might remain on the brink of extinction, further research is urgently required to improve our understanding of the species' skin microbiome and its association to Bd-resistance. The use of culture-independent techniques, like the 16S rRNA amplicon sequencing used in this study, to investigate bacterial communities has led to an important improvement in our understanding of amphibian FIGURE 8 | Heatmap of putative Bd-inhibitory ASVs across life stages of L. vibicarius. Reddish colors indicate increasing abundance, yellowish colors indicate lower abundance, whereas white indicates absence of ASVs. The highest possible taxonomic assignment (maximal at the species level) is shown for these ASVs. The numbers indicate the total number of reads of ASVs assigned to specific taxa. G, genus; F, family. microbial ecology and its interaction with Bd (Rebollar et al., 2016a;Jiménez and Sommer, 2017). However, in the light of the high variability of Bd-inhibitory functions among genera, within bacterial isolates from a single genus, and even at the strain level , further inhibition tests with culture-dependent techniques are necessary (1) to confirm anti-Bd activity of these beneficial bacteria among distinct hosts and (2) to reach more robust conclusions regarding the protective role of these bacteria against Bd-infection. For instance, in three Neotropical frogs (A. callidryas, Dendropsophus ebraccatus, and Craugastor fitzingeri) the majority of the culture-dependent skin bacteria represented in the culture-independent community were Bd-inhibitory (Rebollar et al., 2019). Also, further studies are required to evaluate the degree of protection given by microbiomes in our focal species.
In order to understand the life stage at which the protective skin microbiome is shaped and when shifts occur in abundance of putatively Bd-inhibitory bacteria, we sampled tadpoles, juveniles, and adults of L. vibicarius. We found that life stage is a strong predictor of both alpha diversity (number of observed ASVs, but not Shannon diversity) and beta diversity (Unweighted and Weighted UniFrac distances) of the skin microbiome of L. vibicarius. We also detected differences in the abundance of certain skin bacterial families and genera across life stages via DESeq2 analysis. Our findings clearly support the hypothesis that amphibian life stage influences the skin microbiome (Kueneman et al., 2015). Previous research in amphibians, e.g., Anaxyrus boreas and Eleutherodactylus coqui, has revealed life stage as one of the most important factors structuring microbial communities on amphibian skin (Kueneman et al., 2014(Kueneman et al., , 2015Longo et al., 2015). Within this study, we have mainly observed strong differences between pre-and postmetamorphic stages, indicating a dramatic restructuring of the skin bacterial community after metamorphosis. We found a large number of observed ASVs in tadpoles relative to juveniles and adults, whereas these last two stages did not differ from one another. Also, we observed different community compositions and abundance patterns between the three life stages. Furthermore, effect sizes showed that differences between tadpoles and the two post-metamorphic stages, namely juveniles and adults, are much stronger than between the two later stages. Indeed, the beta diversity metrics suggest that the core microbiomes of juveniles and adults are more similar to each other than to the microbiome of tadpoles. Despite an incomplete understanding of the mechanisms shaping the skin microbiome across development, the variation in the skin microbiome across life stages is probably attributable to a distinct type of skin; keratinized skin tissue after metamorphosis is likely more suitable for distinct bacterial taxa (Kueneman et al., 2015). Moreover, across amphibian development, the skin's defensive glands might shift the production of specific antimicrobial peptides that influence the growth of distinct bacterial taxa (Rollins-Smith, 2009). Further research should explore the microbiome diversity of L. vibicarius on metamorphs, because metamorphosis represents a switchpoint with many drastic physiological and immunological rearrangements, a period during which amphibians might be highly susceptible to pathogens.
In addition to the strong effect of life stage on beta diversity, we have also observed that the study site and year partly explain the variation in distance metrics. These results are congruent with those observed in earlier studies (McKenzie et al., 2011;Longo et al., 2015;Familiar López et al., 2017). Further studies should better explore the role of environment (i.e., study site) and year in order to unravel the relative importance of these covariates in our study species. In particular, some of our sampling sites are under the constant influence of human activities that could either stress hosts and/or modify environmental bacterial taxa that might colonize amphibian skin (Krynak et al., 2016;Huang et al., 2018). These effects could ultimately impact the skin microbiome homeostasis of the amphibians.
Contrary to our expectations, we did not observe a significant increase in the diversity of putatively anti-Bd bacteria after metamorphosis in L. vibicarius, as both pre-and postmetamorphic stages of L. vibicarius hosted a high number of putatively anti-Bd members. The presence of these beneficial bacteria might, at least in part, explain the persistence of this endangered species, now co-existing with Bd well after its presumed extinction. These anti-Bd members can reach high population densities on amphibian skin (Harris et al., 2009). In addition, we observed that a high percentage of these bacterial members are shared among life stages. These findings suggest that tadpoles, which usually have a low susceptibility to Bd and can harbor non-lethal infections, maintain a diverse assemblage of putatively anti-Bd bacterial strains. Tadpoles may obtain these putatively anti-Bd members through horizontal transmission, and subsequently act as reservoir hosts of many of these bacteria for post-metamorphic life stages with keratinized skin tissues vulnerable to Bd infection (Walke et al., 2011;Rebollar et al., 2016c). This is particularly relevant because post-metamorphic frogs in the current study were colonized by many putative anti-Bd bacterial ASVs (e.g., Pseudomonas spp. and Janthinobacterium spp.) Taxonomically diverse assemblages of anti-Bd bacteria are more effective at suppressing Bd growth than individual bacterial isolates (Piovia-Scott et al., 2017;Bell et al., 2018). We observed a lower number of anti-Bd ASVs on the skin of juveniles than on adults, and some bacterial strains were only detected in adults (e.g., Acinetobacter rhizosphaerae and Pseudomonas fragi). Adaptive immune functions such as the production of skin antimicrobial peptides and alkaloids continue to develop after metamorphosis in association with skin shedding (Rollins-Smith, 2009;Kueneman et al., 2015). These factors may influence the Bd-protective bacterial assemblage on juvenile skin and thus explain why juveniles had less putatively anti-Bd ASVs than adults (and may be more susceptible to Bd). Further work is needed to comprehend bacteria/Bd interaction patterns across developmental stages.
The current study extends our knowledge about the skin microbiome in Neotropical montane amphibians, particularly with regard to those species that are Bd-susceptible and have relict populations persisting and recovering after dramatic Bd-associated declines. Despite a limited sample size, our findings point to the presence of a potential Bd-protective microbiome in these endangered species. Our results may also contribute to the development of innovative conservation strategies, such as anti-Bd probiotic treatments, to protect susceptible Neotropical amphibian species from the devastating effects of Bd. Further research should focus on the genes that characterize microbiome function to advance our understanding of amphibian adaptation and persistence in the environment as well as maintenance of a good health status. Given the observed patterns in skin microbiome diversity across developmental stages and Bd-infection status in L. vibicarius, we suggest that the microbiome of this frog is an outstanding model for future research.