key: cord-0802645-b0j0e3mu
authors: Colombo, Federico; Corsiero, Elisa; Lewis, Myles J.; Pitzalis, Costantino
title: Generation of restriction endonucleases barcode map to trace SARS-CoV-2 origin and evolution
date: 2021-06-03
journal: Sci Rep
DOI: 10.1038/s41598-021-91264-6
sha: 1b651d24751383a93b03b1e7419aa21748a94481
doc_id: 802645
cord_uid: b0j0e3mu

Since the first report of SARS-CoV-2 in China in 2019, there has been a huge debate about the origin. In this work, using a different method we aimed to strengthen the observation that no evidence of genetic manipulation has been found by (1) detecting classical restriction site (RS) sequence in human SARS-CoV-2 genomes and (2) comparing them with other recombinant SARS-CoV-like virus created for experimental purposes. Finally, we propose a novel approach consisting in the generation of a restriction endonucleases site map of SARS-CoV-2 and other related coronavirus genomes to be used as a fingerprint to trace the virus evolution.

Since the first report of SARS-CoV-2 in China in 2019, there has been a huge debate about the origin. In this work, using a different method we aimed to strengthen the observation that no evidence of genetic manipulation has been found by (1) 

Coronaviruses (CoVs) goes into the family Coronaviridae causing symptoms primarily in the upper respiratory tracts which range from common cold to severe to fatal illnesses 1 . They have been associated with two major disease outbreaks, the severe acute respiratory syndrome (SARS-CoV, 2002 ) and the Middle East respiratory syndrome (MERS-CoV, 2012) 2 . In December 2019, a new coronavirus (SARS-CoV-2) started to cause viral pneumonia bringing to severe and fatal infection. Although SARS-CoV-2 belongs to the same lineage of CoVs that causes SARS, it is genetically different and it cluster apart exploiting phylogenetic trees 3 . Phylogenetic analysis demonstrated the highly similarity between human SARS-CoV-2 and the sequence isolated from the Bat-Cov-raTG13 4 (97.2% identity) and the Pangolin-SARS-CoV 5 (80% identity), particularly in the receptorbinding-domain (RBD) of the S protein, important to mediate binding to human-receptor-angiotensin-converting-enzyme-2 (hACE2) 6 . The World Health Organization declared a coronavirus disease 2019 (COVID- 19) pandemic in March 2020. Therefore, one of the major discussions around SARS-CoV-2 has been related to its origin with the assumption that SARS-CoV-2 could have been the result of genetic manipulations or spill-over from laboratories studying these viruses. In March 2020, Anderson and colleagues published a detailed analysis showing that SARS-CoV-2 does not derive from a laboratory construct 7 . Although other several coronavirus experts have discredited the hypothesis of a man-made coronavirus [8] [9] [10] [11] [12] , here we aim to present a different method based on the analysis of restriction site (RS) sequences in the genome of SARS-Cov-2 to reconstruct its origin and follow the new variants.

What restriction sites (RS) sequence of the viral genome can say: generation of a restriction endonucleases barcoding map. During the SARS-CoV epidemic outbreak in 2003, a method called reverse genetic to assemble a full-length cDNA of the SARS-CoV-Urbani strain, as a template for manipulation of the viral genome, was published to develop and test candidate vaccines and therapeutics 13 . This resulted in the so-called infectious clone icSARS-CoV containing atypical markers of the wild-type (WT) virus. In particular, several Bgl1 RSs were introduced into the icSARS-CoV cDNA, which can be recognized since mutation are included in the newly formed cDNA. Figure 1A shows the sequence alignment between the WT SARS-CoV-Urbani and the icSARS-CoV. We highlighted the sequence containing the Bgl1 RS used to produce icSARS-CoV.

The newly sequences introduced in the recombinant cDNA of SARS-CoV can be used as markers to follow possible virus laboratory spillage. We analysed natural sequences isolated from four different SARS-CoV (hCoV-19-Italy-Vr/hSARS-CoV-19-Wuhan/hCoV-19 Pangolin/Bat-Cov-raTG13) to look for Bgl1 RS 'marker' www.nature.com/scientificreports/ (GCCNNNN/NGGC). All the genomes did not contain these sites (Table S1 ), in particular, the hSARS-CoV-19-Wuhan and hCoV-19-Italy-Vr. By analysing the sequences of Bat-CoV-raTG13 and the hCoV-19-Pangolin, we observed that only one sequence from SARS-CoV-Urbani (GCC AGC GTGGT) was found in SARS-CoV-2Wu (Wuhan). This is expected since the first part of these two genomes show high similarities.

