Skip to content

Advertisement

  • Research Article
  • Open Access

Evolutionary phylogeography and transmission pattern of echovirus 14: an exploration of spatiotemporal dynamics based on the 26-year acute flaccid paralysis surveillance in Shandong, China

  • 1,
  • 2,
  • 2,
  • 2,
  • 2,
  • 2,
  • 2,
  • 1,
  • 1 and
  • 1, 2Email author
Contributed equally
BMC Genomics201718:48

https://doi.org/10.1186/s12864-016-3418-3

  • Received: 16 September 2016
  • Accepted: 13 December 2016
  • Published:

Abstract

Background

Echovirus 14 (E-14) causes various clinical recognized syndromes, mostly with gastrointestinal syndrome and paralysis. The current study summarized the Shandong E-14 strains isolated from a 26-year acute flaccid paralysis (AFP) surveillance, and elucidated the characterization of phylogenetic and phylogeographic relationships of E-14 worldwide.

Results

As a predominant serotype circulating in AFP surveillance, phylogenetic analysis showed that E-14 exhibited both time and geographic subdivision worldwide. In order to know the evolutionary history and spatial temporal dynamics of E-14, evolutionary phylogeography was reconstructed using BEAST and SPREAD software based on the VP1 sequences. The time of the most recent common ancestor of E-14 was estimated around 85 years and evolved with 9.17 × 10−3 substitutions/site/year. Phylogeographic analysis suggested that two regional transmissions of E-14 were mainly detected, with one located between Europe and Africa countries and the other was in the Asia-Pacific region.

Conclusions

Our study investigates the molecular evolution and phylogeographic of E-14, and brings new insight to the dispersal of E-14 worldwide. Regional transmission was mainly detected and Australia may be responsible for the spread of E-14 in recent years.

Keywords

  • Echovirus 14
  • Acute flaccid paralysis
  • Phylogeny
  • Phylogeography

Background

Enteroviruses (EVs) are common human pathogens belonging to the genus Enterovirus, family Picornaviridae. Most EV infections are usually asymptomatic, but sometimes they are associated with diverse clinical syndromes ranging from minor febrile illness to severe, potentially fatal diseases such as encephalitis, paralysis, myocarditis and neonatal enteroviral sepsis [1, 2]. EVs have traditionally been detected based on virus isolation in cell cultures and subsequently serotyped by serum neutralization test, and were classified into echoviruses (E), coxsackieviruses (CV) group A and B, and polioviruses. Currently, based on the phylogenetic relationships in multiple genome regions, International Committee on Taxonomy of Viruses reclassified EVs into 12 species: EV-A to EV-H, EV-J and rhinovirus A, B and C [3]. Of these species, seven (EV-A to EV-D and rhinovirus A, B, C) of the EV species infect humans while remaining infect non-human animals.

Echovirus 14 (E-14) was first isolated from a hospitalized patient suffering from aseptic meningitis syndrome in 1954 [4] and had been classified as a member of EV-B. In contrast to other common EVs such as EV-A71 and CV-A16, both of which are categorized as EV-A and are major etiological agents for the widely spread hand, foot and mouth disease (HFMD) [57], E-14 is considered relatively rare as it less likely to cause worldwide occurrences of clusters and outbreaks. However, epidemics of E-14 infection had been observed in Europe and Asia in the 1960s [811]. Similarly, association of E-14 with many clinical recognized syndromes such as gastrointestinal syndrome [810], acute serous meningitis [12], paralysis [13], fatal neonatal hepatic necrosis [14], aseptic meningitis [11, 15] has been reported previously. In recent years, E-14 is more than detected in an outbreak or epidemic caused by other pathogens [16, 17] or exhibited sporadic infection in conventional virological surveillance [1820].

As a coastal province, Shandong has a population of about 96 million and major ports that could potentially serve as portals for importation of exogenous viruses. The acute flaccid paralysis (AFP) surveillance, developed for polio eradication program, was launched in Shandong in 1988. Currently, the AFP surveillance has been conducted in all 138 counties of Shandong province and more than 600 sentinel hospitals were included. AFP surveillance can obtain non-polio enterovirus (NPEV) strains as a side benefit and produce baseline data of local NPEV distribution and a genetic overview [21]. Our previous study demonstrated continuous AFP surveillance provided valuable information on the circulation and emergence of different EV types in the context of limited EV surveillance in China [20]. In this present study, we investigated the genetic characterization of E-14 strains isolated from Shandong AFP surveillance during 1988–2013. Phylogenetic and phylogeographic analysis were combined to elucidate the evolution relationship and spatiotemporal dynamics of E-14 worldwide.

Results

Virus isolation

The annual numbers of reported AFP cases (contacts), NPEV and E-14 isolates were illustrated in Fig. 1. From 1988 to 2013, a total of 962 NPEV isolates were identified based on the molecular typing method in Shandong AFP surveillance, with E-14 serotype having the highest number of the total strains. In our study, E-14 was first detected in 1993. As the AFP surveillance and report in China was not so active in the early 1990s, no E-14 isolates were detected for three consecutive years (1994–1996), and numbers of E-14 and NPEV strains increased substantially since 1997. It was worth noting that the proportion of E-14 in contacts is higher than in AFP cases is due to the formation of EV transmission chains, for the rest contacts were probable positive as one infected E-14. The mean constituent rate of E-14 was 7.6% (73/962), and the highest constituent rate reached to 23.1% (9/39) in 2004. Consequently, E-14 under EV-B persisted to be the most commonly isolated serotype circulating in a long term trend in our AFP surveillance.
Fig. 1
Fig. 1

