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

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.


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][3][4][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 F ST ) 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  [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.

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 TCID 50 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 TCID 50 /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].
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.

Sequencing, Assembly and annotation
A total of 1.29 x 10 9 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%.
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. Exposure to Bd elicited far greater transcriptional divergence from controls than did exposure to 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 Interferoninduced very large GTPase 1 proteins (Table 3).
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.

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 postexposure 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 activationinduced 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][51][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][53][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 S1 Text. Re-allocating differentially expressed transcripts in the Bd vs. Ranavirus comparison.