Two Taricha newt species (Taricha granulosa: n = 90, Taricha torosa: n = 40) were sampled for Bd infection, skin microbiome composition, and TTX quantification from 12 distinct locations in Washington, Oregon, and California ( Figure 1 ). Sites were selected to capture a range of TTX levels, based on previous TTX quantification by Hanifin et al. (2008). Sampling occurred in February 2020 for T. torosa, and between May and July 2021 for T. granulosa. One individual T. granulosa did not have a skin punch collected for TTX concentration due to emaciation and health concerns. Mass, snout to vent length (SVL), and sex were also recorded. Sampling was performed in accordance with the following approved permits: Oregon Department of Fish and Wildlife (Permit No. 092-21), Washington Department of Fish and Wildlife (Scientific Collection Permit PIOVIA-SCOTT 20-312), California Department of Fish and Wildlife (Permit No. SC-12430), Washington State University IACUC protocol (ASAF#6339), and University of California, Los Angeles IACUC (2013-011-13 C).

Each individual newt was rinsed twice for 5 seconds with sterile deionized water to reduce transient bacteria (Lauer et al., 2007), and then skin swabs were collected using rayon-tipped sterile applicators for a standard number of strokes on the skin (Roth et al., 2013). The swab was placed into an empty 2 mL microcentrifuge tube for molecular microbiome characterization and Bd infection quantification, and for T. granulosa, a second swab was collected and stored in 20% glycerol for culturing (culturing data are not presented in this study). Samples were stored on dry ice in the field and then transferred to a -80°C freezer within 48 hours of collection (Piovia-Scott et al., 2015).

Skin samples for TTX quantification were obtained from T. granulosa by anesthetizing individuals using buffered 0.1% tricaine methanesulfonate (Ethyl 3-aminobenzoate methanesulfonate salt – MS-222) and collecting a 2.5 mm radius dorsal tissue sample between the pectoral and pelvic girdle using a sterile skin biopsy tool (Hanifin et al., 2004). After sample collection, newts were rinsed with sterile water and received an antiseptic and pain-relieving solution before final release upon anesthesia revival (Downes, 1995). Taricha torosa were sampled for TTX quantification similarly to the methods described for T. granulosa, except a 2 mm radius of dorsolateral tissue was sampled from adults using a sterile skin biopsy tool (Bucciarelli et al., 2014). Tissue samples were stored in 300 µL of 0.1 M aqueous acetic acid and held on dry ice until being transferred to a -80°C freezer for later extraction (Buttimer et al., 2022). All the sampling tools were cleaned and sterilized between individuals and replaced between sites.

Tetrodotoxin concentration was quantified for T. granulosa using a Competitive Inhibition Enzymatic Immunoassay (CIEIA) procedure as in Stokes et al. (2012). Standards were prepared using TTX-citrate (Abcam, Cambridge, United Kingdom) and ranged from 10-500 ng/mL. Values above 500 ng/mL were diluted in 1% PBS-BSA. No samples were below 10 ng/mL in concentration. The average interplate coefficient of variation was 6.84%. Taricha torosa TTX samples were quantified using a high performance liquid chromatography system coupled with fluorescence detection (HPLC-FLD) (Bucciarelli et al., 2014). It has been shown that quantification of TTX using HPLC and CIEIA techniques yield highly correlated results (Williams et al., 2016). Further, studies investigating the TTX levels of newts from the same location in Benton County, Oregon have fallen within the same range, when excluding one newt with 28 mg of TTX that doubled the mean (Stokes et al., 2015; Charles Hanifin personal communication). Although we used different methods for TTX quantification for each species, these methods are comparable to each other (Williams et al., 2016), and the levels of TTX in T. granulosa are generally higher than in T. torosa in other studies (e.g., Johnson et al., 2018). Tetrodotoxin levels of individual newts were assigned a designation of ‘high’ or ‘low’ if the estimated whole newt TTX levels were greater or less than 1 mg, using the procedure outlined in Hanifin et al. (2004) which uses a summation of estimated total skin surface area to predict whole newt TTX levels.

DNA extractions of microbiome skin swabs for T. granulosa and T. torosa were performed using a Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocols for Gram positive bacteria. This included an initial lysozyme incubation step that was conducted by adding 180 μL of a lysozyme buffer (20 mg lysozyme/mL lysis buffer) to each swab sample and incubating at 37°C for 30 minutes. Next, 25 μL proteinase K and 200 µL Buffer AL were added to each sample and incubated at 70°C for 30 minutes. The extracted DNA was used to quantify Bd infection intensity (in ITS gene copy numbers) and characterize skin bacterial communities.

Bd infection intensity was quantified using real-time TaqMan quantitative Polymerase Chain Reaction (qPCR) on the Bio-Rad CFX96 Touch Real-Time PC Detection System following the protocol established in Boyle et al. (2004). Reactions for newt samples were performed in duplicate in 96-well plates with three negative controls containing RNA-spiked water. If detection differed between the duplicate reactions, a third reaction was conducted for verification. Bd infection intensity estimates were determined for each sample using a linear standard curve generated by gBlock ITS (Integrated DNA Technologies, Coralville, Iowa) serial dilutions of 10⁶ through 10¹ performed in triplicate (while standard curve for T. granulosa samples were 10³ through 10⁰ ITS gene copies). Duplicate samples were then averaged and estimated ITS gene copy numbers (i.e., Bd infection intensity) were generated for each individual for analysis.

Newt skin bacterial community composition was characterized using amplicon sequencing of the V4-V5 region of the 16S rRNA gene with primers 515F (barcoded) and 926R for T. granulosa, and the V4 region of the 16S rRNA gene with primers 515F (barcoded) and 806R for T. torosa (Caporaso et al., 2011, 2012; Quince et al., 2011; Parada et al., 2016). While different reverse primers were used, a previous study found little difference between clustering pattern and relative abundance when comparing the primers (Walters et al., 2015). For T. granulosa, triplicate reactions per sample were performed, with each reaction containing 12 μL Qiagen UltraClean PCR-grade H₂O (Qiagen, Hilden, Germany), 10 μL Quantabio AccuStart™ II PCR SuperMix (Quantabio, Beverly, United States), 0.5 μL Illumina Forward primer/barcode 515F (Illumina, San Diego, United States), 0.5 μL Illumina Reverse primer 926R (Illumina, San Diego, United States), and 2 μL sample DNA. Reactions were processed on a Bio-Rad T100™ Thermal Cycler (Bio-Rad Laboratories, Hercules, United States) for 25 μL samples under the following conditions: 1) 94°C for 3 minutes, 2) 94°C for 45 seconds to denature the DNA, 3) 50°C for one minute to anneal primers to the DNA, 4) 72°C for 1.5 minutes to elongate the DNA, 5) Steps 2 to 4 were repeated for 35 cycles, 6) 72°C for 10 minutes, and 7) 4°C hold. Amplicons from triplicate PCR reactions were pooled and analyzed using 1% agarose gel electrophoresis stained with Thermo Fisher™ GreenGlo DNA safe stain (Thermo Fisher Scientific, Waltham, MA, United States) and 1% Tris-Borate-EDTA (TBE). Each sample had an associated negative no-template control used to confirm lack of contamination in reagents. Samples were quantitated using a Qubit™ 4.0 fluorometer (Invitrogen, Carlsbad, CA, United States), then pooled at equimolar concentrations. A Qiagen QIAquick PCR Purification Kit (Qiagen, Hilden, Germany) was used to clean the pooled, multiplexed samples prior to amplicon sequencing. Multiplexed samples were sequenced using Illumina MiSeq (250bp single-end) at the Dana Farber Cancer Institute at Harvard University. For T. torosa, duplicate reactions per sample were performed, with each reaction containing 9.5 μL Millipore Milli-Q water, 12.5 μL Azura 2x HS Taq Red Mix (Azura Genomics, Raynham, MA, United States), 0.5 μL of each primer, and 2 μL sample DNA. Reactions were processed on a Eppendorf Mastercycler Nexus Thermal Cycler (Eppendorf, Hamburg, Germany) using the following conditions: 1) 95°C for 2 minutes, 2) 28 cycles of 95°C for 15 seconds, 50°C for 45 seconds, and 72°C for 90 seconds, 3) 72°C for 10 minutes, and 4) 10°C hold. Amplicons from duplicate reactions were pooled and visualized on a 1.5% agarose gel and then purified and normalized using a Mag-Bind EquiPure Library Normalization Kit (Omega Bio-tek, Norcross, GA, United States). After normalizing, all samples were pooled for library preparation, and the library was quantitated using a QubitTM 4.0 fluorometer (Invitrogen, Carlsbad, CA, United States) and sequenced at the University of Massachusetts, Boston in the Biology Department using an Illumina MiSeq v2 cartridge (300bp single-end, trimmed to 250bp). Sequences were deposited in the NCBI SRA database under Project SUB14754928.

