How much should we sequence? An analysis of the Swiss SARS-CoV-2 surveillance effort

A total of 143,260 sequences were available for a population of 8.7 million inhabitants over the observed period in our data set (cf. Materials and Methods). The median percentage of cases that were sequenced each week was 7% (Fig. 1). In total, this constitutes 3.6% of all cases (n = 3.95 million) for the whole study period.

The first year of the pandemic was characterized by the emergence and persistence of multiple different viral lineages that led to an increase in nucleotide diversity and the lineage diversity index (LDI, as measured by Shannon diversity index, cf. Materials and Methods) (Fig. 1C; Fig. S2 to S6). This diversity on the nucleotide and lineage level was reduced with the first VOC, Alpha, rising to dominance (Fig. S3). All subsequent waves were dominated by a single VOC, with a slow rise in genetic diversity over time and a peak when a new, usually more divergent VOC was introduced (Fig. S3A). In contrast to Alpha, the Pango nomenclature accommodated the diversification of Delta and Omicron by introducing many more sublineages. This leads to a high LDI when all Pango-designated sublineages are considered, but an overall low LDI when these sublineages are consolidated to VOC labels (Fig. S3B). Moreover, as indicated by the arrows in Fig. 1C, Alpha and Delta were already present in Switzerland several weeks before they increased in frequency and truly emerged as the dominant variant. For Omicron, however, this initiation period was much shorter. The pandemic context during the emergence of each VOC was also very different. For example, case numbers were high at the time of introduction of Alpha and Omicron, but low for Delta (Fig. S4). Delta’s rapid growth phase started briefly after a period of decreasing Alpha cases, which coincided with relaxed health and safety measures in Switzerland. This—in addition to its intrinsically higher transmissibility (15)—likely played a role in its rapid increase at the end of June 2021 (calendar weeks 2021-25 onward in Fig. 1C).

To investigate how the observed LDI would have changed with lower sequencing intensity, we downsampled our data set (cf. Materials and Methods).

Naturally, downsampling had an influence on the absolute number of lineages detected (Fig. S6). However, downsampling had only a minimal effect on the LDI during VOC-dominated periods, especially for Alpha and Delta (Fig. S5). The impact of downsampling was stronger during the initial period of the pandemic when the number of available sequences was low. There was also a stronger effect during the later Omicron wave (i.e., once Delta had stopped circulating widely in Switzerland) when downsampling to less than a third of the original effort (less than 50k sequences).

We investigated how the first detection of the VOCs would have changed with a reduced sequencing effort, simulated by downsampling our data set. As expected, the delay of detection increased when fewer sequences were available (Fig. 2A). Surprisingly, downsampling to around a third of the available sequences (50k sequences) only led to a 1-day median delay for Delta and Omicron and no delay for Alpha. In general, we note that Delta was more sensitive to downsampling than Alpha and Omicron, as indicated by the longer tail of the delays upon repeated downsampling.

We reasoned that differences in the impact of downsampling might be due to differences in how steeply each VOC rises shortly after introduction. Thus, we further investigated the speed of introduction of each VOC, defined as how quickly a newly introduced VOC was established and grew in frequency. While the speed of introduction was similar for Alpha and Delta, it was much faster for Omicron (Fig. 3A). We observed that the speed of introduction also depended strongly on downsampling size, used as a proxy for sequencing effort. A higher sequencing effort was most accurate in capturing the respective growth of the emerging VOCs. Interestingly, the estimates of the downsampled data sets converged with increased sampling size toward the value calculated for the complete data set (red line in Fig. 3), suggesting that the actual sequencing effort captured true VOC dynamics accurately.

We then examined how long it takes for a VOC to reach dominance or 50% prevalence (Fig. 3B). In the full data set, Alpha and Delta required 88 and 103 days, respectively, while Omicron required only 31 days. With increased downsampling, the time needed to reach 50% prevalence decreased. In other words, lineages were seen as dominant earlier when fewer sequences were available. Delta was very sensitive for this measure as the spread of data and also the median time in reaching prevalence strongly increased with each downsampling. However, downsampling to around a third (50k sequences) accelerated prevalence only by 8 days median (4 days for Alpha and 1 day for Omicron), recapitulating the original data still well.

