Reproduction of East-African bats may guide risk mitigation for coronavirus spillover.

BACKGROUND: Bats provide important ecosystem services; however, current evidence supports that they host several zoonotic viruses, including species of the Coronaviridae family. If bats in close interaction with humans host and shed coronaviruses with zoonotic potential, such as the Severe Acute Respiratory Syndrome virus, spillover may occur. Therefore, strategies aiming to mitigate potential spillover and disease emergence, while supporting the conservation of bats and their important ecological roles are needed. Past research suggests that coronavirus shedding in bats varies seasonally following their reproductive cycle; however, shedding dynamics have been assessed in only a few species, which does not allow for generalization of findings across bat taxa and geographic regions.
METHODS: To assess the generalizability of coronavirus shedding seasonality, we sampled hundreds of bats belonging to several species with different life history traits across East Africa at different times of the year. We assessed, via Bayesian modeling, the hypothesis that chiropterans, across species and spatial domains, experience seasonal trends in coronavirus shedding as a function of the reproductive cycle.
RESULTS: We found that, beyond spatial, taxonomic, and life history differences, coronavirus shedding is more expected when pups are becoming independent from the dam and that juvenile bats are prone to shed these viruses.
CONCLUSIONS: These findings could guide policy aimed at the prevention of spillover in limited-resource settings, where longitudinal surveillance is not feasible, by identifying high-risk periods for coronavirus shedding. In these periods, contact with bats should be avoided (for example, by impeding or forbidding people access to caves). Our proposed strategy provides an alternative to culling - an ethically questionable practice that may result in higher pathogen levels - and supports the conservation of bats and the delivery of their key ecosystem services.

Background

The order Chiroptera is the second largest order of mammals with more than 1000 identified species [1]. The members of this order, bats, provide important ecosystem services (reviewed in [2, 3]). For example, insectivorous bats can reduce arthropod herbivory [4-6], increase agricultural yields [7], reduce the need for insecticides [8], and prevent large financial losses in agriculture [9-11]. Plant-visiting chiropterans provide pollination and seed-dispersing services (reviewed in [3]), certain nectivorous bats are pollinators of economically important plants [12], and frugivorous bats can be important for reforestation [13]. Finally, cave-roosting bats produce guano, the main energy source in many cave ecosystems [14, 15], and the mining of this product is an income source in poor communities [16].

However, current evidence supports that bats are a natural host of several disease-causing viruses across the globe, including zoonotic viruses, such as rabies virus (Rhabdoviridae, genus Lyssavirus); Hendra and Nipah viruses (Paramyxoviridae, genus Henipavirus); and Marburg and ebolaviruses (Filoviridae, genus Marburgvirus and Ebolavirus, respectively; [17, 18]). Bats are also hosts of several viruses of the Coronaviridae family [19-22]. Molecular evidence suggests that the Severe Acute Respiratory Syndrome Betacoronavirus (SARS-CoV betaCoV) and the Middle-East Respiratory Syndrome betaCoV (MERS-CoV) originated from bats [23-31]. Both viruses emerged in humans during the past two decades, specifically in China (2002) and Saudi Arabia (2012). The SARS-CoV pandemic included 8096 cases in 27 countries with a ~ 10% case fatality [32], while MERS-CoV has affected 2279 people in 27 countries with a case fatality of ~ 35% [33]. Incidental cases of MERS-CoV are still detected mainly in Saudi Arabia [32], and it is thought that camels are important for human infection [34-36].

Today it is known that: i) a high genetic diversity of coronaviruses (CoVs) is present in more than 100 bat species, including viruses related to SARS and MERS CoV [37]; ii) CoVs are prone to move and adapt to new host species [38]; iii) plausibly, all mammalian-adapted CoVs may have originated in bats [20, 21, 38], including a recently emerged highly fatal Alphacoronavirus in piglets [39] and the 229E human CoV [22, 40, 41]; and iv) CoVs found in bats can use human receptors for cell entry [21, 25, 42, 43]. These lines of evidence suggest that future spillover of coronaviruses humans is feasible.

Because CoVs are found in bat species that have adapted to be in close contact with humans, such as the straw-colored fruit bat (Eidolon helvum) and the Brazilian free-tailed bat (Tadarida brasiliensis [44, 45]), high contact “bat-human” interfaces currently exist around the world. If the bats in these interfaces shed CoVs with the ability to infect humans, then opportunities for spillover through direct exposure to feces [37] or the contamination of food are created, as these viruses can remain infectious in the environment for days [46]. Therefore, strategies aiming to mitigate human exposure to CoVs, and thus, the risk of spillover and disease emergence are needed, while supporting the conservation of bats and their important ecological roles.

