A de novo Assembly of the Common Frog (Rana temporaria) Transcriptome and Comparison of Transcription Following Exposure to Ranavirus and Batrachochytrium dendrobatidis

Abstract

Amphibians are experiencing global declines and extinctions, with infectious diseases representing a major factor. In this study we examined the transcriptional response of metamorphic hosts (common frog, Rana temporaria) to the two most important amphibian pathogens: Batrachochytrium dendrobatidis (Bd) and Ranavirus. We found strong up-regulation of a gene involved in the adaptive immune response (AP4S1) at four days post-exposure to both pathogens. We detected a significant transcriptional response to Bd, covering the immune response (innate and adaptive immunity, complement activation, and general inflammatory responses), but relatively little transcriptional response to Ranavirus. This may reflect the higher mortality rates found in wild common frogs infected with Ranavirus as opposed to Bd. These data provide a valuable genomic resource for the amphibians, contribute insight into gene expression changes after pathogen exposure, and suggest potential candidate genes for future host-pathogen research.

About the Authors:

Stephen J. Price

* E-mail: [email protected] (SJP); [email protected] (AGFG)

Affiliations UCL Genetics Institute, University College London, Darwin Building, Gower Street, London, United Kingdom, Institute of Zoology, Zoological Society of London, London, United Kingdom

Trenton W. J. Garner

Affiliation: Institute of Zoology, Zoological Society of London, London, United Kingdom

Francois Balloux

Affiliation: UCL Genetics Institute, University College London, Darwin Building, Gower Street, London, United Kingdom

Chris Ruis

Affiliation: UCL Genetics Institute, University College London, Darwin Building, Gower Street, London, United Kingdom

Konrad H. Paszkiewicz

Affiliation: Wellcome Trust Biomedical Informatics Hub, Biosciences, Geoffrey Pope Building, University of Exeter, Streatham Campus, Exeter, United Kingdom

Karen Moore

Affiliation: Wellcome Trust Biomedical Informatics Hub, Biosciences, Geoffrey Pope Building, University of Exeter, Streatham Campus, Exeter, United Kingdom

Amber G. F. Griffiths

* E-mail: [email protected] (SJP); [email protected] (AGFG)

Current address: FoAM Kernow, Studio E, Jubilee Warehouse, Commercial Road, Penryn TR10 8FG, United Kingdom (Current address)

Affiliation: Environment and Sustainability Institute, University of Exeter, Penryn Campus, Penryn, Cornwall, United Kingdom

Introduction

