key: cord-0720626-k0gehabh
authors: Lai, Alessia; Bergna, Annalisa; Toppo, Stefano; Morganti, Marina; Menzo, Stefano; Ghisetti, Valeria; Bruzzone, Bianca; Codeluppi, Mauro; Fiore, Vito; Rullo, Emmanuele Venanzi; Antonelli, Guido; Sarmati, Loredana; Brindicci, Gaetano; Callegaro, Annapaola; Sagnelli, Caterina; Francisci, Daniela; Vicenti, Ilaria; Miola, Arianna; Tonon, Giovanni; Cirillo, Daniela; Menozzi, Ilaria; Caucci, Sara; Cerutti, Francesco; Orsi, Andrea; Schiavo, Roberta; Babudieri, Sergio; Nunnari, Giuseppe; Mastroianni, Claudio M.; Andreoni, Massimo; Monno, Laura; Guarneri, Davide; Coppola, Nicola; Crisanti, Andrea; Galli, Massimo; Zehender, Gianguglielmo
title: Phylogeography and genomic epidemiology of SARS-CoV-2 in Italy and Europe with newly characterized Italian genomes between February-June 2020
date: 2022-04-06
journal: Sci Rep
DOI: 10.1038/s41598-022-09738-0
sha: 03e4d42553b61b65dfbc7f87d5f85c0de6bb2e9f
doc_id: 720626
cord_uid: k0gehabh

The aims of this study were to characterize new SARS-CoV-2 genomes sampled all over Italy and to reconstruct the origin and the evolutionary dynamics in Italy and Europe between February and June 2020. The cluster analysis showed only small clusters including < 80 Italian isolates, while most of the Italian strains were intermixed in the whole tree. 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 (20B) most probably evolved within Italy and spread from central to south Italian regions, and to European countries. The lineage B.1.1.1 (20D) developed most probably in other European countries entering Italy only in the second half of March and remained localized in Piedmont until June 2020. In conclusion, within the limitations of phylogeographical reconstruction, the estimated ancestral scenario suggests an important 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.

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:// gisan ddata. maps. arcgis. com/ apps/ opsda shboa rd/ index. html#/ bda75 94740 fd402 99423 467b4 8e9ec f6). 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 1 . 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 2 and the related implications such as their potentially increased pathogenicity or transmissibility and, possibly, vaccine escape. Notwithstanding a unique proofreading activity among RNA viruses 3 , 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 4 . 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 necessary 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 5 , Lombardy 6 and Abruzzo 7 . 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 8 or rather the result of genetic drift due to the high level of genetic variation (http:// virol ogical. org/t/ respo nse-to-on-the-origin-and-conti nuing-evolu tion-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.

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.

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 Table 2 shows the most frequent amino acid substitutions stratified by lineage and clade. (Table 3) . Analysis of the international data set. 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 , 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 (Fig. 4 and Supplementary Fig. 2) , the analyses showed that B.1 probably originated in China and spread (Fig. 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 Fig. 3) .