Longitudinal sampling with specific species has shown that the proportion of bats shedding CoVs varies seasonally [47-50] and that fecal CoV-RNA loads can also be heterogeneous over time [51, 52]. If exposure through contact with bat feces is a main pathway for zoonotic CoV spillover to humans but shedding of these pathogens is not uniform over time, then mitigation strategies aiming to prevent bat CoV-shed exposure could be targeted temporally, directed especially at high-risk seasons. Such a strategy could guide policy in limited-resource settings where sampling bats for CoV testing is not feasible and it could support an ethically acceptable management to mitigate spillover risk. However, the few species and locations tested to date do not allow for identification of a potential seasonal shedding pattern to responsibly suggest temporal spillover risk management across species and geography. Therefore, assessment of the CoV dynamics in a broader range of bat species that show different life history traits, as well as in diverse geographic and ecological circumstances, could be extremely useful.

To this end, we evaluated the dynamics of CoV shedding in different bat species sampled in several locations in East Africa at different times of the year. This geographical region has been identified as a hotspot of pathogen emergence [53], where CoV host switching events seem to be higher compared to other areas [22], but, to our knowledge, no study on CoV dynamics in bats has been conducted. Specifically, we hypothesized that bat species exhibit seasonal trends in CoV shedding that are associated with the reproductive season. We assessed this hypothesis by fitting Bayesian statistical multivariable models to evaluate whether CoV shedding in bats is positively associated with the time period when pups are becoming independent from the dam. Beyond the inclusion of several species sampled in different countries at different times, we explicitly identified the reproductive events for each species at the time of sampling and also included other traits, such as the aggregation of individuals at the roost, that may be involved in CoV dynamics.

Methods

Sample collection

Samples (rectal swabs and fresh feces) were collected from bats captured in Uganda, Rwanda, and Tanzania (Fig. 1), between September 2011 and April 2014 with permission from local authorities and under the Institutional Animal Care and Use Committee at the University of California, Davis (protocol number: 16048).

Fig. 1

Unique locations where samples from bats were obtained

Bats were captured in 36 unique locations between latitudes − 0.9739 and − 10.7506 (Fig. 1). These locations were selected because they represented potentially high-risk interfaces for contact between bats and humans, such as areas of land-use change, human-dwellings, ecotourism sites, markets, and places with potential for occupational exposure [54]. Locations in close proximity (Euclidean distance < 20 km) in which sampling was conducted within the same week were considered a single sampling event. The remaining sampling events that occurred in the same location or spatially close to others but conducted in different weeks were considered independent sampling events. As result, we collected samples from 30 unique sampling events.

All captures were conducted using mist nets set in the early morning or at dawn. Individuals were released after sample collection. Samples were handled as previously described [22, 55]. Each sample was immediately transferred to vials containing Viral Transport Media and NucliSens® Lysis Buffer (bioMérieux, Inc., Marcy-I’Étoile, France), which were maintained in liquid nitrogen until storage in a -80C freezer in each country.

Coronavirus detection

RNA was extracted from all samples, and cDNA was prepared as previously described [22, 55]. Two broadly reactive consensus PCR assays targeting non-overlapping fragments of the orf1ab were used to detect both known and novel CoVs [56, 57]. Amplified products of the expected size were cloned and sequenced as described in [22]. Sequences were analyzed and edited using Geneious Prime 2019.1.3 [58].

A sample was considered positive when at least one PCR assay yielded a sequence that matched coronaviruses in GenBank. Coronavirus sequences were classified as belonging to a specific taxa following previously described methodologies [22].

Bats age, sex, and species identification

Bats were categorized as adults and juveniles based on size, and morphological and behavioral characteristics were observed at sample collection. The sex of the bats was also recorded.

Identification of some bat species can be challenging in the field. For this reason, field-identified species were confirmed by DNA barcoding using the cytochrome b or cytochrome oxidase subunit 1 mitochondrial genes [59]. Obtained sequences were compared against sequences in the GenBank and Barcode of Life databases [60, 61]. When possible or necessary, sequences from both genes were used for species identification. A threshold of 97% nucleotide identity was used to confirm the species; sequences with 95–97% nucleotide identity were assigned a “conferre” (cf.) species status, and sequences below 95% nucleotide identity were either classified to the genus level or as unidentified. Sequences with > 97% nucleotide identity to more than one species for either gene, were classified to the genus level unless they clearly clustered with sequences from other animals in the same geographic area.

