A Cross Sectional Sampling Reveals Novel Coronaviruses in Bat Populations of Georgia

Mammal-associated coronaviruses have a long evolutionary history across global bat populations, which makes them prone to be the most likely ancestral origins of coronavirus-associated epidemics and pandemics globally. Limited coronavirus research has occurred at the junction of Europe and Asia, thereby investigations in Georgia are critical to complete the coronavirus diversity map in the region. We conducted a cross-sectional coronavirus survey in bat populations at eight locations of Georgia, from July to October of 2014. We tested 188 anal swab samples, remains of previous pathogen discovery studies, for the presence of coronaviruses using end-point pan-coronavirus RT-PCR assays. Samples positive for a 440 bp amplicon were Sanger sequenced to infer coronavirus subgenus or species through phylogenetic reconstructions. Overall, we found a 24.5% positive rate, with 10.1% for Alphacoronavirus and 14.4% for Betacoronavirus. Albeit R. euryale, R. ferrumequinum, M. blythii and M. emarginatus were found infected with both CoV genera, we could not rule out CoV co-infection due to limitation of the sequencing method used and sample availability. Based on phylogenetic inferences and genetic distances at nucleotide and amino acid levels, we found one putative new subgenus and three new species of Alphacoronavirus, and two new species of Betacoronavirus.


Introduction
Coronavirus (CoV) infection in humans (e.g., 229E, NL63) and other mammals were first reported in the 1960s [1,2]. CoVs are highly pathogenic for livestock, pets, wild animals, as well as birds where they were first described in chickens in the early 1930s [3,4]. During the period 2002-2003, coronaviruses were recognized as zoonotic agents with high pandemic potential after an outbreak of severe acute respiratory syndrome (SARS) that caused 8096 cases and 774 deaths [5,6]. In 2012, the Middle East respiratory syndrome (MERS) emerged as a new public health concern for it spread rapidly to several countries around the globe (mainly the Arabian Peninsula and the Republic of Korea) and for showing Natural vegetation = trees, shrubs and herbaceous cover.
Mist nets and hand nets were utilized depending upon roost type (e.g., caves, buildings, attics). The number of individuals caught per roost and per species was authorized with anticipation by the Ministry of Environment and Natural Resources Protection of Georgia. Bats were individually placed in cotton bags to be transported to the laboratory the day after capture.
Bats were euthanized in a BSL2+ laboratory facility at the NCDC-Lugar Center. Specimens, were subsequently, measured (weight and length), sexed, and morphologically inspected for species identification [20]. Rectal and oral swabs, as well as tissues such as brain, liver, lung, spleen and intestines were collected and stored at −80 • C for subsequent testing. Rectal swabs were the only specimens available to search for coronaviruses for this retrospective study. All tissues were committed to search for other pathogens such as Leshmania spp., Yersinia spp., Leptospira spp., Bartonella spp., Brucella spp., Lyssavirus spp., Hantavirus [17][18][19].

RNA Extraction and PCR Amplification
A total of 188 rectal swabs diluted in 1 mL of virus transport media (VTM) were vortexed and spun down in a refrigerated microcentrifuge for 5 min at 10,000 rpm. For total RNA extraction, 140 µL of the VTM homogenate were placed in a new tube with 560 µL of lysis buffer from the QIAamp Viral RNA Mini Kit (Qiagen, Germantown, MD, USA) according to the manufacturer's instructions. All 188 total RNA extracts were screened with a one-step RT-PCR method described previously that uses degenerated primers targeting a relatively conserved 180 bp region within the RNA-dependent RNA polymerase (RdRp) gene [21]. A second 440 bp non-overlapping fragment was amplified as described by Lelli et al. [21], which together with the 180 bp product would generate an approximately 600 bp sequence that can be used for phylogenetic characterization to infer, genus, subgenus and species with greater confidence [6,9,11,21]. Amplicons were gel purified using QIAGEN MiniElute ® Gel Extraction Kit (250) and subsequently cycle sequenced using the BigDye Terminator Kit, version 3.1 (Applied Biosystems, Foster City, CA, USA) to obtain a consensus sequence. Cycle sequencing products were purified by precipitation with ice-cold 100% isopropanol once, and then with 70% isopropanol. Sequencing was conducted using an Applied Biosystems ABI Prism 3130 XL sequencer with data collection software version 4.0.

