SARS-CoV-2 spillover transmission due to recombination event.

In late 2019, a novel Coronavirus emerged in China. Perceiving the modulating factors of cross-species virus transmission is critical to elucidate the nature of virus emergence. Using bioinformatics tools, we analyzed the mapping of the SARS-CoV-2 genome, modeling of protein structure, and analyze the evolutionary origin of SARS-CoV-2, as well as potential recombination events. Phylogenetic tree analysis shows that SARS-CoV-2 has the closest evolutionary relationship with Bat-SL-CoV-2 (RaTG13) at the scale of the complete virus genome, and less similarity to Pangolin-CoV. However, the Receptor Binding Domain (RBD) of SARS-CoV-2 is almost identical to Pangolin-CoV at the aa level, suggesting that spillover transmission probably occurred directly from pangolins, but not bats. Further recombination analysis revealed the pathway for spillover transmission from Bat-SL-CoV-2 and Pangolin-CoV. Here, we provide evidence for recombination event between Bat-SL-CoV-2 and Pangolin-CoV that resulted in the emergence of SARS-CoV-2. Nevertheless, the role of mutations should be noted as another influencing factor in the continuing evolution and resurgence of novel SARS-CoV-2 variants.

Introduction

Emerging infectious diseases, including Coronaviruses (CoVs), are often the result of a cross-species transmission of viruses from animals to humans. This can happen through several genetic mechanisms including recombination and mutations that give a virus new features and enable the virus to bind and enter into a new host cell with greater efficiency, avoid the immune system, and modifying its virulence (Longdon et al., 2014). An example of this cross species transmission is Severe Acute Respiratory Syndrome (SARS) (China-2002) and Middle East Respiratory Syndrome (MERS) (Saudi Arabia-2012). In the both scenarios, the CoV originated from bats, and amplified in a mammalian host, e.g. Himalayan palm civet (SARS) and dromedary camels (Chan et al., 2020a).

In December 2019, a novel coronavirus was reported in Wuhan, China. The virus was named SARS-CoV-2 and causes Coronavirus disease 2019 (COVID-19) with severe pneumonia symptoms (Guan et al., 2020). In the beginning of 2020, COVID-19 became a national disaster in China, and soon after swept the world. On 11 March 2020, the World Health Organization characterized COVID-19 as a pandemic (Millán-Oñate et al., 2020). In the late 2020, a novel variant of SARS-CoV-2 known as Variant Of Concern (VOC-202012-01) was emerged in England, and very soon became the dominant circulating viral variant in numerous countries around the world (Simmonds, 2020).

The pathogenicity of CoVs differs from each other (Cunha and Opal, 2014). Human CoVs mainly cause mild symptoms (e.g. 229E and OC43), but two Betacoronaviruses cause severe illness and fatality; (i) SARS-CoV infected 8000 humans (774 deaths) (Donnelly et al., 2003), (ii) MERS infected 2494 humans (858 deaths) (Lee et al., 2016). Like all other unknown viruses, very little is known about SARS-CoV-2.

The CoVs have several non-structural proteins (nsps), and structural proteins (sps). The most studied structural protein is the Spike protein (S-protein), which plays a key role in the pathogenicity of SARS-CoV. It is known that SARS-CoV uses human angiotensin-converting enzyme 2 (hACE2) as one of the main receptors (Lai et al., 2020). Mutations in the S-protein of animal SARS-CoVs might change protein structure, especially in RBD, and might make the virus compatible for binding to human cell receptors (Uddin et al., 2020).

Predicting the pathogenicity or transmissibility of a novel agent requires a detailed understanding of multiple factors. A possible tool for forecasting the function of a protein is 3D structure modeling of a protein from its sequences and then uses features of that structure to infer binding domains or other functional characteristics.

Previous studies have demonstrated that bats are the ultimate reservoir hosts for a number of CoVs, including ancestors of SARS-CoV and MERS-CoV (Hu et al., 2015). However, the evolutionary pathways of SARS-CoV-2 elusive. Genetic recombination in different viruses within a species is an important evolutionary process that cause genetic diversity and give new features to an ancestral virus. Studying recombination patterns provide evidence of origin, ecological links, and host range of a recombinant virus (Lam et al., 2013). The occurrence of homologous recombination events has been already detected in coronaviruses (Su et al., 2016). For example, it has been shown that recombination of CoVs in camels caused the emergence of dominant MERS lineage subsequently led to human outbreaks (Zhang et al., 2016).