Another recombinant SARS-CoV was produced in 2007 which derived from fifteen passages of the SARS-CoV-Urbani in BALB/c mouse lungs, therefore it was named Mouse-Adapted (MA)-SARS-CoV 14 . The identity of MA-SARS-CoV compared with the original SARS-CoV-Urbani is 99.97% with only six distinct nucleotides, that cannot be used as markers of this recombinant virus since same mutations are naturally acquired by the WT-SARS-CoV, as demonstrated from the sequences of other isolated SARS-CoV 14 . Both the icSARS-CoV and the MA-SARS-CoV have become the most widely used recombinant viruses to study SARS-like viruses and no specific sequences were found in the human SARS-CoV-2Wu.

In 2008, a consensus sequence called Bat-SCoV (FJ211859) was generated starting from four Bat-SCoVs genomes HKU3-1 (DQ022305), HKU3-2 (DQ084200), HKU3-3 (DQ084199), and RP3 (DQ071615) 15 . The full-length Bat-SCoV infectious clone, generated with the method described by Yount et al. 13 15 ). We observed that these specific sequences were all absent (Fig. 1B, Figure S1 ).

Other markers to identify the origin of SARS-Cov-2. In 2008, Ren et al.showed that SARS-like coronavirus (SL-CoVs) from horseshoe bat, which has a high similarity to SARS-CoV, differed in the N-terminus of the spike protein and particularly in the receptor binding RBS region 16 . Therefore, SL-CoVs were not able to infect hACE2 expressing cells, but only chimeric viruses expressing the spike protein of the SARS-CoV were able to bind the hACE2 which is the functional receptor of SARS-CoV. The authors identified a specific region responsible for the virus entrance into hACE2-expressing cells consisting of a minimal region of less than 200 amino acids. Interestingly, this group showed that chimeric spike proteins, whereby different regions of the SARS-CoV BJ01 (BJ01-S) spike were substituted into the spike of the bat SL-CoV (Rp3), and able to bind the hACE receptor. We generated in silico two of this chimeric spike (CS) sequences (the CS 424-494 and the CS 45-608 ), and then performed a multiple alignments to check similarities between other spikes identified after 2008, including the Bat-Cov-raTG13, the hCoV-19-Pangolin and the human SARS-CoV-2. The similarities of these two chimeric spikes are limited in the RBD of the spike (Fig. 1C ) and in the polybasic cleavage site ( Figure S2 ). Thus, the recombinant spike as possible progenitors of the hSARS-CoV-2 spike sequence can be excluded.

Moreover, we performed a nucleotide blast sequence to find whether these recombinant spikes are found in the recently identified SARS-CoV-2 viruses. As shown in Figure S3 -S5, we observed that, despite high similarities, many gaps (intended as single base mutations) are present between WT viruses and the recombinant spikes.

The turning point arrived in 2013, when Xing-Yi and colleagues published an important paper showing that a WT bat SL-CoV was capable of using hACE2 as an entry receptor, dispelling the observation that no natural SL-SARS-CoV were able to use hACE2. Interestingly, the newly identified bat SL-CoV-WIV1 had high sequence similarity (99.9% identity) to two other identified WT bat coronaviruses, RsSHC014 and RS3367. This study suggested that direct bat-to-human infection is a possible scenario for some bat SL-CoVs. In 2015, Vineet et al. made a recombinant virus between the spike of the bat coronavirus SHC014 and the mouse-adapted SARS-CoV backbone 17 using the well establish reverse genetic approach 13 . According to this method, several Bgl1 RSs were included into the sequence (Table S2) . Moreover, the sequences between the newly mutant SARS-CoV have poor sequence similarity to human SARS-CoV-19-Italy-VR and the SARS-CoV-19-Wuhan (Fig. 1D ).

