Evolutionary dynamics of SARS-CoV-2 in space and time during the first phase of the epidemic in Italy

The aim of this study was the reconstruction of SARS-CoV-2 evolutionary dynamics in time and space in Italy and Europe between February and June 2020. The cluster analysis showed that pure Italian clusters were observed mainly after the lockdown and distancing measures were adopted. Lineage B and B.1 spread between late January and early February 2020, from China to Veneto and Lombardy, respectively. Lineage B.1.1 most probably evolved within Italy and spread from central to south Italian regions, and to European countries. The lineage B.1.1.1 entered Italy only in the second half of March and remained localized in Piedmont until June 2020. In conclusion, the reconstructed ancestral scenario suggests a central role of China and Italy in the widespread diffusion of the D614G variant in Europe in the early phase of the pandemic and more dispersed exchanges involving several European countries from the second half of March 2020.


Introduction
SARS-CoV-2 was first described in Wuhan city, China, likely resulting from adaptation of an animal virus to humans and spread rapidly around the world, causing >170 million documented infections and >3.5 million deaths at the end of May 2021 (https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd4029 9423467b48e9ecf6). Italy was the first European country to experience a major SARS-CoV-2 disease (COVID-19) epidemic, with a first wave of transmission characterized by a relatively high number of deaths starting from February 20, in Codogno, Lombardy. A few weeks later, the first lockdown and other containment measures, such as quarantine of travellers returning from high-risk areas, reduced the number of COVID-19 cases in Italy and prevented the escalation of clusters of community transmission.
The COVID-19 pandemic represents an unprecedented challenge for global public health with the continuous emergence of new genetic variants of the virus [1] and the related implications such as their potentially increased pathogenicity or transmissibility and, possibly, vaccine escape. Notwithstanding a unique proofreading activity among RNA viruses [2], SARS-CoV-2 has been exploring its genetic space due to an exceedingly large number of transmissions and replicative cycles, with an estimated evolutionary rate around 2 mutations per month. Indeed, notable genomic variability can be observed among all viral sequences submitted to the GISAID database, which have been grouped into two main lineages, A and B, each containing a growing number of sub-lineages [3]. Both lineages likely separated early during the Wuhan outbreak, with lineage B now being more widely distributed. In this context, the establishment of surveillance networks at national and international level is mandatory to trace the pandemic and inform the appropriate public health interventions.
Presently no comprehensive data are available to establish the lineage of SARS-CoV-2 strains circulating in Italy and their population dynamics, although regional data have been published for Sardinia [4], Lombardy [5] and Abruzzo [6]. The short time since its identification and the limited number of sequences available in public databases makes it difficult to understand the biological significance of the mutations observed so far, whether they are the product of adaptive selection [7] or rather the result of genetic drift due to the high level of genetic variation (http://virological.org/t/response-to-on-the-origin-and-continuing-evolution-of-sars-cov-2 /418). In addition, as Italy can be considered the first and one of the main incubators for the spread of the epidemic in Europe and in the United States, the analysis of SARS-CoV-2 molecular epidemiology since the first phases of the epidemic in this country is of particular interest for unravelling the first evolutionary steps of the virus outside China and its adaptation to western countries. In this context, the reconstruction of the spatial and temporal dynamics is fundamental to understand the origin and evolution of SARS-CoV-2 from the ancestral strains to the new variants. However, there are several limitations in the phylogeographical inference on SARS-CoV-2 genomes, due to the relatively low evolutionary rate of SARS-CoV-2 in comparison with other RNA viruses, the low sampling rate in several countries in the early epidemic and the rapid dissemination of the infection [2]. Furthermore, a frequent homoplasy [8] , [9], that affects multiple protein sites of the viral genome, and the founder effect played a dominant role in the early evolution of the virus [10] , [11]. For these reasons, a limited number of studies on the SARS-CoV-2 phylogeography have been published, based on maximum parsimony [12], maximum likelihood and Bayesian framework [13] , [2] or phylogenetic network [14]. These studies analysed a limited number of genomes available at the time in the short interval elapsed since the origin of the virus [15] while others highlighted the importance of including travel-related information in the analysis [13].
In this study, viral sequences have been analysed for mutations and phylogeny, in comparison to national and international SARS-CoV-2 genomes, to hypothesize the route of arrival to Italy, the subsequent dispersion and further spread to other countries. Major SARS-CoV-2 infection clusters in Italy were identified and characterized and their role in the international virus spread was assessed by using phylogenetic analyses. In addition, the spatiotemporal SARS-CoV-2 dynamics in Italy was investigated by a relatively new maximum likelihood approach for ancestral character reconstruction, by combining the reconstruction relative to the sampling location with the evolutionary lineage.

Specimen Collection
Sequence and epidemiological data were collected at the centres participating to the collaborative group SCIRE (SARS-CoV-2 Italian Research Enterprise), established at the beginning of the epidemic. SARS-CoV-2 RNA positive samples were collected between 24th February to 18th June 2020 from the respiratory tract of individuals who were either hospitalized or tested within screening programs. Samples were collected in most Italian regions, including Apulia, Campania, Emilia Romagna, Lazio, Liguria, Lombardy, Marche, Piedmont, Sardinia, Sicily, Tuscany, Umbria and Veneto.

SARS-CoV-2 Data Sets
In addition to the 62 sequences previously published by this research group, other 192 sequences were obtained in the present study and combined with all the SARS-CoV-2 genomes collected between 29th January 2020 and 18th June 2020 available from Italy in the GISAID database (October 2020, Supplementary Table 1) to form the Italian data set for further analysis (n=465). Fourteen Chinese strains were added obtaining a global dataset of 479 sequences.
To place the Italian sequences in the context of the international COVID-19 pandemic, an additional dataset was built including European and Chinese (n=52) sequences collected in the same period. Due to the large amount of European sequences available, we included at least 2 strains per country/week (n=858). Identical strains or those with more than 5% of gaps were excluded. For countries with a limited number of sequences, all strains were included. Consequently, the final dataset encompassed 1,375 sequences. SARS-CoV-2 sequences were aligned using MAFFT (https://mafft.cbrc.jp/alignment/server/) and the alignment was manually cropped using BioEdit v. 7.2.6.1 (http://www.mbio.ncsu.edu/bioedit/bioedit.html).

Genetic distance
The MEGA X program was used to evaluate the genetic distance between and within Italian strains on full length genome, with variance estimation performed using 1,000 bootstrap replicates. Amino acid changes were evaluated using MN908947 as the reference sequence (https://www.megasoftware.net/).

Phylogenetic Analysis
SARS-CoV-2 sequences were classified using the Pangolin COVID-19 Lineage Assigner tool v. 2.3.2 (last access 15 April 2021, https://pangolin.cog-uk.io/) and Nextclade v. 0.14.1 (https://clades.nextstrain.org/). The maximum likelihood trees of the two data sets were estimated using IQ-TREE v. 1.6.12 (http://www.iqtree.org/), using the GTR+F+R4 (General time reversible + empirical base frequencies + four number of categories) model selected by the program and 1,000 parametric bootstrap replicates for nodes support. Phylogenetic dating was obtained by the least squares dating method (LSD2) implemented in IQ-TREE with 100 replicates to obtain confidence intervals in node ages [17].
The genome sequence hCoV-19/Wuhan/WH04/2020|EPI_ISL_406801|2020-01-05 was used as an outgroup, as it falls in a basal position with respect to the tree and it results within a reasonable estimate of the time of emergence (time to the most recent common ancestor, tMRCA). Italian clusters (including more than 2 sequences) were identified in the ML tree by Cluster Picker v.1.2.3 using 80% bootstrap support and a mean genetic distance of less than 0.3% as thresholds. Epidemiological characteristics of the identified clusters were further investigated using Cluster Matcher v.1.2 [18] which allows the identification of clusters meeting given criteria. The Italian dataset was also analysed by BEAST v. 1.10 (https://beast.community/) in order to estimate the tMRCAs of the main clades. A previously estimated evolutionary rate of 8×10 −4 substitutions/site [19] and the selected substitution model (GTR+I+G) were used as priors.
An exponential coalescent tree, with priors on population size and growth rate and an uncorrelated relaxed molecular clock model with an underlying lognormal distribution were chosen. We ran a Markov chain Monte Carlo analysis of 300,000 steps sampling every 30,000 steps, and we determined sufficient sampling from the posterior by verifying that all parameters had effective samples sizes of at least 200. Finally, all trees were visualised in FigTree v. 1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).
Each taxon was assigned to its sampling locations character, and we used PastML [20] to reconstruct the ancestral character states and their changes along the trees. We used the MPPA as prediction method (standard settings) and added the character predicted by the joint reconstruction even if it was not selected by the Brier score (option-forced_joint). Additionally, we repeated the PastML analysis for the SARS-CoV-2 lineages. Phylogeographycal reconstruction was conducted using both Italian and European datasets.
In order to avoid ambiguities in the root reconstruction and given that lineages B and B.1 represent monophyletic groups we performed two independent ancestral state reconstructions: one for lineages B.1 and its descendent lineages (B.1.1 and B.1.1.1) and one for lineage B. The same outgroup was used for both analysis (EPI_ISL_406800).

Results
A total of 192 SARS-CoV-2-Italian genomes were newly generated for this study. Travel history was available for 137 (71.3 %) patients. All of them reported no international travel in the two weeks preceding the onset of symptoms. One case of contact with a traveller from Bangladesh was reported. Main patients' information is reported in Table  1.

Genetic distances and mutation analyses
The overall mean p-distance between all the Italian isolates was 3.9 (SE: 0.4) s/10,000 nts corresponding to a mean of 10 Seventeen amino acid substitutions were present in more than 10% of the Italian isolates but only one of them was in the spike protein (D614G). No mutations were observed in the RBD domain in the whole Italian sequence dataset. Only eleven B lineage sequences in the whole dataset, all from Veneto (clade 19A), carried T1543I in orf1a. Overall, the B sequences showed a distinct mutations pattern from those of other lineages, including mutations L3606F, G251V in orf1a and orf3a, respectively. The B.1.1.1 lineage presented additional substitutions in comparison with B.1 and B.1.1 lineages such as T1246I in orf1a in all isolates. Table 2 shows the most frequent amino acid substitutions stratified by lineage and clade Table 2. Aminoacid substitutions found in more than 10% of sequences stratified according to lineage and clade. -    The phylogeography of SARS-CoV-2 identified China as the location of the tree root (

Italian clusters
The phylogenetic analysis by ML of the entire dataset including Italian, European and Chinese genomes showed that the majority of the Italian isolates were dispersed in the entire tree. A total of 80 (out of 465, 17.2%) Italian isolates were included in 22 highly supported clusters (Table 4). Of these, 12 (54.5%) were within the lineage B.1, five (22.7%) were B.1.1/20B, three (13.6%) were B.1.1.1/20D and two (9.1%) were B/19A. All but one B.1 clusters were classified as 20A clade. Cluster #19 was the only exception and included four Italian strains classified as clade 20C (all from Rome), showing a mean tMRCA falling in March 2020. Three clusters (13.6%) were singletons (including only single Italian isolates not linked to other Italian sequences), probably corresponding to sporadic introductions followed by limited circulation, while the remaining 19 clusters encompassed at least two Italian isolates, suggesting a local transmission. Thirteen of these (68.4%) included only Italian strains (suggesting a mainly local circulation of this lineage), while 6 (31.6%) included isolates from other European countries, and one of them (B.1) included also one Chinese genome. The estimate of the clusters tMRCA by ML method confirmed that the first transmission events in Italy dated around the second half of January and early February. Eighteen clusters had a common ancestor dating before the introduction of the containment measures in our country. In particular, B.1/20A clusters predominated (10/14) at earlier time points (before March) while in March other clades (20B, 20C and 20D) prevailed (6/8). Moreover, the mixed and singleton clusters were prevalent at the beginning, while pure Italian clusters were the only clusters observed after the lockdown. The earliest cluster (#1), was lineage B.1/20A, dated back to average 20/01/2020 (CI95% 08/01-24/01/2020) and included only four North Italian strains: one from Lodi, two from Milan (the locations where autochthonous COVID-19 cases were firstly identified in Italy) and one from Piacenza. The first B.1.1 cluster dated back to 10/02/2020 (CI95% 28/01/2020-12/03/2020) and included 3 Italian isolates from Abruzzo. Three B.1.1.1/20D clusters dated back to 02 March (CI95% 22/02/2020-02/03/2020). Only two small Italian clusters supported by significant bootstraps were observed within the ML tree including B/19A isolates. In particular a single pure Italian cluster included 11 genomes from Veneto (province of Padua), characterized by the substitution T1543I in orf1a, not detected in any of the other B/19A genomes in our international dataset.

Phylogeographical analysis in Europe
Combining the ancestral state reconstruction for the location with the lineage (Figure  4), the analyses showed that B.1 probably originated in China and spread to several European countries reaching Italy several times, forming a large cluster which included initially 59 (around the first week of March) and finally 198 genomes, and 6 further independent introductions mainly corresponding to a group of genomes characterized only by the substitution D614G but lacking other substitutions, in particular the P314L in the RdRp identifying the clade 20A (lineage B.1, clade 19A).  (Figure 4). A total of 7 nodes remained undetermined. A separate analysis conducted distinguishing European countries (rather than considering a single generalized group), generally confirmed this scenario and allowed to reconstruct in greater detail the dispersion of the epidemic in the European countries (Supplementary Figure 1).
The analysis of lineage B showed that only 2 nodes remained undetermined between Europe and China. The visualization (Figure 5b) suggested several introductions from China to Italy starting with the end of February. A single cluster corresponding to the previously described cluster#5 was observed, while the other strains apparently represent multiple independent introductions forming small groups of no more than 2 sequences. Two sporadic introductions from Europe were also observed. Unlike the ancestral reconstruction for B.1 lineage, this scenario was different as the migratory flows seem to stop in Italy without further spread. The analysis conducted among European countries (Supplementary Figure 2), highlighted the same ancestral scenario but did not show any introduction from Europe.

Discussion
The present study shows that a few different SARS- CoV-2 lineages (B, B.1 and B.1.1, which largely correspond to clades 19A, 20A and 20B) were the most prevalent in Italy since the beginning of the pandemic, accounting for more than 93.8% of the 465 isolates in the dataset considered. This is in accordance with previous studies of the genomic epidemiology of SARS-CoV-2 in Italy performed by our and other Italian research teams [16] , [5]. Nevertheless, we observed important differences in the distribution of such lineages and clades both in space and in time. Several regions, mostly in Northern Italy, showed a high prevalence of B.1 lineage (clade 20A), while other regions, predominantly in Central/Southern Italy, were characterized by a higher prevalence of B.1.1 (clade 20B). Two regions presented a unique scenario with one highly predominant lineage: Veneto, in North-Eastern Italy, where a high prevalence of lineage B (19A clade) was observed in the early hit area of Padua and Piedmont, in North-Western Italy, where lineage B.1.1.1 (clade 20D) was highly prevalent (73%) as opposed to no cases in the other regions.
Lineages B and B.1, as well as clades 19A and 20A, were largely prevalent at the beginning of the epidemic between January and the first half of March, while other lineages (B. 1.1, B.1.1.1) and clades (20B, 20C and 20D) were more prevalent later in the pandemic. The substitution of distinct SARS-CoV-2 lineages over time was confirmed by the estimation of the main lineages tMRCAs, suggesting an early spread in Italy of lineages B and B.1 in the second half of January 2020, followed later by the emergence of other B.1 descendants (B. 1.1 and B.1.1.1) between late February and March. These results correlate with the consolidated epidemiological data showing that the first introduction of imported COVID-19 cases in Italy occurred at the end of January 2020, after the identification of two Chinese tourists in Rome infected by lineage B SARS-CoV-2, apparently without further spread. The first autochthonous cases of COVID-19 in Italy were documented several weeks later (21 February) when the first Italian transmissions without obvious connections with China were described in Lombardy (Codogno and in other centres of the Lodi's area), where lineage B1 was prevalent, and Veneto (the province of Padua), where dominated lineage B at the beginning. On the contrary, the earliest cases in Central-Southern regions were reported a few days later (in Emilia Romagna and Tuscany on February 24 and in Sicily, Abruzzo, Marche and Apulia on February 26), while the number of cases remained relatively low until at least the second half of March (https://www.epicentro.iss.it/coronavirus/sars-cov-2-dashboard).
The analysis by ancestral character reconstruction of the Italian dataset, assigning each taxon to the region where it was sampled and the main viral lineage to which it belongs, showed two distinct patterns of dispersion of SARS-CoV-2. The first pathway concerns lineage B which was introduced to Veneto giving rise to a cluster that apparently disappeared in that region within the first half of March. The second pattern involved lineage B.1 which seems to have entered Lombardy and spread from there to other Italian regions, mainly in Central (i.e. Marche, Abruzzo) and Northern (i.e. Emilia Romagna, Veneto) Italy. This observation is in agreement with the epidemiological data showing the effective suppression of the SARS-CoV-2 outbreak in Veneto in the early times of the epidemic by a highly effective comprehensive testing and tracing approach and local lockdowns [21]. Central Italy (Abruzzo) seems to represent another important centre of dispersion of the lineage B.1.1 (descendant from B.1) mainly to South in mid-March. The introduction of the SARS-CoV-2 lineage B.1.1.1 occurred in Piedmont in the second half of March and apparently did not spread further.
The cluster analysis performed on the ML tree of the global dataset showed few and small (≤11 isolates) Italian clusters, including 17.2% of total Italian strains, while the majority of them were intermixed in the whole tree, frequently near the clades' root (in particular for clades B.1 and B.1.1). This observation is in contrast with data from other European countries (i.e. Spain, Scotland or UK) describing from hundreds to thousands of phylogenetic clusters. This may be related to several reasons, including the poor sampling in the early stages of the Italian outbreak and the low variability of the virus. Nevertheless, Italy was one of the earliest European countries involved in the pandemic, thus the position of Italian strains near to the root of the tree is not surprising and highlights the central role played by Italy in the early spread of the epidemic. Moreover, most of the earliest clusters, dating before the implementation of the Italian national lockdown (2020-03-11), were frequently "mixed" or singletons, including international isolates, while clusters dating after the lockdown were mainly pure Italian. This could be related to the fact that the social distancing measures "froze" the transmission chain and in turn shut off viral evolution. The blockade or limitation of international travelling likely contributed to halt virus spread.
These observations were more deeply investigated by reconstructing the ancestral scenario from the ML trees obtained with the International data set which, similar to that obtained with Italian data set, showed two well defined phylogeographic patterns. The first pathway is that of lineage B showing a large number of independent introductions indicating multiple importation events, most probably from China (or elsewhere in Asia) to Italy as well as to other European countries, even if sporadic introductions to Italy from other European countries could not be excluded. The second phylogeographic scenario involving lineage B.1, showed initially only a few introductions from Asia to Italy and Europe (more precisely defined as Germany in the more detailed country based analysis) of small clusters corresponding to ancestral B.1/19A isolates, characterized by the substitution D614G in the spike protein in the absence of the P314L in the RdRp. A second introduction to Italy corresponded to the largest B.1 cluster characterized by all the substitutions typical of 20A clade giving rise to smaller clusters dispersed to other European and Asian countries and to a further large Italian cluster, corresponding to the lineage B.1.1 which was then dispersed to Italy and to Europe at multiple times. The biggest European cluster descending from B.1.1, corresponding to the lineage B.1.1.1, spread to Italy only in the second half of March. The phylogeographic analysis based on the sampling countries suggested the important roles played from several other European countries, in particular after the second week of March.
The major limitations of this study are intrinsic to the application of the phylogeographical approach to SARS-CoV-2 to a small number of sequences possibly affected by sampling bias for the characters underrepresented in the data set making difficult to reconstruct the directionality of the spread. These limitations could be overcome by adding more sequences and increasing the signal, which would allow to reduce biases and uncertainties. However, classical Bayesian phylogenetic methods, allowing the joint estimation of tree topology with evolutionary parameters and the character state, is very computationally demanding and time consuming; on the other hand, ML, based on ancestral character reconstruction, can be performed in large trees with thousands of tips in a relatively short time [20]. Another limitation of the study is that travel-related information was not available for all cases. Nevertheless, the absence of international travel among those for whom the information was available is probably due to a low and scattered sampling density restricted mainly to symptomatic patients and a prevalent circulation of the virus in small communities rather than in large cities, during the first phase of the epidemic in Italy. It therefore emphasizes the importance of phylogeographic reconstruction in attempting to formulate hypotheses on the possible flows of the virus in the international context.
In conclusion, a possible scenario was reconstructed by employing an ancestral character method allowing the analysis of a large amount of data. Based on our reconstruction, initial multiple sporadic introductions of B lineage to Italy occurred at least since the second half of January 2020 and remained relatively confined. Subsequently, in the month of February the D614G mutant entered in North Italy rapidly spreading to the rest of Italy and Europe, determining a different epidemiological profile of the Italian epidemic since then sustained only by B.1 lineage and his descendants. B.1.1 apparently emerged from the Italian B.1 cluster, suggesting a local evolution, lineage B.1.1.1 most probably emerged from B.1.1 in other European countries, and was introduced in Piedmont after the Italian national lockdown. Overall, our data suggests a central role of Italy in the exporting of some viral lineages at the beginning of the European epidemic, while subsequently, after mid-March, it was an importing centre from other European countries. The introduction in Italy of the D614G variant with a greater transmissibility and its hidden circulation for weeks before the detection of the first cases in Italy could be responsible for the rapid spread of the epidemic in Northern Italy followed by spread to other Italian regions and possibly to the rest of Europe, similar to what was observed for lineage B.1.7.7, firstly predominating in UK and, subsequently, in many other European (and extra-European) countries (eCDC, rapid risk assessment, 15 February 2021).

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Table S1: Accession IDs and sampling dates of sequences included in the dataset, Figure S1: