Investigating the emergence of a zoonotic virus: phylogenetic analysis of European bat lyssavirus 1 in the UK

Abstract

European bat lyssavirus 1 (EBLV-1, Lyssavirus hamburg) is predominantly detected in serotine bats (Eptesicus serotinus) and is responsible for the majority of bat rabies cases in mainland Europe. A passive bat rabies surveillance scheme detected the virus in a serotine bat in the UK for the first time in October 2018. As of May 2024, 34 cases have been reported, 20 of which involved contact with an animal and 5 reported human contact. We investigated the emergence of EBLV-1 by undertaking comprehensive sequence analysis and Bayesian phylogenetics, based on complete virus genomes of 33 UK sequences and 108 sequences covering six countries in mainland Europe (1968–2023), including 21 French EBLV-1-positive RNA samples sequenced for this study. Sequence analysis revealed extreme similarity among UK EBLV-1 sequences (99.9%–100%), implying a single source of introduction rather than multiple independent introductions. Bayesian analysis revealed that the UK EBLV-1 sequences shared their most recent common ancestor with an EBLV-1 sequence from a serotine bat detected in Brittany, France, in 2001, with an estimated date of divergence of 1997. Within the UK sequences, the earliest divergence was estimated to occur in 2007. This study provides valuable insights into the molecular epidemiology of an emerging zoonotic pathogen and improved understanding of the risks posed to public and animal health.

Introduction

Lyssaviruses (genus Lyssavirus, family Rhabdoviridae) are a group of viruses that cause the disease rabies, a fatal zoonotic infection with nearly 100% mortality once clinical signs develop. Currently, 17 species of lyssavirus have been officially recognized (Walker et al. 2022), and bats are recognized as reservoir hosts of almost all known lyssaviruses.

European bat lyssavirus 1 (EBLV-1, Lyssavirus hamburg) is the causative agent for the majority of bat rabies cases in Europe. Previous phylogenetic analysis of the virus has revealed two distinct subtypes: EBLV-1a and EBLV-1b, with EBLV-1a having geographic distribution along an east–west axis and EBLV-1b a north–south axis, with some partial overlap in Germany, France, Poland, and the Netherlands (Davis et al. 2005, Müller et al. 2007, Smreczak et al. 2009, Mcelhinney et al. 2013, Picard-Meyer et al. 2014). EBLV-1 is primarily associated with serotine bats (Eptesicus serotinus) although detection has also been reported in the congeneric E. isabellinus (Vázquez-Moron et al. 2011). A further six lyssaviruses have been reported in bats in Western Europe: Lyssavirus helsinki (EBLV-2) in Myotis daubentonii (McElhinney et al. 2018); L. bokeloh (Bokeloh bat lyssavirus, BBLV) in M. nattererii (Picard-Meyer et al. 2013, Eggerbauer et al. 2017, Smreczak et al. 2018); Kotalahti bat lyssavirus in M. brandtii (Nokireki et al. 2018); L. caucasicus (West Caucasian bat virus, WCBV) isolated from Miniopterus schreibersii (Leopardi et al. 2021); and L. lleida (Lleida bat lyssavirus, LLEBV), also isolated from M. schreibersii, in Spain and France (Ceballos et al. 2013, Picard‐Meyer et al. 2019), and recently, Divača bat lyssavirus was discovered in a Myotis capaccinii following retrospective lyssavirus surveillance in Slovenia (Černe et al. 2023).

Spillover of bat lyssaviruses into terrestrial mammals and humans is rare but documented cases of cross-species transmission of EBLV-1 have been reported in cats (Dacheux et al. 2009), stone martens (Müller et al. 2004), sheep (Tjørnehøj et al. 2006), and three cases in humans with the most recent occurring in 2020 in France (Kuzmin et al. 2006, Selimov et al. 1989, Regnault et al. 2021). EBLV-1 has also been reported in other bat species, including Pipistrellus nathusii and P. pipistrellus (Schatz et al. 2014), and is suspected in P. pygmaeus (Folly et al. 2021). Virus detections in Pipistrellus bat species are considered likely spillover events based on geographical clustering with serotine populations (Schatz et al. 2014) rather than maintenance host, but as Pipistrellus species are abundant across Europe and the UK, this highlights the continued need to monitor for current and emerging lyssaviruses in bats to effectively assess and manage risk of transmission to the public, wildlife, and companion animals.

