Research Article

SARS-CoV-2 within-host diversity and transmission

See allHide authors and affiliations

Science  09 Mar 2021:
eabg0821
DOI: 10.1126/science.abg0821

Abstract

Extensive global sampling and sequencing of the pandemic virus SARS-CoV-2 have enabled researchers to monitor its spread, and to identify concerning new variants. Two important determinants of variant spread are how frequently they arise within individuals, and how likely they are to be transmitted. To characterize within-host diversity and transmission we deep-sequenced 1313 clinical samples from the UK. SARS-CoV-2 infections are characterized by low levels of within-host diversity when viral loads are high, and a narrow bottleneck at transmission. Most variants are either lost, or occasionally fixed, at the point of transmission, with minimal persistence of shared diversity - patterns which are readily observable on the phylogenetic tree. Our results suggest that transmission-enhancing and/or immune-escape variants are likely to arise infrequently, but could spread rapidly if successfully transmitted.

The ongoing evolution of SARS-CoV-2 has been the topic of considerable interest as the pandemic has unfolded. Clear lineage-defining single nucleotide polymorphisms (SNPs) have emerged (1), enabling tracking of viral spread (2, 3), but also raising concerns that new mutations, or combinations of mutations, may confer selective advantages on the virus, hampering efforts at control. There is compelling evidence that the D614G mutation in the Spike protein, which spread globally during the first year of the pandemic, increases viral transmissibility (46). Current variants of concern include the B.1.1.7. lineage (7, 8), with an estimated transmission advantage of around 50% (9), and the B.1.351 and P.1 lineages (10, 11), which may have decreased sensitivity to natural and/or vaccine acquired immunity (1214). Lineage codes given here are as designated by the Pangolin software (1).

Most analyses have been focused on mutations observed in viral consensus genomes, which represent the dominant variants within infected individuals. Ultimately though, new mutations emerge within individuals, and hence knowledge of the full underlying within-host diversity of the virus at the population level, and how frequently this is transmitted, is important for understanding adaptation and patterns of spread.

The United Kingdom (UK) experienced one of the most severe first waves of infection, with over a thousand independent importation events contributing to substantial viral diversity during this period (15). In this study, we analyzed 1390 SARS-CoV-2 genomes from 1313 nasopharyngeal swabs, sampled predominantly from symptomatic individuals on admission to hospital and from healthcare workers, during the first wave of infection (March - June 2020; Table S1). The dataset comprised samples from 1173 unique individuals, including 41 with samples at 2-4 timepoints, plus 93 anonymous samples, with multiple RNA aliquots from 76/1313 samples resequenced to test for reproducibility. The samples were collected by two geographically separate hospital trusts: Oxford University Hospitals and Basingstoke and North Hampshire Hospital, located 60 km apart. Using veSEQ, an RNA-Seq protocol based on a quantitative targeted enrichment strategy (16), which we previously validated for other viruses (1619), we characterized the full spectrum of within-host diversity in SARS-CoV-2 and analyzed it in the context of the consensus phylogeny.

We observed low levels of intrahost diversity in high viral load samples, with evidence of within-host evolutionary constraint genome-wide, including Spike. Although within-host variants can be observed in multiple individuals in the same phylogenetic cluster, some of whom reside in the same household, most viral variants are either lost, or occasionally fixed, at the point of transmission, with a narrow transmission bottleneck. These results suggest that during early infection, when viral loads are high and transmission is most likely (2022), mutations that increase transmissibility, or potential vaccine- or therapy-escape mutations, may rarely emerge and subsequently transmit. Nonetheless, we identified variants present in multiple individuals that can affect receptor binding or neutralization by antibodies. Since the fitness advantage of escape mutations in highly-vaccinated populations, or those with high levels of natural immunity, could be substantial, and additionally because mutational effects can depend on the genetic background on which they are found, these findings underline the need for continued vigilance and monitoring.

Detection of variants is influenced by viral load

Reliable estimation of variant frequencies requires quantitative sequencing, such that the number of reads is proportional to the amount of corresponding sequence in the sample of interest. The veSEQ protocol has been previously shown to be quantitative for a number of different pathogens (17), including respiratory viruses such as RSV (18). We demonstrated the same quantitative relationship holds for SARS-CoV-2. The number of uniquely mapped sequencing reads we obtained rose log-log linearly with the number of RNA copies in serial dilutions of synthetic RNA controls (Fig. S1A, r2=0.87), and was consequently correlated with cycle threshold (Ct) values of clinical samples (Fig. S1B), indicating that veSEQ reads can be considered a representative sample of viral sequences within the input RNA.

To understand within-host diversity, we quantified the number of iSNVs in the full set of 1390 genomes, testing different thresholds for identifying variants of between 2 and 5% minor allele frequency (MAF). A minimum depth of at least 100 reads was also required to call an iSNV, and all sites with MAF greater than the threshold were included (Fig. 1A).

Fig. 1 Characterization of iSNV frequencies.

(A) Distribution of number of identified iSNV sites in each sample against number of unique mapped reads. The colors represent different minor allele frequency (MAF) thresholds; an iSNV site is identified within a sample if the MAF is greater than the threshold. (B) Distribution of mean MAF in each sample against number of unique mapped reads, with no MAF threshold applied. The black line is the estimated mean value by linear regression, with the green ribbon the 95% confidence interval. (C) Distribution of number of identified iSNV sites, at 3% MAF threshold, when subsampling from high-depth samples. Each color represents a different high-depth sample.