Amphibians are currently undergoing a mass extinction event [1]. Two key pathogens are known to be contributing to amphibian population declines and species extinctions: the fungus Batrachochytrium dendrobatidis (Bd) which causes chytridiomycosis, and the Iridoviridae genus Ranavirus [2–5]. Bd is a non-hyphal zoosporic fungus which causes mortalities on every continent except Antarctica (http://www.bd-maps.net/) and is thought to have caused multiple species extinctions [6]. Ranaviruses are large double-stranded DNA viruses, capable of crossing poikilothermic class boundaries, and implicated in mass die-off events and population declines [2,5,7]. Ranavirus and Bd are both noted for their virulence across a broad range of hosts but previous research on wild and captive animals points to contrasting levels of pathogenicity in European common frogs (Rana temporaria). Common frogs are highly susceptible to Ranavirus infection in the UK [5,8] and Spain [7] but seem relatively resistant to Bd [9]. It is unknown whether this difference in susceptibility reflects differences in the host’s immune response to each pathogen. De novo RNAseq offers an ideal opportunity to further our understanding of the host response to Bd and Ranavirus, allowing the identification of genes and pathways involved in the response to infection.

Bd is thought to kill its host by disrupting the cutaneous integrity and function of amphibian skin [10]. Transcriptome sequencing has suggested that more resistant hosts may increase expression of genes involved in skin structure [11] and dramatic up and down-regulation of pathways relating to collagen, fibrinogen, elastin and keratin have been reported in the skin of adult amphibians experimentally infected with Bd [12]. Enrichment of inflammatory responses when challenged with Bd may be a general response regardless of susceptibility [11,13]. Increased expression of microbial peptides in response to Bd infection has also been identified in several species (Rana sierrae, Rana muscosa adults [12]; Xenopus tropicalis adults [14]), indicating that innate immunity is a component of host defence. Robust adaptive immunogenetic responses to Bd infection have in general not been observed, but components of innate and adaptive immunity have been shown to operate even in species that are highly susceptible [13]. In addition significant up and down regulation of adaptive immune genes (including the Major Histocompatibility Complex (MHC) class I and II in particular) have been shown in experimentally infected adult ranids [12,14] and a comparison of responses to Bd in hosts of varying susceptibility suggested that the ability to escape immunosuppression by mounting T-cell mediated responses may determine resistance [11]. Our understanding of the response of juveniles to Bd remains limited.

To date, there has been one study of the transcriptional response of amphibians to Ranavirus. Ambystoma mexicanum (axolotl) showed a significant immunological response when experimentally challenged with the Ranavirus Ambystoma tigrinum virus (ATV) [15]. Wild ambystomatid salamanders are highly susceptible to ATV, however experimental animals appeared to mount an immune response to infection within 24 hours of exposure. Using spleen tissues processed through microarrays, the authors demonstrated changes in the expression of innate immunity genes and the transcriptional response increased through time [15]. Experimental infection of adult Xenopus laevis with a strain of Ranavirus (FV3) has demonstrated increased expression of pro-inflammatory cytokines e.g. tumour necrosis factor alpha (TNF-α) and interleukin-1β (IL-1β), indicating that (like for Bd) the innate immune response is activated following infection [16]. Tadpoles show weaker and more delayed up-regulation of these genes [17]. Specific MHC class Ia gene supertypes have been found to be associated with infection status of adult wild common frog (R. temporaria) populations, and diseased populations are characterized by more similar supertype frequencies (lower FST) than infected populations, indicating pathogen-driven selection on the MHC [18]. This implies that an adaptive immune response to Ranavirus occurs in R. temporaria, and may be important for survival after infection. While adult Xenopus laevis are able to clear FV3 infections, tadpoles do not appear to mount an adaptive immune response to FV3, and succumb to infection within a month of inoculation [17]. Again, little is known about the immune response of juveniles. In general, metamorphosis is a critical point in amphibian immunity–the adaptive immune response appears to be limited pre-metamorphosis, the innate immune response is transformed at metamorphosis, and during metamorphosis individuals are thought to experience temporary immunosuppression [19].

In order to better understand the host response to Bd and Ranavirus, in this study we (i) used RNAseq to generate an annotated de-novo transcriptome for the common frog (R. temporaria) (ii) conducted comparative expression profiling of the early responses of metamorphic frogs exposed to Bd or Ranavirus relative to control animals, and (iii) identified candidate genes for future studies on the population-level impacts of these pathogens.

Methods

Experimental treatments

Ranavirus (RUK13 isolate; [20]) was cultured at 24ºC in Fathead Minnow cells (FHM, Epithelial-like cells from the posterior anal tissue, obtained from the European Collection of Cell Cultures catalogue number 88102401) in maintenance media (EMEM + 10% FBS + 1% L-Glutamine) and quantified on 96 well flat-bottomed cell culture plates using the TCID50 method [21]. Bd inoculum (Isolate IA 042, BdGPL; [22]) was prepared by culture in mTGhL medium at 18ºC for four days before counting zoospores using a haemocytometer.

R. temporaria eggs were obtained from a private garden pond in Chessington in the UK with the permission of the landowner [5]. This site has a known history of Ranavirus infection, but an unknown history of Bd infection. The eggs were hatched and reared through metamorphosis under controlled conditions, showing good survival and no signs of infection. Metamorphs (n = 45) were moved into an experimental room (18-21ºC, 33–46% humidity, full spectrum UV light) for acclimatisation one week prior to exposure, and were placed in individual boxes and fed on crickets ad libitum.

This work was carried out under Home Office license (Project Licence numbers PPL 80/2214 and PPL 80/2466) and was approved by the Institute of Zoology Ethics Committee and the University of Exeter Ethical Review Board. In total 45 frogs were used—15 animals for each of three exposure treatments (Ranavirus, Bd, Control). Animals were checked, cleaned and given fresh water frequently whilst being mindful of causing them distress through unnecessary handling. Exposure was performed in individual tubes with 29ml aged tap water with 1ml of inoculum (Ranavirus at 1k TCID50/ml or Bd at 100k active zoospores/ml), or 30ml tap water for the control treatment, for four hours. Following exposure, the frogs were returned to their individual boxes and examined daily for signs of Ranavirus infection (oedema of the eye, skin ulceration, bleeding), which would have served as an endpoint for the experiment but none were observed. Frogs were euthanized according to Schedule 1 to the Animals (Scientific Procedures) Act 1986 four days post-exposure by immersion in fresh 5g/L Tricaine methane sulphonate (MS222, Pharmaq Ltd.) solution neutralized with sodium bicarbonate in accordance with Universities Federation for Animal Welfare guidance to ameliorate suffering. We sampled livers because they are large, easily targeted organs which are important for immunity and are a target organ in ranavirus disease [23]. Livers were immediately dissected, preserved individually in RNAlater solution (Sigma Aldrich), and stored at -80ºC.

Sequencing, de novo assembly and abundance estimation

Total RNA was extracted from 5–20μg of liver tissue using the Qiagen RNeasy Mini Kit (Qiagen, Valencia USA) using the standard protocol. Each treatment consisted of 15 frogs; three pools were sequenced per treatment, with five individuals per pool. We therefore prepared nine samples for sequencing; 3 replicates per treatment, each consisting of pooled RNA extracted from the livers of 5 individuals. RNA concentration was estimated using a NanoDrop (NanoDrop, Wilmington USA), and concentrations of individual sample extractions were equalized within pools, resulting in a total pool volume of 80μl (concentrations varied between pools, Supporting Information S1 table). The sample concentrations were not normalized. DNase treatment, library preparation and 100 base pair paired-end sequencing were performed at the Wellcome Trust Biomedical Informatics Hub using an Illumina HiSeq 2500 with two samples per lane.

Raw reads were processed with the fastq-mcf package (ea-utils; https://code.google.com/p/ea-utils/wiki/FastqMcf; [24]) to trim low quality bases and adaptor sequences from the ends of reads and remove short reads as well as those containing non-assigned bases (Ns). The following settings were used: quality threshold of 20, minimum remaining sequence length of 35, minimum identity between adapter sequence and clipped sequence of 85%, no Ns permitted, and a minimum clip length of 3. The results were then evaluated using FastQC [25]. Reads from all samples were combined and processed through the standard Trinity pipeline (r2013-02-25 release) for de novo RNAseq assembly [26]. Assembly computation requirements were reduced by performing in silico normalization on the reads to reduce the total number of reads and remove errors whilst maintaining transcriptome complexity (maximum coverage for reads = 30; minimum kmer coverage for catalogue construction = 2). Isotigs were assembled along Trinity’s three-stage protocol (Inchworm, Chrysalis, Butterfly) with default settings including a minimum isotig length of 200bp. Transcript abundance estimates for each sample were obtained using RSEM (packaged with Trinity) [27].

Functional annotation and comparative transcription rates

The assembled transcriptome was annotated using Trinotate, a suite of programs for functional annotation of transcriptomes that is suitable for use with non-model organisms. Sequences were processed through a pipeline consisting of a homology search (NCBI-BLAST), protein domain identification (HMMER/PFAM), protein signal prediction (singalP/tmHMM), and comparison to EMBL Uniprot eggNOG and the GO Pathways annotation databases [28–37].

The reference assembly was filtered prior to differential expression analysis so that the assembly used for downstream analyses might better reflect transcribed genes and to reduce the number of comparisons undertaken. Filtering was based on the number of mapped reads. Assembled isotigs were retained only if the FPKM (expected fragments per kilobase of transcript per million fragments sequenced [38]) was greater than or equal to one for all three replicates within one or more treatments (FPKM ≥ 1).

We used the CEGMA pipeline to accurately annotate Core Eukaryotic Genes (CEGs; [39]). CEGMA can be used to assess assembly completeness and the impact of downstream filtering steps via the proportion of core genes present. This proportion is calculated relative to a reference set containing 248 of the most highly conserved CEGs and analysis suggests it is a good metric to assess completeness of draft assemblies [40]. We ran CEGMA with default settings on our reference (unfiltered) and FPKM ≥ 1 filtered assemblies to obtain the number of complete [more than 70% of the protein length aligned to isotig(s)] and partial [alignment length is less than 70% but remains higher than a minimum alignment threshold] core genes that these assemblies contained.

Differential expression analysis was performed with edgeR using a Trinity packaged Perl script. Subsets of differentially expressed transcripts were compiled based on log-fold change of one and Benjamini–Hochberg adjusted p-values [41] to control for false discovery rate [FDR<0.05 to obtain the transcript list and FDR<0.10 for downstream Gene Ontology (GO) term enrichment analysis]. All pairwise comparisons were considered; Bd vs. control, Ranavirus vs. control, and Bd vs. Ranavirus. Differentially expressed transcripts in the Bd vs. Ranavirus comparison were allocated to either Bd or Ranavirus (S1 Text).

Analysis of Gene Ontology term enrichment

Bingo (a Cytoscape plugin; [42]) was used to search for GO term enrichment (identifying which GO terms are over or under-represented). Bingo compares the list of differentially expressed products to a user generated list of all genes/transcripts and is therefore of particular use for research on non-model organisms. EdgeR subsets of differentially expressed transcripts at FDR = 0.05 and FDR = 0.10 were analyzed.

Results

Sequencing, Assembly and annotation

A total of 1.29 x 109 reads were generated across all nine samples with a mean quality score of 35.2. All raw data is available through EMBL-EBI Array Express, accession number E-MTAB-3632. In total, 199,602 isotigs were generated with an N50 isotig length of 1,086 base pairs and 30,931 isotigs longer than 1,000bp. There were 134,080,068 bases contained in all isotigs with a GC content of 44%. Assembly summary statistics are shown in Table 1. Trinotate generated a list of 48,263 annotated isotigs from 32,486 disconnected subgraphs. These isotigs hit 17,631 genes when orthologs were included and 11,851 genes when excluding orthologs. The species with the highest number of hits were humans (5,605 genes), mouse (3,800), Xenopus laevis (1,720), rat (1,382), cow (1,154) and Xenopus tropicalis (1,035). In total, 71–76% of reads mapped back to the reference assembly from each set of sample reads (i.e. as part of RSEM analysis); control 73–76%, Ranavirus 71–74%, Bd 72–74%.

[Figure omitted. See PDF.]

Table 1. Trinity assembly summary statistics.

https://doi.org/10.1371/journal.pone.0130500.t001

The reference assembly was almost complete when compared to CEGMA’s 248 CEG set, giving us confidence that our methods have yielded an assembly that is a good approximation to the actual transcriptome. In total, 240 (of 248; 97%) complete CEGs and an additional six partial CEGs (total = 246 of 248; 99%) were found in our reference assembly. Summary statistics measuring the completeness of the assembly (broken down by KOG group) are included in S2 Table.

Some genes were lost through our filtering operation but a large majority of CEGs remained. The FPKM ≥ 1 filtered assembly contained 215 (87%) complete genes out of the 248 most highly conserved CEGs. A further five partial genes were present giving a total of 220 CEGs (89%; S2 Table).

Differential expression

Transcriptional profiles of replicates were consistent within treatments, clustering together within treatments when expression values were compared for each pair of samples (Fig 1). Expression values were more similar within the Ranavirus and control treatments than in the more divergent Bd treatment. The total number of differentially expressed transcripts that we report was affected by the False Discovery Rate (FDR), as well as our procedure for re-allocating some differentially expressed transcripts from the Bd vs. Ranavirus comparisons to the other comparisons (Table 2; Fig 2). Increasing the FDR increased the number of differentially expressed transcripts in our output. Classifying some of the Bd vs. Ranavirus transcripts to Bd or Ranavirus vs. control (S1 Text) also increased the number of transcripts for each pathogen vs. control comparison but reduced the total number of unique differentially expressed transcripts because, for example, some of these transcripts were due to small but opposite effects of the two pathogens compared to controls.

[Figure omitted. See PDF.]

Fig 1. Comparison of transcriptional profiles across all samples.

Heatmap visualizing the hierarchically clustered Spearman correlation matrix resulting from a comparison of the transcript expression values (TMM-normalized FPKM) for each pair of samples; Bd = Batrachochytrium dendrobatidis, Con = control, Rv = Ranavirus.

https://doi.org/10.1371/journal.pone.0130500.g001

[Figure omitted. See PDF.]

Fig 2. Venn diagrams showing distribution of differentially expressed transcripts across comparisons of treatments.

Bd vs. control, Ranavirus vs. control, and Bd vs. Ranavirus at a) FDR < 0.05 & b) FDR < 0.10.

