Molecular Architecture of Early Dissemination and Evolution of the SARS-CoV-2 Virus in Metropolitan Houston, Texas

We sequenced the genomes of 320 SARS-CoV-2 strains from COVID-19 patients in metropolitan Houston, Texas, an ethnically diverse region with seven million residents. These genomes were from the viruses causing infections in the earliest recognized phase of the pandemic affecting Houston. Substantial viral genomic diversity was identified, which we interpret to mean that the virus was introduced into Houston many times independently by individuals who had traveled from different parts of the country and the world. The majority of viruses are apparent progeny of strains derived from Europe and Asia. We found no significant evidence of more virulent viral types, stressing the linkage between severe disease, underlying medical conditions, and perhaps host genetics. We discovered a signal of selection acting on the spike protein, the primary target of massive vaccine efforts worldwide. The data provide a critical resource for assessing virus evolution, the origin of new outbreaks, and the effect of host immune response. Significance COVID-19, the disease caused by the SARS-CoV-2 virus, is a global pandemic. To better understand the first phase of virus spread in metropolitan Houston, Texas, we sequenced the genomes of 320 SARS-CoV-2 strains recovered from COVID-19 patients early in the Houston viral arc. We identified no evidence that a particular strain or its progeny causes more severe disease, underscoring the connection between severe disease, underlying health conditions, and host genetics. Some amino acid replacements in the spike protein suggest positive immune selection is at work in shaping variation in this protein. Our analysis traces the early molecular architecture of SARS-CoV-2 in Houston, and will help us to understand the origin and trajectory of future infection spikes.


[Abstract]
We sequenced the genomes of 320 SARS-CoV-2 strains from COVID-19 patients in metropolitan Houston, Texas, an ethnically diverse region with seven million residents. These genomes were from the viruses causing infections in the earliest recognized phase of the pandemic affecting Houston. Substantial viral genomic diversity was identified, which we interpret to mean that the virus was introduced into Houston many times independently by individuals who had traveled from different parts of the country and the world. The majority of viruses are apparent progeny of strains derived from Europe and Asia. We found no significant evidence of more virulent viral types, stressing the linkage between severe disease, underlying medical conditions, and perhaps host genetics. We discovered a signal of selection acting on the spike protein, the primary target of massive vaccine efforts worldwide. The data provide a critical resource for assessing virus evolution, the origin of new outbreaks, and the effect of host immune response.
To better understand the first phase of virus spread in metropolitan Houston, Texas, we sequenced the genomes of 320 SARS-CoV-2 strains recovered from COVID-19 patients early in the Houston viral arc. We identified no evidence that a particular strain or its progeny causes more severe disease, underscoring the connection between severe disease, underlying health conditions, and host genetics. Some amino acid replacements in the spike protein suggest positive immune selection is at work in shaping variation in this protein. Our analysis traces the early molecular architecture of SARS-CoV-2 in Houston, and will help us to understand the origin and trajectory of future infection spikes. [Introduction] Pandemic disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus is now responsible for massive human morbidity and mortality worldwide (1-4). The virus was first documented to cause severe respiratory infections in Wuhan, China, beginning in late December 2019 (5)(6)(7).
Global dissemination occurred extremely rapidly, and has affected major population centers on many continents, especially in Asia, Europe, and North America (8)(9)(10). In the United States, Seattle and the New York City region have been especially important centers of COVID-19 disease caused by SARS-CoV-2. Similarly, in Seattle and King County, 4,620 positive patients and 303 deaths have been reported as of April 15, 2020 (12).
The Houston metropolitan area is the fourth largest and the most ethnically diverse city in the United States, with a population of approximately 7 million (13). The 2,000-bed Houston Methodist Health System has eight hospitals and serves a large multinational and socioeconomically diverse patient population throughout greater Houston. Although the city is well-known as the energy capital of the world, Houston also has a very large port, and the George Bush International Airport is a major transportation hub for domestic and international flights to Asia, Europe, and Central and South America. Many of these flights provide direct city-to-city service to diverse countries, including global population centers. The first COVID-19 case in metropolitan Houston was reported on March 5, 2020 (14). Community spread was suspected of occurring one week later (14). Many of the first cases in our region were associated with national or international travel, including areas known to have COVID-19 virus outbreaks (14). These facts, coupled with the existence of a central molecular diagnostic laboratory that serves all Houston Methodist hospitals and our very early adoption of a molecular test for the SARS-CoV-2 virus, permitted us to rapidly interrogate genomic variation among strains causing infections in the greater Houston area.
We here report that SARS-CoV-2 was introduced to the Houston metropolitan region many independent times from diverse geographic regions, including Europe, Asia, and South America. The virus spread rapidly and caused disease throughout the metropolitan region. We identified clear genomic signals of person-to-person transmission. Some events were known or inferred based on conventional public health maneuvers, but many were not. In addition, spatialtemporal mapping found evidence of rapid and widespread community dissemination soon after COVID-19 cases were reported in Houston. Analysis of the relationship between distinct genomic viral clades and hospitalization did not reveal significant evidence of more virulent genome types, underscoring the need for a better understanding of the relationship between severe disease, underlying medical conditions, gender, and host genetics.