Newt microbiome sequencing data from each species were processed separately to tailor the filtering parameters and rarefaction depth for the quality and quantity of each species’ sequences. The datasets were then combined to investigate species differences in microbiomes. Because the two newt species had significantly different microbiomes (see Results section), all remaining microbiome analyses were performed on each newt species’ dataset separately.

For T. granulosa, Quantitative Insights Into Microbial Ecology 2 (QIIME2; v.2023.7; Bolyen et al., 2019) was used with a workflow similar to that in Walke et al. (2015). Samples (n = 89) were checked for quality score using visual plots created by FastQC (v.0.12.1. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and visualized with MultiQC (v.1.19: Ewels et al., 2016). Trimmomatic (v.0.39; Bolger et al., 2014) was used to trim sequences with a quality score less than 20 using a sliding window of 5 to 20 with reads less than a 90 minimum length dropped (Bolger et al., 2014). Sequences were imported into QIIME2 and filtered by quality using quality-filter q-score to filter sequences potentially missed by Trimmomatic (Bokulich et al., 2013). After quality filtering, 89 samples containing 9,938,507 total reads were imported (reads per sample ranged from 6 to 926,070). Samples were processed using deblur (https://github.com/biocore/deblur) to correct sequencing errors by clustering into amplicon sequence variants (ASVs) with a 0.005% mean error threshold and 20 read minimum, retaining 88 samples with 984 ASVs (total reads = 1,042,089; reads per sample ranged from 1,135 to 164,401).

Taricha torosa samples (n = 42) were processed similarly using QIIME2 (v.2021.4). Samples were quality filtered and demultiplexed with a total of 916,197 reads generated, and then processed with deblur to further quality filter and cluster reads into ASVs with a 0.005% mean error threshold and 20 read minimum, retaining 40 samples (38 newts, a PCR control, and an extraction control) with 828 ASVs (total reads = 357,000; reads per newt sample ranged from 4,246 to 13,775).

Taxonomy was assigned to the sequences using a sklearn naïve Bayes taxonomic classifier (Pedregosa et al., 2011) trained with the Silva 138 SSU Ref NR 99 (Quast et al., 2012; https://www.arb-silva.de/documentation/release-138/). A multiple sequence alignment was conducted using MAFFT (v.7.525; Katoh et al., 2002), and a 16S rRNA phylogenetic tree (rooted and unrooted) was created using IQ-TREE (v2.3.4; Nguyen et al., 2015). The deblur feature table was filtered using q2-taxa to remove ASVs designated mitochondria, chloroplast, unassigned, and Archaea. For the T. granulosa dataset, infrequently observed ASVs were identified through integration of the decontam plugin (Davis et al., 2018) using quality-control decontam-identify and removed with a threshold of 0.1 using quality-control decontam-remove, removing a total of 874 contaminant reads and 13 contaminant ASVs. The resulting dataset contained 88 samples with 882 ASVs and 1,022,255 total reads. As for the T. torosa dataset, ASVs found in the extraction and PCR controls with greater than 20 reads were removed as contaminants, resulting in 15 ASVs being removed resulting in a total of 813 ASVs and 351,423 total reads.

Samples were rarefied at 2,200 (T. granulosa) and 3,800 (T. torosa) sequences per sample to obtain an equal number of sequences per sample for each species. Rarefaction depth was chosen based on evaluation of alpha rarefaction plots for plateauing levels of diversity with sequencing depth. Samples with fewer reads than the rarefied depth were dropped, retaining 83 (93.25%) of the original 89 T. granulosa samples (875 ASVs and 182,600 reads) and 34 (87.18%) of the original 39 T. torosa samples (744 ASVs and 133,000 reads). QIIME2 was used to calculate ASV frequencies per sample, alpha diversity (Shannon’s diversity index, Faith’s Phylogenetic Diversity, observed ASVs, and Pielou’s evenness), beta diversity (Jaccard similarity, Bray-Curtis dissimilarity, unweighted UniFrac distance, and weighted UniFrac distance), and Principal Coordinates Analysis (PCoA). Two additional samples were excluded from analyses due to only a single individual being sampled at a site, resulting in 32 samples being retained. Data was exported to Rstudio (2023.1: R v.4.3.3. R Core Team, 2023) for further analysis.

To identify if any of our Illumina sequences matched sequences of bacteria known to produce TTX (Vaelli et al., 2020; Dryad Dataset), we compared the representative sequences of the rarefied ASV tables for both newt species using vsearch (Rognes et al., 2016) with a 99% identity threshold. We obtained the total reads associated with known TTX-producing microbes, the total proportion of sample reads associated with known TTX producing microbes, and the overall richness of TTX-producing ASVs found in our dataset. We also compared the ASV tables to TTX concentration to find taxa enriched in high TTX microenvironments using an Analysis of Compositions of Microbiomes with Bias Correction (see below). Similarly, we predicted the proportion of TTX producing bacteria for microbiome data from Vaelli et al. (2020), including skin regions, cloacal, and soil samples using custom scripts that match known TTX producing bacteria to the next generation sequence reads based on 99% sequence similarity in the 16S gene.

Growth inhibition assays were performed to assess if TTX could inhibit the growth of both Bd and Bsal with a dose-response curve analysis. One mg of TTX was dissolved in 1 mL of Millipore Milli-Q water and then sterilized by passing the solution through a 0.22 μm filter. Dilutions were made to test growth inhibition at the following concentrations: 250, 100, 50, 45, 40, 35, 30, 25, 20, and 15 μM. Each concentration was tested in triplicate against Bd (JEL 423) or Bsal (AMFP13/1) by performing a 96-well growth inhibition assay after Bell et al. (2013). Bd was grown on 1% tryptone and Bsal on TGhL agar plates until motile zoospores were released by mature sporangia. Zoospores were harvested by incubating plates with 3 mL 1% tryptone (Bd) or TGhL (Bsal) liquid media for 15 minutes and then normalizing to a concentration of 2 × 10⁶ zoospores/mL. Sample test wells contained 50 μL of Bd or Bsal zoospore solution and 50 μL of the appropriate TTX solution. Positive controls contained 50 μL of zoospore solution and 50 μL of 1% tryptone media for Bd, or TGhL media for Bsal. Two negative controls were added, where one contained 100 μL of appropriate media (as a media blank control to check for assay contamination), and a heat-killed control contained 50 μL of heat-killed zoospore solution (incubated at 60°C for 30 minutes) and 50 μL of the appropriate media. All reactions were set up in six replicates. Growth was measured as optical density at 492 nm on a spectrophotometer on days 0, 3, 5, and 7 for Bd incubated at 19°C, and days 0, 3, 5, 7, and 9 for Bsal at 15°C.

Since our datasets for T. granulosa and T. torosa were collected separately and in different locations, each dataset was analyzed independently. To assess whether TTX concentration varied among locations and between both sexes and species, we used a Kruskal-Wallis test. To address multiple comparisons concerning location, a Dunn’s post-hoc test with a Bonferroni correction was applied to find significant differences. To assess whether Bd infection intensity and prevalence varied among locations and between sexes and species we used a Kruskal-Wallis test, and for multiple comparisons of location we used a Dunn’s post-hoc test with a Bonferroni correction. For analyses investigating Bd infection intensity and intensity across location, only sites with three or more Bd positive individuals were included to meet sample size requirements for statistical validity.

We used Generalized Linear Mixed Models (GLMMs) with a binomial error distribution and a log link function to investigate the effect of TTX on Bd infection status for the overall species datasets and, in a separate set of analyses, Bd positive locations only. The response variable was Bd infection status and the predictor variables were TTX concentration and location; to account for non-independence of individuals from the same location, we used location as a random effect variable in all of our models. To test the significance of predictor variables within the proposed models we used Likelihood Ratio tests. To investigate the association of TTX as a response variable and microbiome alpha diversity as a predictor variable, we used Linear Mixed Models (LMMs) with location as a random variable. We used F-tests with denominator degrees of freedom estimated using Satterthwaite approximations to test for the significance of the model. We also investigated alpha diversity as a function of Bd infection intensity for all samples, Bd positive locations, and Bd positive individuals using LMMs with location as a random effect. We used linear regressions to investigate associations between putative TTX-producing ASVs and TTX concentration using individual ASVs and the relative abundance of all matching putative TTX ASVs as a single predictor for each species.

To decide whether the microbiomes of T. granulosa and T. torosa would be analyzed together as a single dataset or kept separate for analyses, a Permutational Analysis of Variance (PERMANOVA) was conducted within QIIME2 using beta diversity matrices to examine if there was a significant difference in microbiome structure between species. PERMANOVAs were conducted within QIIME2 with beta diversity matrices as the independent variable to examine similarity and dissimilarity with respect to categorical variables (TTX high/low, Bd infection presence/absence, location, sex). For continuous variables (TTX concentration, Bd infection intensity, mass), Mantel tests were conducted on beta diversity matrices within QIIME2 using the diversity beta-correlation function. We conducted PCoA on QIIME2-derived Bray-Curtis matrices for both species using the ggords (V1.18: Beck, 2017) and qiime2R (V0.99.20: Bisanz, 2018) packages and subsequently created ordination plots to visualize clustering among samples. An Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC; Lin and Peddada, 2020) was performed on the rarefied samples in respect to the TTX high classification to assess potential influential community members and their relative enrichment or depreciation in the system using the ancombc2 package (Lin and Peddada, 2020). We characterized the core microbiome for both species using post-rarefaction tables and considered any taxa present in 70% of individuals for a given species and collectively represented 75% of the total frequency to be a core microbiome member (Custer et al., 2023; Yu et al., 2023).

Percent growth of Bd and Bsal in the presence of various TTX concentrations was calculated by 1) taking the slope of the optical density over the incubation days, 2) subtracting the averaged slope of the heat-killed negative controls from the slope of each sample and positive control, 3) dividing the corrected sample slopes by the averaged slope of the corrected positive controls, and 4) multiplying by 100 to get percent growth. Differences in Bd or Bsal percent growth between multiple TTX concentrations and positive controls was determined using a one-way ANOVA and Dunnett’s post-hoc test.