https://doi.org/10.1371/journal.pone.0130500.g002

[Figure omitted. See PDF.]

Table 2. Summary of the total number of differentially expressed transcripts under alternate FDR regimes: FDR <0.05, FDR<0.10, and differentially expressed transcripts from the Bd-Ranavirus comparison at FDR<0.10 allocated to one of the other comparisons (S1 Text).

https://doi.org/10.1371/journal.pone.0130500.t002

Exposure to Bd elicited far greater transcriptional divergence from controls than did exposure to Ranavirus (Table 2, Fig 3; 120 differentially expressed transcripts for Bd compared to 23 for Ranavirus at FDR<0.05). After allocation of the Bd vs. Ranavirus transcripts to one of the pathogen treatments, there were a total of 360 transcripts (294 up-regulated, 66 down-regulated, 122 annotated) for Bd and 65 (35 up-regulated, 30 down-regulated, 16 annotated) for Ranavirus (Fig 4). The overall tendency for up-regulated expression in Bd-exposed animals was not seen in Ranavirus-exposed animals (Fig 4). Amongst annotated transcripts that were differentially expressed after Bd challenge, interferon-induced proteins figure prominently through Interferon-stimulated 20kDa exonuclease-like 2 and multiple versions of Interferon-induced very large GTPase 1 proteins (Table 3).

