BACKGROUND
coxsackievirus  b <dig> belongs to the species enterovirus b  in the genus enterovirus, family picornaviridae.
long-term surveillance of circulating ev serotypes has identified three main outbreak patterns: epidemic, endemic and sporadic  <cit> .
however, most studies of ev have focused on ‘epidemic’ pathogens rather than on ‘sporadic’ or ‘endemic’ pathogens such as cv-b2; thus, no phylogenetic analysis of cv-b <dig> has been reported.
in terms of epidemiology, a virus is considered ‘endemic’ if its population remains constant within a population in a geographic area; it is considered ‘sporadic’ if it appears often but not regularly.
the circulation model of individual serotypes differs by season, geographic location, herd immunity and viral genotypes  <cit> .
thus, analysis history of spatiotemporal transmission and viral population dynamics can reveal trends in their patterns of circulation and succession.
in taiwan, episodes of unusually high cv-b <dig> activity were detected in  <dig>   <dig>   <dig> and  <dig>  and infections are still reported almost annually to the taiwan centers for disease control .
cv-b <dig> has an endemic circulation pattern worldwide where it comprises  <dig>  %– <dig>  % of the annual reported ev rates  .
the peak activity of cv-b <dig>  is rare, but it is commonly identified as one of the top  <dig> most frequently detected ev serotypes globally  <cit> .
cv-b <dig> shares many common properties and clinical features with other cv-bs but they have different outbreak tempos, and overwhelming infection may develop in infants   <cit> .
cv-b infection is also associated with myocarditis, polymyositis and chronic autoimmune diseases such as insulin-dependent diabetes mellitus  <cit> .
moreover, outbreaks of emerging and re-emerging infectious diseases are strongly associated with climate change and human activity  <cit> .
increased cv-b <dig> activity levels and severe infections of young infants have been reported in the usa  and south korea   <cit> .
therefore, identifying the phylodynamics of cv-b <dig> is essential for bridging relationships between successive prevalent serotypes and identifying ev circulation trends.
the bayesian evolutionary analysis sampling trees  program combined various models to infer complex relationships among viral phylodynamics  <cit> .
in addition to genetic variation, viral spatial and temporal components are fundamental issues that require clarification as the dispersal pattern expands.
a clear understanding of the tree topology, transmission history, and demographic dynamics of an endemic pathogen such as cv-b <dig> during viral succession and circulation is needed to identify its prevalence trends.
in evs, the vp <dig> gene sequence encodes the major serotypic epitope, thereby making it the target choice for molecular typing  <cit> .
recombination, which is reportedly a common phenomenon in the ev family, can be identified based on phylogenetic incongruencies in the vp <dig> and 3dp°l gene regions  because 3dpol is located in the last codon region of the ev genome.
to explore the phylogenetic structure of endemic pathogens such as cv-b <dig>  both vp <dig> and 3dpol regions were ananlyzed by the beast program to elucidate its phylodynamic properties and reconstruct its spatiotemporal transmission history.
methods
ethics statement
twenty virus strains were randomly selected from cv-b2-positive viral stock.
these virus strains were collected during 1988– <dig> by the taiwan centers for disease control, taiwan  and two medical centers in southern taiwan, kaohsiung veterans general hospital and kaohsiung medical university hospital.
this study was approved by the ethics committees of both hospitals.
all samples were de-identified and analyzed anonymously.
informed consent was waived because we only conducted experiments on viral isolates obtained from clinically necessary laboratory procedure and no harm to the patients was anticipated because none of the patients’ medical history was studied.
specimens, viral rna extraction, reverse transcription-polymerase chain reaction  and sequencing
among the twenty cv-b <dig> strains, fifteen were isolated from throat swabs and five were isolated from rectal swabs or stool.
forty percent of the infected patients were infant younger than  <dig> year .
the male-female ratio was  <dig> : <dig>  viral rna was extracted using a qiamp viral rna purification kit according to the manufacturer instructions .
the rt-pcr and sequencing were performed as described previously  <cit> .
table  <dig> lists the degenerate primer pairs used for amplifying and sequencing the vp <dig> and partial 3dpol genes.
the purified pcr products were sequenced with an abi prism ready reaction dideoxy terminator cycle sequencing kit .
the obtained sequence data were submitted to genbank under accession numbers ab <dig>  ab <dig>  ab <dig>  ab862115–ab <dig>  lc055763–lc <dig>  lc055768–lc <dig>  and lc057297–lc <dig> table  <dig> primer sets used for amplification and sequencing of coxsackievirus b2
3dpol
3dpol
af: forward primer, r: reverse primer
bnumbering system used for the coxsackievirus b <dig> strain
cprimer designed by the centers for disease control, taiwan
model selection and sequence variation detection
multiple sequence alignments were performed with t-coffee  <cit> .
in addition to the  <dig> cv-b <dig> strains obtained in this study, the sequences available in genbank were also included.
sequences were excluded if they had a nonsense or frame-shift mutation pattern in a single strain.
the final vp <dig> dataset comprised  <dig> sequences and the partial 3dpol dataset comprised  <dig> sequences .
the jmodeltest v  <dig> .
<dig> program  <cit>  was used to select the best-fit model for both datasets according to akaike’s information criterion.
a 4-category generalised time reversible  <cit>  model with a gamma distribution  was used for vp <dig> and 3dpol .
for the beast estimation, eight model compositions were compared: two substitution models   <cit> ), two population models , and two relaxed molecular clock models   <cit> .
the bayes markov chain monte carlo  analysis results were sampled regularly until convergence was reached.
recombination was estimated using the recombination detection program  with the default settings  <cit> .
pairwise comparisons of nucleotide and amino acid  sequences were calculated based on p-distance matrices implemented in the molecular evolutionary genetics analysis version  <dig>  program.
the ratio of non-synonymous/synonymous substitutions rate  was used as an indicator of the selection force acting on the coding sequences and was detected by using single likelihood ancestor counting  method via the datamonkey website  <cit> .
phylodynamic analyses
phylogenetic analyses were performed as described previously  <cit> .
the bmcmc tree analyses were performed using beast v. <dig> .
<dig>  <cit> .
nodal reliability of bmcmc trees were estimated according to posterior probability .
the breakpoint for significant support was set to a pp value >  <dig> .
the beast program was used to co-estimate the nt substitution rate and population growth model in the tracer v. <dig>  program  <cit> .
parameters with effective sample size  values >  <dig> were considered reliable.
briefly, each run of vp <dig> and 3dp°l had a chain length of  <dig> million generations, and sampling was performed once every  <dig> generations.
when the major parameters with ess value > <dig> were considered as candidates for the best model composition, the best model for both the vp <dig> and 3dpol regions was srd06 + uced + con .
all estimation parameters were shown as mean and 95% highest posterior density .
demographic changes over time  were also explored by bsp method  <cit> .
the figtree v. <dig> .
<dig> program was used to construct and visualize the maximum clade credibility tree.
the spread program  <cit>  was used to estimate the route of virus transmission and to calculate the bayes factors  associated with major virus dispersal routes.
a major virus dispersal route was defined as a route in which at least two locations had bf values > <dig>  furthermore, the treestat program in the beast package was used to summarize statistical tests of neutrality  and to calculate three tree-balance statistical measures: cherry count , colless’s tree imbalance , and b <dig>  <cit> .
RESULTS
no recombination events were detected in the vp <dig> or 3dpol region.
high support values  were obtained for four genotypes in the vp <dig> tree  and for five genotypes in the 3dp°l tree  .
both the vp <dig> and 3dpol trees contained spatiotemporally structured clusters, but the vp <dig> and 3dpol trees were clearly incongruent.
based on a genetic discrimination rate of at least 15%, five clusters were obtained in the vp <dig> region .
in the ascent cluster , the prototype ohio strain  and taiwan strain  <dig>  were clustered together in vp <dig> and in 3dpol .
however, the support value was not statistically significant for vp <dig>  in both the vp <dig> and partial 3dpol regions, taiwan strain  <dig> had the same sequence as the prototype ohio strain.
strains gii , giv , and gv  were isolated in asia, where giv was geographically distributed in india.
by contrast, giii was isolated in france  and australia .
the analysis of the historical transmission routes for vp <dig> by using the spread program showed that only the taiwan-south korea route had a bf of  <dig> .
the demographic history determined by bsp showed that the cv-b <dig> viral population had maintained a stable level with only a slight decrease since  <dig>  where the median neτ was  <dig> – <dig>  .
compared with the data reported in year  <dig>  the neτ for cv-b <dig> was much lower than the median neτ values of 4 ×  <dig> reported for ev <dig>  <cit>  and  <dig> for cv-b <dig>  <cit> .fig.
<dig> maximum clade credibility  phylogeny of  <dig> vp <dig> sequences of coxsackievirus b <dig>  for each branch, the thickness indicates the state probability, and the color indicates the most probable location.
support values are indicated on the major nodes.
the genotypes and nucleotide/amino acid similarity within genotypes are shown on the right.
for each strain name, the vp <dig> genotypes are differentiated by colors , whereas the 3dpol genotypes are differentiated by shading .
the branch length is proportional to the evolutionary time, and the scale bar is proportional to calendar time
fig.
<dig> maximum clade credibility  phylogeny of  <dig> 3dpol sequences coxsackievirus b <dig>  for each branch, the thickness indicates the state probability and the color indicates the most probable location.
support values are indicated on the major nodes.
the genotypes and nucleotide/amino acid similarity within genotypes are shown on the right.
for each strain name, the vp <dig> genotypes are differentiated by colors , whereas the 3dpol genotypes are differentiated by shading .
the branch length is proportional to the evolutionary time, and the scale bar is proportional to calendar time
fig.
<dig> bayesian skyline plot of   <dig> vp <dig> sequences and   <dig> 3dpol sequences.
the x-axis is the time scale  and the y-axis is the logarithmic neτ scale .
the thick solid line indicates the median estimates and the shaded area indicates the 95% highest posterior density
phylogenetic comparisons of the vp <dig> topology with those of other evs revealed similar features:  a serotype-based monophyletic origin;  different cv-b <dig> genotypes evolved and diverged sequentially, where each genotype had a different geographic distribution and circulation half-life;  the prevalent lineages were chronologically replaced by emerging different lineages that co-existed at the same time and in the same location and  tree topologies with long internal branches  and intricate short terminal branches  <cit> .
in particular, the long internal branches and short terminal branches indicated a stable pathogen population size  <cit> , which supported bsp estimation of the cv-b <dig> viral population in this study.
five genotypes are depicted in the 3dpol tree .
the ancestral america strains  were clustered in ga–c, especially in gb.
by contrast, most of the taiwan strains were clustered in gd–ge with the other asian strains.
interestingly, gb in the 3dpol tree included two taiwan strains isolated during 2011– <dig> .
this indicated that those strains with gv sequence in vp <dig> region has endemic in taiwan since 2002– <dig>  but the 3dpol gene might have switched from ge to gb, because the 3dpol gene in the latter two isolates  clustered together with us strains in gb.
therefore, changes in the trends in cv-b <dig> in taiwan should be monitored continuously in taiwan.
in contrast to the low  median neτ value estimated by vp <dig> region of cv-b <dig> in this study, the demographic history analysis of the 3dpol region by bsp showed that cv-b <dig> emerged during  <dig>  where it had a larger neτ value of  <dig>  before it began to bottleneck in  <dig>  after peaking in  <dig> the neτ values decreased gradually from  <dig>  to  <dig>  during 1985– <dig> .
only capsid-based analyses reveal that serotypes of evs are monophyletic; however, when the analysis is based on regions outside the capsid region, evs are monophyletic species rather than serotypes  <cit> .
thus, the 3dpol-based estimation of neτ was relevant to ev-b species .
one proposed explanation for this phenomenon is independent evolution of vp <dig> and 3dpol  <cit> , which is consistent with the discordant topologies of vp <dig> and 3dpol determined in the present study.
our results indicate that recombination may occur within discordant strains, and that mutation in the capsid region is a result of immune escape.
however, this study detected high  aa sequence similarity.
based on the p-distances calculated in this study, the estimated maximum differences in the 846-nt  and 282-aa sequences of the vp <dig> gene were  <dig>  % and  <dig>  %, respectively.
the estimated maximum differences in the 441-nt  and 147-aa sequences of the partial 3dpol gene were  <dig>  % and  <dig>  %, respectively.
further, the estimated evolutionary rates of vp <dig> and 3dpol were  <dig>  × 10− <dig>  and  <dig>  × 10− <dig>  substitutions/site/year , respectively.
the vp <dig> mutation rate of cv-b <dig> is within the mutation rate range reported in a previous study of evs   <cit> .
in the present study, the estimated evolution rate of 3dpol  was not only higher than that of vp <dig> but also higher than the estimated evolution rate reported for 3dpol in previous studies of evs   <cit> .
the dn/ds values for the vp <dig> and 3dpol region were  <dig>  and  <dig> , respectively, but neither region contained a positive selection site.
in summary, the low and slightly decreasing neτ value, low evolution rate, and negative selection indicate that the variation in this vp <dig> region was fixed by adaptive fitness.
phylodynamic analysis depicts interaction among epidemics, evolution and selection and reveals the demographic and spatiotemporal transmission signatures  <cit> .
a phylogenetic tree with an imbalanced trunk and short terminal branches suggests acute infection with partial cross-immunity, for example, the phylogenetic tree topology for seasonal influenza virus is ladder-like.
in contrast, the topology of a chronic infection such as hepatitis c virus is star-like with long terminal branches  <cit> .
the tree balance is impacted by selection; therefore, a ladder-like backbone indicates the trajectory of fitness.
table  <dig> summarizes the results of the statistical analyses of tree shape and balance obtained in the present study.
fu and li’s d value is a classic summary statistic for testing neutrality, which indicates the normalize difference between the external branch lengths and total tree length.
as the symmetry of a phylogeny increases, the values for b <dig> and cn are expected to increase, whereas the value of ic is expected to decrease.
the distributions of the mean ic indices among human and zoonotic rna viruses reportedly peak between  <dig>  and  <dig>   <cit> .
the values for the vp <dig> tree were within this range in the present study, but 3dpol had higher tree balance indices compared with vp <dig> table  <dig> summary statistics used in tests of demographic neutrality and tree balance
d
b
b
1
c
cn
d
ic
d
d
b
1
cn
ic
ahpd highest posterior density
b
d fu and li’s d: normalized difference between the external branch lengths and total tree length.
this is a classic summary statistic used for testing neutrality
c
b
<dig> maximum number of nodes between an internal node and the tips of the tree.
a larger value indicates a more balanced tree
d
cn cherry count, the number of internal nodes with exactly two terminal branches.
a larger number indicates a more balanced tree
e
ic colless’s tree imbalance index with a range of  <cit> .
a larger number indicates a more imbalanced tree
a major challenge in modern phylodynamics is qualitatively and quantitatively describing a tree topology to elucidate the epidemiological, evolutionary, and demographic characteristics of a pathogen.
to investigate the evolutionary history of an endemic or sporadic pathogen, a bmcmc sampling framework was used to infer time scaled phylogenies of cv-b <dig> during 1947– <dig>  sampling bias was unavoidable in this study because the incomplete genbank data precluded sampling of all strains isolated in earlier time periods .
an even greater problem was that the sequences tended to accumulate in the same isolation years and locations because cv-b <dig> outbreaks are rare.
on average, cvb <dig> has been implicated in  <dig>  %– <dig>  % of the reported surveillance results for known ev serotypes, and thus the low frequency of cv-b <dig> outbreaks may be related to the stable low neτ values estimates obtained in this study.
furthermore, low herd immunity selection may have an important role in endemic evs because the vp <dig> has high aa sequence similarity and a low rate of evolution.
compare to 3dpol region, the vp <dig> region has a more balanced topology.
CONCLUSIONS
an endemic or sporadic ev, cv-b <dig>  was analyzed in a phylodynamic analysis, and four genotypes in the vp <dig> region and five genotypes in the 3dpol region were determined.
recombination events were not detected in the vp <dig> sequences or 3dpol sequences, but vp <dig> and 3dpol trees were clearly incongruent, thereby indicating the occurrence of recombination events involving the region between vp <dig> and 3dpol.
although the random branching topology observed in this study is typically interpreted as the rapid evolution of an rna virus, the imbalanced tree topology and the high similarity of the vp <dig> sequences indicate that endemic cv-b <dig> is characterized by a low viral population and low herd immune selection.
availability of supporting data
see additional files for supporting data.
additional files
additional file 1:
annual reported rate of coxsackievirus b <dig>  in taiwan.
additional file 2:
global surveillance of coxsackievirus b <dig>
additional file 3:
list of sampled coxsackievirus b <dig>  strains.
additional file 4:
models compared by an extension of air carrier investment model .
abbreviations
aaamino acid
bfbayes factor
beastthe bayesian evolutionary analysis sampling trees
bmcmcbayesian monte carlo markov chain
bspbayesian skyline plot
cneffective sample size
cvcoxsackievirus
eventerovirus
esseffective sample size
gtrgeneralised time reversible
hpdhighest probability density
iccolless’s tree imbalance
mccmaximum clade credibility
neτeffective population size over time
ntnucleotide
ppposterior probability
slacsingle likelihood ancestor counting
srd06shapiro-rambaut-drummond-2006
tmrcatime to most recent common ancestor
uceduncorrelated exponential distribution
hui-wen huang and yao-shen chen contributed equally to this work.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
conception and design of experiments: pyc, khl.
performance of experiments: bcc, ych, lcc, cfw, hjs.
analysis of data: pyc, ysc, jyfc, pll, ycl, hlc.
contribution of reagents/materials/analysis tools: khl, hwh, yys, yct, tsh.
preparation of manuscript: pyc, hwh, bsd.
all authors have read and approved the manuscript.

 coxsackievirus b belongs to the species enterovirus b in the genus enterovirus, family picornaviridae. long-term surveillance of circulating ev serotypes has identified three main outbreak patterns: epidemic, endemic and sporadic. however, most studies of ev have focused on ‘epidemic’ pathogens rather than on ‘sporadic’ or ‘endemic’ pathogens such as cv-b2; thus, no phylogenetic analysis of cv-b has been reported. in terms of epidemiology, a virus is considered ‘endemic’ if its population remains constant within a population in a geographic area; it is considered ‘sporadic’ if it appears often but not regularly. the circulation model of individual serotypes differs by season, geographic location, herd immunity and viral genotypes. thus, analysis history of spatiotemporal transmission and viral population dynamics can reveal trends in their patterns of circulation and succession. in taiwan, episodes of unusually high cv-b activity were detected in and and infections are still reported almost annually to the taiwan centers for disease control. cv-b has an endemic circulation pattern worldwide where it comprises %– % of the annual reported ev rates. the peak activity of cv-b is rare, but it is commonly identified as one of the top most frequently detected ev serotypes globally. cv-b shares many common properties and clinical features with other cv-bs but they have different outbreak tempos, and overwhelming infection may develop in infants. cv-b infection is also associated with myocarditis, polymyositis and chronic autoimmune diseases such as insulin-dependent diabetes mellitus. moreover, outbreaks of emerging and re-emerging infectious diseases are strongly associated with climate change and human activity. increased cv-b activity levels and severe infections of young infants have been reported in the usa and south korea. therefore, identifying the phylodynamics of cv-b is essential for bridging relationships between successive prevalent serotypes and identifying ev circulation trends. the bayesian evolutionary analysis sampling trees program combined various models to infer complex relationships among viral phylodynamics. in addition to genetic variation, viral spatial and temporal components are fundamental issues that require clarification as the dispersal pattern expands. a clear understanding of the tree topology, transmission history, and demographic dynamics of an endemic pathogen such as cv-b during viral succession and circulation is needed to identify its prevalence trends. in evs, the vp gene sequence encodes the major serotypic epitope, thereby making it the target choice for molecular typing. recombination, which is reportedly a common phenomenon in the ev family, can be identified based on phylogenetic incongruencies in the vp and 3dp°l gene regions because 3dpol is located in the last codon region of the ev genome. to explore the phylogenetic structure of endemic pathogens such as cv-b both vp and 3dpol regions were ananlyzed by the beast program to elucidate its phylodynamic properties and reconstruct its spatiotemporal transmission history. methods ethics statement twenty virus strains were randomly selected from cv-b2-positive viral stock. these virus strains were collected during 1988– by the taiwan centers for disease control, taiwan and two medical centers in southern taiwan, kaohsiung veterans general hospital and kaohsiung medical university hospital. this study was approved by the ethics committees of both hospitals. all samples were de-identified and analyzed anonymously. informed consent was waived because we only conducted experiments on viral isolates obtained from clinically necessary laboratory procedure and no harm to the patients was anticipated because none of the patients’ medical history was studied. specimens, viral rna extraction, reverse transcription-polymerase chain reaction and sequencing among the twenty cv-b strains, fifteen were isolated from throat swabs and five were isolated from rectal swabs or stool. forty percent of the infected patients were infant younger than year. the male-female ratio was: viral rna was extracted using a qiamp viral rna purification kit according to the manufacturer instructions. the rt-pcr and sequencing were performed as described previously. table  lists the degenerate primer pairs used for amplifying and sequencing the vp and partial 3dpol genes. the purified pcr products were sequenced with an abi prism ready reaction dideoxy terminator cycle sequencing kit. the obtained sequence data were submitted to genbank under accession numbers ab ab ab ab862115–ab lc055763–lc lc055768–lc and lc057297–lc table primer sets used for amplification and sequencing of coxsackievirus b2 3dpol 3dpol af: forward primer, r: reverse primer bnumbering system used for the coxsackievirus b strain cprimer designed by the centers for disease control, taiwan model selection and sequence variation detection multiple sequence alignments were performed with t-coffee. in addition to the cv-b strains obtained in this study, the sequences available in genbank were also included. sequences were excluded if they had a nonsense or frame-shift mutation pattern in a single strain. the final vp dataset comprised sequences and the partial 3dpol dataset comprised sequences. the jmodeltest v. program was used to select the best-fit model for both datasets according to akaike’s information criterion. a 4-category generalised time reversible model with a gamma distribution was used for vp and 3dpol. for the beast estimation, eight model compositions were compared: two substitution models ), two population models, and two relaxed molecular clock models. the bayes markov chain monte carlo analysis results were sampled regularly until convergence was reached. recombination was estimated using the recombination detection program with the default settings. pairwise comparisons of nucleotide and amino acid sequences were calculated based on p-distance matrices implemented in the molecular evolutionary genetics analysis version program. the ratio of non-synonymous/synonymous substitutions rate was used as an indicator of the selection force acting on the coding sequences and was detected by using single likelihood ancestor counting method via the datamonkey website. phylodynamic analyses phylogenetic analyses were performed as described previously. the bmcmc tree analyses were performed using beast v... nodal reliability of bmcmc trees were estimated according to posterior probability. the breakpoint for significant support was set to a pp value > . the beast program was used to co-estimate the nt substitution rate and population growth model in the tracer v. program. parameters with effective sample size values >  were considered reliable. briefly, each run of vp and 3dp°l had a chain length of million generations, and sampling was performed once every generations. when the major parameters with ess value > were considered as candidates for the best model composition, the best model for both the vp and 3dpol regions was srd06 + uced + con. all estimation parameters were shown as mean and 95% highest posterior density. demographic changes over time were also explored by bsp method. the figtree v.. program was used to construct and visualize the maximum clade credibility tree. the spread program was used to estimate the route of virus transmission and to calculate the bayes factors associated with major virus dispersal routes. a major virus dispersal route was defined as a route in which at least two locations had bf values > furthermore, the treestat program in the beast package was used to summarize statistical tests of neutrality and to calculate three tree-balance statistical measures: cherry count, colless’s tree imbalance, and b.  no recombination events were detected in the vp or 3dpol region. high support values were obtained for four genotypes in the vp tree and for five genotypes in the 3dp°l tree. both the vp and 3dpol trees contained spatiotemporally structured clusters, but the vp and 3dpol trees were clearly incongruent. based on a genetic discrimination rate of at least 15%, five clusters were obtained in the vp region. in the ascent cluster, the prototype ohio strain and taiwan strain were clustered together in vp and in 3dpol. however, the support value was not statistically significant for vp in both the vp and partial 3dpol regions, taiwan strain had the same sequence as the prototype ohio strain. strains gii, giv, and gv were isolated in asia, where giv was geographically distributed in india. by contrast, giii was isolated in france and australia. the analysis of the historical transmission routes for vp by using the spread program showed that only the taiwan-south korea route had a bf of. the demographic history determined by bsp showed that the cv-b viral population had maintained a stable level with only a slight decrease since where the median neτ was –. compared with the data reported in year the neτ for cv-b was much lower than the median neτ values of 4 ×  reported for ev and for cv-b.fig. maximum clade credibility phylogeny of vp sequences of coxsackievirus b for each branch, the thickness indicates the state probability, and the color indicates the most probable location. support values are indicated on the major nodes. the genotypes and nucleotide/amino acid similarity within genotypes are shown on the right. for each strain name, the vp genotypes are differentiated by colors, whereas the 3dpol genotypes are differentiated by shading. the branch length is proportional to the evolutionary time, and the scale bar is proportional to calendar time fig. maximum clade credibility phylogeny of 3dpol sequences coxsackievirus b for each branch, the thickness indicates the state probability and the color indicates the most probable location. support values are indicated on the major nodes. the genotypes and nucleotide/amino acid similarity within genotypes are shown on the right. for each strain name, the vp genotypes are differentiated by colors, whereas the 3dpol genotypes are differentiated by shading. the branch length is proportional to the evolutionary time, and the scale bar is proportional to calendar time fig. bayesian skyline plot of vp sequences and 3dpol sequences. the x-axis is the time scale and the y-axis is the logarithmic neτ scale. the thick solid line indicates the median estimates and the shaded area indicates the 95% highest posterior density phylogenetic comparisons of the vp topology with those of other evs revealed similar features: a serotype-based monophyletic origin; different cv-b genotypes evolved and diverged sequentially, where each genotype had a different geographic distribution and circulation half-life; the prevalent lineages were chronologically replaced by emerging different lineages that co-existed at the same time and in the same location and tree topologies with long internal branches and intricate short terminal branches. in particular, the long internal branches and short terminal branches indicated a stable pathogen population size, which supported bsp estimation of the cv-b viral population in this study. five genotypes are depicted in the 3dpol tree. the ancestral america strains were clustered in ga–c, especially in gb. by contrast, most of the taiwan strains were clustered in gd–ge with the other asian strains. interestingly, gb in the 3dpol tree included two taiwan strains isolated during 2011–. this indicated that those strains with gv sequence in vp region has endemic in taiwan since 2002– but the 3dpol gene might have switched from ge to gb, because the 3dpol gene in the latter two isolates clustered together with us strains in gb. therefore, changes in the trends in cv-b in taiwan should be monitored continuously in taiwan. in contrast to the low median neτ value estimated by vp region of cv-b in this study, the demographic history analysis of the 3dpol region by bsp showed that cv-b emerged during where it had a larger neτ value of before it began to bottleneck in after peaking in the neτ values decreased gradually from to during 1985–. only capsid-based analyses reveal that serotypes of evs are monophyletic; however, when the analysis is based on regions outside the capsid region, evs are monophyletic species rather than serotypes. thus, the 3dpol-based estimation of neτ was relevant to ev-b species. one proposed explanation for this phenomenon is independent evolution of vp and 3dpol, which is consistent with the discordant topologies of vp and 3dpol determined in the present study. our results indicate that recombination may occur within discordant strains, and that mutation in the capsid region is a result of immune escape. however, this study detected high aa sequence similarity. based on the p-distances calculated in this study, the estimated maximum differences in the 846-nt and 282-aa sequences of the vp gene were % and %, respectively. the estimated maximum differences in the 441-nt and 147-aa sequences of the partial 3dpol gene were % and %, respectively. further, the estimated evolutionary rates of vp and 3dpol were  × 10− and  × 10− substitutions/site/year, respectively. the vp mutation rate of cv-b is within the mutation rate range reported in a previous study of evs. in the present study, the estimated evolution rate of 3dpol was not only higher than that of vp but also higher than the estimated evolution rate reported for 3dpol in previous studies of evs. the dn/ds values for the vp and 3dpol region were and, respectively, but neither region contained a positive selection site. in summary, the low and slightly decreasing neτ value, low evolution rate, and negative selection indicate that the variation in this vp region was fixed by adaptive fitness. phylodynamic analysis depicts interaction among epidemics, evolution and selection and reveals the demographic and spatiotemporal transmission signatures. a phylogenetic tree with an imbalanced trunk and short terminal branches suggests acute infection with partial cross-immunity, for example, the phylogenetic tree topology for seasonal influenza virus is ladder-like. in contrast, the topology of a chronic infection such as hepatitis c virus is star-like with long terminal branches. the tree balance is impacted by selection; therefore, a ladder-like backbone indicates the trajectory of fitness. table  summarizes the results of the statistical analyses of tree shape and balance obtained in the present study. fu and li’s d value is a classic summary statistic for testing neutrality, which indicates the normalize difference between the external branch lengths and total tree length. as the symmetry of a phylogeny increases, the values for b and cn are expected to increase, whereas the value of ic is expected to decrease. the distributions of the mean ic indices among human and zoonotic rna viruses reportedly peak between and. the values for the vp tree were within this range in the present study, but 3dpol had higher tree balance indices compared with vp table summary statistics used in tests of demographic neutrality and tree balance d b b 1 c cn d ic d d b 1 cn ic ahpd highest posterior density b d fu and li’s d: normalized difference between the external branch lengths and total tree length. this is a classic summary statistic used for testing neutrality c b maximum number of nodes between an internal node and the tips of the tree. a larger value indicates a more balanced tree d cn cherry count, the number of internal nodes with exactly two terminal branches. a larger number indicates a more balanced tree e ic colless’s tree imbalance index with a range of. a larger number indicates a more imbalanced tree a major challenge in modern phylodynamics is qualitatively and quantitatively describing a tree topology to elucidate the epidemiological, evolutionary, and demographic characteristics of a pathogen. to investigate the evolutionary history of an endemic or sporadic pathogen, a bmcmc sampling framework was used to infer time scaled phylogenies of cv-b during 1947– sampling bias was unavoidable in this study because the incomplete genbank data precluded sampling of all strains isolated in earlier time periods. an even greater problem was that the sequences tended to accumulate in the same isolation years and locations because cv-b outbreaks are rare. on average, cvb has been implicated in %– % of the reported surveillance results for known ev serotypes, and thus the low frequency of cv-b outbreaks may be related to the stable low neτ values estimates obtained in this study. furthermore, low herd immunity selection may have an important role in endemic evs because the vp has high aa sequence similarity and a low rate of evolution. compare to 3dpol region, the vp region has a more balanced topology.  an endemic or sporadic ev, cv-b was analyzed in a phylodynamic analysis, and four genotypes in the vp region and five genotypes in the 3dpol region were determined. recombination events were not detected in the vp sequences or 3dpol sequences, but vp and 3dpol trees were clearly incongruent, thereby indicating the occurrence of recombination events involving the region between vp and 3dpol. although the random branching topology observed in this study is typically interpreted as the rapid evolution of an rna virus, the imbalanced tree topology and the high similarity of the vp sequences indicate that endemic cv-b is characterized by a low viral population and low herd immune selection. availability of supporting data see additional files for supporting data. additional files additional file 1: annual reported rate of coxsackievirus b in taiwan. additional file 2: global surveillance of coxsackievirus b additional file 3: list of sampled coxsackievirus b strains. additional file 4: models compared by an extension of air carrier investment model. abbreviations aaamino acid bfbayes factor beastthe bayesian evolutionary analysis sampling trees bmcmcbayesian monte carlo markov chain bspbayesian skyline plot cneffective sample size cvcoxsackievirus eventerovirus esseffective sample size gtrgeneralised time reversible hpdhighest probability density iccolless’s tree imbalance mccmaximum clade credibility neτeffective population size over time ntnucleotide ppposterior probability slacsingle likelihood ancestor counting srd06shapiro-rambaut-drummond-2006 tmrcatime to most recent common ancestor uceduncorrelated exponential distribution hui-wen huang and yao-shen chen contributed equally to this work. competing interests the authors declare that they have no competing interests. authors’ contributions conception and design of experiments: pyc, khl. performance of experiments: bcc, ych, lcc, cfw, hjs. analysis of data: pyc, ysc, jyfc, pll, ycl, hlc. contribution of reagents/materials/analysis tools: khl, hwh, yys, yct, tsh. preparation of manuscript: pyc, hwh, bsd. all authors have read and approved the manuscript. 
