Transmission patterns of human enterovirus 71 to , from and among European countries , 2003 to 2013

C Hassel 1 , A Mirand 1 2 , A Lukashev 3 , E TerletskaiaLadwig 4 , A Farkas 5 , I Schuffenecker 6 , S Diedrich 7 , HP Huemer 8 , C Archimbaud 1 2 , H Peigue-Lafeuille 1 2 , C Henquell 1 2 , J Bailly 1 2 1. Clermont Université, Université d’Auvergne, EPIE, EA 4843, Clermont-Ferrand, France 2. CHU Clermont-Ferrand, Service de Virologie, Centre National de Référence des Entérovirus et Paréchovirus – Laboratoire associé, Clermont-Ferrand, France 3. Chumakov Institute of Poliomyelitis and Viral Encephalitides, Moscow, Russia 4. Prof. Gisela Enders & Kollegen MVZ GbR and Institute of Virology, Infectious Diseases and Epidemiology, Stuttgart, Germany 5. Division of Virology, National Center for Epidemiology, Budapest, Hungary 6. Laboratoire de Virologie Est des Hospices Civils de Lyon, Centre National de Référence des Entérovirus et Paréchovirus, Bron, France 7. National Reference Center for Poliomyelitis and Enterovirus, Robert Koch Institute, Berlin, Germany 8. Austrian Agency for Health and Food Safety, Vienna, Austria 9. Université Claude Bernard, EA4610 VIRPATH, Lyon, France