The analysis of lineage B showed that only 2 nodes remained undetermined between Europe and China (Supplementary Fig. 4 ). The visualization (Fig. 5) 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 Fig. 5 ), highlighted the same ancestral scenario but did not show any introduction from Europe.

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 6, 9, 10 . 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 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 11 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 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 due to the relatively low evolutionary rate of SARS-CoV-2 in comparison with other RNA viruses, the limited number of sequences, available at the time of the study was performed, possibly affected by sampling bias for the characters underrepresented in several countries, including Italy, in the early epidemic and the rapid dissemination of the infection 3 . Furthermore, a frequent homoplasy 13, 14 , that affects multiple protein sites of the viral genome, and the founder effect played a dominant role in the early evolution of the virus 15, 16 . For these reasons, a limited number of studies on the SARS-CoV-2 phylogeography have been published, based on maximum parsimony 17 , maximum likelihood and Bayesian framework 3, 18 or phylogenetic network 19 . These studies analysed a limited number of genomes available at the time in the short interval elapsed since the origin of the virus 20 while others highlighted the importance of including travel-related information in the analysis 18 .

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 21 . 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. Future studies, employing all the isolates collected in the early phases of the epidemic, many of which became available only after this study began, could provide more information about the origin of the pandemic spread in Italy and Europe.

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).

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 pandemic. 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. All participants gave the written informed consent to the storage of their anonymised data in a protected database. All of the data used in this study were previously anonymized as required by the Italian Data Protection Code (Legislative Decree 196/2003) To place the Italian sequences in the context of the international COVID-19 pandemic, an additional dataset was built including Chinese (n = 52) and European (n = 858) sequences collected in the same period. Due to the large amount of European sequences available, we included at least 2 strains per country/week (Supplementary Table 2 and Supplementary Fig. 6) . Identical strains or those with more than 5% of gaps were excluded. For countries with a limited number of sequences, all strains were included. Only sequences reporting a certain sample collection were included. Consequently, the final dataset encompassed 1,375 sequences. SARS-CoV-2 sequences were aligned using MAFFT (https:// mafft. cbrc. jp/ align ment/ server/) and the alignment was manually cropped using BioEdit v. 7.2.6.1 (https:// bioed it. softw are. infor mer. com/). 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. megas oftwa re. net/).

Assigner tool v. 2.3.2 (last access 15 April 2021, https:// pango lin. cog-uk. io/) and Nextclade v. 0.14.1 (https:// clades. nexts train. 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 22 . 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 23 which allows the identification of clusters meeting given criteria. The Italian dataset was also analysed by BEAST v. 1.10 (https:// beast. commu nity/) in order to estimate the tMRCAs of the main clades. A previously estimated evolutionary rate of 8 × 10 −4 substitutions/site 24 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. Two runs of 150 million interactions were compared to assess convergence and then post-burnin samples pooled to summarize parameter estimates to obtain an effective samples sizes of at least 200. Finally, all trees were visualised in FigTree v. Each taxon was assigned to its sampling locations character, and we used PastML 21 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).

All consensus genomes are being submitted at the GISAID database (https:// www. gisaid. org). All the files used for the analyses are available upon request. 

Genomic characterization and phylogenetic analysis of SARS-COV-2 in Italy

Novel SARS-CoV-2 variants: The pandemics within the pandemic

The emergence of SARS-CoV-2 in Europe and North America

Addendum: A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology

Early Diffusion of SARS-CoV-2 Infection in the Inner Area of the Italian Sardinia Island

Genomic epidemiology of SARS-CoV-2 reveals multiple lineages and early spread of SARS-CoV-2 infections in Lombardy, Italy

Genomic epidemiology of the first wave of SARS-CoV-2 in Italy

Evolutionary analysis and lineage designation of SARS-CoV-2 genomes

SARS-CoV-2 shifting transmission dynamics and hidden reservoirs potentially limit efficacy of public health interventions in Italy

Molecular tracing of SARS-CoV-2 in Italy in the first three months of the epidemic

Suppression of a SARS-CoV-2 outbreak in the Italian municipality of Vo'

Transmission of 2019-nCoV infection from an asymptomatic contact in Germany

Stability of SARS-CoV-2 phylogenies

Emergence of genomic diversity and recurrent mutations in SARS-CoV-2

Quantitative phylogenomic evidence reveals a spatially structured SARS-CoV-2 diversity

A founder effect led early SARS-CoV-2 transmission in Spain

Phylogeography of SARS-CoV-2 pandemic in Spain: A story of multiple introductions, micro-geographic stratification, founder effects, and super-spreaders

Accommodating individual travel history and unsampled diversity in Bayesian phylogeographic inference of SARS-CoV-2

Phylogenetic network analysis of SARS-CoV-2 genomes

Transmission route and introduction of pandemic SARS-CoV-2 between China, Italy, and Spain

A fast likelihood method to reconstruct and visualize ancestral scenarios

IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era

Automated analysis of phylogenetic clusters

Early phylogenetic estimate of the effective reproduction number of SARS-CoV-2

We acknowledge the authors and the originating and submitting laboratories of the GISAID sequences. The research was conducted under a cooperative agreement between Università degli Studi di Milano-Medicina del Lavoro e Clinica delle Malattie Infettive del Dipartimento di Scienze Biomediche e Cliniche "Luigi Sacco", Intesa Sanpaolo and Intesa Sanpaolo Innovation Center.

The authors declare no competing interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-022-09738-0.Correspondence and requests for materials should be addressed to G.Z.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.