If barcoding results for all of the first ten bats tested per species were in agreement with the field identification, we assumed that the field identification for the remaining bats of that species in each country was correct. Otherwise, all of the remaining samples were barcoded to ensure correct speciation.

Bat life history

We recorded when sampled females were pregnant based on abdominal palpation, had attached pups (indicating recent parturition), and were lactating, as well as when juveniles were captured. Therefore, we were able to track pregnancy, lactation, and recent birth pulses. Moreover, we accessed the data in the PanTHERIA [62] and Amniote [63] databases, and we thoroughly reviewed the literature on the biology of the bats species we sampled for latitudes similar to our sampling locations. With the gathering of these information sources we established the timing of the birth pulses, lactation periods, and the weaning of pups for each species. For the details justifying the dates inferred for these three life history events for each species and the corresponding bibliographical references see Additional file 1.

Once the timing of these events was confirmed or inferred, we were able to establish 2 seasons that occur at least once during the year across all observed bat species: i) when juveniles are being weaned and female-pup contact decreases (“Recent weaning” [RW]) and ii) the rest of the year (hereafter “N-RW” for “Not recent weaning”). We chose to evaluate risk of CoV shedding for the first period because past longitudinal studies with microchiropterans in Germany and China found higher CoV-RNA loads approximately 1 month after parturition [51, 52]. Similarly, peaks were found 2 months after the formation of a maternity colony of Myotis macropus [49], which would match a post weaning period for that species. Here, we defined the end of the RW period as 1 month after the last pups were weaned. We assumed that 1 month would provide a reasonable time window for the colony to “clear” the CoV susceptibility status of this period and acquire the susceptibility corresponding to the season(s) when weaning does not occur (N-RW period), if differences actually exist.

Finally, once we determined these two seasons, we categorized each bat sample into one of them depending on the week of the year in which they were taken. Because some species had more than one litter per year, there could be more than one RW period during the year. It is worth noting that we were able to define these periods for those species in our dataset that have synchronized reproduction, whose biology was properly described, and whose taxonomy is generally accepted. When we could not assign a reproductive period to specific bats, this season was imputed (see Methods: Statistical Analysis).

Bat species traits

We characterized specific traits of each bat species studied based on previous scientific literature on pathogen dynamics in bats [51, 64–71]: colony size (small, medium, or large if a typical colony contains one to dozens, hundreds to thousands, or thousands to millions of individuals, respectively); roost type (“closed” if the species has been reported to use caves, mines, roofs, or other confined spaces; or “open” if the bats have been typically reported roosting in the foliage of trees); the aggregation of bats in clusters while roosting (no, yes); and the number of litters per year of the species at equatorial latitudes. References are provided in Additional file 2. We also considered data from PanTHERIA and Amniote [62, 63].

We could not include other species traits, such as the mating strategy (harem or another) and the segregation of females in maternity colonies, because available studies were incomplete or contradictory. We did not include factors, such as multi-species occupancy of the roost, because we did not observe all of the roosts. Further, we did not assess postpartum estrus, as within the study area, it is only known to occur in some Molossid bats [72], of which we only sampled a small number.

Statistical analysis

To statistically estimate the association between the RW and CoV detection we used a Bayesian inference to model the detection of CoVs as a Bernoulli process of the form: where CoV, the detection of CoV in rectal swabs (1 if detected, 0 otherwise) from the ith bat (sample), is assumed to follow a Bernoulli process parameterized by p, the probability of CoV detection on the ith bat. This parameter was related to a set of candidate covariates as:

with S~Normal (0, σ) and Sp~Normal (0, σ).

Here X and Y are binary covariates representing the RW season and juvenile age category. We specifically included these two terms to separate the potential association of the season with CoV detection from the seasonal presence of juveniles. Because it was not feasible to allocate all species in the RW or N-RW seasons based on previous research, we assumed that these unknown reproductive seasons were “Missing at Random” and they were imputed as a function of the latitude at sampling, the day of the year of the sampling event, the number of litters per year of the corresponding species (one litter per year versus more than one litter per year), and the historical precipitation of the month at the sampling event location. The description of the imputation model is provided in Additional file 3.

The terms S and Sp represent the sampling event- and the species-specific intercepts, respectively, because we assumed that bats sampled in the same event and bats belonging to the same species were not independent with respect to CoV detection. No sampling event involved the same bat colony in successive RW and N-RW seasons, therefore, we assumed that CoV detection was not temporally correlated within sampling locations.

