Genetic diversity and spread of COVID-19 in Morocco during the first 3 years, including a unique local variant

The sequence data used in this investigation were extracted from the GISAID EpiCov database. All the sequences analysed in this study are available under the accession number (GISAID identifier): EPI_SET_240225ry (File S1, available in the online version of this article). Supplementary Material is provided alongside the article in PDF and Excel formats. Access to the Supplementary Material is available on Figshare via the following DOI link: https://doi.org/10.6084/m9.figshare.25748970.v1↗ [1].

The results presented in this study are derived from the analysis of 2274 genomic sequences of metadata available on: GISAID Identifier: EPI_SET_240225ry. DOI: 10.55876/gis8.240225ry.

Appearing in 2019 in Wuhan, in the Hubei province, China, COVID-19 is one of the most devastating pandemics of this century, involving millions of cases and deaths worldwide. In February 2020, the first case of COVID-19 was reported in Africa, in Egypt, and the first case in Morocco was reported on 2 March 2020 [2]. Moreover, in Morocco, a national SARS-CoV-2 genomic surveillance consortium has been set up to sequence Moroccan isolates in order to inform the health authorities about circulating variants. As of 6 August 2023, more than 769 million cases of infection and over 6.9 million deaths have been reported worldwide, according to the World Health Organization (WHO) [3]. In Morocco, during the last quarter of 2023, over 1 278 055 confirmed cases and 16 298 deaths were recorded [4].

The SARS-CoV-2, an RNA virus of the coronavirus family responsible for the global COVID-19 pandemic, has a genome length of approximately 29.9 kilobases (kb) [5]. The genome consists of 15 ORFs, with two important polyproteins involved in replication: pp1a and pp1ab (ORF1ab). These polyproteins are cleaved to produce 16 non-structural proteins. Additionally, the genome codes for 11 accessory proteins and four structural proteins, including the Spike glycoprotein, envelope, membrane, and nucleocapsid proteins [57]. Generally, within the SARS-CoV-2 genome, mutations can progressively accumulate over time, in particular, key mutations in structural proteins, of which the Spike glycoprotein is the most closely monitored. Moreover, some mutations in the SARS-CoV-2 genome can significantly impact its fundamental properties. These mutations can increase contagiousness, enhance the severity and virulence of the infection, reduce the effectiveness of drugs and vaccines, allow the virus to evade the immune system, cause reinfections, and escape molecular detection during laboratory diagnostics. Therefore, variants that exhibited the aforementioned characteristics are classified either as variants of concern (VOC), either as variants under monitoring (VUM) or variants of interest (VOI). However, these variants required close monitoring [89].

This study aimed to retrospectively analyse data from sequences derived from Moroccan isolates during the COVID-19 pandemic, in order to better understand the genomic epidemiology of SARS-CoV-2 in Morocco. The main aim of this study was to analyse the genetic diversity of SARS-CoV-2 and identify distinct lineages, as well as to assess their evolution during the pandemic in Morocco, using genomic epidemiology approaches. Furthermore, several key mutations in the functional proteins across different viral lineages were highlighted along with an analysis of the genetic relationships amongst these strains to better understand their evolutionary pathways.

The study conducted will provide access to in-depth data analysis with detailed large-scale analysis. Furthermore, this study will also shed light on the genetic trajectory of SARS-CoV-2, specifically the appearance and distribution of variants during the 3 years after the pandemic's onset in Morocco. This study will extend our understanding of the elements affecting the genetic variation of SARS-CoV-2 in Morocco, such as mutation rates and transmission dynamics. Finally, it will provide a strong basis for understanding the virus's evolution, facilitating the development of more precise preventive measures in future research. However, expanding the analysis by incorporating global GISAID sequences would contribute to a better development of preventative measures.

This study stands out as one of the first to investigate the complete genomic divergence of 2274 sequences of severe acute respiratory syndrome coronavirus 2 in Morocco over the 3 years following the COVID-19 pandemic. Furthermore, this study highlights the impact of SARS-CoV-2 evolution on public health by identifying several distinct mutational signatures that characterize notable variants, thus providing a solid basis for future research.

