Genetic characterization of Tribeč virus and Kemerovo virus, two tick-transmitted human-pathogenic Orbiviruses

Virology 423 (2012) 68–76

Contents lists available at SciVerse ScienceDirect

Virology journal homepage: www.elsevier.com/locate/yviro

Genetic characterization of Tribeč virus and Kemerovo virus, two tick-transmitted human-pathogenic Orbiviruses Meik Dilcher a,⁎, Lekbira Hasib a, Marcus Lechner b, Nicolas Wieseke c, Martin Middendorf c, Manja Marz b, Andrea Koch a, Martin Spiegel a, Gerhard Dobler d, Frank T. Hufert a, Manfred Weidmann a a

Department of Virology, University Medical Center Göttingen, Kreuzbergring 57, D-37075 Göttingen, Germany RNA Bioinformatics Group, Philipps-University Marburg, Marbacher Weg 6, D-35037 Marburg, Germany c Faculty of Mathematics and Informatics, University of Leipzig, Johannisgasse 26, D-04103 Leipzig, Germany d Institute of Microbiology of the Armed Forces, Department of Virology and Rickettsiology, Neuherbergstrasse 1, D-80937 Munich, Germany b

a r t i c l e

i n f o

Article history: Received 26 May 2011 Returned to author for revision 17 September 2011 Accepted 17 November 2011 Available online 20 December 2011 Keywords: Tribeč virus Kemerovo virus Reovirus Orbivirus Co-evolution

a b s t r a c t We determined the complete genome sequences of Tribeč virus (TRBV) and Kemerovo virus (KEMV), two tick-transmitted Orbiviruses that can cause diseases of the central nervous system and that are currently classified into the Great Island virus serogroup. VP2 proteins of TRBV and KEMV show very low sequence similarity to the homologous VP4 protein of tick-transmitted Great Island virus (GIV). The new sequence data support previous serological classification of these Orbiviruses into the Kemerovo serogroup, which is different from the Great Island virus serogroup. Genome segment 9 of TRBV and KEMV encodes several overlapping ORF's in the + 1 reading frame relative to VP6(Hel). A co-phylogenetic analysis indicates a host switch from insect-borne Orbiviruses toward Ixodes species, which is in disagreement with previously published data. © 2011 Elsevier Inc. All rights reserved.

Introduction Tribeč virus (TRBV) was originally isolated from Ixodes ricinus ticks and from the blood of small rodents in the Tribeč mountains (Slovakia) in 1963 (Gresikova et al., 1965; Libikova et al., 1964). Lipovnik virus (LIPV), another closely related Orbivirus was isolated in the same year from I. ricinus ticks in Lipovnik village (Slovakia) (Libikova et al., 1964, 1965). In neutralization tests a strong neutralization index was obtained with antiserum against Kemerovo virus (KEMV, Gresikova et al., 1965), an Orbivirus isolated from I. persulcatus ticks and cerebrospinal fluid of two patients with encephalitis in western Siberia in 1962 (Chumakov et al., 1963). In addition, KEMV was also isolated from the blood of a migrating redstart (Phoenicurus phoenicurus) in Egypt in 1971 (Schmidt and Shope, 1971), which indicates a possible long distance spread of viruses of the Kemerovo complex via birds. TRBV and KEMV caused meningitis in experimentally infected rhesus macaques (Gresikova et al., 1966; Libikova et al., 1970) and mengingo-encephalitis and polyradiculitis in humans have been linked to viruses of the Kemerovo complex in former Czechoslovakia (Libikova et al., 1978); (Malkova et al., 1980) and in the Kemerovo region of Russia (Libikova et al., 1978). In addition, neutralizing antibodies against TRBV and LIPV were found in human sera in Slovakia (Gresikova et al., 1965, 1966; Malkova et al., 1980). Therefore, ⁎ Corresponding author. Fax: + 49 551 3910552. E-mail address: [email protected] (M. Dilcher). 0042-6822/$ – see front matter © 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.virol.2011.11.020

TRBV, LIPV and KEMV could be important causes of diseases of the central nervous system similar to tick-borne encephalitis virus (TBEV) and might be the causative agents of many patients with encephalitis of unknown etiology. This is emphasized by the observation, that in a tick-borne encephalitis focus near Rožňava, former Czechoslovakia, the prevalence of TRBV in I. ricinus was five times that of TBEV (Libikova et al., 1965). There is also evidence that these arboviruses might be involved in the etiology and course of chronic neurological diseases like multiple sclerosis (MS), since antibodies against LIPV were found in CSF of >50% of MS cases (Libikova et al., 1978). A TRBV range from Siberia to central Europe is evident from virus isolates from ticks and antibodies detected in animals (Dobler, 1996; Dobler et al., 2006; Süss and Schrader, 2004). The genus Orbivirus contains 22 serogroups (species) and at least 160 different serotypes (strains) and represents one of 12 genera within the family Reoviridae (Fauquet et al., 2005). Some Orbiviruses are transmitted by insects (midges, flies, mosquitoes) while others are transmitted by ticks. Their genomes consist of 10 segments of double-stranded RNA (dsRNA), which code for at least 10 viral proteins. Seven of these viral proteins are structural and 3 nonstructural. Orbiviruses have no envelope, a double-shelled icosahedral capsid (Attoui et al., 2005) and include pathogenic agents of man (KEMV), domestic animals (Bluetongue virus (BTV) and African horse sickness virus (AHSV)) and wild animals (Epizootic haemorrhagic disease virus (EHDV)) (Gorman, 1979). Type species of the genus is BTV, which is transmitted by Culicoides midges. So far, the insect-borne

M. Dilcher et al. / Virology 423 (2012) 68–76

