The Molecular and Biological Patterns Underlying Sustained SARS-CoV-2 Circulation in the Human Population

Full Text

Introduction

The COVID-19 pandemic was caused by a novel SARS-CoV-2 virus discovered in late 2019 and has caused more than 774 million cases and more than 7 million deaths worldwide [1]. The observed evolution of SARS-CoV-2 has led to the independent emergence of several genetic lineages identified by the World Health Organization at certain periods as variants of concern (VOC) – Alpha, Beta, Gamma, Epsilon, Delta and Omicron [2], which have been associated with upsurges in morbidity and mortality in most countries. Despite the rapid development and introduction of effective prophylactic drugs following the first year of the pandemic [3–6], as the first VOCs emerged, reports of decreased vaccine efficacy began to emerge [7–9]. This became the core reason for the introduction of booster doses [10, 11]. In the case of the Omicron variant, efficacy decreased to such an extent that booster doses with a modified formulation had to be administered [12–15]. However, in the case of the relatively recent XBB variant, even with the use of an updated divalent booster followed by a monovalent vaccine containing the XBB antigen, a decrease in neutralization has been observed [8]. Furthermore, a large amount of information has now been accumulated on the impact of mutations characteristic of dominant variants on the efficacy of therapeutic agents including monoclonal antibodies and antiviral agents [16–18], which has resulted in the loss of efficacy of most monoclonal antibodies (mAbs) against topically circulating variants [19–21]. The emergence of a new virus variant was preceded by the emergence of mutations, which in turn affected changes in the transmissibility and/or pathogenicity of the virus, sensitivity to therapeutic drugs and ability to evade natural or vaccine-induced immunity. It is worth noting that indicators such as viral transmissibility and/or pathogenicity are inextricably linked to the magnitude of viral load. As discussed previously [22], host (role of vaccination or prior infection) and viral factors (SARS-CoV-2 variants) significantly influence viral load dynamics and therefore further influence viral transmission. Despite the accumulated data on the values of some mutations and associated phenotypic traits [22], the factors influencing the dynamics of the epidemiologic process in different periods remain incompletely understood.

This study is based on the results of molecular genetic monitoring of SARS-CoV-2 virus variants, some of which have been previously published [23]. In this study, we assessed both the dynamics of circulating variants and the potential impact of mutations and viral load as the main factors contributing to changes in the profile of circulating lineages, which may be reflected in the observed transition: «increase in mutations ‒ increase in viral load ‒ emergence of a new serotype».

Materials and methods

During the period of Delta variant emergence, nasopharyngeal swabs were collected regularly from different groups of volunteers, including those who were initially infected or vaccinated. Written consent was obtained from all patients in accordance with the order of the Russian Ministry of Health of July 21, 2015, No. 474n «On the procedure for giving informed voluntary consent to medical care within the framework of clinical approbation of methods of prevention, diagnosis, treatment, and rehabilitation, the forms of informed voluntary consent to medical care within the framework of clinical approbation of methods of prevention, diagnosis, treatment, and rehabilitation, and refusal of medical care within the framework of clinical approbation of methods of prevention, diagnosis, treatment, and rehabilitation». All samples were depersonalized before they were received by the research team. The study was approved by the local ethical committee of the N.F. Gamaleya NRCEM (protocol No. 14, September 29, 2021).

Sample collection and viral load testing were performed in two laboratories (hereinafter referred to as Laboratory 1 and Laboratory 2). During the circulation period of Wuhan, Delta and Omicron variants (BA.1.X and BA.2.X) (Laboratory 1), nasopharyngeal swabs were collected in virus transport medium (G00155, GEM, Russia). Total RNA was isolated using the QIAamp Viral RNA Mini Kit (Qiagen, Germany) and a kit for isolation of total RNA from animal and bacterial cells, smears and viruses on columns (RU-250, Biolabmix, Novosibirsk, Russia). Quantitative reverse transcription polymerase chain reaction (PCR) was performed using the SARS-CoV-2 FRT RT-PCR kit (EA-128, N.F. Gamaleya NRCEM, Moscow, Russia). Samples with Ct values < 30 were selected for full-genome sequencing.

