Genetic diversity and ecology of coronaviruses hosted by cave-dwelling bats in Gabon

Little research on coronaviruses has been conducted on wild animals in Africa. Here, we screened a wide range of wild animals collected in six provinces and five caves of Gabon between 2009 and 2015. We collected a total of 1867 animal samples (cave-dwelling bats, rodents, non-human primates and other wild animals). We explored the diversity of CoVs and determined the factors driving the infection of CoVs in wild animals. Based on a nested reverse transcription-polymerase chain reaction, only bats, belonging to the Hipposideros gigas (4/156), Hipposideros cf. ruber (13/262) and Miniopterus inflatus (1/249) species, were found infected with CoVs. We identified alphacoronaviruses in H. gigas and H. cf. ruber and betacoronaviruses in H. gigas. All Alphacoronavirus sequences grouped with Human coronavirus 229E (HCoV-229E). Ecological analyses revealed that CoV infection was significantly found in July and October in H. gigas and in October and November in H. cf ruber. The prevalence in the Faucon cave was significantly higher. Our findings suggest that insectivorous bats harbor potentially zoonotic CoVs; highlight a probable seasonality of the infection in cave-dwelling bats from the North-East of Gabon and pointed to an association between the disturbance of the bats’ habitat by human activities and CoV infection.

