Genomic heterogeneity and clinical characterization of SARS-CoV-2 in Oregon

Cohort Overview

In order to assess the heterogeneity of SARS-CoV-2 genomes across OR, a total of 204 nasopharyngeal swab patient specimens (representing 188 unique patients) were sequenced from diverse clinical sites across OR. A total of 179 of these resulted in complete SARS-CoV-2 genome coverage (average-per-base coverage of 60x or greater) and were included for further study. For a subset of patients, virus from multiple longitudinal swab specimens were sequenced in order to determine whether any novel mutations emerge within an individual over the course of COVID-19. The demographics of the tested patients are reported in Table 1. Specifically, twelve of the patients died during the course of COVID-19, seven required intubation and 54 were admitted for acute or intensive care.

Genomic Diversity of the OR Cohort

Sequence results from the 179 OR specimens revealed a striking diversity of SARS-CoV-2 genomes across Oregon (Figure 1A). Virtually all major worldwide GISAID clades (Shu and McCauley, 2017) are represented in OR genomes, with the largest number in GH and G as well as sizable clusters in S and V-clades. Of note the former major clades have highest representation in Europe, but were also observed as the dominant clades in the spread to South America, Africa and Canada. The S-clade exhibits highest representation in Asia and V-clade represents an early introduction of the virus into Italy and Europe from China (Giovanetti et al., 2020). Genomic variation in these isolates spanned the entirety of the SARS-CoV-2 genome with notable hotspots in ORF1A, ORF3A and Spike (Figure 1B), and variant classes were predominantly missense or synonymous (Figure 1C). A cohort of isolates mapped to S-clade with high similarity to the original WA isolate (one of the earliest introduction events in the United States (Holshue et al., 2020)) as well as other early WA-specific sequences (Bedford et al., 2020). Given the extent of frequent travel between OR and WA, it is perhaps unsurprising that a large number of early isolates bear similarity to WA SARS-CoV-2 sequences. Overall, we observed a shift in clade representation over time, with the early cases dominated primarily by L-, S- and V-clades and later cases predominated by the G-clades (Figure 2A). Of particular note for outbreak tracing, we note that even within clades there is a striking amount of sequence diversity, with 55% of SARS-CoV-2 genomes containing at least one nucleotide variant that was completely unique among the cohort, and an additional 19% that shared an identical SARS-CoV-2 genome with only one other patient (Figure 2B). Only twelve percent of the genomes were common among four or more patients in the cohort.

• Download figure
• Open in new tab

Figure 1. Genomic diversity of SARS-CoV-2 in Oregon.

A) Unrooted phylogenetic tree of major SARS-CoV-2 clades with representative worldwide genomes subsampled from GISAID. OR isolates sequenced in this study are represented by circles. B) Genomic location and histogram of SARS-CoV-2 mutations across the OR isolates. C) Distribution of variant effects across the cohort.

• Download figure
• Open in new tab

Figure 2. SARS-CoV-2 genotype diversity over time in OR.

A) Histogram of clade representation by date for the Oregon SARS-CoV-2 genomes. B) Distribution of SARS-CoV-2 genomes that are unique to a single patient, shared among only two patients, shared among three and four or more patients. C) Phylogenetic tree showing distribution of D614 or G614 Spike genotype. D) Histogram of D614G status of Oregon SARS-CoV-2 genomes over time.

The identification of known functional variants in SARS-CoV-2 is of particular importance, especially given recent data showing increased infectivity and rapid spread of the D614G variant across the United States and worldwide. Our data show that D614G is already the dominant genotype in Oregon, with 69% of sequenced genomes coding the Glycine substitution (Figure 2C). Notably, in Oregon this variant emerged over time as it was present at a significantly higher proportion in the second half of the tested period (Mid-April to May) (p<1⨯10⁻⁸) (Figure 2D). We also tracked prevalence of D614G at specific subregions in OR based on clinic of origin (Figure 3). Notably, D614G prevalence varied dramatically by region, with sites near Hillsboro and Newberg exhibiting 50% or below prevalence of D614G and sites east of the Willamette river exhibiting a very high prevalence of D614G.

• Download figure
• Open in new tab

Figure 3. SARS-CoV-2 genotype diversity over time in OR.

Distribution of SARS-CoV-2 clade (inner pie distribution, see key) as well as D614 or G614 Spike genotype (outer pie distribution) at individual sample isolation sites across the Portland metro region. Number of samples is indicated in the center of each circle.

We also identified two rare mutations in the receptor binding motif (RBM) of the SARS-CoV-2 Spike protein (A475V and P463S). This segment of Spike mediates interaction with ACE2 (Lan et al., 2020; Shang et al., 2020; Wang et al., 2020; Yan et al., 2020). Structural studies indicate that Ala-475 of Spike engages in a hydrogen bond with Ser-19 of ACE2 (Lan et al., 2020; Shang et al., 2020; Wang et al., 2020). A second rare mutation in the RBM of Spike, P463S, is located in a highly-conserved residue that does not directly contact ACE2, but may be involved in maintaining the overall conformation of the RBM (Lan et al., 2020; Shang et al., 2020; Wang et al., 2020; Yan et al., 2020). Given the importance of Spike-ACE2 binding in viral cell entry, follow-up studies on the impact of these variants may be warranted. Mutations with the potential to alter the function of the RBM may be of relevance for infectivity, and might also impact the function of neutralizing antibodies recognizing Spike.

