Results

The most supported molecular clock model was the lognormal relaxed molecular clock, with a Bayes factor with the next most supported model (strict molecular clock) equal to 225, therefore we used this molecular clock prior in all subsequent analyses. Under the SABDSKY model, most parameter estimates were very similar, so we report them jointly unless otherwise specified.
The mutation rate posterior median was 4.46 X 10-3 per year (95% HPD: 3.72 X 10-3 - 5.36 X 10-3 per year), with a coefficient of variance of 0.5 (95% HPD: 0.28-0.9). The root of the tree was estimated to be in 1995.36 (95% HPD: 1995.09 – 1995.58) and the origin in 1995.23 (95% HPD: 1994.92 – 1995.52). The sequence KT006745 (sampled in South Australia in 2009) had a sampled ancestor posterior probability of 38.3% and it was always paired, when sampled ancestor, with the sequence KF494930 (sampled in South Australia in 2010).
Across the four time intervals after 1995, the median estimatedR e had a slightly increasing trend until 2010 followed by a dramatic decline (Figure 1). However, based on the 95% Highest Probability Density (HPD) intervals, no significant difference was detected across the first three epochs in theR e estimate. The range of values forR e in the first five years after the release of the virus was very low (0.76 – 1.4), successively increasing to 0.9 – 1.4 and then declining to below 1 after 2010 (0.2 – 0.9) (Figure 1). The range of values for δ was 0.9 – 5.6 (median=2.6) and the sampling proportion parameter was estimated to be 0.008 (0.0005 – 0.0277) in the SABDSKY model with a common sampling proportion parameter across all the epochs (after 1995), while δ was narrower (0.6 – 2.2, median=1.3) and the sampling proportion showed a declining trend (Figure 2) in the SABDSKY model with a variable sampling proportion parameter. Results from the BDSKY analyses were similar to SABDSKY although with wider 95% HPD (Figure 1), with the exception that δ estimate was much smaller (0.5, 0.17 – 1.8) and the sampling proportion parameter had a much wider 95% HPD (0.02 – 0.9) in the BDSKY model with a common sampling proportion parameter across all the epochs as well as in the BDSKY model with a variable sampling proportion parameter (Figure 2).
The BSP analysis, correctly detected the initial dramatic increase in the virus population size at the beginning of the initial outbreak (Figure 3), despite the paucity of sequence data in the first few years since the release of the virus (Figure S1). The analysis suggests a gradual increase in the virus population size until 2007 and another substantial increase between 2007 and 2009. The virus population size then stabilised in the last five years or so of the analysis. The results of the skyride analysis were similar, and, although they had larger HPD intervals and the second increase in the population size between 2007 and 2009 was not as evident (Figure S2), these are not discussed further.
The analysis with the BDMM did not converge. We further attempted to reduce the time interval for the estimation of R eand fixed δ across the different epochs to reduce the number of parameters to be estimated with no success. The estimation of the sampling proportion and type change events proved particularly difficult and, even with runs in excess of 1.6 X 109 iterations the latter parameter failed to converge.