Genomic surveillance of SARS-CoV-2 in mainland China after ending the zero-COVID policy, December 2022–January 2023

After implementing the dynamic zero-COVID policy for nearly three years, China has recently begun to adjust its COVID-19 prevention and control strategies with the attenuated pathogenicity of omicron subvariants and increasing vaccination coverage, notably by announcing the 10 measures on 7 December 2022. A recent publication in this journal by Jingyi Liang et al. showed that the cumulative SARS-CoV-2 infection rate rose rapidly to 70% within three weeks after the ending of zero-COVID policy in Macao, China. Herein, we report the genomic characteristics and phylodynamic of SARS-CoV-2 complete genome sequences in mainland China after ending zero-COVID policy based on the viral genomic surveillance data deposited in the Global Initiative on Sharing All Influenza Data (GISAID) database.

A total of 8254 complete SARS-CoV-2 genomic sequences from mainland China were extracted from the GISAID database during the period of December 2022 and January 2023 (date of access, February 6, 2023). Sequences were typed using the Phylogenetic Assignment of Named Global Outbreak Lineages (PANGOLIN; version 4.2) web tool. These genomes corresponded to 88 Pango-nomenclature-system-named subvariants, with the dominant lineages of BA.5.2.48 (4881/8254, 59.1%) and BF.7.14 (2223/8254, 26.9%), followed by BA.5.2.49 (525/8254, 6.4%), BA.5.2 (205/8254, 2.5%) and BF.7 (93/8254, 1.1%) (Fig. 1A). The proportion of these dominant lineages were not significantly changed over the 2 months of outbreak (Fig. 1B). While the other 83 lineages accounted for only 4.0%, such as BQ.1.1, BN.1.3 and BN.1.2 (Fig. 1A). What is more, the proportion of these lineages had fallen consistently over the time (Fig. 1B). All these findings demonstrated that the emerging BA.5.2.48 and BF.7.14 were the absolutely dominant drivers of the current COVID-19 outbreak, and the lineages found to be dominant internationally during the same period (i.e., XBB and BQ.1) were quickly cleared and did not prevail in China.

In order to investigate the phylodynamics of BA.5.2.48 and BF.7.14 subvariants, we randomly chose complete genomic sequences of the subvariants from mainland China collected on or before January 31, 2023 in the GISAID database, respectively (date of access, February 16, 2023). To test for temporal structure in the datasets, we performed linear regressions of root-to-tip genetic distances against sampling dates using TempEst (version 1.5.3). The results indicated that both the BA.5.2.48 (n = 350) and BF.7.14 (n = 357) datasets were considered to have sufficient temporal signals for molecular dating after discarding several outliers (Supplemental Fig. S1). Then, Bayesian time-scaled phylogenetic analyses were performed using Markov Chain Monte Carlo (MCMC) methods implemented in BEAST software (version 1.10.4) to estimate the time of the most recent common ancestor (tMRCA), rate of evolution, as well as to infer the effective population size for both BA.5.2.48 and BF.7.14 subvariants. A GTR +Γ substitution model was used along with a strict clock model. A Bayesian skyline plot was used as the tree prior to estimate the median effective population size through time with a 95% highest posterior density (95% HPD). The MCMC chain length was set at 350 million and 300 million states for the BA.5.2.48 and BF.7.14, respectively. Convergence was assessed using Tracer (version 1.7.2), and 10% of states were removed to account for burn-in. Finally, maximum clade credibility (MCC) trees were generated using TreeAnnotator (version 1.10.4).

Our Bayesian phylodynamic analysis revealed that the tMRCA for the BA.5.2.48 and BF.7.14 subvariants were estimated to be 16 July 2022 (95% HPD, 14 June–11 August 2022) and 9 August 2022 (95% HPD, 28 July–15 August 2022) with posterior probability (PP) = 1, respectively (Fig. 2AB). The estimated emerging dates of the subvariants were slight earlier than the first complete genomic sequences collected in GISAID database (15 August 2022 for BA.5.2.48 and BF.7.14). The mean substitution rates of BA.5.2.48 and BF.7.14 were estimated to be 5.353 × 10⁻⁴ (95% HPD, 4.141 × 10⁻⁴ to 6.893 × 10⁻⁴) substitution/site/year and 5.332 × 10⁻⁴ (95% HPD, 4.562 × 10⁻⁴ to 6.129 × 10⁻⁴) substitution/site/year, respectively. The results were in line with a previous study that estimated the average substitution rate of global SARS-CoV-2 (all lineages) was 5.83 × 10⁻⁴ substitution/site/year. Indeed, at the beginning of the pandemic, the evolutionary rate of SARS-CoV-2 was about 1.1 × 10⁻³ substitution /site/year, which means that the new variants present about 2 times slower than the original variants.^,  The Bayesian skyline plot (BSP) showed that the effective population size for both BA.5.2.48 and BF.7.14 underwent expansions during the early phase of the outbreak in November 2022, matching the change of prevention and control policies by announcing the 20 measures during the same period in China. For BA.5.2.48, the BSP analysis revealed one timeframe with a rapid exponential growth from around 7 November to 25 November 2022, continuing at a high average until the end of the sampling time (Fig. 2C). However, the effective population size for BF.7.14 increased earlier than BA.5.2.48 and showed two growth intervals over the time (Fig. 2D). The difference in the growth patterns of the effective population size of these two subvariants might be due to the different viral fitness and the adjustment of COVID-19 control strategies. Besides, which was unexpected, the difference possibly arose from bias due to the limited sample size.

This study provides insights into the genomic epidemiology of SARS-CoV-2 in mainland China after ending the zero-COVID policy. We report that BA.5.2.48 and BF.7.14 were the dominant drivers of the current outbreak. Continuous genomic surveillance of SARS-CoV-2 is crucial for monitoring emerging variants with the potential to cause new COVID-19 outbreaks in China.