Origin and cross-species transmission of bat coronaviruses in China

Taxonomic and geographic sampling

We generated 630 partial sequences (440 nt) of the RNA-dependent RNA polymerase (RdRp) gene from bat rectal swabs collected in China and added 608 bat-CoV and 8 pangolin-CoV sequences from China available in GenBank or GISAID to our datasets (list of GenBank and GISAID accession numbers available in Supplementary Note 1). For each CoV genus, two datasets were created: one including all bat-CoV sequences with known host (host dataset) and one including all bat-CoV sequences with known sampling location at the province level (geographic dataset). To create a geographically discrete partitioning scheme that was more ecologically relevant than administrative borders for our phylogeographic reconstructions, we defined six zoogeographic regions within China by clustering provinces with similar mammalian diversity using hierarchical clustering⁴⁵ (see “Methods”): South western region (SW), Northern region (NO), Central northern region (CN), Central region (CE), Southern region (SO), and Hainan island (HI) (Fig. 1 and Supplementary Fig. 1).

Our host datasets included 701 α-CoV sequences (353 new sequences, including 102 new SADSr-CoV sequences (Rhinacovirus)) from 41 bat species (14 genera, five families) and 528 β-CoV sequences (273 new sequences, including 97 new SARSr-CoV sequences (Sarbecovirus)) from 31 bat species (15 genera, four families) (Supplementary Table 1). Our geographic datasets included 677 α-CoV sequences from six zoogeographic regions (22 provinces) and 503 β-CoV sequences from five zoogeographic regions (21 provinces) (Fig. 1). As some regions or hosts were overrepresented in our datasets, we also created and ran our analyses using a more uniform subset of our sequence data that included ~30 randomly selected sequences per host family or region to mitigate sampling and surveillance intensity bias.

Ancestral hosts and cross-species transmission

We used a Bayesian discrete phylogeographic approach implemented in BEAST⁴⁶ to reconstruct the ancestral host of each node in the phylogenetic tree using bat host family as a discrete character state. The phylogenetic reconstructions for α-CoVs in China suggest an evolutionary origin within rhinolophid and vespertilionid bats (Fig. 2a). The first α-CoV lineage to diverge historically corresponds to the subgenus Rhinacovirus (L1), originating within rhinolophid bats, and includes sequences related to HKU2-CoV and SADS-CoV (Supplementary Fig. 2). Then, several lineages, labeled L2–L7, emerged from vespertilionid bats (Fig. 2a). The subgenus Decacovirus (L2) includes sequences mostly associated with the Rhinolophidae and Hipposideridae and related to HKU10-CoV (Supplementary Fig. 3), while the subgenera Myotacovirus (L3) and Pedacovirus (L5), as well as an unidentified lineage (L4), include CoVs mainly from vespertilionid bats and related to HKU6-CoV, HKU10-CoV, and 512-CoV (Supplementary Figs. 4 and 5). Finally, a well-supported node comprises the subgenera Nyctacovirus (L6) from vespertilionid bats and Minunacovirus (L7) from miniopterid bats, and includes HKU7-CoV, HKU8-CoV, 1A-CoV, and 1B-CoV (Supplementary Fig. 6). These seven α-CoV lineages are mostly associated with a single host family, but each also included several sequences identified from other bat families (Fig. 2a, Supplementary Figs. 2–6, and Supplementary Table 1), suggesting that frequent cross-species transmission events have occurred among bats. Ancestral host reconstructions based on the random data subset, to normalize sampling effort, gave very similar results, with rhinolophids and vespertilionids being the most likely ancestral hosts of most α-CoV lineages too (Supplementary Fig. 7A). However, the topology of the tree based on the random subset was slightly different as the lineage L5 was paraphyletic.

