Head-starting of the federally endangered Houston toad (Anaxyrus houstonensis), that is, the release of egg strands, tadpoles, and metamorphic juveniles produced in captivity into the original breeding ponds, requires assessment of potential threats for the transmission of pathogens from captive to free-ranging toads. We used Illumina-based 16S rRNA V3 amplicon sequencing to investigate the community structure of bacteria from skin lesions of captive Houston toad and habitat (pond) samples. Proteobacteria, alone or together with Actinobacteria and, in some samples, Cyanobacteria represented virtually all reads in tissue lesion samples, whereas pond samples were much more diverse, with Acidobacteria, Actinobacteria, Bacteriodetes, Chloroflexi, Cyanobacteria, Firmicutes, Planctomycetes, Proteobacteria, and Verrucomicrobia present with little variation between samples. If present in lesions, Actinobacteria were largely represented by Mycobacteriaceae, and here mainly by one sequence identical to sequences of members of the Mycobacterium chelonaeabscessus complex. In pond samples, mycobacteria represented only a small portion of the actinobacteria, although at higher diversity with six distinct reads. Sequences for reads obtained from pond samples were identical to those representing the M. chelonaeabscessus complex, a group with Mycobacterium marinum, Mycobacterium kansasii, Mycobacterium avium, a group with Mycobacterium vaccae, Mycobacterium fortuitum, Mycobacterium poriferae, and a group with Mycobacterium elephantis and Mycobacterium celeriflavum, whereas sequences of high similarity were detected for reads related to those of Mycobacterium holsaticum, Mycobacterium pallens, and Mycobacterium obuense, and Mycobacterium goodii. Our results indicated that lesions observed on the Houston toad in captivity are not the result of mycobacteria in every case, and that the presence of mycobacteria in the captive colony does not represent a novel pathogen threat to the wild populations because such bacteria are also seen in the natural pond habitats for the Houston toad.

Habitat destruction and fragmentation by urban and agricultural growth, the introduction of alien species, climate change, environmental pollution, and pathogens, and interactions among these factors have been discussed as underlying reasons for steep declines of amphibian populations, with effects on both abundance and diversity (Bucciarelli et al. 2014; Blaustein et al. 2018; Tornabene et al. 2018). Effects are especially pronounced for habitat specialists because even small changes in environmental characteristics might be detrimental to the populations (Klaus and Noss 2016). The Houston toad (Anaxyrus houstonensis) is a unique Texas endemic, and represents a habitat specialist that occurs in deep, sandy soils under canopy forests in only eight counties in central Texas, US, where it persists at low population levels with evidence of continued declines (Duarte et al. 2014). It is listed on the International Union for Conservation of Nature endangered species list, with habitat loss, droughts, and predation by imported fire ants (Solenopsis invicta) mentioned as contributors to declines (Duarte et al. 2014).

Conservation efforts for endangered amphibian populations often include head-starting, or removal of individuals from at-risk areas, captive propagation, and the subsequent release of offspring into the environment (Clulow et al. 2014; Silla and Byrne 2019). For the Houston toad, head-starting and the development of captive assurance colonies were implemented in 2007, with population assessments, habitat recovery, and restoration introduced as additional conservation efforts conducted by a group of dedicated private landowners, zoos, state and federal agencies, and university research groups. The head-starting efforts for the Houston toad include the collection of subsets of egg strands (less than 50% of any given strand) from natal ponds, propagation of offspring in captivity and the subsequent release of high numbers of eggs, tadpoles, or juveniles back into the natal ponds and adjacent habitat (Forstner and Crump 2011). Some Houston toads are retained from the head-starting effort as the nascent portion of captive assurance colonies in participating entities (i.e., the Houston, Dallas, and Fort Worth Zoos, and the US Fish and Wildlife Service's San Marcos National Fish Hatchery and Technology Center).

High-density propagation in captivity creates stressful conditions for animals, generally increasing susceptibility to diseases and parasites, but also promoting their transmission between individuals. Numerous pathogens on amphibians have been identified (Densmore and Green 2007), with Batrachochytrium dendrobatidis, a zoosporic fungus causing chytridiomycosis currently receiving the most attention (Lips et al. 2003; Fisher and Garner 2020). Viruses of the genus Ranavirus and several bacteria, Aeromonas, Chlamydia and Mycobacteria, are known to contribute to infectious diseases of amphibians (Gray et al. 2009; Hill et al. 2010; Blaustein et al. 2018). Nematodes known to contribute to mortality in captive amphibians have also been identified in holding facilities of captive assurance colonies (Bianchi et al. 2014). Chlamydia pneumonia infections causing neurologic disease and mortality were identified in a breeding colony at the Houston Zoo (Fratzke et al. 2019). Mycobacteria have been detected in skin lesions on many amphibians (Ferreira et al. 2006; Sanchez-Morgado et al. 2009), although they are often considered as opportunistic infections supported by detrimental environmental conditions (Pessier 2002). Recently, acid-fast bacteria, presumably mycobacteria, were also detected in lesions on several individuals of the Houston toad in the captive assurance colony at the Houston Zoo.

The goal of this study was to analyze the bacterial community in tissue lesions from several euthanized Houston toads using next generation sequencing, with special emphasis on the most abundant bacterial genera. Bacterial communities in skin lesions were compared to those from composite water-sediment samples from Houston toad breeding ponds located at the Griffith League Ranch (GLR) in Bastrop County, Texas. Some of these ponds had been used for releases of head-started Houston toads; however, most were native breeding ponds that have never received head-started toads (Forstner and Ahlbrandt 2003).

Our work was approved by the Institutional Animal Care and Use Committee (permit no. 7E1EC3_02) at Texas State University. Tissue samples of skin lesions were obtained from 11 adult Houston toads, that had been hatched in captivity and kept at the Houston Zoo for 6–7 yr as part of a captive assurance colony (Fratzke et al. 2019). The Houston Zoo maintains a 120 m2 Houston toad quarantine facility that is used for captive breeding to support head-starting at the original breeding ponds at GLR. Toads exhibiting skin lesions (Fig. 1) were euthanized between June 2019 and January 2020, and skin lesion tissues were examined for acid-fast bacteria using a modified Kinyoun Stain method (Eng Scientific, Inc., Clifton, New Jersey, USA) at the Houston Zoo. Tissue samples of skin lesions with positive detections were stored frozen, and shipped frozen to Texas State University.

Figure 1

Representative skin lesions found in Houston toads (Anaxyrus houstonensis) that have acid-fast stain-positive bacteria present in aspirates stained with a modified Kinyoun Stain method. Tissue samples of skin lesions were obtained from 11 adult Houston toads that had been hatched in captivity and kept at the Houston Zoo (Houston, Texas, USA) for 6–7 yr as part of a captive assurance colony.

Figure 1

Representative skin lesions found in Houston toads (Anaxyrus houstonensis) that have acid-fast stain-positive bacteria present in aspirates stained with a modified Kinyoun Stain method. Tissue samples of skin lesions were obtained from 11 adult Houston toads that had been hatched in captivity and kept at the Houston Zoo (Houston, Texas, USA) for 6–7 yr as part of a captive assurance colony.

Close modal

Composite water and sediment samples (approx. 1/1, vol/vol) were obtained in sterile 50 mL tubes from 12 ponds (identified as 1–16, without 3, 6, 8, and 13) at GLR in early January 2020, and processed the same day. Wild Houston toads had been found in all of the ponds, and only ponds 2, 5, and 12 were used for head-starting efforts. Water and sediment in samples were mixed by shaking, and 1 mL subsamples transferred to 2 mL cryotubes. After centrifugation at 14,000 × G for 5 min, supernatant was discarded, and DNA was extracted from the remaining pellet (about 300 mg wet weight) using the SurePrepTM Soil DNA Isolation Kit (Fisher Scientific, Houston, Texas, USA), with small modifications (Samant et al. 2012). These modifications included omission of the humic acid removal step, and the bead-beating step was done with 50% of the recommended liquid. The remaining 50% of the liquid was used to re-extract the sample after the bead-beating step. Tissue specimens (100 mg) were treated the same way. The DNA concentrations in extracts were measured with a Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, California, USA).