Unique restriction sequence sites: a novel approach to track the SARS-CoV-2 origin. Exploiting the RS sequences, which are approximately 6-8 base pairs of DNA, as specific markers, we propose an alternative way to trace the SARS-CoV-2 origin. This approach consists in the generation of a RS map of SARS-CoV-2 and the other four related coronavirus genomes. Using the Serial Cloner Restriction Enzyme Library, we generated the RS barcoding map based on the frequency of finding specific RS sequences in the genome. First, we generated a RS barcoding map which was used as genetic fingerprinting of the specific sequence analysed and which easily highlights sequence differences between the genomes. The pattern of the barcode's reconstruction demonstrated high similarity between the coronavirus isolated from the Bat-Cov-raTG13 and the Pangolin, sug- www.nature.com/scientificreports/ gesting a natural evolution and adaptation of the virus. Different sequences of HIV, SARS-CoV and MERS-CoV were used as controls ( Fig. 2A and Figure S6 ). From the full restriction enzyme barcoding map, we identified in the spike (S) gene a sequence of 300 bp that can be used as a barcode to identify the virus and differentiate from others ( Fig. 2B and S7 ). It is also important to note that this method allows small mutations between new variants to be appreciated. This approach is low-cost and does not require full sequencing of the virus genome and extended analyses conducted by bioinformaticians. Indeed, by using a standard PCR reaction to amplify the above mentioned 300 bp spike gene, or simply by using real-time PRC products from swab test, and subsequent sequencing of this region, it is possible to generate an RS barcode that will give us a low-cost system to follow viral mutation and trace it over subsequent years. This would allow more samples to be tracked, especially with the emergence of new variants (English, Brazilian, African, etc.) that may potentially be more infectious and other variants for which vaccines may be less effective 18 . Moreover, this approach can easily be used to discriminate between false negative and false positive which are the reasons for important additional socio-economic disruptions 19 .

Using the data to generate the full barcode map we performed principal component analysis (PCA) to determine whether the observed frequency of RSs is related to the hierarchical distance of the genomes analysed. The PCA plot shows a cluster formed by the Pangolin, the Bat-Cov-raTG13 and human SARS-CoV-2 (5 different sequences including the English, Brazilian, African variant and the isolated sequences in Italy and Wuhan.) (Fig. 2C) . Nearby we find the cluster formed by the bat SARS-CoV related and other 5 SARS-CoV sequences. The HIV genomes and the MERS genomes were used as control and clearly clustered apart showing a greater difference in sequence identity from SARS-CoV virus.

In addition, we focalised on informative RSs to perform hierarchical clustering on the heatmap using Pearson correlation as the distance metric ( Fig. 2D and S8 ). The heatmap confirms that human SARS-CoV-2 (all the variants) and Bat-Cov-raTG13 are closer than hCoV-19 Pangolin and MERS CoV.

Finally, the barcode map of the RSs confirmed the absence of unique sites suggesting that the SARS-CoV-19 is the product of a natural evolutionary process of single base insertions/deletions or/and recombination.

We then focused on the unique RS sequences used to modify the viral genome. In particular, we analysed shared sites between SARS-CoV-19Wu, Bat-Cov-raTG13 and Pangolin-SARS-CoV-19. Only six RS sequences were shared between these genomes and their location does not suggest their use. In the Venn diagram shown in Fig. 3A , there are 12-shared RS. However, they are only six if we consider that some of these enzymes recognize the same sequences. One example is the unique RS sequence recognized by Bsp68I, BtuMI, NruI, RruI found at 319 bp on the Bat-Cov-raTG13 and shifted at 334 bp on the SARS-CoV-19Wu and the Pangolin-SARS-CoV-19. This 15 bp shift is due to single base insertions (Fig. 3B) .

Another example is the unique RS sequence GAG CTC recognized by Ecl136II on the SARS-CoV-19 genome that is located at 15081 bp, while on the Bat-Cov-raTG13 genome we found two of these sequences, one at 15080 bp and the other one at 19768 bp. The latter, if it were to be the result of genetic engineering, would be predicted to produce a gap of 6 bp, while from the local alignment it is clear that a nucleotide substitution occurred from C to T forming the new site ( Figure S9A) .

Finally, the genomic location of these unique sites does not flank specific ORF. Indeed, engineered RSs are typically expected to be at the beginning and at the end of an ORF. Here, all the unique RSs are located inside the ORFs (Figure S9B ), thus not easily editable by conventional genetic engineering.

