Genomic Epidemiology of the SARS-CoV-2 Epidemic in Cyprus from November 2020 to October 2021: The Passage of Waves of Alpha and Delta Variants of Concern

The emergence of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in December 2019 resulted in the coronavirus disease 2019 (COVID-19) pandemic, which has had devastating repercussions for public health. Over the course of this pandemic, the virus has continuously been evolving, resulting in new, more infectious variants that have frequently led to surges of new SARS-CoV-2 infections. In the present study, we performed detailed genetic, phylogenetic, phylodynamic and phylogeographic analyses to examine the SARS-CoV-2 epidemic in Cyprus using 2352 SARS-CoV-2 sequences from infected individuals in Cyprus during November 2020 to October 2021. During this period, a total of 61 different lineages and sublineages were identified, with most falling into three groups: B.1.258 & sublineages, Alpha (B.1.1.7 & Q. sublineages), and Delta (B.1.617.2 & AY. sublineages), each encompassing a set of S gene mutations that primarily confer increased transmissibility as well as immune evasion. Specifically, these lineages were coupled with surges of new infections in Cyprus, resulting in the following: the second wave of SARS-CoV-2 infections in Cyprus, comprising B.1.258 & sublineages, during late autumn 2020/beginning of winter 2021; the third wave, comprising Alpha (B.1.1.7 & Q. sublineages), during spring 2021; and the fourth wave, comprising Delta (B.1.617.2 & AY. sublineages) during summer 2021. Additionally, it was identified that these lineages were primarily imported from and exported to the UK, Greece, and Sweden; many other migration links were also identified, including Switzerland, Denmark, Russia, and Germany. Taken together, the results of this study indicate that the SARS-CoV-2 epidemic in Cyprus was characterized by successive introduction of new lineages from a plethora of countries, resulting in the generation of waves of infection. Overall, this study highlights the importance of investigating the spatiotemporal evolution of the SARS-CoV-2 epidemic in the context of Cyprus, as well as the impact of protective measures placed to mitigate transmission of the virus, providing necessary information to safeguard public health.