For all thresholds, we observed a non-linear relationship between sample viral load (estimated by total unique mapped reads) and the number of detected iSNVs, with the highest number of iSNVs detected at intermediate viral loads (~2000 mapped reads). However, the mean MAF per sample did not vary with viral load when no threshold was applied (p=0.291, linear regression, Fig. 1B). This indicates that as the number of mapped reads decreases, the variance in the observed MAF increases, while the mean stays the same. This effect is at least partially due to the inverse relationship, for the binomial distribution, between the total number of draws and the variance in the proportion of successes observed amongst those draws. In Fig. 1C we demonstrate this effect by downsampling from high-depth samples: the increasing variance associated with sparser sampling causes the number of threshold-crossing iSNVs to increase, until eventually so few reads are sampled that no iSNVs are detected.

This sampling effect of low viral load does not preclude the existence of biological mechanisms also contributing to greater intrahost diversity in low-viral load samples. After the initial peak, viral loads typically decrease as infection progresses (20), while genetic diversity may increase, as observed in other viral infections, such as HIV (23). RNA damage (24) as infection progresses could also contribute to increased observed diversity in low depth samples.

Within-host variant frequencies are reproducible

To calibrate our variant calling and minimise false discovery rates, we compared intrahost single-nucleotide variants (iSNVs) in re-sequenced controls with data for the stock RNA sequenced and provided by the manufacturer (Twist Bioscience) and masked sites vulnerable to in vitro generation of variants (Table S2). We also masked a further 18 sites which were observed to be variant (>3% MAF) in 20 or more high-viral load samples (Table S3, Fig. S3). Most have consistently low MAFs among samples, and some showed evidence of strand bias and/or low reproducibility between technical replicates (Fig. S2), suggesting they are not true genomic variants. Among the excluded sites was 11083, which was observed in 46 samples and is globally ubiquitous in GISAID data. From manual examination of mapped reads in our dataset, this appears to be due to a common mis-calling of a within-host polymorphic deletion upstream at site 11082, occurring in a poly-T homopolymeric stretch. If genuine, this homopolymer stutter may have a structural or regulatory role; however, methodological issues in resolving this difficult-to-map region cannot be ruled out.

Establishing reliable variant calling thresholds for clinical samples, where true variant frequencies are unknown, ideally requires re-sequencing of multiple samples from RNA to test for concordance. Working within the constraints of small volumes of remnant RNA from laboratory testing, we re-sequenced 76 high viral load samples, of which 27 replicate pairs generated sufficient read numbers (>50,000 unique mapped reads) for reliable minor variant detection. Intrahost single-nucleotide variants (iSNVs) with <2% MAF were generally indistinguishable from noise, whereas those >=3% MAF were highly concordant between replicates (Fig. 2A, Fig. S2).

Fig. 2 Comparison of allele frequencies between sequencing replicates of the same sample and multiple time points from the same individual.

(A) Comparison of MAFs from 27 replicate pairs resequenced from RNA, with each point representing a single genomic position in a pair of replicates. The plot represents all MAF frequency comparisons for the 27 samples where both replicates had >50,000 unique mapped reads, limited to genomic sites where the MAF > 0.02 in at least one of the 54 replicates. The green lines are the threshold value of 0.03. (B and C) Comparison of allele frequencies from 41 individuals sampled on different days, with each point representing a genomic position in a pair of samples from the same individual. Each individual is represented by a different color, and for each individual all genomic positions are considered where the MAF >0.03 at either sampling time point, and/or a change in consensus was observed. In all cases, the poly-A tail and sites variable in RNA synthetic controls were excluded, as were sites observed to be variable in over 20 samples at MAF>3% since these are unlikely to represent genomic variants. (C) is an enlargement of the region of (B) near the origin.

Within-host variants vary during infection

In addition, we compared iSNV frequencies, and consensus changes, at different time points for the 41 multiply-sampled individuals, with the duration between sampling ranging between 1 and 20 days apart (median 6 days, Fig. 2B/C). Because viral loads tend to fall as infection progresses, we considered all samples rather limiting to those with >50,000 unique mapped reads. Among the 41 individuals, we observed little concordance in minor variant frequencies across time points within individuals. Our observations, consistent with other studies (2426), suggest a dynamic within-host landscape, but also reflect the inherent stochasticity associated with low viral load samples.

The transmission bottleneck size within households is small

The transmission bottleneck size is a key component in determining the likelihood that new within-host variants will spread in the population (27). Estimating bottleneck size is difficult for SARS-CoV-2, since it requires sufficient genetic diversity to differentiate distinct viruses that may be transmitted in known source-recipient pairs (2831), and confidence that transmission is the cause of variants observed in both source and recipients; inclusion of variants that are not shared by transmission can greatly increase transmission bottleneck size estimations (29). We identified 16 households where two individuals had a first positive sample within two weeks of each other, and assumed direct transmission if the consensus sequences in the individuals had fewer than three differences (hence excluding one household). A further household was excluded since the assumed source individual had no variants above 3% MAF.

Using the exact beta-binomial method (28) we estimated maximum likelihood bottleneck sizes between 1 and 8 among the 14 household transmission pairs (Fig. 3A, Table S4). These observations are consistent with the small bottleneck sizes observed for influenza (3032), and SARS-CoV-2 (3337), but considerably lower than estimates in a recent Austrian study (25). The reasons for the discrepancies are unclear, but could reflect differences in how variants were selected for analysis (37), or how closely observed diversity represents the diversity of virus both available for transmission and successfully transmitted. An association between the route of exposure and the transmission bottleneck has been demonstrated experimentally for influenza (32), and therefore genuine differences in bottleneck sizes in different settings cannot be ruled out.

