Emergence of genomic diversity and recurrent mutations in SARS-CoV-2.

SARS-CoV-2 is a SARS-like coronavirus of likely zoonotic origin first identified in December 2019 in Wuhan, the capital of China's Hubei province. The virus has since spread globally, resulting in the currently ongoing COVID-19 pandemic. The first whole genome sequence was published on January 5 2020, and thousands of genomes have been sequenced since this date. This resource allows unprecedented insights into the past demography of SARS-CoV-2 but also monitoring of how the virus is adapting to its novel human host, providing information to direct drug and vaccine design. We curated a dataset of 7666 public genome assemblies and analysed the emergence of genomic diversity over time. Our results are in line with previous estimates and point to all sequences sharing a common ancestor towards the end of 2019, supporting this as the period when SARS-CoV-2 jumped into its human host. Due to extensive transmission, the genetic diversity of the virus in several countries recapitulates a large fraction of its worldwide genetic diversity. We identify regions of the SARS-CoV-2 genome that have remained largely invariant to date, and others that have already accumulated diversity. By focusing on mutations which have emerged independently multiple times (homoplasies), we identify 198 filtered recurrent mutations in the SARS-CoV-2 genome. Nearly 80% of the recurrent mutations produced non-synonymous changes at the protein level, suggesting possible ongoing adaptation of SARS-CoV-2. Three sites in Orf1ab in the regions encoding Nsp6, Nsp11, Nsp13, and one in the Spike protein are characterised by a particularly large number of recurrent mutations (>15 events) which may signpost convergent evolution and are of particular interest in the context of adaptation of SARS-CoV-2 to the human host. We additionally provide an interactive user-friendly web-application to query the alignment of the 7666 SARS-CoV-2 genomes.

Introduction

On December 31 2019, China notified the World Health Organisation (WHO) about a cluster of pneumonia cases of uene">nkene">nowene">n aetiology iene">n Wuhaene">n, the capital of the Hubei Proviene">nce. The iene">nitial evideene">nce was suggestive of the outbreak beiene">ng associated with a seafood market iene">n Wuhaene">n, which was closed oene">n Jaene">nuary 1 2020. The aetiological ageene">nt was characterised as a SARS-like n class="Species">betacoronavirus, later named SARS-CoV-2, and the first whole genome sequence (Wuhan-HU-1) was deposited on NCBI Genbank on January 5 2020 (Wu et al., 2020). Human-to-human transmission was confirmed on January 14 2020, by which time SARS-CoV-2 had already spread to many countries throughout the world. Further extensive global transmission led to the WHO declaring COVID-19 as a pandemic on March 11 2020.

Coronaviridae comprise a large number of liene">neages that are fouene">nd iene">n a wide raene">nge of mammals aene">nd birds (Shaw et al., 2020), iene">ncludiene">ng the other n class="Species">human zoonotic pathogens SARS-CoV-1 and MERS-COV. The propensity of Betacoronaviridiae to undergo frequent host jumps supports SARS-CoV-2 also being of zoonotic origin. To date, the genetically closest-known lineage is found in horseshoe bats (BatCoV RaTG13) (Zhou et al., 2020). However, this lineage shares 96% identity with SARS-CoV-2, which is not sufficiently high to implicate it as the immediate ancestor of SARS-CoV-2. The zoonotic source of the virus remains unidentified at the date of writing (April 23 2020).

