A novel approach to wildlife transcriptomics provides evidence of disease‐mediated differential expression and changes to the microbiome of amphibian populations

Ranaviruses are responsible for a lethal, emerging infectious disease in amphibians and threaten their populations throughout the world. Despite this, little is known about how amphibian populations respond to ranaviral infection. In the United Kingdom, ranaviruses impact the common frog (Rana temporaria). Extensive public engagement in the study of ranaviruses in the UK has led to the formation of a unique system of field sites containing frog populations of known ranaviral disease history. Within this unique natural field system, we used RNA sequencing (RNA‐Seq) to compare the gene expression profiles of R. temporaria populations with a history of ranaviral disease and those without. We have applied a RNA read‐filtering protocol that incorporates Bloom filters, previously used in clinical settings, to limit the potential for contamination that comes with the use of RNA‐Seq in nonlaboratory systems. We have identified a suite of 407 transcripts that are differentially expressed between populations of different ranaviral disease history. This suite contains genes with functions related to immunity, development, protein transport and olfactory reception among others. A large proportion of potential noncoding RNA transcripts present in our differentially expressed set provide first evidence of a possible role for long noncoding RNA (lncRNA) in amphibian response to viruses. Our read‐filtering approach also removed significantly more bacterial reads from libraries generated from positive disease history populations. Subsequent analysis revealed these bacterial read sets to represent distinct communities of bacterial species, which is suggestive of an interaction between ranavirus and the host microbiome in the wild.

Ranaviruses are large double-stranded DNA viruses that are capable of causing disease in all classes of ectothermic vertebrates (Cunningham et al., 1996;Marschang, 2011;Whittington, Becker, & Dennis, 2010). They are globally distributed and known to be present on all continents except for Antarctica (Duffus et al., 2015).
Ranaviruses have been identified as the cause of amphibian mass mortality events (Cunningham et al., 1996) and have been implicated in the long-term declines and extinction of amphibian populations (Earl & Gray, 2014;Teacher, Cunningham, & Garner, 2010). The susceptibility of amphibians to ranaviruses varies between host species, and a number of ecological risk factors have been identified (Hoverman, Gray, Haislip, & Miller, 2011). There is evidence that the level of genetic variation within a population correlates negatively with susceptibility to ranaviral infection (Gantress, Maniero, Cohen, & Robert, 2003;Pearman & Garner, 2005) and a history of host coevolution with local ranaviral strains also reduces susceptibility within a species (Ridenhour & Storfer, 2008;Schock, Bollinger, & Collins, 2009). This is indicative of an adaptive genetic element to ranaviral immunity, an aspect of host-ranavirus interactions that has received little attention to date.
Most of what is known about the amphibian immune response to ranaviruses comes from an experimental system built around the primitive anuran Xenopus laevis infected with the type species of Ranavirus, Frog Virus 3 (FV3; Hyatt et al., 2000). This system has yielded valuable insights into the pathogenicity and immune evasion strategies of ranaviruses (Grayfer, Andino, Chen, Chinchar, & Robert, 2012) as well as the immune response of X. laevis to FV3 infection (De Jes us Andino, Chen, Li, Grayfer, & Robert, 2012;Maniero, Morales, Gantress, & Robert, 2006;Morales et al., 2010;Robert & Ohta, 2009). Although these findings are fundamental and vital to our understanding of how amphibians may combat ranaviral infection, the diversity of amphibian taxa means there is a clear need to extend this research effort to other amphibian species.
The application of RNA-Seq-based transcriptomics to disease research has revolutionized our ability to understand the processes of and response to diseases, even in nonmodel organisms (Costa, Aprile, Esposito, & Ciccodicola, 2012;Wang, Gerstein, & Snyder, 2009). However, the lack of well-controlled experimental field systems, coupled with the high likelihood of simultaneously sequencing environmental contaminants, has meant that the use of RNA-Seq techniques has remained largely rooted in the laboratory.
Across all taxa and research disciplines, field-based transcriptomic studies are rare (see review by Alvarez, Schrey, & Richards, 2015) and to date all examples of RNA-Seq-based transcriptomics in amphibian disease research have been conducted within laboratorybased study systems (Ellison et al., 2015;Price et al., 2015). As extrapolation of findings from laboratory studies to wild systems is often tenuous, to fully appreciate the true impacts of ranaviral infection we need to understand its impact on the transcript-level response of wild animals (Alvarez et al., 2015).
Frog Virus 3 (FV3) was identified as the causative pathogen behind unusual mortality events seen in UK populations of European common frogs (Rana temporaria) during the late 1980s and early 1990s (Cunningham et al., 1996). The subsequent formation of the Frog Mortality Project (FMP) led to unprecedented public engagement in the study of a wildlife disease (see Price, Garner, Cunningham, Langton, & Nichols, 2016) and the formation of a unique database of frog populations known to have a positive history of ranaviral disease. This database, coupled with a complementary set of populations known never to have experienced a ranavirus-related die-off, has provided a novel framework within which to study the population-level impacts and spread of ranavirus infection in the UK (North, Hodgson, Price, & Griffiths, 2015;Price et al., 2016;Teacher, Garner, & Nichols, 2009a,b;Teacher et al., 2010).
In this study, we combine this controlled field system with a novel approach to the filtering of RNA-Seq read data collected from wild animals to investigate the impact of ranaviral disease history on the transcriptome of wild populations of R. temporaria in the UK.

