ABSTRACT
Lyssaviruses are unsegmented RNA viruses causing rabies. Their vectors belong to the Carnivora and Chiroptera orders. We studied 36 carnivoran and 17 chiropteran lyssaviruses representing the main genotypes and variants. We compared their genes encoding the surface glycoprotein, which is responsible for receptor recognition and membrane fusion. The glycoprotein is the main protecting antigen and bears virulence determinants. Point mutation is the main force in lyssavirus evolution, as Sawyer's test and phylogenetic analysis showed no evidence of recombination. Tests of neutrality indicated a neutral model of evolution, also supported by globally high ratios of synonymous substitutions (dS ) to nonsynonymous substitutions (dN ) (>7). Relative-rate tests suggested similar rates of evolution for all lyssavirus lineages. Therefore, the absence of recombination and similar evolutionary rates make phylogeny-based conclusions reliable. Phylogenetic reconstruction strongly supported the hypothesis that host switching occurred in the history of lyssaviruses. Indeed, lyssaviruses evolved in chiropters long before the emergence of carnivoran rabies, very likely following spillovers from bats. Using dated isolates, the average rate of evolution was estimated to be roughly 4.3 × 10−4 dS /site/year. Consequently, the emergence of carnivoran rabies from chiropteran lyssaviruses was determined to have occurred 888 to 1,459 years ago. Glycoprotein segments accumulating more dN than dS were distinctly detected in carnivoran and chiropteran lyssaviruses. They may have contributed to the adaptation of the virus to the two distinct mammal orders. In carnivoran lyssaviruses they overlapped the main antigenic sites, II and III, whereas in chiropteran lyssaviruses they were located in regions of unknown functions.
Emergence of viral diseases is a permanent threat for animal and public health and may happen whenever environmental conditions are propitious. Virus spillover to new hosts may lead to an emerging disease, provided the virus gains a sufficient fitness. RNA viruses (riboviruses and retroviruses) having a polymerase devoid of a proofreading mechanism are the fastest-evolving organisms (15). They produce a diverse viral population, i.e., quasispecies (14), ready to explore new conditions or escape defense systems. This property makes RNA viruses among the most dangerous pathogens. Indeed, RNA viruses cause two of the six leading infectious killers, i.e., AIDS and measles, and are implicated in two others, i.e., acute respiratory infections and diarrheal diseases (21a). Therefore, understanding their evolution and particularly their emergence may help in fighting them.
The Lyssavirus genus belongs to the Rhabdoviridaefamily of the Mononegavirales order and includes unsegmented RNA viruses causing rabies encephalomyelitis. They are well fitted to vectors belonging mainly to the orders Carnivora andChiroptera. Seven genotypes (GTs) have so far been delineated within the genus (8, 21). These GTs were divided into two immunopathologically and genetically distinct phylogroups (3). Phylogroup I includes two African GTs:Mokola virus, which has been isolated from shrews and cats, although its reservoir remains unknown, and Lagos bat virus, which has been found mainly in frugivorous bats but also in an insectivorous bat. Phylogroup II has five GTs: Duvenhage virus (Africa), European bat lyssavirus 1 (EBLV-1; Europe), EBLV-2 (Europe), Australian bat lyssavirus(Australia), and the classical Rabies virus (RABV; worldwide). Members of the GTs Duvenhage virus, EBLV-1, and EBLV-2 are exclusively found in insectivorous bats, members of the GTAustralian bat lyssavirus are found in both insectivorous and frugivorous bats, and members of the GT RABV are found in carnivores and American bats (insectivorous, frugivorous, and hematophagous) (28). The fact that lyssaviruses are well established in two ecologically distinct mammal orders may very likely be a consequence of successful host switching. Therefore, to phylogenetically investigate the possibility of this host switching during the evolution of the Lyssavirus genus, we assessed the evolutionary forces, rate, and model. We also searched for regions thought to be responsible for host adaptation.
The external viral glycoprotein (G) was appropriate for this study particularly because of its host adaptation potential. The mature G protein (without a signal peptide [SP]) forms a trimer (20) and has an endodomain (ENDO) that interacts with internal proteins (11, 30, 52), a transmembrane (TM) region, and an ectodomain (ECTO) protruding from the viral membrane. The ECTO carries the B- and T-cell antigenic sites (6, 26) and the regions responsible for receptor recognition (27, 48, 50, 51) and membrane fusion (16). It also bears residues important for virulence (9, 12, 33, 34).
MATERIALS AND METHODS
Viruses.Fifty-five isolates studied were of worldwide distribution and collected from various hosts (Table1). Thirty-six of these isolates circulated in carnivores, 17 circulated in chiropters, and 2 corresponded to Mokola virus, which has an unknown vector. Twelve viruses were previously described (3). Eight sequences were retrieved from GenBank, and the G genes of 35 isolates obtained from collaborative laboratories were studied for the first time. Isolates were taken from either the brain of an infected mouse or the brain of a suckling mouse after limited passages.
The 55 isolates studied
Sequencing.RNA extraction from tissue, cDNA synthesis, reverse transcription-PCR amplification, and sequencing were done as described by Sacramento et al. (38). Sequences were obtained from both strands and were carefully verified. Only positive-strand primers were used for cDNA synthesis of the negative-strand viral genome, and consensus sequences were determined without subsequent cloning. This procedure avoided the influence of quasispecies variability on the observed genetic polymorphism. G gene sequencing used 23 primers dispersed between positions 2901 and 5543 of the Pasteur virus strain (PV) genome (49).
Sequence analyses.Sequences were aligned using CLUSTAL W (47), and phylogenetic analyses were done using the PHYLIP phylogeny inference package software (version 3.52c; J. Felsenstein, University of Washington, Seattle [http://evolution.genetics.washington.edu/phylip.html ]) or CLUSTAL W. Sawyer's test was used to look for evidence of recombination (40). For each pair of sequences, the test detects silent polymorphic (or informative) sites and calculates values of four parameters: the sum of the squares of the lengths of condensed fragments [SSCF], the sum of the squares of the lengths of uncondensed fragments [SSUF], the observed maximum length of a condensed fragment [MCF], and the observed maximum length of an uncondensed fragment [MUF]. The length of the condensed fragment or uncondensed fragment is the number (whole-sequence alignment) of silent polymorphic sites or of all sites, respectively, between two adjacent polymorphic sites (pair of sequences). To obtain simulated values for the four parameters, the test does 10,000 random permutations of silent polymorphic sites by respecting their class of degeneracy. The probability of a simulated value being larger than the observed value is estimated, and when the probability is less than 0.05 (or 0.01), recombination is evidenced.
The DnaSP program was used for the Hudson-Kreitman-Aguadé (HKA) (22), Tajima (45), and Fu and Li (19) neutrality tests (37). The HKA test compares the observed polymorphism with that expected by the neutral theory of molecular evolution (24), in which similar levels of polymorphism must be seen between and within species. Tajima's test compares two estimates of nucleotide diversity: the mutation parameter (θ) and the difference between pairs of isolates (π). Fu and Li's test is related to Tajima's and compares the number of singletons to that expected by the neutral theory and can be performed with an outgroup.
Relative-rate tests were done using the K2WuLi (54) and RRTree (36) programs. The first determines the nucleotide distance (all substitutions or only transversions) of an outgroup from two sequences of the tested lineages. The RRTree program performs similarly but accepts amino acid sequences and allows several lineages to be tested. It offers a larger choice of methods to estimate distance and can be weighted with a tree.
The rate of evolution was estimated for each pair of isolates that were within the same (or close) lineage(s) and had been isolated within at least 5 years of each other. The outgroup closest to the lineage was used to estimate the distance by using thedS (32). The evolutionary rate was estimated by dividing the difference between the distances of the two members of the pair by the difference in their years of isolation.
WINA was used to estimate along the G gene the proportions of synonymous substitions (dS ) and nonsynonymous substitutions (dN ) according to the method of Nei and Gojobori (32). The program detects regional differences in accumulations of ds and dn(17). Four levels of significance can be obtained (decreasing order): dN > 2dS and dN > 1, dN > 2dS anddN ≤ 1, dN >dS and dN > 1, anddN > dS anddN ≤ 1.
Nucleotide sequence accession numbers.GenBank accession numbers of the G genes of isolates sequenced for this study are noted in Table 1.
RESULTS
Glycoprotein diversity.In order to analyze the driving forces in Lyssavirus evolution, the G genes of 55 isolates of the main variants and GTs (Table 1) were compared. Twenty isolates were sequenced previously (3) or retrieved from GenBank. Thirty-five isolates were sequenced from the G mRNA start signal to the downstream L mRNA start signal (approximately 2,100 nucleotides). Thus, 44 isolates corresponded to the RABV and 11 were RABV-related viruses (GT2 to -7). The ECTO is the most conserved part of the G protein, showing at least 61% amino acid identity. In particular, all cysteines are conserved among all lyssaviruses, as is the glycosylation site N319. By contrast, the other parts have identities as low as 20.5% (SP plus TM) and 24% (ENDO). The G gene 3′ noncoding region (Ψ) shows significant similarities within, but not between, GTs. All these data outline the importance of the structural and functional constraints on the ECTO in contrast to other parts of the G gene and the considerable genetic diversity of lyssaviruses. A neighbor-joining tree showed that the seven GTs are clearly distinct, and several lineages can be distinguished within GT1 (Fig. 1). Among the 44 RABV isolates, 22 illustrate the cosmopolitan variant that was probably caused by human activity through dog importation (25, 42). This lineage probably originated from the Palearctic region (Europe, Middle East, and North Africa), where it remains almost exclusive to some vectors (dog, red fox, raccoon dog, and wolf). It has been spread worldwide and is maintained mainly in dogs (Africa, Asia, and Latin America) but has also been adapted to wildlife vectors (skunk, gray fox, and coyote in North America). Out of the Palearctic region, the cosmopolitan variant coexists with autochthonous variants present before its importation. At least one autochthonous variant has spread in the Arctic region (arctic and red foxes), two have spread in sub-Saharan Africa (dog and mongoose), three have spread in Asia (dog), and several have spread in the Americas by means of both carnivores (skunk and raccoon) and chiropters (insectivorous, frugivorous, and hematophagous bats). In the last group, the variants are specific for the bat species (31, 41). It is obvious from the phylogenetic tree that carnivoran and chiropteran lyssaviruses are branching apart.
Lyssavirus-rooted phylogenetic tree. The tree was estimated using the neighbor-joining method (39) on the basis of the ECTO nucleotide sequence. Bootstrap values of 1,000 replicates indicate the robustness of the corresponding node. The sequences retrieved from GenBank and those described in the work of Badrane et al. (3) are marked with * and †, respectively. RABV, Rabies virus; EBLV-1, European bat lyssavirus 1; EBLV-2, European bat lyssavirus 2; ABLV,Australian bat lyssavirus; DUVV, Duvenhage virus; LBV, Lagos bat virus; MOKV, Mokola virus.
Evolutionary forces and model.We tested to see if recombination has played a significant role in Lyssavirusevolution. Seventeen isolates have available sequences of their G (this paper) and the nucleoprotein (N) (7) genes. We generated a bootstrapped neighbor-joining tree from each of the two genomic regions separated by almost 2.5 kb. Both trees exhibited very similar topologies with significant bootstrap values (data not shown). In addition, very similar trees were obtained with the N (8) and the G (3) genes in lyssaviruses. Using a different strategy, recombination was further assessed by Sawyer's test. We validated the test with a sample of 10 sequences in which few recombination events were artificially introduced. As shown in Table2, recombination was clearly evidenced in this sample for both polymorphic and informative sites since the probabilities of obtaining larger values for the simulated parameters were equal to zero. With the genuine sequences (no artificial recombination), the probabilities were greater than 11%. The whole of 44 RABV G sequences were tested, and once again the probabilities were high (>10%). All these data provided no evidence of recombination.
Probabilities of Sawyer's test for recombination
Tajima's and Fu and Li's neutrality tests were applied to all sequences from carnivoran and chiropteran lyssaviruses using polymorphic sites or only segregating sites. Neither test rejected the neutral model of molecular evolution. Moreover, the Fu and Li and HKA tests were performed separately on each of the chiropteran and carnivoran set of sequences by taking one sequence from the other group as an outgroup. Similarly, neither test rejected the neutral model (data not shown). These results imply a neutral model inLyssavirus evolution.
Evolutionary rate.K2WuLi and RRTree programs were used for relative-rate tests. In K2WuLi, a G gene nucleotide alignment of representative lyssavirus lineages (9 from chiropters and 10 from carnivores) was used, with Mokola virus as an outgroup and excluding Lagos bat virus (too close to the outgroup). Ninety pairwise relative-rate tests were performed using all substitutions or the transversions alone (Table3). Only four comparisons (2.2%) were significant and concerned two chiropteran GT1 lineages (one from Argentina and one from the United States), which have evolved slightly more slowly than two carnivoran GT1 lineages (French fox and sub-Saharan African dog). Using amino acid sequence, RRTree performed a relative-rate test either weighted or not weighted with a phylogenetic tree and used Chandipura virus (Vesiculovirusgenus) as an outgroup. The test also resulted in nonsignificant differences (data not shown), thus supporting similar rates in all lineages of lyssaviruses.
Frequencies of Z-score intervals from 90 K2WULI relative-rate tests
The evolutionary rate was estimated for six selected pairs of close isolates spanning the whole Lyssavirus genus. ThedS rate was estimated to be 3.1 × 10−4 to 5.5 × 10−4/site/year (Table4). Accordingly, a molecular-clock tree (PHYLIP phylogeny inference package) using corrected distances (23) estimated that the most recent ancestor of the cosmopolitan variant existed 284 to 504 years ago. AdN rate of 1.4 × 10−5 to 2.3 × 10−5/site/year was deduced from the rate ofdS . It allowed us to time distant divergences in which dS were saturated.
Evolutionary rate of dS in theLyssavirus G gene
dN and dS analysis.Although thedS /dN rates are globally high (>7), a delicate analysis of dS anddN was done along the G coding region (1,573 nucleotides). The WINA program estimated dS anddN in pairs of sequences from carnivoran or chiropteran lyssaviruses, and thedN /dS ratio was plotted (Fig. 2). Several segments were identified as having a greater dN thandS . The SP, TM, and ENDO parts were among these segments, as expected from their relatively weak constraints. However, in the ECTO, where the constraints are stronger, two segments with high significance (dN > 2dS and dN ≤ 1) were distinctively detected in carnivoran lyssaviruses (G residues 159 to 193 and 326 to 370; nucleotide positions 475 to 577 and 976 to 1,108, respectively) and chiropteran lyssaviruses (G residues 240 to 255 and 275 to 282; nucleotide positions 718 to 763 and 823 to 844, respectively). Interestingly, in carnivoran lyssaviruses they overlap the two major antigenic sites, II and III (6). In contrast, in chiropteran lyssaviruses the two segments are located around the Western blot-positive epitope in regions of unknown functions. Other segments with lower significance (dN > dS and dN ≤ 1) were found. In chiropteran lyssaviruses three segments (residues 139 to 160, 180 to 200, and 332 to 357; nucleotide positions 415 to 478, 538 to 598, and 994 to 1069, respectively) also overlap the main antigenic sites, II and III, while in carnivoran lyssaviruses one segment (residues 420 to 427; nucleotide positions 1,258 to 1,279) is in the C terminus. These peptides under positive selection are thought to bear host adaptation sites.
Ratios of dN todS along the G gene coding region. We plotted thed S/d Nratio for each pairwise comparison (17) of chiropteran (top graph) and carnivoran (bottom graph) lyssaviruses along their G gene coding regions. Threshold lines of the significance of the ratios are shown at values 1 and 2. A schematic representation of the G gene shows the different domains, SP, TM, ENDO, and ECTO, where antigenic sites are indicated with vertical black boxes. Horizontal black or open boxes represent regions of the chiropteran or carnivoran lyssavirus G gene, respectively, which accumulate significantly mored N thand S. *,d N >d S andd N ≤ 1; **,d N > 2d S andd N ≤ 1.
DISCUSSION
Among the 3,400 taxonomic viral species reported in the virus taxonomy book (30a), 54% are RNA viruses and 43% are DNA viruses. It is noteworthy that, among human viruses, 141 species are RNA viruses while only 32 are DNA viruses (H. Badrane, personal analysis). Of the 42 emerging or reemerging infectious diseases between 1996 and 2000 (21a), 23 (55%) were caused by RNA viruses, showing to what extent they threaten the public health. Understanding the evolution of viruses can be a major step in fighting and controlling viral diseases. We studied lyssavirus evolution by analyzing sequences from the external G protein, a main factor in virus-host interactions. We first assessed the importance of evolutionary forces (point mutations and recombination) and model. We did a phylogenetic analysis to examine the possibility of successful host switching of lyssaviruses from chiropters to carnivores. We also estimated the rate of evolution and analyzed thedN /dS ratio along the G gene to search for putative host adaptation segments.
G sequences were obtained from 35 isolates, and 20 additional sequences were analyzed previously (3) or retrieved from GenBank. So far, our panel of isolates represents the most complete picture of lyssaviruses diversity. Sequence analysis indicated that the Ψ noncoding region is hypervariable, although its maintenance at similar sizes (450 to 514 nucleotides) in all lyssaviruses studied is surprising. This may be explained by the role of noncoding regions in the regulation of transcription of lyssaviruses (18). G proteins exhibited great diversity among lyssaviruses, having as little as 54% amino acid identity. Our data clearly show that point mutation and purifying selection are the major forces in lyssavirus evolution. Phylogenetic analysis of Lyssavirus isolates produced very similar trees on genomic regions separated by 2.5 kb. Moreover, Sawyer's test randomly produced parameters (SSCF, MCF, SSUF, and MUF) with larger values than the observed values (P > 0.1). Hence, both analyses suggest no evidence of recombination inLyssavirus evolution. Although recombination has been described for some DNA and RNA viruses (1), Pringle suggested no recombination for the whole orderMononegavirales (35). Different constraints may prevent recombination (53). The viability of recombined genomes (4) and the replication in the same cell do not seem restrictive. However, the most-limiting processes leading to recombination may be a rare coinfection with very distinct lyssaviruses and/or the constant association of the N protein with the viral genome (35).
As lyssaviruses do not recombine significantly and evolve at similar rates, the tree in Fig. 1 should provide a good estimate of their true phylogeny. It strongly indicates that chiropteran lyssaviruses existed long before carnivoran RABV (Fig. 3), and that successful host switching from chiropters to carnivores has occurred in the history of Lyssavirus. Spillovers of lyssaviruses from bats to carnivores are not unusual, and contemporaneous episodes have been described (10, 29). The phylogeny of lyssaviruses implies at least two ancient spillover events, both within the GT1. One spillover that is predicted to have occurred in North America produced the raccoon RABV lineage (eastern United States) and possibly the closely related skunk lineage (central United States), as was suggested by partial G gene sequencing (data not included). However, a second spillover in an unknown region spread the RABV worldwide in a variety of carnivores from an unknown vector. It produced an RABV responsible for 32,000 human deaths per year (52a), constituting a historical spillover, and made rabies a public health problem. Spillovers of lyssaviruses from chiropters to other animals, mainly carnivores, may have occurred repeatedly in history and still occur (10). However, it is not known why only rare spillovers have succeeded in being maintained through time or why others have been extinguished. The present work give strong evidence that what is today known as carnivoran rabies is indeed a result of a few successful episodes of host switching from bats.
Lyssavirus phylogenetic tree with a molecular clock (PHYLIP phylogeny inference package) derived from nonsynonymous corrected distances (23). Six of the sevenLyssavirus GTs are represented (except GT3). Bold branches distinguish chiropteran Lyssavirus lineages. The geographic locations and vectors of RABV main lineages are indicated. Curved arrows symbolize the two spillover events. Timing estimations of spillovers and of the most recentLyssavirus ancestor are indicated on the scale.
The rate of Lyssavirus evolution was roughly estimated to be 3.1 × 10−4 to 5.5 × 10−4 dS /site/year. This rate date the divergence of the cosmopolitan variant to 284 to 504 years ago, which is in agreement with the intense human movements of the last five centuries. dS are generally neutral and more sensitive to evolutionary time thandN . However, they saturate faster over long periods of evolution. Therefore, to evaluate older divergences, the rate of dN was used. It was deduced from that ofdS to an average of 1.85 × 10−5/site/year. This rate seems low compared to that for Human immunodeficiency virus or Human influenza A virus but similar to that estimated for the G gene ofEbola virus (3.6 × 10−5substitutions/site/year), another member of theMononegavirales (44). The lower rate of molecular evolution in some RNA viruses might be due to a higher fidelity of their polymerase or a lower degree of viral multiplication. A phylogenetic tree was constructed with nonsynonymous corrected distances (23) by assuming a molecular clock (PHYLIP phylogeny inference package) (Fig. 3). It dated the host switching of the RABV from chiropters to carnivores to 888 to 1,459 years ago. This estimation seems questionable, as carnivoran rabies was described back in Mesopotamian civilizations about 4,000 years ago (46). Nevertheless, the type of Lyssavirus responsible for Mesopotamian rabies is unknown, and it is admitted that RABV lineages may be extinct because of historical or environmental events or of the fatality of the disease. Timing the emergence of carnivoran rabies to 4,000 years ago makes the evolution of lyssaviruses extremely slow (4.1 × 10−6 to 6.7 × 10−6 dN /site/year). Hence, the extinction of the Mesopotamian Lyssavirus is the most likely hypothesis, as the neutral model of evolution implies a random genetic drift, which may lead to extinction and loss of polymorphism (24).
Because lyssaviruses migrate from the peripheral to the central nervous systems, they hide from and avoid the immunity of the host. Therefore, it is not surprising that their G genes are globally not subject to selection pressure, as suggested by the neutrality tests. Nevertheless, positive selection may be detected regionally (17) or even at single amino acid sites (43). Hence, several G segments accumulating more dN thandS were distinctively detected in carnivoran and chiropteran lyssaviruses. The conformational main antigenic sites bearing several residues linked to virulence (9, 34) are the putative sites for host adaptation (not including immune system escape), particularly in carnivores and to a lesser extent in chiropters. In addition, two distinct segments of unknown functions are putative sites for adaptation specifically to chiropters.
Insectivorous bats are possible vectors for six LyssavirusGTs and are exclusive to three. It can be speculated that lyssaviruses originated from an insect rhabdovirus, which insectivorous bats contracted from insects. Several arguments are in favor of this speculation. Most of the other Rhabdoviridae genera, exceptNovirhabdovirus, have isolates from insects. Strikingly, three unclassified viruses proposed to belong to theLyssavirus genus have been isolated so far only from insects, i.e., the kotonkan, Obodhiang (5), and Rochambeau (13) viruses. Moreover, Mokola virus, which was also isolated from an insectivorous animal (shrew), was demonstrated to replicate in inoculated Aedes aegypti mosquitoes (2). If this contraction of an insect rhabdovirus by bats really resulted in the emergence of lyssaviruses, it can be timed to the most recent common ancestor of lyssaviruses, about 7,080 to 11,631 years ago (Fig. 3).
ACKNOWLEDGMENTS
We deeply thank the following worldwide collaborators, who provided viral samples: D. Briggs (Kansas State University), D. Lodmell (Rocky Mountain Laboratories, Hamilton, Mont.), A. Fayaz (IPI, Teheran, Iran), B. Swanepoel (NIV, Johannesburg, South Africa), and H. Bourhy (IP, Paris, France). We are grateful to two anonymous reviewers, whose constructive criticism and suggestions improved the manuscript.
H.B. was the recipient of fellowships from the Moroccan Government and French-Moroccan cooperation.
FOOTNOTES
- Received 26 March 2001.
- Accepted 5 June 2001.
- Copyright © 2001 American Society for Microbiology