In Laboratory 2, nasopharyngeal swabs were collected in a virus transport medium (physiological solution or XK-PCR30 transport medium (Jiangsu Xinkang Medical Instrumet Co., Ltd., China)) or GEM transport medium (GEM, Moscow, Russia). During the circulation period of Omicron BA.1.X, BA.2.X, BA.5.X variants, quantitative reverse transcription PCR was performed using the AmpliPrime SARS-CoV-2 DUO kit (NextBio, Moscow, Russia) according to the manufacturer’s instructions. During the circulation period of Omicron BA.5.X, CL.X, XBB.1.X and XBB.1.9.X variants, quantitative reverse transcription PCR was performed using CoV-2 test (TestGen, Ulyanovsk, Russia) according to the manufacturer’s instructions.

Further whole genome sequencing was performed using Ion Torrent (Thermo Fisher Scientific, USA), Illumina (Illumina, USA), and Oxford Nanopore (Oxford, UK) technologies. The sample preparation and analysis procedures have been described in detail previously [21, 23]. The whole genome sequences obtained in this study were uploaded to the GISAID database under the following numbers: EPI_ISL_1710849-1710866, EPI_ISL_2296111-2296286, EPI_ISL_2296288-2296379, EPI_ISL_4572812, EPI_ISL_5334362-5334371, EPI_ISL_5334374-5334389, EPI_ISL_7211325-7211326, EPI_ISL_7263932-7263933, EPI_ISL_9230058-9230062, EPI_ISL_9230064-9230100, EPI_ISL_10627062, EPI_ISL_11864996-11865125, EPI_ISL_11872910, EPI_ISL_421275, EPI_ISL_454732, EPI_ISL_470896-470904, EPI_ISL_572398, EPI_ISL_872628-872643, EPI_ISL_875515, EPI_ISL_1015362, EPI_ISL_1708507-1708509, EPI_ISL_12225322, EPI_ISL_12748381-12748382, EPI_ISL_13431664-13431687, EPI_ISL_14217225-14217226, EPI_ISL_15327072-15327075, EPI_ISL_15858138-15859137, EPI_ISL_15860713-15860737, EPI_ISL_15860739-15860839, EPI_ISL_15860841-15860991, EPI_ISL_15860993-15861048, EPI_ISL_15862338-15863336, EPI_ISL_15863677-15864655, EPI_ISL_15864802-15865776, EPI_ISL_15865821-15866801, EPI_ISL_15867150-15868141, EPI_ISL_15868158-15869145, EPI_ISL_15869218-15870209, EPI_ISL_15871156-15872140, EPI_ISL_15872157-15873150, EPI_ISL_15873159-15874146, EPI_ISL_15874159-15875146, EPI_ISL_15875638-15876623, EPI_ISL_15876640-15877626, EPI_ISL_15879747-15880730, EPI_ISL_15883551-15884536, EPI_ISL_15884833-15885823, EPI_ISL_15885995-15886980. The sequences obtained were also uploaded to the VGARus database.

The information on the circulation of genetic lineages used in this study was supplemented with data obtained from the GISAID database. The obtained data were filtered by the following parameters: Host == “Human”, `Sequence length` >= 27000, `Is complete?` == “TRUE” and Location == “Europe/Russia/Moscow/Moscow” or “Europe/Russia/Moscow” (query date September 18, 2023). Additionally, the retrieved sequences were also filtered by collection date (`Collection date` < “2023-06-01”), resulting in 16,541 records.

The genetic variant information for the remaining sequences was grouped as follows: Wuhan (B.1 + B.1.X), Alpha (B.1.1.7 + Q.X), Beta (B.1.351), Delta (B.1.617.X + AY.X), BA.1.X (B.1.1.529 + BA.1.X), BA.2.X, BA.5.X (BA.5.X + BF.X + BE.X), CL.X, XBB.1.X (without XBB.1.9.X and XBB.1.16.X), XBB.1.9.X, XBB.1.16.X (XBB.1.16.X + FU.X), and «Other XBB»/«Other XBB» (Other XBB lineages + EG.X). Furthermore, the genetic lineages B.1.1.317 and B.1.1.523 were added. The remaining lineages were combined into the «Other»/«Other» group. The obtained data were merged with metadata and visualized in R environment using the dplyr v. 1.1.1.1 [24] and ggplot2 v. 3.4.2 [25] packages.

The Kraskell–Wallis test was used to compare more than two independent groups. If the results were statistically significant, Dunn’s test was applied to determine the difference between groups. A comparison of PCR threshold cycle values was performed among samples where the same quantitative reverse transcription PCR kit was used. For statistical analysis, the ggstatsplot v. 0.12.1 [26] and R v. 4.2.2 packages were used.

Results

Dynamics of SARS-CoV-2 genetic variants in Moscow since the beginning of the COVID-19 pandemic