The analysis of genetic sequence data from pathogens is increasingly recognised as an important tool in infectious disease epidemiology (Rambaut et al., 2008; Greene">nfell et al., 2004). Geene">netic sequeene">nce data sheds light oene">n key epidemiological parameters such as doubliene">ng time of aene">n outbreak/epidemic, recoene">nstructioene">n of traene">nsmissioene">n routes aene">nd the ideene">ntificatioene">n of possible sources aene">nd aene">nimal reservoirs. Additioene">nally, whole-geene">nome sequeene">nce data caene">n iene">nform drug aene">nd vacciene">ne desigene">n. Iene">ndeed, geene">nomic data caene">n be used to ideene">ntify pathogeene">n geene">nes iene">nteractiene">ng with the host aene">nd allows characterisatioene">n of the more evolutioene">nary coene">nstraiene">ned regioene">ns of a pathogeene">n geene">nome, which should be prefereene">ntially targeted to avoid rapid drug aene">nd vacciene">ne escape mutaene">nts.

There are thousands of global SARS-CoV-2 whole-geene">nome sequeene">nces available oene">n the rapid data shariene">ng service hosted by the Global Iene">nitiative oene">n Shariene">ng All Iene">nflueene">nza Data (GISAID; https://www.epin class="Species">cov.org) (Elbe and Buckland-Merrett, 2017; Shu and McCauley, 2017). The extraordinary availability of genomic data during the COVID-19 pandemic has been made possible thanks to a tremendous effort by hundreds of researchers globally depositing SARS-CoV-2 assemblies (Table S1) and the proliferation of close to real time data visualisation and analysis tools including NextStrain (https://nextstrain.org) and CoV-GLUE (http://cov-glue.cvr.gla.ac.uk).

In this work we use this data to analyse the genomic diversity that has emerged in the global population of SARS-CoV-2 siene">nce the begiene">nniene">ng of the n class="Disease">COVID-19 pandemic, based on a download of 7710 assemblies. We focus in particular on mutations that have emerged independently multiple times (homoplasies) as these are likely candidates for ongoing adaptation of SARS-CoV-2 to its novel human host. After filtering, we characterise homoplasies at 198 sites in the SARS-CoV-2 genome. We identify a strong signal of recurrent mutation at nucleotide position 11,083 (Codon 3606 Orf1a), together with two further sites in Orf1ab encoding the non-structural proteins Nsp11 and Nsp13. These, together with a mutation in the Spike protein (21,575, Codon 5), comprise the strongest putative regions under selection in our dataset.

The current distribution of genomic diversity as well as ongoing allele frequency changes both between isolates and along the SARS-CoV-2 geene">nome are publicly available as aene">n opeene">n access aene">nd iene">nteractive web-resource available here:

https://macman123.shinyapps.io/ugi-scov2-aligene">nmeene">nt-screeene">n/.

Material and methods

Data acquisition

7710 SARS-CoV-2 assemblies flagged as “complete (>29,000 bp)”, “high n class="Species">coverage only”, “low coverage excl” were downloaded from the GISAID Initiative EpiCoV platform as of April 19 2020 (11:30 GMT). A full acknowledgements table of those labs which generated and uploaded data is provided in Table S1. Filtering was performed on the downloaded assemblies to exclude those deriving from animals (bat, pangolin), those with more than 1% missing sites, and otherwise spurious assemblies as also listed by nCov-GLUE (http://cov-glue.cvr.gla.ac.uk/#/excludedSeqs). This left a final dataset of 7666 assemblies for downstream analysis. Sequence metadata was obtained from the NextStrain Github repository (https://github.com/nextstrain/ncov/tree/master/data). While results presented here predominately focus on an analysis of the available assemblies as of April 19 2020, equivalent analyses were performed daily from March 24 2020. This allowed tracking of the emergence of genomic variants in public sequence data as assemblies were uploaded during the course of the pandemic.

Multi-sequence alignment and maximum likelihood tree

Assemblies were aligned against the Wuhan-Hu-1 refereene">nce geene">nome (n class="Gene">NC_045512.2, EPI_ISL_402125) using MAFFT (Katoh and Standley, 2013) implemented via the rapid phylodynamic alignment pipeline provided by Augur (https://github.com/nextstrain/augur). Sites in the first 130 bp and last 50 bp of the alignment were masked, as were positions 18,529, 29,849, 29,851 and 29,853, following the protocol also advocated by NextStrain and to account for the fact many putatively artefactual SNPs are located at the beginning and ends of the alignment. Resulting alignments were manually inspected in UGene (http://ugene.net). Subsequently a maximum likelihood phylogenetic tree was built using the Augur tree implementation selecting RAxML as the tree-building method (Kozlov et al., 2019). The resulting phylogeny was viewed and annotated using ggtree (Yu et al., 2017) (Figs. S1-S2). Throughout, site numbering and genome structure are given using Wuhan-Hu-1 (NC_045512.2) as reference.

Phylogenetic dating

The maximum likelihood phylogenetic tree was tested for the presence of significant molecular evolution over the sampling period using the roottotip() function provided in BactDating (Didelot et al., 2018). After confirmation of a significant regression following 1000 random permutations of sampling dates (Fig. S3), temporal calibration of the phylogeny was performed using TreeDater (Volz and Frost, 2017), assuming a strict clock model of evolution, as we do not expect a significant difference in rate variation across lineages at these time scales (Fig. S4). To obtain confidence intervals around each temporal point estimate we conducted a parametric bootstrapping analysis with 50 replicates on the unmasked alignment, keeping the tree topology constant while generating new branch length estimates using a Poisson distribution and running the same model in TreeDater (Volz and Frost, 2017). We also evaluated all currently available estimates for tip-calibration estimates of the tMRCA of SARS-CoV-2 together with rate estimates for other closely related n class="Species">betacoronaviruses (Table 1 , Table S2).

Table 1

Estimates of SARS-CoV-2 time to most recent common ancestor (tMRCA). BCI: Bayesian Credible Interval; HPD: Highest Posterior Density; CI: Confidence Interval. Asterix * denotes non-peer reviewed estimate of tMRCA. ‘N.’ denotes the number of whole genomes analysed.

Reference	N.	Substitution Rate (per site per year)	Estimated tMRCA	Method
Li et al. 2020 (Li et al., 2020)	32	1.0 × 10–3 (95% BCI 1.854 × 10–4, 4.0 × 10–3)	October 15, 2019 (95% BCI May 2, 2019; January 17, 2020)	Rate-informed strict clock model (BEAST v1.8.4)
Li et al. 2020 (Li et al., 2020)	32	1.8266 × 10–3 (95% BCI 7.5813 × 10–4, 3.0883 × 10–3)	December 6, 2019 (95% BCI November 16, 2019; December 21, 2019)	Rate-estimated relaxed clock model (BEAST v1.8.4)
Giovanetti et al. 2020 (Giovanetti et al., 2020)	54	6.58 × 10–3 (95% HPD 5.2 × 10–3, 8.1 × 10–3)	November 25, 2019 (95% CI September 28, 2019; December 21, 2019)	Relaxed clock model (BEAST v1.10.4)
Hill & Rambaut 2020*1	75	0.92 × 10–3 (95% HPD 0.33 × 10–3–1.46 × 10–3)	November 29, 2019 (95% CI October 28, 2019; December 20, 2019)	Unreported clock model (BEAST v1.7.0)
Hill & Rambaut 2020*1	86	0.80 × 10–3 (95% HPD 0.14 × 10–3, 1.31 × 10–3)	November 17, 2019 (95% CI August 27, 2019; December 19, 2019)	Unreported clock model (BEAST v1.7.0)
Hill & Rambaut 2020*1	116	1.04 × 10–3 (95% HPD 0.71 × 10–3, 1.40 × 10–3)	December 3, 2019 (95% CI November 16, 2019; December 17, 2019)	Unreported clock model (BEAST v1.7.0)
Lu et al. 2020* (41)	53	−	November 29, 2019 (95% HPD November 14, 2019; December 13, 2019)	Strict clock model (BEAST v1.10.0)
Duchene et al. 2020*2	47	1.23 × 10–4(95% HPD 5.63 × 10–4, 1.98 × 10–3)	November 19, 2019 (HPD October 21, 2019; December 11, 2019)	Strict clock model (BEAST v1.10)
Duchene et al. 2020*2	47	1.29 × 10–3 (HPD 5.35 × 10–4, 2.15 × 10–3)	November 12, 2019 (HPD September 26, 2019; December 11, 2019)	Relaxed clock model (BEAST v1.10)
Volz et al. 2020*3	53	Model constrained between 7 × 10–4 & 2 × 10–3	December 8, 2019 (95% CI November 21, 2019; December 20, 2019)	Strict clock model (BEAST v2.6.0)
Volz et al. 2020*3	53	Model constrained between 5 × 10–4 & 1.25 × 10–3	December 5, 2019(95% CI November 6, 2019; December 13, 2019)	Maximum Likelihood regression (treedater R package v0.5.0)

1http://virological.org/t/phylodynamic-analysis-of-sars-cov-2-update-2020-03-06/420; 2http://virological.org/t/temporal-signal-and-the-evolutionary-rate-of-2019-n-cov-using-47-genomes-collected-by-feb-01-2020/379; 3https://doi.org/10.25561/77169

Estimates of SARS-CoV-2 time to most receene">nt commoene">n aene">ncestor (tMRCA). BCI: Bayesiaene">n Credible Iene">nterval; n class="Disease">HPD: Highest Posterior Density; CI: Confidence Interval. Asterix * denotes non-peer reviewed estimate of tMRCA. ‘N.’ denotes the number of whole genomes analysed.

1http://virological.org/t/phylodyene">namic-aene">nalysis-of-n class="Species">sars-cov-2-update-2020-03-06/420; 2http://virological.org/t/temporal-signal-and-the-evolutionary-rate-of-2019-n-cov-using-47-genomes-collected-by-feb-01-2020/379; 3https://doi.org/10.25561/77169

Maximum parsimony tree and homoplasy screen

In parallel a Maximum Parsimony tree was built using the fast tree inference and bootstrap approximation offered by MPBoot (Hoang et al., 2018). MPBoot was run on the alignment to reconstruct the Maximum Parsimony tree and to assess branch support following 1000 replicates (−bb 1000). The resulting Maximum Parsimony treefile was used, together with the input alignment, to rapidly identify recurrent mutations (homoplasies) using HomoplasyFinder (Crispell et al., 2019).

HomoplasyFinder provides, for each site, the minimum number of state changes required on the tree to explain the observed character states at the tips, as described by Fitch (Fitch, 1971), and measured via the site specific consistency index. For this analysis all ambiguous sites in the alignment were set to ‘N’. To assess whether aene">ny particular Opeene">n Readiene">ng Frame (ORF) showed evideene">nce of more homoplasies thaene">n expected giveene">n the leene">ngth of the ORF, aene">n empirical distributioene">n was obtaiene">ned by sampliene">ng, with replacemeene">nt, equivaleene">nt leene">ngth wiene">ndows aene">nd recordiene">ng the number of homoplasies detected (Table S3).

HomoplasyFinder identified 1132 homoplasies (1042 excluding masked sites), which were distributed over the SARS-CoV-2 geene">nome (Fig. S5, Table S4). Of these, 40 sites have a derived allele at >1% of the total isolates. However, homoplasies caene">n arise due to coene">nvergeene">nt evolutioene">n (putatively adaptive), recombiene">natioene">n, or via n class="Disease">errors during the processing of sequence data. The latter is particularly problematic here due to the mix of technologies and methods employed by different contributing research groups. We therefore filtered identified homoplasies using a set of thresholds attempting to circumvent this problem (filtering scripts and figures are available at https://github.com/liampshaw/CoV-homoplasy-filtering).

In summary, for each homoplasy we computed the proportion of isolates with the homoplasy p where the nearest neighbouring isolate in the phylogeny also carried the homoplasy (excluding identical sequences). This metric ranges between p  = 0 (all isolates with the homoplasy present as singletons) and p  = 1 (no singletons i.e. clustering of isolates with the homoplasy in the phylogeny). We reasoned that artefactual sequencing homoplasies would tend to show up as singletons, so excluded all homoplasies with p  < 0.1 from further analysis.

To obtain a set of high confidence homoplasies, we then used the following criteria: ≥0.1% isolates in the alignment share the homoplasy (equivalent to >8 isolates), p  > 0.1, and derived allele found in strains sequenced from >1 originating lab and > 1 submitting lab. We also required the proportion of isolates where the homoplasic site was in close proximity to an ambiguous base (± 5 bp) to be zero. The application of these various filters reduced the number of homoplasies to 198 (Table S5). We also plotted the distributions of cophenetic distances between isolates carrying each homoplasy compared to the distribution for all isolates (Fig. S6), and inspected the distribution of all identified homoplasies in the phylogenies from our own analyses and on the phylogenetic visualisation platform provided by NextStraiene">n. Fiene">nally, we examiene">ned whether ambiguous bases were seeene">n more ofteene">n at homoplasic sites thaene">n at raene">ndom bases (excludiene">ng masked sites), which was not the case (Fig. n class="Gene">S7).

To further validate the homoplasy detection method applied to the alignment of the 7666 SARS-CoV-2 geene">nome assemblies, we took advaene">ntage of the geene">nome sequeene">nces for which raw reads were available oene">n the Short Read Arcn class="Disease">hive (SRA). A variant calling pipeline (available at https://github.com/DamienFr/CoV-homoplasy) was used to obtain high-confidence alignments for the 348 (out of 889 as of April 19 2020) SRA genomic datasets both meeting our quality criterions and matching GISAID assemblies. The topology of the Maximum Likelihood phylogeny of these 348 samples was compared to that of the corresponding samples from the GISAID genome assemblies using a Mantel test and the Phytools R package (Revell, 2012) (Figs. S8-S9, see Supplementary text).

As discussed, the GISAID dataset comprises assemblies of variable quality, potentially impairing the detection of genuine homoplasies and/or leading to false positive SNPs due to sequeene">nciene">ng n class="Disease">error or spurious allele assignment during the production of the de novo assembly from raw sequence reads. Therefore, to further assess the detection of homoplasies, we applied HomoplasyFinder to the two datasets comprising the same 348 strains (GISAID and SRA) (Table S6). We detected 19 homoplasies on the dataset originating from the SRA, and 21 on the dataset originating from GISAID assemblies. Of these, 19 were detected in both datasets (Table S7). Using the same filters as for the main dataset (with the exception of the ≥0.1% frequency set to ≥1%), 10 and 11 homoplasies were kept in the SRA dataset and in the GISAID dataset, respectively. Nine sites were detected in both datasets. For sites which failed the filtering thresholds, this was largely due to the low number of studied accessions, which increases the probability of an isolated strain displaying a homoplasy e.g. if n = 2 isolates have a homoplasy, by definition they cannot be nearest neighbours, so p  = 0.

Annotation of variant and homoplasic sites

The alignment was translated to amino acid sequences using SeaView V4 (Gouy et al., 2010). Sites were identified as synonymous or non-synonymous and amino acid changes corresponding to these mutations were retrieved via multiple sequence alignment. We assessed the change in hydrophobicity and charge of amino acid residues arising due to homoplastic non-synonymous mutations using the hydrophobicity scale proposed by Janin (Janin, 1979). The ten most hydrophobic residues on this scale were considered hydrophobic and the rest as hydrophilic. In addition, amino acid residues were either classified as positively charged, negatively charged or neutral at pH 7. The charge of each residue can either increase, decrease or remain the same (neutral mutation) due to mutation (Fig. S10).

Comparison with SARS-CoV-1 and MERS-CoV

SARS-CoV-1 aene">nd n class="Species">MERS-CoV are both zoonotic pathogens related to SARS-CoV-2, which underwent a host jump into the human host previously. We investigated whether the major homoplasies we detect in SARS-CoV-2 affect sites which also underwent recurrent mutations in these related viruses as these adapted to their human host. All Coronaviridae assemblies were downloaded (NCBI TaxID:11118) on April 8 2020 and human associated MERS-CoV and SARS-CoV-1 assemblies extracted. This gave a total of 15 assemblies for SARS-CoV-1 and 255 assemblies for MERS-CoV. Following the same protocol (Augur align) as applied to SARS-CoV-2 assemblies, each species was aligned against the respective RefSeq reference genomes: NC_004718.3 for SARS-CoV-1 and NC_019843.3 for MERS-CoV. This produced alignments of 29,751 bp (187 SNPs) and 30,119 bp (1588 SNPs) respectively.

MPBoot (Hoang et al., 2018) was run on both sets of alignments to reconstruct the maximum parsimony tree and to assess branch support following 1000 replicates (−bb 1000). The resulting maximum parsimony treefiles were used, together with the input alignment, to rapidly identify homoplasies using HomoplasyFinder (Crispell et al., 2019). For SARS-CoV-1 we detected six homoplasies aene">nd for n class="Species">MERS-CoV we detected 350 homoplasies (pre-filtering) (Fig. S11-S12). The distribution of homoplasies was assessed relative to the Genbank annotation files and in the context of the high confidence homoplasies that we detect in SARS-CoV-2.

Results

Emergence of SARS-CoV-2 genomic diversity over time

The 7666 SARS-CoV-2 geene">nomes offer aene">n excelleene">nt geographical aene">nd temporal n class="Species">coverage of the COVID-19 pandemic (Fig. 1a-b). The genomic diversity of the 7666 SARS-CoV-2 genomes is represented as Maximum Likelihood phylogenies in a radial (Fig. 1c) and linear layout (Fig. S1-S2). There is a robust temporal signal in the data, captured by a statistically significant correlation between sampling dates and ‘root-to-tip’ distances for the 7666 SARS-CoV-2 (Fig. S3; R 2 = 0. 20, p < .001). Such positive association between sampling time and evolution is expected to arise in the presence of measurable evolution over the timeframe over which the genetic data was collected. Specifically, more recently sampled strains have accumulated additional mutations in their genome than older ones since their divergence from the Most Recent Common Ancestor (MRCA, root of the tree).

Fig. 1

Global sequencing efforts have contributed hugely to our understanding of the genomic diversity of SARS-CoV-2. a) Viral assemblies available from global regions as of 19/04/2020. b) Cumulative total of viral assemblies uploaded to GISAID included in our analysis. c) Radial Maximum Likelihood phylogeny for 7666 complete SARS-CoV-2 genomes. Colours represent continents where isolates were collected. Green: Asia; Red: Europe; Purple: North America; Orange: Oceania; Dark blue: South America according to metadata annotations available on NextStrain (https://github.com/nextstrain/ncov/tree/master/data). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Global sequencing efforts have contributed hugely to our understanding of the genomic diversity of SARS-CoV-2. a) Viral assemblies available from global regioene">ns as of 19/04/2020. b) Cumulative total of viral assemblies uploaded to GISAID iene">ncluded iene">n our aene">nalysis. c) Radial Maximum Likelihood phylogeene">ny for 7666 complete n class="Species">SARS-CoV-2 genomes. Colours represent continents where isolates were collected. Green: Asia; Red: Europe; Purple: North America; Orange: Oceania; Dark blue: South America according to metadata annotations available on NextStrain (https://github.com/nextstrain/ncov/tree/master/data). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The origin of the regression between sampling dates and ‘root-to-tip’ distances (Fig. S3) provides a cursory point estimate for the time to the MRCA (tMRCA) around late 2019. Using TreeDater (Volz and Frost, 2017), we observe an estimated tMRCA, which corresponds to the start of the COVID-19 epidemic, of 6 October 2019–11 December 2019 (95% CIs) (Fig. S4). These dates for the start of the epidemic are iene">n broad agreemeene">nt with previous estimates performed oene">n smaller subsets of the n class="Disease">COVID-19 genomic data using various computational methods (Table 1), though they should still be taken with some caution. Indeed, the sheer size of the dataset precludes the use of some of the more sophisticated inference methods available.

The SARS-CoV-2 global populatioene">n has accumulated oene">nly moderate geene">netic diversity at this stage of the n class="Disease">COVID-19 pandemic with an average pairwise difference of 9.6 SNPs between any two genomes, providing further support for a relatively recent common ancestor. We estimated a mutation rate underlying the global diversity of SARS-Cov-2 of ~6 × 10−4 nucleotides/genome/year (CI: 4 × 10−4–7 × 10−4) obtained following time calibration of the maximum likelihood phylogeny. This rate is largely unremarkable for an RNA virus (Domingo-Calap et al., 2018; Holmes et al., 2016), despite Coronaviridae having the unusual capacity amongst viruses of proofreading during nucleotide replication, thanks to the non-structural protein nsp14 exonuclease, which excises erroneous nucleotides inserted by their main RNA polymerase nsp12 (Snijder et al., 2003; Minskaia et al., 2006).

Everything is everywhere

Some of the major clades in the maximum likelihood phylogeny (Fig. 1c and Fig. S1) are formed predominantly by strains sampled from the same continent. However, this likely represents a temporal rather than a geographic signal. Indeed, the earliest available strains were collected in Asia, where the COVID-19 paene">ndemic started, followed by exteene">nsive geene">nome sequeene">nciene">ng efforts first iene">n Europe aene">nd theene">n iene">n the USA.

The SARS-CoV-2 geene">nomic diversity fouene">nd iene">n most couene">ntries (with sufficieene">nt sequeene">nces) esseene">ntially recapitulates the global diversity of n class="Disease">COVID-19 from the 7666-genome dataset. Fig. 2 highlights the proportion of the global genetic diversity found in the UK, the USA, Iceland and China. In the UK, the USA and Iceland, the majority of the global genetic diversity of SARS-CoV-2 is recapitulated, with representatives of all major clades present in each of the countries (Fig. 2A-C). The same is true for other countries such as Australia (Fig. S2a).

Fig. 2

Genomic diversity of SARSCoV-2 in the USA, UK, Iceland and China. Strains collected from all four countries are highlighted on the global phylogenetic tree. a) Strains collected in the USA shown in purple. b) Strains from the UK shown in red. c) Strains collected in Iceland shown in red. d) Strains collected in China shown in green. Regional colours match to the global phylogeny shown in Fig. 1c. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Genomic diversity of SARSCoV-2 iene">n the USA, UK, Icelaene">nd aene">nd Chiene">na. Straiene">ns collected from all four couene">ntries are highlighted oene">n the global phylogeene">netic tree. a) Straiene">ns collected iene">n the USA showene">n iene">n purple. b) Straiene">ns from the UK showene">n iene">n red. c) Straiene">ns collected iene">n Icelaene">nd showene">n iene">n red. d) Straiene">ns collected iene">n Chiene">na showene">n iene">n greeene">n. Regioene">nal colours match to the global phylogeene">ny showene">n iene">n Fig. 1c. (For iene">nterpretatioene">n of the refereene">nces to colour iene">n this figure legeene">nd, the reader is referred to the web versioene">n of this article.)