During the COVID-19 pandemic, the year 2020 witnessed a turning point in the evolution of SARS-CoV-2 with the emergence of the D614G mutation [17]. Some mutations in the genome of SARS-CoV-2 have been recognized as determining mutations in variations in transmissibility (more or less contagious), antigenicity variation (immune evasion capacity and/or risk of reinfections), as well as variations in pathogenicity (virulence variation with more or less severe events) [1719]. Mutations in the RBD/RBM region could impact how the virus interacts with the ACE2 receptor, potentially influencing the virus's transmissibility and could assist the virus to escape immune recognition. Moreover, mutations near the furin cleavage site (S1/S2) could accelerate the Spike protein cleavage, enhancing virus infectivity (Fig. 4b) [1719].

In Morocco in 2020, dominant lineages such as, B.1, B.1.528, B.1.177, and B.1.1 with the D614G mutation were prevalent, mainly in clade G. All subsequent SARS-CoV-2 lineages/variants in Morocco are derived from lineage B.1 (clade 20A). These lineages emerged approximately at the same time, indicating multiple introduction events in Morocco, except for the lineage « B.1.528 », which is closely related to clade 20A. This lineage formed a distinct cluster nationally and carried notable mutations in Spike glycoprotein, including D614G, V826L, and A892S (Fig. 4b and Table 5) [20]. These amino acid substitutions, located outside the predicted immunological epitope region, could affect transmissibility and immune evasion [20]. Furthermore, this study highlighted the P1263L mutation, located in the C-terminal region, within B.1.528 sequences. This mutation has also been reported in certain lineages in Europe and Oceania [21]. Interestingly, previous studies have suggested that the P1263L mutation, could potentially alter the structure and function of the SARS-CoV-2 Spike protein, leading to decreased infectivity of the virus [2122]. It should also be noted that non-synonymous mutations such as F2V, V3G, F4L, V6F, and H146Q in the B.1.528 sequences are located within the region of the predicted epitopes [23]. Previous studies have reported that these mutations could contribute to increased transmissibility and pathogenicity of SARS-CoV-2 [23]. Moreover, their location within epitopes could potentially confer to the virus the ability to escape the immune system [2122]. Consequently, all these key mutations in the spike glycoprotein, particularly those associated with the D614G mutation present in the sequences of lineage B.1.528 (Fig. 7b and Table 5), played a crucial role in the emergence of this lineage within clade 20A. These mutations distinguish it from other lineages within the same clade, as well as from lineages of ancestral clades such as 19A, 19B, 20B, and 20C, where the dominant mutation in the Spike protein was sometimes only the D614G mutation (Table 4).

In January 2021, the Alpha variant (B.1.1.7), previously designated as VOC 202012/01, and commonly known as GRY clade (formerly GR/501Y.V1), was identified in Sidi Lahcen, Morocco. Known for its increased transmissibility, Alpha variant led to strengthened health measures in Morocco. The N501Y mutation in the Spike protein, found in the Alpha, Beta, Mu, and Omicron variants (Fig. 4a), is known to enhance virus transmissibility [1724]. In addition, it has been reported that the deletions such as H69del/V70del and Y144del in the NTD region would be involved to immune evasion and increased virulence of SARS-CoV-2 [1824]. The P681H substitution in the S1/S2 cleavage loop, observed in several variants, would favour the processivity of Spike protein cleavage [171925]. Moreover, the coexistence of the N501Y and P681H mutations would increase transmission and disease severity [171925]. However, it is interesting to note that the introduction of the Alpha variant (B.1.1.7) marked a major change in the evolutionary dynamics of SARS-CoV-2 in Morocco. Furthermore, our observations corroborate the findings of a previous study, which reported that at the beginning of 2021, the Alpha variant was dominant in Tunisia, Algeria, and Morocco, while the Eta variant predominated in Egypt, Libya, and Mauritania [26]. In contrast, while Hamzaoui and colleagues observed a frequency of less than 1% for the Alpha variant and an absence of sequences attributed to the Beta variant, our study showed a prevalence of 7.74% for the Alpha variant and 0.31% for the Beta variant [26].

