A SARS-CoV-2 vaccine candidate would likely match all currently circulating variants

To characterize SARS-CoV-2 diversification since the beginning of the epidemic, we aligned 27,977 SARS-CoV-2 genome sequences isolated from infected individuals in 84 countries. The alignment was curated to retain independent sequences that covered over 95% of the ORFs. In addition, because sequences from the United Kingdom constituted 47% of the dataset (n = 12,157), we sampled a representative set of 5,000 UK sequences, yielding a final dataset of 18,514 SARS-CoV-2 genomes (SI Appendix, Fig. S1 and Fig. 1A).

There were 7,559 polymorphic sites (that is, sites where at least one sequence has a change relative to the reference sequence) across the genome (total length: 29,409 nucleotides). Most substitutions were found in a single sequence; only 8.41% (n = 2,474) of the polymorphic sites showed substitutions in two or more sequences (Fig. 1B). Only 11 mutations were found in >5% of sequences, and only 7 were found in >10% of sequences (3 of which were adjacent). The mean pairwise diversity across genomes was 0.025%, ranging between 0.01% for E to 0.11% for N. A phylogenetic tree reconstructed based on all genome sequences reflected the global spread of the virus: Samples from the first 6 wk of the outbreak were collected predominantly from China (Fig. 1C). As the epidemic has progressed, samples have been increasingly obtained across Europe and from the United States (Fig. 1 A and C). The tree shows numerous introductions of different variants across the globe, with introductions from distant locations seeding local epidemics, where infections sometimes went unrecognized for several weeks and allowed wider spread (23). Prior to the severe travel restrictions that were seen in March 2020, intense travel patterns between China, Europe, and the United States allowed transmission of a myriad of variants, which is currently reflected by different lineages in the tree. Yet, the tree topology shows minimal structure, even at the genome level, indicating that SARS-CoV-2 viruses have not diverged significantly since the beginning of the pandemic. To compare how genomes differed from one another, we calculated Hamming distances (which correspond to the number of differences between two genomes) across all pairs of sequences. These Hamming distances showed a narrow distribution, with a median of seven substitutions between two independent genomes, while linked sequences sampled in cruise ships had a median of two substitutions (SI Appendix, Fig. S2). Surprisingly, Hamming distances across genomes sampled in the United States did not show a similar quasi-normal distribution but instead a bimodal distribution, observed despite the large number of sequences compared (n = 5,398). We identified that this bimodal distribution was driven by sequences from Washington State, possibly reflecting separate introductions in that state. Nonetheless, such a bimodal distribution could also indicate a bias in the sampling of sequences (SI Appendix, Fig. S2).

Since the beginning of the pandemic, two mutations across the genome have become consensus: P4715L in ORF1ab (nucleotide 14,143, C to T) and D614G in S (nucleotide 23,403, A to G) (Fig. 1B) (a third consensus mutation, at nucleotide 3,037, is not reported as the site was masked during our sequence-filtering procedure). These mutations were found in 69.3% and 69.4% of sequences, respectively, and are in linkage (Fig. 2B). Given the importance of S for virus entry and as a target of the host neutralizing response, the biologic implications of the D614G mutation are under intense scrutiny (24–28). This mutation was first observed in a sequence from China dated January 24, with seven more sequences sampled until February 8. Then, the D614G mutation was not observed in China until March 13. In contrast, the D614G mutation was introduced in Europe at the end of January (first sequence identified in Germany, dated January 28), and it rapidly became dominant on that continent and at every location where the virus subsequently spread (Fig. 2A). The phylogenetic tree (Fig. 2B) and the distribution of sequences (Fig. 2C) are suggestive of a founder effect. The rapid spread of sequences carrying the D614G mutation implies that the growing viral population should become more homogeneous, that is, the frequency of sequences with shared polymorphisms will increase. We found a median of seven substitutions (based on a comparison of 18,514 sequences) between two independent SARS-CoV-2 genomes (SI Appendix, Fig. S2). Yet, genomes with the D614G mutation showed a median of five substitutions, whereas those with D at position 614 differed by eight substitutions (Fig. 2D).

