Mechanisms for lyssavirus persistence in non-synanthropic bats in Europe: insights from a modeling study

Bats are natural reservoirs of the largest proportion of viral zoonoses among mammals, thus understanding the conditions for pathogen persistence in bats is essential to reduce human risk. Focusing on the European Bat Lyssavirus subtype 1 (EBLV-1), causing rabies disease, we develop a data-driven spatially explicit metapopulation model to investigate EBLV-1 persistence in Myotis myotis and Miniopterus schreibersii bat species in Catalonia. We find that persistence relies on host spatial structure through the migratory nature of M. schreibersii, on cross-species mixing with M. myotis, and on survival of infected animals followed by temporary immunity. The virus would not persist in the single colony of M. myotis. Our study provides for the first time epidemiological estimates for EBLV-1 progression in M. schreibersii. Our approach can be readily adapted to other zoonoses of public health concern where long-range migration and habitat sharing may play an important role.


EBLV-1 metapopulation model design.
We develop a multi-species metapopulation epidemic model [45][46][47] where shelters occupied by bats are represented by patches and migration events between shelters are represented by links connecting different patches. After hibernation in Avenc Davi (AD), M. schreibersii population splits between northern and southern migration routes (Fig. 1a) for mating, birthing and breeding during spring and summer seasons, before reuniting itself in Avenc Davi at the start of the fall (Fig. 1b). M. myotis bats constitute a single colony located in Can Palomeres (CP) year-round where they may get in contact with M. schreibersii throughout spring/summer months. Host mixing and possible local transmission of infection occurs within each populated patch. We propose three models for EBLV-1 infection dynamics in the bat species under study, to account for different hypotheses based on field observations and experimental knowledge. They are all based on a susceptible-exposed-   1, where no infection-induced mortality is considered and immunity wanes with rate ω. ε I is the rate of becoming infective following infection, and μ the recovery rate (b) Compartmental structure for model 2, considering lethal infection to occur with probability ρ, whereas non-lethally exposed individuals (E R ) recover with rate ε R to the permanently immune state. (c) As in (b) for model 3, where immunity wanes with rate ω. Demographic processes in the three diagrams are omitted for clarity. (d) Reproductive numbers R p 0 for M. schreibersii along each patch p of the migration path. The values correspond to the maximum likelihood estimates. The average reproductive number of the metapopulation model, ⟨ ⟩ R 0 , is also shown (black dashed curve).
Seasonality characterizes migrations flows, hosts' birth, and also transmission intensity, as the latter varies upon the degree of bats activity throughout the year. For M. schreibersii, we model it through a patch-dependent variation of the reproductive number R 0 ( Fig. 2d and Methods), a key epidemiological parameter measuring the average number of secondary cases that an infectious host can generate during the infectious period in a fully susceptible population 50 . For M. myotis, we consider a two-step function for R 0 describing hibernation in winter months (as for M. schreibersii in Avenc Davì) and breeding and mating during the rest of the year (as for M. schreibersii in Can Palomeres).
Cross-species transmission between M. schreibersii and M. myotis may occur in Can Palomeres only where the two species share the habitat. We model it with a reduction in transmissibility, 1 , to account for reduced mixing between different species (α = 0 refers to non-mixing conditions).
All parameters are described in Methods and . We compare numerical results across different models and hypotheses by introducing a metapopulation summary measure for M. schreibersii given by the average reproductive number of the metapopulation model across time and patches: EBLV-1 circulation in both species is numerically recovered only in model 1 (temporary immunity and non-lethal infection) with cross-species mixing (Fig. 3, panels a and b). The absence of mixing between M. schreibersii and M. myotis or lethal infection lead instead to negligible or null probability of persistence in one of the species   (Table S6 of the Supplementary Information, for different probabilities leading to the lethal infectious state), respectively. Also, density-dependent transmission would not allow persistence of the pathogen in any of the models explored (Table S6). Persistence probability profiles are very similar in the two host populations in model 1 with cross-species mixing. Virus circulation is maintained for all values of the immunity periods explored, with high persistence ensured at lower transmissibilities if the immunity period is short. Mixing allows persistence also for values of the average metapopulation reproductive number close to or below the critical threshold, favoring viral circulation compared to the non-mixing case (Fig. 3, panels a and c).
We performed a sensitivity analysis to allow for variations in the ecological estimates (bat population sizes, Figs S4 and S5 of the Supplementary Information; starting date of migration events, Fig. S6; duration of migration events, Fig. S7), and in the assumed length of the infectious period ( Fig. S3) yielding no variation in the predicted conditions for EBLV-1 circulation in both species.
Maximum likelihood analysis. We use a maximum likelihood approach to compare the seroprevalence data from the two species 51 with numerical results of model 1 and identify values of the reproductive number and immunity period mostly compatible with observations (see details in the Supplementary Information). Best estimate values for the unknown parameters are ω −1 = 390 (95% CI: 228-772) days and = . R 1 6 CP 0 (95% CI: 1.43-1.84). The latter is associated to a metapopulation average just above the critical threshold, = . ⟨ ⟩ R 1 02 0 (95% CI: 0.91-1.18), whereas the reproductive number is predicted to be largely subcritical in Avenc Davi during the hibernation period, R AD 0 = 0.53 (95% CI: 0.48-0.61) ( Table 1). With R CP 0 set to its maximum likelihood estimate, we find the probability of viral persistence in M. schreibersii to vary strongly, decreasing from 82% to 1% when the immunity period varies within its 95% confidence interval, with a probability equal to 38% for the best estimate ω −1 = 390 days ( Fig. S1 in the Supplementary Information). Such trend in the persistence probability is not substantially altered by variations in the assumed values for cross-species mixing (Fig. S2).

