Field collections

During the field seasons of 2015–2017, samples were collected from 35 sites across six states in the Southern United States: Alabama, Arkansas, Florida, Georgia, Tennessee, and Texas (Fig. 1). We defined “macroscale” patterns of microbial assemblages as those observed in the Southern United States. Sampling emphasis focused on four Tennessee ecoregions (“mesoscale”), incorporating 28 sites across the state (“microscale”), to provide a finer scale understanding of the snake skin microbial assemblage (Fig. 1). All snake individuals encountered were captured and handled with clean nitrile gloves and transient microbes were removed by rinsing with 100 mL of sterile deionized water [28, 29] autoclaved for 2 h [30]. A single skin swab was collected using a sterile rayon-tipped swab (Puritan, VWR cat #10808-146) for later DNA extraction. Each snake was swabbed in a way to standardize “grain size” by taking 15 swab strokes along a 15 cm length of the middle-third body section of the animal, encompassing the ventral, dorsal, and lateral body sites. No snakes were collected or sampled that failed to exceed swabbing length.

Fig. 1 Host species, geography, season, and pathogen presence are predictive of the microbiome across macro- to microspatial scales. Snake sampling sites and presence of O. ophiodiicola are marked with colored circles across the map. Statistically significant explanatory variables from adonis models are listed with the corresponding coefficient of determination (R2) Full size image

Environmental sampling

To assay the environment for O. ophiodiicola, soil samples (n = 48) were procured from the exact point of snake capture, excluding high traffic roads/paths. Soil samples were obtained by aseptically swabbing the soil in concentric circles not exceeding 30 cm in diameter with a rayon-tipped sterile applicator. Samples were immediately stored in dry 2 mL microcentrifuge tubes. To assay the environment for O. ophiodiicola, water samples (n = 29) were collected at the nearest accessible water source (e.g., stream bank adjacent to point capture) in proximity to captured aquatic snakes by aseptically filtering 500 mL of water through a Thermo Scientific Nalgene Analytical Test Funnel (CN 145-2020) with a measured pore size of 0.2 μm. Water was drawn through the filter by a peristaltic pump and the filter paper was aseptically placed into a dry 2 mL centrifuge tube. All three sample types (skin, soil, and water) were immediately stored at −10 °C in a CoolFreeze CF-25 vehicle freezer and transferred directly to a lab freezer at −20 °C until further processing.

DNA extraction and quantitative PCR

DNA was extracted from filter paper (1/4 slice), skin, and soil swab samples using the Qiagen DNeasy PowerSoil HTP 96 kit, per the manufacturer’s standard protocol. A single DNA extraction negative control blank was extracted on each plate of 96 samples or as a single preparation when samples exceeded the 96-well plate. All negative controls were screened using quantitative PCR (qPCR) and/or high-throughput DNA sequencing. To minimize cross-contamination, DNA extraction, qPCR, and library setup were done in separate AirClean Systems AC600 (AirClean Systems, Creedmoor, NC) dedicated to these processes. Each workstation had a dedicated set of pipettes that were routinely autoclaved before experimentation. DNA was concentrated to ≈20 µl using a Thermo Fisher Savant DNA SpeedVac. The molecular presence of O. ophiodiicola in all skin and environmental samples was detected using qPCR. The qPCR assay as described by Bohuski et al. [31] was used to amplify the internal transcribed spacer region of the rRNA gene of O. ophiodiicola. Each 96-well qPCR sample plate was processed with a positive control and two no-template negative control reactions. Sample reactions were run in triplicate and consisted of a 10 μL volume containing 5 μL IDT 2× Primetime MasterMix, 0.4 μL of IDT forward primer (10 μM), 0.4 μL of IDT reverse primer (10 μM), 0.1 μL (20 μM) IDT probe, 2.1 μL PCR-grade water, and 2 μL of DNA template. Thermocycling conditions consisted of a 3 min cycle at 95 °C, followed by 50 cycles of 95 °C for 10 s and 60 °C for 30 s. A positive test was confirmed by an exponential phase appearing at a C t < 39 in each reaction. Samples showing incomplete evidence by testing positive for one or two replicates during the first run were reanalyzed on a separate plate. Samples still showing ambiguous results (one or two positive reactions) after the second round of amplification were considered positive [32,33,34], whereas samples not amplifying in triplicate on the second plate were considered negative. To calculate pathogen load from resulting C t values, a single standard curve was run using a serial dilution of 1 × 1010 − 1 amplicon copies of a synthetic gBlocks (Integrated DNA Technologies) fragment identical to the target qPCR region in Bohuski et al. [31]. The standard curve was used to generate the formula, y = −0.2915 × +11.094, where “x” is average C t value of each sample, and this was used to calculate the log copy number (herein, “pathogen load”) of the ITS region of the rRNA gene of O. ophiodiicola per qPCR reaction. All C t values and resulting pathogen load values are reported in Supplementary Data File 1.

High-throughput sequencing and bioinformatics analyses

A subset of 168 skin swabs, ten soil samples from TN, and ten water samples from TN were selected for high-throughput sequencing. These samples represented a cross-sectional view of samples collected across 35 field sites and all three years. Sequencing was performed according to the Illumina 16S Metagenomic Sequencing Library Preparation protocol on the Illumina MiSeq instrument in two separate sequencing runs. The V4 region of the 16S rDNA was amplified using dual indexed fusion primers as described by Kozich et al. [35]. Samples were cleaned using Ampure XP magnetic beads after both the initial PCR and index PCR step. After amplification and indexing, PCR products were quantified on a Qubit fluorometer 3.0, per the manufacturer’s protocol, and visualized for amplicon size (≈450 bp) on an Agilent 2100 Bioanalyzer according to the DNA 1000 protocol, and then normalized. After library quality control and quantification, the library was loaded on an Illumina MiSeq v3 flow cell and sequenced using a 500-cycle reagent kit (paired-end 2 × 300 reads). A total of 13,908,315 raw data sequence reads were obtained during run one and 10,137,882 in run two. Data were processed per the MiSeq SOP described by Kozich et al. [35] using mothur v1.39.5. After assembling paired-end reads into contigs, sequences were removed from the analysis if they had fewer than 249 bp or >253 bp, contained homopolymers in excess of eight nucleotides, or contained any ambiguous base calls. Of these, unique sequences were aligned to the SILVA v.123 bacterial reference database [36] curated to the V4 region of interest, pre-clustered allowing for two nucleotide differences, and chimeras removed per the mothur UCHIME algorithm [37]. Remaining alignments were classified into taxonomic lineage using classify.seqs at an 80% cutoff value, and lineages identified as Archaea, Eukaryota, chloroplast, mitochondria, and unknown were removed [35]. Sequences passing all quality control metrics were clustered into operational taxonomic units (OTUs) using cluster.split at 97% sequence similarity [38]. Rare OTUs appearing five times or fewer throughout the dataset were removed [39]. In addition, 64 OTUs appearing in the negative control DNA extraction blanks (1,376 reads) were removed from the dataset. A total of 8,654,686 sequence reads were kept after filtering and utilized in all downstream analyses. Ten skin swab samples were removed (n = 178 total samples retained) due to inadequate coverage after data filtering and curation in mothur. To normalize, we ran summary.single in mothur to compute sample coverage, then subsampled at 2,101 sequence reads per sample, and used this dataset in statistical analyses. Code to reproduce the bioinformatics analysis is provided in the Supplementary “mothur code” File.

Statistical analyses—comparison of the skin with environment

All analyses and graphing were conducted in mothur v1.39.5 and R version 3.5.1 using the packages vegan [40], plyr [41], dplyr [42], ggplot2 [43], rcompanion [44], car [45], and gridExtra [46]. To test whether the snake skin microbial assemblages differed from the environment, the skin, soil, and water bacterial OTU beta diversity was calculated using the vegdist function to generate a Bray–Curtis dissimilarity matrix representing sample-to-sample pairwise distances, and these distances were further analyzed using the metaMDS function to generate a non-metric multidimensional scaling (nMDS) ordination. The adonis function was used to perform a permutational multivariate analysis of variance using 999 permutations on the Bray–Curtis dissimilarity matrix, to determine whether skin, soil, or water were explanatory variables for OTU assemblages.

Comparisons of the skin microbial assemblages across spatial extents

As the skin assemblage of snakes differed from the environmental microbial assemblages (see Results below), we used an indicator analysis (indicator values > 30, p < 0.05) to select a subset of 56 OTUs that were descriptive of the variation in the snake skin microbial assemblages across all collected skin samples. To understand the effect of spatial scale on the snake skin microbial assemblages, we subdivided our dataset into three categories for analysis including the following: (1) broad-scale patterns (macroscale) in the Southern United States (n = 158 samples); (2) mesoscale patterns across four geographically distinct Tennessee ecoregions (n = 124 samples); (3) microscale patterns of between site variation within each Tennessee ecoregion (Fig. 1). For the macroscale dataset, all samples were coded as occurring in a geographic region (Alabama, Arkansas, Florida, Georgia, Tennessee, or Texas). Tennessee ecoregion (mesoscale) factors included the blue ridge mountains, interior plateau, ridge and valley, and Southwestern Appalachians. Microscale categories included four sites within the blue ridge mountains (n = 16 skin samples), 13 sites within the interior plateau (n = 73 skin samples), and 8 sites within the Southwestern Appalachians (n = 32 skin samples). Site-based differences within the ridge and valley ecoregion were not analyzed due to a small sample size effect (three sites, three samples). The adonis function was used to perform a permutational multivariate analysis of variance using 999 permutations on the Bray–Curtis dissimilarity matrix of 56 indicator OTUs to test whether the explanatory variables including host species, geographic location, season, year, O. ophiodiicola status (positive/negative), or interactions between these factors were predictive of microbiota diversity at micro-, meso-, and macroscales. We also evaluated alpha diversity as a tentative variable component of snake skin microbial assemblages by comparing inverse Simpson diversity values for the explanatory variables host species and O. ophiodiicola status (positive/negative) using a Kruskal–Wallis test across micro-, meso-, and macroscales.

Host life history and the skin microbiota

“Ecomode” is defined as modal categories (e.g., arboreal, aquatic, fossorial, and terrestrial) of habitat use in species that do not necessarily display convergent morphology [47]. To determine whether snake life history was predictive of microbial assemblages, all snake species were categorized by ecomode including aquatic, arboreal, fossorial, or terrestrial, and an nMDS ordination was run using the metaMDS function on Bray–Curtis dissimilarity matrix values for each micro- to macroscale dataset of the 56 previously described indicator OTUs. The adonis function in vegan was used to perform a permutational multivariate analysis of variance using 999 permutations on the Bray–Curtis dissimilarity matrix of indicator OTUs to test whether the explanatory variables including ecomode, O. ophiodiicola status, or an interaction between these factors was predictive of microbiome diversity at micro-, meso-, and macroscales. We assessed the alpha diversity of ecomode categories using a Kruskal–Wallis test on inverse Simpson diversity values across micro-, meso-, and macroscales.

Pathogen load and the skin microbiota

To understand the structure of the microbiota in the presence of a fungal pathogen across space, we performed indicator analyses (indicator values > 30, p < 0.05) in mothur to select OTUs descriptive of the skin microbial assemblage in the presence/absence of O. ophiodiicola independently across macro-, meso-, and microscales. We conducted a constrained analysis of proximities using the capscale function in vegan to model the continuous explanatory variable, O. ophiodiicola pathogen load (qPCR results), on Bray–Curtis dissimilarity values from microbial assemblages of O. ophiodiicola-positive/negative snakes. Both O. ophiodiicola pathogen load/reaction and OTU relative abundances were cube root transformed before the capscale analysis. We modeled the correlation between O. ophiodiicola pathogen load and microbial assemblages for 16 indicator OTUs at the macroscale, 17 OTUs at the mesoscale, 12 OTUs for the interior plateau (microscale), and 12 OTUs for the Southwestern Appalachian mountain ecoregion (microscale). We did not model the described associations for the ridge and valley ecoregion, as we only collected one O. ophiodiicola positive snake from this location. The complete list of samples processed, taxonomic identifications, and statistical models is presented in Supplementary Tables 1–3. See Grisnik et al. [48] for information about clinical signs for a subset of snakes used in these analyses.