To test whether this site was under selection, we used likelihood-based, phylogenetically informed models that assume branch-specific substitution rates (29) and implemented a sampling strategy to circumvent computational limitations imposed by the large number of sequences. Subsampled alignments (100 times at a 10% sampling fraction) had diversity estimates statistically similar to the complete alignment for each gene (Mann−Whitney U test, P > 0.09; SI Appendix, Fig. S3). In S, only site 614 was estimated to be under diversifying selection in a majority of subsampled alignments (58%); evidence of diversifying selection indicates that genetic diversity increases in the viral population (i.e., there was a higher proportion of mutations causing an amino acid change than not at site 614, or, the nonsynonymous/synonymous substitution rates ratio, dN/dS, was over 1, P < 0.1) (SI Appendix, Fig. S4). Because diversifying selection is often associated with the host adaptive response, we considered whether the D614G mutation coincided with targets of antibody and T cell responses. Site 614 is at the interface between the S1 and S2 subunits and is thus not highly accessible to antibodies (SI Appendix, Fig. S5). Therefore, we predict that antibodies to the native S protein would cross-react with S containing the D614G mutation, in agreement with recent reports (24, 25, 27, 28). Many known neutralizing antibodies target the RBD, yet we found little evidence that mutations could affect binding to the ACE2 receptor, as only five shared mutations were identified at contact sites with the ACE2 receptor, and all were found in 10 or fewer sequences. Of these, one mutation, at position 489, was synonymous and found in three sequences (0.02%). The others were nonsynonymous: G476S (n = 10 sequences, 0.05%), Y453F (n = 5, 0.02%), G446V (n = 3, 0.02%), and A475V (n = 2, 0.01%). To predict the potential immune pressure linked to T cell responses, we developed a T cell immunogenicity index which takes into account the CD8 and CD4 epitope repertoires in the structural proteins of SARS-CoV-2 (S, N, M, E) and the frequency of human leukocyte antigen (HLA) alleles or haplotypes in a given population. We found that sites with mutations, including 614 in S, were not colocalized with T cell epitopes frequently identified in different populations (SI Appendix, Figs. S6 and S7), and there was no significant relationship between the number of mutations and the immunogenicity index (SI Appendix, Fig. S8).

Using phylogenetically informed models (as described above), we identified two sites, residue 614 in S and 13 in N, that were under diversifying selection in a majority of subsampled alignments. For each protein, subsampled alignments tended to have more sites under purifying selection (median = 7.34 ± 4.06% [±SD]) than under diversifying selection (3.10 ± 1.92%) (Mann−Whitney U test, P = 0.057; SI Appendix, Fig. S4) (purifying selection is indicative of a decrease in genetic diversity in the population). Likewise, for each codon separately, the proportion of each phylogeny (i.e., the percentage of total branch length) with dN/dS > 1 was small, indicating diversifying selection was episodic and limited (Fig. 3A). Global measures of dN/dS varied across genes, ranging from 0.35 ± 0.02 (M) to 1.43 ± 0.24 (ORF10), and were significantly lower for structural genes compared to nonstructural genes (Mann−Whitney U test, P = 0.042) (Fig. 3B). Per-lineage nonsynonymous substitution rates were comparable (Student’s t test, P = 0.218) in structural (0.0011 ± 0.021) and nonstructural (0.0012 ± 0.028) genes, although some subsampled alignments showed rates that could be a hundred times higher than the median over all alignments (Fig. 3C). Across structural proteins, mutations were disproportionately neutral: >70.3% of branch length evolved under neutral (or negative) selection for all sites, and over half of all branch length evolved under neutral (or negative) selection for >82.8% of sites (Fig. 3D) (29).

