Tracing SARS-CoV-2 Clusters Across Local-scales Using Genomic Data

Detect Distinct Introduction Events and Identify Locally Circulating Clusters

Our dataset comprised 26,138 SARS-CoV-2 complete genomes, including 9,186 sampled from Houston and 16,952 contextual sequences. We conducted a discrete phylogeographic analysis on the time-calibrated phylogeny to identify distinct SARS-CoV-2 introduction events into Houston. This analysis revealed a total of 1,479 independent introduction events (95% highest posterior density [HPD]: 1,402 to 1,556). Notably, the sizes of resulting circulating clusters were highly skewed (Figure S1). The majority of introductions (909 events, 95% HPD: 853 to 968) resulted in singletons, while a few introductions led to clusters exceeding 2,000 cases. Temporal analysis further revealed that earlier introductions were more likely to result in larger clusters (Figure 2A). Specifically, during EPI Week 17, 55.6% of introductions resulted in clusters larger than 10, 33.3% in clusters smaller than 10, and 11.1% were singletons. By EPI Week 30, the distribution had shiftedsignificantly, with 86.2% resulting in singletons, 13.8% forming clusters smaller than 10, and no clusters larger than 10.

• Download figure
• Open in new tab

Figure 2. Dynamics of Viral Introduction into the Greater Houston.

A. Estimated weekly frequency of viral introductions, with curves colored by the size of the resulting local transmission lineages. Shading indicates the associated 95% Highest Posterior Density (HPD). B. Estimated weekly frequency of viral introductions, with curves colored to distinguish between domestic and international sources. Shading indicates the associated 95% HPD.

We classified viral imports into Houston based on their origin: domestic or international. At the onset of the outbreak, we observed scattered introductions from both sources. After late April, domestic importations rapidly increased and became the dominant type (Figure 2B). There were a total of 1,359 domestic imports (95% HPD: 1,279 to 1,432), significantly outnumbering the 119 international imports (95% HPD: 109 to 132). However, international imports peaked earlier and were associated with larger mean cluster sizes (Figures S2 and S3).

Subtrees extending from the introduction nodes were identified as locally circulating clusters. We selected 181 clusters, each containing five or more isolates, for further analysis of local dispersal. In total, these separate clusters included 7,657 sequences from the Greater Houston area, with the two largest clusters containing 2,198 and 2,031 sequences, respectively.

Phylogeny-Trait Correlation Among Locally Circulating Clusters

We explored the correlations between phylogenetic structures and demographic traits—specifically age and sex—to enhance our understanding of the factors influencing transmission dynamics. Our null hypothesis is that these traits are randomly associated with the phylogenetic structures. A low p-value (p < 0.05) refutes this hypothesis, indicating a strong correlation and suggesting limited viral dispersal between different traits.

We tested the association between 181 locally circulating clusters and demographic traits (Figure 3). In 20 clusters, the age group traits were tightly correlated with the phylogeny (p < 0.05). In contrast, only 6 clusters showed a similar tight correlation for sex traits. Generally, sex groups appear more interspersed within the phylogeny, suggesting that viral transmission is more constrained within age groups than between sex groups.

• Download figure
• Open in new tab

Figure 3. Phylogeny-Trait Correlation Among 181 Locally Circulating Clusters.

Blue circles showed the test result of association between locally circulating clusters and age groups. Green triangles showed the test result against sex groups. Shapes bellowed the red line (p < 0.05) indicated clusters that has tightly correlated with corresponding traits.

Demographic Determinants of Transmission

Locally circulating clusters encompassed 7,657 genomes distributed among diverse age groups: 2,412 young adults (ages 19–35), 2,342 middle-aged adults (ages 36–55), 995 infants and children (ages 0–12), 963 seniors (ages 56 and over), 919 teenagers (ages 13–18), and 26 individuals of unknown age. The sex distribution included 3,954 males, 3,687 females, and 16 individuals of unknown sex.

We jointly estimated a single discrete trait model to all circulating clusters to quantify viral dispersal among age and sex groups. We found young and middle-aged adults were identified as the primary drivers of viral transmission (Figure 4A). The eight most significant transitions included: from young to middle-aged adults at a rate of 4.775 transitions per year, and vice versa at 2.194 transitions per year; from middle-aged adults to infants and children at 1.946 transitions per year; to seniors at 1.567 transitions per year; to teenagers at 1.494 transitions per year; from young adults to seniors at 1.363 transitions per year; to teenagers at 1.245 transitions per year; and to infants and children at 1.174 transitions per year. Each of these transitions was strongly supported by a Bayes Factor exceeding 100. The detailed list for diffusion rate among age groups are available is Table S1.

• Download figure
• Open in new tab

Figure 4. Discrete Trait Diffusion Models of Age Group and Sex.

The age group model (A) and the sex model (B) are presented by circular charts. Chord thickness indicates the magnitude of the transition rate, while color signifies Bayes Factor support. Only transition rates supported by a Bayes Factor greater than 3 (BF > 3) from the discrete trait analysis are displayed. Transmission sources are on the left, transmission sinks are on the right.