[Figure omitted. See PDF.]

Fig 3. Relative expression of differentially expressed (FDR < 0.05) transcripts (rows) across all samples (columns).

Dendrograms show relationships between samples based on expression values (top) and between transcripts based on comparative expression across samples (left). Bd = Batrachochytrium dendrobatidis, Con = control, Rv = Ranavirus.

https://doi.org/10.1371/journal.pone.0130500.g003

[Figure omitted. See PDF.]

Fig 4. Counts of up- and down-regulated transcripts in response to Bd and Ranavirus.

Data are summarized at alternate false-discovery rate levels (FDR < 0.05 & FDR < 0.10) and following re-allocation of transcripts in the Bd vs. Ranavirus comparison.

https://doi.org/10.1371/journal.pone.0130500.g004

[Figure omitted. See PDF.]

Table 3. Differentially expressed transcripts (all comparisons, FDR = 0.05, log-fold-change = 1) that were successfully annotated.

https://doi.org/10.1371/journal.pone.0130500.t003

Five transcripts were significantly differentially expressed in both the Ranavirus vs. control and Bd vs. control sets (FDR<0.05). Three of these were up-regulated after exposure to both pathogens (including the two with the highest logFC change values of all transcripts in both treatments) and two were down-regulated. The AP-4 complex subunit sigma-1 is the only one of these five transcripts that was successfully annotated, and was associated with the highest logFC increase in both treatments (10.9 and 9.74 fold up-regulation in Bd, and Ranavirus treatments respectively). At the FDR threshold of p<0.10, ten transcripts were significantly differentially expressed in both the Bd vs. control and Ranavirus vs. control comparisons. Of these, four were up-regulated in both comparisons and six down-regulated. Four annotations in addition to the AP-4 complex were obtained, including (i) Cholesterol 7-alpha-monooxygenase which may be involved in xenobiotic metabolism, (ii) Acyl-coenzyme A synthetase ACSM3, a mitochondrial gene involved in lipid and/or fatty acid metabolism, (iii) Protein CutA homolog possibly involved in metal ion response, and (iv) Protein FAM136A, another mitochondrial protein. Differentially expressed transcripts with annotation are summarized in Table 3, and the full data is available in S1 File.

GO term enrichment

Enriched GO terms (Table 4; adjusted P-value<0.05) clustered in cell signalling, immunity, inflammation and metabolism. In total, 103 GO terms were significantly enriched in animals challenged with Bd (at 5% level after adjusting for multiple comparisons) compared to the reference set (see Table 4). When parent GO terms were examined, 14 of the 20 differentially expressed transcripts under the “metabolic process” parent GO term (GO:0008152) were up-regulated in Bd compared to controls. All transcripts related to “Biological regulation” (GO:0065007) were up-regulated. Immune system processes (GO:0006958), inflammatory responses (GO:0006954) and response to stimulus (GO:0050896) were also generally up-regulated (driven by complement activation) though there was also a down-regulation of Immunoglobulin J chain. No GO terms were enriched in the Ranavirus challenged animals, and no GO terms (regardless of treatment) were significantly enriched when differentially expressed transcripts in the FDR < 0.05 lists only were used for GO enrichment analyses.

[Figure omitted. See PDF.]

Table 4. Enriched GO terms associated with differentially expressed transcripts in the Bd vs. control comparison (adjusted P-value <0.05).

https://doi.org/10.1371/journal.pone.0130500.t004

Discussion