For Illumina sequencing, 16S rRNA gene fragments were amplified from DNA extracts using primer 515f and barcoded primers 806r, both of which included linker sequences following the instructions from the Earth Microbiome Project (Caporaso et al. 2011, 2012; Gilbert et al. 2014). The PCR was carried out in a 100 µL volume with 1 × Taq Buffer, 1.5 mM MgCl2, 0.2 mM deoxyribonucleotide triphosphates, 0.2 µM primers, 2.5 µg/µL bovine serum albumin, 1 U Taq polymerase (GenScript, Inc., Piscataway, New Jersey, USA), and 1 µL DNA extract. The PCR conditions included an initial denaturation at 94 C for 3 min followed by 35 cycles of 94 C for 45 s, 50 C for 60 s, and 72 C for 90 s, with a final 72 C extension for 10 min.

We cleaned the PCR products using the PureLink® PCR Purification Kit (Invitrogen, Waltham, Massachusetts, USA), and then they were checked and quantified on a 2200 TapeStation (Agilent Technologies, Santa Clara, California, USA) using the Agilent DNA D1000 ScreenTape Analysis (Agilent Technologies). Samples were analyzed on the Illumina MiSeq v3 (Illumina, Inc., San Diego, California, USA) with paired end 2×300 base pair reads using the respective sequencing and index sequence primers.

Primer and adapter-linker sequences were first removed from FASTQ files containing the raw sequence data generated through Illumina sequencing. The Dada2R package version 1.8 (Callahan et al. 2016) was used in R studio (R version 3.4.3; R studio version 1.1.456) to further process the demultiplexed paired reads. Briefly, reads were checked for quality using the R function plotQualityProfile(), followed by filtering and trimming forward and reverse paired reads to 280 and 200 base pair, respectively, with a maximum of two expected errors per read using the R function filterAndTrim() function. Reads were then dereplicated using R function derepFastq(). To infer true sequence variants from unique sequences in the samples, a core sample inference algorithm using R function Dada() was applied. Inferred forward and reverse sequences were merged using R function mergePairs(). Sequences that did not perfectly overlap were considered errors and were removed. An amplicon sequence variant table was constructed using the R function makeSequenceTable(). Chimeras were then removed from the table using the R function removeBimeraDenovo(). A new table with no chimeras was subsequently generated. Taxonomy was then assigned using R function assignTaxonomy() to compare reads to Dada2 formatted reference FASTA files from the National Center for Biotechnology Information RefSeq 16S rRNA database supplemented with sequences from the Ribosomal Database Project (Alishum 2019).

Reads with sequences identifying selected genera were assembled in Geneious 11.1.5 (Bio-matters Ltd., Auckland, New Zealand), aligned by the Geneious alignment tool, and checked in GenBank databases (National Center for Biotechnology Information 2020) using the Basic Alignment Search Tool (Pearson and Lipman 1988). Representative sequences from confirmed strains were added from GenBank databases, and all sequences were realigned. The identities and relationships among the sequences amplified were evaluated using neighbor-joining (NJ) analysis (Saitou and Nei 1987) and maximum likelihood (ML) analysis (Stamatakis 2006) from within Geneious 11.1.5. The neighbor-joining analyses utilized the HKY85 model to correct for substitution bias (Hasegawa et al. 1985). Model parameters for maximum likelihood, which were estimated by the general time reversible model with gamma (Tavaré 1986), were used as input in a ML heuristic search using RAxML (Stamatakis 2014). Bootstrap values (Felsenstein 1985) were estimated from a heuristic search with random stepwise addition sequence for 10,000 NJ and 10,000 ML iterations.