| Ethics statement
This project was approved by the ethics boards of both the University of Exeter and Zoological Society of London and conducted under the project license 80/2466 issued by the UK Home Office.
All field sampling was conducted under the personal Home Office license 30/10730 issued to Lewis Campbell.

| Field site selection and sample collection
Potential field sites were selected from a combined database consisting of the FMP database plus a complimentary database of ranavirus na€ ıve field sites (see Teacher et al., 2010 for more detailed selection criteria). Based on information provided by property owners, the field sites in this database are known to have been either infected by or free from ranaviral disease for at least 20 years. A total of 13 frog populations (seven positive history and six disease free) were recruited.
During the 2014 spring breeding season (February-April), field site owners were asked to capture pairs of adult R. temporaria during amplexus (mating embrace) using nylon hand nets. Captured pairs were placed into individual plastic holding tanks (80 9 50 9 50 cm).
All captured frogs were breeding adults, and none exhibited signs of clinical ranavirus infection (see Cunningham et al., 1996) at the time of sampling. Licensed personnel visited each field site and performed tissue sampling in situ within 24 hr of frog capture. Frogs were rinsed once by immersion into sterile water prior to sampling. The distal portion of the first (inside) digit of a hind limb of each captured frog was clipped using surgical scissors following the application of a topical disinfectant that contained an analgesic (Bactine; WellSpring Pharmaceutical, Florida, USA). Toe clips were immediately placed into individual 1.5 ml microcentrifuge tubes containing 1 ml of RNAlater stabilizing agent (Sigma Aldrich, Missouri, USA). Following sample collection, all animals were released at point of capture.
The number of animals sampled at each site ranged between 2 and 18 and was dependent on the success of field site owners at capturing amplexing pairs.

| RNA extraction and sample pooling
Whole amputated toes (bones and soft tissue) were agitated using a Qiagen High Frequency Tissue Lyser2 (Qiagen, Hilden, Germany) at a frequency of 2,000 Hz for duration of 4 min with lysis buffer and stainless steel lysing beads. RNA was then extracted using NucleoSpin RNA isolation kits from Macherey-Nagel (Duren, Germany) following the manufacturer's instructions.
Approximate concentration of all RNA extracts was assessed using a NanoDrop (NanoDrop, North Carolina, USA). RNA quality of the extracts from the three males and three females with the highest RNA extraction yields from each population was appraised using a BioAnalyser (Agilent Technologies, California, USA). The three positive history and three disease-free field sites with the largest number of high-quality RNA extractions were selected for sequencing ( Figure 1). Due to budget restrictions, samples with RNA integrity scores ≥8 were chosen for pooling. Three males and three females from each of the six sites were pooled, creating six site-specific sample pools of six individuals. To prevent the over-representation of any individuals within a pool, RNA concentrations of each individual within a given pool were normalized to the concentration level of the lowest individual. RNA concentrations were not normalized between pools.

