Experimental infections and tissue harvesting

Four captive-bred adult Atelopus zeteki (41 months postmetamorphosis) and Hylomantis lemur (12 months postmetamorphosis) were used for experimental infections. Prior to experiments, frogs of each species were cohoused, had never been exposed to or tested positive for Bd, and never been treated with antifungals. All individuals appeared healthy and tested negative for Bd at the start of the experiment. Animals were housed individually in plastic shoeboxes with a water dish and paper towel, in a laboratory maintained at 21–22° with a 12:12 hr light:dark photoperiod. All housing materials were replaced every 7 d, water dishes changed every 3 d, and animals were fed crickets and fruit flies sprinkled with Herptivite (Rep-Cal Research Labs) every 3 d.

For experimental infections, we used Bd strain JEL423, cultured (Longcore et al. 1999) and cryopreserved from a field-infected H. lemur, during the 2004 epidemic outbreak at El Copé, Panama (Lips et al. 2006). We inoculated each experimental animal for 10 hr in a 50 ml bath of Bd zoospores (2 × 10⁴ zoopsores/ml). To provide culture-controls of the same passage number (total five passages), the same zoospore stock was used to inoculate four agar plates (1 ml at 1 × 10⁵ zoospores/ml) and cultured at 21–22° for 7 d. This ensured all frogs and plates were inoculated with an equivalent number zoospores (1 million). After culturing, Bd colonies were scraped from plates using sterile instruments, twice syringe-filtered (Life Sciences Acrofilter 10 µm membrane) to obtain sporangia, pelleted by centrifugation, snap-frozen using liquid nitrogen, and stored at −80° until RNA extraction.

Throughout the infection trial, we used a fresh pair of latex powder-free gloves when handling each individual frog. All individuals were assayed for Bd infection status and load at d 8 postinoculation, and thereafter every 3–4 d. We swabbed the abdomen, drink patch, hands, and feet, five times each with a sterile cotton tipped swab, which was then stored in capped tubes containing 30 μl of 70% ethanol (Hyatt et al. 2007). We tested swabs for Bd using PrepMan Ultra extractions and singlicate Taqman quantitative polymerase chain reaction (qPCR) (Boyle et al. 2004). We ran each plate with JEL423 standards of 0.1, 1, 10, 100, and 1000 zoospore genomic equivalents to determine Bd presence and infection intensity. We categorized individuals as Bd-positive when qPCR results showed an infection load ≥1 Bd zoospore genomic equivalent (Kriger et al. 2006). Individuals were monitored daily for clinical signs of chytridiomycosis, and killed by application of 10% Benzocaine to the venter once they lost righting abilities.

Immediately after euthanasia, using sterile instruments, skin from the ventral thigh was harvested from each individual. Skin pieces were placed in Tissue-Tek cryomolds (Sakura Finetek), covered in Tissue-Tek Optimum Cutting Medium (OCT, Sakura Finetek), and immediately snap-frozen in an isopentane/dry ice bath. Molds were stored at −80° until cryosectioning. As time to death due to Benzocaine treatment was rapid (<2 min), and skin tissue was harvested away from site of application, we minimized the potential for the euthanasia method to affect Bd gene expression, and therefore confound comparison with cultured Bd.

The number of Bd zoospores that infect individuals during inoculations is highly variable, causing natural differences in infection intensities and time to death among individuals and species (Carey et al. 2006). In this study, clinical signs of chytridiomycosis developed 15–18 d postinoculation in A. zeteki and 35–51 d in H. lemur. We cannot ascertain that all frogs were at the exact same stages of disease progression at the time of sampling. However, the infected frogs and culture plates used in this study had in common that all were at a “mature” infection stage, colonized with mature zoosporangia and shedding zoospores. At the point of tissue sampling, all frogs carried their highest individual Bd load, and thus all could be considered still within the disease progression phase.

RNA extraction and sequencing