Orbiviruses for which many sequences are available are much better characterized than tick-transmitted Orbiviruses for which few sequences have been described: Broadhaven virus (BRDV, partial) (Moss et al., 1992), Sandy Bay Virus (SBaV, partial (formerly Nugget virus)) (Belhouchet et al., 2010; Doherty et al., 1975; Gorman et al., 1984; Major et al., 2009), St Croix River virus (SCRV, complete genome) (Attoui et al., 2001). Recently, the complete genome sequence of Great Island virus (GIV) as well as segments 1, 2 and 6 of TRBV, KEMV and LIPV have been determined in a classical sequencing approach (Belhouchet et al., 2010). GIV, BRDV and SBaV infect seabirds but not humans. The original geographical distribution of BRDV is Scotland, that of GIV Newfoundland (Canada), whereas SBaV was so far only isolated from ticks on Macquarie Island, a remote subantarctic island (Major et al., 2009). SCRV was isolated from the eggs of I. scapularis, but so far no vertebrate host is known (Attoui et al., 2001). Originally, TRBV, LIPV and KEMV were grouped into the “KEM” subgroup of the Kemerovo serogroup of tick-borne Orbiviruses based on complement fixation data (Gorman et al., 1983). In addition to “KEM” this serogroup also contained the following subgroups: Chenuda “GNU”, Great Island “GI” and Wad Medani “WM”. Later on, KEMV, TRBV and LIPV together with Karagysh virus (KHAV, Skoferts et al., 1972) from Moldova were classified into a newly-formed Kemerovo serogroup of Orbiviruses apart from the Great Island virus serogroup based on limited genome reassortment with members of the “GI” subgroup (Labuda and Nuttall, 2004; Nuttall and Moss, 1989). However, currently the Kemerovo serogroup doesn't exist anymore and TRBV, KEMV and LIPV are grouped together with GIV, BRDV and SBaV into the Great Island virus serogroup which contains 36 different serotypes (Fauquet et al., 2005). So far only partial sequence information was available for ticktransmitted Orbiviruses causing human disease. We decided to sequence the genomes of TRBV and KEMV to investigate differences in host adaptation of insect- and tick-transmitted Orbiviruses and differences between human-pathogenic and non-pathogenic Orbiviruses on a molecular level. Here we present the first complete genome sequence analysis of two tick-transmitted human-pathogenic Orbiviruses. Evolutionary data based on these genome sequences indicate that during the development of Orbiviruses a host switch from insect-borne Orbiviruses toward tick-borne viruses occurred.

69

nucleotides (A). KEMV has 6 conserved bases in the 5′ NCRs of the genome segments and 3 conserved nucleotides in the 3′ NCRs (B). The first and last three nucleotides of all TRBV and KEMV segments are inverted complements. As expected, Proteinortho (Lechner et al., 2011) detected no protein homologues within a specific virus (relaxed settings). We found a significant conservation for VP1, VP3(T2), VP4(CaP), VP5, VP6 (Hel) and VP7(T13). The NS2(ViP) group missed the SCRV homolog and the NS3 group missed SCRV and YUOV proteins. Surprisingly, NS1(TuP) was clearly split into two high scoring clusters which were GIV, KEMV, TRBV, PHSV and YUOV on the one, and AHSV, EHDV, PALV, BTV on the other hand. For VP2 only one homologous group was detected (GIV, KEMV, TRBV, PHSV and PALV). Having a closer look (strict settings), NS1(TuP), NS2(ViP), NS3 and VP2 were split into two characteristic clusters (see Supplementary Table S4). Segment 9 of TRBV and KEMV contains a 27 bp and a 12 bp deletion respectively compared to the sequence of GIV (Fig. 1A). In reading frame +1 relative to the VP6(Hel) reading frame of TRBV, three additional VP6(Hel) overlapping ORFs can be identified via code2aln (Fig. 1B, Stocsits et al., 2005). ORF-Xa, -Xb and -Xc have a length of 288, 231 and 153 nucleotides (nt) possibly coding for proteins of 95, 76 and 50 amino acids (aa), with putative calculated molecular masses of 11.2, 9.1 and 6.1 kDa respectively. Similarly, in KEMV two additional VP6(Hel) overlapping coding sequences, ORF-Xa and -Xb (105 and 456 nt) possibly coding for two proteins (34 and 151 aa, 4.2 and 17.6 kDa) can be identified (Fig. 1B). At 73.3% TRBV ORF-Xa shows the highest homology to KEMV ORF-Xa. TRBV ORF-Xb shares the highest sequence similarity (69.7%) with KEMV ORF-Xb. As described a higher conservation of ORF-X than VP6(Hel) at aa level (Firth, 2008) was observed (Supplementary Table S3). Recently a VP6(Hel) overlapping ORF-2 which encodes a VP6(dBP) protein of 190 aa (22.5 kDa) with a predicted dsRNA binding domain has been identified in segment 9 of GIV (Belhouchet et al., 2010). The GIV VP6(dBP) N terminus shares 55.8% similarity with TRBV ORF-Xa and 58.6% with KEMV ORF-Xa, and the C terminus shares 52.6% similarity with TRBV ORF-Xb when analyzed with the ClustalV multiple alignment algorithm (Fig. 1C). The remaining sequence of GIV VP6(dBP) shares 52.0% similarity with KEMV ORF-Xb (data not shown). Phylogenetic analysis

Results Using the FLAC-method developed for sequencing of dsRNA viruses (Maan et al., 2007), we sequenced segments 3, 4, 6, 7, 8 and 9 as well as parts of segment 2 (around 55%) and segment 1 (around 30%) but no sequences were obtained for segments 5 and 10 of TRBV. The complete genomes of TRBV and of KEMV were then sequenced via pyrosequencing. With this approach we were able to obtain almost full-length sequences of the 10 genome segments of TRBV and KEMV in 1 lane of a 4-lane picotiter plate in a pool of 6 different RL-MID-tagged virus libraries in one 454 titanium sequencing run with 118-fold (TRBV) to 314-fold (KEMV) coverage. For some genome segments the very distal 5′ and 3′ terminal sequences could not be determined. To close these gaps we performed PCR's with an internal primer and a primer complementary to the FLAC-anchor primer (5-15-1) and subsequent cloning and sequencing. General genome comparison Supplementary Table S1 shows the length of the 10 TRBV and KEMV genome segments and the calculated theoretical molecular mass of the deduced encoded viral proteins. In KEMV NS1(TuP) is encoded by segment 5 and VP2 by segment 4, whereas in TRBV its vice versa. Supplementary Table S2 lists the 5′ and 3′ non-coding sequences of the genome segments of TRBV and KEMV. The 5′ NCRs of TRBV contain 5 conserved nucleotides and the 3′ NCRs 4 conserved