The objective of a predictive oversight system is forecasting with a high degree of certainty the pathogenic potential of genomes of newly identified pathogens in comparison to the most closely related agent. Using sequence based prediction; we aim to study the mapping of the SARS-CoV-2 genome, modeling of protein architecture, determination of virulence drivers at the genomic level, as well as SARS-CoV-2 evolutionary origin, potential recombination events.

Material and methods

Dataset and gene mapping

A total number of 85 coronavirus sequences covering four genera, in which 35 were full-length sequences of SARS-CoV-2 from clades; L, V, S, G, GR (GISAID nomenclature) available in GenBank was retrieved. Genomic organization analysis, multiple alignments, open reading frames identification, and amino acid (aa) similarity and substitution were performed using Geneious version 11.0.5 (Biomatters Ltd., Auckland, New Zealand) (Shahhosseini et al., 2021).

3D modeling of S-protein

The homology model for SARS-CoV-2 S-protein and hACE2 were constructed using Oligomeric modeling, which combined interface conservation, structural clustering, and other template features to provide a quaternary structure quality estimate (QSQE) in the SWISS-MODEL workspace. The secondary structure was extracted from the 3D structure using DSSP program (Joosten et al., 2010; Shahhosseini et al., 2020).

Phylogenetic tree construction

In order to construct the phylogenetic tree, a Tamura-Nei genetic distances model and Neighbor-Joining (NJ) method were selected with sorted topologies and 70% threshold. Analyses of the sequences were conducted using Geneious software version 11.0.5. To further assess the precise evolutionary origin of SARS-CoV-2, a split network was made by EqualAngel method using the SplitsTree 4.14.8 software (Huson and Bryant, 2006; Shahhosseini et al., 2016).

Recombination analysis

To analyze the hypothesis whether the SARS-CoV-2 is a result of recombination in parental strains, we used Recombination Detection Program version 4 (RDP4) to detect and characterize recombination events in full-length sequences of SARS-CoV-2 (MN908947). The similarity plots and boot-scanning analyses were performed with Kimura 2-parameter, and window and moving step sizes of 200 and 20 nucleotides respectively. The phylogenetic tree of major and minor parents for recombinant sequence was inferred using UPGMA method in RDP4 software. The number of bootstrap replicates was adjusted on 100 and cutoff percentage on 70% (Chinikar et al., 2016; Lole et al., 1999).

Result

In-depth genome annotation of SARS-CoV-2

The SARS-CoV-2 genome, which encompasses 29,903 nucleotides with 12 open reading frames (ORFs), is compared with genetically relevant Betacoronaviruses in Table 1 . In the SARS-CoV-2 genome, the first two long ORF encompasses two ORFs: ORF1a and ORF1b (21,288 bp) together encodes 7096 aa, which produce non-structural proteins (nsps). The other ORFs encode structural proteins including the S gene (3822 bp) encodes 1273 aa for spike protein (S), the ORF3a gene (825 bp) encodes 275 aa for ORF3a protein, the E gene (225 bp) encodes 75 aa for envelope protein (E), the M gene (666 bp) encodes 222 aa for matrix protein (M), the ORF6 gene (183 bp) encodes 61 aa for ORF6 protein, the ORF7a gene (363 bp) encodes 121 aa for ORF7a protein, the ORF7b gene (129 bp) encodes 43 aa for ORF7b protein, the ORF8 gene (363 bp) encodes 121 aa for ORF8 protein, the N gene (1257 bp) encodes 419 aa for nucleocapsid protein (N), and the ORF10 gene (114 bp) encodes 38 aa for ORF10 protein (Table 1).

Table 1

Genomic organization, nucleotide and amino acid lengths of each ORF of the SARS-CoV-2 compared with genetically relevant Betacoronaviruses.

Virus species (nt length)/Gene Bank Acc. No.	Genomic organization/ Open reading frames
SARS-CoV-2(29,903 bp)/MN908947	ORF1ab266–21,555		S21,563–25,384		ORF325,393–26,220		E26,245–26,472		M26,523–27,191		ORF627,202–27,387		ORF7a27,394–27,759		ORF7b27,756–27,887		ORF827,894–28,259		N28,274–29,533		ORF1029,558–29,674
	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa
	21,288	7096	3822	1273	825	275	225	75	666	222	183	61	363	121	129	43	363	121	1257	419	114	38