Fig. 3 Small transmission bottleneck size within households.

(A) Estimated bottleneck size in 14 households calculated using the exact beta-binomial method described in (28). Bottleneck size for both combinations of potential source and recipient were calculated if the first positives samples from each individual in the household were collected within a week of each other. No estimate was recorded if there are no identified iSNVs >3% MAF in the source individual (household 8), or the two individuals in the household had more than two consensus differences (household 15). The error bars represent the 95% confidence interval determined by the likelihood ratio test. (B) The fate of identified iSNVs within households. Each line links the allele frequency of a given variant in one household member with that in the second member; points and lines are colored by household. Each was identified as an iSNV in at least one individual, but not necessarily (and usually not) both. Where the dates of sample collection differed by at least a week, we also indicate the assumed source and recipient members of the household.

Within-host variants are present in the majority of SARS-CoV-2 samples

To further characterize iSNV sites within individuals, we identified a set of 563 high-confidence iSNV sites that were observed (i) in high-viral load samples with at least 50,000 unique mapped reads (462 samples, 160 from Oxford and 302 from Basingstoke), (ii) at depth of at least 100 reads, (iii) with a MAF of at least 3%, and (iv) not observed to vary in synthetic RNA controls or appear at low frequency in a large number of samples (Table S3). All 1313 samples were included in our analysis under the assumption that by ascertaining on a small set of pre-defined sites, it is less likely that that we include sites which only reach >3% MAF in low viral load samples due to stochastic sampling effects described above.

Among the iSNV sites taken forward for variant analysis, most were only observed in one or two of the 1313 samples (Fig. 4A), but with the majority of samples with >50,000 unique reads (305/462, 66%) harboring at least one iSNV (Fig. 4B). These low levels of SARS-CoV-2 within-host diversity during acute infection are consistent with other reported levels (26, 33), but lower than in some other studies (24, 25), likely reflecting how variants were identified.

Fig. 4 Intra-host variable (iSNV) sites are often found in multiple samples and most samples have at least one iSNV.

(A) Histogram showing the number iSNV sites that are found in N samples. All samples in our dataset are included. (B) Stacked histogram showing the number of samples that have n iSNV sites for all samples with more than 50,000 mapped reads (dark red) and samples with fewer than 50,000 mapped reads (light red). All 563 sites identified for variant analysis were included (see main text), including sites in the 3′ and 5′UTR, but excluding the polyA tail, and excluding the 18 sites variable in 20+ individuals.

Two samples had a particularly high number (15 and 18) of iSNVs, each with high and correlated MAFs consistent with co-infection by two diverse variant haplotypes (38). For one of these samples, laboratory contamination is unlikely since we could not identify any samples that could be the source. We could not distinguish between co-infection and contamination in the other sample since both variant haplotypes within it represent common genotypes in our study.

In general however, the low level of genetic diversity of the virus makes identifying co-infection or contamination, and distinguishing between them, difficult. If sites where a large number of SNPs are present (mutations that distinguish common lineages in our dataset) are only observed to be variant within-host due to co-infection or contamination, then we estimate between ~1 to 2% of samples are potentially affected by co-infection or contamination (Table S2). As a precaution against contamination or batch effects, we sequenced known epidemiologically linked samples in different batches where possible (Fig. S4).

We hypothesized that a proportion of observed within-host variation could be due to co-infection with seasonal coronaviruses, which has been reported in 1-4% of SARS-CoV-2 infections (39, 40). Specifically, closely-matching reads from similar viruses could be mapped to SARS-CoV-2 and appear as mixed base calls. To understand the impact of co-infection, we re-captured and analyzed a random subset of 180 samples spanning the full range of observed SARS-CoV-2 VLs (Ct 14 to 33, median 19.8), using the Castanet multi-pathogen enrichment panel (17), which contains probes for all known human coronaviruses with the exception of SARS-CoV-2. Among the 111 samples that yielded both SARS-CoV-2 and Castanet data, we identified one sample that was also positive for another betacoronavirus, human coronavirus OC43 (Fig. S5). Within the SARS-CoV-2 genome from this sample, which was complete and high-depth, we observed only a single iSNV at position 28580 and no evidence of mixed base calls at any other genomic position. This suggests that even where co-infection is present, it does not impact on the estimation of SARS-CoV-2 within-host diversity in our protocol. However, whether co-infection with OC43 or other coronaviruses exerts a selective pressure on SARS-CoV-2 remains an open question.

Distribution of iSNVs across the genome

We next considered the distribution of the identified high-confidence iSNV sites across the genome. Even excluding the UTR regions, which have a highly elevated density of iSNV sites, there is considerable variability across the genome, with open-reading frames (ORFs) 3a, 7a, and 8, and nucleocapsid (N) showing the highest densities (Table 1). In addition, we calculated dN/dS values, under the assumption that each iSNV was only generated once de novo, and subsequently transmitted (Table 1). Consistent with other studies (24, 33), most areas of the genome appear to be under purifying selection, with dN/dS values less than 1, including the Spike protein, S. Without a full model incorporating within-host evolutionary dynamics and transmission, it is difficult to draw strong conclusions. However, we obtained similar results assuming each iSNV appeared de novo in each individual it was observed (Table S5). These patterns are also broadly consistent with dN/dS values calculated for SNPs among SARS-CoV-2 consensus genomes (41), suggesting evolutionary forces at the within-host level are reflected at the between-host level, at least for within-host variant sites in high viral load samples.