Multiple alignment of the putative protein sequences of TRBV and KEMV with all available Orbivirus protein sequences obtained from GenBank (Benson et al., 2002) or Swiss-Prot (Boutet et al., 2007) showed a high degree of similarity (up to 83%) to sequences of GIV, BRDV and SBaV, all members of the Great Island virus serogroup (Supplementary Table S3). Only the VP2 protein and the VP6(Hel) protein share lower than 50% similarity with the corresponding proteins from GIV. TRBV shares 64%–91% sequence similarity with KEMV and even 96%–100% identity to LIPV in the so far available protein sequences, most probably indicating variants of the same serotype (strain). For phylogenetic analysis of VP2 (Fig. 2B) and T2 proteins (Fig. 3A), TRBV- and KEMV-specific amino acid sequences were aligned to corresponding sequences from GenBank that represented all Orbivirus species via the ClustalW algorithm and NJ-trees were generated with MEGA version 4 (Tamura et al., 2007). In addition, we calculated protein-based NJ- and ML-trees for each protein (see Supplementary Table S6) and for the whole species based on concatenated ClustalW-protein-alignments (Fig. 3B). ML-trees were tested for the best model describing the single protein phylogeny. Surprisingly, nearly all of the single protein trees confirmed the species tree, independently of the eight methods we used. However, some small variants and unclear results should be mentioned. As expected, VP2 was unstable (bootstrap values b30%). Multifurcations for ML-trees of VP6(Hel), VP2 and NS1(TuP) were discovered. NS3 showed for both applications poor quality (bootstrap values b45%).

70

M. Dilcher et al. / Virology 423 (2012) 68–76

Fig. 1. A) Partial (bases 61–240) ClustalW sequence alignment of the nucleotide sequences of GIV segment 9, TRBV segment 9 and KEMV segment 9. Nucleotides that match GIV segment 9 are shaded in gray. A 27 bp deletion in the TRBV segment 9 sequence compared to GIV segment 9 is indicated with a dashed box and a 12 bp deletion in the KEMV segment 9 compared to GIV segment 9 is indicated with a solid rectangle. B) Genome maps of TRBV segment 9, KEMV segment 9 and GIV segment 9 with VP6(Hel), VP6(dBP) and putative overlapping ORF-X sequences as identified via the code2aln software (Stocsits et al., 2005). C) ClustalV sequence alignment of the amino acid sequences of GIV VP6 (dBP), TRBV ORF-Xa and TRBV ORF-Xb. Identical amino acids are shaded in gray.

Co-evolution of viruses and hosts To study the evolutionary dependency between Orbiviruses and their intermediate hosts a co-phylogenetic analysis between the trees of both species was performed using CoRe-PA (Merkle et al., 2010). This revealed a very good overlap between virus and vector evolution (Fig. 4A). Using the adaptive cost evaluation approach of CoRe-PA four most parsimonious reconstructions were found each being optimal

with respect to a certain event cost model. We thereby assumed that all viruses had one ancestor with the same intermediate host. To rate the reconstructions CoRe-PA uses a measurement, which describes how good the cost model fit to the reconstruction. The preferred reconstruction with 5 cospeciations, 4 sortings, 4 duplications and 1 host switching event (Fig. 4B) was rated 0.0032 while the other solutions had much higher values (0.0171, 0.0909 and 0.4). Furthermore the same reconstruction was rated best by using the event costs model

M. Dilcher et al. / Virology 423 (2012) 68–76

71

Fig. 2. A) Neighbor-joining phylogenetic tree of the T2 proteins of different Orbivirus species using the MEGA4 software (Tamura et al., 2007) rooted to SCRV and adapted from Belhouchet et al. (2010). Orbiviruses encoding T2 protein on segment 3 (Culicoides-borne viruses) form a separate group from those encoding T2 protein on segment 2 (mosquito-borne and tick-borne). KEMV, TRBV and LIPV form a distinct clade to GIV and BRDV. Bootstrap-values (1000 replicates) are indicated at the nodes. The scale indicates number of substitutions per site. B) NJ- and ML-total Orbivirus species trees (ten joined protein alignments). NJ-trees with 1000 bootstraps and ML-trees with 100 bootstraps. The scale indicates number of substitutions per site. The dashed rectangles indicate close relationship of mosquito-borne PHSV and YUOV with tick-transmitted Orbiviruses GIV, TRBV and KEMV.

proposed in (Charleston, 1998). The observed host switch occurred in the ancestral virus of KEMV/TRBV/GIV from Culex toward Ixodes hosts, marked green in Fig. 4B. For the co-phylogenetic system a significance test with 1000 randomized host-parasite associations resulted in 99.4% of the solutions having less than 5 cospeciations and none having more (see Supplementary Figure S1). This indicates a significant amount of congruence between the two trees as well as a strong dependency between Orbiviruses and their intermediate host in comparison to what would be expected by chance. Furthermore the randomized

instances were rated significantly worse (see Supplementary Figure S2). To analyze the stability of the observed host switch we computed co-phylogenetic reconstructions with the 22 protein phylogenies as alternative viral histories. In all these scenarios the most parsimonious reconstructions contain a switch from either the Culex (8 of 22), Culicoides (6 of 22) or Culex/Culicoides (8 of 22) branch toward Ixodes (Fig. 4B). In addition a detailed analysis of each single protein phylogeny compared to the species tree was studied. Therefore the cophylogenetic history for each of the 22 protein phylogenies with

72

M. Dilcher et al. / Virology 423 (2012) 68–76

A Accession number TRBV VP2 KEMV VP2 BRDV VP4

HQ266585 HQ266594 -

GIV VP4 SCRV VP3 BTV VP2 AHSV VP2

ADM88596 YP_052944 ACR58459 ACH92678

VP2 homologous protein (kDa) 62 63 63 (Schoehn et al., 1997) 62 74 112 123

VP2 homologous protein (AA) 554 554 n.d. 551 654 961 1056

B

Fig. 3. A) Comparison of the sizes of VP2 homologous proteins of tick-borne (light gray) and insect-borne (dark gray) Orbiviruses. B) Neighbor-joining phylogenetic tree of the VP2 homologous proteins of different Orbivirus species generated with the MEGA4 software (Tamura et al., 2007). KEMV and TRBV form a distinct clade to GIV. Bootstrap-values (1000 replicates) are indicated at the nodes. The scale indicates number of substitutions per site.