Pangolin-CoV(27,213 bp)/MT084071a	ORF1ab266–21,555		S21,635–25,471		ORF325,479–26,310		E 26,334–26,562		M24,001–24,669		ORF627,291–27,477		ORF7a27,483–27,849		ORF7b27,845–27,977		ORF827,983–28,348		N28,363–29,628
	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa		aa
	21,288	7096	3837	1279	831	277	225	75	666	222	183	61	364	121	129	43	363	121	1257	419
Bat-SL-CoV-2(29,855 bp)/MN996532	ORF1ab251–21,537		S21,545–25,354		ORF325,363–26,190		E26,215–26,441		M26,493–27,158		ORF627,169–27,354		ORF7a27,360–27,725		ORF7b27,722–27,853		ORF827,860–28,225		N28,240–29,499
	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa
	21,285	7095	3807	1269	825	275	225	75	663	221	183	61	363	121	129	43	363	121	1257	419
Bat-SL -CoV(29,732 bp)/ MG772934	ORF1ab265–13,389		S21,483–25,220		ORF325,229–26,056		E26,081–26,308		M26,359–27,027		ORF727,038–27,223		ORF827,230–27,595		ORF10b27,730–28,095		N28,110–29,369		ORF1328,120–28,413		ORF1428,570–28,782
	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa
	21,210	7070	3735	1245	825	275	225	75	666	222	183	61	363	121	363	121	1257	419	291	97	210	70
SARS-CoV(29,727 bp)/AY291315	ORF1ab265–21,485		S21,492–25,259		ORF325,260–26,092		ORF3b25,689–26,153		E26,117–26,347		M26,398–27,063		ORF627,074–27,265		ORF727,273–27,641		ORF827,779–28,147		N28,120–29,388		ORF1328,130–28,426
	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa
	21,219	7073	3765	1255	822	274	462	154	228	76	663	221	189	63	366	122	366	122	1266	422	294	98

The submitted sequence of Pangolin-CoV (MT084071) to GenBank contains unread regions. In order to fill the gaps, a consensus sequence was generated from Pangolin-CoV metagenome, NCBI BioProject: PRJNA573298.

Identity of structural proteins

Pairwise identity analysis for ORF1a, ORF1b, S, ORF3, E, M, ORF6, ORF7a, ORF7b, ORF8, N and ORF10 are shown in Table 2 . Given structural proteins, the data showed that SARS-CoV-2 S-protein has highest similarity to Bat-SL-CoV-2 and Pangolin-CoV with 97.4% and 89.7%, respectively. The E-protein of SARS-CoV-2 is 100% identical with Bat-SL-CoV-2, Pangolin-CoV, and Bat-SL-CoV, and is 94.7% identical to SARS-CoV. The identity of SARS-CoV-2 M-protein is the same with Bat-SL-CoV-2, Pangolin-CoV and Bat-SL-CoV with 98.6%, but 89.1% with SARS-CoV. The N-protein of SARS-CoV-2 showed highest identity with Bat-SL-CoV-2 (99%) followed by Pangolin-CoV (97%), Bat-SL-CoV (94.2%), SARS-CoV (90%) (Table 2).

Table 2

Nucleotide and amino acid sequence identities among the genes of SARS-CoV-2 compared with Pangolin-CoV, Bat-SL-CoV-2, Bat-SL-CoV and SARS-CoV.

Virus strains	Gene regions (identity %)
	ORF1a		ORF1b		S		ORF3		E		M		ORF6		ORF7a		ORF7b		ORF8		N		ORF10
	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa	nt	aa
Bat-SL-CoV-2/MN996532	96	98	97.3	99.3	92.8	97.4	96.2	97.7	99.5	100	95.5	98.6	98.3	100	95.5	97	99.2	97.2	96.9	95	96.9	99	99.1	97.3
Pangolin-CoV/MT084071a	89.3	95.8	89.6	99.2	83	89.7	92.3	97	99.1	100	93.3	98.6	95.5	96.4	93.3	97.4	91.4	95.2	92.1	94.2	94.6	97	99.1	97.3
Bat-SL-CoV/MG772934	90.9	95.7	86.1	95.5	74.9	79.8	88.8	91.6	98.6	100	93.4	98.6	95	93.2	89.6	90	95.3	92.7	88.5	94.2	91.1	94.2	100	100
SARS-CoV/AY304488	76	81	86.2	95.6	72.5	75.9	75.5	96.9	93.5	94.7	84.9	89.1	76.5	68.2	84.1	89	86.1	83.8	40.9	16.1	88.0	90	93.1	82.4

