Resurgence of Omicron BA.2 in SARS-CoV-2 infection-naive Hong Kong

Genomic, epidemiologic, and human mobility datasets from Hong Kong

To elucidate the timing and origins of SARS-CoV-2 lineages during the fifth wave in Hong Kong, 116 saliva or nasopharyngeal samples from individual cases between 2 January and 4 February 2022, along with detailed epidemiological records including onset date, report date, and contact history were obtained from the Centre for Health Protection, Hong Kong. This study was conducted under ethical approval from the Institutional Review Board of the University of Hong Kong (UW 20–168). Because samples were collected as part of routine COVID-19 surveillance activities and were de-identified, a waiver of consent was granted. De-identified RT-PCR positive samples were sequenced using the same pipeline as in our recent studies^17,29. Full-genome analysis was conducted at a World Health Organization reference laboratory at the University of Hong Kong (Institutional Review Board no. UW 20–168). QIAamp Viral RNA Mini Kit (Qiagen, Cat. No.: 52906) was used to extract RNA. A number of gene-specific primers (https://github.com/Leo-Poon-Lab/mutations-under-sarscov2-vaccination/blob/main/Source%20Data/) targeting different regions of the viral genome were used to reverse transcribe the extracted RNA. For full-genome amplification, multiple overlapping 2-kb PCRs were performed with LA Taq DNA polymerase (Takara, Cat. No.: RR002M). The QIAquick PCR Purification Kit (Qiagen, Cat. No.: 28106) was used to purify PCR amplicons. DNA Prep (Illumina, Cat. No. 20018704) was used to prepare libraries from purified amplicons obtained from the same specimen. We quantified the libraries using Qubit dsDNA HS Assay Kits (Life Technologies, Cat. No.: Q32851) and sequenced them using Novaseq or iSeq100 sequencers (Illumina). All routine Hong Kong Delta and Omicron sequences deposited in GISAID until 30 April 2022 were also included. In addition, 10 random global (non-HK) sequences and 10 global sequences most similar by pairwise SNP distance to Hong Kong sequences per country per month from November 2021 to April 2022 were included (downloaded on 1 May 2022, Supplementary Data 4) as background to comprehensively and accurately define the monophyletic clade in Hong Kong and possible viral lineage exportations. Finally, reference genomes for each clade were included from GISAID (accessed on 8 May 2022, n = 258).

Pango lineage³⁰ was assigned to each sequence using Pangolin v.4.0.5, data version v1.3³¹. All nucleotide sequences were aligned to reference Wuhan-Hu-1 (GenBank accession MN908947.3), and those shorter than 27,000 nt were discarded. Duplicate sequences were removed, and sites deemed as problematic by other studies were masked (https://github.com/vjlab/omicronwave-hk) prior to phylogenetic analysis. Based on a regression of sample collection dates and root-to-tip genetic distances (from a maximum likelihood (ML) tree constructed in IQ-TREE 2³² and rooted with Wuhan-Hu-1: GenBank accession MN908947.3), sequences that did not deviate more than eight interquartile ranges were considered as high quality and retained for subsequent analysis. As a result, 3317 Hong Kong sequences and 5220 international sequences were included.

Epidemiological trends of confirmed cases, PCR results, and control measures in Hong Kong between January to April 2022 (Fig. 2a) were obtained from Centre for Health Protection (https://www.chp.gov.hk/en/index.html). Given that over 90% of the daily journeys in Hong Kong are made using public transport³³, changes in mobility during January–April 2022 grouped by children, students, adults, and the elderly were obtained from Octopus cards, which are ubiquitously used by the Hong Kong population for daily public transport and small retail payments (https://www.octopus.com.hk/tc/consumer/index.html).

Phylogenetic analysis

Bayesian time-scaled phylogenetic analyses were performed separately for Delta (HK = 126, global = 1426), Omicron BA.1.* (HK = 383, global = 2234), and Omicron BA.2.* (HK = 2807, global = 1361), as they evolved from ancestral SARS-CoV-2 strains independently (Fig. 1a). Molecular clock rates used as priors for the full datasets were estimated from a subset of genomes sampled as evenly as possible across epidemiological weeks (Delta, n = 150; Omicron BA.1.*, n = 181; and Omicron BA.2.*, n = 258) using the HKY + G4 + I substitution model with a strict molecular clock model and an exponential coalescent tree prior for the Omicron lineages and a constant coalescent for Delta. Six independent Markov Chain Monte Carlo (MCMC) chains were each run for 100 million steps, discarding the first 10 million as burn-in and resampling states every 2000 steps.

Lineages resulting from independent introductions into the Hong Kong community were inferred by estimating monophyletic clades from the full datasets using a Bayesian molecular clock phylogenetic analysis pipeline³⁴ implemented in BEAST (v.1.10)³⁵ (commit:d1a45). ML trees with branches scaled to genetic distance in IQ-TREE 2³² and time in TreeTime³⁶ were supplied as priors. Internal branches with less than one substitution were collapsed into polytomies. The analyses were run using a strict clock model with evolutionary rates estimated using the above subsampling datasets (Delta, 5.5 × 10⁻⁴; Omicron BA.1.*, 3.79 × 10⁻⁴; Omicron BA.2.*, 4.0 × 10⁻⁴ substitutions/site/year), the Skygrid population model with weekly grid points and a Laplace root-height prior with mean equal to the time-calibrated tree estimated by TreeTime³⁶ was used, with scale set to 20% of the mean. For each analysis, we ran 40 MCMC chains of 40 million, sampling every 60,000 steps with the first 4 million discarded as burn-in. Model convergence of mixing chains was inspected in Tracer (v.1.7.1)³⁷ to ensure an effective sample size (ESS) of >200 for each parameter. Monophyletic clades in the posterior trees were identified using the R package “NELSI”³⁸. It is notable that SARS-CoV-2 genomes with low variation among transmissions and our epidemiological data showed single introductions led to local outbreaks of HK-BA.2.2, HK-AY.127, and BA.1 (Dance cluster). Global sequences were therefore excluded when defining the three monophyletic clades. The R package “ggtree”³⁹ was used for tree visualization.

Phylogeography of HK-BA.2.2

To infer migration patterns of HK-BA.2.2 in the global context, we used a two-state (HK and non-HK) asymmetric discrete-trait analysis model implemented in BEAST v.10.1.4 with a HKY + G4 + I substitution model, an uncorrelated relaxed molecular clock model (the prior of 4.0 × 10⁻⁴ substitutions/site/year estimated for Omicron BA.2.*) with lognormal rate distribution (UCLN) and an exponential coalescent tree prior. For this analysis, we included the 10 earliest and 10 most recent sequences alongside 125 randomly selected cases from the HK-BA.2.2 monophyletic clade, 23 descendant sequences representing each country and province in mainland China, and two closely related ancestral BA.2.2 sequences (EPI_ISL_13330947 and EPI_ISL_9897214). We removed further outliers using TempEst v.1.5.3⁴⁰ under the premise that there is no major difference between the time signal of the dataset before and after sampling. As contact tracing and confirmatory phylogenetic analysis showed that HK-BA.2.2 virus was first detected in an international traveler arriving on 4 January 2022, an informative Laplace tMRCA of HK-BA.2.2 monophyletic clade prior with a mean (M) of 0.312 and a variance (s) of 0.01 was chosen. Six independent MCMC chains with 40 million states were performed, sampling every 2000 and discarding 10% as burn-in. As a result, 108,000 time-calibrated posterior trees were generated and used as an empirical distribution for the phylogeographic analysis. We combined two independent chains, each run for five million MCMC steps, sampling 1000 steps and discarding 10% as burn-in.

Effective population size (N
_e) and relative case detection rate

For the largest monophyletic clade (HK-BA.2.2, n = 2455) in Hong Kong, the above Bayesian molecular clock phylogenetic analysis pipeline with a strict clock fixed to 5.5 × 10⁻⁴ substitutions/site/year (mean value estimated from relaxed clock rate in phylogeography of HK-BA.2.2) was repeated to estimate changes in the effective population size (N_e) using Skygrid population model. Following Smith et al.⁴¹, by combining N_e and the epidemiological information of conducted tests, we can estimate the dynamics of the relative case detection rate:

subject to

where pop_infected is the number of infections in the population which can be simplified as a constant factor (c, which represents the number of true cases per effective population ‘unit’) times N_e due to their linear correlation. pop is population size (7.4 million) in Hong Kong, r_pos denotes the positivity rate of the PCR tests conducted, and n_t represents the number of tests conducted. Sensitivity sens and specificity spec were set to 1 as the reported COVID-19 cases until 26 February were confirmed twice by PCR tests. However, reducing sens does not change the dynamics of the relative case detection rate, but has an overall increase in the y axis in Fig. 3c.

Effective reproduction number (R
_e)

For improved computational efficiency and tested the effect of subsampling schemes (Supplementary Note and Supplementary Fig. 9) in constructing R_e, we used the sampling schemes recommended by the WHO for practical use in different settings and scenarios^42,43, which included uniform and proportional sampling, to construct three datasets (n = 262, uniform: 20 sequences per week; n = 502, uniform: 40 sequences per week; n = 897, proportional) summarized in Supplementary Fig. 11. A birth-death skyline serial (BDSS) model¹⁸ implemented in BEAST (v.2.6.7)⁴⁴ was used to infer the dynamics of the effective reproduction number (R_e). The HKY + G4 substitution model and a strict clock fixed to 5.5 × 10⁻⁴ substitutions/site/year (mean value estimated from relaxed clock rate in phylogeography of HK-BA.2.2) were used. Given that the BDSS model is affected by biases from sampling proportion (as shown in the sensitivity analysis in the Supplementary Note) and uneven sampling during the sequencing period from January to April, we assume that R_e and the sampling proportion are piecewise constant functions over 16 time intervals, roughly corresponding to weeks between 3 January and 26 April. Specifically, we assume that the sampling proportion per week is 0 before the collection time of the oldest sample, and is given a uniform distribution as prior with an upper bound on the empirical ratio of the number of subsampling sequences per week to the number of weekly reported cases. However, due to extensive sequencing done during the second week from 10 to 17 January (Fig. 1), where very few BA.2.2 cases were reported, the lower bound of sampling proportion prior was set at 0.3 (Supplementary Table 1, the upper and lower bounds of the sampling proportion prior could lead to a higher R_e between 10 and 17 January). A non-informative prior for tOrigin with lower bound set to 1 December 2021 was chosen. A lognormal prior with a mean of 0.0 and a variance (S) of 1.0 was set for R_e. To test the effect of the prior on R_e, we compared different levels of variance S (2 and 3) and found no significant differences to R_e shown in Fig. 3. Given that individuals who test positive in Hong Kong will be isolated, we assumed that there will be no further transmission from these individuals in our analysis. If this assumption is not valid, it could lead to an overestimation of the death rate and consequently the underestimation of R_e. The MCMC runs were performed for at least two independent chains of 100–200 million generations, sampling every 10,000 steps, with at least 10% discarded as burn-in. The R package “bdskytools” (https://github.com/laduplessis/bdskytools) was used to plot changes in R_e over time. The final R_e was selected from the estimation using the uniform subsampling dataset (40 sequences per week), which was better matched with the trend of R_t (Supplementary Fig. 10).

Instantaneous effective reproduction number (R
_t)

We computed R_t based on local cases and those epidemiologically linked to local cases, as defined by the Centre for Health Protection (CHP, https://www.coronavirus.gov.hk/eng/index.html). As SARS-CoV-2 can transmit pre-symptomatically⁴⁵, reconstructing incidence by date of infection provides a more accurate estimate of R_t⁴⁶. Therefore, we reconstructed the epidemic curve by infection date based on confirmation date with the distribution of delay from infection to confirmation using a deconvolution approach²⁸. We conducted the inference in a Bayesian framework and developed Markov chain Monte Carlo algorithms to estimate the posterior distribution of the model parameters and used a bootstrap approach to account for uncertainty associated with deconvolution⁴⁷. As Cori et al.⁴⁶ and Parag et al.⁴⁸ show, R_t measures the average transmissibility over a time window of length τ ending at time t under the assumption that R_t is constant within this time window, where τ is the smoothing parameter. In this study, we take τ = 14, to avoid unstable estimates for time-varying reproduction number. Correspondingly, the estimated R_t would need a few days to move to its true value, but still provide the correct direction of change²⁸.

Estimation of prevalence and incidence

Given the complex dynamics of the fifth wave in Hong Kong, we estimated point prevalence (I) from N_e τ, following a discrete generation model with arbitrary offspring distribution and changing population size⁴⁹. Due to the superspreading dynamics at SARS-CoV-2^21,50, a negative binomial offspring distribution was assumed, for which dispersion parameter (k) controls its shape. Point prevalence (I) can be calculated using the formula below:

subject to:

where N_e is the effective population size scaled by the generation length per year⁵¹ corresponding to N_e τ, where τ denotes generation time, R is the mean number of secondary cases, and k is the dispersion parameter of secondary cases. In this study, we used N_e estimated by a Skygrid coalescent model with 95% confidence interval (CI), τ = 2 and 3 days⁴, R replaced by the estimation of mean R_e using a BDSKY model and k from 0.05 to 0.2 (median = 0.1)²³. Cumulative incidence was calculated by adding the prevalence of each serial interval (2.72 days⁴) together, with the 95% CI restricted to the total population size (7.4 million). Daily incidence was calculated by diminishing cumulative incidence.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.