respect to the viral species tree was reconstructed. All discovered candidates of horizontal gene transfer between the viruses can be found in Supplementary Table S5. In general two different transfers were found in several proteins. One between EEV/AHSV/PALV/EHDV/BTV and PHSV/YUOV (VP5 and NS2(ViP)) and the other between AHSV/ PALV and EHDV/BTV (NS1(TuP), NS2(ViP) and NS3). Discussion Genome analysis of TRBV and KEMV confirms serogroup and serotypes In general, Orbivirus serogroups are defined based on similarities of the proteins T2 and VP7(T13), which form the stable inner-shell or core, while classification of Orbiviruses into serotypes is mainly based on similarities of the outer-shell proteins VP2 and VP5, which interact with neutralizing antibodies. The serogrouping protein T2 (designated VP3(T2) or VP2(T2)) is highly conserved and exhibits T = 2 symmetry (Gouet et al., 1999; Grimes et al., 1998). Different members of an Orbivirus serogroup usually exhibit levels of T2 amino acid sequence identity higher than 91%

(Attoui et al., 2001). The high VP3(T2) homology of TRBV and LIPV (99.6%) leaves KEMV in the serogroup by a tight fit (91.2%) but clearly sets apart GIV and BRDV (82.7%, 76.8%). Therefore, in agreement with observations from Belhouchet et al. (2010) TRBV, LIPV and KEMV form one distinct serogroup apart from the Great Island virus serogroup (GIV, BRDV, Fig. 2A). Based on these amino acid sequence homologies it is possible to group Orbivirus T2 proteins into two main groups, T2 proteins encoded by segment 2 (mosquito- and tick-borne) and T2 proteins encoded by segment 3 (Culicoides-borne) (Attoui et al., 2005, 2009; Belhouchet et al., 2010). The T2 protein of SCRV groups apart from the others, which supports the observation that SCRV is the most divergent Orbivirus (Belhouchet et al., 2010, Fig. 2A). NJand ML-trees based on concatenated ClustalW protein alignments indicate that, as in the case of the T2 protein, mosquito-borne PHSV and YUOV are closer related to the tick-transmitted Orbiviruses GIV, TRBV and KEMV than the Culicoides-borne Orbiviruses (Fig. 2B). VP7(T13), which has T = 13 symmetry is mostly encoded by segment 7 (BTV, GIV, BRDV and SBaV). It is a component of the inner shell, but accessible for antibodies and is therefore used as the major serogroup specific antigen (Fauquet et al., 2005; Schoehn

M. Dilcher et al. / Virology 423 (2012) 68–76

73

Fig. 4. Co-phylogenetic reconstruction of the relationship between Orbiviruses and their corresponding vectors on the basis of concatenated ClustalW-protein-alignments and taxonomic relationships of host organisms based on the phylogenetic protein marker Cox1 using the CoRe-PA program (Merkle et al., 2010). A) Host parasite system of insect hosts and Orbiviruses. B) Co-phylogenetic reconstruction of insect hosts and Orbiviruses with horizontal transfer of viral proteins (CoRe-PA rating: 0.0032).

et al., 1997). TRBV VP7(T13) is encoded by segment 8 whereas segment 7 codes for NS2(ViP). The same pattern can be seen in KEMV and in YUOV (Attoui et al., 2005). TRBV VP7(T13) and KEMV VP7(T13) share 86.3% amino acid sequence similarity. The similarity to GIV and BRDV is with 77.6% and 64.8% significantly lower. The highly variable outer shell serotyping protein VP2 is involved in the attachment of the virus to the host cell surface. In some Orbiviruses the VP2 homologous protein is called VP3 and encoded by segment 3 (YUOV), in others it's called VP4 and encoded by segment 5 (GIV, Supplementary Table S4). So far no VP2 sequence information is available for LIPV. In the case of BTV it was shown that VP2 interacts with neutralizing antibodies and is therefore subjected to antibody selective pressure. Hence, protein sequence variability of VP2 correlates with virus serotype (Mertens et al., 2000). VP2 proteins of tick-transmitted Orbiviruses have only half the size of VP2 proteins

of insect-borne Orbiviruses (Fig. 3A, Schoehn et al., 1997), which indicate that they might posses different mechanisms in cell-attachment and are therefore not absolutely comparable. Interestingly, the VP2 homologous protein VP4 of GIV shows only a low homology of 37.5% to TRBV VP2 (Supplementary Table S3). This difference might give a hint toward the molecular basis of pathogenicity, since the VP2 protein is assumed to interact with the vertebrate cell receptor and GIV infects seabirds, whereas TRBV and KEMV appear to infect rodents and man. The phylogenetic analysis based on VP2 homology (Fig. 3B) therefore supports classification of KEMV and TRBV into a group different from GIV. Additional reading frames in segment 9 The strong selective pressure on the small genomes of RNA viruses leads to compact coding strategies such as overlapping CDSs to

74

M. Dilcher et al. / Virology 423 (2012) 68–76