In May 2021, the Beta variant (B.1.351) of SARS-CoV-2 has been detected in Casablanca, Morocco. This variant is characterized by key mutations in the RBD/RBM (N501Y, E484K, K417N) and NTD (L242del/L244del) regions. Previous studies have suggested that these mutations may contribute to SARS-CoV-2′s immune evasion and increased risk of reinfections [192427]. However, despite concerns about B.1.351 in some countries, the impact of B.1.351 in Morocco has been relatively minor [2829]. Morocco later detected the Eta, Kappa, and Mu variants. Eta, formerly known as G/484K.V3, with the first case identified in Morocco in March 2021, carries the E484K mutation known for its role to circumvent host immune responses [171924]. The Mu variant (B.1.621), previous VOI, was identified in Morocco on 22 May 2021. This variant (B.1.621) contains several notable mutations, including T95I, A262S, D614G, D950N, E484K, N501Y, P681H, R346K, and Y144del, which are likely to increase immune evasion and virus transmissibility [3031]. In general, the Beta, Eta, Kappa, and Mu variants did not have a major impact in Morocco, possibly due to effective health measures or competition with more widespread variants such as Alpha, Delta, and Omicron.

In the last quarter of 2020, the epidemiological dynamics of COVID-19 have evolved substantially. India, in particular, has faced an outbreak of several SARS-CoV-2 variants that have spread to affect many countries. Among these variants, B.1.617.1 and B.1.617.2 emerged, both sharing key mutations such as G142D, L452R, and P681R in the spike glycoprotein, along with other distinct mutations likely to increase virus transmission and promote immune evasion [17,19, 3235]. It is important to note that the Kappa and Delta variants have been labelled under the generic name of «double mutants», due to the presence of key mutations L452R and E484Q (Fig. 4a), which had never been observed together before the emergence of these variants [17,19, 3235]. However, the epidemiological situation in Morocco in 2021 was marked by the emergence of B.1.617.2 variant and its subvariants (B.1.617.2/AY.*), which are phylogenetically related to clades 21A, 21I, and 21J. The first strain of the Delta variant (B.1.617.2) was isolated on 22 April 2021, in Casablanca. Following their introduction, they progressively became the predominant variant, supplanting other variants. Between June and September 2021, B.1.617.2 underwent rapid evolution, giving rise to a major evolutionary node with the emergence of subvariants such as AY.33, AY.112, and AY.122, which became dominant nationwide. Compared to the double mutant B.1.617.1, the B.1.617.2 variant and its subvariants harboured mutations in the NTD domain such as T19R, T29A, G142D, F157del, E156G, R158del, T250I, T299I, and in the RBD domain with mutations L452R and T478K. The SD2 region also revealed mutations such as Q613H, D614G, D950N, as well as an amino acid substitution around the furin cleavage site, P681R (mutation shared by Delta and Kappa). All these mutations differentiated them from other variants by their ability to enhance their transmissibility and their capacity to escape the immune system. Although the Delta variant is devoid of the E484Q mutation, analogous to the immune escape mutation « E484K » found in other variants such as B.1.351 and B.1.525, respectively. The T478K mutation hosted in the Spike glycoprotein could accelerate the replication and transmissibility of the pathogen [1736]. Analogously to the D614G mutation, it is hypothesized that the Q613H mutation may enhance the virulence of SARS-CoV-2 by facilitating the cleavage of the spike protein, thereby enhancing fusion between the viral envelope and the cell membrane [37]. Interestingly, among the mutations of the aforementioned Delta subvariants, mutations such as T29A, T250I, T299I, and Q613H were specific to subvariants AY.33, AY.112, and AY.122. These mutations may confer to these subvariants a selective advantage by enhancing their transmission capacity and their ability to escape the immune system [38]. Moreover, previous studies have reported that the F157del and R158del deletions, as well as the L452R and T478K mutations (Fig. 4a), would confer to Delta and its subvariants the ability to escape neutralizing antibodies [1936]. Consequently, all these characteristics make the Delta variant and its subvariants (AY.*) one of the most virulent variants of SARS-CoV-2 [39]. Furthermore, it should be noted that during the emergence of the Delta variant and its subvariants, the AY.33 and B.1.617.2 lineages predominated in Morocco (Fig. 1a and Table 3), while in Egypt and Algeria, the B.1.617.2 variant was dominant. In contrast, the AY.4, AY.122, and AY.34.1 subvariants were the most dominant in Libya, Tunisia, and Mauritania, respectively [26]. These observations highlight the distinct genetic diversity, reflecting the broader global trends in the evolution of the Delta variant and its subvariants [26].