Illumina-based 16S rRNA V3 amplicon sequencing resulted in a total of 7.6 million effective reads from 23 samples (e.g., tissue lesions from 11 toads, with 104,223 to 944,227 reads per sample, mean±SD: 310,514±222,552 and 12 ponds, with 141,058 to 1,126,442 reads per sample, mean=349,198±309,849). We differentiated reads from sequences to track the input dataset which was then analytically confirmed to a taxonomic sequence identity. Most sequences from tissue samples were identified in the reads as Proteobacteria (range=56.4–100%; mean=81.1±15.8% of all reads), with sequences representing Actinobacteria as second most abundant phylum being much less frequent (range=0–35.8%; mean=9.8±10.6% of all reads) as shown in Table 1. Proteobacteria, alone or together with Actinobacteria and, in some samples, Cyanobacteria, accounted for the majority of all sequences (Fig. 2a), with reads containing additional phyla such as Acidobacteria, Bacteriodetes, Chloroflexi, Firmicutes, Planctomycetes, and Verrucomicrobia as well as several others being present in small numbers, with high variability in abundance among samples (Fig. 2a). In ponds, reads representing Proteobacteria were less abundant than in tissue lesions, with little variation among samples (range=36.8–50.2%; mean=43.1±3.9% of all reads) as shown in Table 2. Little variation among samples was found for percent sequences identifying other phyla, including Actinobacteria with a mean of 4.9±0.8% of all reads, and Acidobacteria, Bacteriodetes, Chloroflexi, Cyanobacteria, Firmicutes, Planctomycetes, and Verrucomicrobia with similar values that all together resembled about 80–85% of all reads (Fig. 2a). The remaining sequences represented a large variety of phyla at abundances generally below 1% of all reads.

Table 1