Since the initial spread of SARS-CoV-2 in Moscow up until February 2021, the predominant strains were characterized as Wuhan (B.1.X), including its diverse genetic sublineages, among which the circulation of a primarily Russia-specific genetic lineage, B.1.1.317, is notable. In winter 2021, the previously circulating variants began to be gradually replaced by the Alpha (B.1.1.7 + Q*) and B.1.1.523 lineages (also characteristic mainly of Russia), which were displaced by the Delta variant (B.1.617.2 + AY*) in May 2021 together with the next COVID-19 incidence increase. By the end of 2021, the Omicron variant (B.1.1.529) succeeded Delta, followed by sequential transition of the genetic lineages BA.1.X, BA.2.X, and BA.5.X. In October 2022, CL.X becomes one of the dominant variants, such a widespread distribution of which has not been reported elsewhere. Along with CL.X, a recombinant variant XBB began circulating in Moscow in November 2022, subsequently displacing the other SARS-CoV-2 genetic lineages (Fig. 1 a).

 

Fig. 1. Genetic variants of SARS-CoV-2 circulating in Moscow since the start of the COVID-19 pandemic. a – dynamics of genetic variants of SARS-CoV-2 circulating in Moscow. The left Y-axis indicates the proportion of genetic lineages, while the right axis represents the number of new COVID-19 cases per 100,000 population; b – changes in the number of mutations in the predominant circulating variants of SARS-CoV-2. The ordinate axis displays the mutation count. The periods (months) are indicated on the abscissa axis.

Рис. 1. Генетические варианты SARS-CoV-2, циркулировавшие в Москве с начала пандемии COVID-19. a – динамика генетических вариантов SARS-CoV-2. Слева по оси ординат показана доля генетических линий, справа – число новых случаев COVID-19 на 100 тыс. населения; б – изменение количества мутаций в основных циркулировавших вариантах SARS-CoV-2. По оси ординат показано количество мутаций. На оси абсцисс отражены периоды (месяцы).

 

Different variants of the virus, which became relatively widespread, had a characteristic profile of mutations.

Consideration of their total number in the S-protein and other genes showed that before the appearance of the Omicron variant there was a gradual increase in their number; the spread of the latter was characterized by a more than twofold increase in the number of nonsynonymous mutations in the SARS-CoV-2 genome (Fig. 1 b). Subsequently, substitutions also continued to appear in the viral genome, but the first Omicron variant, BA.1.X, characterized by 34 mutations in the S-protein, which is 3 more than in the subsequent BA.2.X, and comparable to BA.5.X, falls out of the general series. Despite this fact, the total number of mutations in the genome of the SARS-CoV-2 variant BA.1.X was the lowest among all other Omicron lines and amounted to 60, whereas the last variant considered, XBB.1.9.X, had about 80 amino acid substitutions, of which 41 were in the S-protein.

Mutation profile in the major genetic variants of SARS-CoV-2.

A detailed review of mutations in the selected variants revealed a characteristic mutation profile for the SARS-CoV-2 variants circulating in Moscow. In the RBD of the S-protein, the N501Y mutation was observed for Alpha, E484K and S494P for B.1.1.523. Delta was characterized by two amino acid substitutions (L452R and T478K) in the RBD of the S-protein, and 13 for BA.1.X, and the occurrence of the K417N, N440K, and G446S mutations in some genetic lineages related to this Omicron variant was also observed. The BA.2.X variant is notable for 16 mutations in the RBD, of which only 10 are common with BA.1.X; the dominating variants do not exhibit the R346K and S371L mutations present in BA1.X, and, furthermore, a series of additional substitutions appeared, including S371F, T376A, D405N, R408S, K417N, and N440K. The variants BA.5.X and CL.X had the L452R mutation characteristic of the Delta lineage and F486V, and all further circulating variants were also characterized by amino acid substitutions at position 486. From the end of 2022, the dominant genetic lineage in Moscow was XBB and its subvariants (Fig. 1) with a characteristic minimum of 22 amino acid mutations in RBD, of which 6 were not found in previously circulating variants (Fig. 2). In spring 2023, the main circulating variant was XBB.1.9.X, namely XBB.1.9.1, characterized by the F486P substitution (Fig. 2).

 

Fig. 2. Mutations in the RBD of the spike protein of the predominant SARS-CoV-2 variants. The Y-axis represents SARS-CoV-2 variants, and the X-axis represents amino acid substitutions. Zero values reflect the proportion of mutations tending to zero; empty cells represent the absence of mutations at a specific position.