Annual numbers of AFP cases (contacts), NPEV and E-14 isolates in Shandong AFP surveillance from 1988 to 2013

Homology and phylogenetic analysis

Excluding the same (100% VP1 nucleotide homology) strains and co-infection strain, we used 73 E-14 strains isolated from Shandong province for phylogenetic analysis based on the VP1 sequence. Homologous comparison showed 73 Shandong E-14 isolates had 79.8–99.8% (94.9–100% amino acid sequence identities) VP1 nucleotide sequence identities among themselves and the overall mean p-distance value was 0.142. Phylogenetic analysis was conducted on partial VP1 coding region (nt 2596–2880 on prototype strain Tow) of Shandong isolates with global reference E-14 sequences. Overall, the global E-14 sequences were distributed into two major partitions identified by lineage I and II. Compared with the global E-14 strains, all Shandong E-14 isolates fell into a major phylogenetic lineage (lineage I) that included strains from China (Fig. 2). This is probably due to the more frequent population exchange in domestic than the rest of the world, as the international exchange needs to perform a more cumbersome procedure in China. Therefore, isolates from China studied here existed distinct genetic relationship with other worldwide E-14 strains. We also suggested that E-14 strains from China had the same putative ancestor relative to the strains from other regions of the world.
Fig. 2
Fig. 2

Phylogenetic tree based on partial E-14 VP1 sequences of Shandong isolates and global strains. Strains are color-coded and represented group A-F according to their phylogenetic relationship. Circles indicate strains detected in the present study; triangles indicate strains isolated from other provinces in China. Abbreviations in isolate names are as follows: CHN-SD, Shandong province in China; USA, United States of America; AUS, Australia; IN, India; FRA, France; BAN, Bangladesh; MAD, Madagascar; NL, Netherlands; NIE, Nigeria; CAF, Central African Republic; GRE, Greece; GEO, Georgia; KEN, Kenya; PAK, Pakistan

Within lineage I, all sequences segregated into three distinct genomic groups (A-C). Homologous analysis revealed 79.9–88.7, 80.8–84.7, 79.8–86.3% nucleotide identities between groups A and B, A and C, B and C, respectively. Strains in these three groups showed high correlation with the time of isolation, and a limited temporal overlap was observed among these three groups that circulated in Shandong and other provinces in China. The phylogenetic analysis showed that strains in lineage II mainly from India, Australia, Georgia, Netherlands, France and Bangladesh. For the 26 India strains isolated from 2007 to 2011, 20 strains segregated into group F with strains from Australia, while the other 6 strains isolated from 2008 to 2011 belonged to group E in the phylogenetic tree.

Phylogeographic analysis

In order to estimate the divergence time and substitution rate of E-14 more accurately, the complete VP1 sequence from Shandong province and from global were selected to analyze with the Bayesian MCMC method. Different models were used for data analysis, and all of the results were summarized in Table 1. Considering both the Harmonic Mean and the AICM model selection results, we concluded that the UCLD relaxed clock model was a significantly better fit model than the strict clock model for each of the parameter. Hence, the combination of UCLD relaxed clock model and exponential growth (EG) model were recommended for this study, and the estimated nucleotide substitution rate for global E-14 was 9.17 × 10−3 substitutions/site/year (95% HPD: 7.67 × 10−3–10.80 × 10−3). Estimated the t MRCA of global E-14 strains was August 1928 (95% HPD: February 1910–May 1944), China strains (lineage 1) was June 1972 (95% HPD: February 1965–June 1979), and the strains isolated in other regions of the world (lineage 2) was September 1967 (95% HPD: January 1958–August 1976).
Table 1

Estimators of phylogeographic analysis on E-14 strains under different molecular clock and coalescent model combinations

Parameter

Estimator of the parameter (95% HPD interval) under different molecular combinations

SC & Const

SC & Exp

SC & BSP

UCLD & Log

UCLD & Exp

UCLD & BSP

Harmonic Mean

−18271.839

−18253.871

−20602.113

−18187.626

−18169.821

−20529.023

AICM

36747.72

36686.69

41416.886

36765.743

36742.314

41489.35

t MRCA (global E-14)

88.06 (78.14,98.49)

87.74 (77.86,97.49)

86.29 (77.15,96.05)

85.48 (70.22,103.77)

85.31 (69.64,103.79)

83.02 (68.20,100.25)

t MRCA (Lineage 1)

43.44 (39.22,47.96)

43.52 (39.36,48.18)

42.64 (38.81,46.97)

41.26 (34.69,48.50)

41.46 (34.45,48.79)

40.80 (34.92,47.36)

t MRCA (Lineage 2)

46.87 (41.54,52.31)

47.04 (41.77,52.62)

45.92 (41.07,51.20)

46.32 (37.09,56.15)

46.15 (37.33,55.94)

44.41 (36.14,53.59)

t MRCA (group 1)

