Animal Collection & Husbandry
Following breeding events in March 2016, partial egg masses (n = 7) of L. sphenocephalus were collected from three populations at Edward J. Meeman Biological Field Station (MBFS, The University of Memphis) in Shelby County, Tennessee, USA (35°23’22.66” N, 90°02’15.75” W). Eggs were transported to the MBFS laboratory and kept in aquaria with mesocosm water (described below) with oxygen bubblers until hatching. Upon hatching and reaching free-swimming stage (Gosner 25; [54]), individuals were fed a 3:1 mixture of crushed Kaytee rabbit chow and Tetramin® tropical flakes ad libitum with a 14:10 hour light:dark cycle. After 2 weeks, tadpoles were transferred to randomized outdoor mesocosm tanks (n = 15 per tank) at MBFS. The mesocosms consisted of polyethylene tanks filled with water to a depth of 30.5 cm (~ 613L), 300 g of dried leaf litter (primarily Quercus species), and a 100-ml dose of concentrated zooplankton suspension collected from a local pond [55]. Mesocosm algal and zooplankton communities were given two full weeks to develop before the additions of tadpoles.
At Gosner 42 (emergence of front limbs), individuals were removed from mesocosms and maintained in individual 1.5 L plastic containers with autoclaved sphagnum moss to maintain moisture. Containers were cleaned twice a week (replaced soiled water) and frogs were fed calcium-dusted crickets. Until experimental exposures began, frog were maintained at 19 ºC on a 12:12 light:dark hour photoperiod.
Btk culture preparation
To prepare the bacterial bath used in this study, Bacillus thuringiensis subsp. kurstaki was sub-cultured from Monterey B.t. liquid biocide (Lawn and Garden Products, Fresno, CA) with an inoculating loop and streaked onto 1% tryptone agar plates. After incubating plates at room temperature for 48 hours, visible colonies were transferred to 5 ml tubes of autoclaved 1% tryptone broth. Broth cultures were incubated at 30 °C on an orbital shaker for 28-hrs to reach stationary growth phase (OD600nm = 0.84 ± 0.01, OD targets based on preliminary research). Following incubation, 1 mL of the liquid cultures were centrifuged at 4500 × g for 10 minutes, the supernatant was discarded, and the cells were resuspended in sterile water. This process was repeated 3 times. To enumerate colony forming units (CFU), a 0.1 ml aliquot of 10− 8 diluted stock culture was plated onto 1% tryptone agar with an agar overlay. Following 48-hr incubation at 30 °C, colonies were quantified from 3 plates and averaged. The final concentration of the stock solution used in exposures was 3.4 × 1011 CFU/ml.
Bd preparation
Batrachochytrium dendrobatidis culture was prepared using a strain (FMB 003; [56]) isolated from a local, infected L. sphenocephalus adult in 2010 at Meeman-Shelby State Park, Shelby County, Tennessee, USA. Plates were prepared on 1% tryptone agar from 2 ml of homogenized Bd stock culture, sealed, and incubated at room temperature for 10 days. Following incubation, plates were flooded with 3 ml of aged tap water for 45 minutes to harvest zoospores. The zoospores were pooled from the plates and counted using a hemocytometer and the final concentration used was 2 × 106 zoospores/ml.
Btk and Bd exposures
L. sphenocephalus adults (33 weeks post-metamorphosis) were randomly assigned into four treatment groups: N (negative control; n = 6), Btk (B. thuringiensis kurstaki exposure; n = 8), Bd (B. dendrobatidis exposure; n = 7), BB (Btk + Bd; n = 10). Individuals were rinsed with 30 ml of autoclaved aged tap water to remove transient bacteria and placed into new bleach-sterilized plastic containers (1.5 L) with autoclaved sphagnum moss and 20 ml of autoclaved aged tap water. During the experiment, individuals were removed for weekly cage cleanings, but remained segregated in new zip top plastic bags (1 qt). Crickets were added to the container following cage cleaning.
Individuals in treatment groups Btk and BB were exposed to two inoculation baths of Btk. Individuals in the BB group were later exposed to Bd. The first Btk exposure, on Day 0 (see Fig. 1), was diluted from the stock solution (3.4 × 1011 CFU/ml) to 1.7 × 1011 CFU/ml by adding 5 ml of stock solution and 5 ml of sterile water to sterile 120 ml sample cups with air holes in the lid. The second exposure (Day 7) was diluted from the same stock solution. Due to limited stock volume, 4 ml of Btk stock was combined with 11 ml of sterile water for a final exposure concentration of 9.1 × 1010 CFU/ml. Individuals in the Bd and N groups were exposed to the same volume of sterile water. After 24-hrs in the bath, individuals were placed back into their containers. We used two inoculation baths to mimic the multiple exposures that amphibians experience in nature. Following Btk exposures, frogs were not disturbed for one week, aside from cage cleaning, to allow the potential for Btk colonization.
For Bd exposures, individuals in the Bd and BB groups were exposed to 6 × 105 zoospores/ml of Bd in 30 ml for 24 hours (Day 14; Fig. 1). Frogs in the N and Btk group were exposed to water collected and diluted from 1% tryptone agar plates as a negative control for the residual tryptone media. Mass was measured while individuals were inside of zip top plastic bags with a Pesola Spring scale (20 g ± 0.2) and snout-vent length (SVL) was measured with SPI Dial Calipers (± 0.1 mm).
Sampling and DNA extractions
To assess cutaneous bacterial communities and Bd infection status, samples were taken via sterile cutaneous swabbing. During swab collection (see Fig. 1), each frog was handled with a fresh glove and swabbed 5 times each on the ventral posterior patch, dorsal surface, and ventral side of both hind legs with a sterile cotton swab. A total of 6 swabs were taken for each frog (n = 186) over the experiment at various timepoints. These include: Day 0 (Timepoint 1), Day 3 (Timepoint 2), Day 7 (Timepoint 3), Day 11 (Timepoint 4), Day 23 (Timepoint 5), and Day 29 (Timepoint 6). Swabs were then stored at -20 °C in 1.5 ml sterile microcentrifuge tubes until DNA was extracted using the DNeasy Blood and Tissue DNA Extraction Kit (Qiagen, Valencia, CA USA, animal tissue protocol). DNA was stored at -20 °C until molecular work was conducted.
NGS Library Preparation and Sequencing
Sequencing libraries were generated by selectively amplifying the bacterial 16S (V4) region using a two-step amplification process (following [57]). Briefly, the V4 region of the 16S rRNA gene repeat was amplified using the primer pairs nexF-N[3–6]-515f and nexR-N[3–6]-806r where 515f and 806r [58] are bacteria gene primers, N[3–6] represents four identical primers with the exception of containing a range of ambiguous nucleotides (3–6) mixed to equal molarity to increase nucleotide diversity during sequencing, and nexF and nexR are Nextera forward and reverse sequencing primers. PCRs were conducted in 20 µL reactions using 2 µL extracted DNA, 4 µL 5X Phusion High-fidelity Buffer, 200 µM each dNTP, 1 µM of each forward and reverse primer, 0.2 µL Phusion HotStart II DNA Polymerase (ThermoFisher Scientific; Waltham, MA, USA), and 7.8 µL molecular grade H2O. Primary PCR parameters were 98° C for 30 seconds, 25 cycles of 98° C for 20 seconds, annealing temperature for 30 seconds at 52.5° C, and 72° C for 40 seconds, followed by a final extension at 72° for 10 minutes, all ramp rates were 1° C/second (SimpliAmp Thermal Cycler, Applied Biosystems, Foster City, CA, USA). This resulted in a final 1° PCR construct of nexF-N[3–6]-primer--primer-N[3–6]-nexR. After primary PCR generation, secondary PCR reactions were conducted in 25 µL reactions using forward primers that include the P5-i5-overlap and the reverse primers P7-i7-overlap where P5 and P7 are the Illumina Adaptor sequences, i5 and i7 are 8 bp unique Molecular Identifiers (MIDs – barcode), and the overlap is the partial nexF and nexR sequences that acts as the annealing site for the 2° PCR. The forward and reverse barcoded 2° primers were mixed in a combinatorial fashion to generate unique dual barcoded primers (see Table S1 for primer and MID information) in a working concentration of 10 µM (5 µM for each primer). The 2° PCR reactions contained 2.5 µL of 1° PCR product, 5 µL 5X Phusion High-fidelity Buffer, 200 µM each dNTP, 0.5 µM of each forward and reverse primer, 0.25 µL Phusion HotStart II DNA Polymerase (0.02 U/µL final concentration; ThermoFisher Scientific), and 7.5 µL molecular grade H2O with the PCR parameters of 95° C for 2 minutes, 8 cycles of 95° C for 20 seconds, 50° C for 20 seconds and 72° C for 50 seconds, followed by a final extension at 72° for 10 minutes. This produced the final amplicon constructs of P5-i5-nexF-N[3–6]-primer--primer-N[3–6]-nexR-i7-P7 using a total of 32 cycles.
Secondary PCR products were cleaned using Axygen AxyPrep Mag PCR clean up beads (Axygen Biosciences, Union City, CA, USA) following kit protocol with the modification using a 1:1 bead solution to reaction volume ratio to better select against small fragments and primer-dimers [59]. Cleaned PCR products were quantified using Qubit 3.0 fluorometric assays (dsDNA HS Assay Kit; ThermoFisher Scientific). PCR products were pooled into a library at equal concentrations and sequenced on one Illumina MiSeq (v.3, 300PE) at the Kansas State University Integrated Genomics Facility (Manhattan, KS, USA). Demultiplexing of the raw sequence data using the unique i5 and i7 sequence combinations provided individual paired fastq files. Sequence data is deposited in SRA at NCBI (BioProject xxxxxxxx, BioSample xxxxxxxxx).
Btk reference sequencing and Bd DNA infection confirmation
To determine which operational taxonomic unit (OTU; see below) corresponds to the strain of Btk used in exposures, a representative sequence from Btk was obtained. DNA from Btk cultures were extracted using the DNeasy Blood and Tissue DNA Extraction Kit and PCR was conducted in a 50 µL reactions targeting the partial 16S region (357F and 1100R; [60]), which fully encompasses the V4 region used in our community analyses. PCR parameters consisted of 98° C for 30 seconds, 25 cycles of 98° C for 20 seconds, annealing temperature for 30 seconds at 58° C, and 72° C for 40 seconds, followed by a final extension at 72° for 10 minutes. PCR product was visualized via gel electrophoresis (1% agarose in TBE) and cleaned using DNA Clean & Concentrator kits following protocols (Zymo Research, Irvine, CA, USA and Sanger sequenced (10 ng DNA, and 0.5 µL or each primer [10µM]) at the University of Tennessee Health Science Center (Memphis, TN, USA) and the representative sequence is deposited in GenBank (accession XXXXXX).
To assess Bd infection status, real time quantitative PCR (qPCR; Bio-Rad CFX96 Touch Real-Time PCR Detection System) was used to quantify Bd DNA loads from each sample (following [41]). PCRs were conducted in 20 µL reactions with 900 nM of the PCR primers ITS1-3 Chytr and 5.8S Chytr [61], 240 nM MGB probe, TaqMan Exogenous Internal Positive Control, and 250 nM DyNamo Flash Probe. The qPCR protocol for DyNamo Flash Probe was used to set PCR parameters (Initial denaturation: 95 °C for 7 min, Annealing/Extension: 60 °C for 30 seconds) and included Bd DNA and sterile water for positive and negative controls, respectively. A plasmid standard dilution series (Pisces Molecular, Boulder, CO, USA) was used to quantify zoospore genomic equivalents and all reactions were run in duplicate. We used regression analysis to assess any associations between the Btk abundance (OTU72 - see below) and Bd genomic equivalents.
Bioinformatics
Illumina Sequences were processed using the program mothur (v.1.40; [62]). The forward and reverse sequences were contiged and screened to cull any sequences with ambiguous bases, or greater than 10 homopolymers and merged into a single fasta file and trimmed to eliminate primer sequences. Sequences were aligned against the SILVA (release 132) reference alignment and filtered to exclude non 16S V4 regions. Sequences were preclustered to remove basepair variation due to sequence chemistry errors (using pseudo-single-linkage clustering following [63] as implemented in mothur), screened for chimeras (mothur implemented VSEARCH; [64]), and putative chimeras were culled. Sequences were screened for off-target amplification (non-bacterial in origin) by classifying all sequences using a mothur implemented Naïve Bayesian Classifier [65] against the RDP training set (v.10). Non-target lineages were culled and distance matrices (fully aligned distances not punishing terminal gaps) were generated. Sequences were clustered into operational taxonomic units (OTUs) at a 3% dissimilarity threshold using OptiClust [66]. OTUs with fewer than 10 sequences globally were eliminated to prevent inclusion of potentially spurious OTUs into analyses [67–68]. After all sequence quality control 2,343 OTUs were retained. Negative control (ddH2O) were included throughout extraction and amplicon generation and remained free of visible contamination during amplification, were eliminated completely during these sequence cleanup steps.
Diversity Indices
Diversity estimates based on an iterative subsampling approach using 1000 iterations at a subsampling depth of 5,635 sequences per sample (without replacement - selected to retain all samples) were generated (in mothur) and the average estimates were used for all downstream analyses. These include observed OTU richness (Sobs), the complement of Simpson’s diversity (1-D; 1 - ∑pi2) where pi is the frequency that each OTU occurred in each sample, Simpson’s evenness (ED: (1/D)/Sobs) where D is Simpson’s. Diversity estimators were analyzed to examine if treatments impact microbiome diversity using repeated measures analysis of variance (ANOVA) with Kenward-Roger first order approximations with Kacker-Harville corrections which allows for partial degrees of freedom.
Community Dynamics
To test if bacterial communities differ across treatments, time, and their interactions, a permutational multivariate analysis of variance (PERMANVOA; [69]) of average Bray-Curtis dissimilarity values (iteratively subsampled as above) was conducted in the program R (v.3.3.3) with the package vegan (function adonis with 999 iterations, strata = individual frog to facilitate repeated measures, [70]) and post-hoc multiple comparisons were examined using the package RVAIDeMemoire (function pairwise.perm.manova with FDR corrections, 999 iterations, [71]). To visualize communities, nonmetric multidimensional scaling (NMDS) was conducted (using mothur) using 1000 iterations and optimally resolved to across 4 axes (4D stress = 0.1544). Further, to examine community shifts on individual frogs after perturbation events to test community recovery over time, we used the AWOrD metric (Axis Weighted Ordination Distance; [72]). The AWOrD quantifies distance in ordination space (NMDS here) between any two samples across multiple axes whilst accounting for ordination axes coefficients of determination (scaled by R2 for each axis) and is based on a modified Manhattan distance. If AWOrD values are lower, then two samples are highly similar whereas if values are higher, they are more dissimilar. We calculated AWOrD values for the four solved NMDS axes between sampling timepoint one (T1) and all subsequent sampling timepoints for each individual frog and using two-way ANOVAs, we tested if AWOrD values differed between treatments, across time (T1 vs T2, T1 vs T3, …, T1 vs T6), and treatment by time interactions.
Further, we sought to examine patterns of common core taxa (OTUs) across our experimental framework. To identify core taxa, we compiled a list of OTUs that are present in ≥90% of the frogs in our experiment. Core OTUs were tested using repeated measures ANOVA (relative abundance, logit transformed, sampling effort [time] as a categorical variable) to examine if these core taxa change in abundance across treatments, time, and their interactions with Kenward-Roger first order approximations with Kacker-Harville corrections. Where treatment is a significant effect, Tukey HSD post-hoc tests were conducted to identify how treatments differ. Additionally, to visualize changes in the relative abundances of core taxa over time (continuous), we fit Kernel Smoothing (Loess; [73]) lines using linear local fits, tri-cube weighing, and four iterations to derive best fit lines. Further, to determine if our obtained core taxa putatively inhibit, facilitate, or have no effect on Bd, we compared our core OTUs to sequences obtained from the functional analysis by Woodhams et al. [15] using BLASTn and identified our OTUs as a match if Query Coverage = 100% AND Identity ≥ 99%.
Btk interactions with other taxa
To identify which obtained OTU matched the Btk strain used in our inocuations, we examined representative sequences of each demarcated OTU and compared them (BLASTn) to the Sanger sequence obtained of the isolated Btk used for exposure experiments. OTU72 (hereafter referred to as Btk 72) was identical (100%) to the Btk strain. To investigate if Btk 72 co-associates with (positive correlation) or is mutually exclusive of (negative correlation) core OTUs, Kendall Tau correlations were determined between Btk 72 and all core OTU sequence abundances and alpha levels for significance were adjusted based on Šidàk corrections (adj. α=0.002) for multiple comparison corrections. All statistics were conducted using a combination of R, JMP Pro (v.14), and mothur.