The epidemiological evolution observed in Morocco during the last quarter of 2021 was marked by a transitional period with the emergence of the Omicron variant and its derivatives designated under the GRA nomenclature in GISAID (formerly GR/484A). As of 2022, 82.3% of genomic sequences were assigned to omicron and its subvariants (Table 1). The introduction of the Omicron variant and its subvariants in Morocco resulted in a major epidemiological shift, relegating the previously dominant Delta variant to a marginal position. Among the Omicron strains that have led to substantial infection in Morocco, we note BA.1 and BA.1.1 (displays clade 21K), as well as BA.2 and B.1.1.529 (clade 21L), identified in Morocco in December 2021. However, in June 2022, Morocco rapidly observed a surge in infection cases, with the emergence and concurrent coexistence of subvariants from clade 22B such as BA.5.2.20, BA.5, BA.5.1, BA.5.2.1, and BF.5 (subvariants emerged from Botswana, Hong Kong, and South Africa), succeeding the lineages BA.1 and BA.2, respectively (Fig. 1a and Table 4). Several derivatives of BA.1 subvariant (clade 21K) were also identified, notably the sublineage BA.1.1 (BA.1+R346K), distinguished from BA.1 by the additional presence of mutation R346K. In Morocco, alongside lineage BA.1, lineage BA.2 and its sublineages emerged, identified in February 2022, with a peak reached on April 2022. Lineage BA.1 harbours certain mutations, such as A67V, H69del/V70del, T95I, Y144del, L212I, N211del, ins214EPE, S371L, G446S, G496S, T547K, N856K, L981F, distinguishing it distinctly from lineage BA.2. Conversely, lineage BA.2 is distinguished from BA.1 by specific mutations such as T19I, L24del, P25del/P26del, A27S, V213G, S371F, T376A, D405N, R408S. Our analysis results corroborate the conclusions reported by Balupuri and colleagues in 2023 [40].