in enteroviruses, a clear understanding of the characteristics of tree topology, transmission, and its demographic dynamics in viral succession and circulation are essential for identifying prevalence trends in endemic pathogens such as coxsackievirus b. the balance indices, high similarity, and low evolutionary rate in the vp region indicated mild herd immunity selection in the major capsid region. phylodynamic analysis can reveal demographic trends and herd immunity in endemic pathogens. studies regarding coxsackievirus tend to focus on epidemic outbreaks, an imbalanced topology is considered to be an indication of acute infection with partial cross-immunity. some strains with tree discordance indicated the occurrence of recombination in the region between the vp and 3dpol genes. the similarities of vp and 3dpol gene were %– % and %– %, respectively.
 studies regarding coxsackievirus tend to focus on epidemic outbreaks, an imbalanced topology is considered to be an indication of acute infection with partial cross-immunity. in enteroviruses, a clear understanding of the characteristics of tree topology, transmission, and its demographic dynamics in viral succession and circulation are essential for identifying prevalence trends in endemic pathogens such as coxsackievirus b. this study applied a novel bayesian evolutionary approach to elucidate the phylodynamic characteristics of cv-b a dataset containing vp sequences and a dataset containing partial 3dpol sequencing were analyzed, where each dataset included taiwan sequences isolated during 1988–  four and five genotypes were determined based on the 846-nucleotide vp and 441-nucleotide 3dpol regions, respectively, with spatiotemporally structured topologies in both trees. some strains with tree discordance indicated the occurrence of recombination in the region between the vp and 3dpol genes. the similarities of vp and 3dpol gene were %– % and %– %, respectively. analyses of population dynamics using vp dataset indicated that the endemic cv-b has a small effective population size. the balance indices, high similarity, and low evolutionary rate in the vp region indicated mild herd immunity selection in the major capsid region.  phylodynamic analysis can reveal demographic trends and herd immunity in endemic pathogens. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords coxsackievirus b2phylodynamicspatiotemporal transmissiondemographic historymolecular epidemiologyissue-copyright-statement© the author 2015 