Eria A. Rebollar1*†
Ana Gutiérrez-Preciado2
Cecilia Noecker3
Alexander Eng3
Myra C. Hughey4
Daniel Medina4
Jenifer B. Walke4†
Elhanan Borenstein3,5,6
Roderick V. Jensen4
Lisa K. Belden4,7
Reid N. Harris1,8- 1Department of Biology, James Madison University, Harrisonburg, VA, United States
- 2Unité d'Ecologie, Systématique et Evolution, Université Paris-Sud, Paris, France
- 3Department of Genome Sciences, University of Washington, Seattle, WA, United States
- 4Department of Biological Sciences, Virginia Tech, Blacksburg, VA, United States
- 5Department of Computer Science and Engineering, University of Washington, Seattle, WA, United States
- 6Santa Fe Institute, Santa Fe, NM, United States
- 7Smithsonian Tropical Research Institution, Panama City, Panama
- 8Amphibian Survival Alliance, London, United Kingdom
Skin symbiotic bacteria on amphibians can play a role in protecting their host against pathogens. Chytridiomycosis, the disease caused by Batrachochytrium dendrobatidis, Bd, has caused dramatic population declines and extinctions of amphibians worldwide. Anti-Bd bacteria from amphibian skin have been cultured, and skin bacterial communities have been described through 16S rRNA gene amplicon sequencing. Here, we present a shotgun metagenomic analysis of skin bacterial communities from a Neotropical frog, Craugastor fitzingeri. We sequenced the metagenome of six frogs from two different sites in Panamá: three frogs from Soberanía (Sob), a Bd-endemic site, and three frogs from Serranía del Sapo (Sapo), a Bd-naïve site. We described the taxonomic composition of skin microbiomes and found that Pseudomonas was a major component of these communities. We also identified that Sob communities were enriched in Actinobacteria while Sapo communities were enriched in Gammaproteobacteria. We described gene abundances within the main functional classes and found genes enriched either in Sapo or Sob. We then focused our study on five functional classes of genes: biosynthesis of secondary metabolites, metabolism of terpenoids and polyketides, membrane transport, cellular communication and antimicrobial drug resistance. These gene classes are potentially involved in bacterial communication, bacterial-host and bacterial-pathogen interactions among other functions. We found that C. fitzingeri metagenomes have a wide array of genes that code for secondary metabolites, including antibiotics and bacterial toxins, which may be involved in bacterial communication, but could also have a defensive role against pathogens. Several genes involved in bacterial communication and bacterial-host interactions, such as biofilm formation and bacterial secretion systems were found. We identified specific genes and pathways enriched at the different sites and determined that gene co-occurrence networks differed between sites. Our results suggest that skin microbiomes are composed of distinct bacterial taxa with a wide range of metabolic capabilities involved in bacterial defense and communication. Differences in taxonomic composition and pathway enrichments suggest that skin microbiomes from different sites have unique functional properties. This study strongly supports the need for shotgun metagenomic analyses to describe the functional capacities of skin microbiomes and to tease apart their role in host defense against pathogens.
Introduction
Microbial symbiotic communities are ubiquitous in animals and plants. For decades, most of the studies in animals have focused on insect endosymbionts, which generally involve only a few bacterial species (Dale and Moran, 2006), whereas more recent studies, mainly in humans, have revealed a complex community of microbes coexisting with their host (Turnbaugh et al., 2007; Grice, 2014; Jorth et al., 2014). Many of these symbionts complete a variety of functions for their hosts, such as nutrient acquisition (Dale and Moran, 2006) and pathogen protection (Fraune et al., 2014; Walke and Belden, 2016). Recent studies have determined that the skin of amphibians harbors bacterial communities that are unique to different species (McKenzie et al., 2012; Kueneman et al., 2014; Walke et al., 2014; Belden et al., 2015; Costa et al., 2016; Rebollar et al., 2016; Bletz et al., 2017a) and might play an important role in protecting the host against Batrachochytrium dendrobatidis, Bd (Harris et al., 2009; Becker et al., 2011; Kueneman et al., 2016). This pathogenic fungus causes a skin disease called chytridiomycosis, and it has been considered one of the greatest global threats to amphibian populations (Fisher et al., 2009; Kilpatrick et al., 2010). Field studies have shown that not all species are susceptible to Bd, as some species persist in Bd-endemic regions with no apparent population declines (Crawford et al., 2010; Rebollar et al., 2014; Rodríguez-Brenes et al., 2016). One important component of defense against pathogens that could explain the presence of tolerant and resistant frog species is skin symbiotic bacteria. Hundreds of Bd-inhibitory bacteria have been isolated from the skin of multiple amphibian species from many sites around the world (Antifungal Isolates Database: Woodhams et al., 2015). Moreover, the addition of some of these bacteria (e.g., Janthinobacterium lividum) to captive amphibians protected the hosts from Bd infection and reduced morbidity and mortality (Harris et al., 2009; Becker et al., 2011; Kueneman et al., 2016). In addition, skin bacterial community structure before infection can predict mortality and morbidity after infection (Becker et al., 2015a; Walke et al., 2015).
Even though our knowledge of amphibian skin microbial communities has grown considerably in the past decade, we still lack knowledge of the full range of functional capabilities of these communities and their interactions with their hosts. Most of the studies have focused on describing the skin community via 16S rRNA gene amplicon sequencing (Kueneman et al., 2014; Rebollar et al., 2016; Bletz et al., 2017a) and/or culturing bacteria to characterize their inhibitory capacities against Bd using in vitro challenge assays (Harris et al., 2006; Flechas et al., 2012; Bell et al., 2013; Becker et al., 2015b; Medina et al., 2017). In other cases, functions have been predicted using bacterial data bases (Kueneman et al., 2015; Bletz et al., 2017b) or predictive tools such as PICRUSt (Loudon et al., 2014a). Moreover, changes in skin bacterial community structure associated with Bd infection have been observed in experimental settings (Jani and Briggs, 2014; Becker et al., 2015a; Walke et al., 2015; Longo and Zamudio, 2017); however, we still lack knowledge of the functional changes that may occur in the bacterial community once Bd is present.
Recent studies on several symbiotic systems have identified genes that were originally described in pathogens, which are also important in mutualistic bacterial species, such as bacterial secretion systems (Dale and Moran, 2006; Preston, 2007; Medina and Sachs, 2010) and biofilm formation (Frese et al., 2013; Schmid et al., 2015). In addition, the host can also affect the symbiotic community through the production of molecules, such as antimicrobial peptides (AMPs) (Gallo and Hooper, 2012; Franzenburg et al., 2013). Furthermore, the study of bacterial symbionts on other animal hosts, like amphibians, has revealed symbiotic bacterial interactions that are important in protecting the hosts against pathogens (Walke and Belden, 2016). One mechanism through which bacteria protect their host is through the production of secondary metabolites that have antimicrobial properties (Flórez et al., 2015). In amphibians, some members of the skin microbiota can produce secondary metabolites that inhibit the growth of the fungal pathogen Bd. Examples of these antifungal metabolites are violacein (Brucker et al., 2008b), prodigiosin (Woodhams et al., 2017), tryptophol (Loudon et al., 2014b), indole-3-carboxaldehyde (I3C) (Brucker et al., 2008b) and 2,4-diacetylphloroglucinol (2,4 DAPG) (Brucker et al., 2008a).
Here, we describe the genes involved in bacterial defense and communication that are present in skin microbiomes of the terrestrial Neotropical frog Craugastor fitzingeri using shotgun metagenomics. Our a priori goals were to expand our knowledge of the metabolic capacities that symbiotic bacteria have and to explore the presence of genes that are potentially involved in bacteria-bacteria, bacteria-host and bacteria-pathogen interactions. Previous studies have described the skin bacterial structure of C. fitzingeri across several regions in Panamá (Belden et al., 2015; Rebollar et al., 2016). These studies showed that the skin bacterial community is dominated by the phyla Proteobacteria and Actinobacteria. Specifically, Rebollar et al. (2016) described differences in bacterial OTU relative abundance in frog skin communities between a Bd-endemic site and a Bd-naïve site. The bacterial community structure in the Bd-endemic site was enriched for taxa known to have antifungal properties (e.g., Pseudomonas and members of the Actinobacteria class). These differences may be related to natural selection caused by the presence or absence of Bd in these lowland regions, although other explanations are possible.
In this study, we analyzed samples from three individuals of C. fitzingeri from a Bd-endemic site and three from a Bd-naïve site. These samples were previously analyzed by Rebollar et al. (2016) with 16S rRNA gene amplicon sequencing for bacterial community profiling. We described the taxonomic composition and analyzed the genetic metabolic pathways present in skin microbiomes. Lastly, we compared the bacterial composition and function among frogs from a Bd-endemic site and a Bd-naïve site. We hypothesized that frog skin microbiomes will include genes associated with bacterial communication and bacterial-host interactions, as well as pathways involved in the production of antifungal metabolites and resistance to bacterial toxins, as a result of potential cooperative or competitive interactions within the community. We also hypothesized that taxonomic and functional composition of skin microbiomes will differ between sites with contrasting Bd incidence. A finding of genes associated with bacterial-host interactions would increase our confidence that amphibian skin bacteria are resident species and not transient. Moreover, unraveling the potential functions present in these bacterial communities will advance our knowledge of the interactions occurring among the host, the bacterial symbiotic community and the pathogen, Bd.
Methods
Sample Selection and Molecular Procedures
C. fitzingeri frogs were collected and swabbed from two lowland forests sites in Panamá: Soberanía National Park (Sob), a Bd endemic site, and Serranía del Sapo (Sapo), a Bd naïve site. DNA from these swab samples was extracted and used in a previous study for 16S amplicon sequencing and qPCR detection of Bd (Rebollar et al., 2016), which contains the details of the sample collection and the molecular procedures used to extract the DNA from the swabs. Six DNA samples extracted from this previous study were used to sequence the metagenome as explained in section Shotgun metagenome sequencing. We chose these samples considering their Bd infection status (infected vs. not infected) and collection site (Bd endemic and Bd naïve) (Table S1). Two other samples were sequenced but were not included in further analyses because they could not be properly annotated.
Shotgun Metagenome Sequencing
Six DNA samples (three from Sob and three from Sapo) were used to obtain frog skin shotgun metagenomes (Table S1). The six barcoded samples were randomly distributed on two lanes (three samples on each lane) and were sequenced on HiSeq 2500 (Genomics Sequencing Center, Bioinformatics Institute, Virginia Tech), generating over 500 million, 100 bp, paired end reads (with ~200 bp insert size). The numbers of reads for each of the six samples ranged from 84 million to 124 million reads (Table S1). Metagenomic reads were deposited in the SRA Database (NCBI) with the accession number SRP130893 as part of Bioproject PRJNA429199.
Metagenome Binning for Filtering Out Host Reads
With the purpose of filtering out the host (C. fitzingeri) reads, we assembled the reads for each of the six samples into metagenome contigs using Ray Meta (Boisvert et al., 2012) de novo assembler. Four million reads for each frog sample were then mapped to all the contigs greater than 1,000 bp from all six frogs to determine the relative abundance of each contig. The covariance of these abundances across samples was then used to cluster all contigs with >100 mapped reads into metagenome bins using the PAM k-medoids algorithm in R with k = 10. Contigs in three of these bins or clusters (Clusters 1–3) exhibited similar abundance profiles across the six samples and appeared to be predominantly frog (with consistently low GC content ~42% and best Blastn hits to eukaryotic sequences). We also mapped the reads for each sample to the frog 18S rRNA sequence (assembled from the data) and a database of bacterial 16S rRNA sequences (from www.patricbrc.org) to estimate the relative proportion of eukaryotic (mainly frog) and prokaryotic (mainly bacteria) DNA from each sample. The 18S/16S ratios for the six samples showed good correspondence with the variation in the abundance of the frog DNA clusters. The remaining seven clusters (Cluster 4–10) were composed of contigs belonging to different bacterial groups.
Taxonomic Assignment and Procrustes Analysis
We determined the bacterial species composition of the frog skin samples using Metaphlan (Segata et al., 2012), a method used to profile bacterial communities based on clade-specific marker genes. We determined the relative abundance of bacterial taxa at different taxonomic levels based on the annotated reads that had been previously filtered to eliminate host-derived reads (Table S1).
We used a Procrustes analysis to evaluate how similar the 16S rRNA gene amplicon sequencing data (obtained with QIIME in Rebollar et al., 2016) was to the shotgun metagenome data (obtained with Metaphlan) using the vegan R package (R Core Team, 2015; Oksanen et al., 2016). Procrustes analysis evaluates the congruency between two data sets by the superimposition of principal component analyses (McHardy et al., 2013; Luo et al., 2016). We performed Procrustes analysis on Bray-Curtis distance matrices calculated from the bacterial relative abundance at the genus level in both data sets. We used the PROTEST test function from the vegan R package (R Core Team, 2015; Oksanen et al., 2016), which performs repeated symmetric Procrustes analyses to estimate if the degree of concordance between matrices is greater than expected by random association (Jackson, 1995). Significant p-values below 0.05 indicate that the matrices are more similar than expected by random association.
Read Level Functional Annotation
To annotate the metagenome reads (read level analysis), we assigned reads to the bacterial clusters obtained in section Metagenome binning for filtering out host reads by re-mapping reads from each sample to contigs assembled from the same sample using blastn (Altschul et al., 1990) with an e-value cutoff of 1. This resulted in 8,323,746 reads (1.57% out of the total number of reads) that mapped to contigs with frog-associated cluster assignments, which were discarded.
We assigned KEGG orthology (KO) group (Kanehisa and Goto, 2000) annotations to non-frog-associated reads from the six frog samples using blastx with an e-value cutoff of 1 to map reads to the KEGG gene database downloaded on July 15th, 2013. This produced 192,761,569 (36.36%) reads with best hits to genes with KO annotations (Table S1). We calculated KO abundances using the number of reads assigned to each KO. In the case where a read had best hits to multiple KOs, the read count was evenly distributed among the corresponding KOs. For example, if a read mapped equally well to 3 KOs, each KO received 1/3 of the count. We normalized the abundances of KOs using MUSiCC (Manor and Borenstein, 2015) with the inter-sample correction option, which corrects for biases in quantification using the abundance of universal single-copy KOs.
Descriptive and Comparative Analyses Using Read Level Annotations
To evaluate differences in the functional repertoire among frog skin communities at multiple levels of detail, we aggregated the annotated KO assignments into modules, pathways, and classes based on a custom-curated version of the BRITE hierarchy. We summed the relative abundances of all KOs associated with each pathway, module, or class, following Manor et al. (2016). Pathways and modules were further filtered to verify that downstream analysis considers only bacterial pathways or modules. Specifically, a pathway or module was included in our analysis only if at least 1% of bacterial genomes in KEGG contained at least 1 KO from that pathway or module, and if these bacterial genomes contained at least 5% (20%) of the KOs in the pathway or module, on average.
For each pathway, module or class, normalized abundances from samples from the Sapo region were compared to those from Sob using LefSe Analysis (Linear Discriminant Analysis Effect Size) which incorporates the use of two non-parametric tests (Kruskal-Wallis and Wilcoxon) and a linear discriminant analysis to detect biomarkers or genes that are differentially abundant among the groups tested (Segata et al., 2011).
Contig Level Functional Annotation and Assembly
To determine the taxonomic identity of the genes present on frog skin metagenomes we assembled reads into contigs (contig level analysis), annotated them and determined their taxonomy (Table 1). Reads were assembled into contigs using stringent criteria to facilitate gene prediction. Forward and reverse reads were assembled using MEGAHIT version 1.3.0 (Li et al., 2016) with default parameters, except for a minimum length of 200 bp for the assembled contigs and a starting kmer size of 23 in increasing steps of 10 until reaching a kmer size of 93. Gene prediction was performed on the newly assembled contigs using Prokka (Seemann, 2014). For metabolic and taxonomic classifications of the predicted genes, we used GhostKOALA from KEGG (Kanehisa et al., 2016).
Table 1. Assembly and annotation data of the six skin metagenomes of C. fitzingeri from sites Sapo (N = 3) and Sob (N = 3).
Descriptive and Comparative Analyses Using Contig Level Annotations
Metabolic and taxonomic data inferred from GhostKOALA were visualized in stacked bar charts using ad hoc scripts in R (R Core Team, 2015). We plotted gene counts from five broad functional classes based on KEGG classification: Membrane Transport (MT), Cellular Communication (CC), Metabolism of Terpenoids and Polyketides (MTP), Biosynthesis of Secondary Metabolites (BSM) and Antimicrobial Drug Resistance (ADR). Functional classes or pathways with abundance values < 3 were not plotted. We used LefSe analysis to test for differentially enriched genes (KOs) between sites for all five broad classes and for all the pathways contained within each class (55 pathways in total). Methods on co-occurrence network construction can be found on Supplementary Methods.
Results
Bacterial Composition of C. fitzingeri Skin Microbiomes
We determined the taxonomic composition and relative abundance of the bacteria present in the frog metagenomes from C. fitzingeri. Based on a Metaphlan analysis of the skin metagenome reads, the dominant genera in the Bd-endemic site (Sob) were Pseudomonas (74.43%), Variovorax (4.72%), Sanguibacter (4.65%) and Stenotrophomonas (3.23%), while the Bd-naïve site (Sapo) was dominated by Pseudomonas (40.46%), Acinetobacter (30.6%), Staphylococcus (6.38%) and Delftia (5.15%). A hierarchical clustering analysis based on the relative abundance of bacterial taxa suggested that the community structure was different between the three samples from Sob and the three samples from Sapo (Figure 1A). To determine whether the taxonomic composition obtained with 16S rRNA gene amplicon sequencing (Rebollar et al., 2016) was congruent with the shotgun metagenome approach presented in this study, we performed a Procrustes analysis (Figure 1B). Our results indicate that both methods gave similar results when trying to define the dominant groups present on amphibian skin bacterial communities (PROTEST nperm = 999 p = 0.0069).
Figure 1. (A) Pie charts of the most abundant genera obtained with metaphlan for each frog microbiome sample. An UPGMA on the left shows grouping of similar samples based on the relative abundance of bacterial taxa. (B) Procrustes analysis comparing the relative abundance of bacterial taxa (Bray-Curtis dissimilarity matrices) of the 16S rRNA gene amplicon data (arrowheads) and the shotgun metagenome data (black circles). A p-value of 0.006 (PROTEST test) indicates that the matrices are more similar than expected by random association.
Presence of Main Functional Classes and Gene Relative Abundance in C. fitzingeri Skin Microbiomes
We annotated between 25 and 47% of the reads from the six samples obtained from C. fitzingeri skin (Table S1) and classified them according to main functional classes from KEGG. Skin microbiomes contained most of the functional classes identified in KEGG and had a similar proportion of classes among samples (Figure 2A). The most abundant gene classes were amino acid metabolism (mean = 13.87%, SD = 0.11), carbohydrate metabolism (mean = 12.59%, SD = 0.35), energy metabolism (mean = 9.89%, SD = 0.27), membrane transport (mean = 8.76%, SD = 0.60) and metabolism of cofactors and vitamins (mean = 8.24%, SD = 0.16). When we analyzed the KO relative abundances across samples, these were clustered by site based on Bray Curtis distances (Figure 2B). We found significant differences (p-value < 0.05 and LDA score > 2) between sites at the KO level (Table S2). Specifically, we identified 36 KOs that discriminated between Sob and Sapo sites. Of the 36 KOs, 27 and 9 were enriched in Sob and Sapo, respectively. In Sob, most of the enriched genes were part of functional classes involved in membrane transport (13), cellular communication (5), biosynthesis of secondary metabolites (2) and antimicrobial drug resistance (1). In Sapo, enriched genes were associated with the cell cycle (1), lipid metabolism (1) and xenobiotic degradation metabolism (2). In addition to differences found between sites, there is also individual variation among samples within each site as shown in Figure 2B.
Figure 2. Read level analysis (A) Class abundances normalized with MUSiCC. (B) Heatmap showing the relative abundance of KOs across samples. Rows are individual KOs and columns are frog samples. Dendogram at the bottom indicates clustering of samples based on Bray-Curtis distances. Colors indicate the relative abundance (proportions) of KOs across samples (see color legend on the right hand side of the figure).
Genes Involved in Bacterial Communication, Transport and Defense: Unique and Shared Functional Traits Between Sites
Based on assembled contigs, we specifically decided to explore genes from five functional classes associated with bacterial communication, molecular transport and defense mechanisms: biosynthesis of secondary metabolites (BSM), metabolism of terpenoids and polyketides (MTP), membrane transport (MT), cellular communication (CC) and antimicrobial drug resistance (ADR).
Skin metagenomes from both sites had genes for the five functional classes tested, but these genes were associated with different taxonomic groups depending on the site. As seen in Figure 3, most of the genes within each category belonged to Pseudomonas, Delftia, Azotobacter, and Acinetobacter in Sapo and Pseudomonas, Variovorax, Sanguibacter, Stenotrophomonas, and Microbacterium in Sob. However, no significant differences in the number of genes were found between Sapo and Sob sites for these functional categories (LefSe p-values > 0.05).
Figure 3. Contig level analysis. Gene abundance stacked graphs of five functional categories from KEGG that are involved in bacterial communication and bacterial-host-pathogen interactions. Each category represents the gene abundance for Sapo (N = 3) and Sob (N = 3) sites.
We further analyzed the gene abundance for all pathways within each functional category (55 in total). In the case of the BSM class, 16 out of 25 pathways described in KEGG were present in the skin metagenomes. The most abundant were the biosynthesis of monobactam, prodigiosin, streptomycin and novobiocin (Figure S1). In the case of the MTP class, 16 out of 21 pathways described in KEGG were present in the skin metagenomes. The most abundant were terpenoid backbone biosynthesis, geraniol degradation, polyketide sugar unit biosynthesis and biosynthesis of siderophore group nonribosomal peptides (Figure S2). In the case of classes MT, CC and ADR all pathways were present (Figure S3). In the case of the MT class, and as part of the bacteria secretion system pathway, we found almost the complete Type II and Type VI secretion systems, which are involved in secreting molecules (including toxins) to the external environment (Green and Mecsas, 2016). Also, we identified several complete ABC transporters from all the different types of prokaryotic transporters according to KEGG. Within the CC pathways, several genes involved in biofilm formation were identified, as well as quorum sensing genes that had previously been found in Gammaproteobacteria and Bacilli classes. In the case of ADR pathways, we identified most of the genes involved in beta-lactam and vancomycin resistance. In addition, several genes involved in antimicrobial peptide (AMP) resistance were found including the two-component system PhoQ-PhoP which has been identified as an important component in pathogenic and symbiotic bacteria to adapt to host environments (Clayton et al., 2017).
Of all pathways from the five classes, seven were significantly different between sites (LefSe analysis: LDA score > 2 and p-value < 0.05) (Figure 4). We found a significant enrichment of phenylpropanoid biosynthesis genes (BSM class) in Sob in comparison to Sapo. Genes from this biosynthetic pathway were enriched in most of the abundant genera in Sob (Cellulomonas, Cellulosimicrobium, Sanguibacter, and Microbacterium), which were not as abundant in Sapo, and all are from the Actinobacteria class. We found a significant enrichment of prodigiosin biosynthesis genes in Sapo. Genes from this biosynthetic pathway were enriched in most of the abundant genera in Sapo (Pseudomonas, Agromyces, and Acinetobacter), which are all from the Gammaproteobacteria class. In addition, genes from biofilm formation, AMP resistance, bacterial secretion systems and carotenoid biosynthesis were significantly enriched in Sapo (Figure 4).
Figure 4. Contig level analysis. Gene abundance plots of significant pathways enriched in Sob or Sapo sites within five functional classes from KEGG involved in bacterial communication and bacterial-host-pathogen interactions. Each bar represents the gene abundances for Sapo (N = 3) and Sob (N = 3) sites. Colors indicate the taxonomic assignment.
We obtained gene co-occurrence networks for each site based on the relative abundance of genes from the same five functional classes (BSM, MTP, MT, CC, and ADR) (Figure S4). Both networks had a similar number of nodes (KOs): 298 nodes in Sapo and 280 nodes in Sob (Table S3). Interestingly, the number of significant Spearman correlations in the Sob network was twofold greater than in the Sapo network (Figure S4) and included an equivalent amount of negative (mutual exclusions) and positive ones (co-occurrences). Accordingly, the clustering coefficient and network centralization were higher in Sob (Table S3). Moreover, the 229 nodes (KOs) shared between the Sob and Sapo networks had strikingly distinct degree values (number of connections per node) (Figure S5). Thus, we did not find a significant correlation between the degree values of the genes shared in the two networks (τ = 0.0239, p-value = 0.5987), indicating that connections between nodes are not maintained across sites.
Anti-Bd Secondary Metabolite Pathways Are Present in C. fitzingeri Skin Metagenomes
We looked for genes involved in the production of metabolites that have anti-Bd properties such as violacein (Brucker et al., 2008b), indole-3-carboxaldehyde (I3C) (Brucker et al., 2008b), 2,4-diacetylphloroglucinol (2,4-DAPG) (Brucker et al., 2008a), indole-3-ethanol (tryptophol) (Loudon et al., 2014b) and prodigiosin (Woodhams et al., 2017). We identified 16 genes from pathways involved in the production of anti-Bd metabolites (Table 2).
Table 2. Genes involved in the production of anti-Bd metabolites that are present on frog skin metagenomes.
Specifically, we identified three genes from the prodigiosin biosynthesis pathway (BSM) that were present in both sites and significantly enriched in Sapo (Figure 4). One or more genes from this pathway were present in all the abundant bacterial taxa from both sites except for Erwinia, which was only abundant in Sapo. We also identified the coding gene of the enzyme responsible for producing 2,4 DAPG from the terpenoid backbone biosynthesis pathway (MTP) in both sites (K01641). This gene was found only in Staphyloccocus (abundant in Sapo) and Microbacterium (abundant in Sob). In addition, we identified 12 genes that are part of the tryptophan metabolism from which I3C and tryptophol are produced. Specifically, we detected the enzyme that produces indole-3-acetic acid (K00128), from which I3C can be produced (Stutz, 1958). Most of these genes, including K00128, were present in all the abundant bacterial taxa from both sites with a few exceptions (see Table S4).
Discussion
In this study, we described the taxonomic and functional diversity of C. fitzingeri skin metagenomes, and we compared these microbiomes between two populations, one being from a Bd-endemic site (Sob) and the other from a Bd-naïve site (Sapo). To our knowledge, this is the first study describing the functional diversity of skin microbiomes on adult amphibians using shotgun metagenomics. First, we discuss genes and pathways found in both populations that may be required for bacterial interactions and bacterial-host-pathogen interactions. Then, we discuss differences between populations in genes and pathways. We used two approaches (read-level and contig-level) that showed consistent results, although contig-level analyses allowed us to better describe taxonomic and functional features present on C. fitzingeri skin metagenomes.
We focused on the analysis of five main functional classes that could inform us on the capacity of these bacterial communities to communicate with each other, to interact with the host, and to interact with pathogens like Bd: membrane transport (MT), cell communication (CC), biosynthesis of secondary metabolites (BSM), metabolism of terpenoid and polyketides (MTP) and antimicrobial drug resistance class (ADR). Many of the pathways present in these five functional classes are considered essential for bacteria to communicate with each other and to respond to stimuli from the environment. Examples of these pathways are ABC transporters (MT) (Davidson et al., 2008), quorum sensing (CC) (Rutherford and Bassler, 2012), biofilm formation (CC) (Moons et al., 2009), and secondary metabolite production (BSM and MTP) (Flórez et al., 2015). Here, we found that these essential pathways were prominent in C. fitzingeri skin metagenomes independent of bacterial community structure. Based on these results, it appears that many different bacterial taxa can provide the same functions within these communities, suggesting that functional redundancy may be one of the properties present in these symbiotic communities (Foster et al., 2017). Thus, differences in community structure that have been identified among populations and sites in previous studies, may not always mean these communities differ in function (Belden et al., 2015).
In the case of symbiotic bacteria in other hosts, some genes and pathways within these functional classes are known to be important for bacteria to adapt to the host environment, such as the mammalian gut (Frese et al., 2013). One example is the bacterial secretion system pathways (MT), which were initially described in pathogenic bacteria, but have been recently described in symbiotic bacteria (Green and Mecsas, 2016). The presence of almost the complete set of genes for secretion systems II and VI in C. fitzingeri metagenomes suggest that skin bacteria export molecules (perhaps toxins) through these mechanisms and may in turn exert an effect on their host. In addition to the secretion system, other pathways can be important for bacteria to survive on the host; such is the case of the antimicrobial peptide (AMP) production pathways within the ADR class. In C. fitzingeri metagenomes, we found genes from the AMP production pathway. Symbiotic skin bacteria would be expected to evolve mechanisms to avoid the effects of AMPs to colonize the skin. Specifically, the presence of genes within the AMP resistance pathway suggests that skin microbiomes can adapt to the host environment in part through this mechanism. Many amphibian species produce a vast array of AMPs (Conlon, 2011) but so far, no information on these defense molecules has been published for C. fitzingeri.
Skin metagenomes of C. fitzingeri also have a broad range of secondary metabolite pathways (BSM, MTP and ADR classes). These pathways include the biosynthesis of antibiotics, toxins and aromatic compounds, as well as antibiotic resistance pathways, which may be involved in bacterial competition within the community but could also play a role in protecting the host against pathogens. Indeed, it has been proposed that defense of the host arises as a byproduct of microbial competition (Scheuring and Yu, 2012). One example is the microbial symbionts from sponges and corals which produce terpenoids and polyketides (MTP) and play a protective role (Flórez et al., 2015). In C. fitzingeri metagenomes, we identified genes from the MTP class involved in terpenoid and polyketide biosynthetic pathways, antibiotic biosynthesis pathways like ansamycins and vancomycin from the MTP class and antibiotic resistance pathways within the ADR class.
The MTP class also includes pathways involved in the degradation of several compounds produced mainly by plants such as geraniol, pinene and limonene. These pathways have been previously described in other bacterial genera like Pseudomonas and Rhodococcus (Marmulla and Harder, 2014), which can use these toxic compounds as carbon sources. We found these pathways present in C. fitzingeri metagenomes, which is expected considering that frogs are exposed to many of the same or similar carbon sources in their leaf-litter habitat. Also, terpenoid degradation pathways have been found in other bacterial symbionts (Marmulla and Harder, 2014).
We expected to find genes for the production of anti-Bd metabolites, since C. fitzingeri is apparently resistant or tolerant to Bd infection in the wild. We identified some of the genes involved in the production of prodigiosin, I3C, tryptophol and 2,4 DAPG in both sites (Sapo and Sob). However, proving whether these metabolites are produced on C. fitzingeri skin communities would require additional studies that test the presence of their respective transcripts or the metabolites themselves in vivo.
In this study, we identified clear differences in bacterial community structure between sites that were consistent with previous studies using 16S rRNA gene amplicon sequencing (Rebollar et al., 2016). Some members of these skin communities were shared between sites, but some genera were clearly enriched in either Sapo or Sob. Based on the results obtained here, we consider that the genetic/functional differences between sites may be explained by the unique bacterial genera found in each of the populations that we studied. We expect bacterial communities in both sites to provide a defensive function since potential pathogens are ubiquitous in nature; however, we also expected increased selection for a defensive function in Sob since this is considered a Bd-endemic site (Rebollar et al., 2014).
One clear example of the differences between sites in terms of function and taxonomic composition is the enrichment of genes from the phenylpropanoid biosynthesis pathway at the Sob site (Figure 4). Phenylpropanoids are a diverse family of metabolites that are common plant natural products and play several roles, including resistance to pests in plants (reviewed in Vogt, 2010). In this study, most of the bacteria harboring genes from this biosynthetic pathway in Sob were Sanguibacter, Microbacterium, Cellulomonas, and Cellulosimicrobium, which are all Actinobacteria. The Phenylpropanoid biosynthesis pathway has been previously found in Actinobacteria (Moore et al., 2002), and several bioactive molecules from this family have been identified using metabolomics (Wu et al., 2016). In addition, members of the Actinobacteria class also produce a wide range of antibiotics, and have been identified as symbionts that play a crucial role for the protection of several animal hosts (mainly insects) against pathogens (Flórez et al., 2015). Considering that Sob is a Bd-endemic site, we hypothesize that the enrichment of Actinobacteria in Sob provides C. fitzingeri with protection against the pathogen Bd. Further analyses would be necessary to determine if the presence of Actinobacteria in C. fitzingeri skin indeed plays a role in protecting the host against Bd and which molecular mechanisms are involved. We also identified six metabolic pathways that were enriched in Sapo from the CC, ADR, MT and BSM functional classes (Figure 4) that were mainly explained by genes associated to Acinetobacter, Pseudomonas, and Delftia.
The functional traits enriched in either Sapo or Sob likely reflect distinct interactions within the members of the skin microbiome and potentially different ways to interact with the host. Moreover, the differences found in co-occurrence networks between sites (mainly on the number of connection that nodes have) may be caused by the presence of distinct taxonomic groups that harbor unique genetic repertoires. In particular, the functional repertoire of Actinobacteria in Sob and Gammaproteobacteria in Sapo may cause distinct degree values of the KOs within each network.
We suggest that shotgun metagenomics is a promising tool that could allow a deeper understanding of the functions present in amphibian skin microbiomes. This approach could be used not only in the field but also in experimental settings since it could unveil functional changes in time-series or bacterial manipulation experiments (Davis et al., 2017). In the case of the C. fitzingeri skin microbiome, we have been able to describe important pathways involved in bacterial communication, as well as genes involved in potential bacterial-host-pathogen interactions. However, an important caveat of this study is the sample size (3 microbiomes per site) which may influence our results and may not allow us to detect other significant differences between sites. Thus, in the future, we strongly suggest increasing the sample size to fully describe the functional diversity present in amphibian skin microbiomes.
Ethics Statement
Scientific collection permits were provided by the Panamanian authorities (Autoridad Nacional del Ambiente): permits SE/A-47-12, SEX/A-65- 12, SEX/A-77-12 and SEX/A-89-12. Animal care protocols were approved by the Smithsonian Tropical Research Institute's Animal Care Committee: protocol 2011-1110- 2014 and by Virginia Tech's Animal Care Committee: protocol 11-105- BIOL.
Author Contributions
ER, RH, and LB: contributed to the original idea; ER, RH, LB, RJ, JW, and MH: contributed to the design of the research; MH and DM: carried out the fieldwork; ER, JW, MH, and DM: contributed with laboratory analyses; ER, RJ, AG-P, CN, and AE: performed all data analyses; ER, AG-P, CN, AE, EB, and RH: contributed to analyses interpretation; ER: wrote the manuscript and all authors provided critical feedback.
Funding
This project was funded by the NSF Dimensions in Biodiversity Program: DEB-1136602 to RH and DEB-1136640 to LB.
Conflict of Interest Statement
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.
Acknowledgments
We thank Dr. Santiago Ramírez-Barahona for his input and suggestions, and his help on developing the figures. We thank Dr. Jonathan Friedman for his advice on co-occurrence network construction. We thank Dr. Purificación López-García who allowed us to use her computer resources for some of the bioinformatic analyses.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2018.00466/full#supplementary-material
References
Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2
Becker, M. H., Harris, R. N., Minbiole, K. P. C., Schwantes, C. R., Rollins-Smith, L. A., Reinert, L. K., et al. (2011). Towards a better understanding of the use of probiotics for preventing chytridiomycosis in Panamanian golden frogs. Ecohealth 8, 501–506. doi: 10.1007/s10393-012-0743-0
Becker, M. H., Walke, J. B., Cikanek, S., Savage, A. E., Mattheus, N., Santiago, C. N., et al. (2015a). Composition of symbiotic bacteria predicts survival in Panamanian golden frogs infected with a lethal fungus. Proc. R. Soc. B 282:20142881. doi: 10.1098/rspb.2014.2881
Becker, M. H., Walke, J. B., Murrill, L., Woodhams, D. C., Reinert, L. K., Rollins-Smith, L. A., et al. (2015b). Phylogenetic distribution of symbiotic bacteria from Panamanian amphibians that inhibit growth of the lethal fungal pathogen Batrachochytrium dendrobatidis. Mol. Ecol. 24, 1628–1641. doi: 10.1111/mec.13135
Belden, L. K., Hughey, M. C., Rebollar, E. A., Umile, T. P., Loftus, S. C., Burzynski, E. A., et al. (2015). Panamanian frog species host unique skin bacterial communities. Front. Microbiol. 6:1171. doi: 10.3389/fmicb.2015.01171
Bell, S. C., Alford, R. A., Garland, S., Padilla, G., and Thomas, A. D. (2013). Screening bacterial metabolites for inhibitory effects against Batrachochytrium dendrobatidis using a spectrophotometric assay. Dis. Aquat. Organ. 103, 77–85. doi: 10.3354/dao02560
Bletz, M. C., Archer, H., Harris, R. N., McKenzie, V. J., Rabemananjara, F. C. E., Rakotoarison, A., et al. (2017a). Host ecology rather than host phylogeny drives amphibian skin microbial community structure in the biodiversity hotspot of Madagascar. Front. Microbiol. 8:1530. doi: 10.3389/fmicb.2017.01530
Bletz, M. C., Perl, R. G. B., Bobowski, B. T., Japke, L. M., Tebbe, C. C., Dohrmann, A. B., Bhuju, S., et al. (2017b). Amphibian skin microbiota exhibits temporal variation in community structure but stability of predicted Bd-inhibitory function. ISME J. 11, 1521–1534. doi: 10.1038/ismej.2017.41
Boisvert, S., Raymond, F., Godzaridis, E., Laviolette, F., and Corbeil, J. (2012). Ray Meta: scalable de novo metagenome assembly and profiling. Genome Biol. 13:R122. doi: 10.1186/gb-2012-13-12-r122
Brucker, R. M., Baylor, C. M., Walters, R. L., Lauer, A., Harris, R. N., and Minbiole, K. P. C. (2008a). The identification of 2,4-diacetylphloroglucinol as an antifungal metabolite produced by cutaneous bacteria of the salamander plethodon cinereus. J. Chem. Ecol. 34, 39–43. doi: 10.1007/s10886-007-9352-8
Brucker, R. M., Harris, R. N., Schwantes, C. R., Gallaher, T. N., Flaherty, D. C., Lam, B., et al. (2008b). Amphibian chemical defense: antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus. J. Chem. Ecol. 34, 1422–1429. doi: 10.1007/s10886-008-9555-7
Clayton, A. L., Enomoto, S., Su, Y., and Dale, C. (2017). The regulation of antimicrobial peptide resistance in the transition to insect symbiosis. Mol. Microbiol. 103, 958–972. doi: 10.1111/mmi.13598
Conlon, J. M. (2011). The contribution of skin antimicrobial peptides to the system of innate immunity in anurans. Cell Tissue Res. 343, 201–212. doi: 10.1007/s00441-010-1014-4
Costa, S., Lopes, I., Proença, D. N., Ribeiro, R., and Morais, P. V. (2016). Diversity of cutaneous microbiome of Pelophylax perezi populations inhabiting different environments. Sci. Total Environ. 572, 995–1004. doi: 10.1016/j.scitotenv.2016.07.230
Crawford, A. J., Lips, K. R., and Bermingham, E. (2010). Epidemic disease decimates amphibian abundance, species diversity, and evolutionary history in the highlands of central Panama. Proc. Natl. Acad. Sci. U.S.A. 107, 13777–13782. doi: 10.1073/pnas.0914115107
Dale, C., and Moran, N. A. (2006). Molecular interactions between bacterial symbionts and their hosts. Cell 126, 453–465. doi: 10.1016/j.cell.2006.07.014
Davidson, A. L., Dassa, E., Orelle, C., and Chen, J. (2008). Structure, function, and evolution of bacterial ATP-binding cassette systems. Microbiol. Mol. Biol. Rev. 72, 317–364. doi: 10.1128/MMBR.00031-07
Davis, L. R., Bigler, L., and Woodhams, D. C. (2017). Developmental trajectories of amphibian microbiota: response to bacterial therapy depends on initial community structure. Environ. Microbiol. 19, 1502–1517. doi: 10.1111/1462-2920.13707
Fisher, M. C., Garner, T. W. J., and Walker, S. F. (2009). Global emergence of Batrachochytrium dendrobatidis and amphibian chytridiomycosis in space, time, and host. Annu. Rev. Microbiol. 63, 291–310. doi: 10.1146/annurev.micro.091208.073435
Flechas, S. V., Sarmiento, C., Cárdenas, M. E., Medina, E. M., Restrepo, S., and Amézquita, A. (2012). Surviving Chytridiomycosis: differential anti-Batrachochytrium dendrobatidis activity in bacterial isolates from three lowland species of Atelopus. PLoS ONE 7:e44832. doi: 10.1371/journal.pone.0044832
Flórez, L. V., Biedermann, P. H. W., Engl, T., and Kaltenpoth, M. (2015). Defensive symbioses of animals with prokaryotic and eukaryotic microorganisms. Nat. Prod. Rep. 32, 904–936. doi: 10.1039/C5NP00010F
Foster, K. R., Schluter, J., Coyte, K. Z., and Rakoff-nahoum, S. (2017). The evolution of the host microbiome as an ecosystem on a leash. Nat. Publ. Gr. 548, 43–51. doi: 10.1038/nature23292
Franzenburg, S., Walter, J., Künzel, S., Wang, J., Baines, J. F., Bosch, T. C. G., et al. (2013). Distinct antimicrobial peptide expression determines host species-specific bacterial associations. Proc. Natl. Acad. Sci. U.S.A. 110, E3730–E3738. doi: 10.1073/pnas.1304960110
Fraune, S., Anton-Erxleben, F., Augustin, R., Franzenburg, S., Knop, M., Schröder, K., et al. (2014). Bacteria–bacteria interactions within the microbiota of the ancestral metazoan Hydra contribute to fungal resistance. ISME J. 9, 1543–1556. doi: 10.1038/ismej.2014.239
Frese, S. A., MacKenzie, D. A., Peterson, D. A., Schmaltz, R., Fangman, T., Zhou, Y., et al. (2013). Molecular characterization of host-specific biofilm formation in a vertebrate gut symbiont. PLoS Genet. 9:e1004057. doi: 10.1371/journal.pgen.1004057
Gallo, R. L., and Hooper, L. V (2012). Epithelial antimicrobial defence of the skin and intestine. Nat. Rev. Immunol. 12, 503–516. doi: 10.1038/nri3228
Green, E. R., and Mecsas, J. (2016). Bacterial secretion systems – an overview. Am. Soc. Microbiol. 4, 1–32. doi: 10.1128/microbiolspec.VMBF-0012-2015
Grice, E. A. (2014). The skin microbiome: potential for novel diagnostic and therapeutic approaches to cutaneous disease. Semin. Cutan. Med. Surg. 33, 98–103. doi: 10.12788/j.sder.0087
Harris, R. N., Brucker, R. M., Walke, J. B., Becker, M. H., Schwantes, C. R., Flaherty, D. C., et al. (2009). Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824. doi: 10.1038/ismej.2009.27
Harris, R. N., James, T. Y., Lauer, A., Simon, M. A., and Patel, A. (2006). Amphibian pathogen Batrachochytrium dendrobatidis is inhibited by the cutaneous bacteria of amphibian species. Ecohealth 3, 53–56. doi: 10.1007/s10393-005-0009-1
Jackson, D. A. (1995). Protest - a procrustean randomization test of community environment concordance. Ecoscience 2, 297–303. doi: 10.1080/11956860.1995.11682297
Jani, A. J., and Briggs, C. J. (2014). The pathogen Batrachochytrium dendrobatidis disturbs the frog skin microbiome during a natural epidemic and experimental infection. Proc. Natl. Acad. Sci. U.S.A. 111, E5049–E5058. doi: 10.1073/pnas.1412752111
Jorth, P., Turner, K. H., and Gumus, P. (2014). Metatranscriptomics of the human oral microbiome during Health and Disease. MBio. 5, 1–10. doi: 10.1128/mBio.01012-14
Kanehisa, M., and Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30.
Kanehisa, M., Sato, Y., and Morishima, K. (2016). BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. J. Mol. Biol. 428, 726–731. doi: 10.1016/j.jmb.2015.11.006
Kilpatrick, A. M., Briggs, C. J., and Daszak, P. (2010). The ecology and impact of chytridiomycosis: an emerging disease of amphibians. Trends Ecol. Evol. 25, 109–118. doi: 10.1016/j.tree.2009.07.011
Kueneman, J. G., Parfrey, L. W., Woodhams, D. C., Archer, H. M., Knight, R., and McKenzie, V. J. (2014). The amphibian skin-associated microbiome across species, space and life history stages. Mol. Ecol. 23, 1238–1250. doi: 10.1111/mec.12510
Kueneman, J. G., Woodhams, D. C., Harris, R., Archer, H. M., Knight, R., and Mckenzie, V. J. (2016). Probiotic treatment restores protection against lethal fungal infection lost during amphibian captivity. Proc. R. Soc. B Biol. Sci. 283:20161553. doi: 10.1098/rspb.2016.1553
Kueneman, J. G., Woodhams, D. C., Van Treuren, W., Archer, H. M., Knight, R., and McKenzie, V. J. (2015). Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas). ISME J. 10, 934–944. doi: 10.1038/ismej.2015.168.
Li, D., Luo, R., Liu, C. M., Leung, C. M., Ting, H. F., Sadakane, K., et al. (2016). MEGAHIT v1.0: a fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods 102, 3–11. doi: 10.1016/j.ymeth.2016.02.020
Longo, A. V., and Zamudio, K. R. (2017). Environmental fluctuations and host skin bacteria shift survival advantage between frogs and their fungal pathogen. ISME J. 11, 349–361. doi: 10.1038/ismej.2016.138
Loudon, A. H., Holland, J. A., Umile, T. P., Burzynski, E. A., Minbiole, K. P., and Harris, R. N. (2014b). Interactions between amphibians' symbiotic bacteria cause the production of emergent anti-fungal metabolites. Front. Microbiol. 5:441. doi: 10.3389/fmicb.2014.00441
Loudon, A. H., Woodhams, D. C., Parfrey, L. W., Archer, H., Knight, R., McKenzie, V., et al. (2014a). Microbial community dynamics and effect of environmental microbial reservoirs on red-backed salamanders (Plethodon cinereus). ISME J. 8, 830–840. doi: 10.1038/ismej.2013.200
Luo, G., Fotidis, I. A., and Angelidaki, I. (2016). Comparative analysis of taxonomic, functional, and metabolic patterns of microbiomes from 14 full-scale biogas reactors by metagenomic sequencing and radioisotopic analysis. Biotechnol. Biofuels 9:51. doi: 10.1186/s13068-016-0465-6
Manor, O., and Borenstein, E. (2015). MUSiCC: a marker genes based framework for metagenomic normalization and accurate profiling of gene abundances in the microbiome. Genome Biol. 16:53. doi: 10.1186/s13059-015-0610-8
Manor, O., Levy, R., Pope, C. E., Hayden, H. S., Brittnacher, M. J., Carr, R., et al. (2016). Metagenomic evidence for taxonomic dysbiosis and functional imbalance in the gastrointestinal tracts of children with cystic fibrosis. Sci. Rep. 6:22493. doi: 10.1038/srep22493
Marmulla, R., and Harder, J. (2014). Microbial monoterpene transformations-a review. Front. Microbiol. 5:346. doi: 10.3389/fmicb.2014.00346
McHardy, I. H., Goudarzi, M., Tong, M., Ruegger, P. M., Schwager, E., Weger, J. R., et al. (2013). Integrative analysis of the microbiome and metabolome of the human intestinal mucosal surface reveals exquisite inter-relationships. Microbiome 1:17. doi: 10.1186/2049-2618-1-17
McKenzie, V. J., Bowers, R. M., Fierer, N., Knight, R., and Lauber, C. L. (2012). Co-habiting amphibian species harbor unique skin bacterial communities in wild populations. ISME J. 6, 588–596. doi: 10.1038/ismej.2011.129
Medina, D., Walke, J. B., Gajewski, Z., Becker, M. H., Swartwout, M. C., and Belden, L. K. (2017). Culture media and individual hosts affect the recovery of culturable bacterial diversity from Amphibian skin. Front. Microbiol. 8:1574. doi: 10.3389/fmicb.2017.01574
Medina, M., and Sachs, J. L. (2010). Symbiont genomics, our new tangled bank. Genomics 95, 129–137. doi: 10.1016/j.ygeno.2009.12.004
Moons, P., Michiels, C. W., and Aertsen, A. (2009). Bacterial interactions in biofilms bacterial interactions in biofilms. Crit. Rev. Microbiol. 35, 157–168. doi: 10.1080/10408410902809431
Moore, B. S., Hertweck, C., Hopke, J. N., Izumikawa, M., Kalaitzis, J. A., Nilsen, G., et al. (2002). Plant-like biosynthetic pathways in bacteria: from benzoic acid to chalcone. J. Nat. Prod. 65, 1956–1962. doi: 10.1021/np020230m
Oksanen, J., Guillaume, F. B., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2016). vegan: Community Ecology Package. R package Version 2.4-0. Available online at: https://CRAN.R-project.org/package=vegan
Preston, G. M. (2007). Metropolitan microbes: type III secretion in multihost symbionts. Cell Host Microbe 2, 291–294. doi: 10.1016/j.chom.2007.10.004
R Core Team (2015). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available online at: http://www.R-project.org/
Rebollar, E. A., Hughey, M. C., Harris, R. N., Domangue, R. J., Medina, D., Ibáñez, R., et al. (2014). The lethal fungus Batrachochytrium dendrobatidis is present in lowland tropical forests of Far Eastern Panamá. PLoS ONE 9:e95484. doi: 10.1371/journal.pone.0095484
Rebollar, E. A., Hughey, M. C., Medina, D., Harris, R. N., Ibáñez, R., and Belden, L. K. (2016). Skin bacterial diversity of Panamanian frogs is associated with host susceptibility and presence of Batrachochytrium dendrobatidis. 10, 1682–1695. doi: 10.1038/ismej.2015.234
Rodríguez-Brenes, S., Rodriguez, D., Ibañez, R., and Ryan, M. J. (2016). Spread of amphibian chytrid fungus across lowland populations of Túngara frogs in panamá. PLoS ONE 11:e0155745. doi: 10.1371/journal.pone.0155745
Rutherford, S. T., and Bassler, B. L. (2012). Bacterial quorum sensing: its role in virulence and possibilities for its control. Cold Spring Harb. Perspect. Med. 2, 1–25. doi: 10.1101/cshperspect.a012427
Scheuring, I., and Yu, D. W. (2012). How to assemble a beneficial microbiome in three easy steps. Ecol. Lett. 15, 1300–1307. doi: 10.1111/j.1461-0248.2012.01853.x
Schmid, J., Sieber, V., and Rehm, B. (2015). Bacterial exopolysaccharides: biosynthesis pathways and engineering strategies. Front. Microbiol. 6:496. doi: 10.3389/fmicb.2015.00496
Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069. doi: 10.1093/bioinformatics/btu153
Segata, N., Izard, J., Waldron, L., Gevers, D., Miropolsky, L., Garrett, W. S., et al. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12:R60. doi: 10.1186/gb-2011-12-6-r60
Segata, N., Waldron, L., Ballarini, A., Narasimhan, V., Jousson, O., and Huttenhower, C. (2012). Metagenomic microbial community profiling using unique clade-specific marker genes. Nat. Methods 9, 811–814. doi: 10.1038/nmeth.2066
Stutz, R. E. (1958). Enzymatic formation of indole-3-carboxaldehyde from indole-3-acetic acid. Plant Physiol. 33, 207–212. doi: 10.1104/pp.33.3.207
Turnbaugh, P. J., Ley, R. E., Hamady, M., Fraser-liggett, C., Knight, R., and Gordon, J. I. (2007). The human microbiome project: exploring the microbial part of ourselves in a changing world. Nature 449, 804–810. doi: 10.1038/nature06244
Walke, J. B., Becker, M. H., Loftus, S. C., House, L. L., Cormier, G., Jensen, R. V., et al. (2014). Amphibian skin may select for rare environmental microbes. ISME J. 8, 1–11. doi: 10.1038/ismej.2014.77
Walke, J. B., Becker, M. H., Loftus, S. C., House, L. L., Teotonio, T. L., Minbiole, K. P. C., et al. (2015). Community structure and function of amphibian skin microbes: an experiment with Bullfrogs exposed to a Chytrid fungus. PLoS ONE 10:e0139848. doi: 10.1371/journal.pone.0139848
Walke, J. B., and Belden, L. K. (2016). Harnessing the microbiome to prevent fungal infections: lessons from amphibians. PLoS Pathog. 12, 6–11. doi: 10.1371/journal.ppat.1005796
Woodhams, D. C., Alford, R. A., Antwis, R. E., Archer, H., Becker, M. H., Belden, L. K., et al. (2015). Antifungal isolates database of amphibian skin-associated bacteria and function against emerging fungal pathogens. Ecology 96, 595–595. doi: 10.1890/14-1837.1
Woodhams, D. C., Labumbard, B. C., Barnhart, K. L., Becker, M. H., Bletz, M. C., Escobar, L. A., et al. (2017). Prodigiosin, violacein, and volatile organic compounds produced by widespread cutaneous bacteria of amphibians can inhibit two Batrachochytrium fungal pathogens. Microb. Ecol. doi: 10.1007/s00248-017-1095-7. [Epub ahead of print].
Keywords: skin microbiome, shotgun metagenomics, host-bacteria interactions, amphibians, Batrachochytrium dendrobatidis
Citation: Rebollar EA, Gutiérrez-Preciado A, Noecker C, Eng A, Hughey MC, Medina D, Walke JB, Borenstein E, Jensen RV, Belden LK and Harris RN (2018) The Skin Microbiome of the Neotropical Frog Craugastor fitzingeri: Inferring Potential Bacterial-Host-Pathogen Interactions From Metagenomic Data. Front. Microbiol. 9:466. doi: 10.3389/fmicb.2018.00466
Received: 02 November 2017; Accepted: 28 February 2018;
Published: 20 March 2018.
Edited by:
Sebastian Fraune, Christian-Albrechts-Universität zu Kiel, GermanyReviewed by:
Diogo Neves Proença, University of Coimbra, PortugalLucia Pita, GEOMAR Helmholtz Centre for Ocean Research Kiel, Germany
Copyright © 2018 Rebollar, Gutiérrez-Preciado, Noecker, Eng, Hughey, Medina, Walke, Borenstein, Jensen, Belden and Harris. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Eria A. Rebollar, ea.rebollar@gmail.com
†Present Address: Eria A. Rebollar, Departamento de Ecología Evolutiva, Instituto de Ecología, Universidad Nacional Autónoma de México, Mexico City, Mexico
Jenifer B. Walke, Department of Biology, Eastern Washington University, Cheney, WA, United States