Here, we analysed the peer-review literature of the SARS-related viruses generated in the laboratory over the years used to study the evolution of Coronaviruses and to generate drugs for their treatment. We have demonstrated through the analysis of RS, that SARS-CoV-2 does not contain peculiar RS or other markers that suggest a manipulation deriving from the recombinant viruses known in the literature. Indeed, the use of RS remains today the simplest, fastest and safest way to modify and study recombinant DNA. However, it should be mentioned that other methods have been used for generating recombinant clones, such as the one called transformation-associated recombination (TAR) cloning 20 . This method has proved effective for engineering large viral DNA 21, 22 and has only recently been developed for viruses with large RNA genomes exploiting Saccharomyces cerevisiae [23] [24] [25] . Also, bacterial artificial chromosomes (BACs) were exploited to manipulate the transmissible gastroenteritis coronavirus (TGEV) PUR46-MAD virus 26 and although nowadays other genetic manipulation mechanisms are known that allow no traces to be left, such as the use of the Crispr-cas system 27,28 , these remain more disadvantageous because they require higher technical capacity and higher costs and times. Furthermore, according to our knowledge, in the literature, there are no reports of SARS-CoV-2 coronavirus modifications through these more sophisticated techniques yet.

Finally, we used RS as markers to build a barcode map that could uniquely identify a particular virus. Recently, other authors used mathematical algorithms to find barcodes that identified a particular viral genome and build phylogenetic trees 29,30 while Son and colleagues developed a simple method to discriminate SARS-CoV-2 from SARS-CoV using one-step RT-PCR followed by restriction fragment length polymorphism 31 . Although this method, based on only three restriction enzymes, succeeds in distinguishing between SARS-CoV2 and SARS-CoV, it may not be sensitive to discriminate any mutations outside that region within the virus genome and does not allow new variants to be tracked. Guan and colleagues developed an elegant study based on the identification of a barcode that allows the SARS-CoV-2 genome to be assigned to 5 different clades, however this process requires advanced bioinformatics skills 32 . Here instead we have shown that with our method it is sufficient to sequence a region of 300 bp to build a specific barcode to distinguish the genome of a virus, included the new variants (English, Brazilian and African) and to trace its evolution over time. This would allow us to have useful information quickly and economically during the classic tests performed on swabs. www.nature.com/scientificreports/ Alignments. The sequences were aligned using Serial Cloner, Blastn suite 33 , ClustalW 34 and Jalview 35 .

Generation of restriction enzyme barcode. The restriction enzyme map barcode of each genome was obtained using Serial Cloner library. Using this software each genome was analysed in order to obtain the frequency of each restriction site to occur in that genome. The total frequencies of all the restriction sites present in the library were used to generate the barcode map. The InteractiVenn 36 was used to make Venn diagram.

Genomic distance in bp between restriction enzyme sites. The genomic distance in bp between two or more restriction enzymes sites was calculated with serial cloner and then reported graphically using Prism GraphPad v8.

PCA analyses was performed on the frequencies of the restriction enzymes sites on the different viruses' genomes and plotted by ggbiplot R-studio. Codes available upon request to the authors.

Heatmap and hierarchical clustering. The heatmaps were generated in R-studio by using frequencies of the restriction enzymes sites on the different viruses' genomes. The hierarchical clustering was performed using Pearson correlation as distance metric and Ward D clustering algorithm. Codes available upon request to the authors. www.nature.com/scientificreports/ 300 bp specific region. The 300 bp region was determined analysing the area of major discrepancy (low identity) between genomes, in particular between related genomes. This area was identified inside the spike (S) region. To generate the barcode map of these 300 bp regions we used the same method used for the full-length genomes. Thus, we calculated the frequencies of the restriction sites to generate the heatmap.

Informative sites. As informative sites we chose all those restriction sites that showed strong discrepancy in the cut-off frequency between the various genomes. Thus, to give an example, sites that had a high cut-off frequency in genome A compared to genome B, or sites unique to genome B that are repeatedly frequent in genome A (and vice versa). Then all non-informative sites, designated as those sites equally frequent across genomes, were discarded. In total we selected 104 informative sites here listed: "AatII" "AccBSI" "AcyI" "AfeI" "AloI" "Aor51HI" "AspA2I" "AsuNHI" "AvrII" "AxyI" "BarI" "BbvCI" "BcgI" "BlnI" "BmtI" "BplI" "BsaHI" "Bse21I" "BseYI" "BsiWI" "Bsp19I" "BspOI" "BsrBI" "BssNI" "BstACI" "Bsu36I" "BtgZI" "CchIII" "Cfr9I" "Ecl136II" "Eco32I" "Eco47III" "Eco53kI" "Eco81I" "EcoICRI" "EcoRV" "GdiII" "Hin1I" "Hsp92I" "MbiI" "MreI" "NcoI" "NheI" "NmeAIII" "Pfl23II" "PfoI" "Psp124BI" "PspLI" "PspOMII" "PsrI" "RpaBI" "SacI" "SauI" "SmaI" "SplI" "Sse232I" "SstI" "TspMI" "UcoMSI" "XmaI" "XmaJI" "ZraI" "AasI" "AccIII" "AgeI" "AsiGI" "BsePI" "BshTI" "Bsp13I" "BspEI" "BspMII" "BssHII" "CspAI" "DinI" "DrdI" "DseDI" "EciI" "Eco147I" "EgeI" "FspAI" "KasI" "Kpn2I" "KroI" "KspAI" "McaTI" "Mly113I" "MroI" "MroNI" "NaeI" "NarI" "NgoMIV" "PacI" "PasI" "PauI" "PceI" "PdiI" "PinAI" "PteI" "RceI" "SalI" "SfoI" "SseBI" "SspDI" "StuI". To generate the barcode map of the informative sites we used the method described for the full length and the 300 bp region.