Рис. 2. Мутации в RBD Spike-белка основных доминировавших вариантов SARS-CoV-2. По оси ординат отражены варианты SARS-CoV-2, по оси абсцисс – аминокислотные замены. Нулевые значения отражают долю мутаций, стремящуюся к нулю; пустые ячейки – отсутствие мутаций в данной позиции.

 

A wide variety of amino acid substitutions was observed in the N-terminal domain (NTD) of the S-protein. The genetic lineage B.1.1.523 was characterized by the presence of the F306L mutation. Genetic lineages belonging to the Delta variant had mutations that were not subsequently found in other variants: T19R, E156G, F157del, R158del, and one G142D mutation characteristic of all Omicron lineage variants except BA.1.X. Overall, most BA.1.X mutations were uncharacteristic of subsequent Omicron variants with the exception of deletions 69‒70, which were also present in BA.5.X and CL.X. CL.X, compared to BA.5.X and other Omicron variants, possessed a K150E amino acid substitution; otherwise, variants BA.2.X, BA.5.X, and CL.X had a similar mutation profile. The XBB variants were characterized by V83A, H146Q, Q183E, and V213E mutations not previously seen in other genetic lineages. All variants had a D614G substitution in the C-terminal domain (CTD), which is known to be fixed before the SARS-CoV-2 genetic lineages were divided into VOCs [27]. For Alpha, the largest number of amino acid mutations was observed in the CTD region: P681H, T716I, S982A, and D1118H, not including D614G. B.1.1.523 was characterized by the presence of D839V and T1027I, not further observed in other variants. Delta in CTD had only 2 substitutions besides D614G: P681R near the furin site and D950N. Further, BA.1.X had the highest number of mutations in this region compared to the other Omicron variants: in addition to the characteristic H655Y, N679K, N681H, N764K, D796Y, Q954H, N969K, the first Omicron variant had additional T547K, N856K, and L981F. Overall, little further evolution of the Omicron NTD was observed, but T883I substitution was observed in CL.X and XBB.1.X in a small percentage of cases (29 and 12, respectively) (Fig. S3 a).

The frequency of mutations in non-structural genes was much lower than in the S-protein. The Delta variant was characterized by about 16 mutations in nonstructural genes, and mutations T492I in NSP4 and P323L in RdRp (NSP12) were also characteristic of all subsequent Omicron variants, and G671S in the same gene was characteristic of XBB genetic lineages. The remaining mutations, however, were not further encountered. The first Omicron variant, BA.1.X, had amino acid substitutions NSP3_K38R, NSP3_L1266I, NSP3_A1892T, NSP6_I189V and deletions NSP3_S1265del and NSP6_L105del, which were not found in Omicron variants, as well as a number of other mutations characteristic of Omicron: NS9b_P10S, NS9b_E27-29del, NSP4_T492I, NSP5_P132H, NSP6_S106-107del, and NSP14_I42V. All subsequent Omicron variants compared to BA.1.X possessed NS3a_T223I, NSP1_S135R, NSP3_T24I, NSP3_G489S, NSP4_L264F, NSP4_T327I, NSP6_F108del, NSP13_R392C, and NSP15_T112I mutations. Distinctive features of BA.2.X and XBB were amino acid substitutions in NS6 ‒ D61L and NSP4 ‒ L438F. BA.5.X was characterized by additional mutations NS9b_D16G and NSP13_T127N, which were not found in other variants. CL.X, although being a sublineage of BA.5.X, did not possess these substitutions, but was characterized by the presence of NS8_F41C and NSP13_N102S. XBB variants, in addition to the above mutations, possessed a stop codon at the 8th position of NS8 and NSP1_K47R, while XBB.1.9.X possessed NS9b_I5T, NSP3_G1001S, NSP9_T35I and, in a small percentage of cases, NSP2_D449E (Fig. S3 c).

Regarding other structural genes of SARS-CoV-2, the Delta variant was not characterized by mutations in the E-protein, but substitutions in the matrix protein (M_I82T) and nucleocapsid protein (N_D63G, N_R203M, N_G215C, N_D377Y) were observed, and, as with most other amino acid substitutions, they were not found in Omicron. All Omicron variants circulating in Moscow were characterized by the following mutations: E_T9I, M_Q19E, M_A63T, N_P13L, N_E31-33del, N_R203K, N_G204R, but some genetic lineages had their own peculiarities. For example, BA.1.X carried an additional M_D3G mutation to those listed, the BA.2.X variant and all other Omicron variants considered were N_S413R, BA.5.X and CL.X carried M_D3N. CL.X also possessed the N_A218S substitution, and all XBBs possessed the E_T11A substitution (Fig. S3 b).

