Despite our study focused on the available sequences of HPeV, more accurate estimates could probably be obtained with broader and more homogenous sampling timescale, preferably for all types. Yet, this may be a daunting task because it is difficult to obtain older samples and some of the HPeV types e.g. HPeV2, HPeV4, HPeV5 and HPeV6 appear to be relatively rare (see e.g [17, 18] ). Moreover, a common pitfall on estimating evolutionary rates is its underestimation due to mutational saturation of synonymous sites [41] [42] [43] [44] . By using a gamma distributed substitution model, we assured that rate variation among sites was allowed. Therefore, the effect of possible saturation of synonymous sites was alleviated by permitting a proportion of these sites to change at a higher rate [43] . In addition, we used partitioning in codon positions that allows different codon positions to have different substitution rates (and different amount of rate heterogeneity) (see Methods for details) [31, 45] thus further accommodating rate variation among synonymous and non-synonymous positions.