Skip to main content
Advertisement
  • Loading metrics

Large-Scale Phylogenomic Analysis Reveals the Complex Evolutionary History of Rabies Virus in Multiple Carnivore Hosts

  • Cécile Troupin,

    Affiliation Institut Pasteur, Unit Lyssavirus Dynamics and Host Adaptation, WHO Collaborating Centre for Reference and Research on Rabies, Paris, France

  • Laurent Dacheux,

    Affiliation Institut Pasteur, Unit Lyssavirus Dynamics and Host Adaptation, WHO Collaborating Centre for Reference and Research on Rabies, Paris, France

  • Marion Tanguy,

    Affiliations Institut Pasteur, Unit Lyssavirus Dynamics and Host Adaptation, WHO Collaborating Centre for Reference and Research on Rabies, Paris, France, Institut Pasteur, Genomics Platform, Paris, France

  • Claude Sabeta,

    Affiliation Agricultural Research Council, Onderstepoort Veterinary Institute, OIE Rabies Reference Laboratory, Pretoria, South Africa

  • Hervé Blanc,

    Affiliation Institut Pasteur, Centre National de la Recherche Scientifique UMR 3569, Viral Populations and Pathogenesis Unit, Paris, France

  • Christiane Bouchier,

    Affiliation Institut Pasteur, Genomics Platform, Paris, France

  • Marco Vignuzzi,

    Affiliation Institut Pasteur, Centre National de la Recherche Scientifique UMR 3569, Viral Populations and Pathogenesis Unit, Paris, France

  • Sebastián Duchene,

    Affiliations Marie Bashir Institute for Infectious Diseases and Biosecurity, Charles Perkins Centre, School of Life and Environmental Sciences and Sydney Medical School, The University of Sydney, Sydney, Australia, Centre for Systems Genomics, University of Melbourne, Parkville, Victoria, Australia

  • Edward C. Holmes,

    Affiliation Marie Bashir Institute for Infectious Diseases and Biosecurity, Charles Perkins Centre, School of Life and Environmental Sciences and Sydney Medical School, The University of Sydney, Sydney, Australia

  • Hervé Bourhy

    hbourhy@pasteur.fr

    Affiliation Institut Pasteur, Unit Lyssavirus Dynamics and Host Adaptation, WHO Collaborating Centre for Reference and Research on Rabies, Paris, France

Abstract

The natural evolution of rabies virus (RABV) provides a potent example of multiple host shifts and an important opportunity to determine the mechanisms that underpin viral emergence. Using 321 genome sequences spanning an unprecedented diversity of RABV, we compared evolutionary rates and selection pressures in viruses sampled from multiple primary host shifts that occurred on various continents. Two major phylogenetic groups, bat-related RABV and dog-related RABV, experiencing markedly different evolutionary dynamics were identified. While no correlation between time and genetic divergence was found in bat-related RABV, the evolution of dog-related RABV followed a generally clock-like structure, although with a relatively low evolutionary rate. Subsequent molecular clock dating indicated that dog-related RABV likely underwent a rapid global spread following the intensification of intercontinental trade starting in the 15th century. Strikingly, although dog RABV has jumped to various wildlife species from the order Carnivora, we found no clear evidence that these host-jumping events involved adaptive evolution, with RABV instead characterized by strong purifying selection, suggesting that ecological processes also play an important role in shaping patterns of emergence. However, specific amino acid changes were associated with the parallel emergence of RABV in ferret-badgers in Asia, and some host shifts were associated with increases in evolutionary rate, particularly in the ferret-badger and mongoose, implying that changes in host species can have important impacts on evolutionary dynamics.

Author Summary

Zoonoses account for most recently emerged infectious diseases of humans, although little is known about the evolutionary mechanisms involved in cross-species virus transmission. Understanding the evolutionary patterns and processes that underpin such cross-species transmission is of importance for predicting the spread of zoonotic infections, and hence to their ultimate control. We present a large-scale and detailed reconstruction of the evolutionary history of rabies virus (RABV) in domestic and wildlife animal species. RABV is of particular interest as it is capable of infecting many mammals but, paradoxically, is only maintained in distinct epidemiological cycles associated with animal species from the orders Carnivora and Chiroptera. We show that bat-related RABV and dog-related RABV have experienced very different evolutionary dynamics, and that host jumps are sometimes characterized by significant increases in evolutionary rate. Among Carnivora, the association between RABV and particular host species most likely arose from a combination of the historical human-mediated spread of the virus and jumps into new primary host species. In addition, we show that changes in host species are associated with multiple evolutionary pathways including the occurrence of host-specific parallel evolution. Overall, our data indicate that the establishment of dog-related RABV in new carnivore hosts may only require subtle adaptive evolution.

Introduction

Revealing how viruses jump species boundaries and establish productive infections in new hosts is key to understanding disease emergence. As most recent emerging and re-emerging viruses have RNA genomes [1], it is of central importance to understand the drivers of RNA virus evolution, diversification and cross-species transmission. Clearly, successful virus emergence has diverse causes, likely involving anthropogenic, social and environmental factors [2]. However, the capacity of the viral genome to vary and generate advantageous mutations is also an important element, enabling RNA viruses to exploit new niches, including novel host species, often more rapidly than DNA-based organisms [1, 3, 4]. One important manifestation of RNA virus evolution and diversification is the rate of evolutionary change (i.e. nucleotide substitution), with analyses of how this parameter varies by host species providing important information on the nature of virus-host interactions.

Disease emergence results from complex mechanisms that shape the ability of a virus to be maintained within its primary host species, then be serially transmitted to a new host species and initiate a pathologic process to cause disease [5]. As such, lyssaviruses (family Rhabdoviridae), the causative agents of rabies–an acute and almost invariably fatal encephalomyelitis in humans–represent an informative case study to examine the relationship between virus genetic diversity and disease emergence. In particular, the natural history of these zoonotic viruses provides an excellent model to study how replication in different host species alters the selection pressures that act on virus genomes. Lyssaviruses are single-stranded, negative-sense RNA viruses with a genome size of approximately 12 kb that encodes five proteins: the nucleoprotein (N), the phosphoprotein (P), the matrix protein (M), the glycoprotein (G) and the Large protein or polymerase (L). Currently, the lyssaviruses are classified into 14 species and one tentative species [6]. Like other RNA viruses, lyssaviruses exhibit high rates of mutation due to a lack of proofreading activity in the L protein [7]. Notably, although many mammalian species appear to be susceptible to lyssavirus infection, the virus is only able to establish sustained transmission networks in a relatively small number, indicating that there are major barriers to successful cross-species transmission [811].

One species of lyssavirus, rabies virus (RABV), is present worldwide and circulates in a diverse set of reservoir hosts among the mammalian orders Chiroptera and Carnivora [12]. Its natural evolution provides an illustrative example of multiple host switches, in turn enabling comparative studies of the evolutionary patterns, processes and dynamics associated with host adaptation. Previous studies demonstrated that RABV isolates fall into two major phylogenetic groups; the bat- and the dog-related RABV groups [8, 13, 14]. The ‘bat-related’ RABV group is confined to New World viruses circulating mainly among bats, as well as in some terrestrial carnivores such as skunks and raccoons [1417]. In contrast, the ‘dog-related’ RABV group contains viruses circulating worldwide in dogs, as well as in wildlife carnivores in specific geographic areas such as foxes and raccoon dogs in Europe, foxes in the Middle East, raccoon dogs and ferret-badgers in Asia, skunks, foxes, coyotes and mongooses in the Americas, and mongooses in Africa [14, 16, 1822]. Importantly, dogs are responsible for more than 99% of the human rabies cases worldwide [23] and are likely the main vector for the inter-species transmission of dog-related RABV.