We constructed the model by adding other covariates one-at-a-time to this working model: the C1 categorical variables; and they remained in the model if they were judged to confound the relationship between CoV detection and the reproductive seasons or between age and CoV detection (i.e., their inclusion caused meaningful changes in the Posterior Probability Distributions [PPDs] of the specific reproductive season or age coefficients). Finally, C1 categorical variables could be retained as well if they were marginally associated with CoV detection (the corresponding coefficient PPD did not include zero). To assess the goodness-of-fit of the models, we evaluated the congruence between CoV detection in the data and in the posterior predictive distributions yielded by the models by: i) ages and seasons, ii) age, iii) per season, iv) per age and season, and v) per sampling event.

All models were constructed using “Stan” v. 2.17.0 [73] which was run from “R” v. 3.6.0 [74] through the package “RStan” v. 2.17.3 [75]. Weakly informative priors were assigned for all coefficients: Normal(0, 1.5) for the estimates of α0, β’s, γ, and ρ’s. The σS and σSp had a prior Half-Cauchy (0,5) following previous suggestions [76, 77]. The PPDs were estimated by sampling in parallel from 4 MCMC chains for 4000 iterations following 3000 iteration warm-up for a total of 4000 saved samples for each parameter PPD. Convergence was assessed by the Gelman-Rubin statistic [78] and graphically using trace plots. The code to fit the models is available at https://github.com/dmontecino/East-African-bats-and-CoV-shedding.

Results

Bat samples, age, sex and reproductive seasons

We sampled 753 individuals, all of them aged and successfully identified to belonging to 15 species. Nycteris thebaica, Pipistrellus hesperidus, and Rhinolophus clivosus were assigned the “conferre” status. Hipposiderids were assigned only to the genus level because barcoding did not provide certainty on species identification in line with previously recognized taxonomic difficulties [79-81]; however, the biology of the candidate hipposiderid species is similar (Fig. 2; [82, 83]), and we were able to used them for estimation purposes. We excluded Scotophilus viridis (n = 6) from the analysis because their reproductive traits are unknown, and this species also has taxonomic difficulties for species identification [84]. Therefore, these six individuals were removed, and 747 bats were included in the study.

Fig. 2

Summary of the inferred reproductive periods of the bat species sampled. The red squares show the week of a typical year each bat species was sampled. The yellow, green, and blue polygons show the assigned extent of the birth pulse (yellow), lactation (green), and mating-pregnancy periods (blue). The grey polygons represent the inferred “Recent weaning” period(s) per bat species. In species with unknown lactation length (N. nana and T. persicus), the light green polygons represent the likely overextension of this period not including any other bat in the “Recent weaning” period. The question marks show the period we could not infer the corresponding reproductive activities

We were able to infer the reproductive season for all bats except for Lissonycteris angolensis, Rhinolophus cf. clivosus, and Pipistrellus cf. hesperidus (n = 117, 104 adults and 13 juveniles) that had limited available knowledge on biology and reproductive season. These reproductive seasons were imputed as explained above. This imputation process did not substantively affect the proportion of bats in each reproductive period or the crude CoV detection per age (Additional file 3: Figure S3.2). Moreover, we had partial data for Neoromicia nana and Triaenops persicus lactation period and we assigned one that is likely longer than it would be in reality based on the other species. Even under this overextended period, none of these bats (N. nana and T. persicus) could have being sampled during the RW season, so this knowledge gap did not risk a misclassification (Fig. 2).

In the end, 274 and 356 bats were allocated in the N-RW and RW periods, respectively. The distribution of bats across the reproductive periods per age and sampling event was heterogeneous as expected due to the opportunistically nature of sampling. Indeed, 233 and 273 adult bats were in the N-RW and RW periods, respectively, while 41 and 83 juvenile bats were sampled in these seasons, respectively.

Coronavirus detection and identification

In total, 30.79% bats were positive for CoVs (n = 230). Within the subset of adult bats, 26.22% were positive (n = 160) while 51.09% of the juvenile bats were positive (n = 70). The detection of CoV shedding was variable across seasons and bat species, as well as across sampling events (0, 100, 9.69, and 25.84% for the minimum, maximum, median, and mean detection, respectively). The CoVs found per species are shown in Additional file 4: Table S4.

Species traits

A summary of the roosting and reproductive traits of the bat species sampled is provided in Table 1. All bats except N. nana (n = 9) and E. helvum (n = 315) roosted in “closed” structures, such as caves, abandoned mines, and roofs. Within the group of bats using “open” structures, E. helvum was the only species with CoV positive individuals. Therefore, we did not use this covariate to assess a potential association with CoV shedding.

Table 1

Summary of traits by sampled chiropteran species