Associations Between SARS-CoV-2 Genotype and Patient Clinical Characteristics

Next, we sought to determine whether sequence variants in the isolated SARS-CoV-2 genomes were associated with differential clinical manifestation of COVID-19 in the patients. We looked at associations of SARS-CoV-2 variants with level of clinical care required; this included entry into an Acute Care Unit (ACU), Intensive Care Unit (ICU) or placement on a ventilator (Figure 4). This was done in comparison with other demographic and prior comorbidities as acute care predictors. A list of these characteristics along with a statistical summary across the different clinical categories is provided in Table 1. We note that despite the high prevalence of SARS-CoV-2 D614G in our study cohort, only one of twelve deceased patients had a D614G-positive infection (p<0.001, Fisher’s exact test). Despite this, D614G genotype did not well differentiate between overall admission status in the cohort. In contrast, as expected we see a strong differentiation between the ages of Outpatient vs ACU admitted patients (p =3.72e-10Mann-Whitney U Test), along with current or former smoking status (p = 7.75e-3 G-Test). Comorbidities were the strongest indicator of hospital admission, with all but asthma having a p ≤ 1.79e-3 (G-Test). The p-value for asthma was p = 9.63e-2 however we note the small sample size of this population. Finally, SARS-CoV-2 RT-PCR Ct values are overall higher in the admitted patient group (p = 2.97e-4 Mann-Whitney U Test), suggestive of an inverse relationship between nasal viral titer and disease severity. Applying unsupervised Uniform Manifold Approximation and Projection (UMAP) using the clinical and genomic variables revealed several distinct clusters, however none of these clearly separated outpatient vs admitted patients (Figure 5).

• Download figure
• Open in new tab

Figure 4. Venn diagram of clinical presentation of COVID-19 patients.

The relative frequency of patients requiring different levels of clinical care in the cohort. ACU=Acute care unit. ICU=Intensive care unit.

• Download figure
• Open in new tab

Figure 5. UMAP Projection of EHR and RT-PCR Data.

Clinical outcome of a patient is color-coded while D614G status is indicating by a + marker.

We then asked whether SARS-CoV-2 D614G or other genotypes exhibited utility in predicting hospital admission due to COVID-19. We developed a panel of classifiers based on different algorithmic approaches (logistic regression (LR), random forest (RF) and support vector machines (SVM)). For this we divided the sequenced patients into two populations based on admitted versus outpatient. Classifiers were then built based on patient demographics and prior comorbidities alone, demographics plus viral titer in the swab specimen (RT-PCR Ct value), and finally demographics plus viral titer plus D614G genotype. The results for these models based using LR, RF and SVM are detailed in Figure 6 below. Of note, classification of inpatient versus outpatient was highly accurate based on patient demographics alone using all three approaches, with both the RF and SVM achieving 88% accuracy. Notably, adding viral titer and D614G genotype to the model boosted the accuracy to over 90% for the RF approach, the best performing classifier out of all models tested. From this we conclude that demographic and clinical features have utility for predicting hospital admission in COVID-19 patients, with a modest improvement via the addition of SARS-CoV-2 sample characteristics.

• Download figure
• Open in new tab

Figure 6. Acute versus non-acute classifier using EHR and genomics.

A) Acute versus non-acute modeling using a logistic regression approach. B) Acute versus non-acute using random forest. C) Acute versus non-acute using a support vector machine approach. D) Receiver Operator Characteristic (ROC) curves for the three methods.

To further test possible correlations between genotype and patient outcomes, we assessed the predictive value of other prevalent SARS-CoV-2 variants identified in the cohort. For the remaining variants we calculated the Pearson Correlation, variants that were perfectly correlated with each other |R_p| = 1 were collapsed into a single statistical test (Figure 7). We then tested whether the variant discriminated between outpatients and inpatients using Fisher’s exact test and a logistic regression test (Table 2). We note two variants that show some significance across inpatients. The first variant, a T>C mutation at 3312 were significant at p<0.05 in all statistical tests (however we note the small sample size containing this variant). The second variant, an A>G mutation at 20268 shows >95% significance on Fisher’s Exact test with Benjamini-Hochberg, but has poor significance in both logistic regression and permutation tests. Overall, we see modest association with inpatient status across multiple variants; though larger cohort studies are needed to further evaluate this association.

• Download figure
• Open in new tab

Figure 7. Heatmap Showing the Pearson Correlation among variants that appear with a Frequency > 11 in Viral Genomic Samples.

Patients

All studies were conducted on discarded specimens and associated deidentified clinical information under IRB approval (study protocol 2020000127) with approved waiver of consent.

Specimens