To summarize, we observed that the impact of downsampling on surveillance outcomes was very VOC-dependent (Fig. S7). Downsampling affected Alpha least with regard to the delay in detection and speed of introduction, whereas Delta was affected for both measures. For Omicron, the effect of downsampling on its detection was insignificant, while it was very strong for observing its growth.

These observations mean that differences in the impact of downsampling were not due to differences in growth during the introductory period, as Alpha and Delta displayed a similar speed of introduction for the complete data set but were differently impacted by downsampling. We could observe, however, that Delta had a unique feature. Its introductory period covered the phase in mid-June 2021 (around calendar week 2021-24) when cases (and sequencing) were very low (Fig. 1) and health and safety measures relaxed, giving Delta an additional advantage and leading to its rapid rise. This “gap” might be the reason why different surveillance measures were so strongly affected by downsampling for Delta.

To understand how the transmission information changes when having fewer sequences, we analyzed the effect of downsampling on the detection of clusters.

We first looked at the delay between the detection of the first cluster of each VOC and the first sequence of that VOC. A cluster is defined as a transmission event between three or more people (cf. Materials and Methods). When considering all available sequences, the delay in cluster detection was very VOC-dependent: 49 days for Alpha, 35 days for Delta, and 10 days for Omicron (Fig. 2B). Figure 2B shows the effect of downsampling on the delay to detect clusters. A commonly asked question during the pandemic was whether new cases were due to endemic transmission or new introductions. This is of epidemiological interest as it can affect the backward and forward tracing strategy. Delay upon downsampling was again more pronounced for Delta. When downsampling to 50k sequences, the delay in the detection of the first cluster was much shorter for Alpha (54 days) and Omicron (14 days) but more accentuated for Delta (63 days). The first Delta clusters appeared in May 2021, but it was not until July 2021 that the VOC gained importance, making it difficult to detect those early clusters with a modest sequencing strategy.

The distribution of the normalized density of clusters is shown in Fig. 4 for the whole data set. Upon downsampling to 50k sequences, we observed a similar pattern of peaks and valleys as with the complete data set, with some additional time windows now being without any clusters (e.g., beginning of July 2020 with two additional weeks without active clusters and in May 2022 with three additional weeks without active clusters) (Fig. 4).

Cluster characteristics such as the cluster duration and geographic spread help public health organizations understand the spread of an outbreak, e.g., a superspreading event, and thereby argue on the introduction of countermeasures. Thus, it is important to study how these cluster characteristics change when reducing the sequencing effort. For all the available genomes, a total of 3,014 clusters were found of which 54 were long-distance virus movements (LDVM) (1.8%), defined as events that spread a distance >200 km. The median cluster size included three cases (interquartile range (IQR) = 2), while for LDVM events, the median included five cases (IQR = 4.75). On the other hand, for a downsampling size of 50k sequences, 738 clusters were found, from which 17 were LDVM (2.3%). The median number of cases per cluster was three cases (IQR = 1), while it was four cases (IQR = 2) for LDVM events.

Downsampling had little effect on the duration of clusters, a result consistent for all VOCs. The median duration ranged from 10 days for the complete data set to 14 days for 2,500 sequences (Fig. S8).

Downsampling resulted in a reduction of the maximum distance within cases of a cluster (Fig. S9). Most clusters were found to be localized in a single canton. Clusters spreading up to five cantons were less frequent, and clusters appearing in six or more cantons were sporadic. A small yet significant correlation was found between the number of samples in a cluster and the number of affected cantons (Fig. S10). While the absolute number of captured LDVM decreased upon downsampling, their proportion increased, meaning that LDVM got enriched upon downsampling (Fig. 5).