Chiropteran species	Trait
	Colony size	Roost type	Aggregation in clusters while roosting	Number of litters per year
Chaerephon pumilus	Medium	Closed	No	3b
Coleura afra	Large	Closed	No	2
Eidolon helvum	Large	Open	Yes	1
Hipposideros sp.	Medium	Closed	Yes	1
Lissonycteris angolensis	Medium	Closed	No	2
Mop condylurus	Medium	Closed	Yes	2
Neoromicia nana	Small	Opena	No	1
Nycteris cf. thebaica	Medium	Closed	No	1
Pipistrellus cf. hesperidus	Small	Closed	Yes	1
Rhinolophus cf. clivosus	Small	Closed	No	1
Rousettus aegyptiacus	Large	Closed	Yes	2
Taphozous mauritanus	Small	Closed	No	2
Triaenops persicus	Large	Closed	Yes	1

a Neoromicia nanus roosts in folded banana leaves. Could be considered “Closed” also

b It has been proposed up to 5

The models showed adequate sampling. The 4 Markov Chain Monte Carlo chains converged graphically, whilst all Gelman-Rubin statistics were < 1.004. The selected model had a number of effective samples for each coefficient of at least 1636. The data were properly fitted, as well (Additional file 5: Figure S5.1), although some predictions lacked precision. The PPDs of the fixed coefficients are shown in Additional file 5: Figure S5.2.

The selected model to assess periodic differences in CoV shedding included season and age, species-specific intercepts, and sampling event-specific intercepts. Beyond the species-specific terms, we included a binary categorical covariate equal to 1 for E. helvum and T. persicus and 0 otherwise. We incorporated this term because CoV detection in these species was remarkably higher than the other species. As expected, this fixed effect was correlated with the corresponding species-specific intercepts (the remaining correlations were all low); however, we decided to keep it to assess if the main findings hold even when accounting for the bat species with highest detection. The corresponding means, standard deviations, and 90% HPDI are shown in Table 2.

Table 2

Summary of the posterior probability distributions of the fixed-effects coefficients of the selected model

Covariate	Mean	SD	90% HPDI
Intercept (α0)	−3.21	0.77	−4.375 - -1.944
Recent weaning	1.62	0.68	0.538–2.772
Juvenile age class	0.66	0.26	0.233–1.078
E. helvum – T. persicus	1.33	1.17	−0.511-3.282

SD Standard deviation and 80%, HPDI = 90% high posterior density interval

The coefficients’ PPDs from the selected model indicate an association between age and CoV shedding, with juveniles presenting 1.26–2.94 times higher odds to shed compared to adult bats (90% HPDI). The coefficients’ PPDs also point to an association between the reproductive season and CoV shedding as well, with an estimated odds 1.71–16.00 times higher to shed during the period when pups are being compared to other seasons (90% HPDI). The proportions of CoV shedders estimated by reversing the 90% HPDI of the logits were: 0.02–0.22, 0.09–0.59, 0.01–0.13, and 0.05–0.42 for juveniles during the “N-RW” and “RW” periods, and adults during the “N-RW” and “RW” periods, respectively (90% HPDI). These values refer to bats not belonging to the species E. helvum or T. persicus. Finally, the predicted CoV detections, based on the posterior predictive distributions, were 0.01–0.18, 0.06–0.54, 0.00–0.05, and 0.04–0.36, for these same groups (90% HPDI; Fig. 3, left). In practical terms, these last values imply that juveniles during the “RW” period are, on average, 3.34 times more likely to be detected shedding CoVs compared to juveniles in the “N-RW” period. Adults during the “RW” period are, on average, 3.93 times more likely to be detected shedding CoVs compared to adults in the “N-RW” period. In both seasons, juveniles are, on average, more likely to shed CoVs, than adults.

Fig. 3

a Modeled proportion of chiropterans shedding coronavirus viral particles in two reproductive periods inferred for these bats by age class. b The estimated species-specific intercept coefficients. The lighter colors represent the range where the 90% of the estimated detection values are concentrated (the High Posterior Density Interval). The darker colors show the 51% High Posterior Density Interval

The selected model suggests a higher odds of CoV detection in E. helvum and T. persicus compared to other species. The species-specific intercept terms, once the E. helvum - T. persicus effect is included, suggest no further differences in terms of CoV shedding (Fig. 3, right); however, the estimates are not precise. The sampling event-specific random intercepts suggest that a few specific locations could show differential CoV shedding but that most of them do not explain further variation (Additional file 6).

Discussion