The datasets and codes used and/or analysed during the current study are available from the corresponding author on reasonable request.

(0123456789) Scientific Reports | (2021) 11:11773 |

Coronavirus biology and replication: implications for SARS-CoV-2

The 2019-new coronavirus epidemic: evidence for virus evolution

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

A novel bat coronavirus closely related to SARS-CoV-2 contains natural insertions at the S1/S2 cleavage site of the spike protein

The proximal origin of SARS-CoV-2

Origin and cross-species transmission of bat coronaviruses in China

A genomic perspective on the origin and emergence of SARS-CoV-2

Evolutionary origins of the SARS-CoV-2 sarbecovirus lineage responsible for the COVID-19 pandemic

Probable pangolin origin of SARS-CoV-2 associated with the COVID-19 outbreak

The origin of SARS-CoV-2

Reverse genetics with a full-length infectious cDNA of severe acute respiratory syndrome coronavirus

A mouse-adapted SARS-coronavirus causes disease and mortality in BALB/c mice

Synthetic recombinant bat SARS-like coronavirus is infectious in cultured cells and in mice

Difference in receptor usage between severe acute respiratory syndrome (SARS) coronavirus and SARS-like coronavirus of bat origin

A SARS-like cluster of circulating bat coronaviruses shows potential for human emergence

An action plan for pan-European defence against new SARS-CoV-2 variants

False-positive COVID-19 results: hidden problems and costs

Selective isolation of genomic loci from complex genomes by transformation-associated recombination cloning in the yeast Saccharomyces cerevisiae

Genome-wide engineering of an infectious clone of herpes simplex virus type 1 using synthetic genomics assembly methods

Cloning, assembly, and modification of the primary human cytomegalovirus isolate toledo by yeast-based transformation-associated recombination

Rapid reconstruction of SARS-CoV-2 using a synthetic genomics platform

Infectious RNA transcripts from full-length dengue virus type 2 cDNA clones made in yeast

Rapid one-step construction of a Middle East Respiratory Syndrome (MERS-CoV) infectious clone system by homologous recombination

Engineering the largest RNA virus genome as an infectious bacterial artificial chromosome

Genome engineering using the CRISPR-Cas9 system

A programmable dual-RNA-guided DNA endonuclease in adaptive bacterial immunity

Genetic grouping of SARS-CoV-2 coronavirus sequences using informative subtype markers for pandemic spread visualization

Pitfalls of barcodes in the study of worldwide SARS-CoV-2 variation and phylodynamics

A simple method for detection of a novel coronavirus (SARS-CoV-2) using one-step RT-PCR followed by restriction fragment length polymorphism

A genetic barcode of SARS-CoV-2 for monitoring global distribution of different clades during the COVID-19 pandemic

Basic local alignment search tool

The EMBL-EBI search and sequence analysis tools APIs in 2019

Jalview Version 2-A multiple sequence alignment editor and analysis workbench

InteractiVenn: A web-based tool for the analysis of sets through Venn diagrams

F.C. designed the study, analysed data and wrote the manuscript. E.C. helped the study design, analysed the data, wrote and review the manuscript. M.J.L. gave contribution in the analyses of the data, contributed to the interpretation of the results and reviewed the manuscript. C.P. reviewed the manuscript.

The authors received no specific funding for this work.

The authors declare no competing interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-91264-6.Correspondence and requests for materials should be addressed to F.C. or E.C.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.