While there was only limited evidence of diversification at selected sites, we also assessed whether subpopulations among the globally circulating viral population had become genetically differentiated over time. To do so, we used two measures of population differentiation, the G_ST and D statistics, which characterize changes in allele frequency across populations and can show fitness differences between subpopulations (30–32). Genetic distances between two subpopulations can range between 0 and 1, indicating no and complete differentiation, respectively. We initially compared 30 genomes sampled from the initial outbreak in Wuhan, China, with subsampled alignments of the 18,484 genomes sampled subsequently across the globe. Although distances varied across genes, the median genetic distance between these subpopulations was small for both G_ST (0.0049 ± 0.0047) and D (0.0053 ± 0.0272), indicating little differentiation between the initial outbreak and its global derivatives in the pandemic (Fig. 4A). We then compared subpopulations sampled before and after each consecutive week. Similarly, genetic distances between subpopulations were small for both G_ST (0.0058 ± 0.0096) and D (0.0098 ± 0.0650) and tended to narrow over time rather than diverge (Fig. 4 B and C). Signatures of host adaptation can also be seen in the branching patterns of viral phylogenies. Bursts in transmissibility are emblematic of increases in relative viral fitness and are reflected in imbalances in the phylogeny, which can be estimated at each internal node (SI Appendix, Figs. S9 and S10) (33–35). We estimated phylogenetic η (36, 37) at each internal node of the SARS-CoV-2 phylogeny reconstructed from subsampled (10%) alignments and compared the distribution of estimates through time to phylogenies simulated under models of neutral and positive time-dependent rates (b(t) = be^α(t)). Simulation analyses demonstrated that this metric was robust against sampling fraction (SI Appendix, Fig. S10). The distribution of η in the SARS-CoV-2 phylogenies adhered to expectations of the neutral model and deviated significantly (Student’s t test, P < 0.001) from those of positive time-dependent rates for selection coefficients α ≥ 0.2 (Fig. 4 D and E). Together, the SARS-CoV-2 population and phylogenetic dynamics showed little evidence that the global spread of SARS-CoV-2 was related to viral fitness effects.

Typical vaccine design strategies rely on either 1) selecting sequences sampled from infected individuals or 2) computationally deriving sequences that cover the diversity seen across circulating sequences and are, in theory, optimal compared to an individual isolate (38). Computationally derived sequences include consensus and ancestral sequences, such as the most recent common ancestor (MRCA) of a set of sequences. We inferred the MRCA corresponding to 1) SARS-CoV-2 S sequences sampled from Wuhan within the first month of the epidemic, 2) all currently circulating SARS-CoV-2 sequences, and 3) all SARS-CoV-2 sequences together with closely related sequences sampled from pangolins (n = 6) and a bat. There were 17 mutations between the human MRCA and the human−bat MRCA and 44 mutations between the human MRCA and human−pangolin MRCA. Overall, three segments in S reflected significant variability across species (AA 439 to 445, 482 to 501, and 676 to 690) (Fig. 5 A and C and SI Appendix, Fig. S11) (39). In contrast, when considering only human sequences, SARS-CoV-2 diversity was limited: Both MRCAs (derived from early sequences from Wuhan or from all circulating sequences) were identical to the initial reference sequence Wuhan-Hu-1. Comparing these sequences to the consensus sequence derived from all of the sequences sampled to date, there was only one mutation: D614G (Fig. 5B). Fig. 5D illustrates that mutations found across circulating S sequences were rare: Besides D614G (found in 69.4% of sequences), the next most frequent substitution is found in 1.96% of sequences (synonymous), with sequences sampled from infected individuals, on average, 0.55 mutations away from the consensus sequence (consisting of 0.12 synonymous and 0.43 nonsynonymous mutations). Across the genome, there were, on average, 4.05 nucleotide mutations per individual genome when compared to the consensus, with only P4715L and D614G found in >50% of sequences.