Coronaviruses (CoVs) belonging to the Coronaviridae family are viruses known to infect a wide range of animals and humans. In humans, CoVs are responsible for mild to severe respiratory illnesses such as the Severe Acute Respiratory Syndrome (SARS) and the Middle-East Respiratory Syndrome (MERS). The epidemic of SARS which started in 2002 in China, and spread on various other continents such as North America and Europe, reached a mortality rate of 9% 1 . The other severe epidemic illness (MERS) appeared more recently in the Middle East, and like the SARS, it also spread in other countries in Africa, America and Europe with case fatality rates of 35% (reviewed by de Wit et al. 2 ). Currently the world is facing an ongoing pandemic caused by a new coronavirus (SARS-CoV-2) that emerged in China in December 2019 and caused disease  resulting in 118 326 confirmed cases and 4292 deaths as of March 11, 2020 3 . Since the emergence of these respiratory diseases, CoVs have been considered a real public health problem because of their ability to become epidemic and pandemic 4 . Even though human to human transmission of SARS-CoV and MERS-CoV occurs mainly through nosocomial transmission, some human cases recorded were occasionally the results of a zoonotic transmission 5 . This raised an interest in the identification of animal reservoirs. Indeed, it appears that almost each human CoV have zoonotic origins or otherwise have a close relative that circulate in wild animals (bats) [6][7][8] and domestic animals (camels and cattle) 9,10 . Besides SARS-CoVs and MERS-CoVs, numerous other CoVs have been detected in bats in Africa, Asia, Europe and America and are classified into the genera Alphacoronavirus and Betacoronavirus 11,12 . For SARS-CoV-2, a recent study revealed that the whole genome sequence of SARS-CoV-2 was most closely related with a bat coronavirus detected in bats from Yunnan Province in China 13 . Consequently, bats were defined as reservoirs of ancestral coronaviruses which were transmitted to humans 14,15 . Other wildlife species appear to Molecular characterization of identified CoVs. Fifteen bat nucleotide sequences were obtained and compared to those available in the public database using the algorithm "Blastn" of NCBI BLAST. The results indicated that all the nucleotide sequences from this study matched the deposited CoVs sequences. In order to determine the genetic relationships between the 18 bat CoVs from this study and previously described CoVs, phylogenetic analysis was performed based on 327-bp RdRp truncated sequences (Fig. 2a). The results showed that 12 out of 15 sequences belonged to the Alphacoronavirus (α-CoV) genus (Fig. 2b), whereas the remaining three belonged to the Betacoronavirus (β-CoV) genus (Fig. 2c). Among the 12 sequences of α-CoVs, 7 were obtained from individuals of H. cf. ruber caught in 2009 in the Faucon cave, 2 from other individuals of the same species but caught the same year in the Batouala cave; the other 3 sequences were obtained from 2 H. cf. ruber and 1 H. gigas, all caught in 2010 in the Faucon cave. The 3 sequences of β-CoVs were detected only in H. gigas.
The sequences of Alphacoronavirus were 91.9-100% identical with each other. The 7 sequences from 2009 from H. cf. ruber individuals of the Faucon cave showed 100% identity with each other. Regarding the 2 sequences of Alphacoronavirus from H. cf. ruber of the Batouala cave, the 09GB0761 (MG963198) sequence displayed either 96.5% or 96.8% with the 7 other sequences, while the 09GB0809 (MG963189) sequence showed 100% identity with the sequences from H. cf. ruber of the Faucon cave. Phylogenetic analysis supported this finding (Fig. 2b) and also showed that these 7 sequences clustered with 2 bat α-CoVs sequences (JX174639 and JX174640) from 2 H. cf. ruber from the Faucon cave caught in 2009. These two bat sequences showed 100% identity at nucleotide level with the 7 others. The 3 sequences of α-CoV (10GB0354, 10GB0309 and 10GB0318) from 2010 were 93.2-99.7% identical with each other. The 10GB0309 (MG963199) and 10GB0318 (MG963201) sequences, identical at 99.7%, from individuals of the same species (H. cf. ruber) caught in the Faucon cave, diverged respectively by 6.4% and 6.8% from the 10GB0354 (MG963200) sequence from an H. gigas captured in the same cave.
The 3 sequences of β-CoV, 13GB0214 (MG963186), 13GB0215 (MG963187) and 13GB0273 (MG963188) were obtained from 3 H. gigas caught in 2013 in the Faucon cave. The 13GB0214 and 13GB0215 sequences displayed 99.7% identity at nucleotide level. These two sequences shared a nucleotide identity of 96.5% and 96.9% respectively with the 13GB0273 strain. In addition, these 3 sequences formed a distinct cluster with a bat β-CoV named the Zaria bat coronavirus (ZCoV) (accession number HQ166910), detected in Nigeria in 2008 in a species of bat genetically close to H. gigas, Hipposideros commersoni. The nucleotide identities between ZCoV and 13GB0214, 13GB0215 and 13GB0273 were 96.8%, 96.5% and 97.8%, respectively.
Bats of the species H. cf. ruber and H. gigas positive for alphacoronaviruses shared the same caves with the species Coleura afra and Rousettus aegyptiacus, which were negative for alphacoronaviruses. Moreover, the three H. gigas bats from the Faucon cave caught in 2013 (Table 3), shared the same cave with bats of the species C. afra and M. inflatus in which no Betacoronavirus was detected. Only the H. gigas species was infected with a Betacoronavirus.
High throughput sequencing yielded about 9.9 million reads (Fig. 3). Quality filtering discarded 86 955 (0.87%) sequences. No reads related to coronavirus could be detected with EDGE's taxonomy classification module, but BLAST alignments (98-100% and 90-95%, for respectively query coverage and identity percentage) revealed 2 small coronavirus contigs (210 bp and 271 bp) from one sample. Those contigs have been determined as being parts of replicase polyprotein 1ab. Continued www.nature.com/scientificreports www.nature.com/scientificreports/ Factors driving CoV positivity in bats. The comparison of the proportions of infected bats by category for each variable showed that the differences observed between bat genera (Fig. 4a) were statistically significant (p-value = 0.00002804), as well as the variations within species (p-value = 0.00003354) (Fig. 4b), ages (p-value = 0.0000355) (Fig. 4c), months of capture (p-value = 0.00003423) (Fig. 4e) and sites (p-value = 0.0004276) (Fig. 3f). There was no significant difference between the 2 sexes (Fig. 4d). Furthermore, www.nature.com/scientificreports www.nature.com/scientificreports/ the MCA on all the bats revealed that the most important variables were the month of capture and the site followed by the species. In infected bats, the most important variables were also the month of capture, followed by the site and species (Fig. 5a). Figure 5b and Table 4 show the association between the categories of the variables and the characteristics of the different groups, respectively. According to these results, it appears that CoV infection in H. gigas is correlated with the months of July and October. In addition, infection with CoVs in the H. cf. ruber species can be associated with the months of November but also October. It appears that CoV infection is particularly associated with the Faucon cave.