This genetic diversity of SARS-CoV-2 populatioene">ns circulatiene">ng iene">n differeene">nt couene">ntries poiene">nts to each of these local epidemics haviene">ng beeene">n seeded by a large number of iene">ndepeene">ndeene">nt iene">ntroductioene">ns of the virus. The maiene">n exceptioene">n to this patterene">n is Chiene">na, the source of the iene">nitial outbreak, where oene">nly a fractioene">n of the global diversity is preseene">nt (Fig. 2d). This is also to aene">n exteene">nt the case for Italy (Fig. S2b), which was aene">n early focus of the n class="Disease">COVID-19 pandemic. However, this global dataset includes only 35 SARS-CoV-2 genomes from Italy, so some of the genetic diversity of SARS-CoV-2 strains in circulation likely remains unsampled. The genomic diversity of the global SARS-CoV-2 population being recapitulated in multiple countries points to extensive worldwide transmission of COVID-19, likely from extremely early on in the pandemic.

Genetic diversity along the genome alignment and recurrent mutations (homoplasies)

The SARS-CoV-2 aligene">nmeene">nt caene">n be coene">nsidered as brokeene">n iene">nto a large two-part Opeene">n Readiene">ng Frame (ORF) eene">ncodiene">ng noene">n-structural proteiene">ns, four structure proteiene">ns: n class="Gene">spike (S), envelope (E), membrane (M) and nucleocapsid (N), and a set of small accessory factors (Fig. 3a). There is variation in genetic diversity across the alignment, with polymorphisms often found in neighbouring clusters (Fig. S5). A simple permutation resampling approach suggests that both Orf3a and N exhibit SNPs which fall in the 95th percentile of the empirical distribution (Table S3). However, not all of these sites can be confirmed as true variant positions, due to the lack of accompanying sequence read data. However, we closely inspected those sites that appear to have arisen multiple times following a maximum parsimony tree building step. We identified a large number of putative homoplasies (n = 1042 excluding masked regions), which were filtered to a high confidence cohort of 198 positions (see Methods).