The sequences used in this study were derived from the molecular epidemiological surveillance study performed by the Laboratory of Biotechnology and Molecular Virology at the University of Cyprus (BMV UCY) in collaboration with the Cyprus Ministry of Health and members of the Cypriot Comprehensive Molecular Epidemiological Study on SARS-CoV-2 (COMESSAR) Network. As part of this collaboration, the BMV UCY retrospectively received 2365 whole-genome SARS-CoV-2 sequences from infected individuals in Cyprus from 04 November 2020 to 08 October 2021 from the Cyprus Ministry of Health.
From the initial 2365 SARS-CoV-2 sequences, 13 were excluded due to missing data (collection dates) and poor sequencing quality (too many missing/ambiguous bases) and were therefore not able to be classified into lineages using the Pangolin webtool (https://pangolin.cog-uk.io/, accessed on 7 May 2022) [30]. The remaining 2352 wholegenome SARS-CoV-2 sequences were included in the analysis. The samples of these 2352 sequences originated from 10 facilities in Cyprus where SARS-CoV-2 testing was performed: 1150 sequences from the Cyprus Institute of Neurology and Genetics (CING); 355 sequences from the Limassol General Hospital; 291 sequences from the Nicosia General Hospital; 188 sequences from the Larnaca General Hospital; 157 sequences from the NIPD Genetics; 120 sequences from the S.C.I.N.A. Bioanalysis Sciomedical Centre Ltd.; 39 sequences from Mygene Molecular Diagnostics Ltd.; 31 sequences from Synlab Cyprus; 14 sequences from Archbishop Makarios III Hospital; and 7 sequences from Tymvios Medical Labs. Sequencing of these samples was performed primarily by the laboratory of Eurofins Genomics Sequencing Europe, which sequenced 1989 of the abovementioned samples; the S.C.I.N.A. Bioanalysis Sciomedical Centre Ltd. sequenced 272 samples, and NIPD Genetics sequenced 91. Bioethical approval (EEBK 21.1.04. 43.01) was provided by the Cyprus National Bioethics Committee for analyses of these sequences, and to ensure anonymity, all sequences were coded with a new laboratory code (double coded) upon being received from the BMV UCY lab. The use of the sequences was in accordance with the relevant guidelines and regulations of the Cyprus National Bioethics Committee.

Sample Collection, RNA Extraction and SARS-CoV-2 Real-Time RT-PCR
Nasopharyngeal and/or oropharyngeal swab samples in transport medium were collected or received by each facility for diagnostic purposes. Specifically, Limassol General Hospital and Archbishop Makarios III Hospital used swabs and transport medium from Jiangsu Yuli Medical Instrument, Shengao, China. Larnaca General Hospital, S.C.I.N.A. Bioanalysis Sciomedical Centre Ltd., and Tymvios Medical Labs used swabs and transport medium from Biocomma Limited, Guangdong, China. NIPD Genetics used swabs from Puritan, Guilford, ME, USA and transport medium from Thermo Scientific, Waltham, MA, USA. Mygene Molecular Diagnostics Ltd. used swabs and transport medium from BIOBASE Shandong, China. Synlab Cyprus used swabs and transport medium from Dakewe BioSci, Sacramento, CA, USA. Nicosia General Hospital, which performed both conventional and fully automated protocols, used swabs and transport medium from the following: 1. Bioprepare, Athens, Greece; 2. KANG JIAN Medical Instrument, Taizhou, China; 3. AMV Healthcare Ltd., Nicosia, Cyprus; 4. Copan, Brescia, Italy (only used for the automated protocol). CING received nasopharyngeal swab samples in transport medium from COVID-19 Public Health Clinics (Nicosia, Limassol, Larnaca, Paphos, Famagusta).  RNA samples from SARS-CoV-2-infected individuals were transcribed into complementary DNA (cDNA) using LunaScript RT SuperMix (NEB) (New England Biolabs, Ipswich, MA, USA). Positive and negative controls were included in each run. The cDNA was then used for enrichment of the viral genome with ARTIC V3 [31] primers (or updated versions thereof) covering the full 29.9-kb viral genome. Amplification was carried out using NEBNext Ultra II Q5 Master Mix (New England Biolabs, Ipswich, MA, USA). The AR-TIC amplicons of SARS-CoV-2-infected samples and corresponding positive and negative controls were converted into sequencing libraries, and sample-specific index codes were attached by PCR prior to pooling. The pool was then quantified and analyzed regarding the size of the pool prior to sequencing. Sequencing was performed using an Illumina NovaSeq6000 (Illumina Inc., San Diego, CA, USA) in 150-bp paired-end read mode. All sequencing was performed using original chemistry provided by Illumina (Illumina Inc., San Diego, CA, USA).
Demultiplexed sequencing data were analyzed using an in-house developed pipeline (Eurofins Genomics Covid Pipeline v.1.0 or updated versions thereof). In brief, the primers used for amplicon generation were scanned and removed from the sequencing reads, followed by removal of low-quality bases (<Q30) and reads shorter than 50 nt in length. The remaining high-quality reads were aligned to the SARS-CoV-2 reference genome (isolated Wuhan-Hu-1, NCBI Reference Sequence: NC_045512.2, GenBank: MN908947.3) using the BWA-MEM algorithm [32]. The read alignments were further refined by retaining only full-length primary alignments, followed by removal of redundant alignments using the Gencore algorithm [33]. High-quality cleaned read alignments were then used to generate a final consensus sequence using the Consensusfixer algorithm [34] with at least 15 supporting reads at each position. Base positions that had <15× or zero coverage were filled with Ns. The NGS methodology used by NIPD Genetics and S.C.I.N.A. Bioanalysis Sciomedical Centre Ltd., employed COVIDSeq Assay (Illumina Inc., San Diego, CA USA) with ARTIC V.3 PCR primers, which were described in detail in our previous publication (Section 2.3. Next Generation Sequencing (NGS)) [29]. However, for S.C.I.N.A. Bioanalysis Sciomedical Centre Ltd., the indexing used was IDT for Illumina-PCR Indexes Sets 1-4 (384 Indexes) (part of the Illumina COVIDSeq RUO); the NextSeqTM 550Dx instrument (Illumina Inc., San Diego, USA) was employed, with samples sequenced at a read length of 2 × 151 bp. Additionally, the methodology used by S.C.I.N.A. Bioanalysis Sciomedical Centre Ltd. for demultiplexing, processing and assembling of the raw sequencing reads to make the final consensus involved DRAGEN COVIDSeq Test (RUO) app version 1.3.0 (Illumina Inc., San Diego, CA, USA).

Mutation Calling
Mutations in the sequences in this study were identified using Nextclade Webtool Web version 1.14.1 (https://clades.nextstrain.org/Nextclade accessed on 7 May 2022) [35]. To avoid misinterpretation of mutations from data that may result from sequencing and assembly artifacts, such as an excessive amount of ambiguous nucleotides, missing data, frameshifts and stop codons, the mutation analysis was based on the 1759 sequences that passed a stricter quality control threshold [ To alleviate computational burden, these lineages were analyzed separately, and their delineation was performed according to the same strategy as in our previous study [29]. Specifically, all newly generated near whole genomes were multiple-aligned using the MAFFT v.7.475 multiple sequence alignment program [36]. The alignments were visually inspected and manually edited with the AliView v.1.26 algorithm [37]. A maximum likelihood (ML) tree was then constructed from the edited alignment using IQtree v.2.1.2 software [38], with branch support estimated through the SH-like approximate likelihood ratio test (SH-aLRT) [39] and ultrafast bootstrap (UFB) procedure [40]. An SH-aLRT and UFB threshold of 90 and 100, respectively, was used to identify the lineages of interest.
To infer time-scaled phylogeographic histories, the genomes from the lineages were complemented with publicly available SARS-CoV-2 genomes. To this end, available complete genomes were downloaded from GISAID on 30 March 2022 (9,792,312 sequences) [41,42]. For each lineage, the 25 genomic sequences most similar to each of the newly generated near-complete genomes were selected using National Center for Biotechnology Information Basic Local Alignment Search Tool (NCBI BLASTn v.2.11.0+) [43], and all duplicate hits were removed. After alignment of these sequences with MAFFT v.7.475 [36] and editing with the AliView v.1.26 algorithm [37], a lineage-specific ML tree was produced using IQtree [38]. This served to further reduce the dataset size by replacing sequences sampled in the same country that cluster together with perfect UFB support by a randomly selected sequence thereof. This procedure yielded datasets with 887 (B. datasets precluded integrated analysis to jointly infer epidemic relationships and migration history. We first generated a time-scaled ML tree using IQtree [38] and LSD2 in R (Rlsd2: R-wrapper for LSD2. R package version 1.10) [44]. For the latter, the rate of evolution was set to 0.0008 s/s/y [45], and taxa with z scores > 3 were considered outliers and removed. To account for the branching uncertainty in phylogeographic reconstructions, polytomies were randomly resolved 1000 times, imposing a small nonzero minimum branch length. The resulting collection of bifurcating trees served as an empirical tree distribution [46] to estimate migration history using an asymmetric discrete phylogeographic model that incorporated a model-averaging procedure (the Bayesian stochastic search variable selection procedure [47]) to identify subsets of migration flows that adequately explain the diffusion process [48], as implemented in BEAST v. 1.10 [49]. Only migration links with Bayes factor support ≥ 5 were reported. The expected number and timing of transitions between locations was estimated using stochastic mapping techniques [50]. For B.1.258 & sublineages, the reported country of sampling was used as the location state. The size and

Time-Scaled Phylogenetic Inference
Time-scaled phylogenies from large clades that represent within-Cyprus circulation were estimated through Markov chain Monte Carlo (MCMC) simulation, as implemented in BEAST 1.10 software [49]. To this end, the substitution process was modeled using an HKY85 model with gamma-distributed among-site rate variation [51,52]. A strict clock was applied with a normal distribution with a mean of 0.0008 s/s/y and a standard deviation of 0.001 as the evolutionary rate parameter [45]. The skygrid model [53] was specified as a flexible tree prior, with change points set at 2-week intervals.

Calculations and Figure Information
The data for Figure 1A,B and Figure 2 were obtained from the Cyprus Ministry of Health, the Press and Information Office and the KIOS Research and Innovation Center of Excellence (KIOS CoE) that operates within the University of Cyprus, and were processed to calculate the percent positivity for Figure 1C [7,54,55]. Specifically, the calculations performed for the percent positivity of Figure 1C were performed by dividing the number of positive SARS-CoV-2 cases per week by the total number of SARS-CoV-2 tests per week and then obtaining the percentage [54]. Regarding the calculations for Figure 2, the number of positive SARS-CoV-2 cases per month reported in Cyprus from March 2020 until October 2021 was indicated proportionally to the most prevalent SARS-CoV-2 variants. Additionally, for Figure 2, the illustration of the spike proteins on the colored virions underneath the names of each lineage, as also seen in figures below, was produced by PyMol (Version 2.4.1, Schrödinger, LLC, https://www.pymol.org, accessed on 18 February 2021) and is based on data derived and adapted from Protein Data Bank entry 6XEY [56,57], as well as other sources used to outline spike protein domains [58][59][60][61][62][63][64][65][66].

The Appearance of Lineages and the Waves of SARS-CoV-2 Infection in Cyprus
For this study, 2352 SARS-CoV-2 whole-genome sequences obtained in Cyprus were analyzed for the period from November 2020 to October 2021. The lineage classification analysis revealed a plethora of lineages, with 61 different lineages being identified over the course of this study period (Table 1). For the first three-month period (November 2020-January 2021) (    Figure 1E,F and Figure 2). The Alpha variant (B.1.1.7) (87.69%, 627/715) and Q. sublineages (0.28%, 2/715) became dominant, reaching a prevalence of 87.97% (629/715) ( Table 1). During this period, there was also an increase in percent positivity ( Figure 1B,C), though it was less pronounced with respect to the early phases of the first and second waves of infections. The demise of the third wave coincides with the declining prevalence of the Alpha variant during June 2021, and it was last identified in July 2021 ( Figures 1A and 2). Around this period, the first Delta (B.1.617.2) sequences were detected at low prevalence during the February 2021 to April 2021 period, with 2/715 (0.28%) sequences in April 2021 (Table 1 and Figure 1E Table 1). This period marked the appearance of the fourth wave in Cyprus, which was driven almost uniquely by the Delta variant ( Figure 2). Once more, the upsurge of SARS-CoV-2 infections was accompanied by an increase in percent positivity ( Figure 1B,C). At the start of the fourth wave, some non-Delta variants were detected, including one sequence (0.16%; 1/644) of the Beta variant (B.1.351). Delta and its AY. sublineages gained total dominance in the following months (Table 1 and Figure 1).
In summary, by the end of the sampling period, there were three waves of infections, each characterized by particular dominant groups of lineages. These lineages were (in order of chronological appearance): B.

Spike Protein Mutations of the Most Prevalent Lineages/Variants in Cyprus
The most prevalent lineages that were identified in Cyprus during the period of November 2020 to October 2021, as mentioned above (Section 3.1), were B.  (Tables S1 and S2). The analyses in this section focused on the S protein due to its important role in the transmissibility of the virus, as well as the fact that it is a target for vaccine design and diagnostics [67]. To identify the characteristic and most common mutations of these lineages, sequences were input into Nextclade Webtool [35], and the output was sorted to isolate S protein mutations (deletions and substitutions).
The most common mutations of each lineage are depicted in Figure 3 and tabulated in Table S1. First, during the first wave in Cyprus (Figures 1 and 2), which was described in our previous study [29], B.1.1 (specifically the B.1.1.29 lineage) was the most prevalent, and the sole most common mutation found in the S protein was D614G; thus, it is not included in Figure 3. The most common mutations for the most prevalent lineage in the second wave in Cyprus, B.1.258 and its sublineages, were ∆H69/V70, N439K and D614G. However, it is important to note that the two sublineages of B.  (Table S1).
Last, for the Delta (B.1.617.2) variant and its AY. sublineages, which dominated the fourth wave, the most common mutations were T19R, G142D, ∆E156/F157, R158G, L452R, T478K, D614G, P681R, and D950N ( Figure 3). Note that the mutations ∆E156/F157 and R158G can also be reported as E156G and ∆F157/R158 [68]. Unlike B.1.258 and Alpha (B.1.1.7), for which only a few sublineages were detected, 30 different sublineages of the Delta variant were identified in this dataset (Table S1). Contrary to B.1.258 and Alpha (B.1.1.7) that substantially outnumbered their sublineages, the parental lineage for Delta, B.1.617.2, was not the predominant lineage of this family (Table S1). In fact, the predominant sublineage of Delta was AY.122, at 52.91% (364/688), and that of parental B.1.617.2 was only 5.96% (41/688) of the Delta lineages in this dataset. As a result, analysis of the most common mutations is biased toward better represented sublineages. For example, the T95I S protein mutation is reported to be common in 13 of the 30 Delta AY. sublineages identified in this dataset according to the CoV-Spectrum website (https://cov-spectrum.org/, date accessed 30 May 2022) [23]; in other Delta AY. sublineages, it may be found in lower prevalence (Tables S1 and S2). Similarly, the A222V mutation is commonly found in AY. lineages such as AY.4.2 (Delta Plus) and AY.60; however, as they are represented in lower numbers or not the predominant Delta lineages (Table 1 and Table S1), their mutations were not as commonly represented in the whole Delta family in this dataset (Table S1) [69].
The above shows that as the pandemic progressed, SARS-CoV-2 variants incorporated an increasingly diverse set of mutations ( Figure 4). Lineages that evolved during the early stages of the pandemic (B.1 and its direct descendants) mostly carried the D614G mutation of the S protein [70]. As B.1 became the dominant lineage around the world, essentially all its descendant lineages, including those described in this section, retained the D614G mutation [70]. Other recurrent mutations that appeared later, such as the ∆H69/V70 deletions [71], were not found in all lineages; in this dataset, they are common for B.1.258 and the Alpha variant. However, the ∆H69/V70 deletions do not appear to have evolved from the same ancestor for B.1.258 and the Alpha variant ( Figure 4 and [72]). Indeed, most mutations occurred in the S protein for the Alpha variant, and the entirety of mutations for B.1.258 are concentrated in the S1 subunit, specifically at the start of the N-terminal domain (NTD) and in the receptor-binding domain (RBD) (Figure 3). Similar to the Alpha variant and B.1.258, the Delta variant, which carries approximately the same number of mutations within the S protein as the Alpha variant in this dataset, shows the majority of those mutations concentrated in the S1 subunit (NTD and RBD) (Figures 3 and 4). Furthermore, there are more mutations in the RBD, specifically in the receptor-binding motif (RBM), for Delta (L452R and T478K) than for Alpha (N501Y) or B.1.258 (N439K) (Figure 3). Interestingly, both the Alpha and Delta variant contain a mutation at amino acid P681 in the furin cleavage site that separates the S1 and S2 subunits ( Figure 3). Finally, between B.1.258, the Alpha variant and the Delta variant, only the latter two harbor mutations found in the S2 subunit, namely T716I, S982A and D1118H for Alpha and D950N for Delta ( Figure 3). of mutations for B.1.258 are concentrated in the S1 subunit, specifically at the start of the N-terminal domain (NTD) and in the receptor-binding domain (RBD) (Figure 3). Similar to the Alpha variant and B.1.258, the Delta variant, which carries approximately the same number of mutations within the S protein as the Alpha variant in this dataset, shows the majority of those mutations concentrated in the S1 subunit (NTD and RBD) (Figure 3 and Figure 4). Furthermore, there are more mutations in the RBD, specifically in the receptorbinding motif (RBM), for Delta (L452R and T478K) than for Alpha (N501Y) or B.1.258 (N439K) (Figure 3). Interestingly, both the Alpha and Delta variant contain a mutation at amino acid P681 in the furin cleavage site that separates the S1 and S2 subunits ( Figure 3). Finally, between B.1.258, the Alpha variant and the Delta variant, only the latter two harbor mutations found in the S2 subunit, namely T716I, S982A and D1118H for Alpha and D950N for Delta ( Figure 3).

Timed Migration Histories
The history of spread was reconstructed for lineages that dominated the second to fourth infection waves in Cyprus (Figures 6-9, Table 2). The locations and estimated numbers of import and export events for B.  Figures 7-9, respectively, as derived from data in Table 2.
For B.1.258 and its sublineages, the first import event to Cyprus was previously dated to 07 March 2020 (95%HPD: 17 January 2020-16 April 2020) [29]. However, no B.1.258 & sublineages were detected from June to August 2020. Hence, the likely date of the first introduction of this lineage's variants that spread in Cyprus after the onset of the second wave in September 2020 was also estimated, which indicated that the B.1.258 & sublineages responsible for the second wave first started spreading in Cyprus on approximately 23 August 2020 (95%HPD: 6 August 2020-31 August 2020). In line with previous findings, the majority of B.1.258 and sublineage taxa cluster within a large and nearly perfectly supported clade consisting almost exclusively of genomes sampled from Cyprus [29] ( Figure  6). The demographic history of this clade indicated a period of exponential growth until November 2020, after which a rather stable plateau was reached that persisted through the most recent sampling date [29]. The additional data from this study allowed for refinement of the aforementioned analysis, indicating small growth and decline periods following November 2020 until the final decline of B.1.258 lineages by February-March 2021 ( Figure S1). The fact that this clade encompasses most sampled B.1.258 infections indicates

Timed Migration Histories
The history of spread was reconstructed for lineages that dominated the second to fourth infection waves in Cyprus (Figures 6-9, Table 2). The locations and estimated numbers of import and export events for B.  Figures 7-9, respectively, as derived from data in Table 2.
For B.1.258 and its sublineages, the first import event to Cyprus was previously dated to 07 March 2020 (95%HPD: 17 January 2020-16 April 2020) [29]. However, no B.1.258 & sublineages were detected from June to August 2020. Hence, the likely date of the first introduction of this lineage's variants that spread in Cyprus after the onset of the second wave in September 2020 was also estimated, which indicated that the B.1.258 & sublineages responsible for the second wave first started spreading in Cyprus on approximately 23 August 2020 (95%HPD: 6 August 2020-31 August 2020). In line with previous findings, the majority of B.1.258 and sublineage taxa cluster within a large and nearly perfectly supported clade consisting almost exclusively of genomes sampled from Cyprus [29] ( Figure 6). The demographic history of this clade indicated a period of exponential growth until November 2020, after which a rather stable plateau was reached that persisted through the most recent sampling date [29]. The additional data from this study allowed for refinement of the aforementioned analysis, indicating small growth and decline periods following November 2020 until the final decline of B.  Figure S1). The fact that this clade encompasses most sampled B.1.258 infections indicates that the second wave was to a large extent driven by local transmission and not by frequent importation of infections. This is reflected by the low numbers of import and export events (Table 2), with a combined maximum of less than eight per week ( Figure 10). As the estimated numbers of introductions are sensitive to the sample size, these represent lower boundary estimates. In our reconstructions, we detected no more than two import events per week, except for the second week of January. Conversely, the weekly number of export events was usually higher than the number of weekly importations ( Figure 10). This disparity between the weekly number of import and export events aligns with their imbalance in their overall number of import/export events: for every importation, there were 2.4 exportation events ( Table 2). Although Slovenia was previously found to be a well-supported origin for import into Cyprus, this was no longer the case (i.e., Bayes factor support for a link between Slovenia and Cyprus now at >5). Instead, Sweden and the UK were identified as the only well-supported origins of B.1.258 in Cyprus, each accounting for approximately half of the import events ( Figure 7, Table 2). Overall, this updated analysis of the migration patterns of B.1.258 lineages in Cyprus substantiates and refines our previous results: the UK was previously identified as a major origin location, but the additional data reveal that Sweden has also played an important role. Similarly, the UK, Czech Republic and Denmark were previously found to be well-supported destination locations for migrations out of Cyprus, and this list is now appended with Sweden and Greece ( Figure 7, Table 2).
The first import event of an Alpha lineage (B.1.1.7 & Q. sublineages) was estimated to be approximately 25 November 2020 (95%HPD: 12 November 2020-08 December 2020). This date of import differs only by 6 days from our earlier estimate, which was 01 December 2020 (95%HPD: 08 November 2020-19 December 2020) [29]. Like for B.1.258, there was also a large clade with predominantly Cypriot Alpha taxa (83.5%, 457/547) ( Figure 6). Despite the low support for this clade in the ML tree (SH-aLRT and UFB support values of 73.7 and 12, respectively), its presence by itself is indicative of substantial within-Cyprus spread. This is corroborated by the presence of several large well-supported subclades, as indicated by a posterior support ≥ 0.8 in analysis of this clade in BEAST v.1.10 [49] (Figure S2). In addition to this large clade, two other smaller clades with almost exclusively Cypriot taxa were noted (n = 61/72 and 78/86 taxa) and SH-aLRT/UFB support of 83.8/100 and 73.1/21, respectively. Thus, a total of 596 Cypriot Alpha variants cluster in these clades, which supports the view of substantial within-Cyprus spread of Alpha-lineage infections. The inferred total number of import/export events, however, is higher than for B.1.258 (Table 2, Figure 10), which led to higher levels of weekly import/export events compared to B.1.258 ( Figure 10). Specifically, the sum of the number of import and export events was mostly higher than 5 per week, even reaching 18 per week in the beginning of May 2021 ( Figure 10). The average of the estimated total number of import events was 94.09 (100%), with Greece and the UK being the primary sources of import, accounting for 42.85 (45.54%) and 29.27 (31.11%) of import events, respectively. Sweden completes the list of countries that account for >5% of the import events into Cyprus for this lineage (Figure 8, Table 2). The observed differences between the current and previous results are most likely due to the much wider availability of infection samples due to Alpha, both in Cyprus (10 vs. 900 genomes for both studies) and elsewhere (the GISAID database [41] Figure S3). The absence of taxa from elsewhere that are closely related to the relevant Cypriot taxon indicates that a long branch links the relevant Cypriot taxon to its first common ancestor with taxa from elsewhere and that the branches in the relevant clade are rather long. Because of this, the uncertainty in the location state reconstruction propagates deep in time. Therefore, the date of first introduction was estimated from a subclade with the earliest sampled Delta genomes in Cyprus, while excluding the problematic taxa shown in Figure S3 Table 2). The other source of imports that accounted for >5% of the total number was Southern Asia (namely, Pakistan and Bangladesh), with 38.20 (7.33%) events ( Figure 9, Table 2). Most of these locations were also identified as destinations of exports from Cyprus. Specifically, of an average total of 576.49 (100%) export events, the majority were to Denmark at 141.53 (24.55%). Other notable destinations of export (above 5% of exports) were Sweden at 91.58 (15.89%), the UK at 87.96 (15.26%), Germany at 54.36 (9.43%), Switzerland at 48.35 (8.39%), and Greece at 40.71 (7.06%) (Figure 9, Table 2). Furthermore, the four most important sources of Delta (B.1.617.2 & AY. sublineages) imports into Cyprus accounted for almost equal shares (~16-18% of all imports), which shows that this variant did not predominantly enter Cyprus through one location. Such similarity was not seen for exports.       Countries acting as "sources" or "sinks" for SARS-CoV-2 Alpha (B.1.1.7 & Q. sublineages) transmission are highlighted and labeled, and the estimated average number of migration events is indicated.
To enhance the clarity of the figure, only countries that account for 5% and above of the estimated total average number of import or export events are displayed, except if there were fewer than five countries (       lineages that were denoted as variants of concern (VOC) [16,17]. b "From" indicates the migration events of a Countries/subregions are as denoted by (UN) geographical subregion. d-f Represent average Markov jumps based on the lower and upper bounds of the of the 95% HPD interval migration events towards and from Cyprus. Only links supported by a base factor of at least 5 were considered. g "All" Represents the aggregation of the migration events from each country/subregion. Figure 10. Temporal dynamics of SARS-CoV-2 export and import from and to Cyprus. The width of each column corresponds to one week, and the height corresponds to the mean of the estimated total number of migration events (to/from Cyprus) inferred to have occurred in that week. Exports are shown in gray columns; imports are shown in black columns. The y-axis represents the number of import/export events per week; the x-axis represents time.

Discussion
In the current study, 2352 whole-genome SARS-CoV-2 sequences were analyzed and 61 SARS-CoV-2 lineages were detected in Cyprus during the period of November 2020 to October 2021. The majority of those lineages belong to three groups: B. The period of focus of the present study was November 2020 to October 2021, partly overlapping and largely extending the period of genomic surveillance of our previous work-April 2020 to January 2021 [29]; hence, there was an overlapping period from November 2020 to January 2021. The uninterrupted surveillance allowed for studying the continuous development of the SARS-CoV-2 epidemic in Cyprus, which was initially characterized by an upsurge of B.1.1 lineages (namely, B.1.1.29) (Figures 1 and 2). These lineages harbor the S gene D614G mutation, which results in increased viral load and infectivity, thereby conferring a fitness advantage over lineages without this mutation [73]. Interventions including lockdown and travel restrictions contributed to mitigating the spread of SARS-CoV-2, preventing the collapse of the health and economic systems until they had adjusted to the situation [29,55,74] (Figure 1). The first full lockdown (24 March-03 May 2020), which mainly entailed a curfew and the restriction of movement without permission, was followed by a partial lockdown (04-20 May 2020). The latter mainly entailed a curfew, which was followed by the gradual lifting of travel restrictions (start of travel restrictions 04 April 2020, repatriation program 14 April 2020, gradual lift STAGE-A 09-19 June 2020 and STAGE-B 20 June 2020), with countries being categorized in accordance with their epidemiological profile at the time [74].
After restrictions were alleviated over the summer period, the SARS-CoV-2 burden began to increase again by the start of September 2020 [74] (Figure 1), which coincided with the emergence of B.1.258 and its sublineages (Figures 1 and 2). The first import of Figure 10. Temporal dynamics of SARS-CoV-2 export and import from and to Cyprus. The width of each column corresponds to one week, and the height corresponds to the mean of the estimated total number of migration events (to/from Cyprus) inferred to have occurred in that week. Exports are shown in gray columns; imports are shown in black columns. The y-axis represents the number of import/export events per week; the x-axis represents time.

Discussion
In the current study, 2352 whole-genome SARS-CoV-2 sequences were analyzed and 61 SARS-CoV-2 lineages were detected in Cyprus during the period of November 2020 to October 2021. The majority of those lineages belong to three groups: B. The period of focus of the present study was November 2020 to October 2021, partly overlapping and largely extending the period of genomic surveillance of our previous work-April 2020 to January 2021 [29]; hence, there was an overlapping period from November 2020 to January 2021. The uninterrupted surveillance allowed for studying the continuous development of the SARS-CoV-2 epidemic in Cyprus, which was initially characterized by an upsurge of B.1.1 lineages (namely, B.1.1.29) (Figures 1 and 2). These lineages harbor the S gene D614G mutation, which results in increased viral load and infectivity, thereby conferring a fitness advantage over lineages without this mutation [73]. Interventions including lockdown and travel restrictions contributed to mitigating the spread of SARS-CoV-2, preventing the collapse of the health and economic systems until they had adjusted to the situation [29,55,74] (Figure 1). The first full lockdown (24 March-03 May 2020), which mainly entailed a curfew and the restriction of movement without permission, was followed by a partial lockdown (04-20 May 2020). The latter mainly entailed a curfew, which was followed by the gradual lifting of travel restrictions (start of travel restrictions 04 April 2020, repatriation program 14 April 2020, gradual lift STAGE-A 09-19 June 2020 and STAGE-B 20 June 2020), with countries being categorized in accordance with their epidemiological profile at the time [74].
After restrictions were alleviated over the summer period, the SARS-CoV-2 burden began to increase again by the start of September 2020 [74] (Figure 1), which coincided with the emergence of B.1.258 and its sublineages (Figures 1 and 2). The first import of B.1.258 was estimated to be 07 March 2020 (95%HPD: 17 January 2020-16 April 2020) and traces back to the detection of two B.1.258 infections as early as April-May 2020. No other B.1.258 sequences were identified during June-August 2020, which we believe is linked to the abovementioned measures, lockdowns, enhanced contact tracing, and travel restrictions that helped to mitigate the import and spread of infections (explained in detail in [29]). After these restrictions were alleviated, B.1.258 was estimated to have been reintroduced on approximately 23 August 2020 (95%HPD: 06-31 August 2020) and was first detected approximately a month later. Its prevalence continued to rise until December 2020. The B.1.258 prevalence began to decline in January 2021, and by February 2021, the second wave had ended. It is of note that the second partial lockdown was implemented 23 October 2020-9 January 2021 and was supplemented with governmental decrees that were issued depending on the epidemiological profile of the island, such as the tighter curfew during 13 November 2020-9 January 2021 for two cities in Cyprus (Limassol and Paphos) [74].
Although these interventions clearly could not completely abrogate the spread of SARS-CoV-2 infections, they are likely to have prevented an even higher burden. In turn, a more beneficial situation at the start of the subsequent complete lockdown means that the latter could have led to a quicker reduction in the SARS-CoV-2 burden. The majority of imports and exports involved Sweden and the UK, where B.1.258 was highly prevalent during this period [72,75]. The travel status of both countries was reclassified from category C (increased risk, requiring negative laboratory test results within 72 h prior to traveling, upon arrival, and 14-day isolation with testing at the end of the isolation period) to B (possibly low risk, requiring negative laboratory test results within 72 h prior to traveling from that country) in August 2020 [74]. As persons with a negative test have a small but nonzero probability of carrying SARS-CoV-2 infection [76], an imported lineage may eventually become prevalent.
The successful global spread of B.1.258 lineages is likely linked to acquisition of beneficial mutations. Two notable mutations of B.1.258 are the ∆H69/V70 deletions (Figures 3 and 4). Interestingly, these deletions have recurrently emerged in other lineages, including Alpha (B.1.1.7 & Q. sublineages) [17,19]) as well as B.1.375 and B.1.1.298 [71,72]. ∆H69/V70 are located in the NTD of the S protein, a region with epitopes for neutralizing antibodies and mutational hotspots, and was even proposed to perhaps play a role in virus entry [77][78][79]. The deletions ∆H69/V70 are reported to enhance infectivity through increased spike incorporation into virions, to increase replication, and to even compensate for immune escape mutations that reduce infectivity [71,80]. Additionally, the ∆H69/V70 S gene mutations have been linked to the failure of tests used for detection of SARS-CoV-2 (S gene target failure) [81]. A second mutation, which also independently emerged in other lineages, is N439K. This site is located in the RBD of the S protein, a region of particular importance due to its involvement in binding to the host cell receptor (ACE2), and it serves a prominent target of neutralizing antibodies [82]. The N439K mutation is associated with immune evasion as well as increased affinity for ACE2 binding [80,83]. The fact that such mutations have also been detected in other lineages, some of which have become dominant, highlights the value of complementing investigations of the functional attributes of mutations with studies of their impact on epidemic spread.
The first Alpha variant infection was detected in December 2020 (Figure 1), soon after its estimated first introduction date (25 November 2020, 95%HPD: 12 November 2020-08 December 2020). This lineage was the first to be denoted as a VOC by the WHO [16] and carried a multitude of S gene mutations. Overall, the mutations characterizing the Alpha variant (B.1.1.7 & Q. sublineages) are linked to immune evasion and increased transmission/infectivity. In addition to the ∆H69/V70 deletion it carries in common with B.1.258, other mutations include the following: ∆Y144, N501Y, A570D, D614G, P681H, T716I, S982A, and D1118H (Figures 3 and 4). ∆Y144 is located in the NTD of the S pro-tein and was found to reduce antibody binding affinity, thereby conferring additional immune evasion potential [84]. The N501Y RBD mutation, which has also been identified in three other VOCs (Beta, Gamma and Omicron), is reported to increase ACE2 binding to the S protein RBD, enhancing infectivity [85,86]. The A570D S protein SD1 domain mutation is mainly implicated in increased infectivity [87] and may contribute to immune evasion [80,88]. On the other hand, the mutation P681H, which is adjacent to the furin cleavage site (S1/S2), has no significant impact on viral entry or the spread of the virus between cells, even though it may increase spike cleavage by furin-like proteases [88,89]. The S protein S2 subunit mutation T716I may cause a destabilization effect to favor RBD-up or open states, which is compensated for by substitutions such as D1118H, which stabilize the prefusion spike conformation [80,90]. Finally, S982A is located in the S protein S2 subunit and has been found to enhance viral entry while reducing neutralizing antibody induction [80].
The fact that the Alpha variant quickly became globally dominant after its emergence in the UK in September 2020 [91] is a testament to its competitive advantages with respect to other cocirculating lineages at the time. The highly transmissible nature of the Alpha variant is also evidenced by the fact that a complete lockdown (from 10 January 2021 to 09 May 2021, Figure 1) could not entirely control its spread, as it was able to during the first and second waves. Moreover, Greece (comprising the majority of arrivals from travelers [74,92]), the UK, and Sweden were classified as category B or C travel locations throughout the period of the partial and complete lockdowns, and phylogeographic reconstructions indicate that they were the source of 88% of the imported Alpha lineages in Cyprus ( Figure 8 and Table 2) [74]. Nonetheless, the travel restriction category assignment of countries was highly dynamic [74]. The demise of the Alpha lineage in Cyprus was likely facilitated by increasing SARS-CoV-2 vaccination coverage. Indeed, vaccination began in Cyprus in December 2020/early January 2021, with priority depending on age/risk group, and by 31 May 2021, vaccination coverage had reached approximately 55% for the first dose and 35% for full vaccination in the total adult population [7,74].
Soon after the peak prevalence of Alpha variants (B.  [93,94], became dominant, resulting in the fourth infection wave (Figures 1 and 2). The Delta variant was estimated to have been imported into Cyprus in late March (24 March 2021, 95%HPD: 1 November 2020-19 April 2021), with the first sequences being identified during April 2021 ( Figure 1). Hence, only a few weeks separated the first introduction of B.1.258, Alpha and Delta and their first detection by the genomic surveillance program. In contrast to the decreasing Alpha variant prevalence near the end of the last complete lockdown, the Delta variants increased in relative importance. The number of different well-supported migration links for import to and export from Cyprus is notably higher for Delta than for B.1.258 and Alpha. We believe that to a large extent, this is linked to the combination of Delta's increased infectivity, gradually allowing unrestricted or less restricted travel from more countries, and changes in travel intensity. At the beginning of May (10 May 2021), there were more than 40 countries in the red category, which decreased to 20 by the end of June (28 June 2021), at which point there were also 28 countries in the green category [74]. This was accompanied by an increase in the number of travelers of approximately 40,000 in April 2021 to 100,000 in May 2021 and more than 180,000 in June 2021 [92]. The travel intensity data also coincide with our results showing that Delta imports in Cyprus remained low until June and July 2021, after which there was a gradual rise in import events of this variant ( Figure 10).
The fourth wave occurred during summer of 2021. Thus, climate alone, albeit impactful, cannot explain the timing of increased infection rates [95] and underscores the importance of viral genetic factors. The most commonly found Delta mutations seen in our sequences were T19R, G142D, ∆E156/F157, R158G, L452R, T478K, D614G, P681R, and D950N (Figures 3 and 4). T19R, G142D, ∆E156/F157 and R158G are located within the antigenic supersite of the S protein NTD and affect immune evasion [96,97]. Furthermore, ∆E156/F157 and R158G contribute to increased infectivity [98]. The L452R and T478K S protein RBD mutations are reported to improve transmission and confer immune evasion [73,99,100]. L452R is also reported to reduce the effectiveness of antibodies induced by vaccination or natural infection [101]. P681R, adjacent to the furin cleavage site (S1/S2), improves transmission, confers partial resistance to neutralizing antibodies and enhances pathogenicity [99,102,103]. No obvious impact was reported for the S protein S2 D950N mutation, though it is hypothesized to enhance the fusogenicity of Delta spike [104].
Many of the abovementioned mutations, such as (but not limited to) ∆H69/V70, L452R, T478K, N501Y, and P681H [108], have reoccurred in a number of lineages, and their combination may confer an even more increased transmission potential to an emerging variant. Such was the case for the Omicron VOC, which incorporates over 30 S gene mutations and became the most dominant variant around the world, including Cyprus (Manuscript in preparation), soon after its detection during November 2021 [109]. As such, the Alpha, Delta and Omicron lineages are real-world examples of the principle of survival of the fittest [110].
In conclusion, this study details the ever-evolving and dynamic nature of the SARS-CoV-2 epidemic in Cyprus from November 2020 to October 2021. Most detected lineages belong to three groups that had acquired a set of S gene mutations conferring increased potential for transmission and immune evasion and are related to subsequent periods of increased SARS-CoV-2 spread in Cyprus: B.    Onderzoek- Vlaanderen', G066215N, G0D5117N and G0B9317N).

Institutional Review Board Statement:
The SARS-CoV-2 sequences were analyzed after receiving approval by the Cyprus National Bioethics Committee (EEBK 21.1.04.43.01). The sequences were coded with a laboratory or patient identification number, and then were further coded with a new laboratory identification number to ensure no connection of the sequences could be made to the corresponding study subjects. The collection and use of the samples were in accordance with the relevant guidelines and regulations of the Cyprus National Bioethics Committee.
Informed Consent Statement: Patient consent was waived by the Cyprus National Bioethics Committee (EEBK 21.1.04.43.01) due to the anonymity and double-coding of the sequences.

Data Availability Statement:
The sequences will be available upon publication of the manuscript on the GISAID database [41,42]. Specifically, only the 1759 sequences that were denoted as good quality under "qc.overallStatus" are available to GISAID, to avoid putative misinterpretations of mutations from data that may have derived during sequencing and assembly [35].