In the sex-based model (Figure 4B), we observed that the transition rate from females to males (1.652) was higher than from males to females (0.437). The detailed list for diffusion rate among sex groups are available is Table S2.

Heterogeneous Dynamics of Viral Dispersal in Subregions of Greater Houston

We estimated a jointly fitted single discrete model to reconstruct the viral dispersal history in Greater Houston. This model categorized location traits into 29 groups, including 21 Independent School Districts in Harris County and 8 nearby counties. The joint model (Figure 5A) for these subregions revealed 82 transitions that were decisively supported by Bayes factors (>100). The transition from Fort Bend County to the Houston ISD had the highest transition rate, at 23.780. This was followed by transitions from Houston ISD to Cypress-Fairbanks ISD, with a rate of 10.203, and from Montgomery County to Houston ISD, with a rate of 6.923. The detailed list of rates can be found in Table S3.

• Download figure
• Open in new tab

Figure 5. Distinct Transmission Patterns in Subregions of Greater Houston.

(A) Discrete phylogeographic reconstruction of the dispersal history across 21 independent school districts in Harris County and 8 nearby counties. Subregions on the map are colored based on their associated Source Sink Scores. Arrow thickness indicates the magnitude of transition rates. All transitions shown on the map were decisively supported by Bayes factors (>100) (B) Source Sink Scores (SSS), (C) Local Import Scores (LIS), and (D) Persistence Time (PT) across these subareas. In these bar charts, golden error bars represent the associated 95% Highest Posterior Density (HPD), providing a measure of uncertainty for each score. SSS ranges from 1 (viral source) to -1 (viral sink). LIS ranges from 0 (epidemic is locally maintained) to 1 (epidemic relies on introduction).

We calculated the Source-Sink Scores (SSS) to identify populations as either viral sources or sinks, based on the net viral flow weighted by outbreak size. Using this metric, we ranked subregions from the most dominant sources to sinks (Figure 5B). We identified Houston ISD, with a Source-Sink Score of 0.629 (95% HPD: 0.549 to 0.708), Fort Bend County, with a score of 0.550 (95% HPD: 0.381 to 0.651), and Cypress-Fairbanks ISD, with a score of 0.069 (95% HPD: -0.020 to 0.171), as the primary sources for local dispersal in the Greater Houston area, where Source-Sink Scores were greater than 0. Further, we calculated Local Import Scores (LIS) to assess the relative influence of viral introductions versus local transmission in driving the epidemic (Figure 5C). Viral sources exhibited the lowest Local Import Scores, indicative of strong, locally sustained transmission. In contrast, viral sinks with higher Local Import Scores relied more on external introductions. Additionally, we estimated the median Persistence Time (PT) of viral transmission chains within each subregion (Figure 5D). The primary viral sources demonstrated the longest persistence times, showing more successful local transmission.

SARS-CoV-2 Genomic Dataset

With the support of the Houston Health Department (HHD), we accessed a large dataset of SARS-CoV-2 genomes sampled in Houston (>10,000), along with linked metadata, including ZIP code, age group, and sex. The first reported Delta variant case in Houston occurred in mid-April 2021. Our contextual dataset (non-Houston) was divided into two phases: Phase one included all worldwide sequences available in GISAID (www.gisaid.org) sampled before April 15, while Phase two sampled 1% of worldwide sequences available after April 15 (Figure S4). This sampling scheme balanced the need to include early sequences with the practical limits of handling a rapidly growing dataset.

The combined dataset of Delta sequences from Houston and contextual sources was then aligned to the reference genome (GenBank ID: NC_045512.2) using minimap2 v2.24 (36). We filtered out low-quality sequences with mapped completeness below 93% and trimmed alignments outside the reference coordinates 265:29674, padding with Ns to mask out the 3’ and 5’ UTRs. We then calculated the genetic base-pair differences between the alignments and the reference genome. For samples collected within the same Epi-Week, we excluded sequences with genetic differences greater than 3.0 standard deviations from the mean, aiming to preliminarily filter out those with a poor clock signal. In total, 26,138 alignments passed the filtering criteria. The contextual sequences were categorized as either domestic (excluding Houston) or international. Our dataset included 9,186 sequences from Houston, 5,334 from domestic sources, and 11,618 from international sources (Figure S5).

Time-Scaled Phylogenetic Inference

We inferred the time tree of our dataset in two steps. First, a maximum-likelihood phylogeny was estimated using IQ-TREE 2.3.2 (37) with the default settings. Subsequently, the resulting tree was time-calibrated using BEAST v1.10.5 (38). Given the size of our dataset, Thorney BEAST (13, 14) was applied to significantly reduce computational time by employing an alternative likelihood function. XML files for the BEAST runs were prepared using R scripts. Phylogeny file editing were performed using the tools jclusterfunk v0.0.25 (39) and gotree v0.4.5 (40). To minimize runtime by reducing burn-in, a roughly scaled time tree estimated by TreeTime 0.11.2 (41) was included in the XML file as the starting tree. We executed five chains of 2.5 billion states each, with a burn-in of 1 billion states. Trees were sampled every 7.5 million states, resulting in an empirical tree set size of 1,000. The convergence and mixing of all relevant parameters were inspected using Tracer 1.7 (42) to ensure that their associated effective sample size (ESS) values exceeded 200. One posterior tree (Figure S6) was then randomly selected and used as a fixed time tree for subsequent introduction analysis.