Mean TTX concentration was greater in T. granulosa than T. torosa (χ² = 8.76, p = 0.003). TTX concentration for T. granulosa individuals ranged from 0.061 mg to 6.718 mg, with a mean of 1.94 mg and standard deviation of ± 2.24 (males: mean ± SD = 2.26 ± 2.36; females: mean ± SD = 1.25 ± 1.78), and was significantly different among locations (χ² = 63.405, p < 0.0001; Figure 2 ). For T. granulosa, there were more males than females in our samples (males: n = 61, females: n = 29), but TTX concentration was not significantly different between newt sexes (χ² = 2.61, p = 0.11). TTX concentration for T. torosa ranged from 0.0025 mg to 1.83 mg, with a mean of 0.47 mg and a standard deviation of ± 0.46 (males: mean ± SD = 0.42 ± 037; females: mean ± SD = 0.58 ± 0.63), and was significantly different among locations (χ² = 17.56, p < 0.001; Figure 2 ). For T. torosa, there were also more males than females in our samples (males: n = 26, females: n = 11), and TTX concentration was not significantly different between sexes (χ² = 0.009, p = 0.92).

Batrachochytrium dendrobatidis prevalence was 40.9% overall (51/127 total individuals; T. granulosa = 37/90, 41%; T. torosa = 15/37, 40.54%). Bd infection was present at 70% of sampled locations (7/10 locations overall, T. granulosa = 5/6, T. torosa = 2/4), although one of the five Bd-positive sites for T. granulosa had only one positive individual. Among infected individuals, Bd infection intensity ranged from 45.8 to 23,038.6 ITS copies, with a mean of 2,105.8 ITS copies and standard error of ± 565 (T. granulosa: mean ITS copies ± SD = 1738.3 ± 2899.3: Figure 3A ; T. torosa: mean ITS copies = 1.14 ± 4107.9. Figure 3B ). There were no differences in Bd infection intensity or status between newt species (all samples included: status: χ² = 0.004, p = 0.95, intensity: χ² = 0.26, p = 0.61; Bd-positive samples only: intensity: χ² = 2.33, p = 0.13). Bd infection status was not significantly different between newt sexes for either species when accounting for differences in location (T. granulosa: df = 1, Coef = -0.31, χ² = 0.38, p = 0.54; T. torosa: df = 1, Coef = 2.96, χ² = 1.23, p = 0.27). Bd infection intensity and status were significantly different among locations for both species when all individuals were included (T. granulosa: intensity: χ² = 14.61, p = 0.01, status χ² = 17.05, p < 0.01: Figure 3C ; T. torosa: intensity: χ² = 27.81, p < 0.0001; status: χ² = 27.81, p < 0.001. Figure 3D ). Bd infection intensity was still significantly different among locations when only locations with confirmed Bd infection were included for T. torosa (χ² = 7.04, p < 0.01), but not T. granulosa (χ² = 7.81, p = 0.10). Similarly, when only individuals with confirmed Bd infection were included, Bd infection intensity was significantly different among sites for T. torosa (χ² = 27.81, p < 0.0001), but not T. granulosa (χ² = 6.10, p = 0.19). Bd infection status was significantly different among locations with confirmed Bd infection for both species (T. granulosa: χ² = 9.89, p = 0.04; T. torosa: χ² = 6.60, p = 0.01).