Viral load dynamics of the main circulating variants

Comparative analysis of viral load (equivalent - PCR threshold cycle, Ct) showed that among the Wuhan (B.1.X), Delta (B.1.617.2 + AY.X) and the first Omicron variants (BA.1/BA.2), the highest viral load was characteristic of the Delta variant with a median value of 23.85 (p < 0.0001, Fig. 3 a). Thus, an increase in viral load was observed when the Wuhan variant was replaced by genetic lineages belonging to the Delta variant (within the same “serotype”). The first Omicron variants were characterized by a reduced viral load compared to Delta (median 29.62, p < 0.0001), and further there was also a statistically significant small increase in Ct values for variants BA.2.X and BA.5.X compared to BA.1.X (Fig. 3 b). No statistically significant differences were found between BA.2.X and BA.5.X. The CL.X variant was characterized by a reduced Ct value compared to BA.5.X with median values of 35.13 and 27.27, respectively (p < 0.0001, Fig. 3 c). Furthermore, CL virus variants with the additional T883I mutation showed a higher viral load compared to variants without this mutation (p = 0.001, Fig. S1), but the viral load of XBB.1.X was not affected by the presence of this mutation (Fig. S1). There was a trend towards higher viral load in successive variants BA.5.X, CL.X, XBB.1.X and XBB.1.9.X, where the highest viral load was characterized for XBB.1.9.X (p < 0.01, Fig. 3 c).

 

Fig. 3. Comparative analysis of viral load in predominant SARS-CoV-2 variants in Moscow. a – viral load in Wuhan (B.1.X), Delta (B.1.617.2 + AY.X), and Omicron (BA.1/BA.2.X) variants; b – viral load in Omicron BA.1.X, BA.2.X, and BA.5.X variants; c – viral load in BA.5.X, CL.X, XBB.1.X, and XBB.1.9.X variants. The Y-axis denotes Ct values, while the X-axis represents SARS-CoV-2 variants and the number of samples.

Рис. 3. Сравнение вирусной нагрузки основных доминировавших в Москве вариантов SARS-CoV-2. a – вирусная нагрузка вариантов Ухань (B.1.X), Дельта (B.1.617.2 + AY.X) и Омикрон (BA.1/BA.2.X); б – вирусная нагрузка вариантов Омикрон BA.1.X, BA.2.X и BA.5.X; в – вирусная нагрузка вариантов BA.5.X, CL.X, XBB.1.X и XBB.1.9.X. По оси ординат отражены значения Ct, по оси абсцисс – варианты SARS-CoV-2 и количество образцов.

 

Discussion

During more than 4 years of the COVID-19 pandemic, an unprecedented amount of data, primarily genomic data, was generated that is being used both to understand evolution and to improve prophylactic and therapeutic agents. During this period, virus variants with different phenotypic characteristics, including infectivity, disease severity, and immune evasion, have emerged [28]. A comprehensive analysis of the obtained data on the pathogenesis of SARS-CoV-2 allows us to understand the underlying mechanisms (drivers) of evolution and to recognize the processes that generate this diversity to potentially predict possible future variants of the virus. As the main drivers of the ongoing pandemic, we considered two key aspects – viral load and the profile of mutations in the virus genome, which generally characterizes the observed transition: «increase in mutations – increase in viral load – emergence of a new serotype».

