One of the major questions we addressed in this study was whether bat RABV was re-introduced repeatedly into populations of mesocarnivores in the Flagstaff area during 2001–2009, or perpetuated in their populations after a single introduction in 2001. Our phylogenetic reconstructions demonstrated clearly that multiple introductions of bat viruses from the EF-W1 clade occurred into the skunk and fox populations independently. Moreover, the outbreak of 2001 most likely resulted from two virus incursions, each causing secondary transmission within the skunk population. The differences between viruses from these two clusters are supported by molecular and epidemiological data. Viruses from clade 2001a were collected during January–April from the northeastern part of the Flagstaff area, whereas viruses from clade 2001b were collected during March–July from the southern and western parts 
. An alternative hypothesis, that each rabies case in carnivores could be an individual spill-over infection from a bat, was rejected as well because the sequence similarity of viruses within each outbreak cluster was much greater than the similarity between bat viruses from the EF-W1 lineage, and between viruses from different outbreaks.
Each host shift of bat viruses into carnivores in the Flagstaff outbreaks was transient, likely due to control efforts of state and federal services, including TVR programs for skunks and ORV for gray foxes. However, it is unknown what would have happened if no such significant efforts for outbreak containment had been implemented: whether outbreaks would self-eliminate, or the virus would perpetuate in populations of mesocarnivores. Virus circulation is a balanced equilibrium of multiple components, including genetic and antigenic properties of the pathogen, pathobiology of infection on the individual host level, and ecological properties of the host on the population level 
. Although significant knowledge has been generated during a long history of studies on rabies in wildlife 
, the evolutionary genetics and host shifts of RABV and other emerging viruses are poorly understood 
Putative adaptation at the sequence level, if present, can take place at two stages of host shift. First, it may take place after cross-species transmission, for the virus to adjust to the new host environment (post-shift adaptation). Under this scenario, the virus may undergo a series of active changes in the genome under selective pressures to adapt to the new host, which can be reflected by higher dN/dS values or non-neutral convergent evolution. We performed positive selection analyses on branches that represent the post-host shift evolution. Although the test yielded positive results in M and L genes, they were most likely false positives. In fact, the EF-W1 viruses isolated from mesocarnivores during the Flagstaff outbreaks were closely related and could not be distinguished from those isolated from bats. Only a limited number of amino acid changes occurred at the post-shift stage, and these were found on terminal branches instead of on internal branches, indicating a lack of amino acid fixation after the host shift. Moreover, no convergence was detected in these amino acid changes among the three independent host shift events. Therefore, the RABVs in this study are unlikely to have undergone post-shift adaptation. As we established that each outbreak was caused by a separate virus introduction and continued only for several months, it is not surprising that adaptive changes in viral genomes did not accumulate during such a short time span.
We addressed whether viruses from the EF-W1 lineage that caused Flagstaff outbreaks were different from those circulating in bats historically. The phylogenetic analyses indicated that the historical bat isolates from the EF-W1 lineage were closely related to those causing multiple host-shift events, which suggested that viruses of this lineage have been capable of causing host shifts at least since 1975. Previous presumptive outbreaks in carnivores might have been missed because of the inherited limitations of the passive surveillance system. Alternatively, ecological factors may have favored outbreak conditions and multiple spill-over incidents during 2001, 2004–2005 and 2008–2009.
In general, a virus circulating in a certain reservoir host may already be competent for circulation in a new host (pre-shift adaptation). We examined this scenario by looking for neutral convergence between the EF-W and carnivore RABV lineages, and between the EF-W and DR+TB lineages. A number of sites were discovered to have converged between these lineages/clades. However, it is still unclear whether these convergences rendered a higher chance for successful host shift events. One amino acid substitution which can be considered potentially as a marker of pre-adaptation of bat RABV to cause effective infection in carnivores is the S→T242
in the G ectodomain. The amino acid in this position was significant for pathogenicity in the Nishigahara strain: the parental virus with A242
was pathogenic in a mouse model, whereas a mutant with S242
was attenuated 
. Another recent study demonstrated that A242
increases virus spread between infected cells 
. No studies on other amino acids in this position, including T242
, have been described to date.
The RABV variants associated with carnivores have A242
, and one lineage associated with Mexican skunks (MexSK-1) has T242
. In contrast, the majority of bat RABV variants have S242
. We know only one outbreak, caused by a virus with S242
(of a Myotis
bat origin) in gray foxes in Oregon, during 2009 
. That outbreak was self-eliminating without implementation of any control actions. It is notable that lyssaviruses of other species (which all except Mokola virus are associated with bats) also have S242
, and infect carnivores very infrequently 
. This may suggest that lyssaviruses with S242
are less pathogenic to carnivores.
By comparison, bat viruses with T242
caused at least five outbreaks in carnivores, documented during the last decade, and were transmitted efficiently among skunks, foxes and coatis. However, the significance of T242
for the ability of bat viruses to switch to carnivores should be appreciated with caveats. For example, T242
is present in viruses of the TB lineage. Only one outbreak caused by a virus from this lineage was documented in coatis in Mexico 
. However, RABV of the TB lineage is prevalent in large populations of T. brasiliensis
bats in the southern US, but no outbreaks were documented in mesocarnivores from these areas, despite their frequent chances of exposure to rabid bats. Moreover, the viruses from EF-E1 and EF-E2 lineages broadly distributed in North America have A242
, similarly to the viruses circulating in carnivores. Nevertheless, no outbreaks in carnivores caused by these viruses were documented in the US and Canada. Additional reverse genetics studies coupled with animal experiments in captivity are needed to elucidate the significance of T242
in the RABV G ectodomain for adaptation of bat viruses to circulation in carnivores.
In general, our study resolved several questions on the origin of Flagstaff rabies outbreaks, phylogenetic relationships of viruses involved in these outbreaks, and stability of their genomes over a limited period of time. Concurrently, our data raise questions on different aspects of RABV shifts from bats to carnivores. These include the significance of substitutions detected in viral genomes under different evolutionary models; likelihood of host shift to carnivores for different bat RABV variants; and outcomes of outbreaks in carnivores, caused by bat viruses, if no containment measures are implemented by humans.
In addition, we generated a significant dataset of complete and near-complete genomes of RABV from different lineages which warrant further extensive investigations. Our findings pose important questions on the molecular biology of RABV. For example, what is the significance of A→G10
substitution in the leader region, which previously had been observed only in several non-RABV bat lyssaviruses? This substitution and the following lack of complementarity may be important for virus replication, as the termini are involved in genome encapsidation and polymerase recognition 
. It is remarkable that to date this substitution has been found only in genomes of bat lyssaviruses from different species 
. On the other hand, this finding demonstrates inappropriateness of generation of “complete” viral genomes via RT-PCR with primers complementary to the genome termini, and substituting sequences of the real viral genome termini by sequences of the primers used for their amplification 
In this context, the termini of the SHBRV-18 genome should also be confirmed, as this was the only virus from the lasiurine RABV cluster that had A10
in the leader region. The SHBRV-18 was isolated from a human, presumably infected by L. noctivagans
, and was subjected to a passage in suckling mice before genome sequencing 
. Further phylogenetic refinement demonstrated that SHBRV-18 belongs to the RABV variant associated with P. subflavus
, rather than with L. noctivagans
. Potentially, alteration of the leader sequence could occur during these passages in heterologous models.
Another interesting observation was the insertion of an additional incomplete TTP after the regular TTP of the M gene in one cluster of viruses, recovered from skunks during an outbreak of 2001. Previously, similar insertions were observed in the non-coding region of the N gene, right before the regular TTP in several specimens of European bat lyssavirus, type 1 
. The significance of such insertions is unknown. Perhaps, they appear randomly in RABV genomes as a result of TTP duplication, but are eliminated during further virus passages because of redundancy. At least during the Flagstaff outbreak, such viruses were efficiently transmitted between skunks. These and other observations made in our study warrant further investigations on multiple levels, including evolutionary and functional studies.
The gene sequences, generated during this study, were deposited in GenBank, with accession numbers JQ685892-JQ686013.