late December 2019 3 . Based on the high nucleotide identity of SARS-CoV-2 to a bat coronavirus isolate (96%) 7 , a possible scenario is that SARS-CoV-2 had undergone a period of adaptation in an as yet identified animal host, facilitating its capacity to jump species boundaries and infect humans 3 . The present rapid spread of the virus worldwide, coupled with its associated mortality, raises an important concern of its further potential to adapt to more highly transmissible or virulent forms.

The availability of SARS-CoV-2 genomic sequences concurrent with the present outbreak provides a valuable resource for improving our understanding of viral evolution across location and time. We performed an initial phylogenetic analysis of 749 SARS-CoV-2 genome sequences from late-December 2019 to March 13, 2020 (based on publicly available genome sequences on GISAID) and noted 152 SARS-CoV-2 sequences initially isolated in Europe beginning in February, 2020, which appear to have emerged as a distinct phylogenetic clade. Upon further investigation, we found that these strains are distinguished by a derived missense mutation in the spike protein (S-protein) encoding gene, resulting in an amino acid change from an aspartate to a glycine residue at position 614 (D614G). At the time of our study, the D614G mutation was at particularly high frequency 20/23 (87%) among Italian SARS-CoV-2 sequenced specimens, which was then emerging as the most severely affected country outside of China, with an overall case fatality rate of 7.2% 8 .

During the course of our analysis, the number of available SARS-CoV-2 genomes increased substantially. We subsequently analyzed a total of 2,795 genome sequences of SARS-CoV-2 (Supplementary data Table 1). For those sequences with demographics (65% of sequenced specimens), the male to female ratio was 0.56:0.44, with a median age of 49 years old, and a range from less than 1 to 99 years old. As of 30 March 2020, the D614G clade includes 954 of 1,449 (66%) of European specimens and 1,237 of 2,795 (44%) worldwide sequenced specimens. A comparison against the previous set of genomes collected for our phylogenetic and molecular dating analysis revealed that for samples submitted during the period from March 17-30, 2020, the D614G clade became increasingly prevalent worldwide, expanding from 22 to 42 countries (Fig. 1), as also reported previously 9,10 . The demographic distribution for this mutation, when known, (male to female ratio, 0.56:0.44; median age, 48 years old; age range, less than 1-99 years old) was not significantly different compared to the reported demographics for all sequenced SARS-CoV-2 specimens.

We employed molecular dating to estimate the time of emergence of the D614G clade. Based on a curated set of 442 genomes representing the sequence diversity of SARS-CoV-2 samples available at the time of analysis (30 March 2020), the mean time to most recent common ancestor (tMRCA) was estimated to be on 18 January 2020 (95% highest posterior density (HPD) interval: 10 January-25 January), indicating its relatively recent emergence (Fig. 2). Although the mutation appears to be clade-specific, we noted that D614G also arose another time in a single isolate belonging to a distinct lineage, Wuhan/HBCDC-HB-06/2020 (EPI 412982), collected on 7 February 2020. Given the high degree of nucleotide identity of the D614G clade (~ 99.6%), we expect that future tMRCA estimates will not differ substantially. The mean tMRCA for the root of the tree was estimated to be 13 November 2019 (95% HPD interval: 19 October 2019 to 5 December 2019) which is similar to what has been estimated previously 4 . Multiple factors can influence these estimates as more data become available: increased diversity and number of strains included for analysis, expanded range of sampling dates compared to earlier studies, as well as incomplete/biased sampling of available genomes. Therefore, caution should be used when inferring broad and definite conclusions about the epidemiology of the emerging outbreak in real-time, particularly as the early undocumented stages of SARS-CoV-2 transmission is largely unknown 3 . They should be taken as tentative estimates based on limited sampling, which are subject to change when additional epidemiological information becomes available. For instance, as different countries review and test archived specimens from cases of severe pneumonia or influenza-like illness for SARS-CoV-2, it is expected that additional cases may be identified, such as in France where a patient without travel history to China was identified to have COVID-19 in late December concurrent with the initial reported cases from Wuhan 11 . These retrospective analyses will provide crucial insights into the early transmission dynamics and evolution of SARS-CoV-2 and its rapid global spread.