Fig. 3

Inspection of a major homoplastic site in Orf1ab of SARS-CoV-2 genome (position 11,083). Panel A shows a colour-coded schematic of the SARS-CoV-2 genome annotated as per NC_045512.2 and a plot of all potential homoplastic sites in Orf1ab measured as minimal number of character-state changes on a Maximum Parsimony tree (see Methods). Exemplar homoplasy (denoted with *) has been shown on the radial ML phylogenetic tree in panel B. Panel C shows the distribution of cophenetic distances between isolates carrying the identified homoplasy (red) and the distribution for all isolates (grey), showing that isolates with the homoplasy tend to cluster in the phylogeny. Equivalent figures for other filtered homoplasies are generated as part of the filtering method (see Methods). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Inspection of a major homoplastic site in Orf1ab of n class="Species">SARS-CoV-2 genome (position 11,083). Panel A shows a colour-coded schematic of the SARS-CoV-2 genome annotated as per NC_045512.2 and a plot of all potential homoplastic sites in Orf1ab measured as minimal number of character-state changes on a Maximum Parsimony tree (see Methods). Exemplar homoplasy (denoted with *) has been shown on the radial ML phylogenetic tree in panel B. Panel C shows the distribution of cophenetic distances between isolates carrying the identified homoplasy (red) and the distribution for all isolates (grey), showing that isolates with the homoplasy tend to cluster in the phylogeny. Equivalent figures for other filtered homoplasies are generated as part of the filtering method (see Methods). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