The submitted sequence of Pangolin-CoV (MT084071) to GenBank contains unread regions. In order to fill the gaps, a consensus sequence was generated from Pangolin-CoV metagenome, NCBI BioProject: PRJNA573298.

Molecular structure of Spike protein

The S-protein of SARS-CoV-2 is composed of three identical units of polypeptide (homotrimer structure). To be able to anchor to host cells, the Spike has a RBD that contains 193 amino acids. Among 193 aa, 51 aa are in contact residues, in which the location of five aa known as key aa directly bind with hACE2, as shown in Fig. 1D. Also, the Spike protein possesses glycosylation sites, which has a key role in the penetration of the outer bilayer leaflet of the host cell membrane to initiate cell entry (Fig. 1).

Fig. 1

The aa sequence of the spike is colored based on protein structure (A). The predicted 3D structure of the SARS-CoV-2 spike on the surface of membrane protein. The surface of the spike model is colored based on protein structures (B). The representative shows spike from anchor side to host cells (C). The enlarged RBD with key residues binding to hACE2 (D).

Genetic determinants of virulence in S-proteins

Current sequence prognostication methods might predict the function of an individual protein based on its deduced aa sequences. In order to predict pathogenicity from sequences, we compared 51 aa within the contact residue of RBD in two SARS-CoV-2 variants and SARS-CoV, which showed evolutionary substitutions in 27 aa, in which five aa substitutions (Y455L, L486F, N493Q, D494S, and T501N/Y) occurred in key residues of RBD that are responsible for anchoring with hACE2. Notably, these five aa substitutions in key RBD residues of SARS-CoV-2 are completely identical to the same aa position in Pangolin-CoV (except at position 501 in the new variant of SARS-CoV-2 ‘VOC-202012-01’), but not Bat-SL-CoV-2. In addition, SARS-CoV-2 and Pangolin-CoV has a Glutamic acid insertion at position 484, which is unique and not seen in others compared aa sequences. Among key residues of RBD, only position 505 is strictly conserved in all compared aa sequences, except in Bat-SL-CoV-2 (Fig. 2A).

Fig. 2

The amino acids alignment in RBD (A) and O-linked glycan (B) of SARS-CoV-2 (MN908947), SARS-CoV-2 VOC-202012-01 (GISAID ID: EPI-ISL-601443), Pangolin-CoV (MT084071), Bat-SL-CoV-2 (MN996532), Bat-SL-CoV (MG772934), and SARS-CoV (AY394996) are shown respectively. Key residues in RBD are marked with black boxes (A). Three residues in O-linked glycan are marked with black boxes. The amino acid insertion in SARS-CoV-2 is marked with red box. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In O-linked glycan residues, three aa are critical for O-linked glycosylation attachment in SARS-CoV. Similarly in SARS-CoV-2, at two positions (S673 and S686), no aa substation was observed among aa alignments, but at position 678, Threonine in SARS-CoV2, Bat-SL-CoV-2 and Pangolin-CoV is substituted with Serine in Bat-SL-CoV and SARS-CoV. Notably, SARS-CoV-2 has four extra aa (P/HRRA) in O-linked residues, which is a unique feature in compared to other Sarbecoviruses. In addition, we observed that the acquisition of the polybasic cleavage site (RRAR) is restricted to SARS-CoV-2 at the junction (between aa position 685 and 686) of two subunits of the spike (Fig. 2B).

Phylogenetic tree construction based on full genes

The construction of cladogram revealed the grouping of coronaviruses into four distinct groups (Alphacoronavirus, Betacoronavirus, Deltacoronavirus, and Gammacoronavirus) (Fig. 3A). Further analysis showed that Betacoronaviruses have five distinct clades of viruses classified as Sarbecovirus, Merbecovirus, Embecovirus, Nobecovirus, and Hibecovirus, in which the SARS-CoV-2 fell within Sarbecovirus (Fig. 3B).

Fig. 3