From the pandemic’s onset in 2020 until the early summer of 2023, Moscow witnessed the succession of several genetic variants of SARS-CoV-2, namely B.1.X, Delta (B.1.617.2 + AY.X), and Omicron. Within each of these, certain genetic lineages were dominant at the respective time periods: in 2020. – B.1.1.317, in 2021 – B.1.1.523, Alpha (B.1.1.7 + Q.4) and AY.122 [29], in 2022 – lineages BA.1.X, BA.2.X, BA.5.X, CL.X, in 2023 – recombinant XBB with the respective sublines XBB.1.X and XBB.1.9.X (Fig. 1 а). A consistent feature was a consistent increase in the number of nonsynonymous mutations in the genome of circulating SARS-CoV-2 variants, about half of which were in the S-protein (Fig. 1 b). Subsequent detailed consideration of the molecular-genetic properties of circulating variants of SARS-CoV-2, including the profile of mutations in the SARS-CoV-2 genome and viral load, allowed us to identify certain patterns that could account for the continued spread of new variants of the COVID-19 pathogen and the newly arising rises in morbidity and hospitalizations in Russia and worldwide. Thus, it was found that when switching from variants B.1.X to Delta (B.1.618.2 + AY.X), there was a decrease in viral load by 7.7 PCR cycles, and this phenomenon was previously also observed in other studies [30, 31]. This variant was also characterized by several mutations, including those in the receptor-binding domain of the Spike protein (Fig. 2, Fig. S3), which could potentially account for its increased viral load. In previous studies, the presence of the L452R substitution in the S-protein has been shown to correlate with evasion of neutralizing antibodies [32], increased affinity for the ACE2 receptor, and increased spike stability and viral infectivity, thereby contributing to increased viral replication [33]. Furthermore, mutations E156G, T478K and D614G also increased infectivity and viral affinity for the receptor, while other substitutions such as P681R and D950N promoted higher propagation rate due to more efficient cleavage of S1/S2 at the furin site [34, 35]. In addition to structural proteins, other mutations could also influence the spread of the Delta variant and increase its viral load, respectively. For example, the T492I mutation in NSP4 was shown to have a positive effect on virus replication [36]. Virus variants with the P323L amino acid substitution in the RdRp gene in vitro had a selective advantage over variants without it, and, in addition, the presence of L323 and S671 was associated with more efficient replication at reduced temperature in the upper respiratory tract [37]. Thus, additional mutations characteristic of Delta could directly result in a high viral load, which could determine the faster spread of this variant compared to the previously circulating variants [23, 38, 39].

The first Omicron BA.1/BA.2.X variants, which replaced the Delta variant, were characterized, on the contrary, by a reduced viral load with a higher rate of spread [40]. It is worth noting that Omicron is phylogenetically divergent from previously circulating variants and possesses more than 15 substitutions in the RBD alone (Fig. 2), not counting mutations in other genes (Fig. S3). In addition, this variant is serologically distinct from previously circulating variants, as evidenced by several studies. Whereas vaccination efficacy against Delta remained high [41] and most monoclonal antibodies neutralized the virus [42], in the case of Omicron, a decrease in vaccination efficacy and sensitivity to therapeutic and prophylactic drugs was observed [7, 8]. Thus, it can be assumed that Omicron spread actively at the first stages not due to increased viral load, but due to the fact that the immune system had not encountered such a virus serotype before [43]. Furthermore, mutations in nonstructural proteins (e.g., 203K/204R in the N-protein and ΔSGF/ΔLSG in NSP6) and S-protein (N501Y and H655Y) acquired during the evolution of the Omicron variant may also have contributed to its widespread spread due to its high replication rate and more efficient virus transmission [44–47]. Furthermore, the dramatic increase in the number of mutations in the RBD may be a consequence of their mutual influence, where the appearance of one substitution entailed the appearance of another, as in the case of Q498R and N501Y, which also affected ACE2 binding affinity and infectivity of the virus, respectively [48]. Since the second half of 2022, BA.5.X has become the dominant Omicron variant in Moscow, thereby replacing the previously circulating BA.2.X (Fig. 1 a). Comparative analysis showed that the viral load of genetic lineages BA.2.X and BA.5.X was statistically significantly different from BA.1.X with median Ct values of 25.1, 25.0, and 24.3, respectively, with no significant differences between BA.2.X and BA.5.X (Fig. 3 b). Although the results obtained partially diverge from the data of earlier studies [49–52], they correspond to the dynamics of the epidemic process in Moscow. Thus, with the emergence of the Omicron variant BA.1.X, an increase in the incidence was observed, whereas BA.2.X displaced BA.1.X without an accompanying increase in the number of new COVID-19 cases (Fig. 1 a). The active spread of BA.5.X coincided with the next (6th) rise in incidence, and, theoretically, this was possible not due to an increase in viral load, but due to the appearance of a number of nonsynonymous substitutions in this variant, in particular, two additional mutations in RBD L452R and F486V, the first amino acid substitution being also characteristic of the Delta variant (Fig. 2). It has been shown that due to these mutations there is an evasion of neutralizing antibodies formed after previous infection with the first Omicron variants and vaccination [53, 54]. Furthermore, BA.5.X compared to BA.1/2.X appears to have nonsynonymous substitutions in the ORF9b (D16G) and NSP13 (T127N) genes, which may promote virus spread through evasion from the host immune system [55]. Accordingly, from the Delta variant to BA.1.X and BA.5.X, the accumulation of mutations affecting the ability to evade immunity and penetrate the cell more efficiently was consistent, which allowed these variants to cause several waves of morbidity without a corresponding increase in viral load.