These 198 positions in the SARS-CoV-2 geene">nome aligene">nmeene">nt (0.67% of all sites) were associated with 290 amiene">no acid chaene">nges across all 7666 geene">nomes. Of these amiene">no acid chaene">nges, 232 comprised noene">n-syene">noene">nymous aene">nd 58 comprised syene">noene">nymous mutatioene">ns. Two noene">n-syene">noene">nymous mutatioene">ns iene">nvolved the iene">ntroductioene">n or removal of stop codoene">ns were fouene">nd (*13402Y, *26152G). 53 of the remaiene">niene">ng 101 noene">n-syene">noene">nymous mutatioene">ns iene">nvolved neutral hydrophobicity chaene">nges (Fig. S10a). Iene">n additioene">n, 79 of the remaiene">niene">ng 101 noene">n-syene">noene">nymous mutatioene">ns iene">nvolved neutral chaene">nges (Fig. S10b). Both n class="Gene">Orf1ab and N had a four-fold higher frequency of hydrophilic → hydrophobic mutations than hydrophobic → hydrophilic mutations (Fig. S10). In addition, neutral hydrophobic changes were clearly favoured in the S protein. Lastly, 87 of the remaining 110 non-synonymous mutations involved neutral charge changes.

Amongst the strongest filtered homoplasic sites (>15 change points on the tree), three are found within Orf1ab (nucleotide positioene">ns 11,083, 13,402, 16,887) aene">nd S (21575). We exemplify the stroene">ngest sigene">nal aene">nd our approach usiene">ng positioene">n 11,083 iene">n Fig. 3 aene">nd provide a full list of homoplasic sites, both filtered aene">nd uene">nfiltered, iene">n Tables S4–5. The stroene">ngest hit iene">n terms of the iene">nferred miene">nimum number of chaene">nges required (Fig. 3b-c) at n class="Gene">Orf1ab (11,083, Codon 3606) falls over a region encoding the non-structural protein, Nsp6, and is also observed in our analyses of the SRA dataset (Table S7).

We note that some of the hits also overlap with positions identified as putatively under selection using other approaches (http://virological.org/t/selection-analysis-of-gisaid-sars-cov-2-data/448/3, accessed April 23 2020), with n class="Gene">Orf1ab consistently identified as a region comprising several candidates for non-neutral evolution. Orf1ab is an orthologous gene with other human-associated betacoronaviruses, in particular SARS-CoV-1 and MERS-CoV which both underwent host jumps into humans from likely bat reservoirs (Lau et al., 2005; Memish et al., 2013). We performed an equivalent analysis on human-associated virus assemblies available on the NCBI Virus platform. We identified six putative homoplasic sites within SARS-CoV-1, two occurring within the 3c-like proteinase just upstream of Nsp6 (10,384, 10,793) and a further two homoplasies within Orf1ab at Nsp9 and Nsp13 (Fig. S11). In addition, one homoplasy was identified in the spike protein and one in the membrane protein ORFs.

For MERS-CoV, multiple uene">nfiltered homoplasies were detected, coene">nsisteene">nt with previous observatioene">ns of high recombiene">natioene">n iene">n this species (Dudas aene">nd Rambaut, 2016), though oene">nly oene">ne iene">nvoked more thaene">n a miene">nimum number of 10 chaene">nges oene">n the maximum parsimoene">ny tree (Fig. S12). This correspoene">nded to a a further homoplasy ideene">ntified iene">n n class="Gene">Orf1ab Nsp6 (position 11,631). It is of note that this genomic region coincides with the strongest homoplasy in SARS-CoV-2 which also occurs in the Nsp6 encoding region of Orf1ab. Codon 3606 of Orf1ab shares a leucine residue in MERS-CoV and SARS-CoV-2, though a valine in SARS-CoV. The exact role of these and other homoplasic mutations in human associated betacoronaviruses represents an important area of future work, although it appears that the Orf1ab region may exhibit multiple putatively adapted variants across human betacoronavirus lineages.

The genome alignment of the 7666 SARS-CoV-2 geene">nomes caene">n be queried through aene">n opeene">n access, iene">nteractive web-applicatioene">n (https://macmaene">n123.shiene">nyapps.io/ugi-sn class="Species">cov2-alignment-screen/). It provides users with information on every SNP and homoplasy detected across our global SARS-CoV-2 alignment and allows visual inspection both within the sequence alignment and across the maximum likelihood tree phylogeny. Fig. 3 illustrates some of the functionalities of the web application using position 11083 in the alignment as an example. This particular homoplasy was observed 1126 times across the genomes and requires a minimum of 37 character-site changes to become congruent with the observed SARS-CoV-2 phylogeny (Fig. 3a and b).

Discussion

Pandemics have been affecting humanity for milleene">nnia (Balloux aene">nd vaene">n Dorp, 2017). Over the last ceene">ntury aloene">ne, several global epidemics have claimed millioene">ns of lives, iene">ncludiene">ng the 1957/58 iene">nflueene">nza A (n class="Species">H2N2) pandemic, the sixth (1899–1923) and seventh ‘El Tor’ cholera pandemic (1961–1975), as well as the HIV/AIDS pandemic (1981-today). COVID-19 acts as an unwelcome reminder of the major threat that infectious diseases represent in terms of deaths and disruption.

One positive aspect of the current situation, relative to previous pandemics, is the unprecedented availability of scientific and technological means to face COVID-19. Iene">n particular, the rapid developmeene">nt of drugs aene">nd vacciene">nes has already beguene">n. Moderene">n drug aene">nd vacciene">ne developmeene">nt are largely based oene">n geene">netic eene">ngiene">neeriene">ng aene">nd aene">n uene">nderstaene">ndiene">ng of host-pathogeene">n iene">nteractioene">ns at a molecular level. The mobilisatioene">n to address the n class="Disease">COVID-19 pandemic by scientists worldwide has been remarkable. This includes the feat of the global scientific community who has already produced and publicly shared well over 11,000 complete SARS-CoV-2 genome sequences at the time of writing (April 23 2020), which we have used here with gratitude. Further initatives in the United Kingdom (https://www.cogconsortium.uk/data/) have already to date produced over 10,000 genomes, some of which overlap with those already available on GISAID.

To put these numbers of SARS-CoV-2 geene">nomes iene">n coene">ntext, it is iene">nterestiene">ng to coene">nsider parallels with the 2009 n class="Species">H1N1pdm influenza pandemic, the first epidemic for which genetic sequence data was generated in near-real time (Fraser et al., 2009; Smith et al., 2009). The genetic data available at the time looks staggeringly small in comparison to the amount that has already been generated for SARS-CoV-2 during the early stages of the COVID-19 pandemic. For example, Fraser et al. considered 11 partial hemagglutinin gene sequences two months after the WHO had declared 2009 H1N1pdm influenza a pandemic (Fraser et al., 2009).

This unprecedented genomic resource has already provided strong conclusions about the pandemic. For example, analyses by multiple independent groups place the start of the COVID-19 paene">ndemic towards the eene">nd of 2019 (Table 1). This rules out aene">ny sceene">nario that assumes n class="Species">SARS-CoV-2 may have been in circulation long before it was identified, and hence have already infected large proportions of the population.

Extensive genomic resources for SARS-CoV-1 should iene">n priene">nciple also be key to iene">nformiene">ng oene">n optimal drug aene">nd vacciene">ne desigene">n, particularly wheene">n coupled with kene">nowledge of n class="Species">human proteome and immune interactions (Gordon et al., 2020). Ideally, drugs and vaccines should target relatively invariant, strongly constrained regions of the SARS-CoV-2 genome, to avoid drug resistance and vaccine evasion. Therefore ongoing monitoring of genomic changes in the virus will be essential to gain a better understanding of fundamental host-pathogen interactions that can inform drug and vaccine design.

As most (but not all) pathogens capable of causing epidemic at a pandemic scale, SARS-CoV-2 is iene">n all likelihood of n class="Disease">zoonotic origin. This implies that SARS-CoV-2 may not be fine-tuned to its novel human host. However, it is near-impossible to predict future trajectories for the virulence and transmissibility of horizontally transmitted pathogens (Anderson and May, 1991). It is also possible that the population of SARS-CoV-2 will evolve into different lineages characterised by variable levels of virulence and transmissibility. However, despite existing phylogenetic structure (Rambaut et al., 2020), it is important to stress that there is no evidence for the evolution of distinct phenotypes in SARS-CoV-2 at this stage.

The vast majority of mutations observed so far in SARS-CoV-2 circulatiene">ng iene">n n class="Species">humans are likely neutral (Cagliani et al., 2020; Dearlove et al., 2020) or even deleterious (Nielsen et al., 2020). Homoplasies, such as those we detect here, can arise by product of neutral evolution or as a result of ongoing selection. Of the 198 homoplasies we detect (after applying stringent filters), some proportion are very likely genuine targets of positive selection which signpost to ongoing adaptation of SARS-CoV-2 to its new human host. Indeed, we do observe an enrichment for non-synonymous changes (80%) in our filtered sites. As such, our provided list (Table S5) contains candidates for mutations which may affect the phenotype of SARS-CoV-2 and virus-host interactions and which require ongoing monitoring. Conversely, the finding that 78% of the homoplasic mutations involve no polarity change could still reflect strong evolutionary constraints at these positions (Hughes, 2007; Yampolsky et al., 2005). The remaining non-neutral changes to amino acid properties at homoplasic sites may be enriched in candidates for functionally relevant adaptation and could warrant further experimental investigation.

One of the strongest homoplasies lies at site 11,083 in the SARS-CoV-2 geene">nome iene">n a regioene">n of n class="Gene">Orf1a encoding Nsp6. This site passed our stringent filtering cirteria and was also present in our analysis of the SRA dataset (Table S7). Interestingly, this region overlaps a putative immunogenic peptide predicted to result in both CD4+ and CD8+ T-cell reactivity (Grifoni et al., 2020). More minor homoplasies amongst our top candidates, identified within Orf3a (Table S5), also map to a predicted CD4 T cell epitope. While the immune response to SARS-CoV-2 is poorly understood at this point, key roles for CD4 T cells, which activate B cells for antibody production, and cytotoxic CD8 T cells, which kill virus-infected cells, are known to be important in mediating clearance in respiratory viral infections (Kohlmeier and Woodland, 2009). Of note, we also identify a strong recurrent mutation in nucleotide position 21,575, corresponding to the SARS-CoV-2 spike protein (codon 5). While the spike protein is the known mediator of host-cell entry, our detected homoplasy falls outside of the N-terminal and receptor binding domains.

Our analyses presented here provide a snapshot in time of a rapidly changing situation based on available data. Although we have attempted to filter out homoplasies caused by sequencing error with striene">ngeene">nt thresholds, aene">nd also used available short-read data to validate a subset of homoplasic sites iene">n a smaller dataset, our aene">nalysis nevertheless remaiene">ns reliaene">nt oene">n the uene">nderlyiene">ng quality of the publicly available assemblies. As such, it is possible that some results might be artefactual, aene">nd further iene">nvestigatioene">n will be warraene">nted as additioene">nal raw sequeene">nciene">ng data becomes available.

However, given the crucial importance of identifying potential signatures of adaptation in SARS-CoV-2 for guidiene">ng oene">ngoiene">ng developmeene">nt of vacciene">nes aene">nd treatmeene">nts, we have suggested what we believe to be a plausible approach aene">nd iene">nitial list iene">n order to facilitate future work aene">nd iene">nterpretatioene">n of the observed patterene">ns. More data coene">ntiene">nues to be made available, which will allow oene">ngoiene">ng iene">nvestigatioene">n by ourselves aene">nd others. We believe it is importaene">nt to coene">ntiene">nue to moene">nitor n class="Species">SARS-CoV-2 evolution in this way and to make the results available to the scientific community. In this context, we hope that the interactive web-application we provide will help identify key recurrent mutations in SARS-CoV-2 as they emerge and spread.

Author contributions

L.v.D., and F.B. conceived and designed the study; L.v.D., M.A, D.R L.P.S., C.E.F., L.O., C.J.O., J.P., C.C.S.T., F.A.T.B., and A.T.O analysed data and performed computational analyses; L.v.D., and F.B. wrote the paper with inputs from all co-authors.

Acknowledgments and funding

L.v.D and F.B. acknowledge financial support from the Newtoene">n Fuene">nd UK-Chiene">na n class="Gene">NSFC initiative (grant MR/P007597/1) and the (equipment grant BB/R01356X/1). Computational analyses were performed on UCL Computer Science cluster and the South Green bioinformatics platform hosted on the CIRAD HPC cluster. We thank Jaspal Puri for insights and assistance on the development of the alignment visualisation tool and Nicholas McGranahan and Rachel Rosenthal for their comments on the manuscript. We additionally wish to acknowledge the very large number of scientists in originating and submitting labs who have readily made available SARS-CoV-2 assemblies to the research community.

Declaration of Competing Interest

The authors have no competing interests to declare.