Rabies virus diversification in aerial and terrestrial mammals

Abstract Rabies is a fatal zoonotic infection of the central nervous system of mammals and has been known to humans for millennia. The etiological agent, is a neurotropic RNA virus in the order Mononegavirales, family Rhabdoviridae, genus Lyssavirus. There are currently accepted to be two cycles for rabies transmission: the urban cycle and the sylvatic cycle. The fact that both cycles originated from a common RABV or lyssavirus ancestor and the adaptive divergence that occurred since then as this ancestor virus adapted to a wide range of fitness landscapes represented by reservoir species in the orders Carnivora and Chiroptera led to the emergence of the diverse RABV lineages currently found in the sylvatic and urban cycles. Here we study full genome phylogenies and the time to the most recent common ancestor (TMRCA) of the RABVs in the sylvatic and urban cycles. Results show that there were differences between the nucleotide substitution rates per site per year for the same RABV genes maintained independently in the urban and sylvatic cycles. The results identify the most suitable gene for phylogenetic analysis, heterotachy among RABV genes and the TMRCA for the two cycles.

RABV is one of seventeen members of the genus Lyssavirus (King, A.M.K;Adans, M.J;Carstens, E.B;Lefkowitz 2012), fifteen of which have members of Chiroptera as exclusive reservoirs, showing the importance of this order as a reservoir for the genus (Rupprecht et al., 2017;Hu et al., 2018). This lends support to the hypothesis that various existing RABV lineages had an RABV specific to bats as a common ancestor (Badrane and Tordo, 2001;Hayman et al., 2016;Velasco-Villa et al., 2017) and that this ancestral lineage differentiated into the various lineages in the urban and sylvatic cycles (Oliveira et al., 2010;Streicker et al., 2010Streicker et al., , 2012aKuzmin et al., 2012). However, failure to isolate RABV or antibodies against this virus in bats in the Old World, particularly Asia, where the most ancient representatives of current RABV lineages are believed to be found, would appear to contradict this theory (Kuzmin and Rupprecht, 2007).
There are known to be two rabies transmission cycles: the urban cycle, in which the dog is the main reservoir and transmits the virus to other dogs, other domestic animals and man; and the sylvatic cycle, which is maintained by different terrestrial mammals and chiropterans (Acha andSzyfres, 2003). al., 2012), which could be the basis for classification of rabies cycles into bat-related and dog-related cycles. The bat-related cycle (currently the sylvatic cycle) can be defined as all the RABV lineages maintained in different independent epidemiologic cycles exclusively in the Americas by wild animals, particularly chiropterans, while the dog-related cycle corresponds to all the lineages of rabies circulating in dogs (currently the urban cycle) and wild canids (currently the sylvatic cycle) maintained independently in multiple cycles around the world.
At least 35 RABV lineages have been described in the bat-related cycle, of which 31 are maintained in chiropterans. Of the four remaining lineages, two are maintained in skunks in Mexico and the south of the USA, one in racoons in the USA and one in the common marmoset (C. jacchus) in the northeast of Brazil. Eight major groups of genetic lineages are described in the dog-related cycle, for which the only reservoirs are members of the order Carnivora, particularly the domestic dog and wild canids (Bourhy et al., 2008;Kuzmin et al., 2012). Of these eight groups identified in the dog-related cycle, only the Africa 3 complex, for which the reservoir is the African yellow mongoose, is not maintained in dogs.
In view of the limited number of complete RABV genomes published, the uncertainty surrounding the classification of rabies transmission cycles and the lack of robust evolutionary analyses for RABV, our objective in this study was to generate data and understand the various divergence events for RABV isolated from terrestrial mammals and bats based on estimates of phylogeny, evolution rates and time to the most recent common ancestor (TMRCA) using complete genome sequences.
This work complies with Protocol nº 2030/2010 issued by the "Ethics Committee in the use of animals" of the School of Veterinary Medicine and Animal Science of University of São Paulo.