Toward the end of 2022, CL.X, which is a sublineage of BA.5.1.29, becomes one of the main variants in Moscow [56]. This variant was predominantly distributed in Russia and was characterized by several additional mutations compared to the parental lineage. These mutations include K444N and K150E in the S-protein, A218S in the nucleocapsid protein, F41C in ORF8, and N102S in NSP13 (Fig. S3). Along with the occurrence of the above mutations, an increase in viral load was observed, which may indicate a positive effect of these amino acid substitutions on CL.X replication. However, these mutations were not widespread in the general population of Omicron variant virus: their proportion ranged from 0.05% (for NS8_F41C) to 0.3% (for Spike_K444N) according to GISAID data as of September 7, 2023 [57]. Due to the mutation in RBD, the virus could theoretically avoid neutralizing antibodies [58], which could promote its replication; in addition, they could potentially affect ORF8 and its role in epigenetic regulation [59] or the helicase activity of NSP13 [60]. A detailed review of the genetic structure of this lineage also revealed the T883I mutation, which appeared simultaneously with the onset of CL circulation and was present in about 29% of all collected samples belonging to the CL lineage (Fig. S3 a). The presence of T883I was associated with a higher viral load compared to the variant without substitution (Fig. S1 a), and the proportion of the CL.1.2 variant characterized by this mutation began to increase with the introduction of XBB (Fig. S2). In addition, T883I was also found in XBB.1.X in about 12% of cases, and the proportion of this mutation also began to increase in parallel with the onset of XBB.1.9.X circulation (Fig. S2), but a similar effect in terms of increased viral load was not observed for XBB.1.X + T883I compared with the variant without the substitution (Fig. S1 b). Thus, the mere presence of this mutation did not contribute to an increase in viral load; however, it is possible that the combination of mutations in the S-protein of CL.X, including T883I, could have a positive effect on virus spread, but its proportion in the population was still insufficient.

For the latter dominant variants, XBB.1.X and XBB.1.9.X, a consistent increase in viral load was also observed along with the occurrence of a number of additional amino acid substitutions in both Spike and other proteins. XBB is a recombinant of BJ.1 (sub-variant BA.2.10) and BM.1.1.1.1 (sub-variant BA.2.75) with a recombination point in the S1 of S-protein around 445–460 amino acid positions [61]. Here we observe the first case among all previously dominant SARS-CoV-2 variants of virus adaptation through recombination. The acquired mutations V83A, R346T, L368I and N460K in the S-protein have been shown to result in increased hACE2 binding affinity, fusogenicity and infectivity of the XBB variant virus [62]. These findings are consistent with our results and explain the increased viral load of XBB.1.X compared to BA.5.X and CL.X. However, according to previous studies, the presence of a stop codon at the beginning of the ORF8 gene should have had a negative effect on viral replication with a corresponding decrease in viral copies [59], but our results show the opposite effect. Perhaps the presence of several other substitutions in non-structural genes (Fig. 3) compensates for the absence of full-length NS8. The increased viral load of the XBB.1.9.X variant compared to the previously circulating XBB.1.X (Fig. 1, Fig. 3 c) may have been due to the F486P substitutions in the S-protein present in about 27% of cases, NSP3_G1001S, NSP9_T35I, NS9b_I5T. The presence of F486P was predominantly characteristic of the XBB.1.5 variant, which spread rapidly in several countries due to increased transmissibility [63] and, in agreement with earlier studies, the presence of this substitution correlated with increased affinity to hACE2 [64].

Conclusion

Molecular genetic monitoring of SARS-CoV-2 variants in Moscow since the beginning of the pandemic has allowed not only to characterize the profile of the main genetic lineages of the virus, but also to reveal several regularities that could account for the continued spread of new variants of the virus. Thus, the change of genetic variants of SARS-CoV-2 in Moscow for more than 3 years from the beginning of the pandemic was accompanied by a gradual increase in the number of mutations, where their sharp increase was characteristic of the Omicron variant that emerged at the end of 2021. Furthermore, for several successive variants, for example, from Wuhan to Delta or from BA.5 .X to XBB.1.9.X, there was a consistent increase in viral load along with the appearance of new amino acid substitutions in both the Spike protein and nonstructural proteins. In the remaining cases where this trend was not observed, there was an accumulation of several mutations that promoted evasion of previously established immunity and increased infectivity of the virus, which determined their further spread. These results can be used for modeling and predicting future variants, as they clearly demonstrate the transition: «increase in mutations ‒ increase in viral load ‒ emergence of a new serotype».