Results
Description of metropolitan Houston. Houston, Texas, is located in the southwestern United States, 50 miles inland from the Gulf of Mexico. It is the most ethnically diverse city in the United States (15 Viral genome sequencing and phylogenetic analysis. We sequenced the genomes of SARS-CoV-2 strains dating to the earliest stages of confirmed COVID-19 disease in Houston. Phylogenetic analysis identified the presence of many diverse viral genomes that in the aggregate represent many of the major clades identified to date (16) (Fig. 2). Clades A2a, B, and B1 were the three more abundantly represented phylogenetic groups (Fig. 2). Geospatial and time series. We examined the spatial and temporal mapping of the genome data to investigate community spread (Fig. 3). The figure illustrates evidence of widespread and rapid community dissemination soon after the initial COVID-19 cases were reported in Houston. There is also evidence to suggest there were multiple independent strains introduced into metropolitan Houston, followed by local spread throughout all regions of the community.
Epidemiologically linked patients. We investigated the relationships among some of the genomes that were obtained from patients known to share common epidemiologic associations, such as living in the same household. In all instances, individuals known to be epidemiologically associated had identical or nearly-identical SARS-CoV-2 genomes, a finding consistent with person-toperson transmission of the virus.
Geospatial relationships of genetically similar isolates. We next tested the hypothesis that genetically related viruses were constrained to particular geographic areas of metropolitan Houston. Although in some instances this was the case, an important observation was that most of the individual related subclades were comprised of strains distributed over broad geographic areas (Fig. 4). These findings are consistent with the known propensity of respiratory virus SARS-CoV-2 to spread rapidly from person to person. Patients with shared viral genomes likely constitute epidemiologic clusters, as a consequence of direct transmission to one another, shared transmission through an unknown third party, or via a reticulate network.

Relationship between viral clades, clinical characteristics of infected
patients, and other metadata. It is possible that SARS-CoV-2 genetic subtypes may have different clinical characteristics, analogous to what was thought to have occurred with the Ebola virus (17)(18)(19) and is known for other infectious agents. As an initial examination of this issue in SARS-CoV-2, we tested the hypothesis that patients with disease severe enough to warrant hospitalization were infected with a non-random subset of virus genotypes. We also tested the hypothesis of non-random association between virus clades and disease severity based on in-patient versus out-patient status, the need for mechanical ventilation, and the number of days on a ventilator. We found no apparent simple relationship between viral clades and disease severity using these indicators of disease severity. Similarly, there was no simple relationship between viral clades and other metadata, such as gender, age, or length of stay ( Fig. S1-S7).
Machine learning analysis. Machine learning models were built to predict mortality, length of stay, in-patient status, or ICU admission based on the viral genome sequence. F1 scores (the evaluation metric used in classification algorithms) for most classifiers ranged between 0.5 to 0.6, indicative of classifiers that are performing similarly to random chance. Outcome (lived versus died) was only weakly correlated with age in this data set (PCC = 0.318), and similarly regression models built to predict age and length of stay based on the viral genome sequence also had poor performance, with R 2 values near zero.
Classifiers were also trained to predict host characteristic, gender and ethnic group. The two largest ethnic groups in the patient population were African American and Caucasian, as recorded in the electronic medical record. The binary classification model had an F1 score of 0.67 [0.51-0.83, 95% CI], likely reflecting social networks in early person-to-person transmission. A table of classifier accuracy scores and performance information is provided as Table S1.
Analysis of the nsp12 polymerase gene. The SARS-CoV-2 genome encodes an RNA dependent RNA polymerase (RdRp) used to replicate this RNA virus.
Two amino acid substitutions (Phe476Leu and Val553Leu) in the nsp12 RdRp have been reported to confer significant resistance in vitro to remdesivir, an adenosine analog (20). Remdesivir is inserted into RNA chains by RdRp during replication, resulting in premature termination and inhibition of virus. This compound has shown prophylactic and therapeutic benefit against MERS-CoV experimental infection in rhesus macaques (21). A recent publication describing results from a compassionate use protocol reported that remdesivir may have therapeutic benefit in some hospitalized patients (22). If efficacy is confirmed by a randomized controlled trial, this drug may be be widely used in large numbers of patients worldwide.
To acquire baseline data about allelic variation in the gene encoding nsp12 RdRp, we analyzed the 320 sequenced viral genomes. The analysis identified 12 nonsynonymous single nucleotide polymorphisms (SNPs) in nsp12, resulting in 11 amino acid replacements throughout the protein ( Table 1). The most common amino acid change was Pro322Leu, identified in 224 of the 320 Houston isolates. This amino acid replacement is present in genomes from the A2a clade, and distinguishes the A2a clade from other SARS-CoV-2 clades. The other amino acid changes in nsp12 were mainly present in single isolates from individual Houston strains, and some have been identified in other strains in the global GISAID collection (23,24). A prominent exception was the identification of 29 strains with a Met600Ile polymorphism. All 29 strains were phylogenetically very closely related members of the A2a clade, and they all also had the P322L amino acid replacement characteristic of this clade (Fig. S8). These data indicate that the Met600Ile change is the derived state among strains with the Pro322Leu replacement. Importantly, none of the observed amino acid polymorphisms in nsp12 were located precisely at the two sites associated with remdesivir in vitro resistance (20), and 10 of the 11 polymorphic amino acids are located on the surface of this RdRp. However, importantly, the Ala448Val replacement occurs in an amino acid that is located directly above the nucleotide substrate binding site that is comprised of Lys545, Arg553, and Arg555, as recently shown by structural studies (25,26) (Fig. 5). The Ala448 position is comparable to Val553 (Fig. 5), and a Val553Leu mutation in SARS-CoV was identified to confer resistance to remdesivir (20).

Analysis of the gene encoding the spike protein. SARS-CoV-2 and related
coronavirus such as SARS-CoV and MERS gain entry into susceptible host cells using a densely glycosylated viral-surface molecule knows as spike (S) protein.
The S protein of SARS-CoV-2 virus and its close coronavirus relatives binds directly to the host angiotensin-converting enzyme 2 (ACE2) to enter host cells.
Thus, S protein is a major translational research target, including small molecule inhibitors and extensive vaccine efforts globally (27,28). Analysis of the gene encoding the S protein identified 41 SNPs, including 25 that produce amino acid changes (Table 2, Fig. 6A). Seven of these replacements (A263V, F456L, S967R, T1027I, M1050V, K1157M, and Q1208H) are not represented in the GISAID database as of April 15, 2020. F456 makes contact with the human ACE2 receptor, occupying a pocket formed by ACE2 residues Thr27, Asp30 and Lys31 ( Fig. 6B) (29). The F456L substitution is expected to detrimentally affect binding to ACE2, as the shorter Leu side chain would not fill the pocket well. We note that seven of the amino acid replacements map to the periphery of the S1 subunit N-terminal domain (NTD). Four of these replacements (A27S, T29I, G261V, and A263V) can be mapped to the recently determined cryo-EM structure of the SARS-CoV-2 spike (27), whereas three other amino acid changes (L18F, T22I, and H69Y) are located in flexible NTD loops that could not be modeled (Fig. 6A). The clustering of amino acid replacements to a distinct region of the protein, together with the occurrence of two different amino acid replacements occurring at the same residue is a strong signal of positive selection. Inasmuch as infected patients make antibodies against the NTD region, we favor the idea that host immune selection is a key force contributing to amino acid variation in this region, resulting in an enhanced fitness phenotype of virus variants. Also of note, the D614G replacement was observed in 70% (224 of 320 strains) of the strains sequenced in this study. Residue D614 is located in subdomain 2 (SD-2) of the S protein, and it forms a hydrogen bond and electrostatic interaction with two residues in the S2 subunit of a neighboring protomer (Fig. 6C). Replacement of aspartate with glycine would eliminate both interactions, thereby substantively weakening the contact between the S1 and S2 subunits. We speculate that this weakening produces a more fusogenic spike protein, because S1 must first dissociate from S2 before the S2 subunit can refold and mediate fusion of viral and cell membranes. Stated another way, virus strains with the D614G variant may be better able to enter host cells, potentially resulting in enhanced spread.

Discussion
We discovered substantial genetic diversity of the SARS-CoV-2 viruses causing COVID-19 disease in Houston, Texas. The majority of cases studied are progeny of strains that cause widespread disease in European and Asian countries. We interpret these data as demonstrating that the SARS-CoV-2 virus was introduced into Houston many times independently, likely by individuals who had traveled to different parts of the world. In support of this interpretation, the first cases in metropolitan Houston were associated with a travel history to a known COVID-19 region (14). The data are consistent with the fact that Houston is a large international city characterized by a multi-ethnic population.
The diversity present in our 320 viral genomes contrasts somewhat with data reported recently by Gonzalez-Reiche et al. (30), who studied 84 SARS-CoV-2 isolates causing disease in patients in the New York City region. Those investigators concluded that the vast majority of disease was caused by progeny of strains imported from Europe. Similarly, Bedford et al. (31) reported that much of the COVID-19 disease in the Seattle, Washington area was caused by strains that are progeny of a virus strain recently introduced from China.
The viral diversity present in metropolitan Houston permitted us to test the hypothesis that distinct viral clades were nonrandomly associated with hospitalized COVID-19 patients or disease severity. Although our sample size is relatively small, we found little evidence to support this hypothesis. Clearly, this important matter warrants further study with a larger sample size, and this analysis is currently underway.
We used machine learning classifiers in an attempt to identify SNPs that contribute to increased infection severity or otherwise affect virus-host outcome.
The models could not be trained to accurately predict these outcomes from the viral genome sequence data, which may be due to the small sample size and class imbalance. As such, no particular SNPs were identified that are predictive of disease severity or infection outcome. Classifiers were also trained to search for predictors of host characteristics such as age by decade, gender, and ethnic group. We found that the African American versus Caucasian populations had a predictive signal, and hence potentially significant SNPs, which may be borne out with increased sampling in these two groups. However, examination of the geographic distribution of the viral isolates classified using this model largely coincided with the demographic distribution of ethnic groups in the Houston metropolitan region. As such, the underlying SNPs found by the model may reflect social networks present early in the spread of SARS-CoV-2 in the Houston metropolitan area, rather than distinct viral or human genetic factors.
Remdesivir is a nucleoside analog that has been reported to have activity against MERS-CoV, a coronavirus related to SARS-CoV-2 (21). Recently, Grein et al. (22) reported that this drug shows promise in treating COVID-19 patients in a relatively small compassionate use protocol. Because in vitro resistance of SARS-CoV to remdesivir has been reported to be caused by either of two amino acid replacements in RdRp (Phe476Leu and Val553Leu), we interrogated our data for polymorphisms in the nsp12 gene. Although we identified 11 different inferred amino acid replacements in the 320 genomes analyzed, none of these were located at the two positions associated with resistance in vitro. These findings suggest that if remdesivir proves to be efficacious in COVID-19 patients, and is deployed widely as a treatment, the majority of SARS-CoV-2 strains currently circulating should be susceptible to this drug. However, the Ala448Val polymorphism we identified occurs at an amino acid site that intriguingly is located directly above the nucleotide substrate entry channel and nucleotide binding residues Lys545, Arg553, and Arg555 (20, 25) (Fig. 5). One possibility is that substitution of the smaller alanine residue with the bulkier valine may impose structural constraints for the modified nucleotide analog to bind and thereby disfavor remdesivir binding. This in turn leads to reduced incorporation of remdesivir into the nascent RNA, increased fidelity of RNA synthesis, and thus drug resistance. A similar mechanism has been proposed for a V553L change (20). We also identified one strain with a Lys477Asn replacement in nsp12. This substitution is located very close to a Phe476Leu replacement reported to produce partial resistance to remdesivir in vitro in SARS-CoV patients from 2004, although the amino acid positions are numbered differently in SARS-CoV and SARS-CoV-2. Structural studies have suggested that this amino acid is surfaceexposed, and distant from known key functional elements. Our observed Lys477Asn change is also located in a conserved motif described as a finger domain of RdRp (Fig. 5). One speculative possibility is that Lys477 is involved in binding a yet unidentified cofactor such as nsp7 or nsp8, an interaction that could modify nucleotide binding and/or fidelity at a distance. These data warrant additional study in larger patient cohorts, especially individuals treated with remdesivir.
Analysis of the gene encoding the spike protein identified many new inferred amino acid replacements not present in available databases. These data, coupled with structural information available for S protein, raises the possibility that many of the amino acid variants have functional consequences.
For example, clustering of several amino acid changes to the NTD suggests varying residues in this region bestow a fitness phenotype. Regardless, we are now beginning to acquire critical information about the location and extent of amino acid replacements occurring in the S protein in natural populations of SARS-CoV-2. These data permitted generation of many biomedically relevant hypotheses now under study.

Concluding statement
Our work represents analysis of the largest sample of SARS-CoV-2 genome sequences to date from patients in the southern United States. This investigation was facilitated by the fact that we had rapidly assessed the suitability and performance of the SARS-CoV-2 molecular diagnostic test in January 2020, more than a month before the first COVID-19 patient was diagnosed in Houston. Our large healthcare system has eight hospitals and many clinics located in geographically diverse areas of the city. These facilities serve patients of very diverse ethnicities and socioeconomic status, thus our data likely represent a broad overview of virus diversity causing infections throughout metropolitan Houston. We acknowledge that not every "twig" of the SARS-CoV-2 evolutionary tree in Houston is represented in these data because the samples studied were not comprehensive.
The genomes reported here represent the first strains documented to cause COVID-19 disease in the Houston area. Thus, they are an important resource that will underpin further and our ongoing study of SARS-CoV-2 molecular evolution and dissemination in Houston. As of April 15, 2020, there were 5,602 reported cases of COVID-19 in metropolitan Houston, with the number of cases growing daily. We are now sequencing virus genomes essentially in real time from infected patients, which will permit us to provide important data that can be exploited to facilitate an enhanced public health response to this pandemic. SARS-CoV-2 genome sequence analysis. Consensus viral genome sequences from the Houston isolates were generated using the ARTIC nCoV-2019 bioinformatics pipeline. Publicly available genomes and metadata were acquired through GISAID on April 13, 2020 (16). GISAID sequences containing greater than 1% N characters, and Houston sequences with greater than 5% N characters were removed from consideration. Identical GISAID sequences originating from the same geographic location with the same collection date were also removed from consideration to reduce redundancy. Nucleotide sequence alignments for the combined Houston and GISAID strains were generated using MAFFT version 7.130b with default parameters (32). Sequences were manually curated in JalView (33) to trim the ends and to remove sequences containing spurious inserts. Phylogenetic trees were generated using FastTree with the generalized time-reversible model for nucleotide sequences (34). Trees were parsed, annotated, and visualized using the R packages treeio and ggtree (35)(36)(37). CLC Genomics Workbench (QIAGEN) was used to generate the trees in the BI Filled map visualization. The home address for patients whose isolates were sequenced was matched to a dictionary of addresses downloaded from OpenAddresses.io (https://openaddresses.io/) to obtain the latitude and longitude geocoding data. Because addresses are transcribed and subject to manual error, the fuzzywuzzyR package was used to match the best address in the dictionary.
All address matches were manually reviewed for accuracy. The latitude/longitude coordinates were plotted onto a map using the Microsoft Power BI Desktop Map visualization.
To examine geographic relatedness among genetically similar isolates, the geospatial maps were filtered to four small independent clades. The home address of the patients was again visualized using the Microsoft Power BI Desktop Map visualization.
Machine learning analysis. Machine learning models were trained to predict patient metadata categories including mortality, length of stay, inpatient versus outpatient status, ICU admission, overall outcome, gender, age, and ethnicity from viral sequence data. Models were trained with the consensus whole genome FASTA files by dividing each viral genome into 15-mer oligonucleotide k-mers, which were used as features to train XGBoost models (38) as described previously (39,40). Models were assessed by computing F1 scores for classifiers and R 2 scores for regression models. Unless otherwise stated, data for the first 5folds of a 10-fold cross validation are shown.