We investigated the response of an anuran host (R. temporaria) to the fungal pathogen Bd and the viral pathogen, Ranavirus, and detected a significant transcriptional response to Bd. Enriched GO terms involved the major arms of the immune response (innate and adaptive immunity and complement activation) as well as metabolic processes. Elements of the adaptive immune response were significantly differentially expressed in animals exposed to Bd and those exposed to Ranavirus, and this transcriptional response occurred at only four days post-exposure and before signs of disease consistent with ranavirosis were observed. Despite this, the overall response to Ranavirus was extremely limited. Variation between pools within treatments resulted in low power to detect differential expression when accounting for multiple comparisons (False Discovery Rate). However, this conservative approach allows confidence in the results and likely includes the genes with the largest transcriptional changes.

Perhaps the most notable result is the strong up-regulation of elements of the adaptive immunity, in response to either pathogen. The highest log fold-change for both pathogen treatments was successfully annotated to the AP-4 complex subunit sigma-1 gene (AP4S1). This is a subunit of a protein coat that is involved in targeting proteins from the trans-Golgi network to the endosomal-lysosomal system (http://www.uniprot.org/uniprot/Q9Y587). The trans-Golgi network is the location for the loading of cytokines (cytokines modulate the balance between the humoral and cell-based adaptive immune responses) with signal peptides into vesicles or carriers for delivery to the cell surface or other organelles [43]. The endosomal-lysosomal proteolysis system appears to be crucial in the immune system, in particular by binding class II MHC molecules to create ligands for antigen recognition by the T lymphocyte system. This process of antigen processing is important for immunity to pathogens as well as for the identification of self-peptides [44]. This result indicates that AP4S1 may be a particularly important candidate gene for future studies of amphibian host-pathogen biology. To date, we are unaware of any studies focusing on this gene in this context.

It is also notable that Interferon-stimulated 20 kDa exonuclease-like 2 and multiple versions of an Interferon-induced very large GTPase 1 were significantly up-regulated in the Bd treatment (Table 2). Interferon-induced very large GTPases are thought to mobilise effectors against a broad range of invading pathogens [45]. This result suggests that we have sampled these frogs after an initial innate response and at the point of mobilising other pathways. The result further suggests that these hosts are mounting a major phagocytic response to Bd challenge, an anti-fungal response which has previously been shown to be unimpaired during experimental chytrid infection [46].

We found a much stronger transcriptional response to Bd exposure than to Ranavirus. Our results broadly reflect what is seen in the wild. R. temporaria populations in the UK have undergone serious declines in the face of recurrent FV3-like Ranavirus die-offs [5] and are involved in multi-host mass mortality events associated with infection with a newly emerging Ranavirus lineage in Spain [7]. This host is highly susceptible to Ranavirus and it seems likely that infected individuals are failing to mount an effective immune response. In contrast, common frogs are considered relatively resistant to Bd [9], which may reflect effective immunity. Our observations of a weak transcriptional response to Ranavirus are therefore consistent with the higher pathogenicity of Ranavirus relative to Bd in this host as well as potential immune evasion of Ranavirus (reviewed in [17]).

However, the limited overall response to Ranavirus challenge that we detected was surprising in the light of previous work. Cotter et al. exposed Ambystomatid salamanders to ATV and used a custom microarray to measure host response in spleens at time-points between 24 hours and 6 days [15]. Ambystomatids are naturally infected with ATV in North America and are highly susceptible [47] but immune response (including phagocytosis, cytokine signalling, and complement activation) was detected as early as 24 hours post exposure and transcriptional response increased through the six day experimental period [15]. Differences in methodological approach and host-pathogen system might explain the contrasting findings. Tissue type, dose and exposure route all varied between studies and the host species investigated are highly divergent. The viruses utilised are also divergent. ATV is closely related to fish Ranaviruses and may represent a recent host jump [48]. ATV also seems more specialized for salamanders over other amphibians [49]. On the other hand, FV3-like viruses have been more frequently implicated in multi-host die-offs [2] (and may therefore represent better immune evaders) and are the more derived lineage of amphibian-like ranaviruses.

Previous work has shown increased expression of Immunoglobulin Y and activation-induced cytidine deaminase (AID) after Ranavirus (FV3) exposure in Xenopus laevis, indicating that B-cells are activated, and antibodies to Ranavirus have been detected two to six months post-exposure in this host species [50–52]. As well as this adaptive immune response, there was a rapid up-regulation of pro-inflammatory genes (Arginase 1, IL-1B, TNF-a, largely produced by macrophages), and an increased abundance of macrophages and antimicrobial peptides, which are known to inactivate FV3 in vitro [52–54]. Other studies indicate that we might expect to detect changes in MHC class I gene expression under Ranavirus infection, however it is clear that MHC expression is age-dependent. Pre-metamorphic Xenopus tadpoles do not express MHC class Ia genes, while adults do–in juveniles, the stage used in our experiment, expression is weaker than in adults [55]. We saw evidence for an adaptive immune response being initiated (AP4S1 up-regulation) but little evidence of further responses in our Ranavirus vs. control comparison.

Rosenblum et al. also reported significant transcriptional changes in Interferon-related genes in Bd-exposed animals, although they found many of these were only expressed at 16 days as opposed to 3 days post-exposure [12]. Our data demonstrates that in R. temporaria, these pathways are already active four days post-exposure. Bd may suppress T-cell mediated responses to infection [13] and resistance to Bd may be partly due to an ability to overcome this [11]. Here we have demonstrated the enrichment of lymphocyte and leukocyte mediated immunity in R. temporaria. In addition to immunological responses, Rosenblum et al. reported differential expression in a range of skin integrity, cellular stress, and homeostasis genes—the majority of these detected 16 days post-exposure [12]. Bd disruption of host skin integrity is a key disease process [10] and resistant hosts may mitigate this effect through increased expression of genes contributing to skin structure [11]. Although we sampled livers and not skin, we also observed this type of structural response to Bd through up-regulation of genes involved in fibroblast proliferation as well as up-regulation of actin and Galectin-3.

Previous work has shown that particular MHC class IIB alleles are associated with host resistance to Bd across a range of hosts varying in their susceptibility to the pathogen [56]. Rosenblum et al. reported down-regulation of MHC class I and II in the spleen and up-regulation of MHC class I and II in the skin in R. mucosa experimentally infected with Bd, however these were predominantly shown after 16 days post-exposure and not after only three days post-exposure [12]. We saw no changes in MHC expression in animals exposed to Bd.

A more relaxed false discovery rate threshold is expected to yield an increased number of differentially expressed genes but this effect appeared stronger in our Bd treatment. The relatively large number of genes withstanding the less stringent FDR correction (0.1) relative to the more conservative one (0.05) suggests that a large number of genes were differentially expressed between treatments but that the amplitude of this change was only modest. Our experimental design involved sampling animals at an early time-point, four days after exposure and prior to any observed mortality or observed signs of disease. This design has enabled identification of changes due to pathogen exposure rather than merely identifying differences between healthy and dying hosts. Previous experiments have pointed to reduced transcriptional responses at early time points post-exposure (3 days) relative to late time-points (16 days) in Rana species exposed to Bd, with the majority of transcriptional changes only becoming significant later [12]. Innate immune response to Ranavirus in X. laevis peaks at six days post-exposure, with adaptive immunity still detectable at 2–6 months post-exposure [57]. In contrast, we have demonstrated that certain adaptive host transcriptional responses do occur early for both pathogens examined.

This study demonstrates the utility of using RNAseq with non-model organisms to identify loci that may be important in host responses to pathogens. Amphibians are a highly threatened group, faced with catastrophic declines driven in part by emerging infectious diseases. Whole genome data is currently only available for a single amphibian genus (Xenopus). As such, our transcriptome data for R. temporaria provides a valuable resource for another genus [58], much-needed insight into amphibian immunity, as well as information on specific responses to the two most important multi-host pathogens affecting amphibians. We were also able to identify candidate genes that could serve as markers for understanding the impacts of disease in wild populations and the adaptive potential of populations that are under threat. While this study has provided crucial insights into amphibian gene expression following exposure to pathogens, comparisons across life-history stages, time points post-exposure, and source populations with different previous infection statuses and genetic diversity will all be necessary to gain a more complete picture of the transcriptional responses to amphibian disease. We foresee this line of research offering exciting possibilities for the selection of individuals with disease resistance for captive breeding programmes for conservation.

Supporting Information

[Figure omitted. See PDF.]

S1 Text. Re-allocating differentially expressed transcripts in the Bd vs. Ranavirus comparison.

https://doi.org/10.1371/journal.pone.0130500.s001

(DOCX)

S1 Table. Sample RNA concentrations.

https://doi.org/10.1371/journal.pone.0130500.s002

(DOCX)

S2 Table. CEGMA output summary; results of filtering assembly on CEG coverage.

A) Full assembly, B) Filtered assembly on FPKM> = 1 for all replicates within at least one treatment.