Experimental scenarios.
To assess the impact of several ecological drivers on the probability of persistence of EBLV-1 in both species, we compare our data-driven metapopulation model with a set of experimental scenarios (see Methods).
Discarding yearly seasonality of transmission and keeping a reproductive number constant in space and time leaves the persistence probability almost unaltered (non-seasonal metapopulation, Fig. 4). The essential role of migration is ensured by its northern portion, without which the likelihood of viral maintenance would be strongly reduced for a large set of values of 〈 〉 R 0 , becoming null when 〈 〉 R 0 assumes its maximum likelihood estimate (northern path only vs. southern path only). Finally, considering the breakdown of the Summer refuges patch into smaller subpopulations, i.e. through a higher spatial resolution metapopulation model, would require a slight increase in transmissibility to reach the same persistence values of the reference model (higher resolution).

Discussion
Through a spatially explicit multispecies metapopulation model based on available data from a long-term field survey on M. schreibersii and M. myotis bats in Spanish natural colonies, our study identified the main drivers for EBLV-1 persistence in the ecosystem under study and provide novel numerical evidence informing on previously unknown epidemiological, immunological and ecological factors.
Overall our findings indicate that EBLV-1 persistence relies on host spatial structure through the migration of M. schreibersii bats, on cross-species mixing with M. myotis population, and on a disease progression leading to survival of infected animals followed by temporary immunity. Even a low probability of developing lethal infection together with bats migration and reseeding through cross-species mixing would not be able to sustain the epidemic, regardless of the loss of immunity. While Lyssavirus infections are known to be generally lethal for mammals and for some bat species 52,53 , current knowledge from experimental and natural studies is not sufficient to accurately and satisfactorily define EBLV-1 disease progression in bats in natural conditions 30,33,35,40,42,43 . This is further complicated by biological and regulatory aspects. First, despite recent evidence for E. serotinus bats in France 42 , detecting a subclinical infective state in healthy bats in natural colonies is rather unlikely because of the short duration of the excretion period, and no such evidence exists yet for M. schreibersii or M. myotis. Second, the legal framework protecting European bats 54 makes field studies particularly difficult to implement, as for instance marking bats is forbidden in Europe and special authorizations need to be requested to conduct these studies. Little but precious available field data are providing however an increasing body of evidence pointing to the possibility that various bat species may experience an infective subclinical state after being exposed to EBLV-1, followed by recovery and loss of immunity, with no associated increased mortality 25,30,35,42,44 , in agreement with our model results. A constant survival rate despite recurrent EBLV-1 epidemic cycles was reported for M. myotis 44 , supporting the findings of our model. Additional field data is needed to confirm our numerical predictions on EBLV-1 infection in M. schreibersii.
Two ecological factors emerge as critical drivers for viral persistence: cross-species mixing and host migration. We found that EBLV-1 circulation would not be maintained in a single colony alone in absence of mixing with other species. Multi-species colonies are indeed a phenomenon largely observed in the field that is known to favor virus exchange 33,55 . The importance of cross-species transmission of lyssaviruses was also reported in other ecological settings, identifying phylogenetic distance as the key determinant for cross-species transmission of rabies virus in North American bat species 11,56 . No study of this kind has been performed yet on EBLV-1 hosts species, and while the mixing intensity between M. schreibersii and M. myotis in Can Palomeres is unknown, our persistence estimates remain quite robust against this assumption. Our findings suggest that an increased public ScIentIfIc REPORtS | (2019) 9:537 | DOI:10.1038/s41598-018-36485-y health attention should be focused on caves hosting multiple bat species with targeted surveillance and ecological fieldwork to improve our understanding of the role of species richness in viral exchange. We find that transmission may not depend on host density in the cave. Many bat species are indeed known to form communities (families) that are stable over short terms (daily and nocturnal activities) and long terms (between migrations) thus including members of different generations, and whose sizes are independent on the colony size 57,58 . Given that virus transmission through bites and scratches 41 requires close contact, the regularity of interactions with a limited number of hosts may explain the frequency-dependent transmission selected by our model. This result however seems to depend considerably on the species considered and the roosting behavior 59 . Additional modeling work focusing on different contexts might help better understanding such dependence across different conditions.
The other component critical for stable EBLV-1 circulation is the presence of the migratory species. Movements of infected individuals provide a mechanism for maintaining the chains of transmission through the seeding of epidemics in different patches. The importance of frequent immigration of infected hosts for persistence was also recognized in other settings 49,51,60,61 . Moreover, our findings indicate that the migratory species may contribute to pathogen persistence in species encountered along the migratory path. Given the potential for long-range seasonal movements of M. schreibersii 35,38 , this species may represent a central vector for spatial dispersion of EBLV-1 in southern Europe, where this bat species is abundant, possibly contributing to reconcile the discrepancy observed between the phylogeographical clusters of EBLV-1 variants and the geographic distribution of its common host species E. serotinus 36 .
In the ecological context under study, the northern portion of the migratory path is entirely responsible for viral persistence at the estimated immunity waning, likely because it is composed by a more complex spatial structure including a larger number of patches, thus creating more opportunities for seeding events sustaining coupled but non-synchronous patch epidemics 62 that cannot be otherwise obtained with the southern path only. Increasing spatial resolution and resolving shelters sharing similar ecological and environmental conditions (such as the ones collectively grouped in Summer refuges) would not substantially alter our predictions. These findings are important to inform future field studies minimizing data collection efforts on roosts occupation. Seasonality in mixing between hosts was instead found to have a negligible impact on EBLV-1 maintenance, suggesting that field efforts should be prioritized to provide an accurate characterization of the migration pattern, an important driver for EBLV-1 endemic circulation, instead of hosts' degree of interaction.
The maximum likelihood analysis allows us to provide for the first time previously unidentified parameters characterizing the disease dynamics of EBLV-1 in M. schreibersii. We find that transmission among M. schreibersii bats in Can Palomeres (i.e. under the highest mixing conditions) is similar to what was previously estimated for M. myotis (R 0 = 1.7) in natural colonies in Spain 44 . These conditions of transmissibility are associated to other caves being in close-to-critical or sub-critical conditions for efficient epidemic transmission. Persistence is thus mainly supported by transmission in Can Palomeres, though this cave only is not sufficient to maintain endemicity in M. myotis. The modeling predictions for the M. schreibersii immunity period are in the ballpark of previous empirical estimates of the maximum length of seropositive status observed in M. myotis (2-3 years) 35 and in E. serotinus (4 years) 42 . Individuals are predicted to be immune on average for more than a year, thus hindering virus survival because of the slow replenishment of susceptible hosts. This effect is however counterbalanced by the migration of hosts that occurs on an annual timescale, providing opportunities for seeding events in naïve populations, a mechanism already identified by theoretical works to sustain viral circulation 62 . A rather large confidence interval is obtained for the estimate of the immunity period, indicating that the available serological data are not sufficient to significantly discriminate between approximately 1 year and 2 years of duration of immunity. Also, the likelihood of persistence is predicted to vary quite considerably within this range, as the system is found in the transition between null or negligible persistence and very high probability of viral maintenance. Additional cross-sectional studies in this colony and at higher temporal resolution may help improve our estimates.
The robustness of our model to a set of ecological, epidemiological, and immunological assumptions highlights that, despite the very limited knowledge of the system, our data-driven approach is able to clearly identify the mechanisms underpinning EBLV-1 circulation in the studied populations, contributing to the interpretation of field and laboratory data as well as providing important information for prioritizing and planning further empirical investigations.
Our study provides an example of how models can yield powerful insights for the identification of the main aspects of host ecology and interaction with the pathogen driving the maintenance of the infection in the host population. Though focused on a single ecological and epidemiological context, our model may be applied to other contexts of zoonotic interest through appropriate parameterizations, in order to derive insights into the infectious disease dynamics in bat populations and improve our understanding of the risk of transmission to domestic animals and humans. For instance, through a comparative approach across bat species and lyssavirus types, it may become possible to isolate host-pathogen specific drivers (e.g. the difference in pathogenicity between bats hosting rabies virus 59 and EBLV-1) from commonalities (such as spatial structure and introductions) that would represent the key mechanisms for disease persistence across a range of lyssavirus systems. Understanding those, as well as the impact of possible anthropogenic changes on the infection dynamics in bats (e.g. a cave refuge for bats transformed into a touristic park), is essential for assessing bats capacity to serve as lyssavirus reservoirs in different geographic areas of the world and predicting the associated pathways for spillover events 9 .
Our study may be extended to other pathogens beyond lyssaviruses. While bats that host filoviruses like Ebola and Marburg viruses in Africa display very distinct ecological and behavioral traits compared to the European bat species studied here, recent surveillance work in Europe discovered a novel Ebolavirus-like filovirus in dead M. schreibersii 63 . Empirical evidence suggests that it may be pathogenic for the affected population, in contrast with the asymptomatic circulation of Ebola and Marburg viruses in fruit and insectivorous bats in Africa 64,65 . Given the ScIentIfIc REPORtS | (2019) 9:537 | DOI:10.1038/s41598-018-36485-y large geographic distribution of M. schreibersii, and the critical role of its migratory behavior for pathogen persistence examined here, the discovery of a novel and potentially deadly filovirus in Western Europe is a significant concern. Integrative models may result to be particularly important to explore such scenarios where surveillance data is still rather poor, in order to strategize future empirical data collections. With human-wildlife interactions still difficult to quantify 9 , considerable attention has shifted to animal populations. Spurred by outbreaks of high-impact viral zoonoses causing severe disease in human population, this paradigmatic shift has led to a remarkable increase in wildlife disease surveillance for early recognition of potential threats to humans (e.g. avian influenza in wild birds 66 , coronaviruses in bats 67-69 , etc.). Even after controlling for confounding factors such as investment in surveillance and research, bats were identified as carrying the highest proportion of zoonotic viruses across all mammals 6 . Moreover, hotspots of missing zoonoses were predicted for bats in Latin America and some parts of Asia, suggesting that pathogen diversity in bat species may be even larger than what currently known. Within such diversity, lyssaviruses still represent the only group of pathogens naturally circulating in bats for which a direct causal link between an infected bat and humans is well recorded 19,24,27,[70][71][72][73][74][75][76][77][78][79][80][81][82][83] . Though the largest majority of rabid cases in humans is transmitted by rabid dogs, humans can develop rabies diseases after contact with a bat infected with a lyssavirus, leading to invariably fatal cases. Moreover, accurate human morbidity estimates for infections transmitted by bats are generally hard to achieve because standard diagnostic tests cannot readily resolve the causative lyssavirus 14 . Rabies vaccination is generally recommended in specific contexts to prevent and mitigate occupational hazards, often leading to a heavy financial burden on public health infrastructure 84 .
Our study is affected by some limitations. We did not consider E. serotinus bats in our model, the host species associated with the large majority of EBLV-1 cases detected in Europe 19,20,28 . While present in the region under study, its synanthropic behavior (i.e. living close to humans populations) likely precluded interactions with M. schreibersii and M. myotis, which instead clearly exhibit a non-synanthropic behaviour, roosting in natural caves and in abandoned mines 35 . Moreover, our model considers the two host populations to be closed and isolated in the region under study. While contacts with other non-synanthropic bat species of smaller population sizes may occur in Can Palomeres, our findings show that 2 host species are enough to self-sustain the epidemic given cross-species mixing and 1-species migration, with no need for additional introductions from other sites. This may appear to be in contrast with was found in other contexts, e.g. in a rabies virus model in vampire bats where persistence was largely dependent on immigration of infected individuals 51 . Our model however already accounts for seeding events from one patch to another thanks to the migration of M. schreibersii. Similarly to the study by Blackwood et al., we find indeed that the virus would not be maintained in a single population only.
The uncovered strong spatial component to transmission dynamics may help understanding the observed spatial diffusion of the virus at a larger scale on the continent and across a diverse range of host species, through long-range migration and seeding of local populations. Given the primary role of bats in carrying the largest proportion of zoonotic pathogens across mammals, determining how hosts enable persistence and which host species are critical in a multi-host infectious disease setting is essential for disease prevention and control. Our framework can readily be extended also to other zoonotic viral pathogens of public health concern circulating in spatially dispersed bat populations.