The cladogram of coronaviruses (A) and the phylogenetic NJ tree of Betacoronaviruses (B) were constructed from aligned sequences using tools implemented within the Geneious software. The bootstrap values and number of bootstrap replications were greater than 70% and 1000, respectively. The neighbor-net network based on the full-length sequence of Sarbecoviruses was constructed using the Split tree software (C).

Split tree networking illustrated that the Sarbecovirus clade split into three subclades: Sarbeco-I isolates from bat species in Kenya and Bulgaria, Sarbeco-II divided to three host-specific groups including SARS-CoV isolates from humans, CoV isolates from civets and SARS-like CoV isolates from bats, whilst Sarbeco-III consist (i) SARS-CoV-2 isolated from humans, (ii) SARS-like CoVs isolated from bats, and Pangolin-CoV. All SARS-CoV-2 isolates tightly clustered together with more than 99% similarity (Fig. 3C).

Evidence for recombination

Similarity plotting analysis between SARS-CoV-2 and reference sequences showed recombination event in SARS-CoV-2 (Fig. 4 -A). The recombinant sequence (SARS-CoV-2) is shown as query, while sequences from the major (Bat-SL-CoV-2) and minor parents (Pangolin-CoV) are in purple and green, respectively. Further analysis with Bootscanning indicates recombination breakpoints beginning at nt 22,849 (99% CI: nt 22,711- nt 22,889) and ending at nt 23,093 (99% CI: nt 23,044- nt 23,145). Bat-SL-CoV-2 is the major parent with 98.4% similarity and Pangolin-CoV is the minor parent with 88.2% similarity for recombinant SARS-CoV-2. MaxChi matrix shows the statistically optimal positions of breakpoint pairs. The phylogenetic tree of alternative parents showed close similarity of SARS-CoV-2 (recombinant virus) with Bat-SL-CoV-2 (major parent) at the scale of the complete virus genome, except recombination point that shows highest evolutionary similarity to Pangolin-CoV (minor parent) (Fig. 4 B, C, D, E).

Fig. 4

Recombination analysis by RDP4 of SARS-CoV-2 as query sequence and close strains in Sarbecovirus are shown with Similarity plotting on a window size of 200 nucleotides and moving in steps of 20 nt along the alignment (A). Schematic sequence displays of recombinant region in green. The Bootscan analysis shows the most probable positions of breakpoint pairs, which are enlarged (B). Phylogenetically plausible alternative parents trees inferred in different genomic regions of recombinant sequence based on the major parent (C) and minor parent (D). MaxChi breakpoint matrix. Colors represent chi-squared values for breakpoint; dark red peaks indicate the most probable positions of breakpoint pair (E). The number of bootstrap replicates was adjusted on 100 and cutoff percentage on 70%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion

SARS-CoV-2 is a novel coronavirus and pathogenicity mechanisms are not yet clear. Given the evolution of pathogenicity in viruses, we used computational genomics to solve the controversy on origin and pathogenicity of SARS-CoV-2. Comparison of aa at S-protein showed that SARS-CoV-2 has the highest similarity to Pangolin-CoV than SARS-CoV. In addition, the critical residues binding to hACE2 are 100% identical in SARS-CoV-2 and Pangolin-CoV. This finding suggests that natural selection might happen in pangolin and gave the virus new features for recognition of hACE2. Our sequence-homology data is in accordance with a previous report discussing the low similarity between SARS-CoV S-protein and SARS-CoV-2 S-protein (Yuan et al., 2020). However, the conserved regions in RBD in two SARS-CoV and SARS-CoV-2 is sufficient to expect that SARS-CoV-2 can bind to hACE2 (Ralph et al., 2020).

The 3D modeling of S-protein showed the location of O-linked domain and RBD and its correspondent binding site at hACE2, implying the key role of S-protein in binding to host cells. Importantly, the overall homology of S-protein in SARS-CoV-2 and SARS-CoV are similar, suggesting use of the same cell entry receptor (Zhou et al., 2020). It is shown that at position 493, Gln is compatible for binding to salt bridge between Lys31 and Glu35 (hot spot 31) in hACE2. Also, the Leu455 and Phe486 make a favorable interaction with hot spot 31. The Asn501and Ser494 has favorable interactions with salt bridge between Lys353 and Asp38 (hot spot 353), but interactions are weaker than binding between the counterpart aa in SARS-CoV with hACE2 (Wan et al., 2020). In addition, a Glu484 insertion mutation changes the length of the loop, consequently altering aa positions in the RBD. The new variant of SARS-CoV-2 from England (VOC-202012-01) is defined with two deletions at positions 69–70 and 144. It is demonstrated that deletions (H69 and V70) in the unit S1 of S-protein is likely induce a conformational alteration comparing to other SARS-CoV-2 variants (Xie et al., 2021). The changes in homology in combination with other mutations in the RBD is highly likely enhancing the transmissibility of SARS-CoV-2 (Lauring and Hodcroft, 2021).