Previous phylogenetic analyses have largely been performed on individual genes [1319, 21, 2429] with a few assessing the full-length viral genome [20, 30, 31]. In addition, most of these phylogenetic studies were performed on relatively small numbers of sequences originating from one specific geographical area and/or associated with a specific animal host [20, 22, 30, 32, 33]. Despite these limitations, these studies are consistent in showing that RABV is subject to strong purifying selection [10] coupled to geographical clustering that is occasionally disrupted by human mediated dispersion [13, 34, 35]. Recently, it was shown that nucleotide substitution rates in RABV vary markedly among those viruses infecting bats, such that rates in tropical and subtropical species were markedly higher than those from temporal bat species, perhaps reflecting a combination of host and environmental factors [36]. However, equivalent data for dog-related RABV are lacking. In addition, whether evolutionary rates in RABV vary among wild carnivores and domestic dogs is unknown, although studies in other systems have revealed that rates of RNA virus evolution may differ between wild and domestic animals [37]. Clearly, the large-scale analysis of RABV, particularly comprising full-length genome sequences, is needed to reveal the nature of the selection pressures associated with host switching. That the RABV genome encodes a limited number of proteins that necessarily have multifunctional roles [38], and hence potentially large-scale epistasis, also means that these selection pressures may be complex.

Herein we present the first phylogenetic study of RABV on a genome-wide and global scale, utilizing a data set of 321 whole-genome sequences sampled from 66 countries over a time period of 65 years, with the aim of inferring those evolutionary patterns and processes associated with host-switching. In particular, we compared RABV from wild carnivores and in domestic dogs with respect to selection pressures, evolutionary rates, and the time-scale of their evolutionary history. Importantly, the size of the data set allowed us to reveal any heterogeneity in evolutionary rates among RABV adapted to different primary hosts, and determine the complex evolutionary dynamics of RABV as it adapts to new hosts.

Results

Host and geographical clustering of RABV

A phylogenetic analysis was performed on the (99%) full-length genome sequences of 321 RABV sequences sampled from 66 countries (S1 Fig, S1 Table). Of these viruses, 170 were newly sequenced as part of this study. As expected given the low levels of recombination in RABV, the topology of the maximum likelihood (ML) tree performed on the five concatenated RABV genes (Fig 1) was similar to that obtained for each individual gene (N, P, M, G and L genes) and for the concatenated non-coding sequence (S2 Fig). In particular, two major phylogenetic groups were apparent, corresponding to bat- and dog-related RABVs, each of which can be further subdivided into several major clades. This is consistent with previous analyses of smaller data sets and on individual RABV genes [13, 14, 16, 29].

thumbnail
Fig 1. Maximum likelihood phylogeny of 321 RABV sequences from five concatenated genes.

The major clades of RABV are indicated in boxes. The names of subclades and lineages defined for the Arctic-related, Asian and Cosmopolitan clades are detailed in S1 Table, with corresponding bootstrap values shown for major nodes. The tree is mid-point rooted for clarity only, and shows the division into bat-related RABV including the RAC-SK and bat clades, and dog-related RABV including the Africa-2, Africa-3, Arctic-related, Asian, Indian subcontinent and Cosmopolitan clades.

https://doi.org/10.1371/journal.ppat.1006041.g001

The bat-related group contained two major clades, one including the bat RABVs circulating in the Americas, and the other (RAC-SK) comprising viruses from American skunks and raccoons (Fig 1). In turn, the RAC-SK group contained a number of ‘subclades’ corresponding to Mexican skunks (MeSK-1), North American raccoons (RAC) and South-Central skunks (SCSK) as previously described (S3 Fig) [14, 16, 17, 39].

Similarly, the dog-related group includes six major clades supported by high bootstrap values (S2 Table and Fig 1), and previously identified as the Africa-2, Africa-3, Arctic-related, Asian, Cosmopolitan and Indian subcontinent clades [13]. The phylogenetic analysis based on the five concatenated genes was particularly informative, allowing us to distinguish various subclades and lineages among these six major clades, with some of which are characterized for the first time here (Fig 1 and supplementary text).

Some of these clades and subclades are of particular interest. The SEA2 subclade contains viruses from China and is divided into two lineages, SEA2a and SEA2b, corresponding to isolates from dogs and ferret-badgers, respectively. Subclade SEA5 appears to be specific to RABV circulating in ferret-badgers in Taiwan, an epidemiological cycle that was only identified recently [21, 40, 41]. For the first time, we were also able to fully characterize full-length genome sequences of RABV isolates belonging to the Africa-3 clade (n = 6). These viruses circulating in Southern Africa are monophyletic and phylogenetically distinct from the other major RABV clades, particularly those circulating in Africa [19, 42, 43].

Temporal dynamics and spread of the dog-related RABV

To determine the evolutionary dynamics of RABV, we first determined whether individual data sets contained sufficient temporal structure to undertake detailed molecular clock analyses by performing a regression of root-to-tip genetic distance against the year of sampling. Notably, no correlation between time and genetic divergence was found when the sequences of both bat- and dog-related RABV groups were analyzed together, indicating that there is extensive variation in the rate of RABV evolution among these taxa (and hence that they should not be combined in molecular clock studies) (S4A Fig). In addition, no temporal structure was observed when the sequences of bat-related RABV were analyzed separately, indicating that this subset of viruses is not evolving uniformly (S4B Fig) as noted previously [36]. However, a clear association between genetic distance and time (i.e. a molecular clock) was observed for the dog-related group alone (S4C Fig), allowing us to estimate substitution rates, and hence times to common ancestry, more precisely in this cluster using a Bayesian approach.

The mean rate of evolutionary change in the dog-related RABV was estimated to be 2.44 x 10−4 subs/site/year (95% HPDs of 2.10–2.80 x 10−4 subs/site/year) for the five concatenated genes. Importantly, we were also able to compare the substitution rate of each RABV gene and of the concatenated non-coding regions from the same genomic sequence data set. These estimates varied in the following ascending order: N, L, G, M and P (Fig 2A). However, only the P gene had a nucleotide substitution rate considerably higher than those of N and L genes. As expected, the evolutionary rate in the non-coding regions was significantly higher than those of the coding regions, indicative of weaker selective constraints.

thumbnail
Fig 2. Evolutionary rates of RABV genes in the dog-related group.

(A) Rates of nucleotide substitution per site, per year were estimated for each RABV gene: nucleoprotein (N), phosphoprotein (P), matrix (M), glycoprotein (G) and polymerase (L), for the concatenated non-coding regions (NC) and for the five concatenated RABV genes (5 genes). Both the mean and the 95% highest posterior density (HPD) values on the rate are shown. (B) Substitution rates in the N and G genes of the dog-related group RABV, a sub-set of RABV circulating in mongooses (MG) in Africa-3 clade and in the Caribbean, in ferret-badgers (FB) in Asia, and in dogs in Asia and Africa. Note the different y-axes (rates) in both cases.

https://doi.org/10.1371/journal.ppat.1006041.g002

The estimation of a reliable substitution rate allowed us to determine the mean times to common ancestry (TMRCA) for each RABV clade, subclade and lineage defined above (Fig 3; estimates for different sub-clades shown in S3 Table and discussed in the Supplementary text). For this analysis we utilized the concatenated coding genes as these had the lowest variance. Notably, these TMRCA estimates exhibited less uncertainty than previous studies performed on N and/or G genes alone [13, 25, 27, 29, 43, 44]. Briefly, the TMRCA of the dog-related RABV group was estimated to be approximately between 1308–1510 (95% HPD; mean of 1404). Within this group, the Indian subcontinent clade branched basally and appeared to diversify between 1733–1840 (mean of 1785). The TMRCA of the Asian clade was estimated to be between 1535–1677 (mean of 1604), which is in accordance with other studies [44]. The emergence of the Africa 2 clade was estimated to be between 1750–1852 (mean of 1802), similar to the mean TMRCA found in a previous study conducted on complete N and G genes [27]. The Arctic-related clade appeared between 1725–1815 (mean of 1770), slightly earlier than previously estimated [25], while the Africa-3 clade emerged between 1710–1815 (mean of 1756) in accordance with another study [43]. Finally, for the first time, we estimate that the TMRCA for the Cosmopolitan clade existed between 1687–1773 (mean of 1730).

thumbnail
Fig 3. Maximum clade credibility phylogeny of 248 dog-related RABV utilizing five concatenated genes.

Tip times represent the time (year) of sampling. Bayesian estimates of divergence time are also shown. Upper and lower limits of the 95% highest posterior density (HPD) estimates and the posterior probability values are shown for major nodes.