optimize information content. The VP6(Hel) overlapping ORFs in TRBV and KEMV segment 9 show higher homologies to each other and to the respective regions in the recently discovered VP6(Hel) overlapping ORF-2 encoding VP6(dBP) in GIV than the corresponding VP6(Hel) proteins. Therefore it is very likely that they possess important functions. VP6(Hel) overlapping open reading frames (ORF-X) are also known from other Orbiviruses like BTV, AHSV, PALV, PHSV and YUOV, based on bioinformatic analysis of the segment 9 sequences (Firth, 2008). Interestingly, TRBV and KEMV ORF-Xa seems to be homologous to the N-terminal part of the GIV VP6(dBP) and the ORV-Xb's to the C-terminal part. So far, no functional role could be attributed to ORF-X protein products, but for GIV VP6(dBP) a role in counteracting innate immunity like σ3 of the mammalian Orthoreoviruses has been discussed (Belhouchet et al., 2010). It remains to be seen what roles the ORF-X products may have. Co-evolution of viruses and hosts The co-phylogenetic analysis of Orbiviruses and their vectors indicated that the KEMV/TRBV/GIV ancestor interacted with the ancestor of a Culex species and then suddenly changed its host to some Ixodes species. The stability analysis done with the alternative Orbivirus phylogenies supports the landing host of this switch being Ixodes, but there is an uncertainty about the take-off host being either Culex, Culicoides or the ancestor of both. However, the direction of the switch is clearly toward Ixodes hosts. This is in conflict to the observation of Belhouchet et al. (2010) who describe a host switch from an ancestral tick-borne Orbivirus to a mosquito-borne Orbivirus and that Culicoides-borne viruses were the last to appear. Through the analysis of each single protein horizontal gene transfers could be identified between EEV/AHSV/PALV/EHDV/BTV and PHSV/YUOV (VP5 and NS2 (ViP)) and between AHSV/PALV and EHDV/BTV (NS1(TuP), NS2 (ViP) and NS3). In both transfers the different ancestral viruses shared the same or at least close related intermediate hosts, indicating a chance of exchanging genetic material. Furthermore, this assumption is supported by the lower bootstrap values of the viral species tree at the respective positions. Because Orbiviruses of the Great Island virus serogroup seem to share high sequence similarity with those of the Kemerovo serogroup, but don't seem to cause disease in humans or animals, further studies should be performed to identify the molecular basis for this difference in pathology and virulence between these Orbivirus serogroups. Materials and methods Growth and purification of Tribeč virus and Kemerovo virus TRBV was passaged 4 times in Vero B4 cells maintained in MEM medium and KEMV was passaged 4 times in Vero E6 cells maintained in DMEM medium, both supplemented with 2% fetal bovine serum, 2 mM Glutamine, 10 mM Penicillin and 10 mM Streptomycin in 175 cm 2 tissue culture flasks. Culture supernatants of infected cells were collected when the cells showed 90–100% CPE (app. 1–2 dpi), centrifuged at 700 ×g (2000 rpm) for 10 min then at 2800 ×g (4000 rpm) for 5 min and filtered through a 0.2 μm sterile filter. Twenty milliliters of supernatant was mixed with 1.48 ml 5 M NaCl and 10.8 ml 30% PEG8000 in NTE (10 mM Tris, pH6.5; 1 mM EDTA; 100 mM NaCl), incubated on a shaker for 30 min at 4 °C and subsequently centrifuged at 6000 rpm for 60 min and 4 °C. RNA extraction The dsRNA of TRBV was extracted from virus pellet using a guanidinium isothiocyanate (Trizol) procedure according to instructions

(peqGOLD Trifast, Peqlab, Germany). Each RNA pellet was resuspended in 25 μl RNase-free DEPC-treated water and combined into one tube. The quality of the dsRNA was assessed by 1% Agarose gel electrophoresis (AGE) in TAE buffer containing 0.5 μg/ml of ethidium bromide. Anchor-primer ligation For the amplification of full-length cDNAs the FLAC-method developed by Maan et al., 2007 (Maan et al., 2007) was applied, where a self-complementary ‘anchor-primer’ was ligated to the 3′ ends of the dsRNA genome segments. The ‘anchor-primer’ was ligated to 10 μg of unfractionated viral dsRNA. Electrophoretic separation of ligated RNAs After ligation to the ‘anchor-primer’, the dsRNA genome segments were separated by AGE at app. 5 V/cm in a long chamber with buffer circularisation at 4 °C for app. 7 h. Individual segments or mixtures of segments were recovered from the gel using silica binding methods (RNaid kit, MP Biomedicals, Germany). cDNA first strand synthesis cDNA synthesis was carried out using the Transcriptor High Fidelity cDNA Synthesis Kit from Roche Diagnostics, Germany. 15 to 30 ng of purified adaptor-ligated dsRNA was used. The ‘anchor-primer’ served as a self-priming adaptor primer for the first-strand cDNA synthesis. The cDNA was used directly for PCR amplification or stored at −20 °C. PCR amplification of cDNAs Amplification of cDNAs was performed using primer 5-15-1 (Maan et al., 2007) and PfuUltra-II Fusion HS DNA Polymerase (Stratagene/Agilent, Germany). Amplification was carried out after denaturation at 95 °C for 3 min, followed by 40 cycles of denaturation for 30 s at 95 °C, annealing for 30 s at 58 °C, extension for 20 s to 2.5 min (depending on the length of the amplicon) at 72 °C and a final extension step for 3 min at 72 °C. Cloning and sequencing Amplified genome segments were gel purified (Zymoclean Gel DNA Recovery Kit, Zymo Research Europe, Germany), followed by A-tailing and cloning into pCRII-vectors. Sequencing of plasmids was carried-out with M13-forward and -reverse primer via cycle sequencing (Seqlab, Germany). Pyrosequencing For pyrosequencing of TRBV and KEMV, approx. 20 ml cell culture supernatant from infected vero cells (1 dpi) was first centrifuged for 10 min at 2000 rpm and then for 5 min at 4000 rpm and filtered through 0.2 μm filters, followed by PEG-precipitation (see above) to enrich for virus particles. Subsequently the virus pellet was resuspended in 500 μl PBS. Viral RNA was extracted as described above and treated with DNase to digest cellular DNA's (Turbo DNA free Kit, Ambion, Applied Biosystems, Germany). Double-stranded viral RNA was purified from contaminating single-stranded host RNA by overnight precipitation in 2 M lithium chloride (LiCl) at 4 °C as described elsewhere (Attoui et al., 2000). The quality of the dsRNA was assessed by AGE. To be able to cover the 5′ and 3′ ends of the dsRNA segments we ligated the same anchor-primer sequence used in the FLAC-method to the 3′ ends of the viral RNAs prior to pyrosequencing as

M. Dilcher et al. / Virology 423 (2012) 68–76

described by Maan et al. (2007) and Potgieter et al. (2009). 500 ng purified viral dsRNA was ligated to 500 ng FLAC-anchor primer. The concentration of the adapter-ligated dsRNA was determined via Nanodrop and Qant-iT RiboGreen Assay (Invitrogen, Germany). 60 ng adapter-ligated viral RNA was subsequently amplified and converted to dsDNA with the help of the TransPlex Whole Transciptome Amplification kit (WTA2) from Sigma-Aldrich. After purification with the QIAquick PCR Purification kit (Qiagen, Germany) an additional size exclusion step via Ampure-XP beads (Agencourt, Germany) in a ratio of 1 vol WTA2 product to 0.7 vol Ampure-XP beads was used to remove fragments shorter than 350 bp. 300 ng of the whole genome amplified dsDNA was used for Titanium Shotgun Rapid Library Preparation and pyrosequencing as described in the FLX Titanium Protocol (Roche Applied Sciences, Germany) and in (Margulies et al., 2005), but omitting the DNA fragmentation by nebulization step. In the RL adapter ligation step, RL-MID adapters were used to allow pooling of several samples. Bioinformatics Assembly Assembly of the sequenced TRBV and KEMV genome segments was done with the Genome Sequencer FLX System Software Package version 2.3 (GS De Novo Assembler, GS Reference Mapper) in combination with the commercially available SeqMan Pro Software version 7.2.2 (DNASTAR, Lasergene).