Interestingly, variants BA.1, BA.2, BA.4, and BA.5 share several mutation similarities. The mutations L452R and F486V are specific to BA.4 and BA.5, while H69del/V70del are specific to BA.1, BA.4, and BA.5. However, the Q493R mutation is specific to BA.1 and BA.2 [4142]. Moreover, subvariants BA.4 and BA.5 also display mutation similarities in their Spike protein with previous VOCs. These mutations include the N501Y mutation hosted by the Alpha, Beta, and Mu variants; H69del/V70del hosted by the Alpha variant; K417N observed in sequences of Beta variant; E484K present in sequences of the Beta and Eta variants, while L452R and T478K are observed in sequences of the Delta variant (Fig. 4a). Unlike the Beta variant, which features the E484K mutation, subvariants BA.4 and BA.5 harbour the E484A mutation (Fig. 4a). The mutation profiles described in this study align with findings reported in previous studies [27,4145]. Over time, the BA.2 variant has gradually emerged as the predominant subvariant, relegating the BA.1, BA.1.1, and B.1.1.529 variants to a marginal position (Fig. 1a). Furthermore, a notable descendant of the BA.2 lineage, the subvariant BA.2.12.1 (BA.2+L452Q+S704L), is characterized by an additional mutation, L452Q and S704L, absent in BA.2. Nevertheless, most lineages such as BA.2 and its sublineages, XBB and its derivatives, the Beta, Kappa, Mu, B.1, Delta variants and their derivatives were mostly devoid of the H69del/V70del mutation, unlike some Omicron subvariants, the Alpha and Eta variant (Fig. 4a). These observations correlate with previous research findings [414647]. In addition to the previously mentioned subvariants of Omicron, BA.5 and BA.4 are among those being closely monitored. In the same perspective, the BA.5.2.20 subvariant, a progeny of BA.5, has been observed with the highest proportion of infections within the Moroccan population (Fig. 1 and Table 3). Furthermore, the subvariants BA.2.12.1, BA.4, and BA.5 have also shown a higher rate of infection, in part due to mutations such as S371F, D405N, and R408S, which could allow these subvariants to evade antibodies [4348]. Remarkably, in Morocco, the evolution of SARS-CoV-2 has been characterized by substantial genetic diversity, with a particularly marked predominance of the Omicron variant and its subvariants over time. This evolutionary trend has also been observed in several countries during the pandemic [2649]. In Morocco, Omicron and its subvariants accounted for 59.5% of sequences, while Delta was observed at frequency of 11.30 %. Approximate prevalences of Omicron and Delta variants have been reported elsewhere [26]. Hamzaoui and colleagues reported that in Egypt, 55.3% of sequences were assigned to Omicron and 12.1% to Delta. In addition, Omicron subvariants such as BA.5.2.20 and BA.1 were dominant in Morocco (Fig. 1a), accounting for 11 and 10.1% of sequences, respectively (Table 3). Our observations are consistent with those of Sarkar and colleagues, who reported a predominance of BA.1 in Europe in late 2021 and early 2022 [49]. This trend was also observed in the USA, Brazil, and South Africa, in which BA.1 prevailed over BA.2. Furthermore, in June 2022, a trend of BA.4 and BA.5 subvariants was observed. This observation is consistent with previous studies [2649]. However, unlike some North African countries, such as Egypt, the BA.5.2(13.75%) and B.1.617.2(8.74%) variants were the most prevalent, while in Algeria, BA.5.2(17.27%) and B.1.617.2(12%) were predominant [26]. Nevertheless, despite the notable spread of Omicron worldwide, some pre-Omicron lineages persisted and were dominant in other countries [2649]. Notably, in Libya, B.1.525(44.94%) and A (9.55%) predominated; in Tunisia, AY.122(24.23%) and B.1.1.7(21.67%) were dominant; and in Mauritania, B.1.525(31.03%) and AY.34.1(27.58%) were the most dominant [26]. In general, the peaks of predominance and the heterogeneity of SARS-CoV-2 lineages observed in this study and those reported in previous studies highlight the genetic diversity of SARS-CoV-2 during its worldwide spread [2649].

Furthermore, in Morocco, the BQ.1 lineage (BA.5+R346T+N460K) and one of its sublineages BQ.1.1 (BQ.1+R346T), both descendants of BA.5, which emerged from Botswana/Hong Kong/South Africa, have attracted particular interest and were among the monitored subvariants. These lineages (BQ.1 and BQ.1.1), partially aliased to BA.5.3.1.1.1.1.1.1 were identified in September and November 2022, respectively (File S4). These subvariants stood out by several key mutations, notably R346T, K417N, N440K, K444T, L452R, N460K, F486V, E484A, and N501Y, known for their ability to escape antibodies. In particular, the mutations R346T, K444T, and F486V appear to promote evasion of humoral immune response induced by vaccination [48,5052]. Furthermore, the mutations L452R, N460K, and R346T observed in the Spike glycoprotein would be correlated with the virus's ability to bind more efficiently to ACE2 receptors. Consequently, this could increase the virus's transmissibility or its ability to infect cells. These observations have been reported in previous studies [5052]. It is worth noting that in Morocco, BQ.1.1 and BQ.1, rapidly reached their epidemic peak in December 2022 (Fig. 1a). Similar epidemiological trends have been described elsewhere, such as in the United States and Europe [53]. While BQ.1 and BQ.1.1 were declining in Morocco towards the end of 2022, new subvariants emerged simultaneously, among which was the XBB (22F) lineage and its progeny (Fig. 6a, Table 4). Unlike what was observed in Europe and the United States, cases of XBB emerged with rapid increases. In contrast, our study revealed sporadic circulation of XBB and its progeny [53]. In this regard, the ancestor (XBB) of these new subvariants was characterized by mutations such as D614G, D405N, E484A, F486S, F490S, G446S, K417N, N440K, N460K, P681H, R408S, S371F, S477N, T376A, T478K, and V445P, all known for their ability to increase transmissibility and evade the immune system [445455]. Among the descendants of XBB, the GP.2 subvariant of clade 23C (CH.1.1) stood out due to the presence of the key mutation K444T (XBB+K444T), enhancing the virus's evasive capacity. Furthermore, the mutation G446S observed in XBB and BA.2.75.5 lineage related to BA.2.75 conferred them a resistance to antibodies [50]. Additionally, subvariants such as XBB.1.5 (XBB+F486P), XBB.1.16 (XBB+E180V+F486P), and EG.5.1 (XBB.1+Q52 H + F456L + F486 P / EG.5+Q52 H), harboured mutations that conferred them an evasive capacity against neutralizing antibodies, notably R408S, G446, G339D, D405N, F490S, K417N, and Q498R (Table 4) [55]. It should be noted that EG.5.1 (Eris) from clade 23F, descendant of XBB.1.9.2 (clade 23D), has a mutational profile in the Spike glycoprotein similar to that of XBB.1.5 [55]. However, according to updated data from the WHO on 9 February 2024, these previously mentioned lineages are classified as VOIs currently circulating worldwide [56]. Moreover, subvariants such as XBB.1.9.1 (XBB.1+F486P) and XBB.2.3.11 (XBB+T478R + P521S), derived from the XBB.2.3 lineage (23E) (XBB.1+D253G + F486P + P521S), analysed in this study were also less frequent (Table 4). Nevertheless, these lineages are classified by the WHO as VUMs currently in circulation [56]. Furthermore, this study supports previous observations suggesting that Omicron subvariants such as BA.1.1, BQ.1, BQ.1.1, XBB, and its descendants present distinct characteristics and are associated with variable epidemiological trends worldwide [49]. According to our results, these subvariants which appeared between the end of 2022 and during 2023 (XBB and its derivatives) did not lead to large-scale epidemics, as observed elsewhere [4957]. This observation could be attributed to the implementation of health safety measures and protocols, and the mass vaccination campaign in Morocco aimed at achieving herd immunity in order to curb the spread of SARS-CoV-2.