Chinese β-CoVs likely originated from vespertilionid and rhinolophid bats (Fig. 2b). The maximum clade credibility (MCC) tree was clearly structured into four main lineages: Merbecovirus (lineage C), including MERS-related (MERSr-) CoVs, HKU4-CoV, and HKU5-CoV, and strictly restricted to vespertilionid bats (Supplementary Fig. 8); Nobecovirus (lineage D), originating from pteropodid bats and corresponding to HKU9-CoV (Supplementary Fig. 9); Hibecovirus (lineage E), comprising sequences isolated in hipposiderid bats (Supplementary Fig. 10) and Sarbecovirus (lineage B), including sequences related to HKU3- and SARS-related (SARSr-) CoVs originating in rhinolophid bats (Supplementary Fig. 11). We show that SARS-CoV-2 forms a divergent clade within Sarbecovirus and is most closely related to viruses sampled from Rhinolophus malayanus and R. affinis and from Malayan pangolins (Manis javanica) (Fig. 3). Similar tree topology and ancestral host inference were obtained with the random subset (Supplementary Fig. 7B).

We used a Bayesian stochastic search variable selection (BSSVS) procedure⁴⁷ to identify viral host switches (transmission over evolutionary time) between bat families and genera that occurred along the branches of the MCC annotated tree and calculated Bayesian factor (BF) to estimate the significance of these switches (Fig. 4). We identified nine highly supported (BF > 10) inter-family host switches for α-CoVs and three for β-CoVs (Fig. 4a, b). These results are robust over a range of sample sizes, with seven of these nine switches for α-CoVs and the exact same three host switches for β-CoVs having strong BF support (BF > 10) when analyzing our random subset (Supplementary Tables 2 and 3). To quantify the magnitude of these host switches, we estimated the number of host-switching events (Markov jumps)^48,49 along the significant inter-family switches (Fig. 4c, d) and estimated the rate of inter-family host-switching events per unit of time for each CoV genus. The rate of inter-family host-switching events was five times higher in the evolutionary history of α-CoVs (0.010 host switches/unit time) than β-CoVs (0.002 host switches/unit time) in China. For α-CoVs, host-switching events from the Rhinolophidae and the Miniopteridae were greater than from other bat families, while rhinolophids were the highest donor family for β-CoVs. The Rhinolophidae and the Vespertilionidae for α-CoVs and the Hipposideridae for β-CoVs received the highest numbers of switching events (Fig. 4c, d). When using the random dataset, similar results were obtained for β-CoVs, while rhinolophids were the highest donor family for α-CoVs (Supplementary Tables 4 and 5).

At the genus level, we identified 20 highly supported inter-genus host switches for α-CoVs, 17 of them were also highly significant using the random subset (Fig. 5a and Supplementary Table 6). Sixteen highly supported inter-genus switches were identified for β-CoVs (Fig. 5b). Similar results were obtained for the random β-CoV subset (Supplementary Table 7). Most of the significant cross-genus CoV switches for α-CoVs, 15 of 20 (75%), were between genera in different bat families, while this proportion was only 6 of 16 (37.5%) for β-CoVs. The estimated rate of inter-genus host-switching events (Markov jumps) was similar for α-CoVs (0.014 host switches/unit time) and β-CoVs (0.014 host switches/unit time). For α-CoVs, Rhinolophus and Miniopterus were the greatest donor genera and Rhinolophus was the greatest receiver (Supplementary Table 8). For β-CoVs, Rousettus was the greatest donor and Eonycteris the greatest receiver genus (Supplementary Table 9).

CoV spatiotemporal dispersal in China

We used our Bayesian discrete phylogeographic model with zoogeographic regions as character states to reconstruct the spatiotemporal dynamics of CoV dispersal in China. Eleven and seven highly significant (BF > 10) dispersal routes within China were identified for α- and β-CoVs, respectively (Fig. 6). Seven and five of these dispersal routes, respectively, remained significant when using our random subsets (Supplementary Tables 10 and 11). The Rhinacovirus lineage (L1) that includes HKU2-CoV and SADS-CoV likely originated in the SO region, while all other α-CoV lineages historically arose in SW China and spread to other regions before several dispersal events from SO and NO in all directions (Fig. 6a and Supplementary Fig. 12). A roughly similar pattern of α-CoV dispersal was obtained using the random subset (Supplementary Tables 10 and 12).

The oldest inferred dispersal movements for β-CoVs occurred among the SO and SW regions (Fig. 6b). The SO region was the likely origin of Merbecovirus (lineage C, including HKU4-CoV and HKU5-CoV) and Sarbecovirus subgenera (lineage B, including HKU3-CoV and SARSr-CoV), while the Nobecovirus (lineage D, including HKU9-CoV) and Hibecovirus (lineage E) subgenera originated in SW China (Supplementary Fig. 12). Then, several dispersal movements likely originated from SO and CE (Fig. 6b). More recent southward dispersal from NO was observed. Similar spatiotemporal dispersal patterns were observed using the random subset of β-CoVs (Supplementary Tables 11 and 13).