32.31 (28.95,35.75)

32.36 (29.09,35.90)

31.73 (28.66,34.97)

31.50 (26.51,37.04)

31.46 (26.35,36.96)

31.38 (26.61,36.66)

t MRCA (group 2)

29.85 (26.57,33.27)

29.89 (26.70,33.38)

29.45 (26.36,32.67)

27.41 (22.47,32.85)

27.40 (22.49,32.77)

28.09 (23.16,33.32)

t MRCA (group 3)

29.63 (27.53,31.72)

29.70 (27.70,31.97)

29.53 (27.58,31.60)

28.45 (25.42,31.64)

28.46 (25.26,31.61)

28.85 (25.95,32.21)

t MRCA (group 4)

46.43 (40.44,52.43)

46.62 (40.92,52.85)

45.63 (40.19,51.24)

45.47 (36.29,56.13)

45.30 (35.46,55.05)

43.72 (35.19,53.39)

t MRCA (group 5)

38.91 (33.45,44.89)

39.05 (33.61,44.99)

38.33 (33.06,43.68)

34.62 (27.34,43.15)

34.43 (27.34,42.15)

33.55 (27.35,41.16)

Mean evolutionary rate

8.44 (7.42,9.48)

8.38 (7.39,9.46)

8.57 (7.61,9.58)

9.14 (7.65,10.70)

9.17 (7.67,10.80)

9.30 (7.80,10.80)

The best fitting model combination was highlighted in boldface. The molecular clock models used were the strict molecular clock (SC) and the uncorrelated log-normal distributed (UCLD) relaxed molecular clock. Coalescent models used were the parametric Constant model (Const), Exponential model (Exp), Logistic model (Log) and Bayesian skyline plot model (BSP). 95% HPD 95% highest posterior density, t MRCA the most recent common ancestor, AICM the Akaike information content model selection method

Results of the MCC tree (Fig. 3) constructed by Tree Annotator showed that Shandong E-14 isolates segregated into three groups, which were consist with the phylogenetic tree in Fig. 2. Phylogenetic analysis suggested the origin of global E-14 could be in the USA. Lineage 1 represented strains from China (including all Shandong isolates), and the data indicated that E-14 in China shared the same putative ancestor. Lineage 2 represented the strains from other regions of the world, and the strains were gradually radiated into France, India and Bangladesh. Our results showed that strains in Australia may have also played a dominant role in spreading this lineage.
Fig. 3
Fig. 3

Maximum clade credibility (MCC) tree of the complete VP1 sequences of E-14 strains throughout the world visualized in FigTree. The colors of the branches corresponded to their probable geographic location. The intervals of branch reflect the 95% highest posterior density (HPD) intervals for the branch height. E-14 isolates from Shandong segregated into three groups (1, 2 and 3). Abbreviations of geographic location are shown as described in the legend of Fig. 2

The geographic dispersal of the global E-14 strains was conducted with the partial VP1 sequences and was shown in Fig. 4 (Additional files 1 and 2). Phylogeographic reconstruction was able to identify a single location for the origin of the global E-14 strains in the USA around March 1928. The virus crossed the Atlantic Ocean and arrived in Greece in 1951, and then spread following two distinct routes: i) to the East arrived in Australia around 1958 and ii) to the South arrived in Central African Republic around November 1978. As E-14 strains arrived in Australia, the virus was latent in this region and spread to India, Georgia, Bangladesh and France in about 40 years. Viral strains that from Australia arrived in Guangdong province of China in an early time, and then showed a relatively quick diffusion rate to Shandong and Zhejiang provinces. In general, E-14 strains from the USA and Australia were the most probable source for the introduction of E-14 strains to European countries and the Asia-Pacific region, respectively, with all BF > 19.0.
Fig. 4
Fig. 4