Table 1 iSNVs and dN/dS by gene and over the whole genome.

View this table:

Within-host variant sites are phylogenetically associated

We sought to gain a better understanding of SARS-CoV-2 evolution and determine whether iSNVs could be used to help resolve phylogenies and transmission clusters. For the 1390 genomes in our study, we constructed a phylogeny using the robust procedure outlined by (42) (Fig. 5A). Viral phylogenies are based on the consensus sequence for each sample, with branches indicating differences in the consensus sequence among samples. Given the inferred narrow transmission bottleneck size, we hypothesized that consensus changes on the phylogeny arise due to the emergence of within-host variants which either reach consensus within the individual they emerged in, or that fail to reach consensus but are then transmitted and result in a consensus change in the recipient. In a sufficiently densely sampled population of infected individuals we should therefore be able to observe a phylogenetic association between samples containing iSNVs with branches on the tree leading to a change in consensus at the same locus.

Fig. 5 Consensus phylogeny of all isolates.

In (A), tips are colored by sampling center (Oxford, orange, or Basingstoke, green). The tree scale is in substitutions per site. (B to D) Distribution of samples with iSNVs at three loci. The genomic coordinate (with respect to the Wuhan-Hu-1 reference sequence) appears in the top left. Tree branches are colored by the consensus base at that position, and filled circles indicate samples iSNVs present at minimum 3% for samples with depth of at least 100 at that position, and are colored by the most common minor variant present. For sites 28580 (B) and 20796 (C), an inset panel enlarges the section of the phylogeny where a consensus change is in close proximity to iSNVs with the relevant pair of nucleotides involved. The highlighted samples were prepared in separate batches and the patterns are not due to contamination. (D) Variants at site 21575 (L5F) occur in 14 samples, but with no phylogenetic association with consensus changes at this site. This may represent independent emergence of this variant in multiple individuals. The phylogeny was constructed by maximum likelihood according to the robust procedure outlined by Morel et al. (42).

Of the 563 high-confidence iSNV sites, we identified 153 iSNV sites present in at least two samples and where we also observe differences in the consensus among samples (SNPs). We call these sites iSNV-SNPs. We examined the proximity of tips with the iSNVs to the position of consensus changes (between the two most common bases at the site of the iSNV) on the phylogeny. A highly significant negative association (one-sided Mann-Whitney U-test, p<3x10−16; Fig. S6A) was found between the presence of an iSNV at a given site in a sample and the patristic distance to the nearest example of a consensus change at the same site; that is, intrahost variation clustered on the tree with branches supported by the same variant as consensus. When we tested sites where we had identified at least two iSNVs individually, six showed a significant association after Benjamini-Hochberg correction (p<0.05), reducing to five if only one sample from each individual was included. Repeating this procedure on each of 1000 phylogenetic bootstrap replicates yielded a universal very strong association when taking sites across the whole genome (maximum p=2.46x10−10), while every bootstrapped tree had between one and nine significant iSNV-SNPs (median 7, IQR 5-7).

In Fig. 5B we show the example of site 28580 (significant in 85.8% of bootstrap replicates), with the red clade representing change from the global consensus G to A (a nonsynonymous change D103N in N), and nearby iSNVs occurring, both as minor As in the nodes ancestral to the change branch, and as minor Gs in the branch’s immediate descendants. Based on corresponding epidemiological data, this represents a healthcare-associated cluster with onward transmission to close contacts. In Fig. 5C we give the further example of site 20796 (significant in 98.4% of bootstrap replicates), a synonymous substitution L6843 in ORF1a. Trees for the other significant sites after Benjamini-Hochberg correction appear in Fig. S7. Supporting this relationship between SNPs and iSNVs, we note that in the household transmission pairs we examined, for the five consensus differences where there was sufficient depth, all were within-host variant in one of the two individuals (Fig. 4B).

For the 261 iSNVs that are present in at least two individuals but never reach consensus, we analyzed the association with the phylogeny of each iSNV as a discrete trait, using two statistics: the association index (34) and the mean patristic distance between iSNV tips. After adjustment for multiple testing, no sites showed a p-value less than 0.05 for a phylogeny-iSNV association for either statistic. Similarly, if we simply compare the distance to the nearest iSNV tip amongst iSNV and non-iSNV tips across all 261 iSNV sites, there is also no evidence of phylogenetic association (one-sided Mann-Whitney U-test p~=1, Fig. S6B). Nevertheless, some individual sites do show patterns suggestive of iSNV transmission, with diversity maintained after transmission (22 with p<0.05 before adjustment for multiple testing for at least one of the two statistics; those 9 with p<0.025 are shown in Fig. S7) suggesting we may lack the power to statistically detect some associations. Among the 15 household transmission pairs, we observed only one iSNV shared in two individuals within the same household. This iSNV was unique to these two individuals in our dataset, demonstrating a likely example of transmitted viral diversity (Fig. 3B).

Taken together, our observations suggest the transmission bottleneck can be wide enough to permit co-transmission of multiple genotypes in some instances, but narrow enough that multiple variants do not persist after a small number of subsequent transmissions. In the cases where transmission culminates in a consensus change on the phylogeny, these patterns are readily observable, but in most cases we suggest patterns of co-transmission are drowned out by the high proportion of iSNVs that fail to transmit, or are transmitted but then lost. Analysis of transmission events over multiple generations is needed to fully elucidate these patterns.