The estimated rate of migration events per unit of time along these significant dispersal routes was more than two times higher for α-CoVs (0.026 host switches/unit time) than β-CoVs (0.011 host switches/unit time), and SO was the region involved in the greatest total number of migration events for both α- and β-CoVs. SO had the highest number of outbound and inbound migration events for α-CoVs (Fig. 6c and Supplementary Table 12). For β-CoVs, the highest number of outbound migration events was estimated to be from NO and SO, while SO and SW had the highest numbers of inbound migration events (Fig. 6d and Supplementary Table 13).

Phylogenetic diversity

In order to identify the hotspots of CoV phylogenetic diversity in China and evaluate phylogenetic clustering of CoVs, we calculated the mean phylogenetic distance (MPD) and the mean nearest taxon distance (MNTD) statistics⁵⁰ and their standardized effect size (SES).

We found significant and negative SES MPD values, indicating significant phylogenetic clustering, within all bat families and genera for both α- and β-CoVs, except within the Aselliscus and Tylonycteris for α-CoVs (Fig. 7a, b). Negative and mostly significant SES MNTD values, reflecting phylogenetic structure closer to the tips, were also observed within most bat families and genera for α- and β-CoVs, but we found nonsignificant positive SES MNTD value for vespertilionid bats, and particularly for those in the Pipistrellus genus, for β-CoVs (Fig. 7a, b). In general, we observed lower phylogenetic diversity for β-CoVs than α-CoVs within all bat families and most genera when looking at SES MPD, but the difference in the level of diversity between α- and β-CoVs is less important when looking at SES MNTD (Fig. 7). These results suggest stronger basal clustering (reflected by larger SES MPD values) for β-CoVs than α-CoVs, indicating stronger host structuring effect and phylogenetic conservatism for β-CoVs. Very similar results were obtained with the random subsets for both α- and β-CoVs (Supplementary Tables 14–21).

We found negative and mostly significant values of MPD and MNTD (Fig. 7c and Supplementary Tables 22–25) indicating significant phylogenetic clustering of CoV lineages in bat communities within the same zoogeographic region. However, SES MPD values for α-CoVs in SW were positive (significant for the random subset), indicating a greater evolutionary diversity of CoVs in that region than others (Fig. 7 and Supplementary Tables 22–25). We used a linear regression analysis to assess the relationship between CoV phylogenetic diversity and bat species richness in China and determine if bat richness is a significant predictor of bat-CoV diversity and evolution. α-CoV phylogenetic diversity (MPD) was not significantly correlated to total bat species richness or sampled bat species richness in zoogeographic regions or provinces (Supplementary Table 26). Nonsignificant correlations between bat species richness and β-CoV phylogenetic diversity were also observed at the zoogeographic region level (Supplementary Table 27). However, a significant correlation was observed between sampled bat species richness and β-CoV phylogenetic diversity at the province level (Supplementary Table 27). Similar results were obtained when using the random subsets (Supplementary Tables 26 and 27). These findings suggest that bat host diversity is not the main driver of CoV diversity in China and that other ecological or biogeographic factors may influence this diversity. We observed higher CoV diversity than expected in several southern or central provinces (Hainan, Guangxi, Hunan) given their underlying total or sampled bat diversity (Supplementary Figs. 13 and 14).

We also assessed patterns of CoV phylogenetic turnover/differentiation among Chinese zoogeographic regions and bat host families by measuring the inter-region and inter-host values of MPD (equivalent to a measure of phylogenetic β-diversity) and their SES. We found positive inter-family SES MPD values, except between Pteropodidae and Hipposideridae for α-CoVs and between Rhinolophidae and Hipposideridae for β-CoVs (Fig. 8a, b and Supplementary Tables 28 and 29), suggesting higher phylogenetic differentiation of CoVs among most bat families than among random communities. Our phylo-ordination based on inter-family MPD values indicated that α-CoVs from vespertilionids and miniopterids, and from hipposiderids and pteropodids, as well as β-CoVs from rhinolophids and hipposiderids, were phylogenetically closely related (Fig. 8a, b). We also observed strong phylogenetic turnover between α-CoV strains from rhinolophids and from miniopterids and all other bat families, and between β-CoV strains from vespertilionids and all other bat families (Supplementary Tables 28 and 29). Phylo-ordination among bat genera based on inter-genus MPD confirmed these results and indicated that CoV strains from genera belonging to the same bat family were mostly more closely related to each other than to genera from other families (Fig. 8c, d and Supplementary Tables 30 and 31).