Both SARS-CoV-2 and SARS-CoV has 23 conserved aa at contact residue for binding to hACE2, supporting the idea that hACE2 is also the receptor of SARS-CoV-2. In addition, an in-vitro study to explore a treatment for COVID-19 showed that the production of human recombinant ACE2 can inhibit SARS-CoV-2 infection in early stages, supporting the idea that hACE2 is a receptor for SARS-CoV-2 (Monteil et al., 2020). To explore why SARS-CoV-2 has rapid human-to-human transmission, the binding affinity of SARS-CoV-2 RBD to hACE2 should be studied. The aa substitution in key residues of RBD of SARS-CoV-2 in comparison to SARS-CoV created stronger binding to hACE2, which is around 10–20 fold higher (Wrapp et al., 2020). Hence, the high binding affinity between SARS-CoV-2 RBD and hACE2 makes stronger interaction, which increase chances of virus attachment to host cell for virus entry, and promoting high infectivity. Moreover, it is assumed that decreased virulence of SARS-CoV-2 than SARS or MERS results in a longer time for tissue damage and onset of infection symptoms. It is estimated that the median incubation time for COVID-19 is 5 days (between 1 day to 14 days or longer) (Li et al., 2020), but apparently healthy patients are capable of shedding virus (Lauer et al., 2020). The combination of a long incubation time combined with virus shedding is another factor that contributed in the rapid spread of SARS-CoV-2 causing 2019–21 pandemic. In addition, among 13 non-synonymous mutations identified in SARS-CoV-2 VOC-202012-01, seven mutations (N501Y, A570D, D614G, P681H, T716I, S982A, D1118H) are in the S-gene (Rambaut et al., 2020). Of note, mutation N501Y, where amino acid Asparagine has been replaced with Tyrosine, in the key contact residue of RBD is demonstrated to increase infectivity and virulence in a mouse model (Gu et al., 2020).

The O-linked glycan domain in cleavage site of SARS-CoV-2 is identical to the SARS-CoV at positions Ser673 and Ser686, but position 678 in SARS-CoV-2 is Threonine while the same position in SARS-CoV is Serine. Linkage between Serine or Threonine in S-protein with a glycan has a critical role in viral entry into host cell (Bagdonaite and Wandall, 2018). In addition, we detected a unique aa (Proline) insertion (PRRA) at position 681 in O-linked glycan domain. However, new variant of SARS-CoV-2 (VOC-202012-01) has one non-synonymous mutation (P681H) at this domain. Hence, any changes at this site may play a role in altering pathogenicity capability of SARS-CoV-2. Moreover, it has been shown that introduction of aa into a furin cleavage site at the junction of S1 and S2 in SARS-CoV can mediate membrane fusion and virus infectivity (Andersen et al., 2020; Belouzard et al., 2009). Furthermore, the insertion of a polybasic cleavage site (RRAR) at the junction between unit S1 and S2 of spike can facilitate furin recognition and consequently allowing an efficient cleavage by host cell proteases, which play an important role in viral fusion into host cell membrane (Nao et al., 2017), which might be responsible for the efficient spread of SARS-CoV–2 (Coutard et al., 2020). Similarly, it has been shown that acquisition of a polybasic cleavage site is one of the necessary elements in altering a low-pathogenic avian influenza to a highly pathogenic strain (Stech et al., 2009).

