MYCOBACTERIA IN SKIN LESIONS AND THE HABITAT OF THE ENDANGERED HOUSTON TOAD (ANAXYRUS HOUSTONENSIS) Open Access

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 m² 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.

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 SurePrep^TM 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 MgCl₂, 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.

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).

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).

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.