Bd infection status was not significantly associated with TTX concentration for either newt species when all individuals were included (T. granulosa: df = 1, Coef = 0.24, χ² = 2.16, p = 0.14; T. torosa: df = 1, Coef = 0.07, χ² = 0.0004, p = 0.98 Figure 4 ), nor when only Bd positive locations were included (T. granulosa: df = 6, Coef = 0.17, χ² = 1.42 p = 0.23; T. torosa: df = 1, Coef = 0.20, χ² = 0.07, p = 0.94).

The two Taricha newt species had significantly different microbiome structure and composition ( Figure 5 ; PERMANOVA on Weighted-Unifrac: Pseudo-f = 12.59, p = 0.001); thus, further microbiome analyses were conducted separately for each newt species. The most abundant taxa for T. granulosa species were Gammaproteobacteria (48.6%), Bacteroidia (24.9%), and Alphaproteobacteria (13.3%), while the most abundant taxa for T. torosa were Gammaproteobacteria (61.6%), Bacilli (20.5%), and Bacteroidia (12.3%) ( Figure 5 ). Five ASVs were identified as core members for T. granulosa, but no core ASVs were identified for T. torosa ( Table 1 ).

Skin microbiome alpha diversity for both T. granulosa and T. torosa were significantly different across sites, with the exception of Pielou’s evenness in T. torosa, which was borderline significantly different across sites ( Table 2 ). Alpha diversity was not significantly different between Bd infected and uninfected newts for either species, nor was Bd infection intensity correlated with alpha diversity ( Table 2 ). For TTX levels (high/low), alpha diversity was not different between levels, and this pattern was the same in both newt species ( Table 2 ). For TTX concentration, no alpha diversity metrics were correlated with TTX in T. granulosa, but ASV richness and Faith’s Phylogenetic Diversity were negatively correlated with TTX in T. torosa ( Table 2 ). The ANCOM test identified six ASVs which had an enriched relative abundance in high TTX T. granulosa, but found no microbes that were enriched or depreciated for T. torosa ( Figure 6 ). Female newts had higher Faith’s Phylogenetic Diversity than males for both T. granulosa and T. torosa ( Table 2 ). Mass was not correlated with alpha diversity for either newt species ( Table 2 ).