The phylogenetic tree shows that SARS-CoV-2 has the closest evolutionary relationship with Bat-SL-CoV-2 at the scale of the complete virus genome, and less similarity to Pangolin-CoV. In this context, it can be assumed that the recombination between Bat- and Pangolin-CoVs happened deep in an ancestry of the bat lineage. However, the RBD of Bat-SL-CoV-2 showed low aa identity (78%) to SARS-CoV-2, suggesting that SARS-CoV-2 possibly was not transmitted directly from bats to humans. In contrast, the RBD of Pangolin-CoV is almost identical to SARS-CoV-2 (98%) at the aa level, proposing that mutations in Pangolin gave new features to CoV which enable the virus to bind to hACE2. This finding suggests the Pangolin as the intermediate host for SARS-CoV-2 (Li et al., 2020b). In addition, the role of adaptation pressures should be considered as an influencing factor that may have shaped RBD and mutations induced by host system apparently shaped SARS-CoV-2 genome including the critical amino acids within the RBD (Matyášek and Kovařík, 2020; Simmonds, 2020).

To clarify this contrast that SARS-CoV-2 has highest similarity to Bat-SL-CoV-2 when comparing the complete virus genome, but bats cannot be the direct origin of SARS-CoV-2, we investigated recombination events. Together Bootscanning curves and comparison of phylogenetic tree topology of major and minor parental viruses revealed recombination event between Bat-SL-CoV-2 and Pangolin-CoV with position of breakpoint pairs at RBD. Due to recombination event, the key residues of RBD from Pangolin-CoV integrated to Bat-SL-CoV-2, then zoonotic transmission occurred by a recombinant virus with capability to bind to hACE2. Also, we propose several mutations in the recombinant CoV leading to the emergence of different variants of SARS-CoV-2 in the human community with new features for rapid transmission and severity (Fig. 5 ).

Fig. 5

The proposed cycle for the emergence and resurgence of SARS-CoV-2 variants causing 2019–21 pandemic.

The recombination and phylogenetic analysis support the reliability of our model for SARS-CoV-2 origin and spillover transmission and further suggests that the previously proposed origins of SARS-CoV-2 which identify each bat or pangolin solely may not be accurate. For example, we observed several bias in the previous studies looking for the origin of SARS-CoV-2. The first bias was observed in those articles that identify bat as animal origin for SARS-CoV-2. Such studies argued that SARS-CoV-2 is very likely to be originated from Bat-SL-CoV-2 (RaTG13). However, they ignored Pangolin-CoV in the genomic analysis (Stech et al., 2009). Such approaches lead to this claim to exclude the role of recombination in emergence of SARS-CoV-2 (Stech et al., 2009).

The second bias was observed in articles that compared SARS-CoV-2 with Bat-SL-CoVs (MG772933 and MG772934), and ignored more genetically relevant Bat-SL-CoV-2 and Pangolin-CoV in analysis (Stech et al., 2009). Hence, they predicted that Bat-SL-CoV is the most probable source of SARS-CoV-2. Concurrently, authors could detect the occurrence of recombination events in the S-gene, but due to bias in dataset gathering they suggested recombination is more likely occurring in bat coronaviruses and is not the reason for emergence of SARS-CoV-2 (Chan et al., 2020a; Lv et al., 2020; Zhou et al., 2020). Whereas, another study suggested bat might be the original host of SARS-CoV-2 and an unknown animal might present an intermediate host (Paraskevis et al., 2020). Similarly, another study suggesting SARS-CoV-2 might be a recombinant virus between the Bat-SL-CoV and an origin-unknown CoV (Chan et al., 2020b; Malik et al., 2020; Xiong et al., 2020).

Since the initial stage of COVID-19 pandemic in early 2020, several SARS-CoV-2 variants have emerged in human host due to mutations (Volz et al., 2020). In the SARS-CoV-2 VOC, mutations in the S-gene, particularly N501Y in RBD, were postulated as a causative factor that gave the virus new features for rapid transmission. Consequently, the S-gene and its RBD domain are rapidly and frequently evolving in human hosts (Simmonds, 2020). Collectively, although recombination played a role in the emergence of SARS-CoV-2, convergence has a critical role in the recurrent emergence of novel SARS-CoV-2 variants in different phases of COVID-19 pandemic (Korber et al., 2020).

Conclusions

Here, we provide evidence for recombination event between Bat-SL-CoV-2 and Pangolin-CoV that resulted in the emergence of SARS-CoV-2. Our findings are in accordance with previous findings that viral recombination between different CoVs within animal populations may lead to the emergence of novel zoonotic CoVs that are lethal to humans (Wu et al., 2020).

CRediT authorship contribution statement

NSH, SC designed the study. NSH analyzed data. NSH, SC, GW, GK wrote the manuscript. NSH, GW revised the manuscript.

Declaration of competing interest

No conflict of interests to declare.