Genomic Diversity and Evolution of the Lyssaviruses

1Lyssaviruses (LYSSAV) are RNA viruses with single-stranded, negative-sense genomes of the family Rhabdoviridae that infect a variety of mammals causing rabies-like diseases. Rabies is an ancient disease that may have been reported in the Old World before 2300 B.C.. However, the absence of effective control measures in animal reservoir populations combined with a widespread lack of human access to vaccination means that more than 50,000 people annually die of rabies, particularly in Asia and Africa,. Currently, there are seven recognised genotypes (GT) of LYSSAV defined on the basis of their genetic similarity,: rabies virus (RABV, GT1) responsible for classical rabies in terrestrial mammals globally and in bats on the American continent, as well as the cause of most rabies-related human deaths worldwide; Lagos bat virus (LBV, GT2); Mokola virus (MOKV, GT3); Duvenhage virus (DUVV, GT4); European bat lyssavirus type 1 (EBLV-1, GT5); European bat lyssavirus type 2 (EBLV-2, GT 6); and Australian bat lyssavirus (ABLV, GT7). All genotypes except MOKV (where the host species is unknown) have bat reservoirs, hinting that lyssaviruses originated in these mammals. Additionally, four new lyssavirus genotypes that infect bats in central and southeast Asia have been proposed: Aravan virus, Khujand virus, Irkut virus and West Caucasian Bat virus,. The negative-sense LYSSAV genome encodes five proteins: the nucleoprotein (N), phosphoprotein (P), matrix protein (M), glycoprotein (G) and RNA polymerase (L) in the order 3′-N-P-M-G-L-5′.

2Despite the importance of LYSSAV for human and wildlife populations, the number of complete genome sequences of field isolates of LYSSAV is sparse, with only eight currently available for limited type species –. Herein, we present the first genomic and evolutionary analysis of the seven known genotypes of LYSSAV, therein significantly increasing the extent of available genome sequence data available for these important mammalian pathogens.

3Total RNA (Table 1) was isolated from original specimens or from suckling mice brain after early passage using Tri-Reagent (Euromedex). The only exception was the 8743THA isolate that was adapted on BSR cells (passage 22). For this isolate, total RNA was isolated from infected BSR cells infected at a low multiplicity of infection (0,1). Reverse transcription was performed with random hexamer primer (Roche Boehringer) using Superscript II (Invitrogen) following the manufacturer instructions.