We observed high and positive inter-region SES MPD values between SW/HI and all other regions, suggesting that these two regions host higher endemic diversity (Fig. 9 and Supplementary Tables 32 and 31). Negative inter-region SES MPD values suggested that the phylogenetic turnover among other regions was less important than expected among random communities. Our phylo-ordination among zoogeographic regions also reflected the high phylogenetic turnover and deep evolutionary distinctiveness of both α- and β-CoVs from SW and HI regions (Fig. 9 and Supplementary Tables 32 and 33). Similar results were obtained using the random subset (Supplementary Tables 32 and 33).

Mantel tests

Mantel tests revealed a positive and significant correlation between CoV genetic differentiation (F_ST) and geographic distance matrices, both with and without provinces, including fewer than four viral sequences, for α-CoVs (r = 0.25, p = 0.0097; r = 0.32, p = 0.0196; respectively) and β-CoVs (r = 0.22, p = 0.0095; r = 0.23, p = 0.0336; respectively). We also detected a positive and highly significant correlation between CoV genetic differentiation (F_ST) and their host phylogenetic distance matrices, both with and without genera, including fewer than four viral sequences, for β-CoVs (r = 0.41, p = 0; r = 0.39, p = 0.0012; respectively) but not for α-CoVs (r = −0.13, p = 0.8413; r = 0.02, p = 0.5019; respectively).

Bat sampling