From our findings of the recent emergence of the D614G clade and the increasing number of specimens harboring the mutation identified worldwide, we sought to investigate the potential significance of the mutation on clinical disease severity phenotypes. Unfortunately, only limited clinical disease severity data was available for patients whose sequenced samples were included in our analysis (asymptomatic 64, symptomatic 8, mild symptoms 2, pneumonia 4, hospitalized 163, released 48, recovered/recovering 26, nursing home 3, live 24, intensive care unit 2, deceased 1), preventing us from being able to meaningfully correlate disease severity and genotype from these data. However, using country-wide crude case fatality data for countries from which sequencing data was available, there was no significant correlation between proportion of D614G clade sequences and crude case fatality rate as of 30 March 2020 (Spearman's rank correlation coefficient, r 0.22, 95% CI -0.12 to 0.51). In addition, on analysis of crude case fatality rate by age-group (available for China, Italy, South Korea, Spain, and Canada) there was no significant correlation with proportion of D614G clades in the sequences analysed for these countries (Supplementary Table 2) 8,[12][13][14] . A major limitation of this analysis is that publicly available SARS-CoV-2 genomes are the result of convenience sampling and are not expected to provide an accurate representation of the spatial-demographic distribution of SARS-CoV-2 genotypes. Therefore, one must be cautious about making inferences about severity and transmission of variants 15 from genomic sequence data alone.

The possibility that the D614G mutation may still have a potential impact on the function of the SARS-CoV-2 spike protein could not be excluded. Previous studies of SARS-CoV have shown that the sequential accumulation of mutations in the spike protein increased its affinity to ACE2 and likely impacted its transmission and disease severity during the course of outbreaks in 2002-2004 [16][17][18] . Recently, Phan studied 86 complete SARS-CoV-2 genomes, identifying 42 missense mutations, eight of which occurred in the S-protein gene 2 . However, at the time of the study, the mutation D614G was only found in one sequence from Germany (collected on 28 January, 2020). Modifications in the spike protein are of interest as they might indicate the emergence of a novel strain of SARS-CoV-2 with change in transmissibility or pathogenicity. Therefore, we investigated the potential functional and epidemiological consequences of the D614G mutation with structural modeling of the SARS-CoV-2 spike (S) protein and its interaction with the angiotensin-converting enzyme 2 (ACE2) receptor.

The structure of the SARS-CoV-2 spike (S) protein is shown in Fig. 3. The S protein is a heavily glycosylated trimeric protein that mediates entry to host cells via fusion with ACE2. Recently, Wrapp and colleagues used Cryo-EM to determine the structure of the S protein and analyze its conformational changes during infection 19 . Using their three-dimensional model of the S protein structure, we set out to investigate the effects that a mutation in position 614 might have. First, we analyzed changes in inter-atomic contacts within a radius of 6 Å from position 614 before and after substitution of aspartate by glycine. Notably, four inter-chain destabilizing (i.e., hydrophobic-hydrophilic) contacts are lost with residues of an adjacent chain upon D614G mutation (see Table 1 and Fig. 3C). This suggests that a small repelling interaction between adjacent chains is removed upon this aspartate substitution (see Table 1). However, it is unlikely that this would have a significant effect on recognition and binding to ACE2 given the relative distal position of this mutation with respect to the receptor-binding domain (RBD) (see Fig. 3B), but further analyses would be required to assess whether the D614G mutation has an effect on the way the S protein changes its conformation after interaction with ACE2. Lan and colleagues also showed residues in the RBD act as epitopes for SARS-CoV-2 and mutations can influence antibody binding 20 . Given the important role that glycosylation plays in regulating the function of spike proteins in coronaviruses 21 , we decided to search for potential changes in a glycosylated residue (asparagine) in position 616. As shown in Fig. 3D, this residue and aspartate 614 are oriented in opposite directions, suggesting that substitution by glycine also has a null effect on this interaction. Finally, neither position 614 nor the inter-atomic contacts at positions 854, 859, 860, 861 of the spike protein lie in a polybasic cleavage region which is of importance for SARS-CoV-2 as it has been proposed to activate the protein for membrane fusion 22 . Our results are consistent with the work by Wrapp and colleagues, where they identified nine missense mutations (including D614G) in the spike protein that were thought to be relatively conservative and thus unlikely to affect protein function 19 .

There were important limitations faced in our present analysis which are likely to be a significant hurdle to similar studies in the future. The lack of available clinical metadata prevented our investigation of association between viral clade and disease severity phenotype. Additionally, numbers of sequenced SARS-CoV-2 samples vary greatly between countries and may be subject to potential sampling bias. It is important to note that current country level data on crude case fatality rates and case numbers do not permit robust comparison of clinical phenotype across countries due to significant differences in population demographics, testing protocols, case definitions and implementation of public health measures. Accordingly, it was not possible to draw any conclusions regarding the clinical phenotype of the D614G clade and there is no evidence at this time to suggest this clade is associated with any differences in disease phenotype 23 . As a potential proxy measure for viral load, we reviewed the cycle threshold (Ct) values of the SARS-CoV-2 assay for samples sequenced from Ontario, Canada (total n = 178, with 50 and 128 samples possessing the spike protein 614-D and 614-G amino acids, respectively). In this small selection of samples, we found no significant difference in the average Ct for 614-D and 614-G amino acids (mean 22.35, SD 5.13 and mean 22.32, SD 5.188 respectively, t-test p-value = 0.978) (Fig. 4). We do note, however, that numerous factors affect the measurement and limit the reliability of Ct value measurement, including sample type, quality of swabs, quality of collection method, and time from onset to sample collection 24 . In contrast, another study has found a statistically significant decrease in PCR Ct values associated with the mutated amino acid G variant 10 but its clinical significance remains unclear. In spite of the lack of evidence for phenotypic differences, it is interesting that in a short period of time since its emergence the D614G clade has become widespread all around the world. It is possible that even if the observed mutation does not impact the protein's interaction with ACE2, it may not be completely neutral with respect to viral fitness. For example, given that the molecular weight of glycine is significantly smaller than that of aspartate, the mutation could be advantageous from a cost minimization point of view 25 . It is also possible that the mutation may be in linkage disequilibrium with a selectively advantageous variant impacting another aspect of viral reproduction. For example, most strains of the D614G clade also harbor a mutation (P4715L) in orf1ab and a subclade processes three nucleotide changes in the nucleoprotein (N) gene (GGG to AAC; G204R), which plays diverse roles in virion assembly as well as genome transcription and translation 26,27 . Most recently, it has been suggested that D614G introduces a new elastase cleavage site that may be differentially activated by host genomic mutations thereby facilitating The receptor-binding domain is colored purple and the location of the aspartate residue in position 614 is highlighted in green. (C) Inter-atomic contacts between aspartate 614 (green) in a reference spike monomer (blue) and four residues (pink) in its adjacent spike protein monomer chain (white). These four contacts are destabilizing and create a hydrophilic-hydrophobic repelling effect that is lost upon replacement of aspartate by glycine in the D614G mutation (see Table 1). (D) Spatial distribution of aspartate 614 residue (green) and an adjacent glycosylated asparagine residue in position 616 (orange). The two residues point in opposite directions and thus it is unlikely they share a meaningful interaction. The image (A) was drawn using Affinity Designer (v1.8) (https ://affin ity.serif .com/en-gb/desig ner/). The trimeric and monomeric structures of the Spike protein were generated using Illustrate 19,41 (https ://ccsb.scrip ps.edu/illus trate /) by rendering a protein structure from the Protein Data Bank with ID 6vsb 19 (https ://www.rcsb.org/struc ture/6vsb). The image (B-D) was generated using UCSF Chimera (v1.14) (https ://www.cgl.ucsf.edu/chime ra/) with monomeric protein structure rendered in Chimera 19 . spike processing, and entry into host cells in some host populations but not others 28 . Finally, the global spread of the D614G mutation may have nothing to do with viral biology but may simply be a consequence of the high level of interconnectedness of Europe to the rest of the world. Thus, the emergence of the D614G clade may be explained by a founder event and subsequent clonal expansion in Europe that led to its spreading worldwide. Genomic sequencing and phylogenetic analysis are powerful approaches to track viral evolution during the course of a pandemic and help coordinate the global implementation of strategies for decreasing virus transmission and mitigating global mortality. However, our study highlights that caution is warranted in hastily drawing conclusions from limited observational data. Future work that can leverage clinical outcome data with both viral and human genomic diversity, in vitro and in vivo tests and in silico modeling will greatly enhance our ability to understand if specific SARS-CoV-2 clades are evolving with respect to virulence or transmissibility, and identify genetic variants associated with viral adaptation, transmissibility, and mortality. These data will be essential for effective vaccine development and public health responses to this emergent pandemic.

Sequence data, demographics, maps and statistics. A total of 2,803 genomes of SARS-CoV-2 (excluding non-human sequences specimens) and demographics when available were obtained from GISAID (https ://gisai d.org/) (downloaded on 30 March 2020) to perform the statistical analyses of the D614G mutation (8 sequences were excluded due to ambiguous or unknown nucleotide at this position). The dataset included sequences representing 55 countries. The world geographical maps were built with the geographic information system QGIS (v2.18.21, https ://qgis.org). Global CFR was calculated from data according to Johns Hopkins as of 30 March 2020. CFR for age groups of different countries were extracted from country-specific published data 8,12-14 . Spearman's rank correlation analysis was performed using GraphPad prism (v8.3.1).

and bronchoalveolar lavages) submitted to Public Health Ontario Laboratory was performed as previously described 29 . Specimens with Ct values of 32 or less were selected for whole genome sequencing. Briefly, rRNA depleted RNA was reverse transcribed with SuperScript III Reverse Transcriptase (Invitrogen/Life Technologies, Carlsbad, CA) to cDNA using tagged random primer (5′-GTT TCC CAG TCA CGATA-(N9)-3′) followed by second strand DNA synthesis using Sequenase (Thermo Fisher Scientific, Waltham, MA). Amplification of double stranded cDNA was done with primer 5′-GTT TCC CAG TCA CGATA-3′ using AccuTaq LA DNA Polymerase according to manufacture recommendation. Amplified products were purified using 1.0 × ratio of Agencourt AMPure XP beads (Beckman Coulter) and quantified with Qubit 2.0 Fluorometer by using Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). Nextera XT DNA Library Preparation kit (Illumina) following manufac- turer's instructions was used for Library preparation and indexing. The sequencing library was quantified using Qubit 2.0 (Invitrogen, Waltham, MA, USA) and average size was determined by 4200 TapeStation using the D5000 ScreenTape assay (Agilent Technologies, US). The library was normalized to 4 nMol and sequenced on the Illumina MiSeq platform, using Illumina MiSeq reagent kit v2 (2 × 150 bp), according to the manufacturer's instructions. CLC Genomics Workbench version 8.0.1 (CLC bio, Germantown, MD, USA) was used for both reference assembly and de novo assembly. Additionally, we used protocol published by Artic Network 30 . The method is based on overlapping specific primers producing short amplicon covering entire genome. Briefly, cDNA synthesis and amplicon generation using 2 individual pools of primers was done according to the published protocol 30 . The two pools of amplicons were combined together and cleaned by Agencourt AMPure XP beads (Beckman Coulter) prior to library preparation using Nextera XT (Illumina) following the manufacturer protocol except that half volume of the reagent was used throughout the protocol. Agencourt AMPure XP beads (Beckman Coulter) in ratio of 1.0 were used for the final purification step followed by measurement by Qubit 2.0 and average size determination, using 2,200 TapeStation (Agilent Technologies, USA). All libraries were manually normalized to 4 nMol. Libraries were combined in equal volumes for denaturation and subsequent dilution according to MiSeq protocol recommended by manufacturer.