Sequences were downloaded from GISAID (https://www.gisaid.org/). A full list, along with the originating and submitting laboratories (GISAID_acknowledgment_table_20200518.xls), is available at https://www.hivresearch.org/publication-supplements.

All SARS-CoV-2 sequences available on GISAID as of May 18, 2020 (n = 27,989) were downloaded and deduplicated where possible, and those missing accurate dates (that is, only recording the month and/or year) were removed. Sequences were processed using the Biostrings package (version 2.48.0) in R (49). Sequences known to be linked through direct transmission were removed, and only the sample with the earliest date (chosen at random when multiple samples were taken on the same day) was retained. Sequences were then aligned with Mafft v7.467 using the -addfragments option to align to the reference sequence (Wuhan-Hu1, GISAID accession EPI_ISL_402125) (50). Insertions relative to Wuhan-Hu-1 were removed, and the 5′ and 3′ ends of sequences (where coverage was low) were excised, resulting in an alignment consisting of the 10 ORFs. Any sequences with less than 95% coverage of the ORFs (i.e., >5% gaps) were removed, and 30 homoplasic sites likely due to sequencing artifacts identified by de Maio et al. were masked (https://github.com/W-L/ProblematicSites_SARS-CoV2/blob/master/archived_vcf/problematic_sites_sarsCov2.2020-05-27.vcf).

To identify individual sequences that were much more divergent than expected, given their sampling date, which likely reflected sequencing artifacts rather than evolution, we obtained a tree using FastTree v2.10.1 compiled with double precision under the general time reversible (GTR) model with gamma heterogeneity (51). This tree was rooted at the reference sequence, and root-to-tip regression was performed following TempEst using the ape package in R (52, 53). Outliers were defined as sequences that had studentized residuals greater than 3, and were removed.

Sequences from the United Kingdom corresponded to nearly half of the sequences (n = 12,157/25,671, 47%) of this filtered dataset. To avoid overrepresentation of the UK sequences and bias in subsequent analyses, we investigated the effect of downsampling sequences on the mean Hamming distance and identified the minimum number of sequences required to recover the mean corresponding to the full distribution (SI Appendix, Fig. S1). A subsample of 5,000 sequences satisfied these criteria, and also ensured that there were fewer sequences from the United Kingdom than from the United States (n = 5,398), reflecting the epidemiology. These 5,000 sequences were sampled randomly, with weight proportional to the number of UK sequences collected on that day.

After these filtering steps, the alignment used for subsequent analyses included 18,514 sequences.

The global phylogeny was reconstructed in FastTree v2.10.1 compiled with double precision under the GTR model with gamma heterogeneity (51), and rooted at the reference sequence. The tree was visualized using ggtree in R (54). Lineages were defined using PANGOLIN (Phylogenetic Assignment of Named Global Outbreak LINeages), with lineages with >200 taxa as of the May 19 summary being highlighted in the tree (22) (https://github.com/cov-lineages/lineages). The number of polymorphic sites was calculated as the number of sites which had at least one mutation relative to the reference sequence, Wuhan-Hu-1, ignoring gaps and ambiguities.

For each pair of sequences, we calculated the Hamming distance as the number of sites that are different after removing sites with ambiguities and/or gaps. For computational efficiency, given the size of the alignment, this was implemented in parallel in C++, using Bazel (https://bazel.build/) to build on a Linux system. This implementation is available to download at https://www.hivresearch.org/publication-supplements.

Alignments for each gene were subsampled for sequence and phylogenetic analyses. Each gene alignment was randomly subsampled 100 times per collection date at 5%, 10%, 20%, 30%, and 40%. When fewer than 10 sequences were available for a collection date, all sequences were taken. Median Hamming distances were computed for each set of subsampled alignments. These were bootstrapped 100,000 times, and 95% CIs were estimated and compared to the median Hamming distance for the fully sampled alignment.

Alignments subsampled at 10% 100 times were used to estimate substitution rates. For the set of subsampled alignments for each gene, a mixed-effect likelihood method was used to estimate nonsynonymous (dN) and synonymous (dS) substitution rates globally and at each codon (29). Maximum-likelihood phylogenies were constructed for each alignment using the software IQ-TREE (55) under a best-fit model determined with ModelFinder (56) to prime the dN and dS estimates before branch length optimization. This step serves to expedite the optimization process. Branch length optimization was done with a MG94 model [which is the only model available for this analysis (29)]. The proportion of each phylogeny evolving under neutral (or negative) selection was determined from the mixture density across lineages for each site, assuming different dN and dS along each branch (57). On the same set of subsampled alignments and phylogenies, a fixed-effects likelihood method was used on internal branches to identify sites under pervasive diversifying selection and to estimate global dN/dS (58). Known biases associated with calculating dN/dS on exponentially growing populations (59) were counterbalanced by subsampling phylogenies, as the typical approach to address this bias, which is to ignore terminal branches, would considerably diminish the power of the analysis to detect any significant result. As P values from the fixed-effect likelihood method are uncorrected, results were not averaged over P values; rather, given that P value calculations are conservative for this analysis (58), sites were considered to be under pervasive diversifying selection if their P value was <0.1 in $\ge$50% of alignments, which would account for a typical 5% false discovery rate (58).

Alignments subsampled at 10% 100 times were used to estimate population differentiation. The genetic differentiation of subpopulations within sampled sequences was calculated on each gene separately using Nei’s (30) G_ST. Because comparisons between subpopulations of different sizes can bias genetic differentiation estimates (60), genetic differentiation was also calculated using Jost’s (31) D, which accounts for differences in genetic heterogeneity between subpopulations and is intended to correct for biases in the size of the subpopulations. Both statistics were computed with the mmod package (32) in R (v3.6.1). For each gene, statistics were calculated over 100 bootstrapped samples for each subsampled alignment. Subpopulations were defined in two ways. First, sequences originating from the initial outbreak in the Hubei province (30 sequences) were compared to all other sequences within a subsampled alignment. Second, a 1-wk sliding window was designed to compare all sequences sampled prior to a collection date (subpopulation 1) to all sequences sampled after the same collection date (subpopulation 2). The first collection date for subpopulation 1 was February 14, 2020, the week after the last sequence from the Hubei province was sampled (February 8, 2020), The window was designed to terminate when <30 sequences were available in subpopulation 2.

Time-dependent estimates of phylogenetic diversification were measured by extracting the branches descending from each internal node (above the root) of each phylogeny and calculating the peak height (η) of the spectral density profile of the graph Laplacian of each subtree, which is a measure of the density of branching events (36, 37). The code to perform the analysis is available for download at https://www.hivresearch.org/publication-supplements.

Ancestral S protein sequences were reconstructed from an amino acid alignment of 30 SARS-CoV-2 sequences sampled from the Hubei province, a coronavirus sampled from bat (Yunnan RaTG13), and six SARS-CoV-2-like coronaviruses sampled from pangolins using maximum posterior probability and returning a unique residue at each site assuming a Jones-Taylor-Thornton (JTT) model with gamma heterogeneity (62). The JTT model was the most appropriate model available in the software (62). The bat sequence was retrieved from GenBank, and the pangolin sequences were retrieved from GISAID (63). A sliding window of 10 amino acids (and a step of 1 amino acid) was used to compare the cumulative number of mutations in the human−bat and human−bat−pangolin ancestors with respect to the human ancestral sequence. Median values for each window were compared to a null window (computed as a normal distribution of 10 values with a mean equal to the mean value across the entire S protein, 0.046 mutations) using a one-tailed t test. An alignment including the reconstructed sequences is available at https://www.hivresearch.org/publication-supplements.

CD4+ and CD8+ T cell epitopes were predicted for four SARS-CoV-2 structural proteins: S (accession YP_009724390), N (accession YP_009724397), M (accession YP_009724393), and E (accession YP_009724392). CD4+ T cell epitopes were predicted using a server that predicts binding of peptides to any MHC molecule of known sequence using artificial neural networks, NetMHCIIPan 4.0 (64) with a peptide length of 15. MHC class II HLA alleles of HLA-DQB1, plus the haplotypes of HLA-DPA1-DPB1 and HLA-DQA1-DPB, were selected for predictions if they had frequencies of >1/1,000 in known allele/haplotype distributions (http://17ihiw.org/17th-ihiw-ngs-hla-data/). If multiple peptides had the same core, the peptide with the strongest binding score was selected for analysis. CD8+ T cell epitopes were predicted using NetMHCPan 4.1 (64) with a peptide length of 9. MHC class I HLA alleles of HLA-A, HLA-B, and HLA-C were selected if they were classified as common (frequency ≥ 1/10,000) in any of the populations in the database CIWD 3.0 (Common, Intermediate and Well-Documented HLA Alleles in World Populations) (65). Epitopes predicted as strong binders (with predicted binding affinities below 50 nM) were selected for analyses.

For each site in a predicted epitope, the immunogenicity index was defined as the sum of the frequency of the HLA alleles or haplotypes restricting the corresponding epitope (multiple epitopes can be predicted at a given site in a protein). Total frequencies from CIWD 3.0 were used as the frequencies of the corresponding MHC class I HLA alleles (HLA-A, HLA-B, and HLA-C), and the global frequencies from http://17ihiw.org/17th-ihiw-ngs-hla-data/ were used as the frequencies of the corresponding MHC class II HLA alleles or haplotypes (HLA-DQB1, HLA-DPA1-DPB1, and HLA-DQA1-DPB). This procedure was repeated using the frequencies of MHC alleles or haplotypes in different subpopulations listed in the above HLA frequency dataset.

For comparisons of mean values in normally distributed data, Student’s t test was used. When data were not normal, the Mann−Whitney U test was used. Shapiro−Wilk tests were used to determine normality. Differences in data distributions were estimated using the Kolmogorov−Smirnov test.