75

overlapping ORF sequences were calculated via the code2aln software (Stocsits et al., 2005) with standard parameters. Detected ORFs were cut and cat in succession. Alignment of GIV VP6(dBP) with TRBV ORF-Xa and TRBV ORF-Xb was done with MegAlign version 7.2.2. using the ClustalV algorithm. Co-evolution of viruses and hosts The obtained phylogenetic ML- and NJ-trees of viruses and their intermediate hosts were used for co-phylogenetic analysis. This calculation was performed using the event-based model described in (Charleston, 1998) with co-speciation, duplication, sorting and host switching events. Evaluation was done with CoRe-PA (Merkle et al., 2010) using 10,000 different event cost models. For significance evaluation a randomization test with 1000 instances with randomized host–parasite associations was performed as well as an analysis with the 22 protein phylogenies as an alternative evolutionary history for the Orbiviruses. Additionally a co-phylogenetic analysis of viral proteins compared to the species tree was performed with a modified version of CoRe-PA. We considered only cospeciation, sorting and host switching events, since duplication of protein genes is not expected. An observed “host switch” is interpreted as horizontal gene transfer and a sorting event in this context represents a loss induced by a replacement with the protein of a different virus. We considered all possible reconstructions with maximal number of cospeciations. Acknowledgments

General genome comparison Additional virus genomes were downloaded from GenBank (Benson et al., 2002) and Swiss-Prot (Boutet et al., 2007) and compared to TRBV and KEMV with ClustalW (Thompson et al., 1994) based on protein sequences (all GenBank and Swiss-Prot accession numbers for the used proteins can be found in Supplementary Table S3). In case of partial sequences the shorter sequence was used as reference for percent identity calculation. To test homologous relationships of the single proteins we used Proteinortho v.4 (Lechner et al., 2011) with the settings –selfblast – pairs –cov = 0.5 –e = 1e-10 (strict settings) and with – selfblast – pairs –cov = 0.0 –e = 1e-02 (relaxed settings). Phylogenetic analysis For a convincing phylogenetic relation of the viruses we calculated ML(maximum likelihood)- and NJ(neighbor-joining)-trees (with and without gaps) based on nucleotide and protein sequences of (a) the whole genome (concatenated ClustalW-protein-alignments with and without 5′ and 3′ UTR) (b) each of the ten single proteins and (c) ORF-X overlapping VP6(Hel). NJ-trees were build by ClustalW – tree –bootstrap = 1000. For calculation of the best model of the virus' ML-trees we used ProtTest (Abascal et al., 2005) based on PhyML (Guindon and Gascuel, 2003) with 100 bootstraps. Most proteins were adopted to an LG (Le and Gascuel, 2008) and + I+G+F model, whereas the divergent VP2 was applied to the WAG (Whelan and Goldman, 2001) +I+G+F model and the shifted putative ORF-X protein belonged to an retrovirus model (rtREV (Dimmic et al., 2002) +I+G+F). The obtained model +I+G+F has an invariable fraction of amino acids (Reeves, 1992), different categories of evolution rates (Yang, 1993) and corresponding probabilities (Cao et al., 1998). For calculation of the phylogenetic host tree we used the phylogenetic protein marker Cox1 (Evans et al., 2007; Lin et al., 2009) downloaded from vectorbase (Lawson et al., 2009). In addition, for phylogenetic analysis of the T2 and VP2 proteins, TRBV- and KEMV-specific sequences were aligned to a selection of available corresponding sequences from GenBank that represented all Orbivirus species by using the ClustalW algorithm of commercially available MegAlign software version 7.2.2 (DNAStar, Lasergene) and MEGA version 4 (Tamura et al., 2007). Genome maps and putative

This work was supported by the Federal Ministry of Education and Research (BMBF), grant number 01KI0710, “Research on Zoonotic Infectious Diseases” program, “Emerging arthropode-borne viral infections in Germany: Pathogenesis, diagnostics and surveillance” and the BMBF funded research program “Potential release-oriented biothreat emergency diagnostics (P.R.O.B.E)” for civil security of the German Federal Government as part of the high-tech strategy for Germany. In addition, this work was supported by the German Research Foundation (DFG) through the project “Deep Metazoan Phylogeny” within SPP 1174. Appendix A. Supplementary material Supplementary material to this article can be found online at doi:10.1016/j.virol.2011.11.020. References Abascal, F., Zardoya, R., Posada, D., 2005. ProtTest: selection of best-fit models of protein evolution. Bioinformatics 21 (9), 2104–2105. Attoui, H., Billoir, F., Cantaloube, J.F., Biagini, P., de Micco, P., de Lamballerie, X., 2000. Strategies for the sequence determination of viral dsRNA genomes. J. Virol. Methods 89 (1–2), 147–158. Attoui, H., Stirling, J.M., Munderloh, U.G., Billoir, F., Brookes, S.M., Burroughs, J.N., de Micco, P., Mertens, P.P.C., de Lamballerie, X., 2001. Complete sequence characterization of the genome of the St Croix River virus, a new orbivirus isolated from cells of Ixodes scapularis. J. Gen. Virol. 82, 795–804. Attoui, H., Jaafar, F.M., Belhouchet, M., Aldrovandi, N., Tao, S.J., Chen, B.Q., Liang, G.D., Tesh, R.B., de Micco, P., de Lamballerie, X., 2005. Yunnan orbivirus, a new orbivirus species isolated from Culex tritaeniorhynchus mosquitoes in China. J. Gen. Virol. 86, 3409–3417. Attoui, H., Mendez-Lopez, M.R., Rao, S.J., Hurtado-Alendes, A., Lizaraso-Caparo, F., Jaafar, F.M., Samuel, A.R., Belhouchet, M., Pritchard, L.I., Melville, L., Weir, R.P., Hyatt, A.D., Davis, S.S., Lunt, R., Calisher, C.H., Tesh, R.B., Fujita, R., Mertens, P.P.C., 2009. Peruvian horse sickness virus and Yunnan orbivirus, isolated from vertebrates and mosquitoes in Peru and Australia. Virology 394 (2), 298–310. Belhouchet, M., Jaafar, F.M., Tesh, R.B., Grimes, J., Maan, S., Mertens, P., Attoui, H., 2010. Complete sequence of the Great Island virus and comparison with the T2 and outer-capsid proteins of Kemerovo, Lipovnik and Tribec viruses (genus Orbivirus, family Reoviridae). J. Gen. Virol. 91 (12), 2985–2993. Benson, D.A., Karsch-Mizrachi, I., Lipman, D.J., Ostell, J., Rapp, B.A., Wheeler, D.L., 2002. GenBank. Nucleic Acids Res. 30 (1), 17–20. Boutet, E., Lieberherr, D., Tognolli, M., Schneider, M., Bairoch, A., 2007. UniProtKB/ Swiss-Prot. Methods Mol. Biol. 406, 89–112.

