Identification of novel bat coronaviruses sheds light on the evolutionary origins of SARS-CoV-2 and related viruses

Identification of novel bat coronaviruses

Between May 2019 and November 2020, a total of 283 fecal samples, 109 oral swabs, and 19 urine samples were collected from bats in a tropical botanical garden and adjacent areas in Mengla county, Yunnan province, southern China. The majority of samples were collected from horseshoe bats, comprising Rhinolophus malayanus (n = 88), R. stheno (n = 36), R. sinicus (n = 34), R. siamensis (n = 12), R. pusillus (n = 2), other Rhinolophus sp. (n = 11), and Hipposideros larvatus (n = 59) (Figures 1A and 1B, Table S1). These samples were pooled into 100 libraries (numbered p1 to p100) according to collection date and host species, with each library containing 1 to 11 samples. Meta-transcriptomic (i.e., total RNA) sequencing was performed and coronaviruses contigs were identified in 40 libraries (Table S2). Blastn searches of the de novo assemblies identified 26 long contigs (>23,000 nt in length) that mapped to coronavirus genomes present in 20 libraries, including 9 sarbecoviruses (i.e., from the genus Betacoronavirus) and 17 alphacoronaviruses. The number of read-pairs mapping to these long contigs ranged from 3,433 to 21,498,614, with the average depth ranging from 35.86 to 215,065.00 (Table S3). It should be noted that pool p1 comprising 11 fecal samples from R. malayanus was the same pool previously used to identify the viruses RmYN01 and RmYN02 (Zhou et al., 2020a). The remaining 24 genomes were named in the same manner, in which the first two letters represent an abbreviation of the bat species, YN denotes Yunnan, and the final number is a serial number ranging from 03 to 26. In addition, several short contigs related to SARS-CoV-2 were identified in two other libraries, p7 and p11 (Figure S1 , Table S2).

Further Blastn analyses revealed that four of the seven novel sarbecoviruses identified here (RpYN06, RsYN04, RmYN05, and RmYN08) were related to SARS-CoV-2, while the remaining three (RsYN03, RmYN07, and RsYN09) were more closely related to SARS-CoV. We next designed specific primers and a probe set of quantitative real-time PCR primers (qPCR) (Table S4) that targeted the conserved region of the 1a gene region to detect the presence of the four SARS-CoV-2-related viruses in individual bats (i.e., prior to sample pooling; Figure 1C). Pool p46 only contained a single positive fecal sample, no. 379, collected on May 25, 2020, and the virus was detected with a cycle threshold (Ct) value of 26.97 (Figure 1C). SARS-CoV-2-related viruses were also detected in three (sample nos. 362, 364, and 372) of the six, three (sample nos. 367, 391, and 397) of the eight, and two (sample nos. 448 and 450) of the seven samples in pool nos. p35, 44, and 62, respectively, with Ct values ranging from 26.10 to 32.82 (Figure 1C). Among these, samples 362, 364, 372, and 367 were collected on May 25, 2020, 391 and 397 were collected on June 3, 2020, while both 448 and 450 were collected on July 16, 2020. The 5′ and 3′ termini and the spike gene sequences of the four coronaviruses related to SARS-CoV-2 were verified using individual samples 379, 364, 367, and 450 with 5′ and 3′ RACE (Table S4) and Sanger sequencing. Results from the Sanger sequencing were consistent with those obtained from the meta-transcriptomic sequencing.

Sequence identities between SARS-CoV-2 and related viruses

At the scale of the whole genome, RpYN06 exhibited 94.48% sequence identity to SARS-CoV-2, making it, after RaTG13 (96.10%), the second closest relative of SARS-CoV-2 documented to date (Figure 2 ). However, because of extensive recombination, patterns of sequence similarity vary markedly across the virus genome, and RmYN02 shared 97.18% sequence identity with SARS-CoV-2 in the 1ab open reading frame (ORF), compared to 97.19% for RpYN06. In addition to the ORF1ab, RpYN06 shared the highest nucleotide identities with SARS-CoV-2 in the RdRp (RNA-dependent RNA polymerase; 98.36%), ORF7a (96.72%), ORF8 (97.54%), N (97.70%), and ORF10 (100%) (Figure 2). However, RpYN06 exhibited only 76.33% nucleotide identity to the SARS-CoV-2 spike gene and 60.91% in the receptor binding domain (RBD), in a manner similar to RmYN02, ZC45, ZXC21, and the Thailand coronavirus strains (Figure 2). Excluding the spike gene, the sequence identities of RpYN06, RmYN02, and RaTG13 to SARS-CoV-2 were 97.17%, 96.41%, and 96.49%, respectively.