Identification of specific members of the phyla Proteobacteria and Actinobacteria in skin lesion samples from Houston toads (Anaxyrus houstonensis) kept at the Houston Zoo (Houston, Texas, USA) as part of a captive assurance colony in percent of all reads obtained by Illumina sequencing (Caporaso et al. 2012) and analyzed using the Dada2 R package version 1.8 (Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

Identification of specific members of the phyla Proteobacteria and Actinobacteria in skin lesion samples from Houston toads (Anaxyrus houstonensis) kept at the Houston Zoo (Houston, Texas, USA) as part of a captive assurance colony in percent of all reads obtained by Illumina sequencing (Caporaso et al. 2012) and analyzed using the Dada2 R package version 1.8 (Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).
Identification of specific members of the phyla Proteobacteria and Actinobacteria in skin lesion samples from Houston toads (Anaxyrus houstonensis) kept at the Houston Zoo (Houston, Texas, USA) as part of a captive assurance colony in percent of all reads obtained by Illumina sequencing (Caporaso et al. 2012) and analyzed using the Dada2 R package version 1.8 (Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).
Figure 2

Reads containing sequences identifying major phyla (a), and major families within the two most abundant phyla, the γ-Proteobacteria (b) and Actinobacteria (c) in skin lesion specimens of Houston toad (Anaxyrus houstonensis) and in water-sediment samples from ponds at Griffith League Ranch (GLR), Bastrop, Texas, USA. Reads were generated from 16S rRNA gene amplicons following the instructions from the Earth Microbiome Project (Caporaso et al. 2012) and analyzed on the Illumina MiSeq v3. Raw reads were further analyzed, and the sequences subsequently identified using the Dada2 R package (version 1.8; Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

Figure 2

Reads containing sequences identifying major phyla (a), and major families within the two most abundant phyla, the γ-Proteobacteria (b) and Actinobacteria (c) in skin lesion specimens of Houston toad (Anaxyrus houstonensis) and in water-sediment samples from ponds at Griffith League Ranch (GLR), Bastrop, Texas, USA. Reads were generated from 16S rRNA gene amplicons following the instructions from the Earth Microbiome Project (Caporaso et al. 2012) and analyzed on the Illumina MiSeq v3. Raw reads were further analyzed, and the sequences subsequently identified using the Dada2 R package (version 1.8; Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

Close modal
Table 2

Identification of specific members of the phyla Proteobacteria and Actinobacteria in composite water-sediment samples from Houston toad (Anaxyrus houstonensis) breeding ponds located at the Griffith League Ranch in Bastrop County, Texas, USA, in percent of all reads obtained by Illumina sequencing (Caporaso et al. 2012) and analyzed using the Dada2 R package version 1.8 (Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

Identification of specific members of the phyla Proteobacteria and Actinobacteria in composite water-sediment samples from Houston toad (Anaxyrus houstonensis) breeding ponds located at the Griffith League Ranch in Bastrop County, Texas, USA, in percent of all reads obtained by Illumina sequencing (Caporaso et al. 2012) and analyzed using the Dada2 R package version 1.8 (Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).
Identification of specific members of the phyla Proteobacteria and Actinobacteria in composite water-sediment samples from Houston toad (Anaxyrus houstonensis) breeding ponds located at the Griffith League Ranch in Bastrop County, Texas, USA, in percent of all reads obtained by Illumina sequencing (Caporaso et al. 2012) and analyzed using the Dada2 R package version 1.8 (Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

In tissue specimens, γ-Proteobacteria represented between 80–100% of the Proteobacteria and 30–99.9% of all reads (Table 1). In pond samples, sequences identifying α-, β-, δ-, and γ-Proteobacteria were present in similar percentages for each of the classes, with those representing γ-Proteobacteria making up between 6.8–9.6% of all sequences except for pond 12 with 30.9% (Table 2). In most tissue specimens γ-Proteobacteria were represented by Enterobacteriaceae and Pseudomonadaceae (about 60–100%), with Xanthomonadaceae present in some samples (Fig. 2b). Diversity within the γ-Proteobacteria was much more pronounced in pond samples, with additional families detected in similar community structure patterns among ponds, except for pond 12 where Enterobacteriaceae and Pseudomonadaceae were most prominent (Fig. 2b). Reads identifying Actinobacteria were not detected in three tissue specimens, but were present in the remaining specimens in numbers accounting for 1.2% to 35.8% of all reads (mean=13.2±10.3%). Populations in remaining samples were generally dominated by individual families; for example, the Mycobacteriaceae in seven samples (mean=13.2±10.5%) as shown in Table 1, and the Nocardiaceae in one sample (Fig. 2c). Reads identifying Actinobacteria were also present in pond samples, although at much lower percentages (3.7–6.7% of all reads; mean=4.9±0.8%) as shown in Table 2. Mycobacteriaceae accounted for a mean of 0.22±0.13% of all reads, and thus represented fewer than 10% of the Actinobacteria in very consistent abundance patterns between samples, except for pond 12 where Mycobacteriaceae made up about 70% of the Actinobacteria (Fig. 2c).

The families Pseudomonadaceae, Enterobacteriaceae, and Mycobacteriaceae were generally represented by one to three sequences per sample with patterns relatively consistent in the groups of samples collected from each source (i.e., pond or tissue), but there were differences in the patterns between the sources (Fig. 3). Pseudomonadaceae in specimen samples, for example, were represented by three major sequences (read-1, read-4, read-42), with one of these dominant in a tissue specimen. In pond samples, these sequences were present in very small numbers, except for in samples from pond 12 (read-42). Pond samples were dominated by read-154, with many additional, yet unidentified reads present (Fig. 3a). Enterobacteriaceae in tissue samples were also generally represented by one major sequence (read-37, read-54, or read-135), with read-135 abundant also in most pond samples (Fig. 3b). Pond samples were otherwise dominated by sequences not detected in significant numbers in tissue specimen samples (e.g., read-480), and additional unidentified reads (Fig. 3b).

Figure 3

Sequences identified within the families Pseudomonadaceae (a), Enterobacteriaceae (b), and Mycobacteriaceae (c) in skin lesion specimens of Houston toad (Anaxyrus houstonensis) and in water-sediment samples from ponds at Griffith League Ranch (GLR), Bastrop, Texas, USA. Reads were generated from 16S rRNA gene amplicons following the instructions from the Earth Microbiome Project (Caporaso et al. 2012) and analyzed on the Illumina MiSeq v3. Raw reads were further analyzed, and the sequences subsequently identified using the Dada2 R package (version 1.8; Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

Figure 3

Sequences identified within the families Pseudomonadaceae (a), Enterobacteriaceae (b), and Mycobacteriaceae (c) in skin lesion specimens of Houston toad (Anaxyrus houstonensis) and in water-sediment samples from ponds at Griffith League Ranch (GLR), Bastrop, Texas, USA. Reads were generated from 16S rRNA gene amplicons following the instructions from the Earth Microbiome Project (Caporaso et al. 2012) and analyzed on the Illumina MiSeq v3. Raw reads were further analyzed, and the sequences subsequently identified using the Dada2 R package (version 1.8; Callahan et al. 2016) and a Dada2 formatted reference sequence database (Alishum 2019).

Close modal

Mycobacteriaceae, when present in tissue samples, were represented in all but one sample by the same sequence (read-101; Fig. 3c). In addition to read-101, one specimen (26763) contained an additional sequence (read-5213) with an abundance of about 1% of all sequences representing Mycobacteriaceae (Fig. 3c). Both reads were found in pond samples as well; however, Mycobacteriaceae were generally found at higher diversity in pond samples than in tissue samples, with three to five different sequences per sample (Fig. 3c). Mycobacteriaceae in only two ponds were represented by single sequences, one by read-101 (pond 12) and another by read-1996 (pond 9). Additional sequences found were read-1428, read-1665, read-5213, and read-13792 for a total of six sequences detected in pond samples (Fig. 3c).

All reads representing Pseudomonadaceae matched identical sequences in the databases, generally identifying Pseudomonas spp., with the individual sequence identical to those representing different species (Fig. 4). The same result was obtained for Enterobacteriaceae, with reads identical to sequences from a large variety of species belonging to different genera, including Providentia, Salmonella, Morganella, Escherichia, Klebsiella, or En-terobacteria (Fig. 4). For Mycobacteriaceae, comparative sequence analyses revealed that read-101 was identical to sequences of 16S rRNA gene fragments of Mycobacterium chelonae, Mycobacterium abscessus, Mycobacterium franklinii, and Mycobacterium saopaulense, whereas read-5213 showed high sequence similarity to the 16S rRNA gene fragments of Mycobacterium pallens and Mycobacterium holsaticum (Fig. 5). Read-1428 was identical to sequences of Mycobacterium avium, Mycobacterium marinum, Mycobacterium kansasii, and others (including Mycobacterium tuberculosis), and read-1996 is identical to Mycobacterium fortuitum, Mycobacterium poriferae, Mycobacterium vaccae, and others. Read-13792 is identical to Mycobacterium celeriflavum, Mycobacterium elephantis, and others, and read-1665 has high sequence similarity to Mycobacterium obuense and Mycobacterium goodie (Fig. 5).

Figure 4

Sequence relationships for reads identified as enterobacteria or pseudomonads. Reads were aligned together with representative sequences from confirmed strains by the Geneious alignment tool in Geneious 11.1.5. The identity and relationship among the sequences were evaluated using neighbor joining (NJ) and maximum likelihood (ML) analyses from within Geneious 11.1.5. Numbers above the branches represent the bootstrap values from a NJ bootstrap analysis (10,000 replicates) using the HKY85 correction (Hasegawa et al. 1985), followed by ML bootstrap (10,000 replicates) values, respectively, noted for clades with greater than 45% bootstrap support. All results were plotted on the NJ bootstrap topology.

Figure 4

Sequence relationships for reads identified as enterobacteria or pseudomonads. Reads were aligned together with representative sequences from confirmed strains by the Geneious alignment tool in Geneious 11.1.5. The identity and relationship among the sequences were evaluated using neighbor joining (NJ) and maximum likelihood (ML) analyses from within Geneious 11.1.5. Numbers above the branches represent the bootstrap values from a NJ bootstrap analysis (10,000 replicates) using the HKY85 correction (Hasegawa et al. 1985), followed by ML bootstrap (10,000 replicates) values, respectively, noted for clades with greater than 45% bootstrap support. All results were plotted on the NJ bootstrap topology.

Close modal
Figure 5

Sequence relationships for reads containing sequences identified as mycobacteria. Reads were aligned together with representative sequences from confirmed strains by the Geneious alignment tool in Geneious 11.1.5. The identity and relationship among the sequences were evaluated using neighbor joining (NJ) and maximum likelihood (ML) analyses from within Geneious 11.1.5. Numbers above the branches represent the bootstrap values from a NJ bootstrap analysis (10,000 replicates) using the HKY85 correction (Hasegawa et al. 1985), followed by ML bootstrap (10,000 replicates) values, respectively, noted for clades with greater than 45% bootstrap support. All results were plotted on the NJ bootstrap topology.

Figure 5

Sequence relationships for reads containing sequences identified as mycobacteria. Reads were aligned together with representative sequences from confirmed strains by the Geneious alignment tool in Geneious 11.1.5. The identity and relationship among the sequences were evaluated using neighbor joining (NJ) and maximum likelihood (ML) analyses from within Geneious 11.1.5. Numbers above the branches represent the bootstrap values from a NJ bootstrap analysis (10,000 replicates) using the HKY85 correction (Hasegawa et al. 1985), followed by ML bootstrap (10,000 replicates) values, respectively, noted for clades with greater than 45% bootstrap support. All results were plotted on the NJ bootstrap topology.

Close modal

High-throughput 16S rRNA gene sequencing of tissue lesion specimens from the Houston toad revealed extremely limited diversity, often with two to three sequences representing virtually the entire bacterial community, mainly consisting of Pseudomonadaceae, Enterobacteriaceae, and Mycobacteriaceae. These reads were identical to a variety of bacteria often found as opportunistic pathogens, such as members of the genera Salmonella, Escherichia, Klebsiella, and Pseudomonas or the Mycobacterium chelonae–abscessus complex (O'Rourke and Rosenbaum 2015). These sequences were also present in pond samples, although at abundances two to three orders of magnitude lower than in tissue lesion specimens, and within bacterial communities of much higher diversity. These results indicated that specific environmental microbes might proliferate in skin and tissue lesions of Houston toads, and act as opportunistic pathogens (Good 1985; Falkinham 2002; Primm et al. 2004).

Amphibian skin samples often show Proteobacteria as the most prominent phylum, with relative abundances accounting for 50–100% of all sequences, followed by much smaller numbers of sequences representing Actinobacteria and Firmicutes, and others, such as Cyanobacteria and Bacteriodetes (McKenzie et al. 2012; Walke et al. 2014; Hughey et al. 2017). Contributions of phyla to relative abundances are amphibian-specific and not a reflection of the environmental microbial community, and can vary widely (McKenzie et al. 2012; Kueneman et al. 2014; Walke et al. 2014). Thus, on the phylum level, our analyses of Houston toad tissue lesions were very similar to those often observed for skin microbiomes in amphibians without lesions. A similar statement could be made for pond samples that showed very similar patterns of abundance, distribution, and stability of bacterial phyla to other next-generation sequencing-based studies of bacterial communities in ponds (Fan et al. 2016; Jeong and Ham 2017; Mania et al. 2019).

High-throughput 16S rRNA gene sequencing analyses of microbial communities generally focus on phylogenetic levels between phyla and families, and occasionally include assessments on the genus level (Li et al. 2017; Mashiane et al. 2017). Analyses on lower taxonomic levels such as species or subspecies levels are usually not attempted, because they are often impacted by the limited phylogenetic resolution of sequences of 16S rRNA gene fragments, and the limited number of reads available at this level, as well as by the lack of sufficient sequences for distinct species or subspecies in the databases (Poretsky et al. 2014). In our tissue specimens, individual sequences identifying Proteobacteria might therefore represent different species or even different genera, covering a range from potentially opportunistic pathogens to beneficial symbionts. Consequently, the assignment of a short sequence to a specific bacterial species within the families Pseudomonadaceae and Enterobacteriaceae and the concomitant identification of a specific harmful or beneficial mode of interaction remains highly speculative.

Actinobacteria in our tissue specimens were typically represented by Mycobacteriaceae and here by one sequence, read-101, identifying members of the M. chelonae–abscessus complex, all with identical sequences to read-101, only (Nogueira et al. 2015). The clustering of these short 16S rRNA gene fragments was consistent with results obtained by whole genome analyses (Tortoli et al. 2017). Although individual species could not be separated in our analysis, members of this complex cause zoonotic diseases such as mycobacteriosis in amphibians (Green et al. 2000; Barrows et al. 2017). Sequences representing members of the M. chelonae–abscessus complex were also found in pond samples, as were additional sequences representing other mycobacteria that have been identified in diseased amphibians, such as M. marinum (Ferreira et al. 2006; Haridy et al. 2014), M. avium, or M. fortuitum (Gcebe et al. 2018). Although these species represent opportunistic pathogens, many mycobacteria are saprophytic and thus autochthonous members of terrestrial and aquatic microbial communities (Hruska and Kaevska 2012; Roguet et al. 2016). Although M. chelonae, M. pallens, and M. fortuitum potentially identified in our study as pond residents were not detected in soils (Walsh et al. 2019), they have been identified in high abundance in pond water and sediments, as have many other Mycobacterium species (Hruska and Kaevska 2012; Roguet et al. 2016). Mycobacterium vaccae and M. holsaticum have also been detected in soils (Pontiroli et al. 2013). Thus, mycobacteria in Houston toad specimens, identified as members of the M. chelonae–abscessus complex, are present in pond samples as well, and therefore seem to represent autochthonous members of aquatic microbial communities that infect toads as opportunistic pathogens (Good 1985; Primm et al. 2004; Bland et al. 2005).

Although microscopic observations identified acid-fast bacteria in all tissue specimen, not all tissue specimens contained mycobacteria. Members of the genus Nocardia that were detected in one of our samples, but also other structures such as endospores can display acid-fastness, which could explain our failure to retrieve mycobacterial sequences from these samples. Thus, tissue lesion formation was not dependent on the presence of mycobacteria. When sequences representing mycobacteria were present, they accounted for about 15% of all sequences, with the remaining sequences identifying Pseudomonadaceae and Enterobacteriaceae. Abundance estimates for organisms need to consider copy numbers of rRNA genes that vary profoundly between organisms (Stoddard et al. 2015). Members of the M. chelonae–abscessus complex, the most frequently identified mycobacteria in our tissue specimens, have just one copy, whereas the enterobacterial species Salmonella enterica or Escherichia coli potentially identified in our analyses have seven copies of the 16S rRNA gene (Pei et al. 2010). Abundance estimates correcting numbers of reads for copy numbers of rRNA genes of mycobacteria and enterobacteria could therefore increase the importance of mycobacteria to more than 50% of the bacterial community in tissue lesions, compared to 15% of the sequences in those lesions.

The presence of Houston toads with skin lesions in captive assurance colonies raised concerns regarding the introduction of novel pathogens to the environment. Acid-fast staining of bacteria in these lesions indicated the presence of mycobacteria, and causing speculation on causal relationships between mycobacteria and lesion formation. Our molecular data did not detect mycobacteria in lesions of several individuals, and not support these speculations. The diversity of bacteria in the lesions was limited, with virtually all sequences in each individual lesion represented by members of one family, the Enterobacteriaceae, the Pseudomonadaceae, the Nocardiaceae, or the Mycobacteriaceae, or combinations with one family dominating. The same sequences were detected in natal ponds regardless of head-starting status, although they were present at low abundance in highly diverse bacterial communities. An exception was pond 12 with similar abundance patterns for phyla when compared to all other ponds, but with 10-fold higher abundance of reads representing Pseudomonadaceae and Mycrobacteriaceae, and even higher abundance of reads representing Enterobacteriaceae. The family Enterobacteriaceae includes many potential pathogens and intestinal tract microbes, and thus high abundance might indicate potential small-scale fecal contamination was present in our samples. The possibility of fecal contamination should be addressed in future studies focusing on quantitative polymerase chain reaction-based analysis of members of the M. chelonae–abscessus complex or the genus Salmonella in tissue lesion and pond samples. Our Illumina sequencing data suggest that skin lesions of the captive animals are inhabited by a variety of opportunistic bacteria, with many of these bacterial taxa being normally detected in amphibian microbiomes outside of the occurrence of lesions as well. Their presence in natal ponds demonstrates that potential reintroduction even of diseased animals into their natural habitat will not change the diverse community of microbial taxa present in these habitats.

The authors are indebted to the Graduate College and the Department of Biology at Texas State University for financial support. We worked in collaboration with the Texas Parks and Wildlife Department (TPWD SPR-0102-191), the Houston Zoo, Inc., and the Capitol Area Council of the Boy Scouts of America. We are indebted to the willingness and support all of these groups provided to enable this assessment.

Alishum
A.
2019
.
DADA2 formatted 16S rRNA gene sequences for both bacteria & archaea (version 2).
Accessed April 2021.
Barrows
M,
Koeppel
K,
Michel
A,
Mitchell
E.
2017
.
Mycobacterial arthritis and synovitis in Painted Reed frogs (Hyperolius marmoratus).
J Comp Pathol
156
:
275
280
.
Bianchi
CM,
Johnson
CB,
Howard
LL,
Crump
P.
2014
.
Efficacy of fenbendazole and levamisole treatments in captive Houston toads (Bufo [Anaxyrus] houstonensis).
J Zoo Wildl Med
45
:
564
568
.
Bland
CS,
Ireland
JM,
Lozano
E,
Alvarez
ME,
Primm
TP.
2005
.
Mycobacterial ecology of the Rio Grande.
Appl Environ Microb
71
:
5719
5727
.
Blaustein
AR,
Urbina
J,
Snyder
PW,
Reynolds
E,
Dang
T,
Hoverman
JT,
Han
B,
Olson
DH,
Searle
C,
Hambalek
NM.
2018
.
Effects of emerging infectious diseases on amphibians: A review of experimental studies.
Diversity
10
:
81
.
Bucciarelli
GM,
Blaustein
AR,
Garcia
TS,
Kats
LB.
2014
.
Invasion complexities: The diverse impacts of nonnative species on amphibians.
Copeia
2014
:
611
632
.
Callahan
BJ,
McMurdie
PJ,
Rosen
MJ,
Han
AW,
Johnson
AJA,
Holmes
SP.
2016
.
DADA2: High-resolution sample inference from Illumina amplicon data.
Nat Methods
13
:
581
583
.
Caporaso
JG,
Lauber
CL,
Walters
WA,
Berg-Lyons
D,
Huntley
J,
Fierer
N,
Owens
SM,
Betley
J,
Fraser
L,
Bauer
M,
et al.
2012
.
Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms.
ISME J
6
:
1621
1624
.
Caporaso
JG,
Lauber
CL,
Walters
WA,
Berg-Lyons
D,
Lozupone
CA,
Turnbaugh
PJ,
Fierer
N,
Knight
R.
2011
.
Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample.
Proc Natl Acad Sci U S A
108
(
Suppl 1
):
4516
4522
.
Clulow
J,
Trudeau
VL,
Kouba
AJ.
2014
.
Amphibian declines in the twenty-first century: Why we need assisted reproductive technologies.
In:
Reproductive sciences in animal conservation, Adv Exp Med Biol
753
,
Holt
W,
Brown
J,
Comizzoli
P,
editors.
Springer
,
New York, New York
, pp.
275
316
.
Densmore
CL,
Green
DE.
2007
.
Diseases of amphibians.
ILAR J
48
:
235
254
.
Duarte
A,
Brown
DJ,
Forstner
MRJ.
2014
.
Documenting extinction in real time: Decline of the Houston toad on a primary recovery site.
J Fish Wildl Manag
5
:
363
371
.
Falkinham
JO.
2002
.
Nontuberculous mycobacteria in the environment.
Clin Chest Med
23
:
529
551
.
Fan
LM,
Barry
K,
Hu
GD,
Meng
SL,
Song
C,
Wu
W,
Chen
JZ,
Xu
P.
2016
.
Bacterioplankton community analysis in tilapia ponds by Illumina high-throughput sequencing.
World J Microbiol Biotechnol
32
:
10
.
Felsenstein
J.
1985
.
Confidence limits of phylogenies: An approach using the bootstrap.
Evolution
39
:
783
791
.
Ferreira
R,
de Souza Fonseca
L,
Afonso
AM,
da Silva
MG,
Saad
MH,
Lilenbaum
W.
2006
.
A report of mycobacteriosis caused by Mycobacterium marinum in bullfrogs (Rana catesbeiana).
Vet J
171
:
177
180
.
Fisher
MC,
Garner
TWJ.
2020
.
Chytrid fungi and global amphibian declines.
Nat Rev Microbiol
18
:
332
343
.
Forstner
MRJ,
Ahlbrandt
TL.
2003
.
Abiotic pond characteristics potentially influencing breeding of Houston toads (Bufo houstonensis).
Tex J Sci
55
:
315
322
.
Forstner
MRJ,
Crump
P.
2011
.
Houston toad population supplementation in Texas, USA.
In:
Global reintroduction perspectives: 2011
,
Soorae
PS,
editor.
IUCN/SSC Re-introduction Specialist Group and Environmental Agency
,
Abu Dhabi
, pp.
71
76
.
Fratzke
A,
Howard
LL,
Tocidlowski
ME,
Armién
A,
Oliveira
F,
Ritchie
B,
Berlin
E,
Snook
E.
2019
.
Chlamydia pneumoniae polioencephalomyelitis and ganglionitis in captive Houston toads (Anaxyrus houstonensis).
Vet Pathol
56
:
789
793
.
Gcebe
N,
Michel
AL,
Hlokwe
TM.
2018
.
Non-tuberculous Mycobacterium species causing mycobacteriosis in farmed aquatic animals of South Africa.
BMC Microbiol
18
:
32
.
Gilbert
JA,
Jansson
JK,
Knight
R.
2014
.
The Earth Microbiome Project: Successes and aspirations.
BMC Biol
12
:
69
.
Good
RC.
1985
.
Opportunistic pathogens in the genus Mycobacterium.
Annu Rev Microbiol
39
:
347
369
.
Gray
MJ,
Miller
DL,
Hoverman
JT.
2009
.
Ecology and pathology of amphibian ranaviruses.
Dis Aquat Organ
87
:
243
266
.
Green
SL,
Lifland
BD,
Bouley
DM,
Brown
BA,
Wallace
RJ,
Ferrell
JE.
2000
.
Disease attributed to Mycobacterium chelonae in South African clawed frogs (Xenopus laevis).
Comp Med
50
:
675
679
.
Haridy
M,
Tachikawa
Y,
Yoshida
S,
Tsuyuguchi
K,
Tomita
M,
Maeda
S,
Wada
T,
Ibi
K,
Sakai
H,
Yanai
T.
2014
.
Mycobacterium marinum infection in Japanese Forest Green Tree frogs (Rhacophorus arboreus).
J Comp Pathol
151
:
277
289
.
Hasegawa
M,
Kishino
H,
Yano
T.
1985
.
Dating of human–ape splitting by a molecular clock of mitochondrial DNA.
J Mol Evol
22
:
160
174
.
Hill
WA,
Newman
SJ,
Craig
L,
Carter
C,
Czarra
J,
Brown
JP.
2010
.
Diagnosis of Aeromonas hydrophila, Mycobacterium species, and Batrachochytrium dendrobatidis in an African Clawed frog (Xenopus laevis).
J Am Assoc Lab Anim Sci
49
:
215
220
.
Hruska
K,
Kaevska
M.
2012
.
Mycobacteria in water, soil, plants and air: A review.
Vet Med-Czech
57
:
623
679
.
Hughey
MC,
Pena
JA,
Reyes
R,
Medina
D,
Belden
LK,
Burrowes
PA.
2017
.
Skin bacterial microbiome of a generalist Puerto Rican frog varies along elevation and land use gradients.
PeerJ
5
:
e3688
.
Jeong
CY,
Ham
JH.
2017
.
Comparative analysis of the microbial community in the sediments of two constructed wetlands differentially influenced by the concentrated poultry feeding operations.
J Soil Sediment
17
:
557
566
.
Klaus
JM,
Noss
RF.
2016
.
Specialist and generalist amphibians respond to wetland restoration treatments.
J Wildl Manage
80
:
1106
1119
.
Kueneman
JG,
Parfrey
LW,
Woodhams
DC,
Archer
HM,
Knight
R,
McKenzie
VJ.
2014
.
The amphibian skin-associated microbiome across species, space and life history stages.
Mol Ecol
23
:
1238
1250
.
Li
H,
Liang
T,
Chu
Q,
Xu
F,
Li
Y,
Fu
L,
Zhou
B.
2017
.
Effects of several in-feed antibiotic combinations on the abundance and diversity of fecal microbes in weaned pigs.
Can J Microbiol
63
:
402
410
.
Lips
KR,
Green
DE,
Papendick
R.
2003
.
Chytridiomycosis in wild frogs from Southern Costa Rica.
J Herpetol
37
:
215
218
.
Mania
I,
Gorra
R,
Colombo
N,
Freppaz
M,
Martin
M,
Anesio
AM.
2019
.
Prokaryotic diversity and distribution in different habitats of an alpine rock glacier-pond system.
Microb Ecol
78
:
70
84
.
Mashiane
RA,
Ezeokoli
OT,
Adeleke
RA,
Bezuidenhout
CC.
2017
.
Metagenomic analyses of bacterial endophytes associated with the phyllosphere of a Bt maize cultivar and its isogenic parental line from South Africa.
World J Microbiol Biotechnol
33
:
80
.
McKenzie
VJ,
Bowers
RM,
Fierer
N,
Knight
R,
Lauber
CL.
2012
.
Co-habiting amphibian species harbor unique skin bacterial communities in wild populations.
ISME J
6
:
588
596
.
National Center for Biotechnology Information.
2020
.
Release 236: February 15, 2020.
Accessed March 2020.
Nogueira
CL,
Whipps
CM,
Matsumoto
CK,
Chimara
E,
Droz
S,
Tortoli
E,
de Freitas
D,
Cnockaert
M,
Palomino
JC,
Martin
A,
et al.
2015
.
Mycobacterium saopaulense sp nov., a rapidly growing mycobacterium closely related to members of the Mycobacterium chelonae–Mycobacterium abscessus group.
Int J Syst Evol Microbiol
65
:
4403
4409
.
O'Rourke
DP,
Rosenbaum
MD.
2015
.
Biology and diseases of amphibians.
In:
Laboratory animal medicine
,
Anderson
LC,
Otto
G,
Pritchett-Corning
KR,
Whary
MT,
editors.
Academic Press
,
Cambridge, Massachusetts
, pp.
931
965
.
Pearson
WR,
Lipman
DJ.
1988
.
Improved tools for biological sequence comparison.
Proc Natl Acad Sci U S A
85
:
2444
2448
.
Pei
AY,
Oberdorf
WE,
Nossa
CW,
Agarwal
A,
Chokshi
P,
Gerz
EA,
Jin
Z,
Lee
P,
Yang
L,
Poles
M,
et al.
2010
.
Diversity of 16S rRNA genes within individual prokaryotic genomes.
Appl Environ Microbiol
76
:
3886
3897
.
Pessier
AP.
2002
.
An overview of amphibian skin disease.
Semin Avian Exot Pet Med
11
:
162
174
.
Pontiroli
A,
Khera
TT,
Oakley
BB,
Mason
S,
Dowd
SE,
Travis
ER,
Erenso
G,
Aseffa
A,
Courtenay
O,
Wellington
EM.
2013
.
Prospecting environmental mycobacteria: Combined molecular approaches reveal unprecedented diversity.
PLoS One
8
:
e68648
.
Poretsky
R,
Rodriguez-R
LM,
Luo
C,
Tsementzi
D,
Konstantinidis
KT.
2014
.
Strengths and limitations of 16S rRNA gene amplicon sequencing in revealing temporal microbial community dynamics.
PLoS One
9
:
e93827
.
Primm
TP,
Lucero
CA,
Falkinham
JO.
2004
.
Health impacts of environmental mycobacteria.
Clin Microbiol Rev
17
:
98
106
.
Roguet
A,
Therial
C,
Saad
M,
Boudahmane
L,
Moulin
L,
Lucas
FS.
2016
.
High mycobacterial diversity in recreational lakes.
Antonie Leeuwenhoek
109
:
619
631
.
Saitou
N,
Nei
M.
1987
.
The neighbor-joining method: A new method for reconstructing phylogenetic trees.
Mol Biol Evol
4
:
406
425
.
Samant
S,
Sha
Q,
Iyer
A,
Dhabekar
P,
Hahn
D.
2012
.
Quantification of Frankia in soils using SYBR Green based qPCR.
Syst Appl Microbiol
35
:
191
197
.
Sanchez-Morgado
JM,
Gallagher
A,
Johnson
LK.
2009
.
Mycobacterium gordonae infection in a colony of African clawed frogs (Xenopus tropicalis).
Lab Anim
43
:
300
303
.
Silla
AJ,
Byrne
PG.
2019
.
The role of reproductive technologies in amphibian conservation breeding programs.
Annu Rev Anim Biosci
7
:
499
519
.
Stamatakis
A.
2006
.
RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models.
Bioinformatics
22
:
2688
2690
.
Stamatakis
A.
2014
.
RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies.
Bioinformatics
30
:
1312
1313
.
Stoddard
SF,
Smith
BJ,
Hein
R,
Roller
BRK,
Schmidt
TM.
2015
.
rrnDB: Improved tools for interpreting rRNA gene abundance in Bacteria and Archaea and a new foundation for future development.
Nucleic Acids Res
43
:
D593
D598
.
Tavaré
S.
1986
.
Some probabilistic and statistical problems in the analysis of DNA sequences.
In:
Some mathematical questions in biology—DNA sequence analysis, Lect Math Life Sci
17
,
Miura
RM,
editor.
American Mathematical Society
,
Providence, Rhode Island
, pp.
57
86
.
Tornabene
BJ,
Blaustein
AR,
Briggs
CJ,
Calhoun
DM,
Johnson
PTJ,
McDevitt-Galles
T,
Rohr
JR,
Hoverman
JT.
2018
.
The influence of landscape and environmental factors on ranavirus epidemiology in a California amphibian assemblage.
Freshw Biol
63
:
639
651
.
Tortoli
E,
Fedrizzi
T,
Meehan
CJ,
Trovato
A,
Grottola
A,
Giacobazzi
E,
Serpini
GF,
Tagliazucchi
S,
Fabio
A,
Bettua
C,
et al.
2017
.
The new phylogeny of the genus Mycobacterium: The old and the news.
Infect Genet Evol
56
:
19
25
.
Walke
JB,
Becker
MH,
Loftus
SC,
House
LL,
Cormier
G,
Jensen
RV,
Belden
LK.
2014
.
Amphibian skin may select for rare environmental microbes.
ISME J
8
:
2207
2217
.
Walsh
CM,
Gebert
MJ,
Delgado-Baquerizo
M,
Maestre
FT,
Fierer
N.
2019
.
A global survey of mycobacterial diversity in soil.
Appl Environ Microbiol
85
:
e01181
19
.