Full-genome RT-PCR
Total viral RNA was extracted from infected mouse brain with TRIzol Ô (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions.
For each sample, 15 mL of the extracted RNA was added to the mix for reverse transcription containing 8 mL 5X RT Buffer, 6 mL of the dNTP pool at 10 mM, 5 mL of the "Final" antisense primer at 10 mM , 400 U of Reverse Transcriptase enzyme RevertAid Premium Reverse Transcriptase (Fermentas ), 1 mL of RNAsin (Invitrogen ) and 35 mL of ultra-pure DNAse/RNAse free water to a final volume of 50 mL, carrying out reverse transcription at 42 ºC / 180 minutes followed by a 70 ºC / 15 minute reverse transcriptase inactivation cycle.
Genomic RABV amplicons (circa 12 kb) were obtained by PCR using a pair of primers targeting the leader and trailer regions of the RABV genome , 1 mL of genomic cDNA and GoTaqÔ Long PCR Master Mix 2X (Promega) according to the manufacturer's instructions (Oliveira 2014).
The PCR protocol for the amplification of the complete RABV genome consisted of: 25 mL of GoTaq® Long PCR Master Mix 2X (Promega), 1 mL of the "Início" (sense) and "Final" (antisense) primers at a concentration of 10 mM, 1 mL of c-DNA and 22 mL of ultra-pure DNAse/RNAse free water to a final volume of 50 mL and taken to the thermocycler and subjected to the cycle described in Table 1.

Whole viral genome library preparation and full-genome sequencing
Purified products were quantitated using Quant-IT HS reagents (Invitrogen, Life Technologies, Carlsbad, CA), and approximately 1 ng of each was used in a fragmentation reaction mix using a Nextera XT DNA sample prep kit according to the manufacturer's protocol. Briefly, tagmentation and fragmentation of each product were simultaneously performed by incubation for 5 min at 55°C followed by incubation in neutralizing tagment buffer for 5 min at room temperature. After neutralization of the fragmented DNA, a light 12-cycle PCR was performed with Illumina Ready Mix to 2 Oliveira et al. add Illumina flowcell adaptors, indexes and common adapters for subsequent cluster generation and sequencing. Amplified DNA was then purified using Agencourt AMPure XP beads (Beckman Coulter), which excluded very short library fragments. Following AMPure purification, the quantity of each library was normalized to ensure equally library representation in our pooled samples. Prior to cluster generation, normalized libraries were further quantified by qPCR using the SYBR fast Illumina library quantification kit (KAPA Biosystems) following the instructions of the manufacturer. The qPCR was run on the 7500 Fast Real-Time PCR System (Applied Biosystems). The thermocycling conditions consisted of an initial denaturation step at 95°C for 5 min followed by 35 cycles of [30 s at 95°C and 45 s at 60°C]. The final libraries were pooled at equimolar concentration and diluted to 4 nM. To denature the indexed DNA, 5 mL of the 4 nM library were mixed with 5 mL of 0.2 N fresh NaOH and incubated for 5 min at room temperature. 990 mL of chilled Illumina HT1 buffer was added to the denatured DNA and mixed to make a 20 pM library. After this step, 360 mL of the 20 pM library was multiplexed with 6 mL of 12.5 pM denatured PhiX control to increase sequence diversity and then mixed with 234 mL of chilled HT1 buffer to make a 12 pM sequenceable library. Finally, 600 mL of the prepared library was loaded on an Illumina MiSeq clamshell style cartridge for paired end 250 sequencing. After they had been sequenced in an Illumina MiSeq TM , the reads were assembled without the use of reference genomes (de novo assembly) in CLC Genomics Workbench 6. When the assembled contigs in a genome were smaller than expected based on the size of the sequenced PCR fragment, they were extended by mapping the reads to the reference with 95% similarity and without gap opening or by using de novo assembly with the same similarity criterion in Geneious 6. Both the extended and the de novo assembled contigs had their reads mapped to the reference again and were evaluated in terms of their variability based on the number and quality of reads using the Quality-based Variant Detection function in CLC Genomics Workbench version 5.5.
Low-quality sequences were removed from the analysis, i.e., sequences with a phred Q score lower than 20, using at least 100 reads per site and applying a penalty for homopolymer regions (a quality reduction of 30% for every additional base). After editing, the regions of the RABV sequences corresponding to the Início and Final primers were removed.

The dataset used for the analysis
The dataset used consisted of 138 complete RABV genome sequences from GenBank and the RABV genome sequences obtained in this study ( Table 2). The sequences from GenBank were extracted with a PERL script using the criterion that they should have more than 10,000 nucleotides and that information about year, country and original host should be available. In all, there were 159 sequences, corresponding to 100 in the bat-related RABV cycle and 59 in the dog-related cycle. The region used for the analysis extended from nucleotide 59 to nucleotide 11801 (in relation to the PV strain, accession number M13215.1), i.e., from the first nucleotide in the N gene messenger RNA to the stop codon in the L gene.
The only sequences used were those that did not show any evidence of recombination events when examined by the RDP, GenConv, Chimaera, MaxChi, Bootscan, SiScan and 3Seq methods implemented in RDP4 (beta 4.8 version) (Martin et al., 2010) with a 95% confidence interval and Bonferroni correction for multiple comparisons.

Phylogenetic signal analysis
The phylogenetic signal content for the RABV genome and each RABV gene was investigated with the likelihood mapping algorithm (Strimmer and von Haeseler, 1997) implemented in TREE-PUZZLE v 5.2 (Schmidt et al., 2002), which quantifies the well-resolved phylogenetic quartets in the database.

Heterotachy analysis
Genome sequences were used in the heterotachy analysis with PAML v 4.6 (Yang, 2007) to test the following hypotheses: H 0 -the viral lineages in both RABV cycles are under the same molecular clock, and H A -the lineages in each cycle are under different clocks.
One hundred phylogenetic trees were constructed using the ML (maximum likelihood) method with GARLI v 2.0 (Bazinet et al., 2014), and the tree with the highest value ML was used.
To test the significance of the difference between the likelihood (L) values for each model studied, the likelihood ratio test (LRT) was used according to the equation next. LRT = 2 x (L1 -L2) where L1 is equal to the likelihood of model 1 and L2 the likelihood of model 2. In this case, the degrees of freedom employed were equal to 157, which is equivalent to the number of taxa -2 (N-2). The level of significance established was 0.05, which critical value is 206.390, in a chi-square table.
When the hypothesis of different molecular clocks for the RABV genomes in the two cycles was confirmed, analyses were carried out for each cycle and for each gene to determine the genes or regions of the RABV genome responsible for this phenomenon so that they could then be removed from the joint analyses of the cycles.
Estimates of the nucleotide substitution rate per site per year, estimates of the TMRCA and phylogeny estimates Estimates of the nucleotide substitution rates per site per year and the time to the most recent common ancestor (TMRCA) for genome sequences in which heterotachy was detected were initially made for each of the five genes for each cycle separately. After this analysis, only those genes with similar substitution rates in both cycles were used to estimate the TMRCA and carry out the phylogenetic analysis.   Substitution rates per site per year were calculated in BEAST v 1.7.4 (Drummond et al., 2012). Before the analysis in BEAST, initial substitution rates (m) for use as priors in the Bayesian analyses were estimated with Parth-O-Gen (http://tree.bio.ed.ac.uk/software/tempest/).
Maximum clade credibility (MCC) trees were inferred for each of the five RABV genes using a Monte Carlo Markov Chain (MCMC) approach implemented in BEAST v 1.7.4 (Drummond et al., 2012) and because RABV are under purifying selection (Troupin et al., 2016) we employed the SRD6 model (Shapiro et al., 2006) for coding sequences to avoid underestimation of TMRCA (Wertheim and Pond, 2011) with an uncorrelated relaxed molecular clock and a lognormal distribution. The previously estimated nucleotide substitution rate (m) was used in the model.
The MCMC converged after two independent runs with 50 million generations each. Sampling was performed for every 5000 trees, which was sufficient to obtain a sample of the stationary MCMC. This was examined in Tracer v 1.5 6 Oliveira et al.  Estimates of nucleotide substitution rates per site per year, estimates of TMRCA for the bat-related and dog-related RABV cycles and phylogeny estimates using the concatenated G-L genes MCC trees based on the concatenated G and L genes were inferred using an MCMC approach implemented in BEAST v 1.7.4 (Drummond et al., 2012) and the SRD6 model with an uncorrelated relaxed molecular clock and a lognormal distribution (Drummond et al., 2006a). The previously estimated nucleotide substitution rate (m) was used in the model. The MCMC converged after two independent runs with 50 million generations each. Sampling was performed for every 5000 trees. This was sufficient to obtain a sample of the stationary MCMC, which was examined in Tracer v 1.5 (http://tree.bio.ed.ac.uk/software/). An effective sample size (ESS) of greater than 200 was considered sufficient.

Full-genome sequencing
The 21 DNA samples from amplified RABV genomes were sequenced successfully; after assembly and editing they yielded DNA sequences corresponding to nucleotide 23 to nucleotide 11911 in the fixed PV strain (M13215.1), or the complete genomes without the primer hybridization regions. The sequences were deposited in GenBank under accession numbers KM594023 -KM594043.

Phylogenetic signal analysis
The results of the phylogenetic signal analysis showed that the RABV G gene is the most suitable for use in phylogenetic studies involving the bat-related and dog-related RABV cycles as it provided a better resolution than any of the four other genes or the complete genome (Table 3).

Heterotachy analysis
The hypothesis that the viral lineages in both RABV cycles are under the same molecular clock was less favored than the hypothesis that the lineages in each cycle are under different clocks. This finding was inferred after observing a lower likelihood value for the first model, with a significant difference between the values obtained for each model (LRT Heterotachy was observed in RABV from bat-related and dog-related cycles. Subsequent analyses were therefore carried out separately gene by gene for each cycle in order to identify the genes responsible for this phenomenon so that the phylogenetic analyses for the bat-related and dog-related RABV cycles could be carried out together using the most suitable gene(s) for this purpose.

Estimates of the nucleotide substitution rate per site per year, estimates of the TMRCA and phylogeny estimates
This analysis was first carried out separately for all the five RABV genes for each cycle. The results are shown in Table 4.
In the dog-related cycle, the highest substitution rates were found for the N gene (2.26 E-4). These were approximately twice the rates for the other genes, all of which appear to be evolving at approximately similar rates with overlapping confidence intervals: P (1.17 E-4), M (1.17 E-4), G (1.18 E-4) and L (1.34 E-4). This result suggests that the constraints on nucleotide substitutions per site in the RABV N gene in the dog-related cycle are smaller than those for the P, M, G and L genes.
In the bat-related cycle the mean substitution rates for the P and M genes (2.28 E-4 and 2.49 E-4, respectively) were similar to each other but different from those for the N, G and L genes (1.30 E-4, 1.4 E-4 and 1.15 E-4, respectively), which were also similar to each other. Unlike in the dog-related cycle, in the bat-related cycle the genes with the smallest con-Rabies virus heterotachy 7 straints on nucleotide substitutions would appear to be the P and M genes. The results for intercycle heterotachy showed that the N gene appears to be accumulating nucleotide substitutions in the dog-related cycle twice as fast as in the bat-related cycle (2.26 E-4 and 1.30 E-4, respectively) with non-overlap of confidence intervals. This finding indicates that the constraints on nucleotide substitutions per site in the N gene in RABV circulating in the bat-related cycle are greater than those in the same gene in the dog-related cycle.
The P and M genes in the bat-related cycle had higher substitution rates than in the dog-related cycle with nonoverlap of confidence intervals, indicating that the constraints on nucleotide substitutions in these genes may be greater in the dog-related cycle than in the bat-related cycle.
The L gene exhibited weaker heterotachy, and the higher rates were in the dog-related cycle (1.32 E-4 -1.36 E-4 compared with 1.13 E-4 -1.16 E-4 in the bat-related cycle).
Only the G gene did not exhibit any heterotachy between the cycles, suggesting that nucleotide substitution in both cycles may be subject to the same constraints.
To infer the divergence time between the lineages in the bat-related and dog-related cycles, the analysis was repeated but this time using both cycles together and only the G and L genes concatenated as these not only are evolving at approximately similar rates in both cycles (Table 3) but also provide the best phylogenetic signal for RABV ( Table 2).
The Bayesian inference tree built with the concatenated G and L genes (Figure 1) indicated that RABV separated into bat-related and dog-related cycles around the year 230 CE. Most of the RABV lineages found in these two cycles according to the classification proposed by Bourhy et al. (2008) and Kuzmin et al. (2012), are shown in this figure. 8 Oliveira et al.

Discussion
Although sequencing of complete genomes from different RABV isolates has been performed since 1986 (Tordo et al., 1986), to our knowledge DNA sequencing of complete RABV genomes using a single RT-PCR has not been described to date.
If only one amplification step is carried out to obtain the whole RABV genome sequence, the resulting sequence is more reliable as the use of different amplicons can, even if these are from the same viral sample, generate amplicons of different viral subpopulations for each of the segments being studied (Drummond et al., 2006b;Lauring et al., 2012).
According to the nucleotide substitution rates found for the five RABV genes in the bat-related and dog-related cycles (Table 4), we can conclude that only the G and L genes are accumulating nucleotide substitutions at similar rates per site per year in both cycles and that these genes provide the best phylogenetic signal of all five RABV genes (Table 3). Hence, phylogenetic analyses that use nucleotide substitution rates per site per year to establish the TMRCA for RABV should not use the complete genome but only the genes that appear to be evolving at similar rates in the various lineages in the dog-related and bat-related cycles and that provide the best phylogenetic signal.
The results for the TMRCA for the bat-related and dog-related RABV cycles inferred from the concatenated G and L genes indicate that the divergence occurred around 140 AD (Table 4). In the MCC phylogenetic tree the median value of the node corresponding to the divergence between the two cycles places this divergence in the year 230 AD (Figure 1). Differences of this nature are inherent to Bayesian inferences (Drummond et al., 2005(Drummond et al., , 2012. Using 151 sequences from the complete N gene (25 sequences from the bat-related cycle and 126 from the dogrelated cycle) Bourhy et al. (2008) inferred that the two cycles separated around 1250 AD. In the same study, when they used 74 sequences from the complete G gene (2 from the bat-related cycle and 72 from the dog-related cycle), authors dated this divergence to 1400 AD.
Nevertheless, our results show that although it is the most conserved of the RABV genes, the N gene neither provides the best phylogenetic signal nor has similar nucleotide substitution rates per site per year in both cycles and is therefore not the ideal gene for estimating the TMRCA. This, together with the small number of bat-related lineages and sequences studied by Bourhy et al. (2008) in the bat-related cycle, particularly in the case of the G gene, is the likely cause of the differences in the results for the TMRCA in the two cycles between this study and the study by Bourhy et al., 2008.
Based on the MCC tree built using concatenated G and L genes, the median node dates corresponding to the origin of each cycle are approximately 850 AD for the bat-related cycle (a posterior of 0.99) and 600 AD for the dog-related cycle (a posterior of 0.93). This would indicate that the representatives of RABV in the dog-related cycle are older and have been evolving for longer than those of RABV in the bat-related cycle, although the confidence intervals overlap, which is expected as the lineages in each cycle diverged at about the same time. One of the hypotheses for this difference in the data of origin of each cycle is the representativeness of the RABV sequences used in the analysis, as these may not include the oldest lineages in each cycle.
All the existing RABV lineages maintained in the dog-related cycle shared a common ancestor with the Asia 1/India domestic-dog lineage around 600 AD. Asia 1/India occupies the most basal position in this group, suggesting that this lineage, which circulates in southern India, was the first to diverge and thus the oldest circulating in the dog-related cycle. This would corroborate the findings reported by Bourhy et al. (2008). Troupin et al. (2016), using the five concatenated RABV genes and dog-related lineages, estimated this divergence to have occurred between the years 1308 and 1510 AD. This variation lies within the confidence intervals of the estimates calculated here using the concatenated bat-related and dog-related RABV G and L genes as well as each of the five RABV genes when the TMRCA using only dog-related RABV was assessed (Table 4). The discrepancy regarding the mean values is probably because Troupin et al. (2016) did not include bat-related RABV strains in their analysis and therefore did not take into account the common origin of bat-and dog-related RABV.
Similarly, Velasco-Villa et al. (2017), using the N gene, estimated that dog-related RABV diverged between 1273 and 1562 AD, a result also found in the present study using all genes separately and lineages from the dog-related cycle (Table 4), as bat-related RABV lineages were not used in the study by Velasco-Villa et al. (2017).
The Cosmopolitan complex of lineages, which includes those maintained in canids in Brazil, was estimated to have originated around 1300 AD. The basal group in this complex is a sample isolated in the 1950s in a dog in Israel (RV2324 -KF154998.1). In their study of the phylogeo-graphy of canine rabies around the world, Bourhy et al. (2008) found that the basal group for the Cosmopolitan complex was a sample isolated in a human in the 1970s in Egypt (8692EGY -U22627), which is genetically very similar to the RV2324 isolate.
The isolate from Egypt shares an ancestor with the other lineages in the Cosmopolitan complex, which includes the lineage maintained by wild canids in South America, lineages maintained by the striped skunk and wild canids in North America and some lineages maintained by wild canids in Eurasia and dogs in Africa (Kuzmin et al., 2012).
The lineages maintained in Brazil in independent epidemiologic cycles by the crab-eating fox (Cerdocyon thous) (CT-BR lineage) and domestic dog (Domestic dog-BR) apparently shared a common ancestor and have evolved independently in each of these reservoirs since approximately 1570 AD, corroborating the hypothesis that rabies in wild canids in the Americas is associated with European colonization of the continent from the 14 th century onwards (Baer, 2007).
As already shown by other authors (Bourhy et al., 2008;Kuzmin et al., 2012), in the MCC phylogenetic tree inferred from the concatenated G and L genes (Figure 1), the MexSK-1, SCSK and RAC lineages, which are maintained in the eastern spotted skunk (Spilogale putorius), the striped skunk (Mephitis mephitis) and the raccoon (Procyon lotor), respectively, share a common ancestor with all the other lineages in the bat-related cycle maintained by bats and the common marmoset (C. jacchus). The divergence between these two groups occurred around 850 AD. Although the reservoir of this ancestral RABV remains unknown, some authors have put forward the hypothesis that this original reservoir of the rabies virus in the Americas was a chiropteran (Badrane and Tordo, 2001;Hayman et al., 2016;Velasco-Villa et al., 2017).
In the bat-related cycle, the common ancestor of the complex of lineages maintained in North American carnivores and its diversification into the current lineages was dated in the present study to approximately 1170 AD, while the common ancestor of RABV maintained in chiropterans and its diversification into the many existing lineages was inferred to around 1110 AD. For both estimates, the confidence intervals overlapped.
Among the lineages in the bat-related cycle maintained by bats and C. jacchus, the most basal group consists of two Brazilian lineages, Myotis Brazil (MY-BR) and Eptesicus Brazil 1 (EP1-BR), which shared a common ancestor in around 1280 AD (Figure 1). However, this result should be interpreted with caution, as this basal position was not confirmed in ML analyses performed with the complete genome (Kotait et al., 2019), and many lineages in this cycle were not included in the present study.
The Brazilian lineage C. jacchus (CJ-BR), the only RABV lineage which has a primate as its reservoir (Kotait et al., 2019), was estimated here to have originated around 1600 AD. It shared an ancestor around 1400 AD with RABV lineages maintained by different species of bats in the genus Lasiurus (LB, LC, LS) and the species Lasionycteris noctivagans (LN) and Perimyotis subflavus (PS) in North America (Figure 1). The fact that the LC lineage also circulates in bats in the genus Lasiurus in Brazil (Oliveira et al., 2010) makes this common origin quite plausible.
Although the lineages Nyctinomops Brazil (NY-BR) and Eptesicus Brazil 2 (EP2-BR) have already been described by other authors (Kobayashi et al., 2005(Kobayashi et al., , 2007aBaer, 2007;Oliveira et al., 2010;Albas et al., 2011;Almeida et al., 2011), their probable divergence from a common ancestor around 1550 AD has not been reported in the literature to date.
The close relationship between the lineages T. brasiliensis North America (TB-NA) and D. rotundus (DR) is already known, and our findings suggest that the divergence between the two lineages occurred around 1600 AD. As suggested in other studies, the T. brasiliensis South America (TB-SA) lineage is probably the basal group for the T. brasiliensis North America and D. rotundus lineages (Kobayashi et al., 2005;Oliveira et al., 2010;Kuzmin et al., 2012), and the divergence from the common ancestor that gave rise to the TB-NA and DR lineages probably occurred around 1400 AD.
While the geographic distribution of the RABV lineages can be clearly seen in the dog-related cycle (showing the migration from India to Eurasia and Africa and finally the Americas as a result of the great maritime expeditions in the 15 th century), the geographic isolation is not so readily apparent in the lineages in the bat-related cycle. The clustering of the lineages in this cycle is primarily a result of the specific nature of the reservoirs, although certain populations in some lineages in this cycle have a geographic distribution (Streicker et al., 2010(Streicker et al., , 2012a. Nevertheless, L. cinereus has already been isolated in bats of the genus Lasiurus in North America and Brazil, showing that this lineage circulates between continents (Oliveira et al., 2010;Streicker et al., 2010;Kuzmin et al., 2012;Svitek et al., 2014).
Given the representativeness of the data used in this study, we can consider there to be five host genera (Canis, Vulpes, Nyctereutes, Mephitis and Cynictis) acting in the selective processes involved in the adaptation and selection of representative lineages in the dog-related cycle and at least thirteen (Mephitis, Procyon, Callithrix, Myotis, Eptesicus, Lasiurus, Nyctinomops, Tadarida, Perimyotis, Lasionycteris, Desmodus, Parastrellus and Antrozous) acting in the bat-related cycle. In the latter cycle, species-specific lineages can be found in certain genera (Favoretto et al., 2001;Oliveira et al., 2010;Streicker et al., 2010;Kuzmin et al., 2012).
Like the vast majority of RNA viruses, the RABV RNA-dependent RNA polymerase complex formed by the P and L proteins does not have any proofreading 3' to 5' exonuclease activity. This, together with the large number of progeny and the short interval between RABV life cycles, leads to a high mutation rate with a high probability of neutral mutations getting fixed in mutant subpopulations that maintain themselves by replicating at a lower frequency than the original (master) population in their reservoirs. When this complex composite viral population faced sudden habitat change as a result of, for example, interspecific transmission events, these heterogeneous viral populations may have been responsible for improved adaptability in a new host, leading to intraspecific transmission and the emergence of new RABV lineages or a newly dominant consensus (Streicker et al., 2010(Streicker et al., , 2012aLauring et al., 2012;Borucki et al., 2013).
Different evolutionary patterns resulting from adaptation and evolution in different reservoirs can be reflected in a variety of nucleotide substitution rates in a given site and year for the same gene in RABV distributed between the bat-related and dog-related RABV. Such evolutionary patterns can also lead to differences in the number of sites under selection in these genes as well as to differences in the selection regime acting in these (Lopez et al., 2002).
Each species of RABV reservoir can be considered a specific habitat in which each viral lineage adapts and evolves, and the variety of reservoirs currently found for RABV can be considered different adaptive landscapes (Wright 1932) in which different viral populations can adapt in different ways. Sympatric reservoirs for different species can be considered isolated environments that differ from each other in terms of the viruses for which they are reservoirs even though they occupy the same habitat (Wasik and Turner, 2013).
The separation of RABV into bat-related and dogrelated inferred to have occurred around 200 AD according to the results presented here helped shape specific genetic and phenotypic characteristics and patterns for the viruses in each of these cycles. Nevertheless, this is not taken into account in classifications and descriptions of the epidemiologic cycles of rabies, which are today divided into an urban and a sylvatic cycle (Acha and Szyfres, 2003). Instead, only factors related to reservoirs and the environment are considered, while factors inherent to the etiologic agent are overlooked.