In contrast, RsYN04, RmYN05, and RmYN08 exhibited >99.96% nucleotide identities to each other at the scale of the whole genome. Such strong similarity is indicative of viruses from the same species, even though they were sequenced on different lanes and the samples were collected from different bat species at different times. In addition, they shared relatively low sequence similarity with SARS-CoV-2 across the whole genome (76.5%), particularly in the spike gene, ORF3a, ORF6, ORF7a, ORF7b, and ORF8 with nucleotide identities <70% (Figure 2). Interestingly, when using RsYN04 as the query sequence, the closest hit in the Blastn search was the pangolin-derived coronavirus MP789 (MT121216.1) with 82.9% nucleotide identity. Also of note was that no complete ORF10 was found in RsYN04, RmYN05, RmYN08, and a number of other SARS-CoV-2-related coronaviruses due to premature termination (Figure 2A), consistent with a previous study suggesting that ORF10 is not essential to the replication cycle of SARS-CoV-2 (Pancer et al., 2020).

Evolutionary history of sarbecoviruses

Phylogenetic analysis of full-length genome sequences of representative sarbecoviruses revealed that SARS-CoV-2 was most closely related to RaTG13, while RmYN02 and the Thailand strains formed a slightly more divergent clade. Notably, RpYN06 was placed at the basal position of the clade containing SARS-CoV-2 and its closest relatives from bats and pangolins (Figure 3 A, Table S5). In contrast, RsYN04, RmYN05, and RmYN08 grouped together and clustered with the pangolin-derived viruses from Guangxi, although separated from them by a relatively long branch. Finally, three SARS-CoV-related coronaviruses (RsYN03, RmYN07, and RsYN09) fell within the SARS-CoV lineage, grouping with other bat viruses previously sampled in Yunnan (Figure 3A).

A different topological pattern was observed in the phylogeny of the RdRp (Figure 3B). In particular, RpYN06 grouped with RmYN02 (although with weak bootstrap support), which together formed a clade with RaTG13, the two Cambodian strains, and SARS-CoV-2 (Figure 3B). The two bat-derived strains from Thailand formed a separate lineage. Perhaps more striking was that RsYN04, RmYN05, and RmYN08 now grouped with the Guangdong pangolin viruses (rather than those from Guangxi; Figure 3B). A different pattern again was observed in the phylogeny of the entire ORF1ab (Figure 3C). RpYN06 and RmYN02 now formed a clade and that was the direct sister-group to SARS-CoV-2, with RaTG13 a little more divergent (Figure 3C). In addition, RsYN04, RmYN05, and RmYN08 now clustered with the pangolin-derived strains from Guangxi (Figure 3C), consistent with the complete genome phylogeny.

In the spike gene phylogeny, SARS-CoV-2 and RaTG13 still grouped together, with both pangolin lineages falling as sister groups (Figure 3D). The two Cambodian bat viruses formed a separate and more divergent lineage. Strikingly, RpYN06 exhibited marked phylogenetic movement, this time clustering with two previously described bat viruses from Zhejiang province (ZC45 and ZXC21) whereas the Thailand bat virus clustered closely with RmYN02 (Figure 3D). In addition, RsYN04, RmYN05, and RmYN08 did not fall within the SARS-CoV and SARS-CoV-2 clades, but instead formed a separate and far more divergent lineage (Figure 3D). Finally, in the phylogeny of the RBD region, SARS-CoV-2 clustered with the pangolin viruses from Guangdong with the two Cambodian bat viruses the next most closely related viruses (Figure S2 A). RpYN06 fell within a lineage comprising several bat-derived betacoronaviruses, including ZC45, ZXC21, RsYN09, RsYN03, and RmYN07. As expected from the complete S gene tree, bat viruses RsYN04, RmYN05, and RmYN08 grouped together and formed a lineage characterized by a long branch (Figure S2A). A tanglegram of the representative sarbecoviruses clearly depicted the topological incongruence between the ORF1ab and spike gene phylogenies, particularly the SARS-CoV-2-related viruses, indicative of widespread recombination (Figure S2B).

Molecular characterization of the spike protein of the novel bat sarbecoviruses