https://doi.org/10.1371/journal.ppat.1006041.g003

Variation in evolutionary rates among hosts and selection pressures along the RABV genome

The root-to-tip regression analysis also revealed that different groups of dog-related RABV have seemingly evolved at different rates (S4C Fig), with a number appearing as distinct outliers. Interestingly, these outliers were confined to RABV circulating in mongooses in Africa (Africa-3 clade) and in ferret-badgers in Asia (the SEA5 subclade and SEA2b lineage), suggesting that they might represent species-specific variation. To address this, we compared the evolutionary rates of these clusters to both the entire dog-related RABV group and to subsets of this group representing dog-related viruses circulating in Africa and Asia, and to mongoose viruses circulating in the Caribbean. For this analysis we focused on the N and G genes as they comprise the largest data sets. This analysis revealed that the N gene of those viruses circulating in ferret-badgers in Asia (n = 81) and in mongooses in Africa (n = 47) evolved between 2–4 times more rapidly than those of the whole dog-related group (n = 248), at rates of 7.82 x 10−4 subs/site/year (95% HPD 3.14–12.83 x 10−4 subs/site/year) and 5.88 x 10−4 subs/site/year (95% HPD 3.67–8.11 x 10−4 subs/site/year), respectively (Fig 2B). Importantly, these estimates and their associated uncertainty do not overlap with those for the dog-related group as a whole. This finding is confirmed using smaller subsets of dog-related RABV from more closely related geographically settings in Asia and in Africa (Fig 2B). Interestingly, the rate of RABV evolution in mongooses in Africa is two times higher than that of RABVs from mongooses in the Caribbean (i.e. Puerto Rico, Cuba and Grenada) that belong to the Cosmopolitan clade (Fig 2B). Although less rate variation was observed in the G gene, RABV associated with ferret-badgers in Asia still evolved considerably more rapidly than those obtained with the different subsets of dog-related RABV (Fig 2B). These results were confirmed by using different nucleotide substitution models and a hierarchical phylogenetic model approach (S4 Table) [45, 46].

To determine if the variation in rates of evolutionary change might result from differing selection pressures, we first compared the ratio of nonsynonymous (dN) to synonymous (dS) substitutions per site. This analysis was performed on each of the five RABV genes of the two major RABV groups. For each gene, the dN/dS ratios of the bat- and dog-related groups are very similar (and very low) and followed the same ascending order between genes: N, L, M, G and P genes (Table 1). Furthermore, we explored the number of positively selected sites using several different approaches (SLAC, FUBAR and FEL) [47, 48]. In each of the two major RABV groups, one position was identified as positively selected by at least two of these methods: positions 496 and 484 in the G protein for the bat- and dog-related groups, respectively (Table 1). Interestingly, the dN/dS of the N and G genes for the branches leading to sequences found to be outliers in the analysis of evolutionary rates (Africa-3 clade and ferret-badgers in Asia) were 1.4 to 4.7 times higher than those of dog-related RABV data sets used as controls (S4C Fig and S5 Table), but still relatively low. Together, these results are generally indicative of strong purifying selection among all sites and branches of the RABV phylogeny.

thumbnail
Table 1. Selection pressures in five genes from bat- and dog-related RABVs.

https://doi.org/10.1371/journal.ppat.1006041.t001

To investigate selection pressures in greater detail we utilized a modified MEME analysis that considered internal branches of the tree only (as external branches often contain transient deleterious nonsynonymous substitutions yet to be removed by purifying selection) [49]. Using this MEME-internal analysis, we identified nine positions to be under positive selection (N436, P55, P154, P265, G198, G476, L430, L681, L2091). In addition, position G484 that was identified as positively selected using SLAC, FUBAR and FEL was not significant (at the p<0.05 level) in the MEME-internal analysis (p-value = 0.084).

Finally, it was also clear that specific amino acid substitutions characterized RABV circulating in mongooses in Africa (Africa-3 clade) and in ferret-badgers in Asia (SEA5 subclade and SEA2b lineage) (S6 Table). Four substitutions were specific (i.e. not present in any other dog-related RABV sequences) to mongoose RABV: two in the nucleoprotein, from Asp to Asn at codon position 88 (Asp-N88-Asn) and Leu-N108-Ile, and two in the glycoprotein–Ser-G223-Asn and Pro-G386-Ser. The case of the ferret-badger was more interesting as the host jump to this species from dogs has occurred independently in the SEA5 and SEA2b clades, allowing us to determine whether cross-species transmission in this case is associated with parallel viral evolution. This analysis revealed that two amino acid substitutions were common to all ferret-badger viruses across both clades: Leu-N374-Ser and Lys-L200-Arg. The Leu-N374-Ser substitution is particularly noteworthy as it only occurs in the ferret-badger, this residue is normally highly conserved in RABV, and Leu-to-Ser is a non-conservative amino acid change. Hence, we suspect that Leu-N374-Ser, and perhaps Lys-L200-Arg, facilitate RABV adaptation to ferret-badgers. Notably, neither of these sites was found to be subject to positive selection using the methods employed here (Table 1).

Discussion

The central aim of this study was to determine whether the patterns and processes of RABV evolution vary between viruses sampled from different host species reflect the impact of cross-species transmission. To that end we present the largest phylogenomic analysis of RABVs circulating worldwide performed to date. Although the topology of the RABV phylogeny is similar to those obtained previously [13, 14, 16, 29], it clearly presents a more comprehensive and precise reconstruction of evolutionary history of this virus. In particular, the analysis of the five concatenated genes allowed us to obtain a finer-scale dating of the emergence of the major clades with narrower confidence intervals than obtained previously [13, 25, 27, 29, 43, 44].

RABV undergoes relatively frequent cross-species transmission [8, 11, 13, 18], which provides an opportunity to determine whether host jumping impacts rates of evolutionary change. Notably, we found no correlation between root-to-tip genetic distance and sampling time in the bat-related RABV group, nor when combined with dog-related RABV group, indicating that these viruses have not evolved in a clock-like manner, with substantial rate variation already observed in bat-associated RABV [36]. In contrast, a strong association between genetic divergence and time (i.e. a molecular clock) was observed within the dog-related RABV group, with a mean evolutionary rate of 2.44 x 10−4 subs/site/year (95% HPDs of 2.10–2.80 x 10−4 subs/site/year) for the five concatenated genes. This estimate is evidently more precise than those determined previously [13, 25, 27, 30, 44, 5053].

Despite the relative rate constancy in the dog-related RABV, it was striking that some of the clades or sub-clades have experienced substantially higher rates of nucleotide substitution. In particular, viruses circulating in ferret-badgers in Asia (mainland China and Taiwan) and in mongooses in Africa have evolved at least twice as rapidly as those of the dog-related group. Although there is some uncertainty in these rate estimates, they do not overlap with the estimates for the entire dog-related RABV group. Determining the evolutionary basis to this rate variation is more complex. Changes in evolutionary rate could only be driven either by changes in background mutation rate (which we consider unlikely to differ between dog-related RABV) or, more likely, by changes in the population size and/or incubation time that may vary among different animal hosts [36]. It is also possible that the evolutionary rates estimated here have been impacted by time-dependency, such that they are elevated toward the present (i.e. in closely related sequences sampled recently) due to the presence of transient deleterious mutations that have yet to be removed by purifying selection [54]. However, while this may in part explain the high rate in the recently sampled RABV from ferret-badgers, it is unlikely to explain the higher evolutionary rate in mongoose RABV whose evolutionary history sampled here covers a longer time period. In the case of the ferret-badgers, two amino acid changes (Leu-N374-Ser and Lys-L200-Arg) have evolved in parallel in the two clades associated which is compatible with the occurrence of adaptive evolution, and which have in turn elevated the nucleotide substitution rate. That these two sites were not detected in analyses of dN/dS suggests that these methods may have limitations when identifying adaptive evolution involving limited amounts of amino acid change.

Our analysis also showed that the nucleotide substitution rate varied markedly according to the gene analyzed in the ascending order: N, L, G, M and P. As expected, the two proteins often described as more conserved for RABV—N and L—exhibited the lowest rates, as well as the lowest dN/dS ratios, indicating that they are subject to the strongest purifying selection. Notably, the highest substitution rate and dN/dS was observed in the P protein, perhaps reflecting the weak structural organization of the C-term part of this protein [55, 56].

The presence of relatively constant molecular clock also enabled us to provide a more robust time-scale for the evolution of the principal geographical clusters of dog-related RABV (Fig 3, S3 Table). Accordingly, we estimate that the most recent ancestor of all dog-related RABV dates to between 1308 and 1510. Consequently, any older canid RABV lineages, proposed to have circulated in the Middle-East more than 2000 years ago [57, 58], have not survived to be sampled in the current study. Interestingly, the timing of the most recent ancestor of all dog-related RABV circulating to date coincides with the development of the world's first truly global trade network following the explorations of Columbus, Vasco da Gama and Zheng He, commissioned by the Spanish, the Portuguese and the Chinese Ming Dynasty, respectively. This age of exploration and colonization contributed to the establishment of new long distance commercial practices and transoceanic shipping services between 1450 and 1750 [59]. The concomitant dissemination of RABV during this period, probably by dogs travelling by boats with their owners, therefore provides a powerful example of the early human-mediated dissemination of a zoonotic disease. In addition, all the ancestors of the major clades found circulating today in North and South America, Africa, Asia and Europe originated between 1687 and 1840 at the apogee of this international trade and colonization process [59]. This is further exemplified by the global spread of the Cosmopolitan clade.

A fundamental question in evolutionary virology is how and why some viruses are seemingly better able to jump species boundaries than others. A compelling theory is that the more closely related the host species in questions, the greater the chance of successful transmission [9, 60, 61]. However, it is unclear how strictly this theory holds for RABV [11], and our results confirm species jumps of RABV among animal species of the order Chiroptera, and from bats to striped skunks (Mephitis mephitis) [14, 62, 63]. In addition, there is also clearly a geographic component to cross-species transmission as bat-related RABVs are only found in the Americas. More notably, our study clearly confirms that although spill-over infections from wildlife species to dog take place, species jumps involving dog-related RABVs generally occur from dogs to wildlife species of the order Carnivora; not only to the family Canidae (dog, red fox, raccoon dog), but also to more distant species belonging to the families Mustelidae (ferret-badger), Herpestidae (mongoose) and Mephetidae (skunk) (S5 Fig) [13, 18, 42, 43]. These changes in primary animal host species have occurred independently in different localities and at different times during RABV evolution. Further, some carnivore species, notably skunks, are infected by RABV of both dog and bat origin [14, 16, 39].

Revealing the respective roles of genetic drift and the selection of advantageous mutations in shaping the genetic diversity of RABV, particularly during host shifts, is a central evolutionary question. There is currently no definitive data on whether dog-related RABV emergence requires active adaptive evolution (i.e. positive selection) to in a new host species, or whether it is largely a chance process involving ecological factors facilitating the transmission of a viral strain with the pre-existing necessary genetic characteristics [64]; the latter has been proposed for the repeated outbreaks of bat-related RABV in striped skunks and gray foxes in Arizona [14] and of gray foxes due to skunk-associated RABV in California [65]. Our analysis showed that the dog-related RABV group is subject to strong purifying selection, and when positive selection did occur on internal branches of the phylogenetic tree it was not obviously associated with host jumping. As noted above, however, the failure to detect positive selection in the case of ferret-badger RABV despite the occurrence of parallel evolution suggests that these methods may suffer from false-negatives.

Successful cross-species transmission is a complex ecological and evolutionary process, beginning with exposure and contact between the two species, followed by the successful infection of the new host species, and potentially host-adaptive evolution to enable long-term sustained transmission [66, 67]. However, due to complex interactions among the five viral proteins and with their cellular counterparts, including epistasis [68], it is often difficult to clearly determine which mutations are advantageous or fixed by genetic drift. Moreover, some mutations in the RABV P protein can improve the modulation of the innate immune response of the host but reduce replication efficiency [69]. That two amino acid changes have evolved in parallel in the ferret-badger alone suggests that they have played a role in host adaptation. Further, it is possible that some of the other amino acid substitutions that define individual viral clades associated with different host species represent host-adaptive sites that have not been identified as positively selected through simple analyses of dN/dS. Clearly, additional large-scale analyses of RABV based on full-length genome sequences, extending that presented here, followed by linked experimental studies including generation of mutant RABVs by reverse genetics and phenotypic testing, are needed to reveal the nature of complex evolutionary processes that occur during host switching.

In conclusion, RABV is capable of infecting many mammals but paradoxically is maintained in distinct epidemiological cycles associated with animals almost exclusively from the orders Carnivora and Chiroptera. This strict association between RABV and host-species most likely arose from a combination of historical human-mediated spread of RABV and jumps into new primary host species. These data also suggest that the establishment of dog-related RABV in new carnivore hosts may only require subtle adaptive evolution as demonstrated by parallel evolution in the ferret-badger. Evidently, along with more defined analyses of individual mutations, additional studies are needed to determine the role played by the frequency of exposure, animal host behavior, density of the recipient species, duration of incubation and optimum infectious doses in cross-species transmission.

Materials & Methods

Samples

A total of 321 complete genome sequences of RABV isolates were analysed, originating from a wide variety of host species and collected in 66 countries between 1950 and 2015. Details of these isolates are described in S1 Table and S1 Fig. Among these genome sequences, 170 came from the archives of the World Health Organization Collaborative Center for Reference and Research on Rabies, or from the National Reference Centre for Rabies, both located at Institut Pasteur, Paris, France. These samples were newly sequenced as part of this study. These data were combined with 151 full-length genome sequences extracted from GenBank and selected to be representative of the overall phylogenetic diversity of RABV.

RNA extraction and next-generation sequencing

Total RNA was extracted using Trizol (Ambion) according to the manufacturer’s instructions from primary brain samples or after an amplification passage on suckling mouse brain. RNA was then reverse transcribed using Superscript III reverse transcriptase with random hexamers (Invitrogen) according to manufacturer’s instructions. The complete viral genome (excluding the 3’ and 5’ extremities, corresponding to the leader and the trailer regions, respectively) of 160 isolates was amplified with six overlapping PCR fragments by using the Phusion polymerase (ThermoFisher). Details of primers are given in S7 Table. After electrophoresis, each PCR fragment was independently purified using the NucleoSpin Gel and PCR clean-up kit (Macherey-Nagel) and quantified using Picogreen dsDNA quantification kit (Invitrogen). For each sample, all six PCR fragments were pooled with equimolar proportions to obtain 500 ng of dsDNA.

Different protocols were used for the preparation of libraries and next-generation sequencing on Illumina platforms (NextSeq 500, HiSeq2000, HiSeq2500 or MiSeq platforms), depending on the isolates considered (details provided in S1 Table). Briefly, three different protocols were used: (i) dsDNA was fragmented by ultrasound with Bioruptor (Diagenode), libraries were prepared using NEXTflex PCR-Free DNA-Seq kit (Bioo Scientific), and then sequenced using an 100 or 150 nucleotides single-end strategy on the HiSeq2500 platform or a 2 x 300 nucleotides paired-end strategy on the MiSeq platform, (ii) dsDNA was fragmented by NEBNext dsDNA fragmentase (New England Biolabs), libraries were prepared using NEBNext Ultra DNA Library Prep kit (New England Biolabs) and sequenced using an 100 nucleotides single-end strategy on the NextSeq500 platform, and (iii) dsDNA libraries were constructed using Nextera XT kit (Illumina) and sequenced using a 2 x 150 nucleotides paired-end strategy on the NextSeq500 platform. For nine remaining isolates (S1 Table), the viral RNAs were reverse transcribed using Superscript III reverse transcriptase (Invitrogen) and then amplified using the whole-transcription amplification (WTA) protocol (QuantiTect Whole Transcriptome kit; Qiagen) as previously described [70]. dsDNA was fragmented by ultrasound, libraries were prepared using TruSeq protocol (Illumina) and sequenced using an 100 nucleotides single-end strategy on the HiSeq2000 platform. Finally, the sequence of 09035FRA was determined using a shotgun base approach [31].

Genome sequence analyses