A total of 2,803 SARS-CoV-2 genomes were obtained from GISAID (downloaded on March 30th 2020). Initial quality filtering was performed by iteratively removing redundant sequences with 100% nucleotide identity, having less than 2,940 informative positions, and lacking informative sampling dates. The resulting 2,017 curated sequences were aligned using mafft 31 (v7.450: -retree 2 -maxiterate 2 -auto settings), and unaligned ends and gap positions from the resulting alignment were removed using Geneious (vR10.2.6, https ://www.genei ous.com). A Maximum-likelihood tree was constructed using IQTree 32 (v1.6.1: -model HKY) and used for identification and removal of highly divergent sequences with Treetime 33 (clock setting). Following previous recommendations (adapted from https ://virol ogica l.org/t/issue s-with-sars-cov- 2-seque ncing -data/473) 34 to account for regions which might potentially be the result of hypervariability or sequencing artifacts, alignment positions showing significant homoplasy were first identified using Treetime (homoplasy setting). Top-10 significant homoplastic positions were merged from three separate analyses, using all sequences from the initial alignment, and sequence subsets generated using Illumina (1,450 samples) or Nanopore (341 samples) sequencing technologies. Next, recombinant clusters of mutations were identified with ConalFrameML 35 . Identified positions were masked using a custom-written python script, and further filtered for sequences with 100% nucleotide identity, resulting in an alignment of 1,225 sequences. Finally, to reduce sampling bias only one sample per country was kept for date of sampling, leading to 442 sequences. Molecular dating analysis was performed by the Markov chain Monte Carlo (MCMC) method implemented by Bayesian Evolutionary Analysis on Sampling Trees (BEAST) 1.10.4 36 . A HKY85 nucleotide substitution model was used. Fixed and relaxed molecular clock models were fitted, and constant population size was compared with exponential population growth. All models were run using default priors except for the exponential growth rate (Laplace distribution) in which scale was set to 100. The chain length was set to 100 million states and burn-in of 10 million. Convergence was checked with Tracer 1.7.1 37 . Based on ESS values, estimates for tMRCA of the European Clade spike protein mutation and root age are reported under the strict molecular clock and constant population growth (ESSs > 900). The resulting coalescent tree was generated using TreeAnnotator 38 and visualized using ggtree package 39 in R version 3.6.1 (https ://www.r-proje ct.org/).

Visualization, analysis and in silico mutations of protein structures were done in UCSF Chimera 18 . First, we downloaded the molecular structure of the spike protein from the Protein Data Bank (PDB). This structure corresponds to that resolved by Wrapp and colleagues 19 and deposited in PDB with identifier number 6vsb. We replaced the aspartate residue in position 614 by glycine using the "Rotamers" function in Chimera with default parameters. Then, inter-atomic contacts of both residues at position 614 were derived with the CSU package 19 freely available at https ://oca.weizm ann.ac.il/oca-bin/lpccs u.

Ethics statement. Specimens were collected from patients by submitters and sent to PHOL for testing as part of routine clinical service. These data are also used for routine laboratory surveillance, which is a mandate of Public Health Ontario. Therefore, consultation with our organization's privacy office or ethics committee was not required. To protect patient privacy and confidentiality, data are reported in an anonymized format.