If coronavirus shedding by bats follows temporal patterns that are generalizable across species and locations, then mitigation strategies targeting the prevention of human exposure and potential spillover could be directed toward high-risk periods, through mechanisms that can also support bat-human coexistence and the provision of bat ecosystem services. Previous research has focused on viral identification in specific locations and in few species [47-50], resulting in a limited representation of viral dynamics in association with few ecological settings, biological traits, and reproductive strategies. Additionally, few studies of coronavirus shedding patterns have employed statistical models and, in consequence, the potential complex web of factors and causal relationships that may determine this process has not been fully explored. Here, we aimed to address these issues by statistically modeling coronavirus RNA detection in several bat species, captured at different times and locations in East Africa and involving different ecological contexts and life histories. Using data from several hundreds of bats, we found that, beyond spatial, taxonomic, and life history differences; i) the odds of coronavirus shedding is higher during the period when pups are being weaned (up to a month after the lactation period is over), and ii) juvenile bats have higher odds to shed these viruses. Moreover, the ratios of predicted detections per bat category (age and reproductive season) suggest that juvenile bats during the recent weaning period have relatively higher shedding compared to bats out of this period no matter their age. Caution must be taken with these ratios because we used a logit link and our data had high proportion of CoV shedding in specific groups. However, our results are consistent, and they are in agreement with previous research conducted in a restricted number of species and locations.

Similar seasonality of coronavirus shedding has been observed in Germany, Australia, Thailand, China, and Ghana (West Africa). In the specific species involved in these previous studies, higher coronavirus shedding and viral loads were detected weeks after the birth pulse [47–52, 85]. Further, and consistent with our results, detection of higher levels of coronavirus in juveniles has been reported in micro- as well megabats from Africa, Asia, Europe, and North and South America [22, 48, 50, 86–89].

It has been proposed that the increased detection of coronaviruses after the birth pulse is attributable to the waning of passively-received maternal antibodies in juveniles [51]. This idea has been frequently cited; however, we are not aware of any longitudinal age-specific coronavirus seroprevalence study in bats. Such studies are important to understand the drivers of pathogen persistence and spillover risk, and in consequence, to ethically manage and prevent bat pathogen exposure. Nevertheless, this kind of research is difficult to conduct due to logistical challenges, our questionable ability to obtain statistically-representative samples across age groups, cross-reactivity of serological assays, and the difficulties to differentiate serodynamics derived from closed-population processes from those caused by migratory movements. Although extrapolations for antibody dynamics across viruses and species are not simple [90], bat serodynamics for Hendra virus are congruent with the increased detection of coronaviruses after the birth pulse. Pups passively receive maternal Hendra virus antibodies which decline after the first month of age up to 6 months of age [91-94]. This decline would lead to a period in which young bats tend to be more susceptible to infection, become infected, and then shed virus. Consistently, coronavirus shedding peaked weeks after the birth pulse in a German and a Chinese species [51, 52] and immunologically naïve bats shed higher coronavirus loads [30]. Over time, as young bats clear Hendra virus infection, they become seropositive again [91-94]. Concordantly, capture-mark-recapture studies support the clearance of coronaviruses in infected bats [47, 95], which would become seropositive. However, young bats may not reach adult seroprevalence levels until they are older than a year, as occurs with Hendra virus [91, 92, 94]. Therefore, the population of juvenile bats would remain comparatively more susceptible to viral infection and shedding beyond the period immediately after weaning. Age-specific henipavirus seroprevalence in African E. helvum is in agreement with the serodynamics described for Hendra virus [90].

Additionally, coronavirus transmission may be favored by high colony density created by the birth pulse, as previously proposed [85, 95], and then the seasonal influx of susceptible juveniles could accelerate viral spread across the entire colony, including adult bats. Indeed, adult Myotis macropus in an Australian colony showed a peak of coronavirus detection after the birth pulse [85]. The peak of coronavirus detection for two Hipposideros species and Nycteris cf. gambiensis sampled in Ghana occurred during the months that encompassed the birth pulse and nursing after accounting for the age of the sampled individuals [50]. Higher coronavirus infection has also been reported in lactating females [85, 88], which overlap with the period of pup weaning and decay of maternally-derived immunity; however, the opposite has also been found [48, 50, 85].