Introduction
The results of infections by enterovirus 71 (EV-71) range from the absence of symptoms to acute manifestations, including hand, foot, and mouth disease (HFMD) as well as neurological conditions such as acute meningitis, encephalitis, and poliomyelitis-like disease [1].The outbreaks caused by this virus since the late 1990s pose serious threats to public health in countries of eastern and south-east Asia (http://www.wpro.who.int/emerging_diseases/HFMD/en/) because a number of infections are involved in cardiopulmonary failure and neurological diseases which, in infancy (<1 yearold), can lead to death [2].During the HFMD epidemics in China between 2010 and 2012, there were 1,737 laboratory-confirmed deaths, of which 93% were associated with EV-71 infections [3].The occurrence of fatalities in Asia, notably in Cambodia, China, Malaysia, and Vietnam does not necessarily indicate that the EV-71 variants and strains circulating in these countries have a greater virulence [4].Rather, this probably reflects the high total number of infections (HFMD cases and asymptomatic infections).The overall disease burden caused by EV-71 led health authorities to develop surveillance systems in parts of eastern and south-east Asia.The World Health Organization provides a guide to clinical management and outbreak prevention of EV-71 infection (http://www.wpro.who.int/emerging_diseases/documents/HFMDGuidance/en/), for which five inactivated EV-71 vaccine candidates are being evaluated in clinical trials [5].
Unlike for some areas in eastern and south-east Asia, there is currently no particular epidemiological surveillance of EV-71 infections in the European Union (EU).The virus was involved in outbreaks of neurological diseases in the 1970s in Sweden (1973), Bulgaria (1975; 44 fatalities), Hungary (1978; 47 fatalities), and France (1979) [6][7][8][9].Sero-epidemiological surveys conducted on samples from 1999 to 2007 in Germany showed that ca 12% of children aged < 5 years, and ca 40 to 60% of adults develop neutralising antibodies against EV-71, which suggests that mild and mostly undiagnosed EV-71 infections occur in the general population [10,11].We and others showed that acute EV-71 infections are regularly investigated in patients admitted to hospital in a number of European countries, such as Austria, France, Germany and the Netherlands [12][13][14].Severe manifestations have been occasionally reported in France [15,16], and most documented EV-71 infections in EU Member States included sporadic cases of febrile illness and acute meningitis [17,18].
EV-71 is one serotype among a hundred human enteroviruses (EVs) that belongs to the Picornaviridae family.Virus strains can be divided into six genogroups designated A to F [19] but since the mid-1960s only genogroups B and C have been reported in outbreaks and individual cases of infection.The virus strains of the two major genogroups are classified into subgenogroups B0 to B5 and C1 to C5 on the basis of genetic relationships [20].Subgenogroups B4, B5, and C4 The chronogram trees were inferred with the 1D gene encoding the capsid viral protein 1 sequence datasets for subgenogroups B4/B5 (panel A), C1 (panel B), C2 (panel C), and C4 (panel D).The tree topologies show the genetic relationships between taxa sampled over the periods indicated on the x-axis (calendar years).The phylogenetic relationships were inferred with a Bayesian method using a relaxed molecular-clock model.For clarity, the sequence names are not included in the tree.Asterisks indicate key nodes with posterior probability (pp) density values > 0.90.
Each branch tip represents a sampled virus sequence.The branches in the genealogies of each subgenogroup were coloured to investigate relationships between phylogenetic clustering and the geographical origins of taxa.The geographical areas where the virus strains were sampled are indicated by different colours: blue, Europe (including European Union Member States, Iceland, Norway and Switzerland); light blue, Azerbaijan, Georgia, Kazakhstan, Kyrgyzstan Russia and Ukraine; green, North America; and red, Australia and eastern and south-east Asia.
are mainly restricted to Asian countries while subgenogroups C1 and C2 are chiefly found in Europe [20].Many aspects of the epidemiological and evolutionary dynamics of circulating virus strains remain unknown.
In particular, how virus transmission occurs over time, across space, and among genogroups has not been extensively studied in geographical areas other than eastern and south-east Asia [21].We investigated the evolutionary dynamics of EV-71 to determine the origin of virus strains sampled in Europe and their relationships with viruses reported elsewhere in the world.
To this aim, we used large sequence datasets of virus isolates sampled worldwide since the mid-1980s, which were mainly representative of EV-71 subgenogroups circulating in Europe and eastern and southeast Asia.We estimated the genetic diversity with the 1D gene encoding the capsid viral protein 1 (hereafter designated 1D VP1 ) over time, across geographical areas, and among subgenogroups, something that, to our knowledge, has never been done on such a large scale before.

Figure 2
Plots of genetic diversity a over time (calendar years) reconstructed for the enterovirus 71 subgenogroups C1 (panel A), C2 (panel B), and C1/C2 (panel C), 1985-2013 AUS: Australia; AUT: Austria; DEU: Germany; ESP: Spain; FRA: France; GBR: United Kingdom; HUN: Hungary; JPN: Japan; MYS: Malaysia; NLD: the Netherlands; NOR: Norway; SGP: Singapore; THA: Thailand; TWN: Taiwan; USA: United states.a Global genetic diversity over time was estimated with gene 1D encoding the capsid viral protein 1 (1D VP1 ) for each subgenogroup.Genetic diversity is estimated with the Bayesian skyline plot model and is expressed as log10Net (Ne: effective size of virus population; t: generation time).Genetic diversity reflects the effective number of infections averaged over time under the assumption of a neutral evolutionary process.On panels A and B the geographical areas where the virus strains were detected frequently are indicated above the peaks of genetic diversity.On panel C, plots estimated for the EV-71 subgenogroups C1 and C2 are superimposed on the same time scale.The epidemic peaks are numbered.The sporadic introductions of EV-71/C4 strains in European Union Member States (as described in this study and [34]) and Canada are indicated with open triangles.The introduction events in Russia are shown with full triangles.The introduction events of the EV-71 subgenogroup B5 in Denmark [30] and France [31] are shown with diamonds.

Data collection and compilation of nucleotide sequence datasets
The sequences determined were analysed with sequences obtained from GenBank.The 1D VP1 sequence datasets were constructed by collating all GenBank entries (as of July 2013 for the 1D VP1 locus) including a sequence of the VP1 capsid protein for any human isolate of EV-71.Entries reporting nt sequences of < 891 nt were discarded and only sequences with fully specified dates (year) and countries of origin were used.The sequence datasets were constructed with BioEdit v.7.2.5 software (http://www.mbio.ncsu.edu/bioedit/bioedit.html)and were compiled with all sequences determined in our laboratory.The EV-71 1D VP1 gene sequences were distributed into five datasets corresponding to subgenogroups B4/B5 (n = 217 sequences), C1 (n = 280), C2 (n = 322), and C4 (n = 675).On the basis of earlier phylogenetic data indicating that subgenogroups B4 and B5 had a common ancestor [21], we analysed jointly their sequences to increase the sample size.To investigate the EV-71/C4 subgenogroup, an initial dataset of 775 complete sequences was downsized to 675 by removing all sequences but one in genetic clusters containing multiple sequences with ≥ 99.5% nt identity with one another.

Coalescent estimation of divergence dates and evolutionary rates
Genealogical trees of EV-71 subgenogroups were reconstructed with the Bayesian evolutionary analysis by sampling trees (BEAST) v1.7.5 programme [23].Uncorrelated lognormal prior distributions of substitution rates among lineages were used for the molecular clock model [24].The evolutionary history was reconstructed with the substitution general time reversible (GTR) model.The phylogenetic parameters were coestimated in the different analyses by a Markov chain Monte Carlo (MCMC) process involving 100 x 10

Reconstruction of the demographic history and geographical spread of enterovirus subgenogroups
The Bayesian skyline plot model (BSP) does not assume a pre-defined model of demography [25].It infers a demographic parameter representing 'virus population size', also referred to as relative genetic diversity, using the coalescent theory and temporal information of the molecular data (virus collection date).Genetic diversity is the product of the effective size of the virus population (N e ) and the generation time (t), and reflects the effective number of infections averaged over time under the assumption of a neutral evolutionary process.In our study, this model was used to investigate possible variations in virus population size over time for subgenogroups B4, B5, C1, C2, and C4 to describe past epidemiological events on a continuous time scale using large sequence datasets for the 1D VP1 gene.EV-71 infections are characterised by epidemic transmission in a number of countries, a high rate of asymptomatic or pauci-symptomatic infections, and short-lived acute infections in susceptible individuals.Accordingly, the skyline plot model appears epidemiologically better suited than the other demographic models (constant population size, exponential growth, logistic growth, and expansion growth) for modelling the spread of EV-71 populations and thus for reconstructing the diffusion of HFMD over time and across geographical regions.The BSP analyses were done with a lognormal distribution (piecewise-constant population size model) to account for variation in substitution rates among the phylogenies.The countries where the virus strains were collected were used as discrete character states (or traits) in the phylogeographical analysis to estimate changes of geographical locations on the phylogenetic trees [26].The geographical locations of the ancestral nodes were co-estimated with the phylogeny by using a discrete phylogeographical diffusion model (symmetrical substitution model).Bayesian stochastic search variable selection was also used for identifying the statistically significant transition rates between locations.TreeAnnotator software was used to produce the MCC tree, with the branches coloured by traits.The MCC tree was obtained from a Bayesian MCMC analysis (discrete phylogeographical model) and shows for more clarity the most relevant parts of tree topologies (time is indicated as calendar years).The FigTree programme was used to display a number of information.
Branches and circles at the tree nodes are coloured according to the geographical location that had the highest probability.The size of each circle represents the location probability.Line width indicates the posterior probability of the corresponding lineage (thick lines indicate high posterior probability values).The tree topologies were used to determine the phylogenetic patterns (numbered sequentially) showing the most probable virus spread events.The main features of the most probable virus spread events (i.e.geographical location, location probability, and spread direction) are indicated in the figure.

Phylogenies of five enterovirus 71 subgenogroups estimated with datasets of the 1D gene encoding the capsid viral protein 1
We used the large number of 1D VP1 sequences determined in circulating strains to obtain a comprehensive comparison of virus transmission among EV-71 subgenogroups.The country of virus isolation was checked for each sequence to determine the distribution in four main geographical areas: (i) eastern and south-east Asia as well as Australia, (ii) Azerbaijan, Georgia, Kazakhstan, Kyrgyzstan, Russia and Ukraine (iii) Europe (whereby data originated from EU Member States, as well as Iceland, Norway and Switzerland) and (iv) North America, as observed in the genealogies.The phylogenies displayed a 'ladder-like' shape mainly characterised by a long trunk with short side branches, evidence of rapid extinction of lineages over time (Figure 1).These topologies arose from a combined effect of temporal sampling in various countries and rapid coalescence.A number of long side branches were also inferred in the C1, C2, and C4 trees, which indicated transmission of lineages over time and the persistence of these lineages alongside more recent ones.Of the subgenogroups examined, all but one (i.e.subgenogroup C1) were inferred to arise in the 1990s (Figures 1A, 1C, 1D).The time to most recent common ancestor (TMRCA) of subgenogroup EV-71/C1 was estimated in 1983.5 (95% HPD: 1982.2−1984.4) (Figure 1B).Most B4, B5 and C4 lineages were restricted to eastern and south-east Asia while subgenogroups C1 and C2 originated from this geographical area as well as European countries and Russia over the sampling periods.After 2000, several phylogenetically independent EV-71 C2 lineages were inferred to have circulated in Europe, during three main periods in 2007, 2010, and 2013.The C4 phylogeny also included a few lineages in EU Member States and Russia, which suggests that a number of persons in these countries sporadically acquired the virus from individuals elsewhere in the world, most probably eastern and south-east Asia.

Genealogy-based population dynamics of enterovirus 71 subgenogroups
The demographic histories showed variations between the five subgenogroups.The pattern estimated by our coalescent-based reconstruction for subgenogroup C4 clearly indicated an annual series of increases in genetic diversity interspersed with genetic bottlenecks.This pattern was caused by the large HFMD epidemics recently observed in China [3] and confirmed the usefulness of the BSP model for investigating EV-71 transmission.
We provide the detailed data obtained for the most frequent subgenogroups (C1 and C2) in Europe (Figure 2).The pattern estimated for subgenogroup C1 (Figure 2A) depicted successive variations in virus population size.This demographic pattern was consistent with a few outbreaks over time and across different geographical areas as documented by epidemiological data [20,27].This pattern may also indicate uninterrupted global circulation of the virus.
The pattern reconstructed for EV-71/C2 disclosed two distinct rises in virus population size (Figure 2B).The subgenogroup experienced a first sharp increase in 1998 that coincided with the occurrence of a large HFMD epidemic in Taiwan [28].The virus population size then slowly decreased over several years.The second exponential rise in virus population size (in 2007) was of a high magnitude.This major increase in virus genetic diversity coincided with the occurrence of an outbreak in the Netherlands [29] and a concomitant rise in virus circulation in Austria, France, and Germany [12].The transmission of the virus remained at a sustained high level after the 2007 epidemic, yet the different lineages sampled in 2010 and 2013 (Figure 1C) were not associated with increases in virus population size.
A combination of the C1 and C2 patterns showed that the decrease in the EV-71/C1 population size in 1996 was concurrent with the initial occurrence of EV-71/ C2 (Figure 2C).The 1998 increase in virus population size (Taiwan epidemic, subgenogroup C2, peak 1) coincided with a lower genetic diversity of EV-71/C1.The marked increase in C2 genetic diversity in 2007 (peak 5) also occurred in Europe during a period of low C1 diversity.The C1 pattern showed three yearly ascending rises (peaks 2, 3, and 4) during a long period of low genetic diversity of the C2 virus.The population size of both C1 and C2 subgenogroups remained low over the years 2004 to 2006 and 2008 to 2013, periods during which the EV-71/C4 subgenogroup occurred in several European countries (Austria, Croatia, France, Germany, and Hungary) and Canada (see below).EV-71/B5 caused an outbreak in Denmark (2007) [30] and a sporadic infection in France (2013) [31].

Spread of enterovirus 71/C1 between distinct countries
Spread of EV-71 between different countries was investigated with the MCC trees of the above analyses.A phylogenetic pattern indicative of a probable virus spread event between two nodes was defined as follows: (i) the nodes exhibited pp > 0.9, (ii) the inferred location probabilities were > 0.7 at both nodes, and (iii) the difference between the TMRCA values or the 95% HPD interval values of nodes were in a range of one year.With this conservative approach, the MCC tree topology inferred for subgenogroup C1 indicated a total of 13 consistent virus transmission events (Figure 3A).Eight phylogenetic patterns (numbered 03-05, 07-09, 11 and 12 in Figure 3A) indicated that the virus was transmitted between persons within various European countries.The United Kingdom, with a probability of 0.98, was likely the country of origin in three inferred virus spread events (numbered 03, 08, and 12, Figure 3A) to France (probabilities of 0.95) and one to Spain (0.99).France was inferred as the origin of two spread events, respectively to Germany and the United Kingdom.The EV-71/C1 tree also showed that virus spread occurred from Austria (location probability, 0.83) to Germany and from the Netherlands (0.84) to Finland.Three other phylogenetic patterns (numbered 06, 10, and 13 in Figure 3A) were indicative of long distance virus spread from a EU Member State to either Australia or Japan (n = 2 events).In the spread event number 13 (Figure 3A), the probability of the United Kingdom being the source country was estimated to be 0.56 against 0.28 for France.Two consistent phylogenetic patterns showed virus spread from Malaysia (location probability ≥ 0.94) respectively to Singapore and to Thailand.
A schematic representation of all virus spread events described above is shown in Figure 4.

Spread of enterovirus 71/C2 between distinct countries
In the MCC tree inferred for subgenogroup C2 (Figure 3B), 12 phylogenetic patterns indicated that the virus spread between distinct European countries, Georgia, and Russia.The most frequent likely countries of origin of infected persons in the inferred virus spread events were the Netherlands (n = 5 events; location probability range: 0.87-1) and France (n = 5; 0.74-0.98)followed by Finland and Germany.The 2007 epidemic lineage in the Netherlands was involved in only one virus spread event (Figure 3; number 7).Three virus spread events were estimated between distant countries: from Finland to Georgia and from France to Russia.An additional event of spread of the virus from Japan (location probability: 0.85) to Spain was also inferred in the genealogy before 1995 (not shown in Figure 3B).All the virus spread events occurring after 1998 are shown in Figure 4.

Spread of enterovirus 71/C4 between distinct countries
We analysed separately a sequence subset corresponding to the lineage highlighted in Figure 1D.The distributions of geographical location probabilities were indicated for the most relevant nodes in the MCC tree inferred for this sequence subset and were used to investigate virus spread between countries (Figure 5).The phylogeny pattern indicated that between 2001 and 2002 a C4 virus strain spread from China to Japan, and from there to the EU in early 2003.The virus was disseminated in different EU Member States between 2003 and 2004.It was sampled in Croatia in 2005 and the phylogenetic analysis estimated a probability of 0.92 that the source country was Germany.The same virus strain also moved from the EU to Canada (sampling year 2006).The phylogenetic analysis suggested that the virus was spread from Germany with a probability of 0.72.

Discussion
The recent reporting of sporadic circulation of EV-71 strains of the 'Asian subgenogroups' B5 [30,31] and C4 [32][33][34] in European countries makes it important to determine the evolutionary origins and spatiotemporal spread of EV-71 infections in geographical areas other than eastern and south-east Asia.There is a large body of literature on the molecular epidemiology of EV-71 [35] but, to our knowledge, the spread of this virus has not yet been analysed explicitly with recently developed phylogeographical models [26].In this report, we therefore applied demographic and phylogeographical analyses to compare the population structure and the dissemination patterns of EV-71 subgenogroups in Europe and Asia.We found that the EV-71 subgenogroups exhibit distinct phylodynamic patterns and that these patterns differed from those of past epidemics in frequency and geographical location.The transmission history of subgenogroup C2 was characterised by fewer increases in genetic diversity than in subgenogroups C1 and C4.For instance, the genetic diversity pattern of subgenogroup C2 shows clearly the two outbreaks (Taiwan, 1998 and the Netherlands, 2007) documented in the literature [28,29].The variations in genetic diversity over the study period also reflect gaps in virus sampling over time and surveillance differences between countries.In Europe, EV-71 sequences are determined solely in patients admitted to hospitals and in a few countries only, whereas in a number of Asian countries real-time epidemiological data are obtained from sentinel surveillance of HFMD cases.In this study, we used only full-length EV-71 1D VP1 sequences because phylogenetic trees based on short-length sequences did not provide a high level of statistical confidence.Partial nt sequences were notably determined for virus genotyping and ranged over different parts of the 1D VP1 gene.Accordingly, we had to discard a certain amount of sequence data, some from Europe, which might have been helpful in phylogenetic analysis of virus spread.
We show how transmission of EV-71 strains (subgenogroups C1, C2, and C4) occurs in Europe as discrete and temporally defined virus introductions, occasionally followed by limited local spread.The only exception to this pattern was the epidemic expansion in the Netherlands in 2007.Interestingly, we found limited phylogenetic evidence that the Dutch outbreak was the source of a large spread to other European countries.The phylogenetic patterns also show that European countries may experience multiple virus introduction events within the same year.This dissemination mode was observed within particular countries (e.g.France, Germany) for both subgenogroups C1 and C2 but was more clearly seen during three waves of infections in 2007, 2010, and 2013 (subgenogroup C2).The epidemiological and biological factors involved are unknown but the occurrence of these infection waves is consistent with the hypothesis that the sustained circulation of an EV-71 strain throughout Europe depends on the proportion of susceptible hosts in different countries.The data may alternatively indicate that the immunity elicited by the C1 and C2 infections is cross-protective, as suggested by earlier studies that identified common epitopes and suggested that human immune sera among virus strains of different EV-71 subgenogroups have cross-neutralization properties [36].In this respect, we also provide phylogenetic evidence that the spread of the C2 virus across European countries in 2007 happened at a time when C1 infections had been transmitted at low rates for at least three years.
It is particularly noteworthy that the transmission of EV-71 strains C1 and C2 in Europe is mainly dependent on the frequency of virus spread events between neighbouring countries.This has been previously described for EV-71 [37] and for another EV, coxsackievirus B5 [38].The present study also indicates that the long-term survival of EV-71 (C1 and C2) depends on continued virus transmission between individuals across larger geographical areas, notably Russia and Asia.We also reconstructed a consistent transmission chain caused by a C4 virus strain throughout Europe (Austria, Croatia, France, Germany, and Hungary), which suggests that the virus persisted between 2003 and 2005.The chain was brought about by a virus from Japan which eventually reached Canada between 2005 and 2006.Transmission of the C1 and C2 infections was low during the whole period in Europe.Similarly, the recent C4 infections in France (2012), and Russia (2013) arose during a transmission trough of both C1 and C2 infections.Thus, transcontinental transmission events between individuals should be considered in the global epidemiology of the virus: they provide the epidemiological bases for explaining the long-term survival of lineages despite abrupt extinctions that drastically reduce virus diversity after the occurrence of outbreaks in particular locations.
Earlier hypotheses proposed that yearly epidemics of EV-71 in Japan arose from spread of the virus from neighbouring countries [37].We suggest that virus strains are sustained by complex migration dynamics involving Europe and Russia, and possibly other geographical regions.This argument is lent weight by a recent study reporting the occurrence of EV-71 C4 in Denmark in 2012 and 2013 [39].Our Bayesian analysis of the partial 1D VP1 sequences showed that the C4 viruses isolated in Denmark were distributed in two distinct phylogenetic clusters, which suggests independent virus sources (data not shown).Thus, the frequency of epidemics in Asian countries, the involvement of EV-71 in neurological conditions, the recent isolation of C4 variant strains in Europe, and the spread patterns of the virus between distant and neighbouring countries underline a need for enhanced surveillance of EV-71 infections in Europe using common enterovirus genotyping methods.

Figure 3
Figure 3 Spatial and temporal distribution of enterovirus 71 lineages showing the geographical dissemination of subgenogroups C1 (panel A) and C2 (panel B)

Figure 5
Figure 5Bayesian analyses of the spatial and temporal spread of an enterovirus 71 subgenogroupC4 lineage, 2002-2006

Methods Virus samples, gene amplification, and nucleotide sequencing
.ox.ac.uk/Evolve/Software.html) was used to check for convergence and mixing (operator effective sample size > 200).The trees estimated with the MCMC procedure were sampled to obtain a final 20,000 trees.Mean estimates and 95% highest probability density (HPD) intervals calculated for each operator were compiled by analysing the output files obtained from the BEAST programme with the Tracer v.1.5programme.Maximum clade credibility (MCC) trees were calculated with the TreeAnnotator v.1.7.5 programme (http://evolve.zoo.ox.ac.uk/Evolve/Software.html) and topological support was assessed by estimating the values of the posterior probability (pp) density of each node.