All reads were pre-processed to remove low-quality or artifactual bases. Library adapters, PCR primers used for amplification of the genome, and base pairs occurring at 5’ and 3’ ends with a Phred quality score <25 were trimmed using AlienTrimmer as implemented in Galaxy [7174] (https://research.pasteur.fr/en/tool/pasteur-galaxy-platform/). Reads with lengths of less than half of the original read after these pre-processing steps or those containing >20% of bp with a Phred score of <25 were discarded. The filtered reads were then mapped to complete genome sequences specific for each RABV clade obtained from GenBank using the CLC Genomics Assembly Cell (http://www.clcbio.com/products/clc-assembly-cell/) implemented in Galaxy. The majority nucleotide (>50%) at each position with a minimum of coverage of 200 was used to generate the consensus sequence.

All consensus sequences were manually inspected for accuracy, such as the presence of intact open reading frames, using BioEdit (http://www.mbio.ncsu.edu/bioedit/bioedit.html). A sequence alignment of the 170 newly sequenced genomes combined with the 151 complete genome sequences from GenBank was constructed using ClustalW2 with default parameters [75] (http://www.ebi.ac.uk/Tools/msa/clustalw2/) implemented in Galaxy and manually adjusted when necessary. Sequence alignments of individual RABV genes (N, P, M, G and L genes) and concatenated non-coding regions (from the stop codon in N to the initiation codon of L) were also generated. All the full-length genome sequences generated in the present study have been submitted to GenBank (S1 Table).

Phylogenetic analysis

We used jModelTest2 [76, 77] to determine the best-fit model of nucleotide substitution according to the Bayesian Information Criterion. This revealed that the general time reversible model with proportion of invariable sites plus gamma-distributed rate heterogeneity (GTR+I+Γ4) was optimal for all the RABV data sets compiled here. Phylogenetic trees using the different data sets (i.e. individual genes, concatenated genes or non-coding regions) were then estimated using the maximum likelihood (ML) method available in PhyML 3.0 [78] utilizing SPR branch-swapping. The robustness of individual nodes on the phylogeny was estimated using 1,000 bootstrap replicates for the five concatenated gene data set, and using the approximate likelihood ratio test (aLRT) with SH-like supports for each individual RABV gene as well as the concatenated non-coding region data set [79].

Estimates of RABV evolutionary dynamics and time-scale

To determine the degree of clock-like structure in each data set we employed root-to-tip linear regression as available in the TempEst program [80]. For those data sets with sufficient phylogenetic structure we then inferred a maximum clade credibility (MCC) tree using the Bayesian Markov chain Monte Carlo (MCMC) method available in the BEAST v1.8 package [81] by incorporating information on sampling time (year) of the dog-related RABV group (isolates for which the date of sampling was unavailable and vaccine strains were excluded). Posterior probability values provided an assessment of the degree of support for each node on the tree. This analysis utilized the GTR+I+Γ4 model of nucleotide substitution, a relaxed (uncorrelated log-normal) molecular clock and the constant population size model as a coalescent prior. Ten independent MCMC analyses were run for 100 million steps and sampled every 10,000 states. The log and tree files of each MCMC chains were combined using Logcombiner v1.8.2 (http://tree.bio.ed.ac.uk/software/beast/), with a burn-in of 10%. The convergence of each parameter in this combined file was checked using TRACER v1.6 (http://tree.bio.ed.ac.uk/software/tracer/) and indicated by an effective sample size >200. The MCC tree was obtained using TreeAnnotator v1.8.2 (http://tree.bio.ed.ac.uk/software/beast/). Additional analyses were performed utilizing the GMRF Bayesian Skyride [82] and Bayesian SkyGrid [83] demographic models, and gave similar results.

Based on the BEAST analysis, we also estimated the rate of nucleotide substitution per site, per year (see below) and the time of most recent common ancestor (TRMCA) for host-specific clusters of sequences. The degree of statistical uncertainty in each parameter estimate was given by the 95% highest posterior density (HPD) values.

The root-to-tip regression analysis performed on the 248 sequences of the dog-related RABV group revealed a number of clear outlier taxa characterized by anomalously high evolutionary rates (S4C Fig). These outliers belong to three clades or sub-clades: the Africa-3 clade that is specific to mongooses in Southern Africa, and the SEA5 and SEA2b subclades that are confined to viruses from ferret-badgers in Taiwan and China, respectively. To further assess if there are considerable differences in evolutionary rate in these clades, we performed additional analyses on the N and the G proteins for which relatively large numbers of sequences were available on GenBank. We therefore collected from GenBank an additional 41 N and 26 G sequences from the Africa-3 clade, and an additional 72 N and 71 G sequences from the ferret-badger in Taiwan and China (S8 Table). These data sets were compared to the N and G sequences of the dog-related RABV group (n = 248) and two RABV subsets corresponding to viruses circulating in dogs in Asia (n = 51) and in Africa (n = 46). As the Africa-3 clade is specific to the mongoose, we also estimated the evolutionary rate of RABV circulating in mongooses in the Caribbean region, for which we constructed a data set of 64 N sequences (no G sequences were available). As the ferret-badger data set was small, covered a relatively short time-range, and comprised two groups sampled during different time periods, it was unfortunately impossible to analyse the evolutionary dynamics in these two groups separately.

Estimates of nucleotide substitution rate of each data set were performed using BEAST as described above. Preliminary analysis on the N and G gene data sets using different nucleotide substitutions models (GTR+I+Γ4 or GTR+I), strict or relaxed (uncorrelated log-normal) molecular clocks, constant population size or Bayesian skyline coalescent priors gave similar results. Therefore, all analyses were performed using the GTR+I+Γ4 substitution model, a relaxed (uncorrelated log-normal) molecular clock, and a constant population size.

Finally, to assess the robustness of our rate estimates we also utilized and hierarchical phylogenetic models [46]. This analysis considered the lineages of the N and G genes defined previously (i.e. those viruses circulating in ferret-badgers, in mongooses in Southern Africa and the Caribbean, and in dogs in Africa and Asia) which we treated as data partitions. To be as robust as possible we used two substitution models–SRD06 [45] and GTR+I+Γ4. For the SRD06 model we specified hyperprior distributions to govern κ (the relative rate of transitions to transversions) and the shape parameter, α, of the Γ-distribution among the first and second, and third codon positions. In the case of the GTR+I+Γ4 model we linked each of the six rate parameters of the substitution matrix, α, and the proportion of invariable sites (I). Importantly, for both substitution models we set separate uncorrelated lognormal relaxed clock models and constant-size coalescent tree priors for each of the partitions, which is appropriate because they involve different taxa.

Analysis of selection pressures

To reveal the selection pressures acting on the RABV genome we compared the numbers of nonsynonymous (dN) and synonymous (dS) substitutions per site for the different RABV genes and phylogenetic clusters using the Single Likelihood Ancestor Counting (SLAC), Fixed Effect Likelihood (FEL), the internal branch Mixed Effects of Model Evolution (MEME-internal; Kosakovsky Pond SL, personal communication) and the Fast Unbiased Bayesian Approximation (FUBAR) models [4749]. Only codon positions with a p-value < 0.05 for the SLAC, FEL and MEME models and with a posterior of probability > 0.95 for the FUBAR method were considered as containing evidence for positive selection. For each data set and gene, the best-fit model of nucleotide substitution model was determined using the model selection tool available on the DATAMONKEY server [84, 85].

Supporting Information

S1 Text. Description of each dog-RABV clade, subclade and lineage identified in Fig 1 and their TRMCA estimates.

https://doi.org/10.1371/journal.ppat.1006041.s001

(DOCX)

S1 Fig. Geographic distribution of the 321 RABV isolates analyzed in the study.

(A) Triangles and dots represent the bat- and dog-related RABV, respectively, with sizes proportional to the number of isolates as indicated in the legend. (B) The different colors represent the major clades in the bat- and dog-related RABV groups.

https://doi.org/10.1371/journal.ppat.1006041.s002

(PDF)

S2 Fig. Comparison of the maximum likelihood phylogenies of 321 RABV sequences representing the N, P, M, G, L genes and concatenated non-coding regions.

The ML trees are mid-point rooted and aLRT values are shown for each clade, subclade and lineage named according to Fig 1 and the nucleoprotein (A), phosphoprotein (B), matrix (C), glycoprotein (D), polymerase (E) genes and concatenated non-coding regions (F) are shown separately.

https://doi.org/10.1371/journal.ppat.1006041.s003

(PDF)

S3 Fig. Maximum likelihood phylogeny of 67 bat-related RABV sequences representing five concatenated genes.

The tree is mid-point rooted with bootstrap values shown for each major subclades. The subclade names are given according to Kuzmin et al., (2012) [14]. Subclade abbreviations: SCSK–South-Central skunk; RAC–North-American Raccoon; MexSK-1 –Mexican skunk, variant 1; EF-W1 and EF-W2 –Eptesicus fuscus, in western USA; MYu–Myotis yumanensis; LX–Lasiurus xanthinus; LS–Lasiurus seminolus; LC–Lasiurus cinereus; LB–Lasiurus borealis; PS–Perimyotis subflavus; LN–Lasionycteris noctivagans; LI–Lasiurus intermedius; TB–Tadarida brasiliensis; DR–Desmodus rotundus; MYsp–Myotis spp; PH–Parastrellus hesperus; AP–Antrozous pallidus; EF–Eptesicus fuscus, in eastern and central USA.

https://doi.org/10.1371/journal.ppat.1006041.s004

(PDF)

S4 Fig. Root-to-tip regression of genetic distances against the year of sampling for five concatenated RABV genes.

The root-to-tip regressions were obtained using TempEst [80], on (A) a combined bat- and dog-related RABV data set (n = 315), (B) bat-related RABV data set (n = 67), and (C) dog-related RABV data set (n = 248) (isolates for which the date of sampling were unavailable and vaccine strains were excluded). The red circle indicates a number of outlier strains characterized by anomalously high rates (employed an arbitrary cut-off). The inferred rate of nucleotide substitution corresponding to the slope and the correlation coefficient are also indicated.

https://doi.org/10.1371/journal.ppat.1006041.s005

(PDF)

S5 Fig. Maximum likelihood tree of 321 RABV from the five concatenated genes labelled by host species.

Tip names are colored according to the isolation species of each virus, dog in black, bat in grey, ferret-badger in magenta, mongoose in orange, human in red, other carnivores in green, herbivores and/or omnivores in blue and vaccine strains in purple. The major clades of RABV are indicated in boxes like in Fig 1. The names of subclades and lineages defined for the Arctic-related, Asian and Cosmopolitan clades are detailed in S1 Table. The tree is mid-point rooted for clarity only.

https://doi.org/10.1371/journal.ppat.1006041.s006

(PDF)

S1 Table. List of viruses used in the full-length genome analyses.

https://doi.org/10.1371/journal.ppat.1006041.s007

(DOCX)

S2 Table. Bootstrap values of nodes corresponding to all clades, subclades and lineages of the dog-related group defined in Fig 1.

https://doi.org/10.1371/journal.ppat.1006041.s008

(DOCX)

S3 Table. Evolutionary characteristics of dog-related RABV group clades, subclades and lineages.

https://doi.org/10.1371/journal.ppat.1006041.s009

(DOCX)

S4 Table. Substitution rates of the N and G genes among different host species in the dog-related RABV group.

https://doi.org/10.1371/journal.ppat.1006041.s010

(DOCX)

S5 Table. Selection pressures in the N and G genes among different host species in the dog-related RABV group.

https://doi.org/10.1371/journal.ppat.1006041.s011

(DOCX)

S6 Table. Amino acid substitutions specific to mongoose-related RABV (Africa-3 clade) or ferret-badger-related RABV (SEA5 subclade and the SEA2b lineage).

https://doi.org/10.1371/journal.ppat.1006041.s012

(DOCX)

S7 Table. List of primers used in this study.

https://doi.org/10.1371/journal.ppat.1006041.s013

(DOCX)

S8 Table. List of additional nucleotide and glycoprotein sequences used to estimate the evolution rates among different hosts.

https://doi.org/10.1371/journal.ppat.1006041.s014

(DOCX)

Acknowledgments

We thank Laurence Ma and Magali Tichit from the Genomics Platform and Andreea Alexandru from the Mutualized Platfom for Microbiology for technical assistance on the NGS sequencing. We thank Alexis Criscuolo and Julien Guglielmini from the Mutualized Platfom for Microbiology of the Pasteur International Bioresources Network (PIBnet) and Nizar Fawal for their advices in bio-informatics analyses. We thank Sergei L. Kosakovsky Pond for his help on the MEME-internal analysis.

Author Contributions

  1. Conceived and designed the experiments: CT LD ECH HBo.
  2. Performed the experiments: CT MT HBl CB ECH.
  3. Analyzed the data: CT LD MT SD ECH HBo.
  4. Contributed reagents/materials/analysis tools: CT LD CS HBl CB MV HBo.
  5. Wrote the paper: CT LD CS CB MV SD ECH HBo.

References

  1. 1. Holmes EC. Evolutionary history and phylogeography of human viruses. Annual review of microbiology. 2008;62:307–28. pmid:18785840
  2. 2. Kuiken T, Fouchier R, Rimmelzwaan G, Osterhaus A. Emerging viral infections in a rapidly changing world. Current opinion in biotechnology. 2003;14(6):641–6. pmid:14662395
  3. 3. Cleaveland S, Haydon DT, Taylor L. Overviews of pathogen emergence: which pathogens emerge, when and why? Current topics in microbiology and immunology. 2007;315:85–111. pmid:17848062
  4. 4. Cleaveland S, Laurenson MK, Taylor LH. Diseases of humans and their domestic mammals: pathogen characteristics, host range and the risk of emergence. Philosophical transactions of the Royal Society of London Series B, Biological sciences. 2001;356(1411):991–9. pmid:11516377
  5. 5. Childs JE, Richt JA, Mackenzie JS. Introduction: conceptualizing and partitioning the emergence process of zoonotic viruses from wildlife to humans. Current topics in microbiology and immunology. 2007;315:1–31. pmid:17848058
  6. 6. (ICTV) ICoToV. Virus Taxonomy: 2015 Release 2016. http://www.ictvonline.org/virustaxonomy.asp.
  7. 7. Steinhauer DA, Holland JJ. Rapid evolution of RNA viruses. Annual review of microbiology. 1987;41:409–33. pmid:3318675
  8. 8. Badrane H, Tordo N. Host switching in Lyssavirus history from the Chiroptera to the Carnivora orders. Journal of virology. 2001;75(17):8096–104. pmid:11483755
  9. 9. Faria NR, Suchard MA, Rambaut A, Streicker DG, Lemey P. Simultaneously reconstructing viral cross-species transmission history and identifying the underlying constraints. Philosophical transactions of the Royal Society of London Series B, Biological sciences. 2013;368(1614):20120196. pmid:23382420
  10. 10. Holmes EC, Woelk CH, Kassis R, Bourhy H. Genetic constraints and the adaptive evolution of rabies virus in nature. Virology. 2002;292(2):247–57. pmid:11878928
  11. 11. Rupprecht CE, Turmelle A, Kuzmin IV. A perspective on lyssavirus emergence and perpetuation. Curr Opin Virol. 2011;1(6):662–70. pmid:22440925
  12. 12. Nel LH, Markotter W. Lyssaviruses. Critical reviews in microbiology. 2007;33(4):301–24. pmid:18033596
  13. 13. Bourhy H, Reynes JM, Dunham EJ, Dacheux L, Larrous F, Huong VT, et al. The origin and phylogeography of dog rabies virus. The Journal of general virology. 2008;89(Pt 11):2673–81. pmid:18931062
  14. 14. Kuzmin IV, Shi M, Orciari LA, Yager PA, Velasco-Villa A, Kuzmina NA, et al. Molecular inferences suggest multiple host shifts of rabies viruses from bats to mesocarnivores in Arizona during 2001–2009. PLoS pathogens. 2012;8(6):e1002786. pmid:22737076
  15. 15. Biek R, Henderson JC, Waller LA, Rupprecht CE, Real LA. A high-resolution genetic signature of demographic and spatial expansion in epizootic rabies virus. Proceedings of the National Academy of Sciences of the United States of America. 2007;104(19):7993–8. pmid:17470818
  16. 16. Davis R, Nadin-Davis SA, Moore M, Hanlon C. Genetic characterization and phylogenetic analysis of skunk-associated rabies viruses in North America with special emphasis on the central plains. Virus research. 2013;174(1–2):27–36. pmid:23524137
  17. 17. Kuzmina NA, Lemey P, Kuzmin IV, Mayes BC, Ellison JA, Orciari LA, et al. The phylogeography and spatiotemporal spread of south-central skunk rabies virus. PloS one. 2013;8(12):e82348. pmid:24312657
  18. 18. Bourhy H, Kissi B, Audry L, Smreczak M, Sadkowska-Todys M, Kulonen K, et al. Ecology and evolution of rabies virus in Europe. The Journal of general virology. 1999;80 (Pt 10):2545–57.
  19. 19. Nel LH, Sabeta CT, von Teichman B, Jaftha JB, Rupprecht CE, Bingham J. Mongoose rabies in southern Africa: a re-evaluation based on molecular epidemiology. Virus research. 2005;109(2):165–73. pmid:15763147
  20. 20. Oem JK, Kim SH, Kim YH, Lee MH, Lee KK. Complete genome sequences of three rabies viruses isolated from rabid raccoon dogs and a cow in Korea. Virus genes. 2013;47(3):563–8. pmid:23975690
  21. 21. Tsai KJ, Hsu WC, Chuang WC, Chang JC, Tu YC, Tsai HJ, et al. Emergence of a sylvatic enzootic formosan ferret badger-associated rabies in Taiwan and the geographical separation of two phylogenetic groups of rabies viruses. Veterinary microbiology. 2016;182:28–34. pmid:26711025
  22. 22. Zhao J, Liu Y, Zhang S, Zhang F, Wang Y, Mi L, et al. Molecular characterization of three ferret badger (Melogale moschata) rabies virus isolates from Jiangxi province, China. Archives of virology. 2014;159(8):2059–67. pmid:24643334
  23. 23. Hampson K, Coudeville L, Lembo T, Sambo M, Kieffer A, Attlan M, et al. Estimating the global burden of endemic canine rabies. PLoS neglected tropical diseases. 2015;9(4):e0003709. pmid:25881058
  24. 24. Nadin-Davis SA, Abdel-Malik M, Armstrong J, Wandeler AI. Lyssavirus P gene characterisation provides insights into the phylogeny of the genus and identifies structural similarities and diversity within the encoded phosphoprotein. Virology. 2002;298(2):286–305. pmid:12127791
  25. 25. Pant GR, Lavenir R, Wong FY, Certoma A, Larrous F, Bhatta DR, et al. Recent emergence and spread of an Arctic-related phylogenetic lineage of rabies virus in Nepal. PLoS neglected tropical diseases. 2013;7(11):e2560. pmid:24278494
  26. 26. Streicker DG, Altizer SM, Velasco-Villa A, Rupprecht CE. Variable evolutionary routes to host establishment across repeated rabies virus host shifts among bats. Proceedings of the National Academy of Sciences of the United States of America. 2012;109(48):19715–20. pmid:23150575
  27. 27. Talbi C, Holmes EC, de Benedictis P, Faye O, Nakoune E, Gamatie D, et al. Evolutionary history and dynamics of dog rabies virus in western and central Africa. The Journal of general virology. 2009;90(Pt 4):783–91. pmid:19264663
  28. 28. Zhao J, Zhang S, Liu Y, Zhang F, Hu R. Complete Genome Sequence of a Rabies Virus Isolate from a Ferret Badger (Melogale moschata) in Jiangxi, China. Genome Announc. 2013;1(3).
  29. 29. Zieger U, Marston DA, Sharma R, Chikweto A, Tiwari K, Sayyid M, et al. The phylogeography of rabies in Grenada, West Indies, and implications for control. PLoS neglected tropical diseases. 2014;8(10):e3251. pmid:25330178
  30. 30. Brunker K, Marston DA, Horton DL, Cleaveland S, Fooks AR, Kazwala R, et al. Elucidating the phylogynamics of endemic rabies virus in eastern Africa using whole-genome sequencing. Virus Evolution. 2015;1(1):vev011. pmid:27774283
  31. 31. Delmas O, Holmes EC, Talbi C, Larrous F, Dacheux L, Bouchier C, et al. Genomic diversity and evolution of the lyssaviruses. PloS one. 2008;3(4):e2057. pmid:18446239
  32. 32. Tang HB, Lu ZL, Zhong YZ, He XX, Zhong TZ, Pan Y, et al. Characterization of the biological properties and complete genome sequence analysis of a cattle-derived rabies virus isolate from the Guangxi province of southern China. Virus genes. 2014;49(3):417–27. pmid:25142164
  33. 33. Zhang J, Zhang HL, Tao XY, Li H, Tang Q, Jiang XY, et al. The full-length genome analysis of a street rabies virus strain isolated in Yunnan province of China. Virologica Sinica. 2012;27(3):204–13. pmid:22684475
  34. 34. Talbi C, Lemey P, Suchard MA, Abdelatif E, Elharrak M, Nourlil J, et al. Phylodynamics and human-mediated dispersal of a zoonotic virus. PLoS pathogens. 2010;6(10):e1001166. pmid:21060816
  35. 35. Bourhy H, Nakoune E, Hall M, Nouvellet P, Lepelletier A, Talbi C, et al. Revealing the Micro-scale Signature of Endemic Zoonotic Disease Transmission in an African Urban Setting. PLoS pathogens. 2016;12(4):e1005525. pmid:27058957
  36. 36. Streicker DG, Lemey P, Velasco-Villa A, Rupprecht CE. Rates of viral evolution are linked to host geography in bat rabies. PLoS pathogens. 2012;8(5):e1002720. pmid:22615575
  37. 37. Fourment M, Holmes EC. Avian influenza virus exhibits distinct evolutionary dynamics in wild birds and poultry. BMC evolutionary biology. 2015;15:120. pmid:26111936
  38. 38. Fooks AR, Banyard AC, Horton DL, Johnson N, McElhinney LM, Jackson AC. Current status of rabies and prospects for elimination. Lancet. 2014;384(9951):1389–99. pmid:24828901
  39. 39. Velasco-Villa A, Orciari LA, Souza V, Juarez-Islas V, Gomez-Sierra M, Castillo A, et al. Molecular epizootiology of rabies associated with terrestrial carnivores in Mexico. Virus research. 2005;111(1):13–27. pmid:15896399
  40. 40. Chang JC, Tsai KJ, Hsu WC, Tu YC, Chuang WC, Chang CY, et al. Rabies Virus Infection in Ferret Badgers (Melogale moschata subaurantiaca) in Taiwan: A Retrospective Study. Journal of wildlife diseases. 2015;51(4):923–8. pmid:26267459
  41. 41. Chiou HY, Hsieh CH, Jeng CR, Chan FT, Wang HY, Pang VF. Molecular characterization of cryptically circulating rabies virus from ferret badgers, Taiwan. Emerging infectious diseases. 2014;20(5):790–8. pmid:24751120
  42. 42. Davis PL, Rambaut A, Bourhy H, Holmes EC. The evolutionary dynamics of canid and mongoose rabies virus in Southern Africa. Archives of virology. 2007;152(7):1251–8. pmid:17401615
  43. 43. Van Zyl N, Markotter W, Nel LH. Evolutionary history of African mongoose rabies. Virus research. 2010;150(1–2):93–102. pmid:20214938
  44. 44. Gong W, Jiang Y, Za Y, Zeng Z, Shao M, Fan J, et al. Temporal and spatial dynamics of rabies viruses in China and Southeast Asia. Virus research. 2010;150(1–2):111–8. pmid:20214936
  45. 45. Shapiro B, Rambaut A, Drummond AJ. Choosing appropriate substitution models for the phylogenetic analysis of protein-coding sequences. Molecular biology and evolution. 2006;23(1):7–9. pmid:16177232
  46. 46. Suchard MA, Kitchen CM, Sinsheimer JS, Weiss RE. Hierarchical phylogenetic models for analyzing multipartite sequence data. Systematic biology. 2003;52(5):649–64. pmid:14530132
  47. 47. Kosakovsky Pond SL, Frost SD. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Molecular biology and evolution. 2005;22(5):1208–22. pmid:15703242
  48. 48. Murrell B, Moola S, Mabona A, Weighill T, Sheward D, Kosakovsky Pond SL, et al. FUBAR: a fast, unconstrained bayesian approximation for inferring selection. Molecular biology and evolution. 2013;30(5):1196–205. pmid:23420840
  49. 49. Murrell B, Wertheim JO, Moola S, Weighill T, Scheffler K, Kosakovsky Pond SL. Detecting individual sites subject to episodic diversifying selection. PLoS genetics. 2012;8(7):e1002764. pmid:22807683
  50. 50. Guo Z, Tao X, Yin C, Han N, Yu J, Li H, et al. National borders effectively halt the spread of rabies: the current rabies epidemic in China is dislocated from cases in neighboring countries. PLoS neglected tropical diseases. 2013;7(1):e2039. pmid:23383359
  51. 51. Meng S, Sun Y, Wu X, Tang J, Xu G, Lei Y, et al. Evolutionary dynamics of rabies viruses highlights the importance of China rabies transmission in Asia. Virology. 2011;410(2):403–9. pmid:21195445
  52. 52. Ming P, Du J, Tang Q, Yan J, Nadin-Davis SA, Li H, et al. Molecular characterization of the complete genome of a street rabies virus isolated in China. Virus research. 2009;143(1):6–14. pmid:19463716
  53. 53. Wu H, Wang L, Tao X, Li H, Rayner S, Liang G, et al. Genetic diversity and molecular evolution of the rabies virus matrix protein gene in China. Infection, genetics and evolution: journal of molecular epidemiology and evolutionary genetics in infectious diseases. 2013;16:248–53. pmid:23453987
  54. 54. Wertheim JO, Kosakovsky Pond SL. Purifying selection can obscure the ancient age of viral lineages. Molecular biology and evolution. 2011;28(12):3355–65. pmid:21705379
  55. 55. Assenberg R, Delmas O, Ren J, Vidalain PO, Verma A, Larrous F, et al. Structure of the nucleoprotein binding domain of Mokola virus phosphoprotein. Journal of virology. 2010;84(2):1089–96. pmid:19906936
  56. 56. Mavrakis M, McCarthy AA, Roche S, Blondel D, Ruigrok RW. Structure and function of the C-terminal domain of the polymerase cofactor of rabies virus. J Mol Biol. 2004;343(4):819–31. pmid:15476803
  57. 57. Steele JH, Fernandez PJ. History of rabies and global aspects. In: Baer GM, editor. The Natural History of Rabies. 2. Boca Raton, USA: CRC Press; 1991. p. 1–26.
  58. 58. Theodorides J. Histoire de la Rage. Paris: Masson; 1986. 289 p.
  59. 59. History APW. Key Concept 4.1 Globalizing Networks of Communication and Exchange. http://apworldipedia.com/index.php?title=Key_Concept_4.1_Globalizing_Networks_of_Communication_and_Exchange.
  60. 60. Streicker DG, Turmelle AS, Vonhof MJ, Kuzmin IV, McCracken GF, Rupprecht CE. Host phylogeny constrains cross-species emergence and establishment of rabies virus in bats. Science. 2010;329(5992):676–9. pmid:20689015
  61. 61. Villarreal LP, Defilippis VR, Gottlieb KA. Acute and persistent viral life strategies and their relationship to emerging diseases. Virology. 2000;272(1):1–6. pmid:10873743
  62. 62. Engeman RM, Christensen KL, Pipas MJ, Bergman DL. Population monitoring in support of a rabies vaccination program for skunks in Arizona. Journal of wildlife diseases. 2003;39(3):746–50. pmid:14567243
  63. 63. Leslie MJ, Messenger S, Rohde RE, Smith J, Cheshier R, Hanlon C, et al. Bat-associated rabies virus in Skunks. Emerging infectious diseases. 2006;12(8):1274–7. pmid:16965714
  64. 64. Mollentze N, Biek R, Streicker DG. The role of viral evolution in rabies host shifts and emergence. Curr Opin Virol. 2014;8:68–72. pmid:25064563
  65. 65. Borucki MK, Chen-Harris H, Lao V, Vanier G, Wadford DA, Messenger S, et al. Ultra-deep sequencing of intra-host rabies virus populations during cross-species transmission. PLoS neglected tropical diseases. 2013;7(11):e2555. pmid:24278493
  66. 66. Parrish CR, Holmes EC, Morens DM, Park EC, Burke DS, Calisher CH, et al. Cross-species virus transmission and the emergence of new epidemic diseases. Microbiol Mol Biol Rev. 2008;72(3):457–70. pmid:18772285
  67. 67. Woolhouse ME, Haydon DT, Antia R. Emerging pathogens: the epidemiology and evolution of species jumps. Trends Ecol Evol. 2005;20(5):238–44. pmid:16701375
  68. 68. Sanjuan R, Cuevas JM, Moya A, Elena SF. Epistasis and the adaptability of an RNA virus. Genetics. 2005;170(3):1001–8. pmid:15879507
  69. 69. Wiltzer L, Okada K, Yamaoka S, Larrous F, Kuusisto HV, Sugiyama M, et al. Interaction of rabies virus P-protein with STAT proteins is critical to lethal rabies disease. J Infect Dis. 2014;209(11):1744–53. pmid:24367042
  70. 70. Dacheux L, Berthet N, Dissard G, Holmes EC, Delmas O, Larrous F, et al. Application of broad-spectrum resequencing microarray for genotyping rhabdoviruses. Journal of virology. 2010;84(18):9557–74. pmid:20610710
  71. 71. Criscuolo A, Brisse S. AlienTrimmer: a tool to quickly and accurately trim off multiple short contaminant sequences from high-throughput sequencing reads. Genomics. 2013;102(5–6):500–6. pmid:23912058
  72. 72. Goecks J, Nekrutenko A, Taylor J, Galaxy T. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome biology. 2010;11(8):R86. pmid:20738864
  73. 73. Blankenberg D, Von Kuster G, Coraor N, Ananda G, Lazarus R, Mangan M, et al. Galaxy: a web-based genome analysis tool for experimentalists. Current protocols in molecular biology / edited by Ausubel Frederick M [et al]. 2010;Chapter 19:Unit 19 0 1–21.
  74. 74. Giardine B, Riemer C, Hardison RC, Burhans R, Elnitski L, Shah P, et al. Galaxy: a platform for interactive large-scale genome analysis. Genome research. 2005;15(10):1451–5. pmid:16169926
  75. 75. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23(21):2947–8. pmid:17846036
  76. 76. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.
  77. 77. Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Systematic biology. 2003;52(5):696–704. pmid:14530136
  78. 78. Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Systematic biology. 2010;59(3):307–21. pmid:20525638
  79. 79. Anisimova M, Gascuel O. Approximate likelihood-ratio test for branches: A fast, accurate, and powerful alternative. Systematic biology. 2006;55(4):539–52. pmid:16785212
  80. 80. Rambaut A, Lam TT, L.M. C, Pybus OG. Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). Virus Evolution. 2016;2.
  81. 81. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC evolutionary biology. 2007;7:214. pmid:17996036
  82. 82. Minin VN, Bloomquist EW, Suchard MA. Smooth skyride through a rough skyline: Bayesian coalescent-based inference of population dynamics. Molecular biology and evolution. 2008;25(7):1459–71. pmid:18408232
  83. 83. Gill MS, Lemey P, Faria NR, Rambaut A, Shapiro B, Suchard MA. Improving Bayesian population dynamics inference: a coalescent-based model for multiple loci. Molecular biology and evolution. 2013;30(3):713–24. pmid:23180580
  84. 84. Delport W, Poon AF, Frost SD, Kosakovsky Pond SL. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics. 2010;26(19):2455–7. pmid:20671151
  85. 85. Pond SL, Frost SD. Datamonkey: rapid detection of selective pressure on individual sites of codon alignments. Bioinformatics. 2005;21(10):2531–3. pmid:15713735