In practical terms, public health managers could anticipate high risk periods for coronavirus shedding to target interventions. Assuming that higher spillover risk is a function of higher viral shedding [67] and that all coronaviruses with zoonotic potential behave ecologically similarly to coronaviruses detected in this study, managers could target the prevention of human-bat direct (consumption) or indirect (bat droppings) contact specifically during the high-risk season: around and just after weaning, the timing of observable juveniles or individuals smaller than adults. For the species and interfaces defined herein, those management periods have now been determined (Fig. 2). For others, direct observation of bats at high-risk transmission interfaces could be used to identify time periods when non-adult sized bats are present. However, observation of dependent pups is not always easy [96]. Of course, for specific species, birth pulses and lactation seasons could also be used to more precisely establish high risk periods similar to the methods we used here, including a combination of direct observation, reports from previous literature, and consultation with knowledgeable bat biologists. Our proposed risk-driven strategy i) is evidence-based, as it builds upon coronavirus shedding patterns observed across several chiropteran species present around the world; ii) does not require the advanced laboratory capacity often lacking in resource-restricted settings where intense bat-human interfaces usually occur; iii) is a good alternative to the ideal but expensive and resource-intensive longitudinal surveys; and iv) it may prevent the exposure to viruses belonging to other taxa whose observed bat shedding dynamics resemble our findings for coronaviruses (e.g. paramyxoviruses [97]),

The lower coronavirus detections in African emballonurids (C. afra and T. mauritianus) and the higher detections in E. helvum, African hipposiderids (Hipposideros sp. and T. persicus), and R. cf. clivosus that we found are consistent with previous reports [22, 44, 50, 98] and should be considered by managers when providing risk-based spillover prevention strategies. Moreover, SARS-like coronaviruses in Africa have been found in hipposiderid, rhinolophid, and molossid bats [22, 57, 98], and MERS-related coronaviruses have been found in vespertillionid bats [22, 28, 29]. Therefore, it seems reasonable to prioritize the identification of birth pulses and lactation seasons, and thus determine high-risk periods of coronavirus shedding, for these bat families. Interestingly, E. helvum roost in tight clusters that can contain hundreds of individuals [99], similar to T. persicus. On the other hand, emballonurids, showing the lowest crude detection levels across families (represented by Coleura afra and Taphozous mauritanus here) tend not to cluster while roosting [100, 101]. We did find an association between coronavirus shedding and whether the species typically aggregate in clusters while roosting when the variable “E. helvum – T. persicus” was not considered, but we chose a different model not including this term because we did not directly assess bat roosts and our categorization may oversimplify the continuum from mostly solitary roosting (e.g., Neromocia nana) to common tight aggregations of bats (e.g., Mops condylurus). Using this categorization could be misleading, as some species differentially cluster while roosting depending on temperature, colony size, colony type (e.g., maternity colony versus not a maternity colony), and season [72]. Future studies should consider the roosting habits of bats, as this trait could further support risk-based management to prevent or reduce human exposure.

The risk-driven strategy we propose provides a contact-reduction alternative that is ethically favorable compared to often-employed measures, such as culling or other reactive measures, that ensue when the public becomes aware of a health threat without a suggested practical option to reduce their risk for exposure. In addition to ethical concerns and being logistically difficult and expensive, culling has failed to reduce disease in wild populations and can result in even higher pathogen levels. For example, “badger culling can make no meaningful contribution to cattle tuberculosis control in Britain” [102]. A culling program to reduce Echinococcus multilocularis prevalence in red foxes (Vulpes vulpes) resulted in an increase of infection [103]. This strategy has also failed to control rabies in canids around the globe [104]. Similar results have been observed in bats. Culling failed to reduce rabies seroprevalence in Desmodus rotundus in Perú and could have increased the levels of exposure to the virus [105, 106]. In Argentina, the extermination of bats changed the direction of spread of rabies in livestock but did not prevent its advancement [107]. In Uganda, miners exterminated a colony of Rousettus aegyptiacus bats after an outbreak of Marburg virus in 2007 that involved 4 miners in close contact with these bats. Five years later, a new outbreak occurred in miners from the same mine. The second time, Marburg virus RNA was detected in a higher proportion in the R. aegyptiacus that recolonized the mine (13.3%, n = 400; [108]) compared to RNA detection before culling in this cave (5.1%, n = 611; [109]) and other caves in Uganda (2.5%, n = 1622; [67]) and Gabon (4.8%, n = 187; [110]) where culling has never been reported.