To minimize RNA degradation, all sectioning, staining, and laser-capture microdissection (LCM) steps were performed as quickly as possible. Prior to cryosectioning, PEN membrane slides (Zeiss) were sterilized by dipping in RNase Zap (Ambion) for 30 sec, rinsed in nuclease-free water, and dried for 30 min under UV lights. Tissue molds were allowed to equilibrate to −20° within the cryostat for ∼15 min prior to sectioning. Approximately five 8 µm skin sections were mounted on each prechilled slide. Slides were fixed in ice-cold 100% ethanol for 5 min, placed individually in sterile 50 ml falcon tubes containing desiccant, and stored at −80° (for a maximum of 4 wk) until LCM. Cryostat blades were changed between each sample.

Immediately before LCM, slides were removed from −80° and placed in ice-cold nuclease-free water for 30 sec to remove OCT. To aid rapid identification of highly infected skin areas, slides were placed for 30 sec in Calcofluor white (1%, Sigma-Aldrich)—a chitin-binding fluorescent stain. Slides were then rinsed in ice-cold nuclease-free water, and dehydrated for 30 sec in 70%, 75%, 95%, and 1 min in 100% ice-cold ethanol. Stained and dehydrated slides were immediately microdissected using a Zeiss PALM MicroBeam laser capture system with a 430 nm filter. For each slide, we dissected as many infected areas of the stratum corneum and the immediately underlying epidermal cells as possible within 2 hr. Previous studies of RNA integrity during LCM indicate that dehydrated slides can be worked on for up to 4 hr without significant loss of RNA quantity or quality (Grover et al. 2012). Microdissected skin areas were collected in 500 µl adhesive-capped tubes (Zeiss). After collection, tissue was lysed by addition of 120 µl of Qiagen RLT buffer with β-mercaptoethanol (10 µl in 1 ml buffer), vortexed for 30 sec, and incubated at room temperature in an upside-down position for 30 min. Tubes were then centrifuged for 5 min, and stored at −80° until RNA extraction.

For each individual, three lysates were pooled, and total RNA extracted using RNeasy micro kits (Qiagen) and eluted in 15 µl nuclease-free water. Due to the extremely low quantities of RNA (typically undetectable using high-sensitivity RNA Qubit assays; <0.2 pg/µl), whole transcriptome amplification was performed prior to Illumina sequencing library preparation. cDNA reverse transcription and amplification was performed using Clontech SMARTer Ultra Low RNA Amplification kit using 5 µl eluted total RNA, and the minimum number of PCR cycling steps (nine cycles). Amplified cDNA was quantified using Qubit high-sensitivity DNA assays. Sequencing libraries were prepared and uniquely indexed using the Illumina Nextera XT kit with 100–150 pg input cDNA. Libraries were assessed for quality and fragment size range using Bioanalyzer 2100 assays. SMARTseq whole transcriptome amplification is demonstrated to have negligible systematic biases, particularly in multi-cell preparations as used here (Ramsköld et al. 2012). However, identical RNA extraction and mRNA amplification protocols were used for cultured Bd to avoid potential biases between the culture and in vivo Bd transcriptomes.

Libraries were pooled and sequenced using Illumina NextSequation 75 bp runs. Pilot sequencing studies indicated that six frog/Bd samples per lane (∼60 million reads per sample) would provide sufficient sequencing depth to fully capture the transcriptomes of both host and Bd. Cultured Bd libraries were sequenced at a sequencing depth equivalent to 16 samples per lane (20 million reads per sample) because only Bd gene expression would be captured in these samples (and thus would not be dominated by host transcripts).

Gene expression analyses

After Illumina standard quality control filtering (all reads flagged as low quality were removed), read quality for each sample was visualized using FastQC version 0.10.0 (Andrews 2010), checking for per base sequence quality, GC content, and overrepresented sequences (e.g., contamination). We used Trimmomatic version 0.27 (Lohse et al. 2012) to: (1) trim Illumina adapter sequence, (2) trim the 5′ and/or 3′ end of reads where quality score dropped below Q20, (3) trim anywhere within each read where a 5 bp window dropped below Q20, and (4) discard any trimmed reads <36 bp long. This ensured only the highest quality reads were used for subsequent analyses. For each A. zeteki sample, reads were separated into those derived from frog and Bd by first mapping all reads to the published JEL423 transcriptome (Batrachochytrium dendrobatidis Sequencing Project, Broad Institute) using bowtie v2.2.4 with default parameters (Langmead et al. 2009). Unmapped reads were then mapped to a previously assembled A. zeteki transcriptome (Ellison et al. 2014). To further ensure read identity, we removed all reads that mapped to both host and pathogen transcriptomes. Gene expression raw read counts were obtained for both host and pathogen using RSEM v1.2.2 (Li and Dewey 2011). As H. lemur has no available transcriptome, reads were first mapped to the Bd transcriptome, and then unmapped reads were pooled from all samples and de novo assembled into H. lemur transcripts using previously published pipelines utilizing the Trinity assembly package (Haas et al. 2013; Ellison et al. 2014, 2015).