At the six amino acid positions deemed critical for binding to the human angiotensin-converting enzyme 2 (hACE2) receptor, SARS-CoV-2 and the three bat-derived viruses identified here (RsYN04, RmYN05, and RmYN08) shared L455 and Y505. In contrast, despite being a closer overall relative, RpYN06 only possessed one identical amino acid with SARS-CoV-2: Y505 (Figure 4 A). At the S1/S2 cleavage site of the spike gene, none of the four SARS-CoV-2-related viruses reported here possessed a similar insertion/deletion (indel) pattern as SARS-CoV-2 (Garry et al., 2021) such that there was no furin cleavage site in these viruses (Figure 4A). Interestingly, however, the recently sampled bat virus from Thailand possessed a PVA three amino acid insertion at this site, similar to the PAA insertion found in RmYN02. The QTQTNS motif, proposed as the landing site for furin and other proteases, was similarly not observed in the newly described viruses (Figure 4A). In addition, two indel events have been identified in the RBD of many bat-associated coronaviruses (Holmes et al., 2021), and RpYN06 was characterized by indel patterns identical to those of ZC45 and ZXC21 (Figure 4A). There were no indel events in SARS-CoV-2 and the pangolin-derived coronaviruses in the RBD, and RsYN04, RmYN05, and RmYN08 possessed one unique indel event different to those observed in other sarbecoviruses (Figure 4A). In a similar manner to other bat-derived coronaviruses, the four novel SARS-CoV-2-related viruses possessed several indel events in the N-terminal domain, while RsYN04, RmYN05, and RmYN08 were again characterized by a unique indel pattern (Figure S3 A). Notably, RpYN06, ZC45, ZXC21, and the Guangdong pangolin virus shared the same indel pattern, with RpYN06 exhibiting high amino acid identity to these viruses in the N-terminal domain (85.29% to 99.02% at the amino acid level; Figure S3B).

We predicted and compared the three-dimensional structures of RpYN06, RsYN04, and SARS-CoV-2 using homology modeling (Figures 4B-4D). In a similar manner to RmYN02 (Zhou et al., 2020a), the RBD of RpYN06 had two shorter loops than those observed in SARS-CoV-2, while RsYN04 had only one shorter loop (Figure 4D). In addition, near the S1/S2 cleavage sites, the conformational loop of RpYN06 and RsYN04 were different from those of SARS-CoV-2 (Figures 4B and 4C). Notably, RsYN04 exhibited greater amino acid identity (71.28%) and shared more structural similarity with the SARS-CoV-2 RBD than RpYN06 (63.08%). Importantly, the conformational variations caused by these amino acid substitutions and deletions were speculated to interfere with the binding of RpYN06 and RsYN04 RBD to hACE2 (Figure 4D). However, RsYN04 exhibited lower structural similarity with SARS-CoV-2 in the N-terminal domain (NTD) (36.32% amino acid identity; Figure 4C, black arrowheads) than RpYN06 (60.77% amino acid identity).

To determine to what extent the deletions in the RBD region might interfere with the binding of the RpYN06 and RsYN04 RBDs to hACE2, RBDs from SARS-CoV-2, RsYN04, and RpYN06 were prepared and analyzed using flow cytometry and surface plasmon resonance (SPR) assays. The SARS-CoV-2 RBD that readily binds to the BHK-21 cells expressing hACE2 was used as a positive control (Figures 4E, 4H, and S4A–S4C). The RsYN04 RBD caused the fluorescence shift to the BHK-21 cells with the expression of hACE2, whereas the RpYN06 RBD did not (Figures 4F, 4G, and S4A–S4C). The equilibrium dissociation constant (K _D) of SARS-CoV-2 RBD binding to hACE2 was calculated to be 11.37 ± 0.57 nM (Figures 4H and S4 D), similar to previous results (Wu et al., 2020b). Consistent with the results from flow cytometry, no detectable binding was observed between the RpYN06 RBD and hACE2, even when the concentration of the RpYN06 RBD was increased to 100 μM (Figures 4I and S4D). However, RsYN04 RBD associated with hACE2, with the K _D of 15.87 ± 3.12 μM (Figures 4J and S4D).

Phylogenetic analysis of novel bat alphacoronaviruses

As well as betacoronaviruses, we identified 17 novel bat alphacoronaviruses. Phylogenetic analyses of the full-length genomes (Figure 5 A), the RdRp genes (Figure 5B), and ORF1ab (Figure S2C) of these 17 viruses and background representatives were consistent, with all trees revealing that the viruses newly identified here fell within four established subgenera: Decacovirus (n = 12), Pedacovirus (n = 1), Myotacovirus (n = 1), and Rhinacovirus (n = 2) (Figure 5). Of note were MlYN15 and RsYN25 isolated from Myotis laniger and R. stheno bats and closely related to swine acute diarrhea syndrome coronavirus (SADS-CoV) (Figure 5; Zhou et al., 2018) sharing nucleotide identities 87.55%–87.61%. In addition, HlYN18, isolated from a Hipposideros larvatus bat, fell within the subgenus Pedacovirus and was close to the porcine epidemic diarrhea virus (PEDV) lineage (Figure 5). Notably, the virus CpYN11 (isolated from Chaerephon plicatus) clustered with WA3607 (GenBank: MK472070; isolated from a bat from Australia), which together might represent an unclassified subgenus (Figure 5). Finally, RsYN14, RmYN17, McYN19, and RmYN24, although isolated from different bat species and sequenced on different lanes, were almost identical (with nucleotide identity >99.98% to each other) and might represent a novel species of subgenus Decacovirus.