Culling can also cause demographic changes, leading to a higher proportion of juvenile individuals. This change may occur because of a disproportionate cull of older individuals; the potential increase in survival of pups at lower population densities, followed by higher recruitment of juvenile females into reproductive age [111]; the hypothetical increase of young dispersers immigrating from neighboring colonies into culled, less dense, and better resourced colonies [105, 112]; or by causing compensatory reproduction [113]. This last possibility may have not been studied in bats but seems unlikely due to their high conception rates and usual litter size of one. Examples of younger populations after culling have been reported in the red deer (Cervus elaphus), racoon (Procyon lotor), American mink (Mustela vison), and Australian brushtail possums (Trichosurus vulpecula), among others [114-117]. As our results and past research consistently show higher viral shedding and detection in young individuals, activities leading to a younger bat population are not advisable for viral spillover management. Similar results are expected when fruit bats are culled based on being categorized as “agricultural pests”; therefore, this kind of management may create higher risk of viral exposure to the human population.

Virological, ecological, and epidemiological research on bats over the last 15 years has helped to identify chiropterans as hosts of zoonotic viruses and to document that human-driven environmental change, human behavior, and human-to-human transmission are the key drivers for the creation of bat-human interfaces, spillover, and epidemics of emergent viruses, respectively [118-120]. In the context of the current biodiversity and bat conservation crisis [121, 122], we must not omit these facts when attempting to effectively, and responsibly frame and communicate disease risks associated with bats. Realistic, data-based risk communication is of paramount importance to avoid framing bats as a threat to humans and to support bat conservation given their important ecological roles [123, 124]. With this background, it seems a proper time for the scientific community studying “bat-associated” viruses to move the conversation from bat spillover risk assessments to the planning of pro-biodiversity and subsequently pro-ecosystem strategies aiming to mitigate spillover risk. Science is valued not only for the diagnosis of problems but because it finds solutions to them. Here, we have attempted to aid the progress of scientific and management dialogue by proposing, not only a management strategy to limit potential coronavirus spillover, but one that is context- and logistically-grounded and pro-conservation, promoting the delivery of the key ecosystem services provided by bats.

Conclusions

Data from hundreds of bats collected in East-Africa show that coronavirus shedding is expected to be more frequent when pups are becoming independent from the dam, independently of the age of the bats, their species, their location, and their life histories; however, the odds of shedding do differ by species. These results can guide temporal-based mitigation strategies to prevent bat-associated coronavirus exposure using non-lethal methods in limited-resource settings, where longitudinal surveillance is not feasible, by identifying high-risk periods for coronavirus shedding when contact with bats should be avoided.

Additional file 1: Summary of the inferred start date of the birth pulse, the end of the lactation period, and the start date of the mating period per microchiropteran species.

Additional file 2: Bibliographic references for the traits of the bat species included in the study.

Additional file 3: Imputation of the reproductive season to those bats whose biology is insufficiently known. Model methods, results, and literature cited. Table S3. Summary of the coefficients’ posterior probability distributions of the selected model for the imputation of the reproductive seasons of bats missing this data. Figure S3.1. The discrepancy between the observed (inferred) period and the period predicted by the selected imputation model considering 5 thousand posterior predictions of the observed period. Values of zero indicate matching, (observed period - predicted period). Figure S3.2. Left: the proportion of bats in each period after the imputation for those bats with the non-inferred season. Each line connects the proportion of individuals per reproductive season in each of the 5,000 Markov Chain Monte Carlo iteration. Right: the distribution of the crude coronavirus detection per reproductive season across the 5,000 Markov Chain Monte Carlo sampling iterations after imputing the periods when un-inferred. The black boxplots show the distribution of the coronavirus detection per period, while the light and dark colored boxes above and below show the interquartile detection in non-adults and adult bats, respectively, per period.

Additional file 4: Table S4. Summary of the Alpha- and Betacoronaviruses (alphaCoV and betaCoV, respectively) found in the microbats tested.

Additional file 5: Results of the model to assess the association between CoV shedding and the “Recent weaning” season in Eidolon helvum and microbats. Figure S5.1. Posterior predictive distributions of coronavirus positive bats (histograms) and the observed coronavirus positive bats (vertical lines). The histograms show the distribution of four thousand predictions of detection in the sampled bats. The lines show the observed detection. A) Coronavirus detection across all bats. B) Coronavirus detection across age categories: adults (light blue) and non-adults (yellow). C) Coronavirus detection across the reproductive seasons: “Not recent weaning” (light blue) and “Recent weaning” (yellow). D) Coronavirus detection across the reproductive seasons in the non-adult bats: “Not recent weaning” (light blue) and “Recent weaning” (yellow). E) Coronavirus detection across the life history seasons in the adult individuals: “Not recent weaning” (light blue) and “Recent weaning” (yellow). Figure S5.2. The density of the coefficients’ posterior probability distributions of the selected model. RW refers to the “Recent weaning” period (versus “Not recent weaning” season).

Additional file 6: Sampling events -specific random intercepts values.