Variants occurring repeatedly, but without phylogenetic association, could indicate sites under selection in distinct individuals (43). Of particular note are variants we observed at three sites in S: 21575 (L5F), 22899 (G446V) and 24198 (A879V), with G446V lying within the receptor binding domain (RBD). The minor variant F5 was observed in 14 samples, and represented SNPs in 8 samples, but did not have phylogenetic association in our iSNV-SNP analysis (p = 0.771 before multiple testing adjustment, Fig. 5D). This L5F mutation has been shown to increase infectivity in vitro (44), and has previously been identified as a potential site subject to selection (45). This variant has repeatedly been observed in global samples, including as minority variant, but appears to be increasing in frequency slowly if at all, suggesting it is only advantageous within a small subset of individuals, with the variant either ‘reverting’ in subsequent infections (as seen in HIV (46)), or failing to transmit at all. Similarly, we observed the minor variants V446 and V879 in 4 and 6 individuals respectively. Both variants have previously been shown to reduce sensitivity to convalescent sera in vitro (44), and V446 strongly reduces binding of one of the antibodies (REGN10987) in the REGN-Cov2 antibody cocktail (47), suggesting these may represent antibody escape mutations. We did not observe N501Y, or E484K, both mutations of concern, in any of our samples (48).

Concluding remarks

We uncovered a consistent and reproducible pattern of within-host SARS-CoV-2 diversity in a large dataset of over 1000 individuals, with iSNV sites showing strong phylogenetic clustering patterns if they are also associated with a change in the consensus variant at the same site. However, most samples harbored few intrahost variants, and estimated transmission bottleneck sizes were very small, with maximum likelihood estimates between 1 and 8 among household transmission pairs, meaning that if mutations do arise they will be prone to loss at the point of transmission. The dense sampling and deep sequencing of SARS-CoV-2 has enabled us to witness ‘evolution-in-action’, with variants generated in one individual, if transmitted, leading to a change in consensus and fixation in subsequently infected individuals. This suggests within-host variants could be used, at least in some instances, to help better resolve patterns of transmission in a background of low consensus diversity.

Our observations indicate the within-host emergence of vaccine- and therapeutic-escape mutations is likely to be relatively rare, at least during early infection when viral loads are high. However, even in the absence of vaccine- or therapeutic selection pressure, potential host-adaptive mutations are observable with sufficient frequency that even a rare transmission event, combined with narrow bottleneck size, could result in rapid spread. Here, we identified 30 nonsynonymous minor variants in Spike that are present in multiple individuals (Table S2). Two of these (G446V and A879V) have previously been shown to escape antibody binding (44), and a third, L5F, has been shown to increase viral infectivity (44). We suggest that commonly occurring iSNVs, along with variants known to affect transmissibility, severity of infection, or immune responses, should be investigated and monitored, particularly as vaccines and therapeutics are rolled out more widely.

The emergence of new variants of concern, including B.1.1.7, B1.351 and P.1, underscores the need for continued vigilance. A leading hypothesis is that these variants, characterized by a large number of nonsynonymous mutations, originated within individuals with long durations of infection where the virus was subject to prolonged immune pressure (7, 8), potentially facilitated by the within-host emergence of deletions (49). However, the presence of multiple mutations on the same genetic background is not a necessary prerequisite for a new variant to be worrying. The single D614G Spike mutation spread globally after it emerged during the early stages of the pandemic, likely due to a transmission advantage (50). The potential for mutations including N439K and E484K, which may enable the virus to evade host-immune responses (47, 51), to emerge on the highly-transmissible B.1.1.7 background is also troubling, particularly as population immunity builds due to natural infection and vaccination.

Our work demonstrates that an essential requirement for incorporating intrahost variants in any analysis is an understanding of the observed intrahost diversity in the context of the methods used to produce the deep sequencing data. Throughout, we aimed to minimize sequencing artefacts and sample contamination where possible. Moreover, our results emphasize the power of open data, large and rigorously controlled datasets, and the importance of integrating genomic, clinical, and epidemiological information, to gain in depth understanding of SARS-CoV-2 as the pandemic unfolds.

Materials and methods

RNA extraction