| Sequencing
Library preparation and strand-specific sequencing were performed by the Biomedical Informatics Hub at the University of Exeter. The amount of ribosomal RNA present within each library was minimized using the poly A capture method. Reverse transcription to cDNA was performed using the Superscript II reverse transcription kit (Invitrogen, Massachusetts, USA) and a combination of random hexamer and oligo dT primers. Each library was indexed with a different sequencing barcode, and all six were run simultaneously over the same three lanes of an Illumina Hiseq flow cell on an Illumina Hiseq 2000 machine using v3 chemistry, generating 100 base pair pairedend reads.

| Read quality control
Adapter sequence and low-quality bases were trimmed from the ends of reads in the raw read data sets of each of the six populations independently, using the fastq-mcf utility from the ea-utils package (Aronesty, 2013). Short reads and those containing nonassigned bases (Ns) were also removed from read sets at this stage.
The following quality parameters were used: quality threshold of 20; minimum remaining sequence length of 35 bases; minimum identity between adapter sequence and sequence to be trimmed of 85%; no Ns permitted, and a minimum trim length of 3 bases. The read sets were then evaluated using FastQC (Andrews 2010).

| Environmental contaminant screening
To remove any reads that may belong to potential environmental contaminants, we performed a novel and systematic filtering step using the software package BIOBLOOM TOOLS (version 2.0.13b. ;Chu et al., 2014). BioBloom tools use Bloom filters, which are a memory efficient, probabilistic data structure that tests set membership, to screen data sets (Bloom, 1970). We produced Bloom filters using comprehensive sets of publicly available genomes for a suite of contaminant groups that may be present in urban freshwater environments. These were (in alphabetical order) Algae, Archaea, Bacteria, Fungi, Plants and Viruses. Genomes were downloaded from ENSEMBL (Yates et al., 2016). The program ntCard (Mohamadi, Khan, & Birol, 2016) was used to estimate the number of distinct kmers (subsequences of length k) for each set of contaminant genomes and estimate the number of elements to be inserted into each  Bloom filter (Table S1). All Bloom filters were created with a target false-positive rate (FPR) of 2%. Read sets for each site were run independently against each filter, and any reads that hit a contaminant filter were removed (threshold read categorization score of 0.05). To retain as many potential frog reads as possible, we performed a conservative step using a further Bloom filter constructed from the genomes of three frog species: Lithobates (Rana) catesbeianus (Hammond et al., 2017); Nanoranaparkeri (Sun et al., 2015) and Xenopustropicalis (Hellsten et al., 2010). Any reads that hit against contaminant filters but which also matched this conservation filter were returned to their respective read sets.