Bat oral and rectal swabs and fecal pellets were collected from 2010 to 2015 in numerous Chinese provinces (Anhui, Beijing, Guangdong, Guangxi, Guizhou, Hainan, Henan, Hubei, Hunan, Jiangxi, Macau, Shanxi, Sichuan, Yunnan, and Zhejiang). Fecal pellets were collected from tarps placed below bat colonies. Bats were captured using mist nets at their roost site or feeding areas. Each captured bat was stored into a cotton bag, all sampling was non-lethal and bats were released at the site of capture immediately after sample collection. A wing punch was also collected for barcoding purpose. Bat-handling methods were approved by Tufts University IACUC committee (proposal #G2017-32) and Wuhan Institute of Virology Chinese Academy of Sciences IACUC committee (proposal WIVA05201705). Samples were stored in viral transport medium at −80 °C directly after collection.

RNA extraction and PCR screening

RNA was extracted from 200 μl swab rectal samples or fecal pellets with the High Pure Viral RNA Kit (Roche) following the manufacturer’s instructions. RNA was eluted in 50 μl elution buffer and stored at −80 °C. A one-step hemi-nested reverse transcription-PCR (Invitrogen) was used to detect CoV RNA using a set of primers targeting a 440-nt fragment of the RdRp gene and optimized for bat-CoV detection (CoV-FWD3: GGTTGGGAYTAYCCHAARTGTGA; CoV-RVS3: CCATCATCASWYRAATCATCATA; CoV-FWD4/Bat: GAYTAYCCHAARTGTGAYAGAGC)⁸². For the first round PCR, the amplification was performed as follows: 50 °C for 30 min, 94 °C for 2 min, followed by 40 cycles consisting of 94 °C for 20 s, 50 °C for 30 s, 68 °C for 30 s, and a final extension step at 68 °C for 5 min. For the second round PCR, the amplification was performed as follows: 94 °C for 2 min, followed by 40 cycles consisting of 94 °C for 20 s, 59 °C for 30 s, 72 °C for 30 s, and a final extension step at 72 °C for 7 min. PCR products were gel purified and sequenced with an ABI Prism 3730 DNA analyzer (Applied Biosystems, USA). PCR products with low concentration or bad sequencing quality were cloned into pGEM-T Easy Vector (Promega) for sequencing. Positive results detected in bat genera that were not known to harbor a specific CoV lineage previously were repeated a second time (PCR + sequencing) as a confirmation. Species identifications from the field were also confirmed and re-confirmed by cytochrome (cytb) DNA barcoding using DNA extracted from the feces or swabs⁸³. Only viral detection and barcoding results confirmed at least twice were included in this study.

Sequence data

We also added bat-CoV RdRp sequences from China available in GenBank to our dataset. All sequences for which sampling year and host or sampling location information was available either in GenBank metadata or in the original publication were included (as of March 15, 2018). Our final datasets include 630 sequences generated for this study and 616 sequences from GenBank or GISAID (list of GenBank, China National Genomics Data Center and GISAID accession numbers available in Supplementary Note 1, and Supplementary Tables 34, 35, and 36). Nucleotide sequences were aligned using MUSCLE and trimmed to 360 base pair length to reduce the proportion of missing data in the alignments. All phylogenetic analyses were performed on both the complete data and random subset, and for α- and β-CoVs separately.

Defining zoogeographic regions in China

Hierachical clustering was used to define zoogeographic regions within China by clustering provinces with similar mammalian diversity⁴⁵. Hierarchical cluster analysis classifies several objects into small groups based on similarities between them. To do this, we created a presence/absence matrix of all extant terrestrial mammals present in China using data from the IUCN spatial database⁸⁴ and generated a cluster dendrogram using the function hclust with average method of the R package stats. Hong Kong and Macau were included within the neighboring Guangdong province. We then visually identified geographically contiguous clusters of provinces for which CoV sequences are available (Fig. 1 and Supplementary Fig. 1).

We identified six zoogeographic regions within China based on the similarity of the mammal community in these provinces: SW (Yunnan province), NO (Xizang, Gansu, Jilin, Anhui, Henan, Shandong, Shaanxi, Hebei, and Shanxi provinces and Beijing municipality), CN (Sichuan and Hubei provinces), CE (Guangxi, Guizhou, Hunan, Jiangxi, and Zhejiang provinces), SO (Guangdong and Fujian provinces, Hong Kong, Macau, and Taiwan), and HI. Hunan and Jiangxi, clustering with the SO provinces in our dendrogram, were included within the central region to create a geographically contiguous Central cluster (Supplementary Fig. 1). These six zoogeographic regions are very similar to the biogeographic regions traditionally recognized in China⁸⁵. The three β-CoV sequences from HI were included in the SO region to avoid creating a cluster with a very small number of sequences.

Model selection and phylogenetic analysis

Bayesian phylogenetic analysis was performed in BEAST 1.8.4⁴⁶. Sampling years were used as tip dates. Preliminary analysis were run to select the best-fitting combination of substitution models (HKY/GTR), codon partition scheme, molecular clock (strict/lognormal uncorrelated relaxed clock), and coalescent models (constant population size/exponential growth/GMRF Bayesian Skyride). Model combinations were compared and the best-fitting model was selected using a modified Akaike information criterion implemented in Tracer 1.6⁸⁶. We also used TEMPEST⁸⁷ to assess the temporal structure within our α- and β-CoV datasets. TEMPEST showed that both datasets did not contain sufficient temporal information to accurately estimate substitution rates or time to the most recent common ancestor. Therefore, we used a fixed substitution rate of 1.0 for all our BEAST analysis.

All subsequent BEAST analysis were performed under the best-fitting model, including an HKY substitution model with two codons partitions ((1 + 2), 3), a strict molecular clock and a constant population size coalescent model. Each analysis was run for 2.5 × 10⁸ generations, with sampling every 2 × 10⁴ steps. All BEAST computations were performed on the CIPRES Science Getaway Portal⁸⁸. Convergence of the chain was assessed in Tracer so that the effective sample size (ESS) of all parameters was >200 after removing at least 10% of the chain as burn-in.

Ancestral state reconstruction and transition rates

A Bayesian discrete phylogeographic approach implemented in BEAST 1.8.4 was used to reconstruct the ancestral state of each node in the phylogenetic tree for three discrete traits: host family, host genus, and zoogeographic region. An asymmetric trait substitution model was applied. These analyses were performed for each trait on the complete dataset and random subsets. MCC tree annotated with discrete traits were generated in TreeAnnotator and visualized using the software SpreaD3⁸⁹.

For each analysis, a BSSVS was applied to estimate the significance of pairwise switches between trait states using BF as a measure of statistical significance⁴⁷. BFs were computed in SpreaD3. BF support was interpreted according to Jeffreys in 1961⁹⁰ (BF > 3: substantial support, BF > 10: strong support, BF > 30: very strong support, BF > 100: decisive support) and only strongly supported transitions were presented in most figures, following a strategy used in other studies^91,92. We also estimated the count of state switching events (Markov jumps)^48,49 along the branches of the phylogenetic tree globally (for the three discrete traits) and for each strongly supported (BF > 10) transition between character states (for bat families and ecoregions only). Convergence of the MCMC runs was confirmed using Tracer. The rate of state switching events per unit of time was estimated for each CoV genus by dividing the total estimated number of state switching events by the total branch length of the MCC tree.

To assess the phylogenetic relationships among SARS-CoV-2 and other CoVs from the Sarbecovirus subgenus, we also reconstructed an MCC tree in BEAST 1.8.4 and median‐joining network in Network 10.0⁹³, including all Sarbecovirus sequences, two sequences of SARS-CoV-2 isolated in humans (GenBank accession numbers: MN908947 and MN975262), one sequence of SARS-CoV (GenBank accession number: NC_004718), eight sequences from Malayan pangolins (M. javanica) (GISAID accession numbers: EPI_ISL_410538-410544, EPI_ISL_410721) and one from Rhinolophus malayanus (GISAID accession number: EPI_ISL_412977) (Supplementary Note 1 and Supplementary Table 36).

Phylogenetic diversity

The MPD and the MNTD statistics⁵⁰ and their SES were calculated for each zoogeographic region, bat family, and genus using the R package picante⁹⁴. MPD measures the MPD among all pairs of CoVs within a host or a region. It reflects phylogenetic structuring across the whole phylogenetic tree and assesses the overall divergence of CoV lineages in a community. MNTD is the mean distance between each CoV and its nearest phylogenetic neighbor in a host or region, and therefore it reflects the phylogenetic structuring closer to the tips and shows how locally clustered taxa are. SES MPD and SES MNTD values correspond to the difference between the phylogenetic distances in the observed communities versus null communities. Low and negative SES values denote phylogenetic clustering, high and positive values indicate phylogenetic over-dispersion, while values close to 0 show random dispersion. The SES values were calculated by building null communities by randomly reshuffling tip labels 1000 times along the entire phylogeny. Phylogenetic diversity computations were performed on both the complete dataset and random subset for each trait. A linear regression analysis was performed in R to assess the correlation between CoV phylogenetic diversity (MPD) and bat species richness in China. Total species richness per province or region was estimated using data from the IUCN spatial database, while sampled species richness corresponds to the number of bat species sampled and tested for CoV per province or region in our datasets.

The inter-region and inter-host values of MPD (equivalent to phylogenetic β diversity), corresponding to the MPD among all pairs of CoVs from two distinct hosts or regions, and their SES were estimated using the function comdist of the R package phylocomr⁹⁵. The matrices of inter-region and inter-host MPD were used to cluster zoogeographic regions and bat hosts in a dendrogram according to their evolutionary similarity (phylo-ordination) using the function hclust with complete linkage method of the R package stats (R core team). These computations were performed on both the complete dataset and random subset.

Mantel tests and isolation by distance

Mantel tests performed in ARLEQUIN 3.5⁹⁶ were used to compare the matrix of viral genetic differentiation (F_ST) to matrices of host phylogenetic distance and geographic distance in order to evaluate the role of geographic isolation and host phylogeny in shaping CoV population structure. The correlation between these matrices was assessed using 10,000 permutations. To gain more resolution into the process of evolutionary diversification, these analyses were also performed at the host genus and province levels. To calculate phylogenetic distances among bat genera, we reconstructed a phylogenetic tree, including a single sequence for all bat species included in our dataset. Pairwise patristic distances among tips were computed using the function distTips in the R package adephylo⁹⁷. We then averaged all distances across genera to create a matrix of pairwise distances among bat genera. Pairwise Euclidian distances were measured between province centroids and log transformed. Mantel tests were performed with and without genera and provinces, including <4 viral sequences to assess the impact of low sample size on our results.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.