Introduction Analysis

We performed a preliminary discrete phylogeographic analysis on the fixed time tree obtained in the previous step to identify descent clusters likely arose from independent introductions. Tips of the tree were assigned as either ‘Houston’, ‘Domestic’, or ‘International’. XML files for BEAST runs were generated using R scripts. We executed five chains of 20 million states each, with a burn-in of 4 million states. Trees then were sampled every 160,000 states, resulting in an empirical tree set size of 500. The convergence and mixing aspects of all relevant parameters were inspected using Tracer 1.7 (42).

We considered an introduction event when a node was assigned the location ‘Houston’, while its parent node was labeled as non-Houston (either ‘Domestic’ or ‘International’). We defined locally circulating clusters as the subtree extending downward from the introduction node. Introduction time was identified as the midpoint on the branch connecting the introduction node and its parent. We summarized the weekly count of introduction events from the posterior tree sets, providing 95% highest posterior density estimates for these counts. Phylogeny reading and manipulation, including summarizing cluster size, were performed using ape 5.8 (43), treeio 1.20.2 (44) and ggtree 1.14.6 (45) packages.

We selected a representative tree (Figure S7) from the posterior tree set that matched the posterior median for total importations. From this tree, we extracted 181 distinct locally circulating clusters with five or more sequences. In each cluster, tips sampled outside Houston were pruned using the ape 5.8 package (43).

Tip-Trait Association Test

For each locally circulating cluster, we used the Association Index (29, 46) to quantify phylogeny-trait correlations. To assess the significance of these correlations, we generated null distributions of the Association Index by randomizing trait assignments on the tips 1,000 times, enabling us to perform the tip-trait association test. All scripts for this test were bundled into an R package named TTAT, which is publicly available on GitHub at https://github.com/leke-lyu/TTAT. This package takes various clades and trait data as input, using metrics such as the Association Index and parsimony score to perform the test. The p-values from the tests are displayed in a scatter plot of sex group vs. age group using ggplot2 (47).

Jointly Fitted Discrete Trait Model

Discrete trait analysis models the evolution of discrete states on a phylogeny, similar to sequence evolution. Assuming all circulating clusters shared underlying characteristics, we jointly estimated a single transition rate matrix (28, 48, 49) to describe transitions between traits—location, sex, and age groups—with the aim of reconstructing geographic dispersal patterns and identifying demographic determinants of transmission. XML files for BEAST runs were created using R scripts. For each trait, we ran five chains of 100 million states each, with a burn-in period of 20 million states. Trees were subsequently sampled every 800,000 states, ensuring each cluster achieved an empirical tree set size of 500. Convergence and mixing of all relevant parameters were inspected using Tracer.

Over half of the Greater Houston population resides in Harris County (50). In our phylogeographic model, we divided Harris County into 21 Independent School Districts (ISDs) and combined these with the eight surrounding counties, resulting in a total of 29 trait categories. When translating ZIP codes (the sampling location record) into these subregions, some ZIP codes spanned multiple ISDs; for instance, parts of one ZIP code might fall in both ISD M and ISD N. We treated these intersections as ambiguous traits, allowing the model to interpret the trait as either M or N. Additionally, 16 sequences lacked associated sex records, and 26 sequences lacked age data, which we also categorized as ambiguous traits.

Posterior Processing

The joint estimation procedure reconstructed the ancestral states for 181 clusters. Given these posterior tree sets, we estimate the following epidemiological metrics:

1. Source Sink Score (51): This metric, ranging between -1 and 1, measures net viral exports, weighted by outbreak size. A score approaching 1 suggests that the region primarily functions as a viral source. Conversely, a score nearing -1 indicates that the region predominantly acts as a viral sink.

2. Local Import Score (51): Ranging from 0 to 1, this metric measures the fraction of introductions relative to the total count of new cases in a region. A score near 1 indicates that importations predominate, while a score near 0 suggests that local transmissions dominate, indicating that the epidemic is primarily sustained locally.

3. Persistence Time: This metric measures how long a viral lineage circulates in a region. It is calculated by tracing the number of days it takes for a tip to move from its sampled location, going backward up the phylogeny until the node location differs from the tip location (34, 52).

Given the size of our dataset, we employed a divide and conquer strategy to efficiently manage the input data. The empirical tree sets for each local cluster were divided into separate files, with each file containing a tree corresponding to a unique state. Tree files shared the same state were read into ‘phylo’ objects using the treeio 1.20.2 package (44). These ‘phylo’ objects were then converted into structured data frames using the tidytree 0.4.6 package (53), facilitating the easier estimation of epidemiological metrics.