Discussion
We studied 1867 animals from 39 species including rodents, bats, primates and other wildlife species, from 6 provinces (Estuaire, Haut-Ogooué, Ogooué-Ivindo, Ogooué-Lolo, Ngounié and Woleu-Ntem) in Gabon. This study identified coronaviruses (CoVs) in 18 bats belonging to the species H. cf. ruber, H. gigas and M. inflatus from two caves located in the Ogooue-Ivindo province located northeast of Gabon. The nucleotides sequences obtained showed that these CoVs belonged to the Alphacoronavirus and Betacoronavirus genera. In a previous   www.nature.com/scientificreports www.nature.com/scientificreports/ study, the identification of CoVs was already reported in H. cf. ruber in Gabon 17 . Taken together, these results indicated that genetically distinct coronaviruses co-circulate among bats in the Ogooue-Ivindo province. These findings supported that bats carry a great genetic diversity of CoVs and Woo et al. 19 suggested that bats were ideal hosts for both alphacoronaviruses and betacoronaviruses and could play an important role in the ecology and evolution of coronaviruses. Although no coronavirus was identified in other animal species, alpha-and betacoronaviruses were already found naturally in various wild mammals such as civets 5 , buffalos 20 and rodents [21][22][23] . In their studies in rodents from China, Lau et al. 21 and Wang et al. 22 analyzed 725 and 1465 rodents, respectively, whereas we analyzed only 494 rodents. The lack of detection of CoVs in our study could be linked to the small sample size, unlike the aforementioned studies. However, in their study, Ge et al. 23 found 23 infected rodents out of the 177 studied. Unfortunately the rodents in our study come mainly from urban localities, rodent captures should also be considered in rural forest areas around caves inhabited by bats. Moreover, experimental infections have been described in non-human primates 24,25 but no natural infection of CoV in NHP have yet been reported. For other wild animals, the small size of samples investigated by species could explain the lack of CoVs found in this study.
Alphacoronaviruses and betacoronaviruses were detected only in H. cf. ruber and H. gigas, 2 among the 5 species tested. These results from the Ogooue-Ivindo province suggest that these 2 species could be both natural reservoirs of and specific for the alpha-and the betacoronavirus, respectively. A close co-evolutionary association could be suspected between CoVs, particularly betacoronaviruses, and bats of this family 26 . Interestingly, in our study, we detected a specific Hipposideridae Betacoronavirus clustering with a Betacoronavirus previously detected in Thaïland in 2007 26 . It was suggested that phylogenetic relationships between betacoronaviruses are mainly driven by the host phylogeny 26 . However, in a study conducted in Africa, Quan et al. 27 suggested that there was no strict association between bat species and betacoronaviruses. Besides this study, previous studies on CoV host species in Asia have suggested that there was no strict association between bat species and alphacoronaviruses 28,29 . Indeed, many studies have described coronaviruses in a wide range of bat species 30,31 . Nevertheless, despite the report of alphacoronaviruses in numerous species of Miniopteridae and Vespertillionidae, Ar Gouilh et al. 32 suggested a strict association between Alphacoronavirus EPI6 and Rhinolophidae, a familly usually associated with betacoronaviruses.
No CoVs were detected in chiropteran species sharing the same caves with the positive species (H. cf. ruber and H. gigas), suggesting a specificity of CoVs for these species or an absence of contact or interaction between some species. It was observed, for example, in the Faucon cave, that colonies of H. cf. ruber and H. gigas were found separately, in specific areas of the cave, the same with the other species (M. inflatus and C. afra), except for R. aegyptiacus which was not present in this cave. This spatial segregation could limit direct contact between species and represent a behavioral barrier to interspecific transmission, as previously described for lentiviruses and felids 33,34 and more recently for bats and coronaviruses 32 .
A strain of Alphacoronavirus was identified in 7 H. cf. ruber nesting in the same cave (Faucon cave) and caught the same year (2009). Likewise, 2 other H. cf. ruber from the Batouala cave but caught in 2010 were infected by the same alphacoronavirus infecting the 7 H. cf. ruber from the Faucon cave. Corman et al. 14 had already identified a closely related strain of Alphacoronavirus in 2 H. cf. ruber from the Faucon cave in 2009. Altogether, these findings first suggest the transmission of CoVs between individuals of the same species through close contacts probably increased by social and roosting behavior, large colony size and reduced space in these caves 26 ; this strain of Alphacoronavirus could be enzootic in H. cf. ruber in this area of Gabon and roosts harboring a large number of bats nesting together promotes viral diffusion in the colony 29,35 , and second, local movements of bats individuals in the region imputable to seasonal/reproductive or opportunistic roosting behaviors, as the two caves being only 67 km apart. The migratory behavior of some bats provides an opportunity for pathogens to spread over long distances, as reported for flying foxes of the species Eidolon helvum able to migrate more than 1000 km 36 38 . In Ghana, Pfefferle et al. 6 identified a strain of Alphacoronavirus in H. cf. ruber closely related to the human coronavirus 229E. Our finding argues that close relatives of the human coronavirus 229E exist in African bats 14 . The phylogenetic proximity of the CoV strains identified in this study with the human CoV-229E strain suggests a risk of emergence of this virus in humans in Africa. Furthermore, the 10GB0354 alphacoronavirus strain (MG963188), from H. gigas bats caught in 2010 in the Faucon cave, was phylogenetically close (93.8% nucleotide identity) to the Alpaca CoV strain detected in the Vicugna pacos species, commonly known as alpaca (JQ410000). Maganga et al. 17 and Corman et al. 14 already described two Gabonese CoV strains from hipposiderid bats caught in 2009 in the Faucon cave closely related to Alpaca CoV. Indeed, this coronavirus was described in 2007 during the epidemic of severe respiratory disorders associated with abortions in alpacas. It could be the direct progeny of its common ancestor with HCoV-229E, sharing with the latter a nucleotide www.nature.com/scientificreports www.nature.com/scientificreports/ identity of 92.2% 39,40 , on 440 bases of the nsp12 (polymerase). The phylogenetic link between these 2 CoVs could originate from an inter-species transmission, as a result of the opportunistic evolution of their common ancestor. As a result of genetic mutations and recombinations, the Alpaca CoV would have passed from bats to alpaca and vice versa. Phylogenetic analysis also revealed that two betacoronaviruses from 2 H. gigas, caught in 2013 in the Faucon cave, were phylogenetically close to a bat betacoronavirus named the Zaria bat coronavirus (ZCoV), detected in Nigeria in 2008 in a species of chiroptera genetically close to H. gigas, H. commersoni 27 , suggesting that phylogenetic proximity between host species may promote inter-species transmission 41 .
BLAST succeded in finding 2 contigs linked to coronavirus while no reads from this species could be detected by all EDGE's taxonomy tools. To identify coronavirus reads having being assembled into those 2 contigs, further analyses using BLAST against NCBI RefSeq virus database were carried, allowing for the identification of 7    www.nature.com/scientificreports www.nature.com/scientificreports/ coronavirus reads. Since BLAST proved to be more sensitive, the same was done for the 2 other samples in order to verify that no coronavirus reads could have been missed. No more reads related to coronavirus were found. Moreover, this analysis overall showed that viral reads were extremely under-represented. One could argue that not finding the targeted virus might be caused by low DNA concentration in sequencing libraries. However, even in the sample where the highest concentration has been recorded, no coronavirus reads were found. Furthermore despite having achieved a satisfying quality sequencing, the biased amplification (majority species will be over-represented) and the low viral load explain why coronavirus was identified with such weak coverage.
Knowledge on the diversity of coronaviruses in bats and the factors promoting this diversity is important to prevent the risks of emergence in humans. Based on the analysis of data on life history traits of bats collected in the field, our findings suggest that the diversity of bat coronaviruses was associated with the seasonality and the site of capture, as well as the species of bats. It appeared that CoV infection in H. gigas correlated with the months of July and October, while CoV infection in H. cf. ruber was associated with the months of October and November.
In tropical ecosystems, some bat species seem to breed year-round, whereas others have annual or biannual birth peaks 42 . In Gabon, the parturitions in H. cf. ruber occur synchronously but some colonies give birth in March and others in October. In H. gigas from the Republic of Congo, a neighboring country, the young are born at the beginning of the rainy season. In Gabon, the months of October and November correspond to the beginning of the rainy season. A study demonstrated a high prevalence of CoVs in Chinese bats during the beginning of the rainy season 43 . The same finding was reported for the shedding of astroviruses in insectivorous bats in Asia (Borneo) 42 . Moreover, these months also correspond to the parturition season for H. gigas and H. cf. ruber. Indeed, lactating females of H. cf. ruber and juveniles of the same species were observed in October in northeastern Gabon. Regarding H. gigas, these observations are more difficult, Brosset and Saint Girons 44 reported that the seasonal events of mating and parturition in H. gigas were confined in large caves where these large bats are entirely dependent to raise their young. According to Plowright et al. 45 , reproductive stress was an important driver of Hendra virus antibody prevalence. The prevalence of antibodies against the Hendra virus was higher in pregnant and lactating female of Pteropus scapulatus contrary to non-reproductive females. Different studies showed that parturition and lactation were important risk factors for CoVs shedding in insectivorous bats 46,47 . These two physiological states are energetically costly and thus lead to a depression of the immune system.
Moreover, our findings suggested that CoV infection was particularly associated with the Faucon cave. Of all the caves studied, the Faucon cave was the largest, with many cavities sheltering a diversity of bat species living in sympatry. This cave is regularly visited by villagers to hunt bats for human consumption. This hunting pressure could cause habitat disturbance and stress to these animals. Although Seltmann et al. 42 found no association between habitat disturbance and the detection rate of CoVs in insectivorous bats, numerous studies in bats reported that pathogen prevalence was associated with habitat disturbance 48 . Many authors argued that habitat disturbance in some bats species may cause chronic stress and consequently immunosuppression 42,49-51 , increasing the susceptibility of animals to acquire and shed viruses.
In summary, the results of our study highlight a genetic diversity of CoVs in insectivorous bats in northeastern Gabon and show that the beginning of the rainy season and mainly the parturition season constitute a period of risk for CoV infection of insectivorous bats H. cf. ruber and H. gigas. These results also suggest an association between the disturbance of the bats' habitat by human activities and the infection of these animals with these viruses. This study highlights a probable seasonality of the infection of insectivorous cave bats of the north-east of Gabon by CoVs. To conclude, a longitudinal monitoring should be carried out at different months of the year in order to determine the reproduction periods of the different bat species and therefore the periods favorable to viral shedding. The effect of the habitat disturbance on viral shedding should also be studied further. Bats harbor a large number of potentially zoonotic CoVs and pose a threat to public health for populations that are in contact with reservoir species of these viruses. The knowledge of CoVs ecology in bats is essential to prevent emergence in human populations.