https://doi.org/10.1371/journal.pone.0130500.s003

(DOCX)

S1 Fig. Expression pattern of re-allocated Bd vs. Ranavirus transcripts.

Annotated transcripts from the Bd vs. Ranavirus (FDR<0.10) comparison that were also found in either Bd vs. control, Ranavirus vs. control or both prior to FDR filtering (protein name of best blast hit given).

https://doi.org/10.1371/journal.pone.0130500.s004

(DOCX)

S1 File. Full list of differentially transcribed genes

https://doi.org/10.1371/journal.pone.0130500.s005

(CSV)

Acknowledgments

This work was carried out under Home Office license (Project Licence numbers PPL 80/2214 and PPL 80/2466) and was subject to review by the Institute of Zoology Ethics Committee and the University of Exeter Ethical Review Board.

We thank the Exeter Sequencing Service, the National Institute for Health Research University College London Hospitals Biomedical Research Centre, and Computational core facilities at the University of Exeter. We also thank Chris Durrant for help in preparing samples for sequencing and two anonymous reviewers whose comments improved the manuscript.

Author Contributions

Conceived and designed the experiments: AGFG TWJG FB. Performed the experiments: CR AGFG SJP TWJG. Analyzed the data: KHP KM CR SJP FB. Contributed reagents/materials/analysis tools: KHP KM. Wrote the paper: SJP AGFG.

Citation: Price SJ, Garner TWJ, Balloux F, Ruis C, Paszkiewicz KH, Moore K, et al. (2015) A de novo Assembly of the Common Frog (Rana temporaria) Transcriptome and Comparison of Transcription Following Exposure to Ranavirus and Batrachochytrium dendrobatidis. PLoS ONE 10(6): e0130500. https://doi.org/10.1371/journal.pone.0130500

References

1. Wake DB, Vredenburg VT. Are we in the midst of the sixth mass extinction? A view from the world of amphibians. Proc Natl Acad Sci U S A. 2008;105: 11466–11473. pmid:18695221

2. Green DE, Converse KA, Schrader AK. Epizootiology of sixty-four amphibian morbidity and mortality events in the USA, 1996–2001. Ann N Y Acad Sci. 2002;969: 323–339. pmid:12381613

3. Pounds JA, Bustamante MR, Coloma LA, Consuegra JA, Fogden MPL, Foster PN, et al. Widespread amphibian extinctions from epidemic disease driven by global warming. Nature. 2006;439: 161–167. pmid:16407945