Newts sampled from different locations had significantly different microbiome structure; this was the case for both newt species and for all beta diversity metrics ( Figures 7A, D ; Table 3 ). No newt skin microbiome structures were correlated with TTX concentrations in either species ( Figures 7B, E ; Table 3 ); however skin microbiome structure was correlated with TTX concentration classified as high and low for T. granulosa, but not T. torosa ( Table 3 ). On the other hand, newt skin microbiome structure was correlated with Bd infection intensity in T. torosa, but not in T. granulosa ( Figures 7C, F ; Table 3 ). T. torosa newts also had significantly different microbiomes depending on Bd infection status (present/absent) ( Table 3 ). Male and female newts had similar microbiome structure in both species, and lastly, body mass was correlated with at least one beta diversity metric in both species ( Table 3 ).

Matching our Illumina sequences to the Vaelli et al. (2020) dataset of TTX-producing bacteria identified four Pseudomonas spp. found on T. granulosa, and three Pseudomonas spp., one Aeromonas sp., and one Shewanella sp. found on T. torosa ( Table 4 ). These putative TTX-producing bacteria represent an overall relative abundance of 6.05% and 1.32% on T. granulosa and T. torosa, respectively. TTX concentration in T. granulosa was negatively correlated with the relative abundance of putative TTX-producing bacteria (t-value = - 2.85, p < 0.01), but was not significantly correlated with the richness of TTX-producing ASVs (t-value = -1.44, p = 0.15). For the four matching Pseudomonas ASVs, only one exhibited a negative correlation with TTX concentration (TX174011: Coef = -0.06, t-value = -2.47, p = 0.02), while the other three matches exhibited no correlation with TTX concentration ( Table 5 ). In T. torosa, TTX concentration was not significantly correlated with either relative abundance of putative TTX-producing bacteria (t-value = -0.52 p = 0.61) or number of TTX-producing ASVs, richness (t-value = -0.49, p = 0.62). For the six matching ASVs, none exhibited a significant correlation with TTX concentration ( Table 5 ).

Overall, there were low levels of predicted TTX-producing bacteria in the microbiomes from Vaelli et al. (2020) samples ( Supplementary Figure S1 ). The soil samples had nearly no predicted TTX-producing bacteria compared to that of Taricha body site samples.

TTX did not inhibit the growth of Bd (F_10,41 = 1.248; P = 0.291) nor Bsal (F_10,41 = 1.807; P = 0.090) in the growth inhibition assays. From the Dunnett’s Post Hoc test, no concentration of TTX significantly affected the growth of Bd and Bsal ( Supplementary Figure S2 ; Supplementary Table S2 ).