Residual RNA from COVID-19 RT-qPCR-based testing was obtained from Oxford University Hospitals (‘Oxford’), extracted on the QIASymphony platform with QIAsymphony DSP Virus/Pathogen Kit (QIAGEN), and from Basingstoke and North Hampshire Hospital (‘Basingstoke’), extracted with one of: Maxwell RSC Viral total nucleic acid kit (Promega); Reliaprep blood gDNA miniprep system (Promega); or Prepito NA body fluid kit (PerkinElmer). An internal extraction control was added to the lysis buffer prior to extraction to act as a control for extraction efficiency (genesig qRT-PCR kit, #Z-Path-2019-nCoV in Basingstoke, MS2 bacteriophage (52) in Oxford). The #Z-Path-2019-nCoV control is a linear, synthetic RNA target based on sequence from the rat ptprn2 gene, which has no sequence similarity with SARS-CoV-2 (GENESIG PrimerDesign pers. comm, 6 April 2020). The MS2 RNA likewise has no SARS-CoV-2 similarity (52). Neither control RNA interfered with sequencing.

Targeted metagenomic sequencing

Samples with suspected epidemiological linkage, where this information was available prior to sequencing, were processed in different batches. Sequencing libraries were constructed from remnant volume of nucleic acid after clinical testing, ranging from 5 to 45 μl (median 30μl) for each sample depending on the available amount of eluate. These volumes represented 1-15% of the original specimen (swab). Libraries were generated following the veSEQ protocol (16) with some modifications. Briefly, unique dual indexed (UDI) libraries for Illumina sequencing were constructed using the SMARTer Stranded Total RNA-Seq Kit v2—Pico Input Mammalian (Takara Bio USA, California, US) with no fragmentation of the RNA. An equal volume of library from each sample was pooled for capture. Size selection was performed on the captured pool to eliminate fragments shorter than 400nt, which otherwise may be preferentially amplified and sequenced. Target enrichment of SARS-CoV-2 libraries in the pool was obtained through a custom xGen Lockdown Probes panel (IDT, Coralville, USA), using the SeqCap EZ Accessory Kits v2 and SeqCap Hybridization and Wash Kit (Roche, Madison, US) for hybridization of the probes and removal of unbound DNA. Following 12 cycles of PCR for post-capture amplification, the final product was purified using Agencourt AMPure XP (Beckman Coulter, California, US). Sequencing was performed on the Illumina MiSeq (batches 1-2) or NovaSeq 6000 (batches 3-27) platform (Illumina, California, US) at the Oxford Genomics Centre (OGC), generating 150bp or 250bp paired-end reads.

Quantification controls

A dilution series of in vitro transcribed SARS-CoV-2 RNA (Twist Synthetic SARS-CoV-2 RNA Control 1 (MT007544.1), Twist Bioscience) was included in every capture pool of 90 samples starting from batch 3, and sequenced alongside the clinical samples. Control RNA was serially diluted into Universal Human Reference RNA (UHRR) to a final concentration of SARS-CoV-2 RNA of 500,000, 50,000, 5,000, 500, 100 and 0 copies/reaction. From this we produced a standard curve demonstrating linear association between viral load and read depth (Fig. S1). For an experiment comparing iSNV presence with and without probe capture, we additionally sequenced two replicates of the Twist RNA control without capture, diluted into UHRR to give an expected concentration of 50,000 copies per reaction.

As an additional validation step, we compared intrahost single-nucleotide variants (iSNVs) in re-sequenced controls with data for the stock RNA sequenced and provided by the manufacturer (Twist Bioscience). Six well-defined iSNVs, which were present in the manufacturer’s data and presumably arose during in vitro transcription, were also recovered by our protocol (Fig. S8). In addition, we identified 112 sites that appeared vulnerable to low-frequency intrahost variation in vitro (Table S3), possibly as a result of structural variation along the genome or interaction with the sequencing protocol. We blacklisted vulnerable sites from further analysis.

In-run controls

In addition to the synthetic RNA standards described above, each batch included a non-SARS-CoV-2 in-run control consisting of purified in vitro transcribed HIV RNA from clone p92BR025.8, obtained from the National Institute for Biological Standards and Control (NIBSC) (53). For batches 1 and 2, which were sequenced prior to synthetic RNA becoming available, we included negative buffer controls. As additional negative controls, we sequenced 6 matched clinical samples from non-COVID-19 patients, distributed across different sequencing runs; none contained any SARS-CoV-2 reads.

Minimizing risk of index misassignment

All samples had unique dual indexing (UDI) to prevent cross-detection of reads in the same pool. We used the in-run HIV RNA controls to estimate index misassignment, as this provided a sequence-distinct source of RNA: <3 SARS-CoV-2 reads were detected in any HIV control (median 0) and <10 HIV reads were detected in any SARS-CoV-2 control (median 0), suggesting that index misassignment, if present, occurred at extremely low levels.

Bioinformatics processing

Demultiplexed sequence read pairs were classified by Kraken v2 (54) using a custom database containing the human genome (GRCh38 build) and the full RefSeq set of bacterial and viral genomes (pulled May 2020). Sequences identified as either human or bacterial were removed using filter_keep_reads.py from the Castanet (17) workflow (55). Remaining reads, comprised of viral and unclassified reads, were trimmed in two stages: first to remove the random hexamer primers from the forward read and SMARTer TSO from the reverse read, and then to remove Illumina adapter sequences using Trimmomatic v0.36 (56), with the ILLUMINACLIP options set to “2:10:7:1:true MINLEN:80”. Trimmed reads were mapped to the SARS-CoV-2 RefSeq genome of isolate Wuhan-Hu-1 (NC_045512.2), using shiver (57) v1.5.7, with either smalt (58) or bowtie2 (59) as the mapper. Both mappers generated comparable results; smalt was used for the final analysis. Only properly paired reads with insert size under 2000 and with at least 70% sequence identity to the reference were retained. For analysis of consensus genomes, consensus calls required a minimum of 2 uniquely mapped (deduplicated) reads per position, equivalent to >15 raw reads per position. Analysis of within-host diversity was restricted only to positions with minimum raw depth of 100, except when examining diversity within presumed recipients of transmissions in the bottleneck analysis. Minor allele frequencies were computed at every position using shiver (57) (tools/AnalysePileup.py), with the default settings of no BAQ and maximum pileup depth of 1000000. Lineages were assigned by the Pangolin web server (60) using the determined consensus genome for each sequenced sample.

Alignment

Oxford and Basingstoke samples were selected if the consensus sequence (inferred from unique mapped reads) consisted of no more than 25% N characters. As an alignment to the reference sequence was already performed in shiver, no further alignment was necessary. To place these data into the global phylogenetic context and help resolve ancestry, a collection of non-UK consensus sequences from the GISAID database (61) were included in the set of sequences to be aligned. All GISAID (62) sequences were downloaded from the database on the 26th April 2020 and filtered to remove sequences that were less than 29800 base pairs in length, were more than 1% Ns, or were from the United Kingdom. The remaining sequences were clustered using CD-HIT-EST (63) using a similarity threshold of 0.995, and then one sequence per cluster picked. The resulting set, along with the reference genome Wuhan-Hu-1 (RefSeq ID NC_045512), were aligned using MAFFT (64), with some manual improvement of the algorithmic alignment and removal of problematic sequences performed as a post-processing step. Indels with respect to Wuhan-Hu-1 in both the Oxford/Basingstoke and GISAID alignments were deleted, resulting in two alignments of 29903 nucleotides that could be readily combined.

Demonstration of the effect of read downsampling

To demonstrate the effect of read depth on estimated iSNV counts, we selected the 30 samples with the highest total number of mapped reads, and, for each, chose a variety of downsampling fractions and removed all but that proportion of called bases from consideration. We then determined, for each sample and fraction, the number of iSNVs that would be identified at a threshold of 3% minor allele frequency at a minimum depth of 100, if only that fraction of called bases were available to us.

Transmission bottleneck analysis

Sixteen potential transmission pairs were identified by shared address (household) and first positive sample within two weeks. If samples from the two individuals in the household differed by fewer than three consensus differences (15 households), direct transmission was assumed. Apart from one genome position in household 6, and one in household 12, all sites associated with a consensus difference within a household were within-host variable in at least one member of the household pair, lending support to assumption of direct transmission (the exceptions are associated with low-read samples). Household 15 had six consensus differences, and was therefore excluded from our bottleneck analysis, although we note that for all six positions, the site was within-host variable in one or other individual. This pattern is inconsistent with direct transmission, but may represent transmission from a common source. Where the first samples for each individual in the household were more than one week apart, we assumed the earlier sampled individual was the donor, otherwise we considered both possible directions of transmission. Where individuals had more than one sample, or replicate sequences from the same sample, we used the sample/replicate with the highest number of mapped reads.

Bottleneck size was calculated using the exact beta-binomial method described in (28). Since most samples in the analysis had <50,000 mapped reads, we considered all sites in the genome, including sites in the 3′ and 5′ UTR, but excluding the poly-A tail (positions 29865-29903), the 18 “highly shared” sites, and those identified from the synthetic controls. All sites >3% MAF and more than 100 reads in the assumed source individual were used in the analysis. In the recipient all reads at these sites were considered, with an error threshold of 0.5% MAF. Following (28), 95% confidence intervals were calculated using a likelihood ratio test. No estimate was recorded for household 8 since there were no identified iSNVs >3% in the donor.

Calculation of dN/dS

The total number of synonymous and nonsynonymous substitutions in the SARS-CoV-2 genome was estimated using the first method of (65) applied to the coding regions of the Wuhan-Hu-1 reference sequence. Overlapping reading frames were accounted for such that a substitution was considered nonsynonymous overall if it was nonsynonymous in either frame.

We took two approaches to this calculation, first by counting all iSNVs individually, and second by counting only unique iSNVs. In the latter case, where we detected iSNVs with different base changes at the same position, we included only the most frequent. The results of the former are the basis for Table 1, while those of the latter appear in Table S5.

The dN/dS ratio for iSNVs over a genomic region G was then calculated as:pGipNTGN/pGipSTGS where ipN is the fraction of iSNVs at p that are nonsynonymous, or 0 if there are no iSNVs at p, TGN the total number of potential nonsynonymous substitutions in G, and the denominator replaces N with S to represent synonymous substitutions.95% confidence intervals for these estimates were obtained using the likelihood ratio test.

Phylogenetics

Phylogenetic reconstruction was performed on the alignment consisting of the 1390 consensus sequences, along with the GISAID set and the Wuhan-Hu-1 reference sequence. We followed the recommendations of Morel et al. (42) whereby 100 separate maximum likelihood phylogenies were generated using RAxML-NG (66) and the GTR+G substitution model, such that each reconstruction used a different random starting parsimony tree. The final phylogeny was then obtained from this set using majority rule. This final tree was rooted with respect to the reference sequence, and then that and all GISAID isolates were pruned.

To identify homoplasic sites, we selected sites that changed state more than once along the tree, after inferring the states at internal nodes using ancestral state reconstruction as implemented in ClonalFrameML (67) and rooting the tree using the reference genome NC_045512.

The recommendations of Morel et al. do not easily lend themselves to fast bootstrapping, so in order to explore phylogenetic uncertainty we performed an additional phylogenetic reconstruction on the same alignment analysis using the ultrafast bootstrap procedure in IQ-TREE (68). 1000 bootstrap replicates were used.

Phylogenetic association of iSNVs and SNPs

Where an iSNV corresponded to a consensus SNP (by the base pair involved, not simply the site), we performed ancestral state reconstruction on the consensus trees using ClonalFrameML (67) to identify all branches upon which that substitution was involved. Tips derived from the same clinical sample were then pruned until only one (the one with the highest overall depth) remained. We then, for each tip in the tree, calculated the patristic distance from that tip to the midpoint of the closest one of these branches, and used a one-tailed Mann-Whitney U-test to test for association between the iSNV existing in a sample and this distance. Multiple testing was controlled for using the Benjamini-Hochberg adjustment. As a sensitivity analysis, this was repeated such that all but one tip per infected individual, rather than per clinical sample, were pruned. These analyses were done both on an individual site level and across all sites of interest.

To confirm that the associations we observed here were unaffected by phylogenetic uncertainty, we used the set of 1000 IQ-TREE bootstraps. We repeated the Mann-Whitney U tests above for each of these 1000 trees.

Phylogenetic association of iSNVs at consensus invariant positions

For the remaining iSNVs, we calculated the extent of association with the consensus phylogeny by treating the presence of an iSNV as a discrete character and calculating the association index, and the mean patristic distance between iSNV tips. Once again the consensus tree was pruned such that tips corresponding to samples with read depth <100 at the position and all but one tip coming from the same individual were removed. A null distribution was generated by permuting the tip labels of this tree 10,000 times, and a one-sided permutation test p-value calculated. Multiple testing was adjusted for as above. In addition, for each tip in the phylogeny at each site of interest, we calculated the minimum patristic distance to a different tip corresponding to an iSNV, and used the Mann-Whitney U-test again to compare the distribution of these distances between iSNV and non-iSNV tips.

Supplementary Materials

science.sciencemag.org/cgi/content/full/science.abg0821/DC1

Figs. S1 to S8

Tables S1 to S5

List of OVSG members

List of COG-UK consortium names and affiliations

MDAR Reproducibility Checklist

https://creativecommons.org/licenses/by/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

References and Notes

Acknowledgments: We thank Robert Esnouf, Adam Huffman, and the BMRC Research Computing team for unfailing assistance with computational infrastructure. We also thank Benjamin Carpenter and James Docker for assistance in the laboratory, and Lorne Lonie, Maria Lopopolo, Chris Allen, John Broxholme, Angela Lee and the WHG high-throughput genomics team (Oxford Genomics Centre) for sequencing and quality control. The HIV clone p92BR025.8 was obtained through the Centre For AIDS Reagents from Drs Beatrice Hahn and Feng Gao, and the UNAIDS Virus Network (courtesy of the NIH AIDS Research and Reference Reagent Program). Ethics: The COVID-19 Genomics UK (COG-UK) consortium study protocol was approved by the Public Health England Research Ethics and Governance Group (reference: R&D NR0195) on the 8th of April 2020. For seasonal coronavirus screening, samples were collected with consent to assay for infectious causes of respiratory disease from patients admitted to Hampshire Hospitals NHS Foundation Trust to aid diagnosis and outbreak management and inform public health surveillance. Funding: We gratefully acknowledge the UK COVID-19 Genomics Consortium (COG UK) for funding. COG-UK is supported by funding from the Medical Research Council (MRC) part of UK Research & Innovation (UKRI), the National Institute of Health Research (NIHR) and Genome Research Limited, operating as the Wellcome Sanger Institute. The research was also supported by the Wellcome Core Award Grant Number 203141/Z/16/Z with additional funding from the NIHR Oxford Biomedical Research Centre. The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health. KAL and MAA were supported by The Wellcome Trust and The Royal Society (107652/Z/15/Z to KAL and 220171/Z/20/Z to MAA). MH, LF, MdC, GMC, NO, LAD, DB, CF and TG are supported by Li Ka Shing Foundation funding awarded to CF. PS is supported by a Wellcome Investigator Award (WT103767MA). JAT was supported by Wellcome Core Award (203141/Z/16/Z). CEM is supported by the Fleming Fund at the Department of Health and Social Care, UK, the Wellcome Trust (209142/Z/17/Z), and the Bill and Melinda Gates Foundation (OPP1176062). DWE is a Robertson Fellow and an NIHR Oxford BRC Senior Fellow. Author contributions: Conceptualization: TG, KAL, MH, LF, CF; Methodology: TG, KAL, MH, MdC, AT, D Bonsall; Software: TG, MH; Validation: TG, KAL, MH, PS; Formal analysis: MH, KAL; Investigation: KAL, MH, LF, TG, MdC, AT, GM-C; Resources: MA, ELW, NM, JL, SK, MM, RW, GV, AJ, NO, SMN, MAA, CEM, TEAP, DWE, RS, D Buck, AG, JAT, OVSG Analysis Group, PS, JS, SA, AdaSF, COGUK; Data Curation: DWE, NM, TG; Writing – original draft preparation: KAL, MH, TG; Writing – review and editing: KAL, MH, TG, CF, D Bonsall, LA-D, PS, JAT, COGUK; Visualization: MH, TG, KAL; Supervision: TG, D Bonsall, CF, ECT, TRC, JAT; Project administration: D Bonsall, TG, AT, AG, LA-D, D Buck; Funding acquisition: D Bonsall, COGUK, CF. Competing interests: DWE declares personal fees from Gilead outside the submitted work. MA is on the advisory board of Prenetics. All other authors declare no competing interests. Data and materials availability: All genomic data has been made publicly available as part of the COVID-19 Genomics UK (COG-UK) Consortium (69) via GISAID (62) and via the European Nucleotide Archive (ENA) study PRJEB37886. All other data are available in the main text, supplementary materials, or (70) (including the full alignment of consensus sequences and inferred tree). This work is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/. This license does not apply to figures/photos/artwork or other content included in the article that is credited to a third party; obtain authorization from the rights holder before using such material.

Stay Connected to Science

Navigate This Article