All sequences (consensus sequences) shared with SPSP (https://www.spsp.ch) until the 5th of August 2022 were used (n = 143,260) in this study. Genomes sequenced with the explicit purpose of outbreak investigation were removed. The earliest genome collection date was 24 February 2020, the latest 1 August 2022. Genomes with incorrect collection dates (obvious data entry errors) were excluded. Quality control was done by the submitting laboratories (only genomes passing quality control of the individual labs are intended to be shared with SPSP). Pango lineage assignments were performed on SPSP in regular intervals and correspond in the data set to the latest Pangolin version available on 5 August 2022 (22). Additionally, Swiss COVID-19 case data were taken from the FOPH dashboard (https://www.covid19.admin.ch/, last accessed 20 September 2022).

The complete data set was downsampled to sizes ranging from 2,500 to 100,000 using random uniform sampling and 100 iterations. VOC labels (Alpha, Delta, and Omicron) were given to all genomes with assigned Pango sublineages associated with a given VOC. The introduction date was determined as the earliest collection date of a given VOC in the data set. The speed of introduction was calculated as the slope of a linear regression model of the log of exponential growth phase of a given VOC. The end of the exponential growth phase was determined based on the complete data set as 2021-03-01 (week 2021-09) for Alpha, 2021-08-01 (week 2021-30) for Delta, and 2022-01-01 (week 2022-01) for Omicron.

The lineage diversity was calculated as the Shannon diversity index (23) using the R package vegan (24). It accounts for both lineage richness and evenness, i.e., the number of lineages as well as their proportion. For example, the diversity in two samples with the same number of lineages would be greater in the sample with even proportions for each lineage compared to the sample with a few dominant and many rare lineages.

The nucleotide diversity was calculated for whole genome alignments [produced with mafft v7.505 (25) of all sequences per calendar week using the R package pegas (26).

A three-step method was used to identify sequence clusters. First, sequences with 0 single-nucleotide polymorphism (SNP) distance (including missing regions) that were collected maximally 20 days apart were grouped to obtain mutually exclusive clusters. These early clusters should have a minimum of three sequences. At this point, if a SNP difference occurred in a region that was missing in an otherwise identical sequence, this was counted as a SNP difference of 1, and the two sequences would be placed in different clusters. To overcome this and account for the possible mutations within a transmission chain, 0-SNP-clusters were expanded by adding sequences at 1 SNP distance. This resulted in multiple non-mutually exclusive clusters with overlapping sequences. Finally, the expanded clusters sharing sequences were merged using hierarchical clustering to create mutually exclusive superclusters. To ensure the quality of the superclusters, superclusters containing sequences collected more than 30 days apart were removed from the study. This time window was selected based on the estimated mutation rate of SARS-CoV-2 of 1.1  ×  10⁻³ subs/site/year (i.e., around two to three mutations per month) (27). To maximize the number of superclusters, the duration of 0-SNP-clusters was studied by comparing the number of sequences, cluster counts, and percentage of superclusters spanning more than 30 days for different duration thresholds (Fig. S11). A threshold of 20 days for the 0-SNP-clusters was found to maximize the number of sequences analyzed and the number of clusters found and to minimize the resulting number of superclusters exceeding 30 days.

The delay in VOC cluster detection was computed in relation to the first detection of a VOC sequence. The effect of downsampling on the detection of the first VOC clusters was analyzed by finding the distribution of the first cluster detection dates in the downsampled data sets. The downsample size of 50k samples was considered the smallest size in which the VOC detection dates could be retrieved without major delays (i.e., 1-day median delay for Delta and Omicron and no delay for Alpha; cf. Results on the first detection of VOCs upon downsampling and Fig. 2A). Therefore, this size was selected as the downsample size for further analysis.

Clusters containing the same sequences in the complete data set and a downsampled data set (50k sequences) were matched. Discrepancies in cluster matching can arise because (i) large clusters lasting longer than our 30-day threshold are discarded but are found in part in the downsampled data set or (ii) larger clusters from the complete data set were split into several clusters of shorter duration in the downsampled data set given that the common sequences are not in the sample anymore.

The distance between the cantons affected by clusters was computed as the Euclidean distance in kilometers between the capital of each canton.

Clusters were determined for the complete data set as well as the 50k downsampled data set (see above). A cluster was considered active during the period between the isolation of its first and last sequences. Given that the sequencing effort was not consistent during the study period, only considering the absolute counts of clusters per day would introduce a bias toward the periods with more sequencing effort. Therefore, the absolute number of clusters active each day was normalized by the number of sequences isolated in a time window of ±15 days (30 days being the estimated maximum duration of a cluster).