4Long-range PCR products were obtained using ExTaq (Takara) and specific primers (Table S1) using manufacturer recommendations. For sequence determination we used a shotgun base approach called LoPPS (Long PCR Product Sequencing),. 3′ genomic ends were generated by RACE protocol using a 5′ phosphorylated reverse complementary T7 primer. T7 cDNAs were further used for heminested-PCR with ExTaq using T7 and two strain specific primers designed in the N coding region (supplementary Table 1). To determine the 5′ sequence of the genomic RNA we used a 5′RACE version 2.0 kit from Invitrogen following manufacturer instructions. The PCR products (5′ or 3′ RACE) were then purified on gel using Qiaquick gel extraction kit (Qiagen) and cloned in PCR 2.1 TOPO T/A (Invitrogen) for sequencing. Each position of the consensus nucleotide sequence was determined from at least three independent sequences. All consensus sequences obtained using Sequencher 4.7 (Gene Codes) software were aligned using ClustalX 1.83.1. The untranslated regions were further aligned manually using the SE-AL program (http://tree.bio.ed.ac.uk/). GenBank accession numbers for the sequences newly acquired here are designated EU293108-EU293121.

5Phylogenetic analysis of LYSSAV genomes was based on a multiple alignment of concatenated coding region sequences (12105 nt). A maximum likelihood (ML) phylogenetic analysis of these data was undertaken using PAUP* employing the best-fit GTR+I+Γ4 model of nucleotide substitution inferred by ModelTest. To determine the extent of support for different groupings on the tree a bootstrap resampling analysis was undertaken employing 1000 replicate neighbor-joining trees estimated under the ML substitution model.

6In total we determined 14 new complete genome sequences of field isolates representing six (GT1, GT2, GT3, GT4, GT5 and GT6) of the seven genotypes of LYSSAV, with complete genome from GT4 obtained for the first time. These genomes were combined with eight genomes described previously (with the exception of one Australian bat lyssavirus for which leader and trailer sequences are unavailable). Eight field isolates of viruses isolated from humans, canids and bats were chosen as representative of the diversity of GT1. Two vaccine strains (SAD-B19 and PV) were included in all sequence comparisons but not in the phylogenetic analysis. Our study also represent the first analysis of the intrinsic genetic diversity of GT2, GT3, GT4, GT5 and GT6 based on full length genomes.

7All genomes have the same structural organization although their lengths varied between 11918 nt. (GT7) and 12016 nt. (GT2) (Table 2). The predicted size of the coding regions is similar among genotypes, with the M protein identical in length across all genotypes and the P protein the most variable,,. As observed in other RNA viruses, all genotypes show a bias toward G+C richness, with the lowest G+C content observed in GT2 and the highest in GT1 (Table 2). All genomes have a polycistronic genome organization surrounded by untranslated regions (Table S2) similar to that already described,. The extent of genetic diversity, reflected in percentage identity, varies within and among proteins (Figure 1), in the order N>L>M>G>P (95.2, 94.2, 92.3, 85.8, 81.5% amino acid identity, respectively). A similar pattern was previously observed using more limited data sets. This same order was also observed in terms of overall selection pressure, measured as the mean ratio of nonsynonymous (dN) to synonymous substitutions (dS) per site (dN/dS), estimated using the maximum likelihood SLAC (Single Likelihood Ancestor Counting; http://www.datamonkey.org/) method: N = 0.048; L = 0.055; M = 0.078; G = 0.119; P = 0.187. This approximately four-fold difference in mean dN/dS reflects major differences in selective constraint among proteins. This trend was also reflected in previous analyses of full length genomes of vaccine strains and through partial gene comparisons,.

8Our study represents the largest analysis of the 3′ and 5′ UTR of the lyssavirus genomes undertaken to date. The 3′ UTR comprises 70 nt and includes the leader regions potentially transcribed into the leader RNA. The 5′UTR region comprises 86–145 nt and contains the trailer regions of size 68–69 nt. Both the 3′ and 5′ UTR have conserved signals that play a role to modulate replication and transcription (Figure S1). Our data also reveals a strict complementarity limited to the 9 terminal nucleotides as well as nucleotide positions 14 and 16 from both ends of the genome,,.

9There have been several attempts to estimate the evolutionary relationships among lyssaviruses, with most utilizing only one or two genes,,,, –. We therefore undertook a phylogenetic analysis of 22 genomes representative of the seven genotypes of LYSSAV based on a multiple alignment of concatenated coding sequences. Our phylogenetic analysis reveals the separation of LYSSAV into two major branches previously defined as different ‘phylogroups’ and 7 component lineages defined as genotypes,. Phylogroup 1 comprised GT1, 4, 5, 6 and 7, while phylogroup 2 contains only GT2 and GT3 (Figure 2). Notably, phylogroup 2 contains viruses of sampled exclusively from Africa – LBV and MOKV – while a third African genotype (DUVV) is found within phylogroup 1,. Also of note was the observation that although GT5 and GT6 both circulate in European insectivorous bats, the former is more closely related to the African GT4 viruses. Hence, there has clearly been an independent origin of genotypes 5 and 6 in European bats, as previously documented in analyses of the N and G genes in isolation. Finally, that bats appear as the principle host species across such a large phylogeographic range indicates that the association between lyssaviruses and bats is likely to be the ancestral condition (with a secondary loss of bat transmission in GT3), such that the movement of bats is likely to be responsible for the global dissemination of these viruses.

10Notably, our study represents the first analysis of the genetic diversity of four complete genomes of GT2 and GT4, both of which are African in origin. While little variation is seen within GT4, the degree of divergence among the two GT2 isolates (Lagos bat virus – LBV) is striking (23.7% and 12.1% at the nucleotide and amino acid levels, respectively) and greater than that seen within any other genotype. Hence, although 0406SEN and 8619NGA are related according to the arbitrary classification system based on nucleotide identity between N coding regions (80.3% between 0406SEN and 8619NGA compared to a cut-off of 80%),, this classification system will likely need to be revised as expanded surveys of LYSSAV in Africa (this study and ) and in Eurasia,, reveal greater genetic diversity. More fundamentally, if, as we suggest, complete genomes represent the best tools for genotyping, we propose that 0406SEN should constitute a new GT8 different from GT2 (8619NGA) and that the genotype division should be set at 76.4 to 81.6% nucleotide identity at coding sequences for all five viral proteins. Such a cut-off would provide more discriminatory power than systems that utilize the N gene in isolation (Table 3).

11Finally, we suggest that the phylogenetic methods used here – based on a realistic model of nucleotide substitution, a robust phylogenetic method, and rigorous bootstrap resampling – represent a more powerful method of lyssavirus classification than those based on pairwise genetic diversity alone, particularly as they account for any lineage-specific rate variation that will compromise all distance-based approaches used to date. This method has also been proposed for HIV to try to standardize viral classification confirming the interest of this method for viral classification.