Funding. The study was funded by the Ministry of Health of the Russian Federation (Project No. 123031400022-0, titled Investigation of SARS-CoV-2 variability in relation to the biological risks associated with the reduced efficacy of therapy and prevention used during the COVID-19 pandemic).

Acknowledgement. We are grateful to all GISAID submitted and originating laboratories for access to SARS-CoV-2 sequencing data. We are grateful to Anastasia A. Zakharova and Timofey A. Remizov for their technical assistance in the implementation of the project.

Conflict of interest. The authors declare no apparent or potential conflicts of interest related to the publication of this article.

Ethics approval. The study was conducted with the informed consent of the patients. The research protocol was approved by the Ethics Committee of the National Research Centre for Epidemiology and Microbiology Named after Honorary Academician N.F. Gamaleya (Protocol No. 14 dated 29 September 2021).

About the authors

Daria D. Kustova

National Research Centre for Epidemiology and Microbiology Named after Honorary Academician N.F. Gamaleya of the Ministry of Health of the Russian Federation; Federal State Budgetary Educational Institution of Higher Education Lomonosov Moscow State University

Email: kustovadaria@yandex.ru
ORCID iD: 0000-0002-8382-275X

junior researcher, Laboratory of mechanisms of population variability of pathogenic microorganisms; PhD student, Department of Virology, Faculty of Biology

Россия, 123098, Moscow; 119991, Moscow

Andrei A. Pochtovyi

National Research Centre for Epidemiology and Microbiology Named after Honorary Academician N.F. Gamaleya of the Ministry of Health of the Russian Federation; Federal State Budgetary Educational Institution of Higher Education Lomonosov Moscow State University; I.M. Sechenov First Moscow State Medical University of the Ministry of Health of the Russian Federation (Sechenov University)

Author for correspondence.
Email: a.pochtovyy@gamaleya.org
ORCID iD: 0000-0003-1107-9351

Cand. Sci. (Biol.), senior researcher, Laboratory of mechanisms of population variability of pathogenic microorganisms

Россия, 123098, Moscow; 119991, Moscow; 119435, Moscow

Olga G. Shpakova

Moscow Healthcare Department

Email: shpakovaog@dcli.ru

Head of the laboratory of the Moscow Scientific and Practical Center for Laboratory Research

Россия, 127006, Moscow

Irina A. Shtinova

Moscow Healthcare Department

Email: shtinovaia@dcli.ru

Head of Laboratory Center of the Moscow Scientific and Practical Center for Laboratory Research

Россия, 127006, Moscow

Nadezhda A. Kuznetsova

National Research Centre for Epidemiology and Microbiology Named after Honorary Academician N.F. Gamaleya of the Ministry of Health of the Russian Federation

Email: nadyakuznetsova0@yandex.ru
ORCID iD: 0000-0002-7399-7628

Cand. Sci. (Biol.), senior researcher, Laboratory of mechanisms of population variability of pathogenic microorganisms

Россия, 123098, Moscow

Denis A. Kleimenov

National Research Centre for Epidemiology and Microbiology Named after Honorary Academician N.F. Gamaleya of the Ministry of Health of the Russian Federation

Email: mne10000let@yandex.ru
ORCID iD: 0000-0001-9422-7238

Cand. Sci. (Med.), Head, Laboratory of translational biomedicine

Россия, 123098, Moscow

Andrey G. Komarov

Moscow Healthcare Department

Email: komarovag@zdrav.mos.ru
ORCID iD: 0009-0000-8597-7125

Cand. Sci. (Med.), Head of the Moscow Scientific and Practical Center for Laboratory Research

Россия, 127006, Moscow

Vladimir A. Gushchin

National Research Centre for Epidemiology and Microbiology Named after Honorary Academician N.F. Gamaleya of the Ministry of Health of the Russian Federation; Federal State Budgetary Educational Institution of Higher Education Lomonosov Moscow State University; I.M. Sechenov First Moscow State Medical University of the Ministry of Health of the Russian Federation (Sechenov University)

Email: wowaniada@yandex.ru
ORCID iD: 0000-0002-9397-3762

Dr. Sci. (Biology), Head, Laboratory of mechanisms of population variability of pathogenic microorganisms, Reference center for coronavirus infection; senior researcher, Department of virology, Biological faculty

Россия, 123098, Moscow; 119991, Moscow; 119435, Moscow

References

Supplementary files