Epidemiology and analysis of SARS-CoV-2 Omicron subvariants BA.1 and 2 in Taiwan

Prevalence of Omicron subvariants between December 2021 and January 2023 in Taiwan

BA.1 and BA.2 and their sublineages entered Taiwan in December 2021 and January 2022, respectively (Supplementary Table S1 and Fig. 1). These two Omicron lineages did not cause COVID-19 outbreaks until March and April 2022 (Supplementary Table S1 and Fig. 2). Most COVID-19 cases were imported from abroad in the nonoutbreak period, and autochthonous COVID-19 cases began to outnumber imported COVID-19 cases in April 2022. BA.1 lineages dominated in January and February, and BA.2 lineages gradually replaced BA.1 lineages beginning in February (Fig. 1b). COVID-19 cases correlated with BA.1 lineages were reported primarily in northern Taiwan, and BA.2 lineages spread widely in southern Taiwan. BA.2 lineages dominated between March and June, as they did worldwide. Further results revealed that BA.2.3.7 dominated between April and July 2022. BA.2.3.7 accounted for 78.7%, 86.1%, 86.8%, and 89.2% of all isolated SARS-CoV-2 infections between April and July. BA.5 and BA.4 lineages entered Taiwan in June and July 2022, respectively. BA.5 lineages replaced the BA.2 lineages rapidly and composed 54.5% to 100% of sequenced SAR-CoV-2 between August and October 2022. BA.5 lineages were gradually replaced with other Omicron lineages (e.g., BF.7, BN.1, BQ.1, and some other non-BA.5 lineages) from November 2022 to January 2023 (Fig. 1b). The geographical distribution of autochthonous cases between December 2021 and January 2023 is shown in Supplementary Fig. S1. Most of the cases were distributed in New Taipei City, Taichung City, and Kaohsiung City, which are the top three cities in terms of total population, and Taoyuan City, where the largest international airport is located. Although the sample sizes of sequences deposited in GISAID were small between August 2022 and January 2023, the trend of sequence replacement was similar to the data released by the Central Epidemic Command Center (CECC) in Taiwan. The Taiwan CECC started to release sequencing results weekly from mid-October 2022. BA.5 lineages overwhelmingly dominated between October and December 2022 (82–99%). BA.5 lineages made up 62–55% of cases during most of January 2023, declined to 35% in the last week of January and were replaced by the BA.2.75 lineage¹⁴.

Source of data: GISAID https://gisaid.org/.

Monthly distribution of different variants of concern and clade replacement over time between December 2021 and January 2023 in Taiwan. (a) Distribution of the Delta variant, Omicron variant, and their subvariants. (b) Clade replacement during the emergence of the Omicron variant and its subvariants in Taiwan.

Source of data: The web-based notifiable disease surveillance system maintained by the Taiwan Centers for Disease Control (CDC) https://nidss.cdc.gov.tw/nndss/disease?id=19CoV.

Monthly COVID-19 data between December 2021 and January 2023 in Taiwan.

Detection and analysis of SARS-CoV-2 identified in this study

We collected nasopharyngeal swabs from the patients suspected with COVID-19 in KMUH from January 2022 to September 2022. Among the collected nasopharyngeal swabs, SARS-CoV-2 genomic RNA was detected in ninety-one swab-VTM samples by using real-time RT‒PCR. We evaluated the integrity of the SARS-CoV-2 genome in these swab-VTM by using virus culture. We observed a CPE in 42 out of 91 swab-VTM samples inoculated into Vero E6 cells either in cells inoculated with original swab-VTM or in cells with blind passage. We noticed that the CPE was observed only in the swab-VTM samples with a cycle threshold (Ct) value between 16.3 and 23. The presence of SARS-CoV-2 in the culture medium was confirmed by using real-time RT-PCR. The results suggested all CPEs observed were caused by the presence of SAR-CoV-2 in the culture medium. The swab-VTM samples that caused Vero E6 cells to produce a CPE were sent for next-generation sequencing for either the full-length genome or spike gene. The lineage and clade assignments of the sequencing results were analyzed by using the Pangolin COVID-19 Lineage Assigner¹⁶ and Nextclade¹⁷. The results revealed that 2 sequences were assigned to BA.1 lineages, 36 sequences were assigned to BA.2 lineages, and 4 sequences were assigned to BA.5 lineages. These results were basically consistent with those internationally at that time. We focused on the BA.1 lineages and BA.2 lineages in the later analysis. These sequences were uploaded to GenBank and GISAID; the demographic characteristics are listed in Supplementary Table S2. No patients died from Omicron variant infection in our study.