Geographic distribution and inferred dynamics of global E-14 strains. The map was reconstructed using the ArcGIS (http://www.esri.com/), and was identical to the original image created by the SPREAD and GoogleTM Earth. Arrows indicate the inferred routes of spread of global E-14. The number represented the time when E-14 arrived. Abbreviations of geographic location are shown as described in the legend of Fig. 2

Discussion

Compared with many EV serotypes which can lead to serious illness, E-14 is considered relatively rare as it less likely to cause worldwide occurrences of clusters or outbreaks. To the best of our knowledge, no detailed molecular epidemiology study was carried out based on a large number of E-14 sequences. Owning to the launch of the Global Poliomyelitis Eradication Initiative, we had accumulated a significant number of E-14 sequences associated with AFP surveillance which was developed for casting a wider net for poliovirus detection in Shandong since 1988. For the detection of circulating EVs, AFP surveillance has also provided a major source of information in China, as specialized EV surveillance based on other approaches or populations was still very limited. In addition, studies conducted in other countries [18, 2227] also described the NPEV serotypes identified from the AFP surveillance, and the information had attracted more attention of the scholars worldwide.

Although E-14 may cause mild to severe neurological disorders as reported from Russia [12], France [9, 10], Italy [8] and Japan [11, 13, 15], unfortunately, no clinical information, molecular epidemiology and disease correlation of this serotype have been reported. This study, which was performed a comprehensive molecular epidemiologic analysis of Shandong E-14 isolates, was feasible because of the large number of AFP database during the last 26 years in Shandong, China. During 1990–2013, a total of 10322 AFP cases and contacts were reported, and 958 NPEV isolates were obtained from stool specimens in Shandong AFP surveillance [20]. A total of 53 NPEV serotypes were identified, serotypes such as E-6, E-14 and EV-A71 were the most commonly reported. It should also be notified that no EV-D68 strains were detected in our AFP surveillance, although this EV serotype was associated with AFP and other neurological complications. The proportion of E-14 isolated in this study was high (7.6%) as compared to other identified NPEV serotypes [20]. Additionally, E-14 was also found a prevalent EV serotype with AFP surveillance in other regions of the world [18, 19, 28]. Consequently, all these findings suggest E-14 as the possible causative agent of paralysis, and the molecular epidemiologic analysis of E-14 is still necessary.

EVs are well known for their high genome plasticity, due to both high mutation and recombination rates. In this study, we analyzed the global E-14 genetic diversity, of 73 isolates from Shandong and 54 strains isolated in twelve different countries. The phylogenetic tree (Fig. 2) revealed that E-14 strains clustered in both time and regional distributions. For E-14 isolated in China, strains in group C and B were from 1993 to 2006 and 2000 to 2010, respectively and there were 10 and 6 years had passed since the last isolation in these two groups. The isolation years of group A spanned from 1993 to 2014, revealing a long time circulation of this group in Shandong region. Moreover, it is obvious that E-14 isolates belonging to group B and C may have been disappeared after circulating for a period of time and was replaced by group A. Viruses in group A may have become the dominant group circulating in China. Alternatively, the information on Chinese E-14 isolates is limited, and more sequences are needed to determine domestic E-14 phylogenetic relationships.

We used a phylogenetic molecular clock approach to further infer the historical background and geographic diffusion of global E-14 strains. The MCC tree (Fig. 3) constructed was congruent with the results observed in phylogenetic analysis (Fig. 2). It was perceived by Bayesian MCMC evolutionary analysis that E-14 strains were circulated for about 80 years in the world. Moreover, E-14 strains isolated in China may be the product of geographic segregation and has evolved independently in China for a long time. Two observations suggest that the common ancestor of global E-14 lineages existed in the USA. First, USA was the first country reported to detect E-14 [11] and second, the USA strain is phylogenetically more diverse within the phylogenetic tree (Fig. 2) and the MCC tree (Fig. 3) than those from elsewhere. However, these observations may reflect differences in AFP surveillance intensity among countries and more information is required before we can conclude that E-14 was introduced to the USA multiple times from other regions. Interestingly, our results indicate that circulation of E-14 in Australia is very dynamic, as the virus likely to be the reservoir for evolutionary origin of E-14 in the world. The results that the nucleotide substitution rate of E-14 VP1 gene was similar to other EV serotypes [2934] also suggested the evolutionary pressure acting on E-14 existed no significant difference with other serotype strains.

A particular focus of this study was to determine the phylogeographic relationship between E-14 collected from Shandong province and compare these with those from other countries. To our knowledge, the current study is the first applying the Bayesian phylogeographic reconstruction approach to gain insight into the emergence and spread of E-14 strains worldwide [35]. The phylogeographic analysis revealed that the transmission of E-14 in China likely originated from the Australia, and then the endemic circulation was observed without spreading to foreign countries. However, E-14 strains in China and Australia have no direct evolutionary relationship according to the phylogenetic analysis. The reason why the phylogenetic and phylogeographic analysis was inconsistent between E-14 strains in China and Australia remains unclear, and analysis of more sequences with greater detail information will yield more precise diffusion process. In addition to spread from Australia to India, E-14 strains seem to have been reintroduced into Australia from India, and the same scene occurred among Shandong province with Guangdong and Yunnan provinces in China. We believe that large-scale patterns in human mobility between these locations will provide a more plausible explanation on E-14 emergence and introduction. Overall, as a well documented origin from the USA, E-14 in this area has become silent over the past few decades. Australia may be responsible for the transmission of E-14 in recent years, and two regional transmissions were mainly detected, with one located between Europe and Africa countries and the other was in the Asia-Pacific region.

The results in this study are subject to a selection bias. Sequences from GenBank are limited, especially for the complete VP1 genes of E-14. In addition, the first E-14 sequence represents a strain isolated in 1954 in the USA, whereas the next sequences can be obtained from GenBank originated in 1980. Therefore, partial VP1 sequences were used for the ananlysis of the phylogenetic relationships among global E-14 strains. Due to availability of limited gene information of global E-14, further collection and analysis, together with detailed epidemiological information is required to understand the spatiotemporal dynamics of E-14 worldwide. All these indicate that EV surveillance in most countries remains to be strengthened and the systematic surveillance are necessary in investigating the EV circulation and understanding the social burden of EV infection.

Conclusions

In this study, we collected the complete VP1 sequences of E-14 based on the 26-year acute flaccid paralysis surveillance in Shandong province, and inferred the historical background and geographic diffusion of global E-14 strains. The estimated nucleotide substitution rate for global E-14 was 9.17 × 10−3 substitutions/site/year, and the estimated t MRCA was August 1928. Our research confirmed that the origin of E-14 could be in the USA, however, regional transmissions were detected with one located between Europe and Africa countries and the other was in the Asia-Pacific region. Australia may be responsible for the spread of E-14 in recent years.

Methods

Virus isolation

Stool samples were collected from children aged <15 years who were presenting acute flaccid type of paralysis and the contacts were collected, and processed according to standard procedures recommended by the WHO [36]. Contacts were defined as high-risk AFP cases with (1) <5 years of age, (2) <3 doses of oral poliovirus vaccine (OPV) immunization or unknown OPV history, (3) inadequate stool specimen or (4) diagnosis of clinical poliomyelitis. The AFP surveillance has been ongoing without any changes since 1990 in line with the WHO recommended standards. Three cell lines, human rhabdomyosarcoma (RD) cell, human epidermoid carcinoma (Hep-2) cell and mouse cell expressing the gene for the human cellular receptor for poliovirus (L20B) were used for virus isolation [36]. All cell lines were gifted from the WHO Global Poliovirus Specialized Laboratory in the USA and were all originally purchased from the American Type Culture Collection. A total of 200μl of chloroform-treated stool solution was added to each of the cell culture tubes, and the inoculated cells were examined daily. After 7 days, the tubes were frozen, thawed, and re-passaged and another 7-day examination was performed. Infected cell cultures were harvested and used for further examination until complete cytopathic effect (CPE) was obtained. To avoid cross contamination, cell tubes of normal L20B, RD and Hep-2 cells served as negative controls. Isolates from RD or Hep-2 cell lines were re-passaged to L20B cell line, and were designated as NPEV if no CPE was observed.

Nucleotide sequencing and molecular typing

Total RNA was extracted from 140μl of the infected cell culture using QIAamp viral RNA mini kit (Qiagen, Valencia, CA, USA) according to the manufacturers’ recommended procedure. Primer pairs 008/013 [37] and 187/011 [38] were used to amplify VP1 gene and the combination of the two sequences yielded the entire VP1 coding region. The reverse transcription-PCR (RT-PCR) was performed using Access RT-PCR System (Promega, USA) according to the manufacturers’ procedures. In order to avoid cross contamination, a RT-PCR reaction using the RNA extracted from normal RD cell served as a blank control, and a negative control containing all the components of the reaction except for the template was also included.

PCR amplicons were visualized in agarose gels, and positive products were purified and sequenced with the BigDye Terminator v3.0 Cycle Sequencing kit (Applied Biosystems, Foster City, CA). The nucleotide sequences were analyzed by ABI 3130 genetic analyzer (Applied Biosystems, Hitachi, Japan). The PCR products were sequenced in both directions to avoid possible ambiguous nucleotides. Molecular typing based on VP1 sequences was performed using BLAST, obtained online from the National Center for Biotechnology Information (NCBI). Isolates showing >75% nucleotide sequence identities with the E-14 Tow prototype strain were designated relative serotype [37].

Homologous comparison and phylogenetic analysis

VP1 sequences of E-14 isolated in other countries were blasted and collected by using the nucleotide program in the NCBI. All of the references were downloaded with the strain name, GenBank number, geographical region and collection date. Nucleotide sequence alignments were carried out by BioEdit software v7.0.5.3 [39]. Phylogenetic trees were constructed by using MEGA v6.0 [40], using the neighbor-joining method after estimation of genetic distance using the Kimura two-parameter method [41]. A bootstrapping test was performed with 1,000 duplicates, and the transition/transversion rate was set at 2.0. To the sequences detected in worldwide, we chose one sequence to represent sequences with similarity over 99 percentage and identical isolation years. Recombination can sometimes bias evolutionary relationship when constructing phylogenetic trees [42], hence recombination of VP1 sequences was excluded based on the analysis of SimPlot.

Bayesian Markov Chain Monte Carlo (MCMC) method was used to analyze the evolution rate and molecular clock phylogeny of global E-14 strains containing Shandong isolates with the BEAST software package version 1.7.5 [43], and the time of the most recent common ancestor (t MRCA ) with 95% highest posterior density (HPD) of global E-14 and Shandong isolates were estimated (Additional file 3). The calibration point was the year that each strain was isolated. We used the general time reversible (GTR) nucleotide substitution model with gamma-distributed rates of variable among sites, which were identified as the best fitting model by jModelTest v2.1.7 on the basis of Akaike Information Criterion (AIC) [44]. Runs were performed using the constant size, under the strict clock model and the uncorrelated log-normal distributed (UCLD) relaxed molecular clock model. Multiple combinations of molecular clock and coalescent models were explored, and to select the best fitting model we used posterior-simulation Akaike information content (AICM) [45] and Harmonic Mean [46] model selection methods described by Raftery and Baele, respectively. Bayesian MCMC analyses were run with a chain of 30 million steps and sampled every 1,000 steps. Convergence of parameters was identified by TRACER v1.5 with the effective sample size (ESS) exceeding 200. The Maximum Clade Credibility (MCC) tree was calculated by TreeAnnotator with a burn-in period of 3,000 and then visualized in FigTree v1.4.2.

Phylogeographic analysis

In order to infer the spatiotemporal dynamics of E-14, the Bayesian Stochastic Search Variable Selection (BSSVS) was used to provide evidence for statistically supported diffusion between state variables under BEAST v1.7.5. The results of BSSVS was summarized using SPREAD v1.0.6 [47], a keyhole markup language (KML) file was generated to identify the major routes of geographic diffusion (Additional file 4). Bayes factor (BF) test was used to select the most probable diffusion process. To visualize the geographic dispersal of E-14, the KML file was imported to Google™ Earth to produce a graphical animation of the estimated spatiotemporal pathways of global E-14.

Nucleotide sequence accession numbers

The VP1 sequences of E-14 isolates described in this study were deposited in the GenBank database under accession numbers KX840391-KX840454. All sequences data that support the findings of this study have been deposited in TreeBASE and are accessible through the accession number 20306 (http://purl.org/phylo/treebase/phylows/study/TB2:S20306).

Abbreviations

AFP: 

Acute flaccid paralysis

AIC: 

Akaike information criterion

AICM: 

Akaike information content

BF: 

Bayes factor

BSSVS: 

Bayesian stochastic search variable selection

CPE: 

Cytopathic effect

CV: 

Coxsackieviruses

E-14: 

Echovirus 14

EG: 

Exponential growth

ESS: 

Effective sample size

EVs: 

Enteroviruses

GTR: 

General time reversible

Hep-2: 

Human epidermoid carcinoma

HFMD: 

Hand, foot and mouth disease

HPD: 

Highest posterior density

KML: 

Keyhole markup language

L20B: 

Mouse cell expressing the gene for the human cellular receptor for poliovirus

MCC: 

Maximum Clade Credibility

MCMC: 

Bayesian Markov Chain Monte Carlo

NCBI: 

National Center for Biotechnology Information

NPEV: 

Non-polio enterovirus

OPV: 

Oral poliovirus vaccine

RD: 

Human rhabdomyosarcoma

RT-PCR: 

Reverse transcription-PCR

t MRCA

The time of the most recent common ancestor

UCLD: 

Uncorrelated log-normal distributed

Declarations

Acknowledgements

We would like to thank the staff of EPI division in Shandong Center for Disease Control and Prevention for their contribution in data and sample collection. Additionally, we would like to thank the Environmental Systems Research Institute for our nonprofit application of ArcGIS software.

Funding

This study was supported by a grant from the National Natural Science Foundation of China (81302481) and a grant from Taishan Scholar Program of Shandong Province (ts.201511105).

Availability of data and materials

The datasets supporting the results of this article are included within the article and its Additional files 1, 2, 3 and 4. Newly generated VP1 sequences of Shandong E-14 isolates have been submitted to GenBank, under the following accession numbers: KX840391-KX840454. All sequences data that support the findings of this study have been deposited in TreeBASE and are accessible through the accession number 20306 (http://purl.org/phylo/treebase/phylows/study/TB2:S20306).

Authors’ contributions

PC, YL (Yan Li) and AX conceived the study and drafted the paper, ZT, HW, XL and YL (Yao Liu) gathered and analyzed the data, and SW, NZ and PW participated in RNA extraction. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

All protocols for this study were reviewed and approved by Ethics Review Committee of the Shandong Center for Disease Control and Prevention (No. 2016043), and the study was carried out in accordance with the principles of the Declaration of Helsinki. Written informed consents for the use of their clinical samples were obtained from all subjects (the legal guardians of the patients and contacts).

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Department of Epidemiology, School of Public Health, Shandong University, No. 44 Wenhuaxi Road, Jinan, 250012, People’s Republic of China
(2)
Shandong Provincial Key Laboratory of Infectious Disease Control and Prevention; Shandong Center for Disease Control and Prevention, No. 16992 Jingshi Road, Jinan, 250014, People’s Republic of China

References

  1. Enteroviruses: polioviruses, coxsackieviruses, echoviruses, and newer enteroviruses. Fields Virology. 5th ed. Philadelphia: Lippincott Williams & Wilkins; 2007.Google Scholar
  2. Khetsuriani N, Lamonte-Fowlkes A, Oberst S, Pallansch MA, Centers for Disease Control and Prevention. Enterovirus surveillance-United States, 1970–2005. MMWR Surveill Summ. 2006;55:1–20.PubMedGoogle Scholar
  3. Gao C, Ding Y, Zhou P, Feng J, Qian B, Lin Z, et al. Serological detection and analysis of anti-VP1 responses against various enteroviruses (EV) (EV-A, EV-B and EV-C) in Chinese individuals. Sci Rep. 2016;6:21979.View ArticlePubMedPubMed CentralGoogle Scholar
  4. Melnick JL. Problems associated with viral identification and classification in 1956. Ann N Y Acad Sci. 1957;67:363–82.View ArticlePubMedGoogle Scholar
  5. Ho M, Chen ER, Hsu KH, Twu SJ, Chen KT, Tsai SF, et al. An epidemic of enterovirus 71 infection in Taiwan. Taiwan Enterovirus Epidemic Working Group. N Engl J Med. 1999;341:929-35.Google Scholar
  6. Zhao K, Han X, Wang G, Hu W, Zhang W, Yu XF. Circulating coxsackievirus A16 identified as recombinant type A human enterovirus. China Emerg Infect Dis. 2011;17:1537–40.PubMedGoogle Scholar
  7. Chang J, Li J, Liu X, Liu G, Yang J, Wei W, et al. Broad protection with an inactivated vaccine against primary-isolated lethal enterovirus 71 infection in newborn mice. BMC Microbiol. 2015;15:139.View ArticlePubMedPubMed CentralGoogle Scholar
  8. Grosso E, Bergamini F. Epidemic infantile gastrointestinal syndrome caused by enteroviruses (ECHO-14, ECHO-9 and poliovirus of type 3). Boll IstSieroter Milan. 1960;39:495–509.Google Scholar
  9. Lepine P, Samaille J, Maurin J, Dubois O, Carre MC. Isolation of the ECHO 14 virus in the course of a nursery epidemic of gastroenteritis. Ann Inst Pasteur (Paris). 1960;99:161–6.Google Scholar
  10. De SanctisMonaldi T, Benedetto A, Formato FG. Epidemic of “acute enteric disease” associated with infection by “Echo 14” virus. Helv Paediatr Acta. 1966;21:239–49.Google Scholar
  11. Hinuma Y, Murai Y, Nakao T. Two outbreaks of echovirus 14 infection:a possible interference with oral poliovirus vaccine and a probable association with aseptic meningitis. J Hyg (Lond). 1965;63:277–84.View ArticleGoogle Scholar
  12. Galko NV, Vashukova SS. Antigenic variant of ECHO 14, the causative agent of acute serous meningitis. Vopr Virusol. 1980;3:345–9.Google Scholar
  13. Wako H, Yoshida S, Kawana R, Kaneko M. Case of bilateral ascending paralysis caused by ECHO-14 virus. Iqaku To Seibutsuqaku. 1965;71:265–9.Google Scholar
  14. Hughes JR, Wilfert CM, Moore M, Benirschke K, de Hoyos-Guevara E. Echovirus 14 infection associated with fatal neonatal hepatic necrosis. Am J Dis Child. 1972;123:61–7.PubMedGoogle Scholar
  15. Nakao T, Miura R, Takagi M, Inagaki T, Abe S, Oyanagi K, et al. The prevalence of aseptic meningitis caused by Echo-14 virus. Nihon Shonika Gakkai Zasshi. 1964;68:963–6.PubMedGoogle Scholar
  16. Krasota A, Loginovskih N, Ivanova O, Lipskaya G. Direct identification of Enteroviruses in cerebrospinal fluid of patients with suspected meningitis by nested PCR amplification. Viruses. 2016;8:1–10.View ArticleGoogle Scholar
  17. Dumaidi K, Frantzidou F, Papa A, Diza E, Antoniadis A. Enterovirus meningitis in Greece from 2003-2005: diagnosis, CSF laboratory findings, and clinical manifestations. J Clin Lab Anal. 2006;20:177–83.View ArticlePubMedGoogle Scholar
  18. Rao CD, Yergolkar P, Shankarappa KS. Antigenic diversity of enteroviruses associated with nonpolio acute flaccid paralysis, India, 2007-2009. Emerg Infect Dis. 2012;18:1833–40.View ArticlePubMedPubMed CentralGoogle Scholar
  19. Apostol LN, Suzuki A, Bautista A, Galang H, Paladin FJ, Fuji N, et al. Detection of non-polio enteroviruses from 17 years of virological surveillance of acute flaccid paralysis in the Philippines. J Med Virol. 2012;84:624–31.View ArticlePubMedPubMed CentralGoogle Scholar
  20. Tao Z, Wang H, Liu Y, Li Y, Jiang P, Liu G, et al. Non-Polio Enteroviruses from acute flaccid paralysis surveillance in Shandong Province, China, 1988–2013. Sci Rep. 2014;4:6167.View ArticlePubMedPubMed CentralGoogle Scholar
  21. Bingjun T, Yoshida H, Yan W, Lin L, Tsuji T, Shimizu H, et al. Molecular typing and epidemiology of non-polio enteroviruses isolated from Yunnan province, the People’s Republic of China. J Med Virol. 2008;80:670–9.View ArticlePubMedGoogle Scholar
  22. Khetsuriani N, Lamonte-Fowlkes A, Oberst S, Pallansch MA. Enterovirus surveillance-United States, 1970-2005[J]. MMWR Surveill Summ. 2006;55:1–20.PubMedGoogle Scholar
  23. Centers for Disease Control and Prevention. Nonpolio enterovirus and human parechovirus surveillance-United States, 2006–2008[J]. MMWR Morb Mortal Wkly Rep. 2010;59:1577–80.Google Scholar
  24. Ortner B, Huang CW, Schmid D, Mutz I, Wewalka G, Allerberger F, et al. Epidemiology of enterovirus types causing neurological disease in Austria 1999-2007: detection of clusters of echovirus 30 and enterovirus 71 and analysis of prevalent genotypes. J Med Virol. 2009;81:317–24.View ArticlePubMedGoogle Scholar
  25. Oyero OG, Adu FD, Ayukekbong JA. Molecular characterization of diverse species enterovirus-B types from children with acute flaccid paralysis and asymptomatic children in Nigeria. Virus Res. 2014;189:189–93.View ArticlePubMedGoogle Scholar
  26. Tan CY, Ninove L, Gaudart J, Nougairede A, Zandotti C, Thirion-Perrier L, et al. A retrospective overview of enterovirus infection diagnosis and molecular epidemiology in the public hospitals of Marseille, France (1985-2005). PLoS One. 2011;18:e18022.View ArticleGoogle Scholar
  27. Kim H, Kang B, Hwang S, Lee SW, Cheon DS, Kim K, et al. Clinical and enterovirus findings associated with acute flaccid paralysis in the Republic of Korea during the recent decade. J Med Virol. 2014;86:1584–9.View ArticlePubMedGoogle Scholar
  28. Oyero OG, Adu FD. Non-polio enteroviruses serotypes circulating in Nigeria. Afr J Med Med Sci. 2010;39(Suppl):201–8.PubMedGoogle Scholar
  29. Tao Z, Wang H, Li Y, Xu A, Zhang Y, Song L, et al. Cocirculation of two transmission lineages of echovirus 6 in jinan, china, as revealed by environmental surveillance and sequence analysis. Appl Environ Microbiol. 2011;77:3786–92.View ArticlePubMedPubMed CentralGoogle Scholar
  30. Bailly JL, Mirand A, Henquell C, Archimbaud C, Chambon M, Regagnon C, et al. Repeated genomic transfers from echovirus 30 to echovirus 6 lineages indicate co-divergence between co-circulating populations of the two human enterovirus serotypes. Infect Genet Evol. 2011;11:276–89.View ArticlePubMedGoogle Scholar
  31. Kyriakopoulou Z, Bletsa M, Tsakogiannis D, Dimitriou TG, Amoutzias GD, Gartzonika C, et al. Molecular epidemiology and evolutionary dynamics of Echovirus 3 serotype. Infect Genet Evol. 2015;32:305–12.View ArticlePubMedGoogle Scholar
  32. McWilliam Leitch EC, Cabrerizo M, Cardosa J, Harvala H, Ivanova OE, Kroes AC, et al. Evolutionary dynamics and temporal/geographical correlates of recombination in the human enterovirus echovirus types 9, 11, and 30. J Virol. 2010;84:9292–9300.View ArticlePubMedPubMed CentralGoogle Scholar
  33. McWilliam Leitch EC, Cabrerizo M, Cardosa J, Harvala H, Ivanova OE, Koike S, et al. The association of recombination events in the founding and emergence of subgenogroup evolutionary lineages of human enterovirus 71. J Virol. 2012;86:2676–85.View ArticlePubMedPubMed CentralGoogle Scholar
  34. Abdelkhalek I, Seghier M, Yahia AB, Touzi H, Meddeb Z, Triki H, et al. Molecular epidemiology of coxsackievirus type B1. Arch Virol. 2015;160:2815–21.View ArticlePubMedGoogle Scholar
  35. Alfonso-Morales A, Martínez-Pérez O, Dolz R, Valle R, Perera CL, Bertran K, et al. Spatiotemporal phylogenetic analysis and molecular characterization of infectious bursal disease viruses based on the VP2 hyper-variable region. PLoS One. 2013;8:e65999.View ArticlePubMedPubMed CentralGoogle Scholar
  36. World Health Organization. Isolation and identification of polioviruses. WHO Polio Laboratory Manual. 4th ed. Geneva: World Health Organization; 2004.Google Scholar
  37. Oberste MS, Maher K, Kilpatrick DR, Pallansch MA. Molecular evolution of the human enteroviruses: correlation of serotype with VP1 sequence and application to picornavirus classification. J Virol. 1999;73:1941–8.PubMedPubMed CentralGoogle Scholar
  38. Oberste MS, Maher K, Flemister MR, Marchetti G, Kilpatrick DR, Pallansch MA. Comparison of classic and molecular approaches for the identification of untypeableenteroviruses. J Clin Microbiol. 2000;38:1170–4.PubMedPubMed CentralGoogle Scholar
  39. Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.Google Scholar
  40. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6:Molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.View ArticlePubMedPubMed CentralGoogle Scholar
  41. Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111–20.View ArticlePubMedGoogle Scholar
  42. Schierup MH, Hein J. Consequences of recombination on traditional phylogenetic analysis. Genetics. 2000;156:879–91.PubMedPubMed CentralGoogle Scholar
  43. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.View ArticlePubMedPubMed CentralGoogle Scholar
  44. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772.View ArticlePubMedPubMed CentralGoogle Scholar
  45. Raftery A, Newton M, Satagopan J, Krivitsky P. Estimating the integrated likelihood via posterior simulation using the harmonic mean identity. In: Bernardo JM, Bayarri MJ, Berger JO, editors. Bayesian Statistics. Oxford: Oxford University Press; 2007. p. 1–45.Google Scholar
  46. Newton MA, Raftery AE. Approximating Bayesian inference with the weighted likelihood bootstrap. J R Stat Soc B Methodol. 1994;56:3–48.Google Scholar
  47. Bielejec F, Rambaut A, Suchard MA, Lemey P. SPREAD: spatial phylogenetic reconstruction of evolutionary dynamics. Bioinformatics. 2011;27:2910–2.View ArticlePubMedPubMed CentralGoogle Scholar

Copyright

© The Author(s). 2017

Advertisement