Comparative genomics of an emerging Amphibian virus 1 2 3 Brendan Epstein* and Andrew ...

50 Pages·2015·0.93 MB·English

G3: Genes|Genomes|Genetics Early Online, published on November 3, 2015 as doi:10.1534/g3.115.023762
1 Comparative genomics of an emerging Amphibian virus
2
3
4 Brendan Epstein* and Andrew Storfer*
5
6 *School of Biological Sciences
7 Washington State University
8 PO BOX 644236
9 Pullman, WA 99164
10
11
12
13 Data in public repositories
14 The data set generated for this article is available in GenBank: BioProject PRJNA257291.
15
1
© The Author(s) 2013. Published by the Genetics Society of America.
16 Running title: ATV Ranavirus genomics
17
18 Keywords: Ranavirus, Ambystoma tigrinum virus, range expansion
19
20 Corresponding Author:
21 Andrew Storfer, [email protected], 509-335-7922
22 School of Biological Sciences
23 Washington State University
24 PO BOX 644236
25 Pullman, WA 99164
26
2
27 Abstract
28 Ranaviruses, a genus of the Iridoviridae, are large double-stranded DNA viruses that
29 infect cold-blooded vertebrates world-wide. Ranaviruses have caused severe epizootics in
30 commercial frog and fish populations and are currently classified as notifiable pathogens in
31 international trade. Previous work shows that a ranavirus that infects tiger salamanders
32 throughout Western North America (Ambystoma tigrinum virus, or ATV) is in high prevalence
33 among salamanders in the fishing bait trade. Bait ATV strains have elevated virulence and are
34 transported long distance by humans, providing widespread opportunities for pathogen pollution.
35 We sequenced the genomes of 15 strains of ATV collected from tiger salamanders across western
36 North America and performed phylogenetic and population genomic analyses and tests for
37 recombination. We find that ATV forms a monophyletic clade within the rest of the Ranaviruses
38 and that it likely emerged within the last several thousand years, before human activities
39 influenced its spread . We also identify several genes under strong positive selection, some of
s
40 which appear to be involved in viral virulence and/or host immune evasion. In addition, we
41 provide support for the pathogen pollution hypothesis with evidence of recombination among
42 ATV strains and potential bait-endemic strain recombination.
43
44 Introduction
45 Emerging infectious diseases are increasingly appreciated as a leading health concern for
46 humans, wildlife, and economically important agricultural populations (Daszak et al. 2000;
47 Smith et al. 2009). Indeed, pathogens are now listed as a leading cause of species’ declines and
48 extinctions (De Castro and Bolker 2005; Smith et al. 2006). Ranaviruses, a genus of the
3
49 Iridoviridae, are globally-distributed pathogens of amphibians, reptiles and commercial fish
50 species (Chinchar 2002; Miller et al. 2011). These large, double-stranded DNA viruses are
51 considered emerging due to increases in incidence and geographic range over the last 30 years
52 (Daszak et al. 2000; Chinchar 2002; Miller et al. 2011).
53 Ranaviruses are now classified as notifiable pathogens in international trade because o
f
54 their effects on commercial and wildlife populations (Schloegel et al. 2010). Pathogen pollution
55 is of particular concern, whereby non-native ranavirus strains are introduced into host
56 populations with which they have no evolutionary history and potentially leading to large scale
57 epizootics (Cunningham et al. 2003). Phylogenetic analyses provide several lines of evidence for
58 host switching events among the ten completely sequenced Ranavirus genomes (Jancovich et al.
59 2010; Abrams et al. 2013). These analyses suggest that TFV (tiger frog virus) and GIV (grouper
60 iridovirus) are likely strains of the geographically-distant FV3 (frog virus 3) and SGIV
61 (Singapore grouper iridovirus), respectively (Chinchar et al. 2011). In addition, strains isolated
62 from different vertebrate classes, such as STIV (soft-shelled turtle iridovirus; Reptilia) and FV3
63 (Amphibia) are very similar in genome organization and, like RGV (Lei et al. 2012), both have
64 truncated versions of the viral homolog of eukaryotic translational initiation factor 2α (vIF-2α
65 [Huang et al. 2009; Tan et al. 2004]). Additional phylogenetic analyses, combined with
66 comparisons of genomic organization suggest that the most recent common ancestor of all
67 ranaviruses was a strain that infected fish (Jancovich et al. 2010; Chinchar et al. 2011). Evidence
68 of a fish-amphibian host switch comes from the strong co-linearity of EHNV (Epizootic
69 haematopoietic necrosis virus; isolated from rainbow trout) and ATV (Ambystoma tigrinum virus,
70 isolated from salamanders). A possible mechanism for host switching among vertebrate classes
4
71 was inferred from evidence of positive selection on 12 ranavirus genes; six of these genes were
72 apparently acquired as a result of host switches (Abrams et al. 2013).
73 The broad host range and propensity for host switching among ranavirus strains has led to
74 widespread concerns surrounding the potential effects of pathogen pollution. Two ranavirus
75 strains that infect amphibians have received considerable attention in this regard. FV3 and FV3-
76 like strains, which have a broad host range, are commonly found when international shipments of
77 bullfrogs are tested using PCR-based methods (Schloegel et al. 2009). In addition, ATV - a
78 ranavirus that infects tiger salamanders across the North American cordillera (Jancovich et al.
79 1997; Jancovich et al. 2003; Jancovich et al. 2005) - has been found to have > 80% prevalence in
80 tiger salamanders in bait shops across Arizona, Colorado and New Mexico (Picco and Collins
81 2008). ATV causes seasonally recurrent epidemics in primarily larval tiger salamander
82 populations that inhabit natural and human-made ponds throughout western North America
83 ([(Jancovich et al. 1997; Jancovich et al. 2005 ; Brunner et al. 2004]; Figure 1). Phylogeographic
84 analyses suggest that long-distance ATV dispersal follows highways that are known to be used
85 for transporting bait (Jancovich et al. 2005), and angler surveys confirmed release of unused live
86 bait salamanders into lakes and ponds (Picco and Collins 2008). Thus, extensive opportunities
87 exist for pathogen pollution by ATV.
88 A reference genome of ATV has been sequenced (Jancovich et al. 2003), and
89 phylogenetic analysis of ATV strains (using DNA sequences comprising about 2% of the
90 genome) throughout the western United States suggests monophyly (Jancovich et al. 2005;
91 Storfer et al. 2007). However, ATV appears to evolve differently in different host populations.
92 That is, of nine genes (ORFs) that have a putative function in host immune evasion and/or viral
5
93 virulence (Jancovich et al. 2003), two were under positive selection among tiger salamander host
94 populations (Ridenhour and Storfer 2008). Indeed, there is also evidence of coevolution of ATV
95 and tiger salamander hosts as comparative phylogenetic analyses show strong concordance
96 between host and virus phylogenetic trees (Storfer et al. 2007). Thus, it is important to
97 understand genomic variation among ATV strains to help assess the potential effects of pathogen
98 pollution as coevolutionary dynamics of hosts and pathogens could be disrupted by introduction
99 of non-native bait strains.
100 Here, we sequenced the complete genomes of 15 ATV native and bait-associated strains
101 across a broad geographic area to assess genomic diversity, identify genes that are experiencing
102 strong positive selection, and to test for evidence of recombination. Recombination has been
103 reported to occur among ranaviruses (e.g. Abrams et al. (2013)), but has not been verified or
104 quantified. A previous study showed that an ATV strain isolated from a bait shop had
105 significantly greater virulence than endemic, coevolved strains (Storfer et al. 2007). Thus,
106 recombination among non-native and native ATV strains could enhance virulence, leading to
107 severe epizootics and possibly local extirpation. We also conducted phylogenetic analyses of
108 ATV in comparison to other Ranavirus genomic sequences to provide a rigorous test of
109 monophyly versus a pattern of non-monophyly resulting from repeated host-switches, as well as
110 timing of the emergence of ATV to assess human involvement in its geographic spread.
111
112 Materials and Methods
113 Strain collection
114 The genome sequences of eleven previously sequenced Ranavirus strains were obtained
6
115 from Genbank: Ambystoma tigrinum virus (ATV), Andrias davidianus ranavirus (ADRV),
116 common midwife toad virus (CMTV), tiger frog virus (TFV), frog virus 3 (FV3), Rana grylio
117 virus (RGV), soft-shelled turtle iridovirus (STIV), european sheatfish virus (ESV), epizootic
118 haematopoietic virus (EHNV), grouper iridovirus (GIV), and Singapore grouper iridovirus
119 (SGIV) (Table 1). The genome sequences of 15 additional Ranavirus isolates collected from
120 North America between 1997 and 2003 are described for the first time here (Table 1). Reads and
121 assemblies for these 15 strains have been deposited at NCBI in BioProject PRJNA257291.
122
123 DNA Extraction
124 ATV strains were isolated from infected tiger salamanders at each of the 15 localities by
125 first freezing the entire salamander at -70°C. Then, each salamander was homogenized in Eagle’s
126 Minimum Essential Medium (MEM) and slurry was filtered in preparation for viral cell culture.
127 Cell culture protocols for replication of ATV have been standardized (Jancovich et al. 1997;
128 Jancovich et al. 2001); the 15 virus isolates to be used herein all produced cytopathology in
129 cultured Epithelioma papilloma cyprini fish cells (Fijan et al. 1983) within 48-72 hours. Virus
130 strains were triple plaque-purified and cultivated with no more than three total passages and
131 equal number of passages per strain, because serial passages through clonal cell lines could result
132 in viral evolution, such as adaptation to the cell line and consequent attenuation in the host (see
133 Ebert (1999)). Each strain was replicated to approximately 100 ml in volume in MEM and viral
134 DNA extraction followed a modified sucrose-gradient protocol (Thurber et al. 2009). In brief,
135 DNA was extracted by lysing the host cells with three freeze-thraw cycles, then incubating with
136 DNase, layering the DNA onto a 25% sucrose pad, and spinning at 16,000 rpm for eight hours in
7
137 an ultracentrifuge. The DNA was then purified by treatment with proteinase K and RNase
138 followed by a phenol extraction.
139
140 Sequencing, assembly, and annotation
141 Libraries were prepared using standard protocols and ATV strains were sequenced at the
142 University of Idaho Genomics Resources Core using 300 bp paired-end reads on an Illumina
143 MiSeq, resulting in a mean of 670,000 reads per strain (Table 1). The reads were cleaned with
144 seqyclean (available from https://bitbucket.org/izhbannikov/seqyclean) using ESTs from
145 Pimephales promelas as contaminants, then overlapping read pairs were merged using Flash
146 (v1.2.8; Magoč and Salzberg 2011). A reference-guided assembly was performed using ARC
147 (Hunter et al. 2015), with both ATV and FV3 as the reference genome in separate analyses,
148 followed by manual finishing. For most strains, the best assembly (most continuous) was
149 obtained using SPAdes (v3.0.0; Bankevich et al. 2012) as the assembler, except for ORV and
150 JMH (see Table 1), for which we used Newbler (Marguiles et al. 2005). We obtained assistance
151 with the ORV and JMH assemblies from the University of Idaho Genomics Core. To correct
152 errors in the assembly and remove any reference bias, we aligned the reads for each strain back
153 to its assembly. Any position not covered by ≥ 10 uniquely aligning reads, and with a single base
154 call supported by at least 80% of the uniquely aligning reads with a quality score ≥ 30, and with
155 support on reads aligning to both strands was changed to “N.” If the aligned reads supported a
156 different base than the assembly, the assembly was changed. We relaxed the requirement for
157 alignment to both strands for the first and last 100 bp of the assembly, and we trimmed any
158 remaining terminal ambiguities. The final fold coverage of the assemblies ranged from
8
159 approximately 350 to 2300.
160 Although relatively long, paired-end reads and very high coverage should result in high
161 accuracy, we verified the assembled sequence at seven loci, including two repetitive regions,
162 using 50, 833 bp of Sanger sequence data from 14 strains (no Sanger data was available for
163 DOT312) (Ridenhour and Storfer 2008). These data included two variable repeat regions.
164 We annotated the genomes using the Rapid Annotation Transfer Tool (RATT; Otto et al.
165 2011), with ATV, FV3, and CMTV as the reference genomes, ATG as the start codon, TGA,
166 TAA, and TAG as stop codons, no splice site correction, and no pseudogene correction. We used
167 a custom python script to remove potentially redundant annotations (we picked the longest ORF
168 for annotations with the same stop codon). We also removed 20-30 ORFs that were split into
169 multiple pieces, because none of the amphibian-like ranavirus (ALRV) sequences on from
170 Genbank have introns or split features. BLASTP (v2.2.28; Camacho et al. 2009) was used to
171 search Swissprot (downloaded 9 April 2014; Boeckmann et al. 2003) to find more descriptive
172 annotations for the genes.
173 Sets of orthologous genes were identified using blastclust (Camacho et al. 2009) using a
174 threshold of 80% amino acid identity across 80% of the length of both sequences; to construct
175 the genealogy with GIV strains, we required only 60% identity across 60% of the length of both
176 sequences. The protein sequences of orthologous genes were aligned using muscle (Edgar 2004),
177 and the protein sequence alignment was then used to construct the nucleotide alignment. We also
178 aligned the nucleotide sequences separately from the protein sequences and considered indels
179 that were not a multiple of three base pairs long to be potential frameshifts. To avoid frameshift
180 mutations being erroneously treated as non-synonymous mutations, alignments in genes with
9
181 potential frameshifts were trimmed to the last codon before the frameshift before running
182 phylogenetic or population genomic analyses.
183
184 Dotplots
185 Pairwise genome comparisons were conducted using the dnadiff program from MUMmer
186 (Kurtz et al. 2004), and dotplots were drawn using custom R (R Core Team 2014) code. Dotplots
187 are used for visualizing structural differences between pairs of genomes: the axes are position in
188 each genome, and a dot is placed at every location where the two genomes have similar
189 sequence. Identical genomes appear as a diagonal line, while inversions appear as changes in the
190 orientation of matches.
191
192 Genealogical analyses
193 Genealogical analyses were conducted on the concatenated nucleotide sequences of the
194 52 genes that were found in all ALRV strains using BEAST (v1.8.0; Drummond et al. 2012) and
195 MrBayes (v3.2.2; Ronquist et al. 2012). Based on model-selection tests conducted with
196 jModelTest (v2.1.4; Darriba et al. 2012) to determine the most likely model of molecular
197 evolution, we used the general time-reversible model with a gamma shape parameter (GTR+G)
198 and four site categories. In BEAST we used a lognormal relaxed clock, empirical base
199 frequencies, and a Bayesian skyline coalescent tree prior, on data partitioned by codon position.
200 The tip date option was used to provide dating information; the tip dates were entered as shown
201 in Table 1 with a precision (the accuracy of the collection date) of 0.5, except for EHNV, ESV,
202 and FV3, for which we entered a precision of 2.0, and TFV for which we entered a precision of
10