76

M. Dilcher et al. / Virology 423 (2012) 68–76

Cao, Y., Janke, A., Waddell, P.J., Westerman, M., Takenaka, O., Murata, S., Okada, N., Paabo, S., Hasegawa, M., 1998. Conflict among individual mitochondrial proteins in resolving the phylogeny of eutherian orders. J. Mol. Evol. 47 (3), 307–322. Charleston, M.A., 1998. Jungles: a new solution to the host/parasite phylogeny reconciliation problem. Math. Biosci. 149 (2), 191–223. Chumakov, M.P., Ernek, E., Kozuch, O., Sarmanov, E., Sergeeva, G.I., Tapupere, V.O., Bychkova, M.B., Karpovich, L.G., Libikova, H., Rehacek, J., Mayer, V., 1963. Report on Isolation from Ixodes Persulcatus ticks and from patients in Western Siberia of a virus differing from agent of tick-borne encephalitis. Acta Virol. 7 (1), 82–83. Dimmic, M.W., Rest, J.S., Mindell, D.P., Goldstein, R.A., 2002. rtREV: an amino acid substitution matrix for inference of retrovirus and reverse transcriptase phylogeny. J. Mol. Evol. 55 (1), 65–73. Dobler, G., 1996. Arboviruses causing neurological disorders in the central nervous system. Arch. Virol. 33–40. Dobler, G., Wolfel, R., Schmuser, H., Essbauer, S., Pfeffer, M., 2006. Seroprevalence of tick-borne and mosquito-borne arboviruses in European brown hares in Northern and Western Germany. Int. J. Med. Microbiol. 296, 80–83. Doherty, R.L., Carley, J.G., Murray, M.D., Main, A.J., Kay, B.H., Domrow, R., 1975. Isolation of Arboviruses (Kemerovo-Group, Sakhalin-Group) from Ixodes-Uriae Collected at Macquarie Island, Southern Ocean. Am.J.Trop. Med. Hyg. 24 (3), 521–526. Evans, K.M., Wortley, A.H., Mann, D.G., 2007. An assessment of potential diatom “barcode” genes (cox1, rbcL, 18S and ITS rDNA) and their effectiveness in determining relationships in Sellaphora (Bacillariophyta). Protist 158 (3), 349–364. Fauquet, C.M., Mayo, M.A., Maniloff, J., Desselberger, U., Ball, L.A., 2005. Virus Taxonomy. Classification and Nomenclature of Viruses. Eighth Report of the International Committee on the Taxonomy of Viruses. Elsevier Academic Press, San Diego, USA VIII. Firth, A.E., 2008. Bioinformatic analysis suggests that the Orbivirus VP6 cistron encodes an overlapping gene. Virol. J. 5 (48). Gorman, B.M., 1979. Variation in Orbiviruses. J. Gen. Virol. 44, 1–15 (Jul). Gorman, B.M., Taylor, J., Walker, P.J., 1983. Orbiviruses. In: Joklik, W.K. (Ed.), The Reoviridae. Plenum, New York, pp. 287–357. Gorman, B.M., Taylor, J., Morton, H.C., Melzer, A.J., Young, P.R., 1984. Characterization of nugget virus, a serotype of the Kemerovo group of Orbiviruses. Aust. J. Exp. Biol. Med. Sci. 62, 101–115 (Feb). Gouet, P., Diprose, J.M., Grimes, J.M., Malby, R., Burroughs, J.N., Zientara, S., Stuart, D.I., Mertens, P.P.C., 1999. The highly ordered double-stranded RNA genome of bluetongue virus revealed by crystallography. Cell 97 (4), 481–490. Gresikova, M., Nosek, J., Kozuch, O., Ernek, E., Lichard, M., 1965. Study on ecology of Tribec virus. Acta Virol. 9 (1), 83–88. Gresikova, M., Rajcani, J., Hruzik, J., 1966. Pathogenicity of Tribec virus for Macaca Rhesus monkeys and white mice. Acta Virol. 10 (5), 420–424. Grimes, J.M., Burroughs, J.N., Gouet, P., Diprose, J.M., Malby, R., Zientara, S., Mertens, P.P.C., Stuart, D.I., 1998. The atomic structure of the bluetongue virus core. Nature 395 (6701), 470–478. Guindon, S., Gascuel, O., 2003. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol. 52 (5), 696–704. Labuda, M., Nuttall, P.A., 2004. Tick-borne viruses. Parasitology 129, S221–S245. Lawson, D., Arensburger, P., Atkinson, P., Besansky, N.J., Bruggner, R.V., Butler, R., Campbell, K.S., Christophides, G.K., Christley, S., Dialynas, E., Hammond, M., Hill, C.A., Konopinski, N., Lobo, N.F., MacCallum, R.M., Madey, G., Megy, K., Meyer, J., Redmond, S., Severson, D.W., Stinson, E.O., Topalis, P., Birney, E., Gelbart, W.M., Kafatos, F.C., Louis, C., Collins, F.H., 2009. VectorBase: a data resource for invertebrate vector genomics. Nucleic Acids Res. 37, D583–D587. Le, S.Q., Gascuel, O., 2008. An improved general amino acid replacement matrix. Mol. Biol. Evol. 25 (7), 1307–1320. Lechner, M., Findeiss, S., Steiner, L., Marz, M., Stadler, P.F., Prohaska, S.J., 2011. ProteinOrtho: detection of (co-)orthologs in large-scale analysis. BMC Bioinforma. 12 (124). Libikova, H., Rehacek, J., Ernek, E., Gresikov, M., Somogyio, J., Kozuch, O., 1964. Cytopathic viruses isolated from Ixodes ricinus ticks in Czechoslovakia. Acta Virol. 8 (1), 96–106. Libikova, H., Rehacek, J., Somogyio, J., 1965. Viruses related to Kemerovo virus in Ixodes ricinus ticks in Czechoslovakia. Acta Virol. 9 (1), 76–82. Libikova, H., Tesarova, J., Rajcani, J., 1970. Experimental infection of monkeys with Kemerovo virus. Acta Virol. 14 (1), 64–69.