4. Skerratt LF, Berger L, Speare R, Cashins S, McDonald KR, Phillott AD, et al. Spread of chytridiomycosis has caused the rapid global decline and extinction of frogs. Ecohealth. 2007;4: 125–134.

5. Teacher AGF, Cunningham AA, Garner TWJ. Assessing the long-term impact of Ranavirus infection in wild common frog populations. Anim Conserv. 2010;13: 514–522.

6. Kriger KM, Hero J-M. Chytridiomycosis, Amphibian Extinctions, and Lessons for the Prevention of Future Panzootics. Ecohealth. 2009;6: 6–10. pmid:19421815

7. Price SJ, Garner TWJ, Nichols RA, Balloux F, Ayres C, Mora-Cabello de Alba A, et al. Collapse of Amphibian Communities Due to an Introduced Ranavirus. Curr Biol. 2014;24: 2586–2591. pmid:25438946

8. Cunningham AA, Langton TES, Bennett PM, Lewin JF, Drury SEN, Gough RE, et al. Pathological and microbiological findings from incidents of unusual mortality of the common frog (Rana temporaria). Philos Trans R Soc Lond Ser B-Biol Sci. 1996;351: 1539–1557.

9. Balaz V, Voeroes J, Civis P, Vojar J, Hettyey A, Sos E, et al. Assessing Risk and Guidance on Monitoring of Batrachochytrium dendrobatidis in Europe through Identification of Taxonomic Selectivity of Infection. Conserv Biol. 2014;28: 213–223. pmid:24033675

10. Voyles J, Young S, Berger L, Campbell C, Voyles WF, Dinudom A, et al. Pathogenesis of Chytridiomycosis, a Cause of Catastrophic Amphibian Declines. Science. 2009;326: 582–585. pmid:19900897

11. Ellison AR, Tunstall T, DiRenzo GV, Hughey MC, Rebollar EA, Belden LK, et al. More than skin deep: functional genomic basis for resistance to amphibian chytridiomycosis. Genome Biol Evol. 2014; evu285. 10.1093/gbe/evu285

12. Rosenblum EB, Poorten TJ, Settles M, Murdoch GK. Only skin deep: shared genetic response to the deadly chytrid fungus in susceptible frog species. Mol Ecol. 2012;21: 3110–3120. pmid:22332717

13. Ellison AR, Savage AE, DiRenzo GV, Langhammer P, Lips KR, Zamudio KR. Fighting a losing battle: vigorous immune response countered by pathogen suppression of host defenses in the chytridiomycosis-susceptible frog Atelopus zeteki. G3 Bethesda Md. 2014;4: 127575hesda

14. Ribas L, Li M-S, Doddington BJ, Robert J, Seidel JA, Kroll JS, et al. Expression Profiling the Temperature-Dependent Amphibian Response to Infection by Batrachochytrium dendrobatidis. Plos One. 2009;4: e8408. pmid:20027316

15. Cotter JD, Storfer A, Page RB, Beachy CK, Voss SR. Transcriptional response of Mexican axolotls to Ambystoma tigrinum virus (ATV) infection. Bmc Genomics. 2008;9: 493. pmid:18937860

16. Morales HD, Abramowitz L, Gertz J, Sowa J, Vogel A, Robert J. Innate Immune Responses and Permissiveness to Ranavirus Infection of Peritoneal Leukocytes in the Frog Xenopus laevis. J Virol. 2010;84: 4912–4922. pmid:20200236

17. Grayfer L, Andino FDJ, Chen G, Chinchar GV, Robert J. Immune Evasion Strategies of Ranaviruses and Innate Immune Responses to These Emerging Pathogens. Viruses-Basel. 2012;4: 1075–1092.

18. Teacher AGF, Garner TWJ, Nichols RA. Evidence for Directional Selection at a Novel Major Histocompatibility Class I Marker in Wild Common Frogs (Rana temporaria) Exposed to a Viral Pathogen (Ranavirus). Plos One. 2009;4: e4616. pmid:19240796

19. Rollins-Smith LA. Metamorphosis and the amphibian immune system. Immunol Rev. 1998;166: 221–230. pmid:9914915

20. Cunningham AA, Hyatt AD, Russell P, Bennett PM. Emerging epidemic diseases of frogs in Britain are dependent on the source of ranavirus agent and the route of exposure. Epidemiol Infect. 2007;135: 1200–1212. pmid:17181914

21. Reed LJ, Muench H. A Simple Method of Estimating Fifty Per Cent Endpoints,. Am J Epidemiol. 1938;27: 493–497.

22. Farrer RA, Weinert LA, Bielby J, Garner TWJ, Balloux F, Clare F, et al. Multiple emergences of genetically diverse amphibian-infecting chytrids include a globalized hypervirulent recombinant lineage. Proc Natl Acad Sci U S A. 2011;108: 18732–18736. pmid:22065772

23. Miller D, Gray M, Storfer A. Ecopathology of Ranaviruses Infecting Amphibians. Viruses. 2011;3: 235151ng Am pmid:22163349

24. Aronesty E. ea-utils : “Command-line tools for processing biological sequencing data” [Internet]. 2011. Available: http://code.google.com/p/ea-utils

25. Babraham Bioinformatics. FastQC A Quality Control tool for High Throughput Sequence Data [Internet]. 2015. Available: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

26. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8: 1494–1512. pmid:23845962

27. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12: 323. pmid:21816040

28. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29: 644–U130. pmid:21572440