| Exploratory investigation of bacterial communities
We used Kraken (Wood & Salzberg, 2014 Jost (2006). Effective number of species provides a more interpretable measure of differences in community diversities than standard entropies (Shannon, Simpson, etc.) and can be calculated from any of them (Jost, 2006).
Its use circumvents the largely subjective choice between more commonly used indices of diversity, allowing for easier comparison and repeatability of diversity estimates between studies (Jost, 2006). We compared these values between disease status groups using a linear model fitted in R. The effective number of species value of each population was fitted as the response variable and disease status of population as the sole predictor. We compared the identified bacterial community composition using a combination of multidimensional scaling (MDS) ordination plots in PHYLOSEQ. We first took observed abundance of individual species into account by ordinating our populations using Bray-Curtis dissimilarity distance measures based on count data of each identified species. Second, we negate the impact of the observed abundance of individual bacterial species on the distances between our sampled populations using the Jaccard index computed using binary presence/absence data of each identified species.
We used DESEQ2 (version 1.10.1, alpha = 0.05, Love, Huber, & Anders, 2014), which corrects for multiple hypothesis testing and discrepancy in sampling effort (read set size), to identify bacterial species that were differentially abundant between positive history and disease-free populations.

| Transcriptome assembly
Our transcriptome was assembled using a multistep workflow; first, the transcriptome assembler TRANS-ABYSS (version 1.5.5; Robertson et al., 2010) was used to produce multiple assemblies for each of our six libraries. It has previously been noted that an assembly pipeline incorporating multiple k-mer lengths accounts for differences in sequence coverage due to variation in transcript abundance within a sample (Robertson et al., 2010). Our strand-specific assemblies were therefore produced over a range of k-mer lengths (k = 24, 32, 40, 48, 56, 64) using the following additional parameters: minimum base quality at read ends = 10, minimum base quality throughout a read = 10, minimum output sequence length of 500 bases and indel size tolerance = 2.
Second, each field site's varying k-mer length assemblies were merged independently using the merge function of trans-ABySS (using the same parameters as initial construction, and percentage identity between sequences to be amalgamated = 95%) to create a total of six assemblies (one per sampled population).
Lastly, a pan-sample transcriptome (representative of all six sample pools) was produced using CD-HIT (strand-specific alignments, sequence identity of contigs to be amalgamated of 80%; Li & Godzik, 2006) to collapse individual library assemblies down to a common reference transcriptome. CD-HIT clusters sequences based on similarity and retains the longest member of each cluster. We assessed the completeness of our final assembly using CEGMA, which reports a proxy for assembly completeness based on the presence of a set of highly conserved "core eukaryotic genes" (CEGs, default parameters; version 2.5, Parra, Bradnam, Ning, Keane, & Korf, 2009).

| Functional annotation
The BLASTX (nucleotide to protein alignment) program of the BLAST+ software package (-strand "plus", Camacho et al., 2009) was used to align our assembled transcripts to a database of reference vertebrate proteins. Our reference database was constructed by compiling the UNIPROT (Bateman et al., 2015) submission databases for vertebrates, human, rodent and mammalian proteins as these gene annotation repositories are maintained independently. Transcripts were annotated as their best alignments based on the alignment criteria of ≥25% identity and ≥25% coverage.
To test the validity of our assembly, transcripts that did not hit against any annotations from the UNIPROT database were then We then used the TransDecoder program (in strand-specific mode with default parameters; Haas et al., 2014) to search for putative protein-coding regions (open reading frame; ORF) within each of our entire set of assembled transcripts.

| Differential expression, permutation and GO term analysis
The read set for each library was individually aligned to our final pan-sample transcriptome assembly using BWA-MEM (version 0.7.12-r1039, Li, 2013 (Table S2). The one exception was for reads identified as having a bacterial origin. Read sets from disease-free populations contained considerably lower percentages of putative bacterial reads (Table 1).
This discrepancy was shown to be statistically significant based on comparison of bacterial read counts per disease history group (  T A B L E 1 Library preparation, sequencing and bacterial read-filtering/classification results. RNA libraries consisted of pools of six individuals (three males and three females) from a given population. The concentration of each individual within a pool was equalized to the concentration of the lowest individual within that pool. The concentration between RNA pools was not equalized. FV3 = the number of reads per population that were mapped to the genome of FV3. Reads generated = the total number of reads generated for that library. Putative bacterial reads = the number of reads that hit the bacterial Bloom filter but not the conservative frog Bloom filter. % total reads = percentage of the total reads generated identified as potentially bacterial in origin. % reads classified = the number of putative bacterial reads that were then classified to a species ID by Kraken. Effective number of species = the exponent of Shannon diversity index      (Table 3). In addition to those annotated to proteins, a further 255 of our differentially expressed transcripts were present in previously assembled frog transcriptomes from other studies. Only 12 of these 255 transcripts were predicted to contain complete ORFs.
No GO categories were significantly enriched in either our up-or down-regulated differentially expressed sets when p-values were adjusted for multiple testing (Benjamini Hochberg correction for false discovery rate).

| Differential expression
Although there was no obvious clustering of expression profiles between our six sampled populations across the whole of our transcriptome, we identified a relatively small suite of 407 transcripts, both up-and down-regulated, that showed strong clustering within disease history classes. Below we will discuss differentially expressed, annotated genes that are potentially most interesting given what is currently known about ranavirus immunology and pathology. This includes increased production of specific antigen receptors (Grayfer et al., 2014) and the production of immunoglobulins that protect against reinfection (Maniero et al., 2006). Two genes that can be linked to such responses also appear in our up-regulated gene list: immunoglobulin lambda variable 3-1, which has the second highest fold change of any up-regulated gene, and endogenous retrovirus group K member 18 (ERVK-18).

| Up-regulated genes
Ranaviruses are known to cause damage to the epithelial mem-

| Down-regulated genes
The largest fold change seen in the down-regulated gene set is docking protein 4 (DOK4). DOK4 is involved in the regulation of T-cell-mediated immune response. This is potentially surprising as the importance of T-cells in clearing FV3 infection in X. laevis has been noted (Robert et al., 2005). This reduction in transcript levels of DOK4 could suggest reduced T-cell-mediated immune capability in R. temporaria populations with a positive history of ranavirus infection. However, observations in X. laevis show that a reduction in T-cell response can also represent increased immune competency and specificity of response (Morales & Robert, 2007).
Expression of an antimicrobial peptide (AMP), Preprofallaxidin-7, was also down-regulated in positive disease history populations. This is perhaps unexpected as the importance of amphibian antimicrobial peptides in the deactivation of ranaviruses has been demonstrated (Chinchar, Wang, Murti, Carey, & Rollins-Smith, 2001). However, the production of AMPs by amphibians is known to contribute to the regulation of commensal microbial communities present on their skin (Walke et al., 2014). This change in expression of an AMP may therefore be contributory to the potential differences in the skin microbial communities of frogs based of ranaviral disease history that we present herein.
Little is known about how ranaviruses impact the development of the nervous system, but it is known that the brain is an important The importance of each of the genes discussed above in relation to ranaviral disease history remains to be verified. However, given that the immunological response of R. temporaria to ranaviruses is very poorly characterized, the genes listed herein, along with others presented in Table 3, are good candidates for further, targeted investigation.  (Cunningham et al., 2007;Teacher et al., 2010).
The RNA pooling approach that we used for this experiment reduced the amount of interindividual variation in gene expression that could be detected, resulting in a bias towards the detection of only those transcripts that are most significantly differentially expressed between the pools (Alvarez et al., 2015). This may explain the relatively low number of differentially expressed transcripts when compared to disease-challenged amphibians of other species and of R. temporaria itself to other pathogens (Cotter, Storfer, Page, Beachy, & Voss, 2008;Ellison et al., 2015;Price et al., 2015).
However, R. temporaria experimentally exposed to ranavirus by Price et al. (2015) showed almost no response at the transcript level suggesting that a potentially weak transcript response to even acute ranaviral infection may be a true reflection of the immune capabilities of R. temporaria and thus contributory to their high susceptibility to disease. However, due to limitations of the experimental design, it should be noted that it is possible that the fold changes in baseline expression of control vs. exposed animals observed by Price et al. (2015) were not large enough to reveal differential gene expression.
Additionally, transcript response in amphibian immunity is known to be tissue specific (e.g., Ellison et al., 2015), so the possibility that both the present study and previous studies have missed transcriptlevel responses to ranaviral infection in other tissues or organs can also not be dismissed.
The diversity of function in our differentially expressed, annotated gene set and our relatively low numbers of functionally annotated transcripts (compared to the number of contigs in our assembled transcriptome) is potentially responsible for the lack of significant enrichment of any one functional pathway. However, functional enrichment was recorded in response to ranaviral infection in the Mexican axolotl (Ambystomamexicanum; Cotter et al., 2008); this may be explained by the use of more target-specific microarray techniques, yielding differentially expressed genes from fewer functional pathways, all of which were functionally annotated.

| Suggestion of a role for long noncoding RNA
Our annotation pipeline was only able to identify those protein-coding genes that have been discovered and functionally evaluated in other vertebrates. However, aside from our 42 annotated proteincoding genes an additional 255 of our 407 differentially expressed transcripts were found to be present in previously constructed de CAMPBELL ET AL.
| 1421 novo amphibian transcriptomes. Only 12 of these 255 transcripts were predicted to have protein-coding potential. The lack of ORF and length (>200 bp) of these transcripts present the possibility that some of the remaining 243 may be long noncoding RNAs (lncRNA).
It has long been acknowledged that noncoding RNA has numerous important regulatory roles in vertebrates (see Pang et al., 2005) and recent work has been shedding light on their importance to vertebrate immune surveillance and response (Carpenter et al., 2013;Guttman et al., 2009). Differential expression of noncoding RNAs has been reported in response to viral infections in several species (Chicken, Ahanda et al., 2009;Mice, Peng, Gralinski, Armour, & Ferris, 2010;Humans, Yin et al., 2013), but their expression in amphibian disease responses is unevaluated. This result suggests that lncRNA may form a large proportion of the differentially expressed transcripts between populations of varying ranaviral disease history.

| Novel approach to wildlife transcriptomics
Harnessing the potential of RNA-Seq to examine the mRNA transcript profiles of wild animals requires that contamination be minimized. However, traditional approaches to the removal of contaminant read data, consisting of limited alignment protocols against the genomes of humans and Escherichia coli bacteria (e.g., Yang, Qi, Bi, Fu, & Access, 2012), rarely encompass the possible range of contaminants of a wild-derived RNA sample. Our adoption of RNA-Seq read-filtering techniques previously utilized in clinical genomics studies (e.g., Burk et al., 2017;Jeraldo et al., 2015Jeraldo et al., , 2016Strong et al., 2014;Zheng et al., 2016) has allowed us to overcome the challenge of contamination in wild-collected samples, facilitating the identification and separation of potential contaminant reads from those of the organism of interest.
In the present study, we found that read sets from populations with a positive history of ranaviral infection contained significantly more bacterial reads than those from disease-free populations. The use of only one sterile wash prior to sampling each animal means that the detection of transient microbes cannot be entirely excluded, although the use of skin tissue should mean that most of the RNA extracted from each sample is directly associated with the host or their commensal microbiota. The detectable bacterial species diversity in samples from positive disease history populations was significantly lower than disease-free populations, with commensal communities of frogs from positive disease history populations being dominated by a single species, Bacillus subtilis. The bacterial community composition between the populations of different disease history was also shown to be distinct (Figure 3). When distances between each population's microbial community structure were computed based on presence/absence data, rather than species abundance counts, we found high similarity between disease-free populations but more divergence between the microbial communities of positive disease history populations (Figure 3). This suggests that between disease-free and positive disease history populations, both the species assemblages and the abundance of represented bacterial species vary but that similarities between the microbial communities from positive disease history samples are primarily driven by the sheer abundance of the most dominant detectable taxa (i.e., Bacillus subtilis).
In total, five bacterial OTUs were identified as differentially abundant between disease history groups. Bacillus subtilis was found to be two orders of magnitude more abundant in positive history vs. disease-free populations. A well-known probiotic, B. subtilis has been shown to increase immune function of fish and reduce the impact of infectious diseases in aquaculture (Aly, Abdel-Galil Ahmed, Abdel-Aziz Ghareeb, & Mohamed, 2008;Liu, Chiu, Wang, & Cheng, 2012;Ran et al., 2012). Although B. subtilis is present in many fish care products, no such products have been used in the ponds sampled in this study. This presents the possibility that the huge overabundance (compared to disease-free populations) of this potentially beneficial bacterium could represent an adaptive change in host microbiome that may infer greater resistance to ranaviral disease.
However, ranaviruses have previously been shown to coinfect with several other micro-organisms in both wild (Landsberg et al., 2013;Warne, La Bumbard, La Grange, Vredenburg, & Catenazzi, 2016;Whitfield et al., 2013) and captive amphibian populations (Miller et al., 2008). An alternative explanation is therefore that these relatively more abundant bacterial species are exploiting the potentially immune-compromised hosts from populations with a positive history of ranaviral disease to establish themselves at higher levels on amphibian skin.
A growing body of work has examined links between the amphibian skin microbiome and Batrachochytrium dendrobatidis (Bd), which are now known to be intimately linked. Infection with Bd has been shown to alter the structure of the commensal bacterial communities on amphibian skin (Jani & Briggs, 2014;Woodhams et al., 2014), and several species of bacteria isolated from amphibian skin have been shown to inhibit Bd in the laboratory (Antwis, Preziosi, Harrison, & Garner, 2015;Park, Park, Collingwood, St-hilaire, & Sheridan, 2014). In fact, the structure of amphibian microbial communities has been shown to be an effective predictor of disease risk associated with Bd (Becker et al., 2015;Jani, Knapp, & Briggs, 2017;Woodhams et al., 2014), and it is postulated that the relative immunity of R. temporaria to Bd infection may be due to the microbiota present on their skin (Woodhams et al., 2014). Laboratory studies have also demonstrated that a more diverse skin microbiome increases survival of captive-reared R. temporaria experimentally infected with ranavirus (Harrison, Price, Hopkins, Leung, Sergeant, & Garner, 2017). Although incidental, our findings suggest that a relationship may also exist between the skin microbial communities of R.
temporaria and ranaviruses in the wild. However, exactly how a host's microbiome is impacted by a pathogen depends on its preinfection state (Becker et al., 2015;Harrison et al., 2017;Jani & Briggs, 2014;Longo & Zamudio, 2017). With no knowledge of the structure of the amphibian skin microbiomes at positive disease history populations prior to disease outbreak, it is not possible to know whether potential changes in the microbiome structure are a cause or consequence of disease history.
It is also important to note that these findings are based on reads generated as a by-product of shotgun RNA sequencing. This means that the read counts mapping to each species may not represent those species which are most abundant in the skin microbial communities in terms of biomass but those species which are expressing some portion of their genome at the highest levels. For example, it may not be that B. subtilis is necessarily two orders of magnitude more abundant in positive disease history populations but that the B. subtilis that are present are expressing a portion of their genome at a level two orders of magnitude higher than in diseasefree populations. Given the perfect correlation between positive ranaviral disease history and high abundance of B. subtilis, we have no way of concluding whether the transcript changes discussed previously are due to a positive history of acute ranaviral infection or to the abundance or transcriptional activity of B. subtilis, and further experimental work is warranted to this end.
Elucidating the nature and direction of this potential relationship along with its underlying mechanisms will require the application of several complementary approaches but is no doubt worthy of the required research effort. As our finding relies on tissue samples of a small number of individuals, a more comprehensive study of the interaction between ranaviruses and the structure of amphibian microbiomes in the wild would prove valuable. The use of a more traditional method of community classification such as 16S rRNA amplicon sequencing would demonstrate whether the over abundance of B.
subtilis observed here is due to increased biomass of that species or potentially elevated levels of gene expression. Additionally, investigations of the interaction between ranaviruses and B. subtilis both in vivo and in vitro might prove important as understanding potentially inhibitory interactions between commensal bacteria and ranaviruses may lead to innovative use of probiotics in disease mitigation strategies, similar to those already observed in response to Bd (e.g., Woodhams et al., 2014) or in aquaculture (e.g., Aly et al., 2008).
Finally, although none of the animals sampled for this study were symptomatic of ranaviral infection, we found that sample pools from both positive disease history and disease-free populations contained similar numbers of reads that mapped to the genome of FV3. The population monitoring of our field site owners has previously been verified by diagnostic sampling of R. temporaria collected from most of the populations which were sampled for this study (Teacher et al., 2010); therefore, we have no reason to doubt the disease history classification of each population. Due to the possible detection of transient microbial diversity from the environment, this result may suggest that ranaviruses are present even at those populations where acute disease has not been observed. This may mean that ranaviruses are more ubiquitous than currently thought and that other factors dictate whether disease occurs within a given population. The use of new environmental DNA sampling techniques (i.e., Hall, Crespi, Goldberg, & Brunner, 2016) will make it possible to explore how universally distributed ranavirus virions are compared to incidences of known disease outbreak and to investigate potential relationships between species diversity within an environment (microbial and otherwise) and the occurrence of ranaviral disease.

| CONCLUSION AN D FUTUR E DIRE CTION S
Our results suggest that a disease history of ranavirus may alter the gene expression profile of wild R. temporaria populations but that such changes may not be associated with heightened ability to cope with the presence of a ranavirus. We provide a set of 42 candidate genes which form a good basis from which to begin a more extensive study of the response of R. temporaria to ranaviral infection.
We also found many differentially expressed transcripts which could not be currently annotated as vertebrate protein-coding genes and which appear to possess no ORF. This fact, along with their length, suggests that they could represent long noncoding RNAs.
There is growing recognition of the role that lncRNA plays in the response of various vertebrates to viral infections, so a more concerted effort to evaluate the expression of lncRNAs in amphibian response to disease is warranted.
The incorporation of a novel read-filtering technique has provided us with evidence of a relationship between ranaviruses and the bacterial communities inhabiting amphibian skin in the wild. The relationship between pathogens and amphibian disease is a burgeoning area of research in relation to fungal pathogens. Therefore, further targeted investigation into how ranaviruses and amphibian microbiomes interact is timely, and the implications of such research for amphibian conservation are potentially profound.