Libikova, H., Heinz, F., Ujhazyova, D., Stunzner, D., 1978. Orbiviruses of Kemerovo complex and neurological diseases. Med. Microbiol. Immunol. 166 (1–4), 255–263. Lin, S., Zhang, H., Hou, Y.B., Zhuang, Y.Y., Miranda, L., 2009. High-level diversity of dinoflagellates in the natural environment, revealed by assessment of mitochondrial cox1 and cob genes for dinoflagellate DNA barcoding. Appl. Environ. Microbiol. 75 (5), 1279–1290. Maan, S., Rao, S.J., Maan, N.S., Anthony, S.J., Attoui, H., Samuel, A.R., Paul, P., Mertens, C., 2007. Rapid cDNA synthesis and sequencing techniques for the genetic study of bluetongue and other dsRNA viruses. J. Virol. Methods 143 (2), 132–139. Major, L., Linn, M.L., Slade, R.W., Schroder, W.A., Hyatt, A.D., et al., 2009. Ticks associated with Macquarie Island penguins carry arboviruses from four genera. PLoS One 4 (2), e4375. Malkova, D., Holubova, J., Kolman, J.M., Marhoul, Z., Hanzal, F., Kulkova, H., Markvart, K., Simkova, L., 1980. Antibodies against some arboviruses in persons with various neuropathies. Acta Virol. 24 (4), 298. Margulies, M., Egholm, M., Altman, W.E., Attiya, S., Bader, J.S., Bemben, L.A., Berka, J., Braverman, M.S., Chen, Y.J., Chen, Z.T., Dewell, S.B., Du, L., Fierro, J.M., Gomes, X.V., Godwin, B.C., He, W., Helgesen, S., Ho, C.H., Irzyk, G.P., Jando, S.C., Alenquer, M.L.I., Jarvie, T.P., Jirage, K.B., Kim, J.B., Knight, J.R., Lanza, J.R., Leamon, J.H., Lefkowitz, S.M., Lei, M., Li, J., Lohman, K.L., Lu, H., Makhijani, V.B., McDade, K.E., McKenna, M.P., Myers, E.W., Nickerson, E., Nobile, J.R., Plant, R., Puc, B.P., Ronan, M.T., Roth, G.T., Sarkis, G.J., Simons, J.F., Simpson, J.W., Srinivasan, M., Tartaro, K.R., Tomasz, A., Vogt, K.A., Volkmer, G.A., Wang, S.H., Wang, Y., Weiner, M.P., Yu, P.G., Begley, R.F., Rothberg, J.M., 2005. Genome sequencing in microfabricated high-density picolitre reactors. Nature 437 (7057), 376–380. Merkle, D., Middendorf, M., Wieseke, N., 2010. A parameter-adaptive dynamic programming approach for inferring cophylogenies. BMC Bioinforma. 11 (Suppl. I), S60. Mertens, P.P.C., Arella, M., Attoui, H., Belloncik, S., Bergoin, M., Boccardo, G., Booth, T.F., Chiu, W., Diprose, J.M., Duncan, R., et al., 2000. Family Reoviridae. Virus Taxonomy. Seventh Report of the International Committee on Taxonomy of Viruses, pp. 395–480. Moss, S.R., Jones, L.D., Nuttall, P.A., 1992. Comparison of the major structural core proteins of Tick-borne and Culicoides-borne Orbiviruses. J. Gen. Virol. 73, 2585–2590. Nuttall, P.A., Moss, S.R., 1989. Genetic reassortment indicates a new grouping for Tickborne Orbiviruses. Virology 171 (1), 156–161. Potgieter, A.C., Page, N.A., Liebenberg, J., Wright, I.M., Landt, O., van Dijk, A.A., 2009. Improved strategies for sequence-independent amplification and sequencing of viral double-stranded RNA genomes. J. Gen. Virol. 90, 1423–1432. Reeves, J.H., 1992. Heterogeneity in the substitution process of amino-acid sites of proteins coded for by mitochondrial-DNA. J. Mol. Evol. 35 (1), 17–31. Schmidt, J.R., Shope, R.E., 1971. Kemerovo virus from a migrating common redstart of Eurasia. Acta Virol. 15 (1), 112. Schoehn, G., Moss, S.R., Nuttall, P.A., Hewat, E.A., 1997. Structure of Broadhaven virus by cryoelectron microscopy: correlation of structural and antigenic properties of Broadhaven virus and bluetongue virus outer capsid proteins. Virology 235 (2), 191–200. Skoferts, P., Yarovoi, P.I., Gaidamov, S., Korchmar, N.D., Obukhova, V.R., Melnikov, E., Klisenko, G.A., 1972. Isolation in Moldavian Ssr of a Kemerovo group arbovirus from Ixodes-ricinus ticks. Acta Virol. 16 (4), 362. Stocsits, R.R., Hofacker, I.L., Fried, C., Stadler, P.F., 2005. Multiple sequence alignments of partially coding nucleic acid sequences. BMC Bioinforma. 6 (160). Süss, J., Schrader, C., 2004. Tick-borne human pathogenic microorganisms found in Europe and those considered non-pathogenig Part I: ticks and viruses. Bundesgesundheitslatt 47, 392–404. Tamura, K., Dudley, J., Nei, M., Kumar, S., 2007. MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0. Mol. Biol. Evol. 24 (8), 1596–1599. Thompson, J.D., Higgins, D.G., Gibson, T.J., 1994. Clustal-W — improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22 (22), 4673–4680. Whelan, S., Goldman, N., 2001. A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach. Mol. Biol. Evol. 18 (5), 691–699. Yang, Z.H., 1993. Maximum-likelihood-estimation of phylogeny from DNA-sequences when substitution rates differ over sites. Mol. Biol. Evol. 10 (6), 1396–1401.