Omicron-BA.1 Dispersion Rates in Mexico Varied According to the Regional Epidemic Patterns and the Diversity of Local Delta Subvariants

Purpose: The Omicron subvariant BA.1 of SARS-CoV-2 was first detected in November 2021 and quickly spread worldwide, displacing the Delta variant. In this work, a characterization of the spread of this variant in Mexico is presented. Methods: The time to fixation of BA.1, the diversity of Delta sublineages, the population density, and the level of virus circulation during the inter-wave interval were determined to analyze differences in BA.1 spread. Results: BA.1 began spreading during the first week of December 2021 and became dominant in the next three weeks, causing the fourth COVID-19 epidemiological surge in Mexico. Unlike previous variants, BA.1 did not exhibit a geographically distinct circulation pattern. However, a regional difference in the speed of the replacement of the Delta variant was observed. Conclusions: Viral diversity and the relative abundance of the virus in a particular area around the time of the introduction of a new lineage seem to have influenced the spread dynamics, in addition to population density. Nonetheless, if there is a significant difference in the fitness of the variants, or if the time allowed for the competition is sufficiently long, it seems the fitter virus will eventually become dominant, as observed in the eventual dominance of the BA.1.x variant in Mexico.


Introduction
SARS-CoV-2 emerged in late 2019 and caused a global pandemic. During this time, multiple lineages of this virus have emerged. In late 2020, the appearance of viral variants associated with increased transmissibility or virulence, detrimental changes in epidemiology or clinical disease presentation, or decreased effectiveness of diagnostics, vaccines, or therapeutic methods were reported [1]. These variants were subsequently classified as Variants of Concern (VOCs) by the World Health Organization (WHO) and have become regionally or globally dominant. So far, there have been five VOCs, dubbed Alpha [2], Beta [3], Gamma [4], Delta [5], and Omicron [6]. The latter was declared a VOC on November 26, 2021, and was initially designated B.1.1.529 in the Pangolin classification. Omicron was detected a few days prior in South Africa and Botswana and made the headlines for its unprecedented ability to rapidly outcompete other circulating variants [6]. By the beginning of December 2021, it had spread to several European countries, causing epidemic surges, and soon became dominant worldwide [7].
As sequences accumulated, phylogenetic analyses showed that the Omicron variant consisted of two main separate clades, designated as BA.1 and BA.2 in the Pangolin nomenclature system. BA.1, labeled clade 21 K by Nextclade, was the first of the two clades to spread, and it caused epidemic outbreaks throughout the world between December 2021 and January 2022. As of March 30, 2022, the BA.1 lineage had diversified into fifty-four sublineages (hereafter referred to as BA.1.x). Each one has a distinct mutation pattern and differences in geographical distribution, although no distinctive phenotypic variation has been reported for any of them.
BA.1 is characterized by forty-four nonsynonymous substitutions, ten synonymous changes, and seven deletions (ranging between one and three amino acids), compared to the reference strain. This large number of mutations had no precedence in all the described SARS-CoV-2 variants and posed questions regarding its possible origin. Notably, the gene coding for the spike protein (S), which mediates virus cell entry via ACE2 recognition, harbors thirty nonsynonymous substitutions, three deletions, and one insertion. The accumulation of genetic changes in S has been related to an increase in viral transmission rate compared to the original strain [8,9]. In the case of Omicron, it also exhibits a transmission advantage against previous VOCs [10]. For instance, in Denmark and Italy, it has been estimated that BA.1 had an effective reproduction number (Rt) around two to three times greater than the Delta variant [11,12]. The large number of mutations in S resulted in an enhanced antigenic escape, reflected by the high number of infections among the vaccinated population and a sharp decline in the neutralization titers of reported monoclonal antibodies, as well as in the sera of vaccinated subjects [13,14].
Additionally, Omicron BA.1.x transmission advantage seems to be associated with a preference towards a TMPRSS2-independent endosomal entry pathway rather than membrane fusion at the cell surface, which has been associated with a preferential infection of the upper respiratory tract resulting in a reduction in virulence, as described in clinical reports from multiple countries [15][16][17][18]. Furthermore, in hamster and murine infection models, Omicron exhibited reduced virulence with respect to previous variants [19,20], and in vitro assays demonstrated a limited activation of the NF-κB pathway, which may lead to lower inflammation [21].
By the first week of December 2021, when the Omicron variant began expanding in Mexico, the country had experienced three epidemic surges (peaks). The third outbreak, which occurred during the summer of 2021, was caused by the Delta variant, resulting in the largest number of cases up to that date. However, since vaccination of the most vulnerable population had already begun, in the third surge, the number of hospitalizations and deaths was lower than during the second epidemic peak [22]. Early in December, vaccination efforts reached 60% of the adult population with a complete scheme, and a booster dose campaign had started, aimed at people over 60 years of age in light of the imminent arrival of the Omicron variant.
Given the large number of variants that have emerged from the SARS-CoV-2 virus, global genomic surveillance is required to monitor its evolution and to detect new variants that may impact clinical outcomes. This study examines the entry and spread of Omicron BA.1.x, which caused the fourth pandemic wave in Mexico, describing the geographical and temporal dispersion patterns and the factors contributing to its varying spread rates in Mexico.