To compare Bd gene expression between culture and the two host species, differential expression analyses were performed using DESeq2, which normalizes raw read counts by library depth (Love et al. 2014). Genes were deemed significantly differentially expressed using a false discovery rate (FDR) corrected p-value of 0.05. Due to the potentially variable degree of transcriptome coverage resulting from the within-host samples, only genes that had detectable expression in at least two samples of all three conditions (culture, A. zeteki host, and H. lemur host) were considered for analysis (8268 of the 8819 JEL423 reference transcripts). Functional annotation [gene ontology (GO) mapping and InterPro scans] of the JEL423 reference transcriptome was performed using Blast2Go (Conesa et al. 2005), using default settings (Supplemental Material, File S4). GO enrichment tests (with FDR correction) were carried out on each group of differentially expressed genes via Blast2Go to detect significantly over-represented biological processes and molecular functions. Similarity of per-sample Bd gene expression profiles was visualized using hierarchical clustering of Poisson distances of gene expression values using DESeq2.

Expressed single nucleotide variant (eSNV) analyses

Bd sequence reads for each treatment (culture, A. zeteki or H. lemur) were pooled, and single nucleotide variants (SNV) were called using samtools mpileup using default settings (Li et al. 2009). RNAseq data do not allow us to discriminate whether SNV allele frequency differences are due to allele-specific expression, or underlying genotype frequency differences in pooled populations (Konczal et al. 2014); therefore, we refer only to “expressed genotypes.” We used the program Popoolation2 (Kofler et al. 2011a) to calculate expressed SNV (eSNV) allele frequencies within each treatment pool. We applied the conservative parameters of minimum mapping quality = 20, minimum count = 6, and minimum coverage = 20 (Kofler et al. 2011b), and also tested a higher minimum allele count of 10 to ensure the results were robust (Figure S1). To further filter out potential sequencing errors, only eSNVs with a minor allele frequency of ≥0.1 were used in subsequent analyses. These criteria allowed investigation of differences in eSNV allele frequencies between treatment populations, but not detection of newly acquired SNVs, as each polymorphic site had to be detected in all populations (i.e., in culture and both host species) to be considered further. Significant differences in eSNV allele frequencies were determined using two methods. First, Fisher’s exact tests (and FDR-corrected p-values) were used to estimate the significance of eSNV allele frequency differences between treatment pools. Second, each sample was treated as a separate, replicate pool of Bd and tested in pairs using the Cochran-Mantel-Haenszel test implemented in Popoolation2. Only biallelic eSNVs found to be significant in both analyses were considered further. Custom Perl scripts were used to determine if SNVs were synonymous or nonsynonymous. Blast2Go enrichment analyses were used to test for GO term enrichment in sets of genes containing eSNVs found to have significant differences in allele frequencies between treatment groups.

Data availability

All sequence data are available at the NCBI Short Read Archive with the accession number SRP064950. File S1 contains details of all differential gene expression tests. File S2 contains all eSNVs with significant differences in frequencies and their locations. File S3 contains details of the eSNVs matched to previous resequencing studies. File S4 contains the GO annotation of the JEL423 reference transcriptome generated by Blast2Go. File S5 contains the custom perl script used to determine synonymous and nonsynonymous eSNVs. Table S1 and Table S2 contain the top 20 differentially expressed Bd genes of all comparisons between culture plates and amphibian hosts. Figure S1 shows sample separation by eSNV frequencies using two different minimum allele count thresholds and when considering only SNVs found in previous resequencing studies. A Venn diagram of the overlap of genes found to have significant differential expression and eSNV frequencies is provided in Figure S2.