Although the phylogenetic trees of the spike gene (Figure S2D) and protein sequences (Figure S2E) were topologically similar to those of the full-length genome, RdRp, and ORF1ab, a number of notable differences were apparent and indicative of past recombination events. First, CpYN11 clustered with HKU8 rather than WA3607 in the spike gene tree where they formed a separate lineage. Second, the topology of the subgenus Decacovirus in the spike gene tree was different to those observed in other gene regions. Finally, the two viruses belonging to the subgenus Tegacovirus were placed into the subgenera Pedacovirus (GenBank: NC_028806) and a separate lineage (GenBank: DQ848678), respectively.

Ecological modeling of the distribution of Rhinolophus species in Asia

To better understand the ecology of bat coronaviruses, we modeled the distribution of 49 Rhinolophus species in Asia using the collated distribution data and several ecological measures (Figures 6 and S5 ). The models performed well with a mean area under curve (AUC) of 0.96 for training and 0.92 for testing, and all training AUCs were above 0.88. Continentality (reflecting the difference between continental and marine climates) was, on average, the most important factor, contributing an average of 14.91% (based on permutation importance), followed by temperature seasonality at 11.7% average contribution, mean diurnal temperature range at 5.69%, and annual potential evapotranspiration at 5.38%. Three additional ecological factors also contributed more than 5% on average: minimum precipitation at 5.25%, potential evapotranspiration seasonality at 5.17%, and Emberger’s pluviothermic quotient (a measure of climate type) at 5%. The next most important factor was the distance to bedrock (an indicator of potential caves and rock outcrops) at 4.46%. Thus, local climate, especially factors that influence diet availability across the year, is seemingly key to determining bat species distributions across the region.

Although we could not accurately model diversity for Indonesia because of limited recently available data and likely high endemism, mainland Southeast Asia was well mapped (Figures 6 and S5). Most of mainland Southeast Asia’s remaining tropical forests showed a high diversity of rhinolophid bats, with a maximum of 23 species estimated to exist concurrently (Figure 6A). Rhinolophid hotspots occurred in forests throughout much of mainland Southeast Asia, with the largest contiguous hotspots extending from South Laos and Vietnam to Southern China (Figure 6A). Hotspots were also identified in the Hengduan mountains, and some parts of northern Myanmar and Nagaland in India (Figure 6A).

Interestingly, R. affinis (Figure 6B) and R. pusillus (Figure 6C) were widely distributed in Southeast Asia and southern China, and most bat species shared hotspots in Cambodia and peninsula Thailand. Several rhinolophid species extended their ranges northward into southern China reflecting the presence of forest (R. affinis and R. pusillus), whereas the geographic range of R. malayanus only just reached southern China (Figures 6D–6F). Ecological drivers for these species unsurprisingly showed some differences. Specifically, R. affinis was also influenced by temperature seasonality (16.59%), followed by Emberger’s pluviothermic quotient and mean diurnal range (8.79 and 8.7%), while R. malayanus (a smaller species) was mainly influenced by annual potential evapotranspiration mean (33.79%) and seasonality (14.57%). R. pusillus was influenced by temperature seasonality (12.44%) and continentality (9%), and R. shameli was largely influenced by annual potential evapotranspiration seasonality (34.81%) followed by annual evapotranspiration (9.79%). Overall, these factors control the range limits and food availability for these bat species.

It should be noted that the ecological modeling identified several other rhinolophid species with wide geographic distributions: R. huanensis, R. lepidus, R. luctus, R. macrotis, R. marshalli, R. microglobosus, R. pearsoni, R. rouxii, R. stheno, R. thomasi, and R. yunnanensis (Figure S5). Notably, R. stheno was found to host both SARS-CoV-2 and SARS-CoV-like coronaviruses.

Limitations of the Study

This study presents the identification of four SARS-CoV-2 related coronaviruses in bats, including one virus displaying high sequence identity to SARS-CoV-2 in most genomic regions. However, the direct evolutionary progenitor of SARS-CoV-2 remains unclear, and our sampling only considers a small number of bat species from a restricted geographic region. In addition, we did not consider other potential mammalian hosts for SARS-CoV-2 related coronaviruses that may also shed light on virus origins.

Key resources table

Resource availability

Experimental model and subject details

A total of 23 different bat species were tested in this study (Table S1). Samples were collected between May 2019 and November 2020 from Mengla County, Yunnan Province in southern China (101.271563 E, 21.918897 N; 101.220091 E, 21.593202 N and 101.297471 E, 21.920934 N). The Xishuangbanna Tropical Botanical Garden has an ethics committee that provided permission for trapping and bat surveys within this study.

Method details

Quantification and statistical analysis