Epidemiological Information
National epidemiological data derived from the governmental SARS-CoV-2 surveillance program in Mexico is made available by the General Epidemiological Directorate (Dirección General de Epidemiología) through the Secretaría de Salud website (https: //www.gob.mx/salud/documentos/datos-abiertos-152127), which was accessed on 24 April 2022. Confirmed positive cases of COVID-19 and the associated deaths were collated per day, from 25 March 2020 to 8 April 2022, for the epidemiological analyses. The information was plotted at national and regional levels using a 7-day rolling average on a sliding window with a step of one day. Epidemiological data processing and graphs were created using R v.4.2.1 [23].

Sequences Data Source
All sequences and metadata from the Omicron variant available in the Global Initiative on Sharing All Influenza Data (GISAID) database (accessed on 18 May 2022) were downloaded with a submission dateline on 31 March 2022. In total, the set had 13,562 sequences from Mexico and 3,191,009 from the rest of the world. Regarding Mexican sequences, 5752 (42.4%) were generated by the National Institute for Genomic Medicine (INMEGEN), 4097 (30.2%) by the Mexican Consortium for Genomic Surveillance (CoViGen-Mex), 2851 (20.0%) by the Institute for Diagnostics and Epidemiological Reference (InDRE), and the remaining (6.4%) by other public and private institutions. The information on these Mexican genomes is included in Table S1 (Supplementary Materials). Regarding the CoViGen-Mex sequences, samples were taken and sequenced, then consensus genomes were generated as described in Taboada et al., 2022 [24].
To evaluate the displacement of the Delta variant by BA.1.x lineages, sequences and metadata from this variant prior to and during the transition period (November 2021 to January 2022) were also retrieved from GISAID; in total, 5489 Mexican sequences of the Delta variant were obtained (Table S2), with 2590 being generated by CoViGen-Mex.
All sequences and metadata are available through GISAID EPI_SET_220927gw and genome sequences from CoViGEn-Mex were also being deposited in GenBank (Tables S1 and S2).

Lineage and Mutation Frequencies
Lineage clade assignment was carried out using Pangolin PLEARN v1.8 [25] and mutation annotation was done using the NextClade tool v.2.0.0-beta.3 alignment [26]. Variant calling was used to determine the genetic diversity and frequency of BA.1.x sublineages in Mexico and the rest of the world. Regarding mutation analysis, all sequences containing more than 5% ambiguous nucleotides (Ns) were excluded from the analyses, resulting in 12,304 genomes from Mexico. to ensure that every state and month were represented. Additionally, 740 sequences from other countries were included with a distribution based on the relative prevalence of the three sublineages on each continent. Sequences were aligned using MAFFT v.7 with the addfragments option [27], and a maximum likelihood tree was built using iqtree2 [28] and then time-scaled with LSD2 [29]. The resulting phylogeny was visualized in R using ggtree [30]. Similarly, a haplotype network was built, since haplotypes more accurately reflect epidemiological clustering over short periods of time. Thus, using the same alignment employed for the tree, PopART v1.7 [31] was used to generate a Templeton-Crandall-Sing (TCS) haplotype network [32]. A Spearman correlation analysis was performed to analyze differences in lineage distribution between regions. Heatmaps were generated to show the abundance of lineages throughout time, using the pandas v.1.3.5, seaborn v.0.11.0, and matplotlib v.3.5.2 libraries of Python v.3.7.6. Finally, using the Shannon's index, the diversity of Delta lineages prior to the entry of BA.1.x in each region was determined. The significance of any differences in the Shannon diversity between regions was evaluated using the Wilcoxon signed-rank test and its associated p-value with an alpha = 0.05.

BA.1.x Entry into Mexico
According to official records from the Mexican government, the first SARS-CoV-2positive case was registered in Mexico on 28    Notably, the fourth wave (W5-Dec2021 to W08-Feb2022) accumulated more cases than the previous waves in a shorter period, with maximum daily records 8.5, 3.9, and 3.2 times greater than those registered during 1st, 2nd, and 3rd waves, respectively. In contrast, reported deaths attributed to COVID-19 have decreased since the second wave, Notably, the fourth wave (W5-Dec2021 to W08-Feb2022) accumulated more cases than the previous waves in a shorter period, with maximum daily records 8.5, 3.9, and 3.2 times greater than those registered during 1st, 2nd, and 3rd waves, respectively. In contrast, reported deaths attributed to COVID-19 have decreased since the second wave, having a case fatality rate of 11.4%, 8.7%, 4.1%, and 1.2% for waves 1st to 4th, respectively. This difference might be due to an increase in the population immunity, as a consequence of the protection acquired from prior natural SARS-CoV-2 infections and the vaccination campaign, which had covered most adults by the end of the third wave (complete scheme) and included a booster dose for the elderly by late 2021. eages circulated in Mexico with a prevalence higher than 1% (Table 1); however, the most abundant three accounted for 91.5% of all sequenced genomes. As shown in Table 1, BA.1.1 accounted for more than 50% of the sequenced genomes, while BA.1 and BA.1.15 had a prevalence of around 17% each. In contrast, some unique sublineages circulated predominantly in Mexico during the Delta wave compared to the rest of the world [18]; these three Omicron sublineages circulated at roughly the same frequencies in the USA. Interestingly, BA.1.1 accounted for around 50% of the sequences in North and South America, except for Canada and Brazil, while it exhibited a low prevalence in the rest of the world. On the other hand, BA.1 was more prevalent in Canada and Brazil (about 38%) than in Mexico and the USA (12% and 16%, respectively). Whole-genome mutation analysis was carried out using only high-quality genomes (having <3% ambiguous bases) of the six BA.1.x lineages with a prevalence higher than 1% in Mexico (Table 1). We identified 7144 distinct nucleotide changes, but only 110 (1.5%) were found in more than 1% of the genomes. Regarding nonsynonymous changes, 3939 were identified, 84 (2.1%) of which had a prevalence higher than 1%. These changes include the 46 amino acid signature substitutions and 7 amino acid deletions (between 1 to 3 aa in length) previously reported in BA.1.x, which were identified in over 90% of the Mexican isolates, as well as a 3-aa insertion that was present in over 40% of sequences (Table S3) Only two previously unreported amino acid substitutions were identified in these lineages, with a frequency greater than 5%.
A phylogenetic reconstruction and a haplotype network were built to determine the presence of BA.1.x Mexican haplotypes among the circulating Omicron viruses ( Figure S3). To this end, only the three most prevalent sublineages, BA.1, BA.1.1, and BA.1.15, were included in the analysis. Interestingly, these sublineages did not form distinct clusters in the network and were highly polytomic in the phylogeny ( Figure S3), probably due to their low genetic diversity. The Mexican sequences were interspersed with those from abroad, similarly pointing to the limited diversification of these lineages. Moreover, when analyzed by geographical region, Mexico showed no difference in sublineage genetic diversity ( Figure S3). This limited sequence diversity could be related to the rapid global and national spread of the BA.1.x variant.

BA.1.x Entry in Mexico Was Asynchronous
A series of maps were built to characterize the dynamics of the BA.1.x spread in Mexico. Figure 2 shows that at epidemiological week 49, all states were dominated by the Delta variant (blue). However, the following week, some states showed a larger proportion of BA.1.x sequences (red). During the last week of December 2021, most states in Central and South Mexico and many in the north had a majority of BA.1.x sequences. By the first week of 2022, the national prevalence of BA.1.x exceeded 50%, and during January its prevalence continued to rise. Interestingly, in some northern states, the replacement of the Delta variant was slower than in others.

BA.1.x Entry in Mexico Was Asynchronous
A series of maps were built to characterize the dynamics of the BA.1.x spread in Mexico. Figure 2 shows that at epidemiological week 49, all states were dominated by the Delta variant (blue). However, the following week, some states showed a larger proportion of BA.1.x sequences (red). During the last week of December 2021, most states in Central and South Mexico and many in the north had a majority of BA.1.x sequences. By the first week of 2022, the national prevalence of BA.1.x exceeded 50%, and during January its prevalence continued to rise. Interestingly, in some northern states, the replacement of the Delta variant was slower than in others. To further characterize the replacement of Delta by BA.1.x throughout the country, we compared four geographical regions, which were determined by considering the similarity of weekly lineage abundance in the respective states, as described in the Methods section. The weekly relative frequency of BA.1.x lineages was plotted against the epidemiological week to compare their dynamics of dispersion and prevalence by region ( Figure  3A). For clarity, the number of weeks required to reach a given relative frequency in each area was plotted, beginning with the first week of detection ( Figure 3B). A limitation of this analysis was the low number of sequences that were available from the SE region during the initial weeks of the study, so that by the time it was first detected, its prevalence had already ramped up to 50%. To overcome this limitation, the analysis was carried out by comparing the number of weeks it took BA.1.x to increase their prevalence from 50% to 100% in each region. Figure 3B shows that the increase from 50% to 100% required four and six weeks in the SE and CS regions, respectively, but instead took eight weeks in the N and CN regions. In addition to the shorter period it took BA.1.x to dominate the CS and SE, the variant was detected earlier in these regions ( Figure 3A). To further characterize the replacement of Delta by BA.1.x throughout the country, we compared four geographical regions, which were determined by considering the similarity of weekly lineage abundance in the respective states, as described in the Methods section. The weekly relative frequency of BA.1.x lineages was plotted against the epidemiological week to compare their dynamics of dispersion and prevalence by region ( Figure 3A). For clarity, the number of weeks required to reach a given relative frequency in each area was plotted, beginning with the first week of detection ( Figure 3B). A limitation of this analysis was the low number of sequences that were available from the SE region during the initial weeks of the study, so that by the time it was first detected, its prevalence had already ramped up to 50%. To overcome this limitation, the analysis was carried out by comparing the number of weeks it took BA.1.x to increase their prevalence from 50% to 100% in each region. Figure 3B shows that the increase from 50% to 100% required four and six weeks in the SE and CS regions, respectively, but instead took eight weeks in the N and CN regions. In addition to the shorter period it took BA.1.x to dominate the CS and SE, the variant was detected earlier in these regions ( Figure 3A).

Factors Affecting BA.1.x Dispersion Dynamics by Region
To understand the differences in the dynamics of BA.1.x dispersion, we analyzed and plotted the sublineage diversity of Delta and BA.1.x per Mexican region from W39-Oct2021 to W05-Feb2022 ( Figure 4). Interestingly, the diversity of Delta sublineages in the weeks prior to the entry of BA.1.x also differed by region (W44-Nov to W49-Dec). In the N region, the observed diversity was higher, with six lineages reaching at least 10% of prevalence. In comparison, the CN and SE regions had four and three lineages above 10%, respectively, whereas the CS region exhibited the lowest diversity, clearly dominated by a single sublineage. These results were in agreement with the differences between Shannon diversity indices, which considered the abundance of Delta sublineages from W39-Oct to W48-Nov 2021 to be significantly higher in the N (H = 2.

Factors Affecting BA.1.x Dispersion Dynamics by Region
To understand the differences in the dynamics of BA.1.x dispersion, we analyzed and plotted the sublineage diversity of Delta and BA.1.x per Mexican region from W39-Oct2021 to W05-Feb2022 ( Figure 4). Interestingly, the diversity of Delta sublineages in the weeks prior to the entry of BA.1.x also differed by region (W44-Nov to W49-Dec). In the N region, the observed diversity was higher, with six lineages reaching at least 10% of prevalence. In comparison, the CN and SE regions had four and three lineages above 10%, respectively, whereas the CS region exhibited the lowest diversity, clearly dominated by a single sublineage. These results were in agreement with the differences between Shannon diversity indices, which considered the abundance of Delta sublineages from W39-Oct to W48-Nov 2021 to be significantly higher in the N (H = 2.  A SARS-CoV-2 high lineage diversity has been associated with population density in Mexico [9]. This factor, rather than lineage diversity per se, could explain the differences in BA.1.x dispersion. Accordingly, the N region of Mexico has the lowest population density and the slowest BA.1.x dispersion (see Figure S4). However, the region with the second lowest population density and the fastest spread of BA.1.x was the SE region, suggesting that population density alone cannot determine lineage diversity nor explain disparities in BA.1.x spread rate.
Notably, in the CS and the SE, the proportion of Delta sublineages that circulated during the fourth wave [24] was practically the same as that of the inter-wave (W42-Oct to W50-Dec) period. In contrast, in the N region, and to a lesser extent in the CN, the introduction and growth of additional Delta sublineages occurred during the inter-wave interval, such as AY.3, AY.103, and AY.113 (Figure 4).
To investigate other factors contributing to the variations in the BA.1.x spread, aside from the diversity of Delta lineages before the onset of BA.1.x and population density in the regions, the epidemic curves of these regions were examined ( Figure 5). The epidemic curve of the N region shows an active inter-wave interval ( Figure 5A). The maximum number of cases in N during the Delta wave was 1880 per day, and during the inter-wave period, cases persisted at an average of 977 per day, corresponding to 52% of the Delta wave's peak. Additionally, infections in this region began to increase during epidemiological W47-Nov and peaked at 1170 daily cases during W48-Nov. Throughout this period, all sequences identified in the N region corresponded to the Delta variant (see Figure 4A). In contrast, in the CN region, SARS-CoV-2 circulation during the Delta-BA.1.x inter-wave period was roughly 13% of the Delta peak, whereas it was 10% and 6% in the CS and SE regions, respectively ( Figure 5B-D).
Additionally, in the N region, the BA.1.x peak was 5.9 times higher than Delta's, while in the CN region it was 3.3 times higher. The increase in the CS and SE regions was 2.9 and 2.5, respectively, resulting in a smaller disparity. Taking all data together, a higher Delta peak (relative to the BA.1.x peak), lesser virus circulation during the inter-wave period, and/or a lower Delta variant diversity were associated with a faster replacement of Delta by BA.1.   period, cases persisted at an average of 977 per day, corresponding to 52% of the Delta wave's peak. Additionally, infections in this region began to increase during epidemiological W47-Nov and peaked at 1170 daily cases during W48-Nov. Throughout this period, all sequences identified in the N region corresponded to the Delta variant (see Figure 4A). In contrast, in the CN region, SARS-CoV-2 circulation during the Delta-BA.1.x inter-wave period was roughly 13% of the Delta peak, whereas it was 10% and 6% in the CS and SE regions, respectively ( Figure 5B-D).

Discussion
In this study, we analyzed the SARS-CoV-2 genome sequences isolated in Mexico around the time of the introduction of the Omicron-BA.1 variant to the country to better understand its spread during the pandemic's fourth wave in Mexico. The first national case was reported by late November 2021, quickly becoming the dominant variant by mid-December. The first significant observation was the rapid increase in the number of people infected, which represented the highest peak in the country thus far during the pandemic, occurring in mid-January 2022 with more than 62,000 daily infections. Although this increase coincided with Christmas and New Year celebrations, which may have had an impact on the number of cases [33], the second wave that was caused by variant B.1.1.519, which began in December 2020, had a significantly lower peak, with a maximum of 21,000 cases per day [9].
In contrast, and as was reported in other countries, this variant had a relatively minor impact on hospitalizations and deaths. This phenomenon has been attributed to high vaccination rates and a general reduction in the susceptible population, despite the ability of the BA.1 virus to evade neutralization, and to Omicron's potentially lower virulence [7,[34][35][36].
The spread of BA.1.x sublineages in Mexico did not result in the acquisition of specific new mutations nor the clustering of sequences by sublineage or geographical region, in contrast to previous reports on the spread of Alpha and Delta VOCs, which described the emergence of local specific mutations and sublineages, as well as a solid geographical association with those changes [24,37]. The low diversity of Omicron sublineages was probably associated with its unprecedentedly rapid spread across the country. This obser-vation is consistent with a previous report describing BA.1.x entry and spread in Mexico City [38], and with descriptions from other countries. For instance, in an analysis of South American Omicron sequences, the authors divided them into clusters using phylogenetic methods; however, these clusters contained sequences from different sublineages and countries, suggesting an overall low diversity [39]. Additionally, in Finland, a phylogenetic analysis of BA.1.x sequences reported multiple singletons and low clustering of Finnish sequences [40].
Despite the reported transmission advantage of BA.1 over Delta [6,7,12], this variant propagated more slowly in the N and CN of Mexico than in the CS and SE regions. The slower spread rate may be attributed to a higher diversity of Delta sublineages and epidemic activity during the BA.1.x inter-wave interval. Notably, in the N region, the number of active SARS-CoV-2 cases reported during the Delta-BA.1.x inter-wave period was the highest in the country, representing 52% of the maximum peak of the Delta wave in that region, indicating continued transmission during the fall of 2021. The active circulation of the virus may have allowed the introduction of new Delta sublineages, sustained by the low population density and a low proportion of people with prior Delta infection. Consequently, this region exhibited the highest diversity of Delta sublineages before the onset of Omicron.
Conversely, CS and SE had the lowest number of cases in their Delta-BA.1.x inter-wave period and the smallest differences in the ratio between the peaks of the third and fourth waves. Additionally, these regions had the lowest diversity of Delta lineages at the onset of Omicron. Furthermore, in contrast to the N region, no change was observed between the sublineages circulating previously during the Delta peak and in the inter-wave period [24]. The CN region seems to be a transition zone between the N and the CS and SE regions, as it shares a slower replacement and a higher diversity of Delta with the N region. Still, its epidemic curve was comparable to that of the CS and SE regions. To our knowledge, the association between the epidemiological and genomic profile of the viruses in a region with the rate of Omicron dispersion has not been reported.
Geographical and prevalence differences in the circulation of distinct lineages in Mexico have been previously reported [9,24,37]. The prevalence of B.1.1.519, which dominated the center and south of the country between February and May 2021, was lower in the north [37]. Conversely, in the north of Mexico, the Alpha variant circulated at a higher prevalence than in the rest of the country during roughly the same period. Whereas the Alpha variant had a proven transmission advantage and managed to displace other lineages in many countries, the B.1.1.519 lineage was only dominant in Mexico, suggesting that its fitness advantage was less than that of Alpha. This observation indicates that fitness alone is insufficient to predict the spread of a SARS-CoV-2 variant, which is, in fact, a complex phenomenon. In this study, Delta lineage diversity seems to have influenced the rate of spread of Omicron-BA.1.x in Mexico. However, the sublineages that circulated in these regions could also play an important role in this regionally delayed replacement. Sublineages AY.3, AY.103, and AY.113 circulated at higher frequencies in the N and CN, which may have interfered with the arrival of BA.1.x. In particular, AY.3 and AY.103 have been associated with breakthrough infections in vaccinated people, with AY.3 having caused asymptomatic cases with high viral loads [41] and AY.103 having a breakthrough rate of 41%, which is higher than other Delta sublineages [42]. Therefore, these viruses could continue circulating during the inter-wave period and possibly delay the entry of BA.1. A molecular dynamic simulation has shown that AY.3 may exhibit a higher affinity for ACE2 than Omicron-BA.1 [43], further suggesting that some Delta sublineages may have higher fitness and thus be maintained in the population longer.

Conclusions
The genetic homogeneous dispersion of BA.1.x lineages in Mexico contrasts with the accumulation of specific mutations reported in variants of earlier waves. However, it was found that the rate of spread of the BA.1.x lineages varied by region. This variation was