Methods
Migration. Following hibernation in Avenc Davì (AD), M. schreibersii population splits between northern and southern migration routes (Fig. 1b). On the northern route, M. schreibersii reach Castanya (C) for mating, and then they progressively start migrating to Can Palomeres (CP) for mating and birthing. An important fraction of the population (estimated 60%) continues the migration further north to other refuges composed of breeding or summer colonies ("Summer refuges", SR, in Fig. 1a), where they stay approximately all summer. During the same period, the remaining bats stay in Can Palomeres where they share the refuge with M. myotis. Once summer is over, the entire M. schreibersii colony returns to Avenc Davì following the northern route in the opposite direction: from Summer refuges to Can Palomeres, to Castanya, to Avenc Daví for hibernation to conclude the annual migration. Bats following the southern route from Avenc Daví reach a set of caves near the coast ("Other caves", OC) and return to Avenc Daví for hibernation, reuniting with the bats following the northern route. Migration rates are set to estimates from field data (Table S1 of the Supplementary Information) 38 , and a sensitivity analysis on starting date of migration events and on their duration was performed.

EBLV-1 infection and vital dynamics.
The three models considered are all a variation of a standard SEIR compartmental model, considering transmission in each patch and possible cross-species mixing in Can Palomeres. As an example, we provide here the force of infection for M. schreibersii in Can Palomeres in model 1, whereas full details for all models are reported in the Supplementary Information: Numerical simulations and persistence analysis. We perform discrete stochastic numerical simulations of the EBLV-1 transmission in the host populations to account for the discrete nature of hosts and for stochastic extinction events that may be favored by small host population sizes. Time is considered to be discrete with a daily timescale. The epidemic is seeded with 100 infected M. schreibersii bats and 10 infected M. myotis bats in the hibernation period, and different initial conditions are explored for sensitivity. For each model and under each hypothesis considered, we ran 10 3 stochastic simulations starting from the same initial conditions and reaching the endemic equilibrium. Simulations provide at each time step the number of M. schreibersii and M. myotis in each compartment in each cave, and the number of M. schreibersii migrating from one cave to another. The metapopulation framework is implemented in C++, and technical details for simulations are reported in the Supplementary Information.

Experimental scenarios.
To evaluate the role of seasonality in transmission, we build a non-seasonal metapopulation epidemic model with the same spatial structure of the data-driven metapopulation model (i.e. patches, demographics and migration dynamics of Fig. 1), but with no variation in the transmissibility associated to the caves. The reproductive number is assumed to be constant in time and space, and equal to ⟨ ⟩ R 0 . To identify the portions of the migration path that are mostly relevant for disease persistence, we consider a northern path only metapopulation model considering only bats following the path including Avenc Davi, Castanya, Can Palomeres and Summer refuges, and a southern path only metapopulation model considering instead only bats following the path from Avenc Davi to Other caves and return. Time and rates of migration events remain as in the full migration path.
To assess the role of spatial resolution in the identification of patches and associated mixing opportunities, we consider a higher resolution metapopulation model where all caves in the Summer refuges are independently considered as patches (corresponding migration flow estimates are provided in the Supplementary Information).
Finally, we test density-dependent transmission rates for EBLV-1 dynamics in both species, alternative to the frequency-dependent assumption considered in Eq. (1), to explore different mixing conditions.