Substitution rates and dates of ancient divergence were estimated with Bayesian MCMC in BEAST version 1.5.3 [35] , with collection times in years used as calibration points in the clock model. The youngest strain was collected in 2008, which sets this year as the origin for past time estimates. Each dataset was evaluated individually for best fitting substitution model, which ranged from HKY+C 4 +I to GTR+C 4 +I. However, analyses performed under GTR family models neither converged nor mixed well, possibly due to an insufficiency of data to estimate these highly parametric substitution models. Hence, the simpler, less parameter rich, HKY+C 4 +I model was used throughout the BEAST investigation. We tested the impact of using a GTR model by running an analysis for 20610 6 generations. Estimates for the parameters of interest were largely concordant (data not shown), albeit the analyses returned very low ESS and much wider confidence intervals. Pairwise comparisons of Bayes factors calculated in Tracer selected the uncorrelated lognormally distributed relaxed-clock (UCLN) and the Bayesian Skyline coalescence model [36] as the best fitting clock and demographic models following the procedure in Hon et al. [37] . We defined two partitions that separated first and second positions from third codon positions. For each analysis, four independent MCMC chains were run for 20610 6 generations and their log output combined with 10% burn-in samples discarded. Tracer was used to determine the degree of mixing, shape of the probability density distribution, median and highest posterior density regions at 95% (HPD) for the relevant parameters. The modes and parameters of the posterior distributions were estimated using the distribution fitting software EasyFit 5.3 (MathWave Technology). For all analyzed parameters, we modeled the posterior distributions with gamma distributions. The analytical framework of the BEAST analyses is presented in Figure 1 and the details are explained below.