Since 1986, a passive surveillance programme at the Animal and Plant Health Agency (APHA, Surrey, UK) has screened bat carcasses submitted by the public for lyssaviruses. Two lyssaviruses are currently known to occur in the UK: EBLV-2, routinely detected in M. daubentonii since 1996, and EBLV-1, which was detected for the first time in E. serotinus in October 2018, followed by 33 further cases (as of May 2024). Previously, the only data supporting the earlier presence of EBLV-1 in the UK were a single seropositive serotine sampled during an active surveillance programme in Southern England in 2004 (Harris et al. 2009), but no viral sequences could be obtained from the bat. The English Channel has been previously suggested not to be a substantial barrier to gene flow in Daubenton’s bats (Atterby et al. 2010) and theorized as a transmission pathway for EBLV-2, despite no cases of EBLV-2 being reported in France (Smith et al. 2006). Similarly, for serotines it was suggested that the Channel restricted some gene flow but was not a substantial barrier with evidence suggesting frequent serotine movements between Southern England and mainland Europe (Moussy et al. 2015). Thus, EBLV-1 was predicted to eventually emerge in the UK, although the likelihood of disease establishment and spread was uncertain due to the previously reported fragmented and marginal nature of the serotine populations at the edge of their range in Southern England. Yet now regular yearly detection of EBLV-1 may indicate that the virus has become established in UK serotine bat populations. Alternatively, recent observations may be due to successive introductions of virus from continental Europe. Whether EBLV-1 is established and sustaining in the UK or being repeatedly introduced, both may be the result of a changed pattern of bat behaviour, for instance, range expansion due to climate change or land use changes (Rebelo et al. 2010, Gili et al. 2020).

In this study, we endeavour to investigate the emergence of EBLV-1 in the UK and the potential implications of recently repeated detections. We aim to determine whether detections in the UK stem from a single incursion with onward transmission, which may lead to disease establishment in UK serotine populations and associated ongoing public and animal health risks, or if it is the result of multiple independent introductions. Additionally, we aim to estimate when the virus may have diverged from its presumed mainland European ancestor, providing valuable insight into the potential timeframe of the virus circulating in the UK serotine population prior to its 2018 detection. To achieve this, we assembled a panel of virus sequences from 33 UK EBLV-1 cases (2018–23) and compared these with previously unpublished French EBLV-1 sequences, and published EBLV-1 sequences from across mainland Europe, to explore the rate of adaptation and identify the location and date of a putative most recent common ancestor.

Materials and methods

Sample collection

This study utilized a comprehensive dataset of 141 EBLV-1 whole-genome viral sequences (EBLV-1a: n = 37 and EBLV-1b: n = 104), sampled across seven European countries (the UK, France, Spain, Germany, the Netherlands, Denmark, and Russia) over a 55-year period (1968–2023). Within the dataset, 87 sequences were obtained from GenBank and 54 derived for this study through the routine bat passive surveillance scheme at the APHA, UK (n = 33) and at the ANSES Laboratory for Rabies and Wildlife, Malzéville, Nancy, France (n = 21). See Supplementary Table S1 for full details.

At the APHA (UK), 354 serotine bats (identified morphologically) were received by the APHA passive surveillance scheme between 1987 and May 2024. Of these, 34 were confirmed to have EBLV-1 infection, with the first detection occurring in October 2018. The majority of infections were in serotines found in the county of Dorset in Southern England (n = 24), with some found in Somerset (n = 9) and a recent case found in Wiltshire (n = 1) (Table 1). All cases were a mixture of adults, juveniles, males, and females. Additionally, 20 of the cases reported animal contact and five mentioned human contact.

At ANSES (France), 551 serotine bats were received by the ANSES passive surveillance scheme between 2001 and 2022, 85 of which were identified as EBLV-1 positive. A subset of 21 cases were sequenced to be used as part of this study analysis (Table 2).

Total RNA preparation

RNA was prepared from EBLV-1-positive clinical brain samples for the UK cases and ANSES-135181 using TRIzol, following the manufacturer’s instructions. RNA was prepared from the remaining French EBLV-1-positive clinical brain samples using the iPrep PureLink Virus Kit (ThermoFischer, France), following the manufacturer’s instructions. All were screened for lyssavirus using the SYBR Green reverse-transcription polymerase chain reaction (RT-PCR) assay, with additional details described earlier (Marston et al. 2019, Folly et al. 2021).

In preparation for next-generation sequencing (NGS), samples underwent one of two pre-treatment methodologies (Supplementary Table S2). Samples were either depleted of host DNA via column, using the RNeasy Mini Kit (Qiagen) following the manufacturer’s instructions, or underwent sequence-independent, single-primer amplification (SISPA) for which protocol details can be found in Supplementary Material S3.

Illumina sequencing and full-genome sequence generation