In addition to the key mutations observed in the Spike glycoprotein of different lineages of SARS-CoV-2, some mutations such as P13L and R203K/G204R located in the nucleocapsid protein have been identified in notable proportions among the 2274 genomic sequences analysed (Fig. 4a). Therefore, these mutations are known to be associated with increased infectivity, structural stability, and virulence of the virus, which could confer on it a better ability to escape the immune response [5859]. In general, our analysis revealed a prevalence of 71% of R203K/G204R mutations in our isolates, significantly higher than the proportion of 55% reported by Ahmad and his colleagues [60]. These mutations are observed at high frequency in the omicron and its subvariants, as well as in the Alpha variant and the B.1.1 lineage (Fig. 4a).

However, this study had several limitations. This study included both high and low coverage genomic sequences. Although low-coverage sequences can provide additional information on the genomic diversity of the virus, this can also pose methodological challenges and limitations in the interpretation of results, with possible under-representation of certain viral strains or genome regions. In addition, in-depth documentation of mutations in structural and non-structural proteins could provide a better understanding of the genomic diversity of Moroccan isolates of SARS-CoV-2 and their impact on the dynamics of the virus.

This study provided a detailed analysis of the genomic epidemiology and genetic diversity of SARS-CoV-2 lineages identified in Morocco during the 3 years of the pandemic, enabling a better understanding of the evolution and phylogenetic relationships among different lineages. Several lineages identified in Morocco were closely related to those observed worldwide, except for lineage B.1.528, before their local spread, highlighting the impact of human mobility on the introduction and spread of these lineages during the pandemic. Viral dynamics in Morocco, characterized by a predominance of Alpha, Delta, Omicron variants, and their subvariants, reflected global trends in their evolution. However, the epidemiological trends of some Delta and Omicron subvariants showed variable patterns compared to those observed in other countries. Additionally, several key mutations identified within the lineages analysed were correlated with variations in transmissibility, pathogenicity and antigenicity, which could have contributed to affecting vaccine efficacy and pandemic management. However, the set-up of the SARS-CoV-2 genomic surveillance consortium in Morocco and vaccination campaigns have contributed to control and reduce infection rates and severe forms of COVID-19, thus mitigating the impact of infections at national level.