29. Finn RD, Clements J, Eddy SR. HMMER web server: interactive sequence similarity searching. Nucleic Acids Res. 2011;39: W29–W37. pmid:21593126

30. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40: D290–D301. pmid:22127870

31. Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat Methods. 2011;8: 785–786. pmid:21959131

32. Krogh A, Larsson B, von Heijne G, Sonnhammer ELL. Predicting transmembrane protein topology with a hidden Markov model: Application to complete genomes. J Mol Biol. 2001;305: 567–580. pmid:11152613

33. Altschul S, Gish W, Miller W, Myers E, Lipman D. Basic Local Alignment Search Tool. J Mol Biol. 1990;215: 403–410. pmid:2231712

34. Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012;40: D109–D114. pmid:22080510

35. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25: 25–29. pmid:10802651

36. Powell S, Szklarczyk D, Trachana K, Roth A, Kuhn M, Muller J, et al. eggNOG v3.0: orthologous groups covering 1133 organisms at 41 different taxonomic ranges. Nucleic Acids Res. 2012;40: D284–D289. pmid:22096231

37. Lagesen K, Hallin P, Rodland EA, Staerfeldt H-H, Rognes T, Ussery DW. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 2007;35: 3100–3108. pmid:17452365

38. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28: 511–U174. pmid:20436464

39. Parra G, Bradnam K, Korf I. CEGMA: a pipeline to accurately annotate core genes in eukaryotic genornes. Bioinformatics. 2007;23: 1061–1067. pmid:17332020

40. Parra G, Bradnam K, Ning Z, Keane T, Korf I. Assessing the gene space in draft genomes. Nucleic Acids Res. 2009;37: 289–297. pmid:19042974

41. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995; 28999995

42. Maere S, Heymans K, Kuiper M. BiNGO: a Cytoscape plugin to assess overrepresentation of Gene Ontology categories in Biological Networks. Bioinformatics. 2005;21: 3448–3449. pmid:15972284

43. Lacy P, Stow JL. Cytokine release from innate immune cells: association with diverse membrane trafficking pathways. Blood. 2011;118: 9–18. pmid:21562044

44. Watts C. The endosome-lysosome pathway and information generation in the immune system. Biochim Biophys Acta-Proteins Proteomics. 2012;1824: 14–21.

45. Kim B-H, Shenoy AR, Kumar P, Bradfield CJ, MacMicking JD. IFN-Inducible GTPases in Host Cell Defense. Cell Host Microbe. 2012;12: 432–444. pmid:23084913

46. Fites JS, Ramsey JP, Holden WM, Collier SP, Sutherland DM, Reinert LK, et al. The Invasive Chytrid Fungus of Amphibians Paralyzes Lymphocyte Responses. Science. 2013;342: 366–369. pmid:24136969

47. Jancovich JK, Davidson EW, Morado JF, Jacobs BL, Collins JP. Isolation of a lethal virus from the endangered tiger salamander Ambystoma tigrinum stebbinsi. Dis Aquat Organ. 1997;31: 161–167.

48. Jancovich JK, Bremont M, Touchman JW, Jacobs BL. Evidence for Multiple Recent Host Species Shifts among the Ranaviruses (Family Iridoviridae). J Virol. 2010;84: 2636–2647. pmid:20042506

49. Schock DM, Bollinger TK, Chinchar VG, Jancovich JK, Collins JP. Experimental evidence that amphibian ranaviruses are multi-host pathogens. Copeia. 2008;2008: 133–143.

50. Chinchar VG, Wang J, Murti G, Carey C, Rollins-Smith L. Inactivation of frog virus 3 and channel catfish virus by esculentin-2P and ranatuerin-2P, two antimicrobial peptides isolated from frog skin. Virology. 2001;288: 351–357. pmid:11601906

51. Chinchar VG, Bryan L, Silphadaung U, Noga E, Wade D, Rollins-Smith L. Inactivation of viruses infecting ectothermic animals by amphibian and piscine antimicrobial peptides. Virology. 2004;323: 268–275. pmid:15193922

52. Chen G, Robert J. Antiviral Immunity in Amphibians. Viruses-Basel. 2011;3: 2065–2086.

53. Marr S, Morales H, Bottaro A, Cooper M, Flajnik M, Robert J. Localization and differential expression of activation-induced cytidine deaminase in the amphibian Xenopus upon antigen stimulation and during early development. J Immunol. 2007;179: 6783–6789. pmid:17982068

54. Panchanathan V, Chaudhri G, Karupiah G. Antiviral protection following immunization correlates with humoral but not cell-mediated immunity. Immunol Cell Biol. 2010;88: 461–467. pmid:20066003

55. Flajnik M, Hsu E, Kaufman J, Dupasquier L. Changes in the Immune-System During Metamorphosis of Xenopus. Immunol Today. 1987;8: 58–64. pmid:25291685

56. Savage AE, Zamudio KR. MHC genotypes associate with resistance to a frog-killing fungus. Proc Natl Acad Sci U S A. 2011;108: 16705–16710. pmid:21949385

57. Savage AE, Kiemnec-Tyburczy KM, Ellison AR, Fleischer RC, Zamudio KR. Conservation and divergence in the frog immunome: Pyrosequencing and de novo assembly of immune tissue transcriptomes. Gene. 2014;542: 98–108. pmid:24680726

58. Savage AE, Kiemnec-Tyburczy KM, Ellison AR, Fleischer RC, Zamudio KR. Conservation and divergence in the frog immunome: Pyrosequencing and de novo assembly of immune tissue transcriptomes. Gene. 2014;542: 98–108. pmid:24680726