Samples were sequenced using an Illumina NextSeq platform with a Nextera XT (Illumina) kit used for library preparation. For each sample, reads were mapped to an initial EBLV-1 reference sequence (KP241939), and duplicate reads were removed to generate a final consensus sequence using a custom reference-guided alignment script (https://github.com/AMPByrne/WGS). Subsequently, reads were remapped to the closest consensus sequence derived in this study (ANSES14) to maximize sequence coverage and depth. The coverage and read depth were inspected in Tablet (v1.17.08.17), and the details are reported in Supplementary Table S2. Generated consensus sequences were aligned and compared in MegAlign 15 to identify base changes in the UK sequences. Subsequently, the genome sequences obtained in this study were aligned to previously published EBLV-1 genomes in MegAlign 15 to compare sequence distances across the dataset.

Bayesian phylogenetic analysis

The EBLV-1 dataset (n = 141) was analysed using TempEST (Rambaut et al. 2016) by performing a linear regression of root-to-tip divergences as a function of sampling time to ensure the data had sufficient temporal signal. The whole-genome sequences were then aligned in Mafft v7.445 and imported into BEAST v2.6.3 to carry out Bayesian phylogenetic analysis, allowing for the estimation of dates of divergence. The best-fitting nucleotide substitution model was determined to be a general time reversible model with gamma-distributed variation (GTR+G) using jModelTest (Darriba et al. 2012) and the Bayesian Information Criterion value.

To determine the most appropriate clock model and population prior for the dataset, models with varying clock models and population prior were investigated using BEAST v2.6.3 (Suchard et al. 2018), and the 95% highest probability density (HPD) values were compared. Strict, relaxed log-normal, and relaxed exponential clock models were compared, and coalescent constant, coalescent exponential, and Bayesian skyline population priors were compared. For each model, 1 × 10⁸ Markov chain Monte Carlo (MCMC) generations were used. Log files were analysed using Tracer v1.7.1 to ensure effective sample size (ESS) values were beyond the threshold (>200). Models where ESS <200 were discarded.

TreeAnnotator v2.6.3 was used to generate a maximum credibility tree (MCC) with a 10% burn-in, which was then visualized and annotated using FigTree v1.4.4 to include posterior probability values, time scales, and analysis of node ages. All models with appropriate ESS values also had overlapping estimated dates of divergence (95% HPD), and therefore choosing the best-fitting model was arbitrary. To analyse the evolutionary rate and in keeping with recommended parameters (Hughes 2008), a relaxed log-normal clock model with coalescent constant population prior was used for the final phylogenetic analysis.

Additionally, the dataset was aligned using MEGA v6.06, and a phylogenetic tree scaled by the number of substitutions per site was constructed via the maximum likelihood method, based on the GTR model, and using 1000 bootstrap replicates.

Analysis of viral heterogeneity

Viral heterogeneity (VH), defined as the number of genetic variations within a virus population, was compared between UK and ANSES sequences. A custom script (https://github.com/JaviNunez/ViralVariantAnalysis) was used to compare raw reads against the consensus sequence and generate the heterogeneity index (H index), as previously described (Marston et al. 2017). Here, it was used as a measure for exploring diversity within samples and comparing UK and ANSES sequences to determine if samples in either country are undergoing more adaptation (higher H index value). For consistency, when comparing the H index across samples, only those that underwent the column pretreatment were included in this analysis.

Analysis of selection pressure

The presence of selection pressure and the possibility of positive selection within the EBLV-1 genome were investigated by comparing nonsynonymous (d_N) and synonymous (d_S) substitutions per site across the virus genes for both subtypes separately and together and for the UK sequences alone. Several models were implemented using the Datamonkey server (Weaver et al. 2018, Kosakovsky Pond et al. 2020): single-likelihood ancestor counting (SLAC), mixed-effects model of evolution (MEME), fixed-effects likelihood (FEL), and fast unbiased Bayesian approximation (FUBAR) (Kosakovsky Pond and Frost 2005, Murrell et al. 2012, 2013). Only sites with a P value of <.05 for SLAC, MEME, and FEL models and with a posterior probability of >.95 for the FUBAR model were considered as showing evidence of positive selection.

The d_N/d_S ratios for samples of different sources: bat brain (n = 101), mouse brain (n = 23), and cell culture supernatant (n = 12), were also determined and compared to assess whether sample type or culture history has an effect. Samples from unknown, cat, or human sources (n = 5) were grouped together, and the d_N/d_S ratio also assessed and compared.

Evolutionary rate

The evolutionary rate of EBLV-1 viral sequences was obtained through BEAST analysis. The same BEAST model used in the phylogenetic analysis was used to determine the evolutionary rate of the entire EBLV-1 dataset (n=141). To determine the evolutionary rate of just the UK sequences and ANSES14 (n=34), first JModelTest was re-run on a separate alignment of the 34 whole genome sequences. The Hasegawa, Kishino and Yano model was determined to be the most suitable substitution model with a kappa value of 8.8587. Alongside this, parameters in BEAST were set as a log-normal clock model, coalescent constant population prior, and 1 × 10⁸ MCMC. Tracer v1.7.1 was used to confirm that ESS >200 and then obtain the ‘ucldMean’ parameter which indicates the average number of substitutions per site per year (subs/site/year). Figtree v1.4.4 was used to observe individual branch rates of substitutions.

A summary of the various tools used here to investigate the incursion of EBLV-1 into the UK can be found in Supplementary Table S4.

Results

EBLV-1 full-genome sequence generation

The EBLV-1 genome is ∼12 kb and encodes five structural proteins, which are, in order, the nucleoprotein (N), phosphoprotein (P), matrix protein (M), glycoprotein (G), and RNA polymerase gene (L). Whole-genome viral sequences were constructed for all 21 ANSES samples and 33 of the UK EBLV-1 cases (Supplementary Table S2); a genome of adequate coverage could not be constructed for Bat23-807.

EBLV-1 full-genome sequence analysis reveals the UK EBLV-1 sequences are most closely related to EBLV-1 from Brittany

Alignment of the UK whole genomes revealed 99.9%–100% identity. Alignment of the entire dataset (n = 141) then showed the UK sequences to be most closely related to five isolates from Brittany with 99.5% to >99.9% identity. For comparison, the UK sequences displayed 97.4%–99.2% identity to remaining sequences in the EBLV-1b subtype and 95.7%–96% identity to sequences in the EBLV-1a subtype.

UK EBLV-1 sequences are highly related with seven nonsynonymous changes identified across the entire genome

Across the UK sequences, there were 51 single nucleotide polymorphisms (SNPs) and one insertion observed. The breakdown of SNPs across the genome was as follows: four occurred in the N gene; four in the P gene; eight in the M gene; eight in the G gene; nine in the G–L intergenic region, where the insertion was also observed; and 18 in the L gene. Seven of the SNPs resulted in a nonsynonymous codon change, which, along with the changes in the intergenic region, are described in Table 3. The remaining SNPs were synonymous (Supplementary Table S5). Three UK sequences (Bat18-762, Bat19-177, and Bat19-949) had no SNPs. When aligned with the most closely related EBLV-1 sequence from Europe, the ANSES14 sequence had four base changes not shared with any of the UK sequences, therefore not detailed here.

Including those sequences where no SNPs were observed, the average number of observed SNPs within each year group was calculated and an increase in SNPs was found: the average number in 2018 was 0.5; the average number in 2019 was 0.67; the average number in 2020 was 2; the average number in 2021 was 3; the average number in 2022 was 4.2; and the average number in 2023 was 5.6. Bat23-736 and Bat23-1138 had the highest number of SNPs (n = 8), both sequences had one unique SNP and the remaining seven were shared with other UK EBLV-1 sequences.

Variant frequencies were also investigated for each SNP. Of the 51 observed, all except three had a frequency of >98% at their position. The three exceptions were two synonymous SNPs at position 4928 of Bat18-791 and position 10824 of Susp22-01, which had a ratio of 70:30 C:A and 80:20 T:C, respectively (Supplementary Table S5), and one nonsynonymous SNP at position 2027 of Bat23-736, which had a ratio of 90:10 C:A (Table 3).

Phylogenetic analysis and divergence dates

The data were determined to have sufficient temporal signals to proceed with the molecular clock analysis and estimate the time to most recent common ancestor (TMRCA) following a linear regression using TempEST. Bayesian phylogenetic analysis was performed on the whole-genome sequences of all 141 EBLV-1 viral sequences (Fig. 1) and an MCC tree was produced. There are two major phylogenetic clusters corresponding to EBLV-1a (n = 37) and EBLV-1b (n = 104), which can each be further divided into several clusters, consistent with previous work. To maintain consistency, the same grouping system as Troupin et al. (2017) was used, naming A1 and A2 for EBLV-1a subtypes and B1 to B7 for EBLV-1b subtypes, with the addition of splitting the B4 cluster into B4a and B4b. Case locations were also mapped where accurate geographical descriptors were known (Fig. 2). The EBLV-1a subtype consists of two main clusters: A1, containing sequences from Denmark, Germany, and Russia, and A2, which exhibits a more specific distribution with all eight sequences originating from southern France. Clusters within the EBLV-1b subtype appear to be loosely defined by specific regions although some overlap occurs, particularly in north-western France (Brittany)—B1: 18 isolates from north-eastern France, 2 from central France, and 1 from Brittany; B2: 2 isolates from northern France; B3: 5 isolates from the Netherlands; B4a: 33 isolates detected in the UK and 5 originating from north-western France (Brittany); B4b: 10 isolates from north-western and central France; B5: 3 isolates in Spain; B6: 2 isolates from central France; and B7: 17 isolates from central France and 1 isolate from Brittany. MF187812, originating from north-east France, and LT839609, from Germany, remained in their own independent lineages.

Bayesian inference revealed that the UK sequences were most closely related to ANSES14, detected in Plouguin (Brittany, France), sampled in 2001. Followed by MF187816, MF187845, ANSES-BM0211, and ANSES-135181 obtained from Plounéour-Ménez in 2000, Guingamp in 2010, Clémont in 2020, and Bringolo in 2022, respectively, all located within the Brittany region of France. The TMRCA of the B4a and B4b clusters was estimated to be 1895 (95% HPD 1868–1921). The TMRCA of the most closely related France sequence (ANSES14) and UK EBLV-1 sequences was estimated to be 1997 (95% HPD 1994–2000).

The TMRCA of only the UK sequences was estimated to be 2007 (95% HPD 2002–2012), after which the sequences divide into two main clusters, only one of which had a strong posterior probability (1.0) and an estimated divergence time of 2013 (95% HPD 2009–2017). For the majority of the UK clusters, there was no clear pattern across location, sex, or age. One exception was the cluster of cases Bat23-279, Bat23-41, Bat23-1067, and Bat21-993, which all originated from Somerset. These sequences shared their MRCA with two cases from Dorset (Bat23-870 and Bat23-997) with an estimated TMRCA of 2014 (95% HPD 2010–2018). Across the rest of the UK sequences, clusters with a strong posterior probability (1.0) were estimated to have emerged between 2017 and 2022 (95% HPD 2014–2022).

A maximum likelihood tree scaled by the number of substitutions per site is provided in Supplementary Fig. S6.

SISPA pretreatment inflates H index values

Initially, all 2022–23 UK samples were pretreated by SISPA. However, subsequent VH analysis revealed that the H index values had inflated to ∼10× the values generated for the 2018–21 cases, which were not SISPA pretreated (Supplementary Table S7). The phylogenetic and divergence dating analyses were not affected. From visualizing the genomes in Tablet (v1.17.08.17), it appeared that SISPA was increasing the number of individual variant bases across genomes, thus inflating the H index. Therefore, there was a slight effect on SNP analysis in genome regions of low coverage (<2 bp). Where the material was available, the 2022–23 samples were re-extracted and pretreated with a host depletion step via column as described earlier. NGS was then performed again. In subsequent VH analysis, the H index values had decreased to be within the range of the 2018–21 samples.

NGS could not be repeated for five UK and four France samples (Supplementary Table S2), and therefore the sequences generated following SISPA pretreatment were utilized within this study. However, due to inflated H index values, they were excluded from VH analysis, and coverage of any observed SNPs was also checked to ensure it was adequate and there was a consensus within the read depth for that SNP.

VH index reveals EBLV-1 heterogeneity is low across all sequences including the UK samples

The H index value ranged from 234 to 386, with a mean of 286, for the UK sequences and 254 to 860, with a mean of 335, for the ANSES sequences. ANSES19 had the highest H index at 860, and excluding this sequence, the range reduced to 254–378, with a mean of 303. It is currently unclear why ANSES19 has an increased H index value, and further work would be needed to investigate this in detail. Overall, the similar distribution of indices suggests that the analysed UK EBLV-1 sequences and the majority of ANSES EBLV-1 sequences have broadly similar levels of sequence variation and are under similar pressures to adapt.

Selection pressure analysis indicates that EBLV-1 is not under selective pressure

The SLAC method was used to obtain the d_N/d_S ratio for each EBLV-1 gene. The overall ratio across each gene for each EBLV-1 subtype and for the subtypes together was very low (range 0.027–0.231). Across analyses of the entire dataset and the two EBLV-1 subtypes separately, the d_N/d_S ratios followed the same ascending order across genes of N-L-G-P-M. Overall, seven sites were indicated to be under positive selection by SLAC, FEL, MEME, or FUBAR (Table 4), with two of these (site 244 in the G gene and site 1373 in the L gene) being confirmed by multiple analytical methods. All genes except the N and M gene observed at least one site being under positive selection by at least one of the models.

When the EBLV-1 subtypes were analysed separately, two sites remained identified as being under positive selection for EBLV-1a and four sites for EBLV-1b, although only one was confirmed by multiple methods. When UK samples were considered independently, only site 171 of the P gene was identified as being under positive selection by SLAC. This was identified across the entire dataset by FEL. The ascending order of the d_N/d_S ratio also differed (N-G-P-L-M); however, the index range was broadly similar, and the UK sequences do not appear to be under increased selection pressures.

The d_N/d_S ratios for samples derived from bat brain, mouse brain, supernatant, or a miscellaneous source were also very low, with a range of 0.021–0.168 (Supplementary Table S8). Across the sample types, six sites were identified as being under positive selection for bat brain (n = 3), mouse brain (n = 1), and miscellaneous source (n = 2). Two of these were confirmed by multiple analytical methods. The four sites identified within the bat brain and mouse brain samples were also identified across the dataset and some within the EBLV-1 subtypes; based on the repeated identification, it is unlikely these sites are under positive selection solely because of their sample type.

Evolutionary rate

The mean rate of evolutionary change for the whole-genome EBLV-1 dataset (n = 141) was 4.17 × 10⁻⁵ subs/site/year (95% HPD 3.47 × 10⁻⁵–4.99 × 10⁻⁵). Separate analysis of the whole-genome UK and ANSES14 sequences (n = 34) revealed a highly similar mean rate of evolutionary change of 6.92 × 10⁻⁵ subs/site/year (95% HPD 3.75 × 10⁻⁵–1.08 × 10⁻⁴).

Discussion

The comparative description of virus sequences from 33 cases of EBLV-1 in serotines in the UK, alongside the analysis of French virus sequences, permits the first epidemiological description of this incursion since its detection in 2018 (Folly et al. 2021). This allows us to investigate the introduction of the virus and the epidemiological progress of disease. While the ecology of serotine bats is not well understood in the UK, we establish a plausible description of bat behaviour, which crucially permits us to contextualize the inference drawn from the molecular investigation of these cases.

Whole-genome sequences were constructed for 33 of the 34 EBLV-1-positive cases identified in serotine bats in the UK between the first detection of the virus in 2018 and May 2024. Here we determine a striking similarity in virus sequences, suggesting that the UK sequences share an epidemiologically recent common ancestor (i.e. short transmission chain between case infections). UK whole-genome sequences exhibit 99.9%–100% identity. In comparison, the sequences had 99.5% to >99.9% identity to five closely related sequences from Brittany and 97.4%–99.2% identity to the remaining EBLV-1b subtype sequences. Based on the high similarity of the UK sequences, combined with the diversity of sequences observed in Brittany (Fig. 2), current evidence suggests that the emergence of EBLV-1 in the southern counties of the UK follows a single introduction event. An understanding of the occurrence and frequency of serotine movements across the channel would be necessary to provide a more robust conclusion.

Understanding the host ecology is critical to understanding the epidemiological process and the observed history of EBLV-1 in UK serotines. Serotines in a community mix across a network of shared roosts in their summer landscapes, including philopatric reuse of a traditional nursery site by breeding females and their pups/juveniles. This network can be widely distributed with Catto et al. (1996) identifying roosts as far apart as 10 km that appeared to be used by the same community. Using this network of roosts will also be immature and non-breeding adult females; likely to make periodic use of the nursery roost as well as other roosts known to the community (Catto et al. 1996, O’Shea et al. 2011, Chauvenet et al. 2014, Webber et al. 2016). Generally, individual roost residency times are likely to be relatively short with the community continually resorting into fission–fusion groups across its network of sites (Kerth et al. 2011, Linton and Macdonald 2019), although the use of a nursery site by reproductive females and pups is more regular. Consequently, during the active season (spring/summer), contact between community members is likely to be frequent, especially for females and young, even if the complete community rarely occupies the same site at the same time (Webber et al. 2016). Daily foraging by these bats can occur over extended areas, and combined with the distribution of the roost network, this will produce a geographical footprint within which bats might be grounded as disease develops. In England, this maximum foraging radius is reported to vary between 7.4 (Robinson and Stebbings 1997) and 11.5 km (Catto et al. 1996). Little is known of serotine hibernation behaviour in the UK, but if we assume they are consistent with the species in Europe, they are likely to use natural caves and underground features with the closest extensive karstic cave system found to the west and north of Dorset, in Somerset.

Based on the cases reported here, the progress of disease appears to conform to our understanding of the host ecology. First, disease in juveniles is evidence of transmission within summer communities, likely within a nursery setting based on the seasonally early cases in Bat18-791, Bat22-880, and Susp22-03 (Table 1), but probably true for all cases in juveniles. Furthermore, the intimate clustering of bats in their roosts provides an opportunity for an infected individual to quickly infect many of its roost mates, resulting in a close epidemiological relationship through short transmission chains or sequential infections by a common host. As long as disease progression is not rapid, this behaviour would help to explain the exceptional similarity of the UK sequences despite their geographical separation. For the majority of the UK EBLV-1 cases, the locations were intermingled; notably, the single case from Wiltshire was ∼35 km from the Dorset case it was clustered with (Bat22-881). Despite the distances between these various cases, which may be unexplained by foraging behaviour associated with a single serotine community, these appear to be infections produced in a close epidemiological association. Furthermore, only four Somerset cases were phylogenetically closely associated with each other (Bat23-279, Bat23-41, Bat23-1067, and Bat21-993), and the rest were dispersed among Dorset cases. It is plausible that this cluster of Somerset cases is an indication of EBLV-1 being present and persisting within a Somerset serotine roost. However, the dispersal of the other viral sequences from Somerset cases suggests that they may have acquired infection in Dorset, possibly in proximity to the bats they share base changes with, and bat behaviour has resulted in their geographically separated discovery.

Assuming bats retain their normal movement behaviour as rabies infection progresses, the early cluster of cases in eastern Dorset potentially represents infection within the focal summer community of serotines, with cases occurring within the putative community footprint. Subsequent cases in Dorset represent the spread to neighbouring summer colonies, with more recent autumnal recoveries from Somerset representing either further spread or alternatively a signature of a western seasonal movement of bats from colonies in Dorset to hibernacula in the west. Knowledge of the ecology of serotine bats in the UK at their very northern range edge is too uncertain to confidently corroborate this narrative. However, molecular evidence from this study suggests intimate/direct links in transmission between all cases, and our description of bat ecology explains current observations.

Direct evidence of serotine bats crossing the English Channel is yet to be shown; however, circumstantial evidence suggests that bats can cross the channel (Moussy et al. 2015) and can therefore bring across diseases. How regularly this occurs will have a direct impact on the likelihood of EBLV-1 being brought into the UK. Estimating when EBLV-1 was introduced into UK serotines is important to anchor discussions of epidemiology (establishment and spread) and for evaluating the sensitivity of the UK passive bat surveillance scheme, previously suggested to be too insensitive to provide timely detection of lyssavirus disease (Folly et al. 2021). Our molecular clock analysis estimated the closest ancestor of the UK and France EBLV-1 sequences circulated in 1997 (HPD 95% 1994–2000) and the ancestor of the UK only sequences circulated in 2007 (HPD 95% 2002–2012). It is important to note that these dates cannot differentiate the geographic locations of these ancestral bats (between France and the UK), but along with epidemiological data can help inform our interpretation.

In 2004, during a capture-sample-release survey of healthy bats, a single juvenile male bat from Sussex was reported with a detectable anti-EBLV-1 titre (Harris et al. 2009). However, no viral sequences were recovered from the bat and the roost was generally well sampled (51 bats across 3 years), leaving uncertainty about virus or disease presence at that time. A direct relationship between that bat and the cases reported here would imply continuous cryptic circulation of EBLV-1 in UK serotines in Southern England for at least 14 years (2004–18), without a single positive case being identified via the passive surveillance scheme during that time despite consistent surveillance efforts. We cannot precisely determine if the 2004 cases and current cases are linked or determine exactly when EBLV-1 arrived in the UK. The ancestral viruses of the current UK strains may have been circulating undetected in France, the UK, or another unknown location or have arrived in the UK from multiple concurrent introductions from the same source. However, the exclusive detection of UK EBLV-1 cases from 2018 despite consistent surveillance efforts since 1986, combined with the viral sequences being extremely similar and closely linked geographically and epidemiologically, provides strong evidence to suggest that the virus was introduced to the UK during a single transmission event in the recent past.

The source of the UK virus appears to be Brittany, in north-western France, clustering with a series of cases, ANSES14, MF187816, MF187845, ANSES-BM0211, and ANSES-135181 (Fig. 1). ANSES14, detected in Plouguin in 2001, was determined to be the most likely to share a common ancestor with UK sequences, proposing a transmission route of ≈309 km between Brittany and Dorset. Given that the English Channel does not appear to be a substantial barrier to gene flow (Atterby et al. 2010, Razgour et al. 2014, Moussy et al. 2015), and high gene flow among serotines in mainland Europe implies serotines can travel long distances (Bogdanowicz et al. 2013), a cross-channel route of transmission is possible, although a natural direct route between Plouguin and Dorset implies a continuous 280-km sea crossing. The shortest crossing of the English Channel is ≈44 km between Calais and Dover and would seem a less risky route for a bat to take but would include longer cross-country journeys between Dover and Dorset (≈233 km), and Calais and Brittany (≈535 km). Additionally, this route may be more directly associated with sequences within the B2 cluster, which are in northern France (Fig. 2). Alternatively, an anthropogenically assisted move is plausible with regular cross-Channel ferry routes between France (Roscoff, Saint-Malo, and Cherbourg) and England (Plymouth, Poole, and Portsmouth), and even a route between Spain (Santander and Bilbao) and England (Plymouth and Portsmouth) which closely passes Brittany. All routes implicate some cross-country journeys. Dorset, as the primary location of EBLV-1 cases within the UK, is very close to a ferry terminal at Poole, which serves ports at Saint-Malo in Brittany and Cherbourg in neighbouring Normandy which does suggest a likely pathway where the origin of EBLV-1 in the UK is Dorset and its ancestor originated from Brittany. We note that there may have been EBLV-1-positive serotines in France or the UK that were missed by passive surveillance. Indeed, a lack of EBLV-1 samples from across northern France means we cannot rule out the possibility of the source coming from another northern location.

Situating this study within the context of other EBLV-1 studies, there is a general agreement between analyses. For instance, our TMRCA estimations of the MRCA of EBLV-1 and the individual subtypes are ∼15–25 years later than previous estimates (Troupin et al. 2017) but with overlapping 95% HPD. Likewise, the clustering structure of our MCC tree (Fig. 1) was consistent with previous similar works (Troupin et al. 2017, Calvelage et al. 2021). The low d_N/d_S ratios are indicative of strong selective constraints and are consistent with Troupin et al. (2017), including the ascending order of the d_N/d_S ratios across genes being N-L-G-P-M; the EBLV-1b subtype exhibited lower d_N/d_S values than EBLV-1a; the M protein ratio halved between the two subtypes; and the N gene had the lowest ratio in all analyses, which is consistent with it being the most conserved gene in the Lyssavirus genome (Bourhy et al. 1993). The mean substitution rates generated for the whole-genome analysis are consistent with previous estimates (Troupin et al. 2017) and provide further evidence that EBLV-1 is under strong selective constraints. Finally, the H index of the UK and France EBLV-1 sequences was investigated here under the assumption that the UK serotine populations are naïve to EBLV-1 infection, in comparison to France populations, so a transitional increase in the H index may be observed, as previously shown for classical rabies virus in a host-switching event from domestic dogs to foxes (Marston et al. 2017). Analysis presented here additionally confirms that EBLV-1 is not under selective pressure to adapt to UK serotines, and the relatively high number of EBLV-1-positive serotines identified suggests that EBLV-1 is establishing within this naïve serotine population and spreading across the landscape.

None of the inference here would be possible without the bat disease passive surveillance scheme and the effort and engagement of those members of public who submit bat carcasses for testing. Especially important here are the networks of bat carers and rehabilitators, as well as private veterinary surgeons, who receive grounded bats into their care and subsequently submit many for testing, especially in the southern counties of the UK. Without their diligence and skill, this work would not be possible. Their contribution also highlights the need for the surveillance scheme as it is the principal tool to detect new diseases of bats, which may be of policy concern (e.g. public health risk). These workers, as well as health care professionals, will need the earliest and best description of a new risk to public health produced by exotic or emerging bat-borne diseases. Fortunately, between government initiatives, as well as a very active nongovernmental organization and stakeholder community, substantial guidance about the symptoms, risks, and safe working practices for members of the public and bat carers is available. While beyond the scope of the study here, evaluating the apparent sensitivity of the current scheme, using this incursion of exotic disease as a measure, appears to be a worthwhile next step to ensure that the scheme could detect a new disease of policy concern in a timely fashion.

This study identifies several gaps in information necessary to draw more confident inferences and highlights the requirement for interdisciplinary work to address them. Virus transmission between the UK and mainland Europe requires a better understanding. Factors such as climate change can change bat population structures (Rebelo et al. 2010, Piksa and Gubała 2021) and, alongside seasonality of bat activity, can impact viral evolution (Streicker et al. 2012). Therefore, reassessment of the serotine population structure is necessary to understand how well the virus may transmit within UK serotine populations, along with continued reassessment of the evolutionary rate. The use of bat radio tracking or radio telemetry, such as the Motus Wildlife Tracking System (Taylor et al. 2017), would also be beneficial towards understanding how frequently serotine and other bat species are crossing the English Channel, thus providing insight into the possibility of other lyssaviruses in mainland Europe emerging in the UK. Secondly, we suggest an active surveillance programme to specifically describe epizootiological parameters necessary for the accurate description of disease progress and risks to people, including assessing seroprevalence and rates of virus excretion.

Conclusion

Our study describes in detail an incursion of a zoonotic exotic virus to bats in the UK. We conclude that the most likely origin of the UK EBLV-1 in this study was an infected serotine from Brittany in a single introductory event. Currently, we cannot determine if the introductory event shortly preceded the first detection of disease in 2018, or if the virus was already present in the UK at a level too low to be detected by passive surveillance. Our findings highlight the potential that EBLV-1 is established in the UK serotine population and support the need for further research activity, including enhancement of the current passive surveillance scheme. Principally, this is to best provide further information to quantify the scale of the risk to people and companion animals but ultimately would help determine if EBLV-1 in the UK will spread further. The emergence of a known zoonotic virus in a new region, and the potential for changes in bat behaviour to have influenced this, has implications for disease emergence more broadly. This study provides a framework for investigating the incursion of an emerging disease within a country. The technology and analyses utilized here can be applied to investigate the emergence of future lyssaviruses and other bat-borne diseases of policy concern.

Acknowledgements

The authors would like to thank the APHA Viral Zoonoses team and the APHA central sequencing unit for technical support and advice, as well as Javier Nunez-Garcia and Alex Byrne for scripts used in analysing sequences and technical help in utilizing these. The authors would also like to thank the Bat Conservation Trust for their continued support of the surveillance scheme and bat workers/rehabilitators and members of the public who submit bats for routine testing.

Author contributions

M.E.G., G.W., D.A.M., D.L.H., and L.M.M. assisted with the conceptualization; M.E.G. assisted with the methodology; M.E.G. assisted with the formal analysis; M.E.G. assisted with the investigation; M.E.G., R.W., E.P.-M., and A.S. assisted with the resources; M.E.G. assisted with the data curation; M.E.G. assisted with the writing—original draft preparation; M.E.G., G.W., R.W., E.P.-M., A.S., J.N.A., D.A.M., D.L.H., and L.M.M. assisted with the writing—review and editing; M.E.G., G.W., and D.L.H. assisted with the visualization; M.E.G. assisted with the project administration; L.M.M. and G.W. assisted with the funding acquisition; and M.E.G., G.W., J.N.A., D.L.H., and L.M.M assisted with the supervision. All authors have read and agreed to the published version of the manuscript.

Supplementary data

Supplementary data is available at VEVOLU online.

Conflict of interest:

None declared.

Funding

This study was supported by the Department of Environment, Food and Rural Affairs (Defra), the Scottish Government, and the Welsh Government, through grant nos SV3500, SE0554, SE0570, and SE0433, as well as supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 871029. M.E.G. is supported by a PhD studentship under SE0433.

Data availability

All full-genome sequences generated as part of this study have been submitted to GenBank with accession IDs OR888933–OR888959 and PP942897–PP942924 (Supplementary Table S1).

References