Sequence Analyses and Phylogenetic Reconstructions
Excellent quality sequences were obtained in only 9 out of the 46 positive hits (20%) for the 180 bp RT-PCR screening amplicon. Thus, consistent concatenation with their 440 bp counterparts to generate a~600 bp sequence was not possible. Nonetheless, we were able to obtain high-quality sequences ranging from 385 to 440 bp in all 46 CoV positive samples using the 440 bp RT-PCR amplicon. Thus, we trimmed all 46 sequences to 411 bp that corresponded to the longest sequence length we could obtain with GenBank accession numbers OL791325 to OL791370. This 411 bp fragment was used to search for homolog sequences ranging from 65 to 98% nucleotide identity, with 100% query coverage, in the GenBank database using the nucleotide-NCBI-BLAST and MOLEBLAST tools [22]. Gathered GenBank reference sequences were then used to reconstruct preliminary phylogenetic trees for initial genus-level identification (data not shown). Moreover, we used review articles [6,10,11,[23][24][25][26][27] to identify CoV reference sequences representing all currently ICTV recognized subgenera and species pertaining to Alphacoronavirus and Betacoronavirus [9]. We retrieved 164 reference sequences encompassing 14 subgenera and 19 species pertaining to Alphacoronavirus, and 106 reference sequences comprising 5 subgenera and 14 species within Betacoronavirus. All sequences were aligned with the 46 Georgian CoV positive sequences using MUSCLE [28] and trimmed to 411 bp using BioEdit [29]. Taxa were subsequently divided into two datasets, one for BetaCoV that encompassed 106 BetaCoV reference sequences, 27 CoV sequences obtained from Georgian bats and 5 taxa representing some AlphaCov species, while the second data set comprised 164 AlphaCoV reference sequences, 19 sequences from Georgian bats and 4 BetaCoV taxa encompassing different species. MEGA7 [30], was used to determine the most appropriate model of molecular evolution for these combined datasets. The GTR + G + I substitution model was chosen, out of 24 models tested, for our phylogenetic reconstructions based on the Bayesian information criterion (BIC), Table S1 [30]. A Bayesian phylogenetic reconstruction using MrBayes V 3.2 software [31], was run on these two datasets with two independent runs under the GTR + G + I substitution model with 15 million generations, 4 Markov chains each run and the sampling for tree parameters every 1000 generations to assess branch support and/or calculation of Bayesian posterior probabilities. Using the pairwise alignment tool available in BioEdit [29], we calculated both the average nucleotide identity (ANI) and the average amino acid identity (AAI) by comparing one on one sequences of the Georgian CoV clusters (separately, from both the AlphaCoV and BetaCoV alignments) with taxa of their closest subgenera, species, as well as with unidentified closest CoVs neighbors observed in the phylogenetic reconstruction. The subgenera and species demarcation criterion suggested in the current literature were also used to assign Georgian CoV clusters to these taxonomical classification levels [6,9,11,23,27].
To better appreciate the extent of the geographic distribution of AlphaCoV and Beta-CoV found in bat species of Georgia with those bat species harboring highly similar CoVs through Eurasia, we mapped their natural geographic distribution using species range data from the International Union for Conservation of Nature IUCN red list of threatened species [32][33][34][35][36][37][38][39][40][41][42][43][44]. Restricting ranges to the geography currently occupied by each species according to expert's records (i.e., extant-resident). Scientific names and synonymies were corroborated following the standards of the Integrated Taxonomic Information System. Bat species were then grouped based on the genus of the coronavirus detected (i.e., AlphaCoV and BetaCoV). To create the maps, administrative boundaries were generated from maps downloaded from DIVA-GIS. Spatial data were handled and displayed using ArcGIS 10.8 (ESRI 2021 and R R Core Team 2021) [45,46].
A total of 183 bats were collected at 7 locations with different degrees of agricultural perturbation across Georgia, where we found a 25% (46/183) CoV positivity rate ( Table 1). The eighth collection site at Saadamio cave, village Saadamio of Senaki municipality, presented an unperturbed landscape dominated by trees, shrubs embedded in a mosaic of herbaceous cover, where we only caught 5 CoV negative individuals ( Table 1). The overall CoV positive rate for all eight collection sites was 24.5% (46/188), (Table 1, Figure 1).
were generated from maps downloaded from DIVA-GIS. Spatial data were handled and displayed using ArcGIS 10.8 (ESRI 2021 and R R Core Team 2021) [45,46].
A total of 183 bats were collected at 7 locations with different degrees of agricultural perturbation across Georgia, where we found a 25% (46/183) CoV positivity rate ( Table  1). The eighth collection site at Saadamio cave, village Saadamio of Senaki municipality, presented an unperturbed landscape dominated by trees, shrubs embedded in a mosaic of herbaceous cover, where we only caught 5 CoV negative individuals ( Table 1). The overall CoV positive rate for all eight collection sites was 24.5% (46/188), (Table 1, Figure  1). Collected bats belong to two families Rhinolophidae (n = 84) that encompassed a single genus with three species namely, Rhinolophus euryale (n = 40), R. ferrumequinum (n = 39) and R. blasii (n = 5), and Vespertilionidae (n = 104) that comprised 5 genera with a total of 8 species (Table 2). The overall CoV positivity rates within the Rhinolophidae and Vespertilionidae were 25% and 24%, respectively (Table 2). However, Rhinolophus spp. presented higher CoV positive rates for BetaCoV (19%) than for AlphaCoV (6%) (with exception of Rhinolophus blasii that was negative for both CoV genera). In contrast, Miniopterus Collected bats belong to two families Rhinolophidae (n = 84) that encompassed a single genus with three species namely, Rhinolophus euryale (n = 40), R. ferrumequinum (n = 39) and R. blasii (n = 5), and Vespertilionidae (n = 104) that comprised 5 genera with a total of 8 species ( Table 2). The overall CoV positivity rates within the Rhinolophidae and Vespertilionidae were 25% and 24%, respectively (Table 2). However, Rhinolophus spp. presented higher CoV positive rates for BetaCoV (19%) than for AlphaCoV (6%) (with exception of Rhinolophus blasii that was negative for both CoV genera). In contrast, Miniopterus schreibersii exclusively presented AlphaCoV, while Myotis blythii (35% overall positivity rate, with 15% AlphaCoV and 19% BetaCoV) presented similar positivity rates to those observed in Rhinolophidae bats ( Table 2, Figures 2 and 3). The remaining Vespertilionidae species (Eptesicus serotinus, Pipistrellus pygmaeus, Myotis mystacinus, M. alcathoe and Nyctalus leisleri) were CoV negative along this transversal sampling (Table 2, Figure 1). Notably, CoVs were not detected in all bat species in which less than 10 individuals were collected. Albeit different individuals of the same species were found infected with AlphaCoV and BetaCoV within the same locations ( Figure 1, Table 2 Table S2).  The 19 Georgian AlphaCoVs grouped as two unclassified subgenera (n = 12 taxa), and as two unclassified species within two subgenera (n = 7), (Figure 2). The largest cluster encompassing 10 sequences (locate at the top of the tree, marked with a bar in color dark red) collected from cave Gliana, village Kumistavi of Tskaltubo municipality corresponded to Georgian bats of three species, Myotis blythii (n = 7), Miniopterus schreibersii (n = 2), and Rhinolophus euryale (n = 1). This cluster formed a monophyletic clade with unclassified CoV sequences found in Myotis myotis bats from Hungary, Germany, Spain and Italy presenting ANI and AAI values of 97.1% and 99.5%, respectively ( Figure 2, and Table S2). Other CoVs sequences found in a Myotis daubentonii bat from China, an unidentified bat from Korea and a Myotis sp. from Hong Kong were also identified as closely related to this cluster of Georgia sequences. However, they presented lower overall ANI and AAI values around 87.1% and 96.1%, respectively. At a more ancestral node (0.94 posterior probability) this cluster of unclassified AlphaCoV, mainly associated with vespertilionid bats from Europe and Asia shared a common ancestry with Pedacovirus I, Pedacovirus and Colacovirus subgenera with more distant overall ANI and AAI values around 77.6% and 86.6% (Figure 2).
We identified a Georgian CoV sequence (GE_CoV10) obtained from a Myotis emarginatus bat as sister taxon of Rhinolophus ferrumequinum alphacoronavirus HuB-2013 (with an ANI of 94.7 and an AAI of 100%), which falls within the Decacovirus clade ( Figure 2). The fourth AlphaCo cluster (comprising 6 sequences) was the most diverse considering the number of bats species it contained, Miniopterus schreibersii (n = 3) from cave Gliana at village Kumistavi of Tskaltubo municipality, Rhinolophus ferrumequinum (n = 1) from the managed reserve at Gardabani, R euryale (n = 1) from cave Gliana at village Kumistavi of Tskaltubo municipality, and Myotis blythii (n = 1) from cave Gliana at Tshaltubo. This Georgia cluster was monophyletic with sequences pertaining to Myotis ricketti alpha coronavirus Sax-2011 and Miniopterus bat coronavirus HKU8 bat CoV, as well as with two CoV obtained from Molossus rufus bats from Brazil, all occupying the Myotacovirus clade (Table S1, Figure 2). The closest relatives (albeit they grouped inside a polytomy with a low support value 0.57) to this cluster of Georgian bat CoVs were sequences obtained from Myothis blythii from Kazakhstan (91% ANI and 100% AAI), Miniopterus schreibersii (90% ANI and 99.3% AAI) from Luxemburg and Spain, and Myotis emarginatus (89.8% ANI and 97.8% AAI). These Georgian sequences, together with all other sequences inside the polytomy shared a common ancestor with CoV sequences obtained from Molossus rufus from Brazil (also pertaining to Myotacovirus), with ANI (77.8%) and AAI (83.2%) values observed (Figure 2).
All CoV positive bat species that were found in Georgia presented an overlapping geographic distribution across most of southern Europe, a portion of North Africa and across central Asia. Only the distributions of Myotis blythii and Rhinolophus ferrumequinum's extended to parts of Southeast Asia with the latter having the widest distribution in the region ( Figure 4). Myotis ricketti (renamed as Myotis pilosus) whose natural distribution seems to be restricted to Southeast Asia appears to be the primary host of Myotis ricketti alpha coronavirus Sax-2011. However, close relatives of this AlphaCoV species were found circulating in Georgia in Miniopterus schreibersii, Rhinolophus ferrumequinum and Myotis blythii, suggesting that broadly distributed bats species such as Rhinolophus ferrumequinum and Myotis blythii could have dispersed this CoV species to Central Asia and Western Europe (Figure 4).

Discussion
The genetic characterization of mammal-associated CoV using a region spanning 411 bp within the RdRp was robust enough to identify novel CoVs circulating in bat populations across Georgia, albeit no reliable evolutionary relationships among subgen-

Discussion
The genetic characterization of mammal-associated CoV using a region spanning 411 bp within the RdRp was robust enough to identify novel CoVs circulating in bat populations across Georgia, albeit no reliable evolutionary relationships among subgenera and species could be established consistent with previous investigations using RdRp fragments shorter than 806 nucleotides [6,9,23,27]. Overall, our analyses based on phylogenetic inference, as well as considering the ANI and AAI criteria proposed recently [6,9,27], demonstrated the circulation of one putative unclassified subgenus closely related to the Pedacovirus holotypes and the Colacovirus subgenus. Three unclassified species two of which are within the Decacovirus, and the other one within the Myotacovirus. Additionally, we found two putative unclassified species within Sarbecovirus and Merbecovirus.
Wilkinson et al.; provided robust evidence that strong monophyletic groups demarcating a subgenus and species should have Bayesian posterior probabilities higher than 0.9 [6]. The genera level trees constructed for this investigation depicted with high support values all ICTV recognized subgenera and species and were consistent with previous investigations where robust analyses were conducted [6,9,11,27]. Interestingly, most of our Georgian CoV clades grouped consistently within CoVs associated with the same bat species or genus across Europe, Asia, Africa or the Americas, which supports the hypothesis that CoVs have strong epizootiological associations with their bat hosts spanning most of their natural geographic distributions [9,16,23,[47][48][49][50][51][52]. In some instances, there were multiple bat species with the same viruses such as in the SARS-like clade where Rhinolophus spp. and Myotis blythii shared closely related viruses [11,51,53]. Similarly, two different clusters within AlphaCoVs, one that was embedded within Myotacovirus clade from Georgia was dominated by sequences obtained from Miniopterus schreibersii and Rhinolophus spp. However, this group contained a CoV species mainly associated with Myotis ricketti, while highly similar CoV sequences have been obtained from Myotis blythii in Kazakhstan and Italy and other parts of Europe, suggesting Myotis spp. could be the main reservoir hosts for this CoV [49,[53][54][55]. Another Georgia clade predominantly containing sequences from Myotis blythii, may represent a yet unclassified Myotis blythii alphacoronavirus subgenus, which was found also infecting Rhinolophus euryale and Miniopterus schreibersii bats. Interestingly, the cross-sectional sampling we undertook could not identify an infection pattern among sedentary (Rhinolophus ferrumequinum, R. euryale, R. blasii, Myotis emarginatus, M. alcathoe and Eptesicus serotinus), sedentary or short-distance migrants (M. blythii, M. mystacinus), short-distance seasonal migrant (Miniopterus schreibersii), partially migrant (Pipistrellus pygmaeus) and long-distance migrant (Nyctalus leisleri) bat species, that would intuitively render seasonal and long-distance migrant bat species with a greater diversity of CoVs [39,40]. Although Miniopterus schreibersii bats have been found infected either with alpha or beta CoVs throughout Europe, they are predominantly infected with AlphaCoV in South East Asia [47,51,56], as we noted in Georgia. All together, these results indicate relatively frequent CoV spillover infections among bat species sharing roosts highlighting a broad susceptibility of multiple bats species to CoV infection that could potentially lead to host shifts [23,27]. Alternatively, these results may suggest that there might be multiple bat hosts species associated with the circulation maintenance of each viral genus acting as host-community complexes. Bat species with broad natural ranges could be playing a central role in the gradual dissemination of alpha and beta CoVs to other susceptible bats species across Eurasia, making CoV infection persistence and its dissemination dynamics a rather complex issue across bat populations of World.
Wilkinson and colleagues proposed an ANI threshold of 77.6% for different subgenera within alpha CoVs and 71.7 for BetaCoVs with high confidence levels [6]. Conversely, Geldenhuys and colleagues report pairwise amino acid average divergence thresholds (that we converted to identity values using the following formula, 100-divergence value = identity) that demarcate species to subgenus with AAI values lower than 92.4% and from subgenus to genus with AAI values lower than 85.3% [27]. Thereby, all our inferences to characterize Georgian CoV to subgenus and species levels fell into expected ranges [6,27].
Our results are not surprising since Georgia contains fauna native to both Europe and Asia and is relatively close to the Middle East [15]. Its geographic location at the middle of this great bio-diversity corridor, explains the outstanding diversity of CoV found, despite its relatively small territory (69,700 km 2 ), as well as predicts a similar CoV richness for all the unexplored countries in the surrounding region [16]. We noticed that all Georgian bat species found positive with CoVs had broad natural geographic distributions across Eurasia, parts of Africa and the Americas. Strikingly, Georgian CoVs share recent common ancestry with CoV found circulating in other widely distributed bats species across Eurasia and the Americas. This empirical observation suggests mammal-associated CoVs may undergo complex transmission dynamics across sympatric bat communities, perhaps highlighting the susceptibility of multiple bat species to CoV infection. Recent modeling efforts have assessed the risk for the presence of bat-associated zoonoses globally [56]. Some of these results suggest Western Asia (especially the region comprising Georgia) is a hot spot for the presence of BetaCoV with pandemic potential, outside of South East Asia Western Europe and Central Asia [15,16].
The great CoV diversity found in single caves and the close proximity among positive locations suggests CoV co-infections are likely to occur in bat populations across the region. However, the sequencing approach we used and the amount of sample available hampered the detection of CoV co-infections in our sample set. In addition, the poor sequencing results obtained in the 180 bp pan-coronavirus RT-PCR could be due to a poor amplicon yield associated to low primer affinity (due to more than two sequence mismatches) along the highly variable region of the RdRp gene this 180 bp RT-PCR targets. This is clearly indicated by the 7 degenerated positions in the primer set used to amplify this 180 bp region [21]. Thus, this short 180 bp amplicon is practically useless for further CoV phylogenetic characterization. For future research we would only recommend the use of the 440 bp long amplicon (that anneals along a much more conserved region within the RdRp), which could be used for screening and phylogenetic characterization when a limited amount of sample is available [21]. Moreover, the sequencing of such an amplicon by next generation platforms may allow the detection of CoV co-infections in single bats.

Conclusions
CoVs found in Georgian bat populations were strikingly diverse. Some of these Georgian CoVs were closely related to those that caused human pandemics with epicenters in China (SARS) and the Middle East (MERS). We predominantly described new bat CoV species within previously described subgenera within Alphacoronavirus and Betacoronavirus, with one likely novel, yet unclassified, a subgenus of Alphacoronavirus associated with Myotis blythii. Our results provide further insight into the global genetic structure of CoV associated with bats in Eurasia and may inform the overall regional CoV diversity, critical to assess the global epidemiology of emerging diseases.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/v14010072/s1, Tables: Bat species captured and associated coronaviruses in the present study. The GTR + G + I substitution modeling and phylogenetic recon-structions based on the Bayesian information criterion (BIC). Tables S1-S4: Shows all the metadata associated to each of the samples collected in Georgia, which include collection date, location's name, sample number, identified bat species, taxon's name, morphometric characteristics corresponding to each bat collected. Putative CoV genus identified according to the phylogenetic inference, BLAST best match and overall nucleotide identity over a 411 bp fragment. Red color rows indicate taxa identified as positive for alphacoronavirus RNA and blue color rows indicate taxa identified as positive for beta CoV RNA via RT-PCR.  Data Availability Statement: The data that have been used in this study will be available from the corresponding author upon requests.