Discarded nasopharyngeal swab specimens previously assessed as positive for SARS-CoV-2 by RT-PCR using either a Cobas SARS-CoV-2 RT-PCR assay on a Cobas 6800 system (Roche Diagnostics) or a BioGX SARS-CoV-2 assay on a BD MAX system (Becton Dickinson) were selected for sequencing. Specimens were suspended in 3 mL of Viral Transport Medium (VTM) (Hanks’ Balanced Salt Solution [HBSS] supplemented with 2% heat-inactivated fetal bovine serum [FBS], 100 µg/mL Gentamicin, and 0.5 µg/mL Amphotericin B), phosphate-buffered saline (PBS) or Copan Universal Transport Medium (UTM). Specimens were collected over a period of March to May 2020.

Nucleic Acid Purification

A mixture of total nucleic acids (NA) was isolated from 200 µL of swab-inoculated media using the QIAsymphony DSP Virus/Pathogen (Mini or Midi) kits and a QIAsymphony SP instrument (Qiagen). NA purification took place following the addition of carrier (poly-A) RNA and an internal control RNA, according to the manufacturer’s instructions. Purified NA were eluted in 85 µL of Elution Buffer AVE (Qiagen; RNase-free water with 0.04% NaN₃).

Sequencing Library Preparation

Sequencing libraries were prepared using Illuminia TruSeq RNA Library Prep for Enrichment kits (Illumina), and SuperScript II Reverse Transcriptase (ThermoFisher Scientific). Given the challenge of quantifying SARS-CoV-2 RNA within the purified NA, 8.5 µL of QIAsymphony eluate was used as starting material for library preparation. When samples possessed a SARS-CoV-2 threshold cycle (C_t) value of lower than 10.2, a 10-fold dilution of purified NA was performed prior to library preparation.

Hybrid Capture

Indexed cDNA libraries were pooled in groups of eight; each pool contained a total of 1.5 µg. Such pools were constituted of samples having similar SARS-CoV-2 C_t values (cobas 6800, Target 1, Roche). SARS-CoV-2 cDNAs were enriched using a hybrid-capture oligonucleotide panel covering the SARS-CoV-2 genome (GenBank Accession MN908947.3; 29,903 base pairs [bp]; Twist Bioscience) and a commercial hybrid-capture kit (Twist Bioscience). This approach yielded adequate (>60x average per base) coverage for clinical specimens with an initial RT-PCR cycle threshold of 30 or below (Figure S1).

DNA Sequencing

SARS-CoV-2-enriched libraries were sequenced on a NovaSeq 6000 or iSeq 100 instrument (Illumina) using paired-end reads of 150 bp each.

Sequencing Data Analysis

Raw sequencing reads were mapped against the reference sequence MN908947.3 using BWA-MEM v0.7.12-r1039 (Li and Durbin, 2009). Mapping statistics were calculated via samtools flagstat (samtools v1.10) (Li et al., 2009). Coverage was assessed by generating a coverage plot for each sample illustrating percent coverage across the genome via deepTools plotCoverage v3.4.3 (Ramirez et al., 2014). Samples with 60x or less coverage at the 90% percentile were excluded from downstream analysis. Average coverage for the 179 samples passing filter was approximately 3000x. Variant calling was performed using Freebayes v1.3.2-40-gcce27fc with standard filters and ploidy of 1. Variant calls were further filtered via bcftools vcffilter (bcftools v1.10.2-61) with parameters requiring a quality score > 3000, minimum read depth > 100, and minor allele frequency >= 0.5. Remaining variants passing filter were annotated via SnpEff v4.3 (Cingolani et al., 2012). The resulting variant calling file was used in conjunction with reference genome MN908947.3 to construct a consensus genome for each sample via bcftools consensus. Consensus sequences were manually assessed for variant call accuracy using Geneious Prime v2020.1.2 to confirm accuracy before submission to GISAID. SARS-CoV-2 genotype diversity was assessed by extracting variant calls from annotated VCFs using a custom Python script and generating plots using custom R scripts.

Variance Significance

Variants from Variant Call Files (vcfs) are mapped into an N_s × M indicator matrix is where N_s is the number of samples while M is the number of unique variants in the sample set. Initially, N_s = 181 and M = 233. We filter out all variants that do not appear in at least 11 samples. On the remaining variants we calculate the Pearson Correlation matrix and find any groups of variants that are perfected correlated or anticorrelated. From these groups we filter out all but the first variant. Significance (p values) is first calculated using Fisher’s Exact Test for two-tailed distributions. To correct for multiple-hypothesis tests we adjust the p values using the Benjamini-Hochberg correction with a family-wise error rate of 0.1. For further validation, we apply a permutation correction (with number of permutations= 20000) along with Fisher’s Exact Test to correct the p value. Finally, we calculate a final set of p values using a logistic regression model.

Clinical Prediction and Analysis

Patient’s age and Ct values have been normalized to be [0,1] by divind by either the max Ct value or dividing by 100 for the patient’s age. Current and former smoking status was merged into a single categorical variable. Finally, we removed any categorical variable whose total sample size < 5. For the UMAP projection we set N_neighbors = 15, min_distance = 0.1, and used the Euclidean distance metric. All predictive models (LR, SVM, RF) where trained using the Scikit-Image library.