Methods
Ethical approval. All experiments were performed in accordance with relevant guidelines and regulations. www.nature.com/scientificreports www.nature.com/scientificreports/ liver, spleen and gut, were then kept in liquid nitrogen until their delivery at the CIRMF (Centre International de Recherches Médicales de Franceville) and stored at −80 °C for further laboratory investigations.
In the Ogooué-Lolo province, bats were caught in the Kessipoughou cave (0°86 S 12°77 E); and in the Haut-Ogooué province, bats were captured in the Djibilong cave (1°36 S 13°46 E), located north of Franceville. For bats, sex, stage of maturity or age (adult, subadult, juvenile) and capture season were recorded. Rodents were sampled in 2012 and 2013, using Tomahawk and Sherman traps as previously described 54 . The traps were set in houses and their vicinities in three provinces of Gabon: Estuaire (Libreville and Owendo), Ogooué-Ivindo (Makokou) and Haut-Ogooué (Franceville and Leconi). After autopsy, organs were collected in bats and rodents, frozen and then stored at −80 °C until analysis. Bat and rodent species were identified by trained field biologists.
Feces collection. Fecal samples were collected from NHPs, bats and wild animals. Fresh bat droppings were collected on plastic film arranged below roost sites in the caves 55 . Additionally, fecal pellets were directly collected in bags in which the bats were put individually after capture. Fecal samples from NHPs were collected in the field in the Ogooué-Ivindo province (Lope National Park) in 2015. Feces from other wild animals were sampled in 2015 in the Haut-Ogooué (Lekedi Park in Bakoumba) and Ogooué-Ivindo (Lope National Park and Mopia) provinces. Fresh fecal samples were preserved in approximately 10 ml of RNAlater (Ambion) or frozen and transferred to the CIRMF laboratory, where they were stored at −80 °C until analysis.
Specimen preparation and RNA extraction. Total RNAs were extracted from the intestines of bats and rodents coming from bushmeat, and from the feces of bats, PNHs and other wild animals. Around 100 mg of intestine was pooled according to the species and crushed in 600 µl of cold phosphate buffer saline (PBS) (Biological Diagnostic Supplies Ltd) in a ball-mill tissue grinder (Geno/Grinder 2000, Spex Centripep). Total RNA was extracted from 100 µl of homogenate supernatant using a Biorobot EZ1 and the EZ1 RNA tissue mini kit (Qiagen) according to the manufacturer's guidelines in a BSL3. For feces preserved in RNAlater solution, RNAlater was taken out before feces was suspended in PBS and centrifuged at 1500 rpm for 5 min. RNAs were then extracted from 400 µl of supernatant, using the EZ1 RNA Viral Mini kit (Qiagen) according to the manufacturer's recommendations.
The RNA was then tested for the detection of CoVs using a nested reverse transcription-polymerase chain reaction (nRT-PCR) 56  Hight-throughput sequencing. The hight-throughput sequencing was performed using a MiSeq machine.
Phylogenetic analyses. The partial RdRp sequences obtained were cleaned and the contigs assembled using ChromasPro 1.5. The nucleotide sequences were compared to those available in GenBank using the "Blastn" algorithm of NCBI BLAST. Multiple sequence alignments were performed using the ClustalW algorithm implemented in the MEGA 7 software 59 . Phylogenetic analyses were performed by maximum likelihood algorithm using 1,000 bootstraps replications with PhyML 60 and also under a Bayesian statistical framework implemented in BEAST (version 1.8.3) 61 using the general time reversible model of substitution, with a gamma distribution and a proportion of invariant sites (GTR + I + G).
Nucleotide sequence accession numbers. All the sequences obtained in this study have been deposited in Genbank under the accession numbers MG963186-MG963189 and MG963191-MG963201.