To understand the phylogenetic relation of the SARS-CoV-2 variants identified in this study with those in public SARS-CoV-2 databases, such as GISAID, GenBank, COVID-19 Genomics UK Consortium, and the China National Center for Bioinformation, we analyzed the 13 full-length sequences identified in this study with Ultrafast Sample placement on Existing tRee (UShER, https://genome.ucsc.edu/cgi-bin/hgPhyloPlace) first. The analysis results suggested that the 13 full-length genomes identified in this study, including those for the BA.1, BA.1.1, BA.2, BA.2.3.7, and BA.2.64 lineages, were clustered in 11 subtrees in the existing tree with 14,359,598 SAR-CoV-2 genomes. The other 25 spike gene sequences were merged with selected sequences from 11 subtrees to reconstruct a single phylogenetic tree by using IQ-TREE 2.2.0 with their spike genes as described in our previous studies^12,15. There were 257 spike sequences, including 251 Omicron BA.1/BA.2 sequences and 6 reference sequences (1 wild type, 1 Alpha variant, 2 Delta variants, and 2 Omicron BA.5.2), in the reconstructed phylogenetic tree. The results revealed that approximately 60% of BA.1, BA.2, and their sublineages in this study were located among sequences identified in Taiwan, and some others were closer to the sequences identified in other countries (Fig. 3). In general, autochthonous COVID-19 cases in Taiwan were possibly initiated by the importation of SARS-CoV-2 by travelers returning to Taiwan or international travelers visiting Taiwan.

Phylogenetic analysis of Omicron subvariants identified in this study. UShER enables real-time sequence placement for the SARS-CoV-2 pandemic using an existing phylogenetic tree generated by the sarscov2phylo pipeline, containing 14,359,598 genomes from GISAID, GenBank, COG-UK, and CNCB (March 31, 2023). Some nearest neighbor SARS-CoV-2 sequences already in the existing phylogenetic tree were selected for further analysis. A total of 251 Omicron BA.1/BA.2 spike gene sequences, including those selected from the UShER results and spike gene sequences identified in this study, were analyzed by using the maximum likelihood and fits of 484 different nucleotide substitution models. The results suggested that K3Pu + F + R2 was the best-fitting model with the lowest Bayesian information criterion score of 11,355.955 among the models tested. Tree topology was automatically computed to estimate maximum likelihood values. The optimal log-likelihood for this computation was − 5089.245. There was a total of 3,180 positions in the final dataset. An original tree is displayed using FigTree v1.4.4 with a color indicator for bootstrap values and a scale bar. Viruses are shown as the Virus name| EPI_ISL ID in GISAID. The BA.1 lineage includes BA.1, BA.1.1, and BA.1.14. The BA.2 lineage includes BA.2, BA.2.3.7, BA.2.12.1, BA.2.9, and BA.2.64. The Omicron subvariants identified in this study are highlighted in purple (38 strains). Asterisk, BA.5.2; filled diamond, B.1.617.2.43 (Delta); filled square, B.1.1.7 (Alpha); filled circle, Wild-type. Wuhan/WIV04/2019 was used as the root of this tree.

We analyzed genomic variations, including single nucleotide variation (SNV) and insertion/deletion (Indel), of the 38 SARS-CoV-2 BA1 and BA.2 lineages identified in this study by comparing their sequences to the reference sequence Wuhan-Hu-1/2019 (MN908947), and we focused on the spike gene. The results revealed that the BA.1 and BA.2 lineages acquired additional genetic variations in the spike gene as they spread and evolved when compared with the original BA.1 and BA.2 lineages. These genetic variations resulted in nonsynonymous codon variations, such as L5F (5.3%), A163V (2.6%), S255F (2.6%), R346K (5.3%), A879S (2.6%), G1251V (13.2%), and K97E (78.9%) (Fig. 4). We found that spike R346K appeared only in BA.1.1, and spike K97E and G1251V were characteristic in BA.2.3.7. In addition, some reverse mutations were observed.

Genetic variation frequencies in BA.1, BA.2, and their sublineages isolated in this study.

Spike gene variations in BA.1 lineages and BA.2 lineages identified in Taiwan

We downloaded full-length genomes of BA.1 lineages and BA.2 lineages identified between December 2021 and January 2023 in Taiwan from GISAID. The quality of SARS-CoV-2 genomes was analyzed by using Nextclade, and only high-quality sequences were preserved for further analysis. The spike gene in 82 out of 93 BA.1 lineages, including BA.1.1, B.1.15, BA.1.16, BA.1.17, and BA.1.20 sublineages, and in 2010 out of 2069 BA.2 lineages, including BA.2.10, BA.2.3, BA.2.3.7, BA.2.68, and BA.2.75.2 sublineages, were aligned and compared with the spike gene of the reference sequence Wuhan-Hu-1/2019 (MN908947). The genetic Indels and SNVs resulting in amino acid (AA) deletions and nonsynonymous codon variations are shown in Table 1. The results revealed that the BA.1 sublineages and BA.2 sublineages acquired additional gene deletions and nonsynonymous codon variations in the spike gene compared with the original BA.1 and BA.2 lineages. Additional AA deletion in the spike protein was found at a low frequency in the BA.1 and BA.2 lineages (0.05% to 2.4%). We searched the genetic deletions found in this study in GISAID. We found 215 Omicron sequences with spike 61-73del (82.8% were in BA.1 sublineages) and 336 Omicron sequences with spike 81-95del (24.5% were in BA.1 sublineages). There was no spike 81-95del/166-180del and only 33 Omicron sequences with spike 361-376del/380-397del (3.3% were in BA.1 lineages, 36.4% were in BA.2 sublineages, and 15.2% were in BA.5.2 lineages) in the database. There were 45 Omicron sequences (including BA.1, BA.2, and BA.5 lineages) with similar spike AA deletions found in BA.2 lineages in this study. These genetic events were rare since the numbers of Omicron sequences deposited in GISAID were over 7.9 million at the time of analysis. In the context of nonsynonymous codon variations, the results revealed that additional spike R346K (40.2%) was found exclusively in BA.1 sublineages, and additional spike K97E (84.8%) and G1251V (59.5%) were found exclusively in BA.2 sublineages. Further analysis revealed that spike R346K was exclusively found in the BA.1.1 sublineage (100%). In addition, spike K97E (100%) and G1251V (71%) were characteristic mutations observed in the BA.2.3.7 sublineage (Table 2).

The CFR and severity of COVID-19 with Omicron subvariants

The accumulated CFR gradually decreased after BA.1 and BA.2 entered Taiwan in late 2021 and early 2022 with the accumulation of confirmed COVID-19 cases with a limited number of deaths between December 2021 and February 2022. During the outbreak related to the BA.2 sublineages, the accumulated CFR decreased drastically (Supplementary Table S1, Fig. 1, and Supplementary Fig. S2). The accumulated CFR dropped drastically after Omicron subvariants, mainly BA.2.3.7, entered Taiwan after April 2022. The CFR was 0.17% in Taiwan, which was lower than that worldwide (0.25%) between April 2022 and January 2023, and the Omicron subvariants, regardless of sublineage, dominated.

Both the accumulated and monthly CFR in Taiwan dropped rapidly after the outbreak in April 2022 to a level lower than the global CFR (Table 3). The CFRs for BA.2.x were 84.6% (monthly CFR 0.35%) and 91.7% (monthly CFR 0.31%) in April 2022 and May 2022, respectively. The dominant subvariants in Taiwan were BA.2.3.7 (86.1%, monthly CFR 0.21%) and BA.2.x (86.8%, monthly CFR 0.22%) in May 2022 and June 2022, respectively. The CFR for BA.2.x was 1.4–1.7-fold greater worldwide than that of BA.2.3.7 in Taiwan.

The COVID-19 vaccination program started in late March 2021. National vaccination rates rose rapidly in response to the COVID-19 outbreak related to the Alpha variant beginning in May 2021. The approved COVID-19 vaccines in Taiwan are AstraZeneca ChAdOx1-S (March 2021), Moderna mRNA-1273 (June 2021), BioNTech BNT162b2 (August 2021), Spikevax Bivalent Original/Omicron BA.1 (September 2022), and Spikevax Bivalent Original and Omicron BA.4/BA.5 (September 2022). Basically, the CFR was inversely correlated with the overall vaccination rate from then on (Supplementary Fig. S1).

According to the data released by the Taiwan CECC, moderate and severe COVID-19 cases increased rapidly from April 2022 after the outbreak related to BA.2 sublineages. COVID-19 patients over 70 years of age were at higher risk for moderate to severe disease and worsened after the Omicron-related outbreak, and cases increased from 60.8% in April 2022 to 71.6% in July 2022 (Supplementary Table S3). Among moderate and severe COVID-19 patients, 91% of patients had at least one chronic disease, and 91% of patients were > 60 years old. The rates for moderate and severe COVID-19 were 0.20–0.27% and 0.03–0.21%, respectively, between January 2022 and January 2023. We observed an upward trend in severe COVID-19 rates between April 2022 and July 2022 (Fig. 5). We also noted that the nonsynonymous codon variations in the spike protein resulted from genetic mutation increased during the prevalence of Omicron subvariants in Taiwan. The results suggested that AA variations might play a role in disease severity. Taken these results together, the severity of Omicron-related COVID-19 was associated with older age, the presence of chronic diseases, and the accumulation of AA variations in the spike protein.

Monthly data on genetic variation frequencies in the spike protein, illness, and COVID-19 cases between December 2021 and January 2023 in Taiwan. The percentage of illness is the original data multiplied by 100 for presentation. Data source: GISAID https://gisaid.org/ and Taiwan CDC https://nidss.cdc.gov.tw/nndss/disease?id=19CoV.