 Methodology article
 Open access
 Published:
Inferring demographic parameters in bacterial genomic data using Bayesian and hybrid phylogenetic methods
BMC Evolutionary Biology volumeÂ 18, ArticleÂ number:Â 95 (2018)
Abstract
Background
Recent developments in sequencing technologies make it possible to obtain genome sequences from a large number of isolates in a very short time. Bayesian phylogenetic approaches can take advantage of these data by simultaneously inferring the phylogenetic tree, evolutionary timescale, and demographic parameters (such as population growth rates), while naturally integrating uncertainty in all parameters. Despite their desirable properties, Bayesian approaches can be computationally intensive, hindering their use for outbreak investigations involving genome data for a large numbers of pathogen isolates. An alternative to using full Bayesian inference is to use a hybrid approach, where the phylogenetic tree and evolutionary timescale are estimated first using maximum likelihood. Under this hybrid approach, demographic parameters are inferred from estimated trees instead of the sequence data, using maximum likelihood, Bayesian inference, or approximate Bayesian computation. This can vastly reduce the computational burden, but has the disadvantage of ignoring the uncertainty in the phylogenetic tree and evolutionary timescale.
Results
We compared the performance of a fully Bayesian and a hybrid method by analysing six wholegenome SNP data sets from a range of bacteria and simulations. The estimates from the two methods were very similar, suggesting that the hybrid method is a valid alternative for very large datasets. However, we also found that congruence between these methods is contingent on the presence of strong temporal structure in the data (i.e. clocklike behaviour), which is typically verified using a daterandomisation test in a Bayesian framework. To reduce the computational burden of this Bayesian test we implemented a daterandomisation test using a rapid maximum likelihood method, which has similar performance to its Bayesian counterpart.
Conclusions
Hybrid approaches can produce reliable inferences of evolutionary timescales and phylodynamic parameters in a fraction of the time required for fully Bayesian analyses. As such, they are a valuable alternative in outbreak studies involving a large number of isolates.
Background
Genomic data are increasingly used to investigate infectious disease outbreaks caused by microbial pathogens. Recent developments in sequencing technologies have made it possible to obtain data for a very large number of samples, at low cost and within a very short timeframe. Phylogenetic methods can make use of these data to infer their evolutionary dynamics, known as phylodynamic inference. For example, genome data obtained during the first months of the 2013â€“2016 Ebola virus epidemic were used to determine the time of origin of the outbreak and the basic reproductive number (R_{ 0 }) of the circulating strains [1, 2]. Some of the key requirements for these inferences are that the data must have sufficient genetic diversity and that they should be a representative sample of the circulating strains.
Serially sampled data are particularly useful because their sampling times can be used to calibrate the molecular clock. This consists of calculating the rate of evolution, which is the amount of genetic change that has accumulated per unit of time. The rate of evolution is key to infer an evolutionary timescale, typically represented by a phylogenetic tree where the branch lengths correspond to time, known as a chronogram. Some methods assume that the rate of evolution is constant over time, known as a strict molecular clock, but popular Bayesian implementations, such as that in BEAST [3, 4], include relaxedclock models that use a statistical distribution to describe rate variation across time and lineages (reviewed in [5]). Phylodynamic models can be used to estimate the epidemic growth rate (r), R_{ 0 }, and other parameters [6, 7]. Importantly, these models describe the expectation of the distribution of node times in the chronogram. As such, inferences drawn from phylodynamic models rely on accurate estimates of evolutionary rates and timescales. A number of statistical methods are available to assess the robustness of inferences of evolutionary rates and timescales; those that are most widely used are implemented under a Bayesian framework (reviewed in [8]).
Bayesian phylogenetic approaches allow sophisticated evolutionary models to be specified. For example, the evolution of a pathogen during an outbreak can be defined as an exponentially growing population with considerable evolutionary rate variation among lineages; which can be modelled by specifying a nucleotide substitution model, a relaxedclock model and an exponentialgrowth tree prior. The parameters for all these models are obtained simultaneously and their estimates correspond to posterior probability distributions, such that their uncertainty is a natural byproduct of the analysis. Bayesian methods require specifying a prior distribution for all parameters. Although specifying a prior distribution is not trivial for some parameters, their influence can be assessed by comparing them to the posterior. An advantage of specifying prior distributions is that it is possible to include previous knowledge about the data. As a case in point, a known probability of sampling can be represented with a prior distribution in birthdeath models [9].
Whilst Bayesian phylogenetic methods have many desirable properties, analysing large genomic data sets under complex models is often computationally prohibitive (e.g. [10, 11]). An alternative to full Bayesian methods is to conduct the analysis in several steps. In this hybrid approach the phylogenetic tree, evolutionary rates and timescales, and demographic parameters are estimated separately.
Phylogenetic trees can be rapidly estimated using various maximum likelihood implementations [12,13,14,15]. These methods assume a substitution model, but not a molecularclock or demographic model, such that the branch lengths of the trees represent the expected number of substitutions per site, and are known as phylograms.
Next, phylograms can be used to estimate evolutionary rates and chronograms, for example, using a recently developed molecular clock method based on leastsquares optimisation, called LSD (Least Squares Dating) [16]. LSD is more computationally tractable than Bayesian molecularclock methods, such that it is feasible to analyse genomic data sets with thousands of samples. Although LSD assumes a strict molecular clock, its accuracy is frequently similar to that obtained using more sophisticated Bayesian clock models [17]. Other nonBayesian molecularclock methods have also been developed recently with the purpose of analysing large genomic data sets [18,19,20].
Finally, a range of tools are available to infer phylodynamic parameters from a chronogram, such as that obtained using LSD. For example: TreePar uses maximum likelihood to fit birthdeath and skyline models [21]; BEAST2 [4] and RevBayes [22] can fit a range of birthdeath, coalescent, and Skyline models using Bayesian inference [7]; and approximate Bayesian computation (ABC) approaches that use tree summary statistics have recently been developed to fit phylogenetic epidemiological models [23, 24]. The main disadvantage of these approaches over those that are fully Bayesian is that the estimates are based on a single tree, such that uncertainties in tree topology, branch lengths, and evolutionary rates are ignored. A potential solution is to repeat the analysis using nonparametric bootstrap replicates, but combining the different sources of uncertainty under this framework is not trivial.
Here, we compare the following two methods to infer evolutionary rates and timescales, and demographic parameters:

(i)
The fully Bayesian method, implemented in BEAST2, to simultaneously infer the phylogenetic tree, evolutionary timescales and phylodynamic parameters;

(ii)
The hybrid method: phylogram inference using maximum likelihood in PhyML v3.1 [14], chronogram inference using LSD v0.3, and estimation of phylodynamics parameters in BEAST2 using Bayesian inference.
To compare the performance of these two methods, we analysed previously published whole genome SNP bacterial data sets of Mycobacterium tuberculosis Lineage 2 [25], Vibrio cholerae [26], Shigella dysenteriae type 1 [11], and Staphylococcus aureus ST239 [27]. Because these data sets have small numbers of samples (nâ€‰=â€‰63 for M. tuberculosis, nâ€‰=â€‰122 for V. cholerae, nâ€‰=â€‰121 for S. dysenteriae, and nâ€‰=â€‰74 for S. aureus) their analyses are computationally tractable using both approaches. We also demonstrate the unique potential of the hybrid approach by analysing two genomic data sets with larger numbers of sequences, which have been difficult to analyse using a fully Bayesian approach; a global sample of S. dysenteriae type 1 (nâ€‰=â€‰329) and S. dysenteriae type 1 lineage IV (nâ€‰=â€‰208) [11]. Finally, we validated the performance of the hybrid approach using a simulation experiment.
Results
Estimates of evolutionary rates and timescales
We compared estimates of rates and evolutionary timescales using the full Bayesian approach in BEAST2 and LSD. Because our data consist of SNPs, we used ascertainment bias correction by specifying the number of constant sites from the core genome. In BEAST2 we used both the strict and the uncorrelated lognormal (UCLN [28]) clock models. We investigated the degree of rate variation among lineages by inspecting the coefficient of rate variation, estimated in the UCLN model. This parameter is the standard deviation of branch rates divided by the mean rate. The data are considered to display clocklike behaviour if the distribution for this parameter abuts zero. Therefore, we used this parameter to select the clock model in BEAST2 for each data set, as suggested in previous studies [29, 30]. The M. tuberculosis data set was the only data set to support a strict clock over the UCLN model, whereas the remaining data sets favoured the UCLN model (Fig.Â 1). We set uniform prior distributions for the clock rate, the growth rate (r) and the scaled population size (Î¦). In the context of pathogen evolution, r determines the speed of spread of the pathogen in the host population, while Î¦ is proportional to the infected host population size at present.
The estimates of evolutionary rates and timescales from these different methods were largely congruent (Fig.Â 1). In all four cases, the 95% credible intervals for the evolutionary rate and age of the root node obtained with BEAST2 overlapped with the 95% confidence intervals obtained for the same parameters with LSD (Fig.Â 1). However, we observed some differences in the mean evolutionary rate estimates, with the estimates from BEAST2 consistently producing higher values than those from LSD. The largest difference in mean rate estimates was observed in M. tuberculosis, with a mean rate of 9.37â€‰Ã—â€‰10^{âˆ’â€‰8} (95% credible interval: 4.25â€‰Ã—â€‰10^{âˆ’â€‰8} â€“ 1.73â€‰Ã—â€‰10^{âˆ’â€‰7}) using BEAST2, and 1.10â€‰Ã—â€‰10^{âˆ’â€‰8} (95% confidence interval: 1.00â€‰Ã—â€‰10^{âˆ’â€‰10} â€“ 2.02â€‰Ã—â€‰10^{âˆ’â€‰7}) in LSD (see Fig.Â 1). In contrast we found more congruent mean rate estimates in the V. cholerae data set, with estimates of 7.20â€‰Ã—â€‰10^{âˆ’â€‰7} (95% credible interval: 5.87â€‰Ã—â€‰10^{âˆ’â€‰7} â€“ 8.65â€‰Ã—â€‰10^{âˆ’â€‰7}) for the BEAST2 and 6.76â€‰Ã—â€‰10^{âˆ’â€‰7} (95% confidence interval: 5.76â€‰Ã—â€‰10^{âˆ’â€‰7} â€“ 8.89â€‰Ã—â€‰10^{âˆ’â€‰7}) for LSD. The differences in estimates of the rootnode age were similar, with the largest difference in the mean rootnode age found in S. aureus ST239 (mean rootnode age of 1958 for BEAST2 and 1949 for LSD) (Fig.Â 1). In most cases, the estimates from BEAST2 were more uncertain with credible intervals that were wider than the confidence intervals from LSD. We investigated two aspects of phylogenetic data that can affect estimates of evolutionary rates; the topological uncertainty and the degree of clocklike variation. We found that the maximum likelihood trees were highly supported, according to local likelihood ratio tests (aLRT) [31] (which ranges from 0 to 1, for low to high branch support, respectively). The median aLRT values across nodes were 0.9 for M. tuberculosis, 0.83 for V. cholerae, 0.99 for S. dysenteriae type 1, and 0.92 for S. aureus.
Assessing temporal structure using a daterandomisation test
We assessed the reliability of our estimates of evolutionary rate and timescales by conducting a daterandomisation test [32, 33]. The motivation of this test is similar to that of roottotip regressions implemented in TempEst [34]. That is, to determine whether there is sufficient sampling in the data. However, roottotip regressions should be interpreted for visual inspection, as opposed to daterandomisations, which are a formal statistical test. The date randomisation test consists in repeating the analysis several times after randomising the sampling dates. The resulting rate estimates correspond to the expected values if there is no association between sampling times and genetic divergence. The data are considered to have strong temporal structure if the rate estimate obtained using the correct sampling times is not contained within the range of values from the randomisations. In a Bayesian context, 10 to 20 randomisations appear to be sufficient [33, 35]. We conducted this test in BEAST2 using 20 randomisations and in LSD using 100 randomisations (Fig.Â 2). Interestingly, the results from both tests were congruent, and consistent with visualisations of clocklike behaviour of the data using roottotip regressions (AdditionalÂ fileÂ 1: Figure S1). The M. tuberculosis data set had no temporal structure with either method (Fig.Â 2): the credible interval of the Bayesian estimate with the correct sampling times overlapped with those from all of the randomisations; using LSD, the estimate with the correct sampling times was around the lower threshold in the program, at 1.00â€‰Ã—â€‰10^{âˆ’â€‰10} subs/site/year, which also corresponds to the value obtained for most of the randomisations. The other data sets showed strong temporal structure with both daterandomisation tests: the Bayesian credible intervals using the correct sampling times did not overlap with those from any of the randomisations, and the estimates from LSD using the correct sampling times were not contained within the distributions of the 100 randomisations (Fig.Â 2).
Inference of phylodynamic parameters
We analysed the data sets using the exponentialgrowth coalescent model in BEAST2, which has two parameters, r and Î¦. Because these are compound parameters, they cannot be interpreted in an absolute scale without additional information about the size of the infected host population at present [36]. In most cases, the posterior distributions of both parameters were very similar when using either BEAST2 or the hybrid approach, with similar means and uncertainties (Fig.Â 3). Although the intervals overlapped in V. cholerae, S. dysenteriae, and S. aureus, the mode of the posterior distribution of Î¦ was higher when using the hybrid approach. The posterior distributions of r were almost identical across methods for the three data sets with temporal signal (Fig.Â 3). The uncertainty in estimates of this parameter did not include 0, except in the case of V. cholerae, suggesting that most of these bacterial data sets were undergoing population growth. Interestingly, the M. tuberculosis data set, which had no temporal structure, was the only data set to display large differences in estimates among the methods (Fig.Â 3).
Application: analysing large data sets using the hybrid approach
Having demonstrated good performance of the hybrid approach on small data sets with strong temporal signal, we applied it to analyse two published genomewide SNP data sets whose sample size was prohibitively large to analyse under a full Bayesian framework in the original publication. These data sets consisted of: (i) 329 samples of S. dysenteriae type 1 from [11], which included BEAST2 analysis of a subset of 125 samples; and (ii) 208 samples of lineage IV of S. dysenteriae type 1, which was represented by 61 samples in the BEAST2 analysis in the same study [11]. These three data sets displayed strong temporal structure according to the daterandomisation test in LSD, with rate estimates that were not contained within the range of estimates from 100 daterandomisations (Fig.Â 4). The evolutionary rate estimates from LSD were 5.93â€‰Ã—â€‰10^{âˆ’â€‰7} (95% confidence interval: 3.65â€‰Ã—â€‰10^{âˆ’â€‰7}  1.65â€‰Ã—â€‰10^{âˆ’â€‰6}) subs/site/year for S. dysenteriae type 1, and 7.04â€‰Ã—â€‰10^{âˆ’â€‰7} (95% confidence interval: 3.92â€‰Ã—â€‰10^{âˆ’â€‰7}  1.54â€‰Ã—â€‰10^{âˆ’â€‰6}) subs/site/year for S. dysenteriae type 1 Lineage IV (Fig.Â 4). Interestingly, the estimate of r for S. dysenteriae type 1 lineage IV was over an order of magnitude higher than that for the global data set of this bacterium, with a mean of 2.00â€‰Ã—â€‰10^{âˆ’â€‰2} for lineage IV compared with 3.40â€‰Ã—â€‰10^{âˆ’â€‰3} for the global data set. Importantly, the posterior distributions of r for these three data sets did not include zero, indicating epidemic growth (Fig.Â 4).
Validation using simulations
Although our empirical analyses suggest that the hybrid and the full Bayesian method can produce largely congruent results, it is unclear whether the methods are accurate. That is, whether they can recover the true parameter estimates. To investigate this, we conducted a simulation experiment. We simulated 100 whole genome data sets using similar parameters to those we inferred for our S. dysenteriae data set. We extracted the SNPs from the synthetic genomes and analysed them using the hybrid and full Bayesian methods, with the same settings that we used for the empirical data. Our daterandomisations in LSD indicated that all of these data sets had temporal structure, with pvalues of 0.00. The estimates for the age of the rootnode from both methods were very similar. However, it is important to note that our hybrid method uses a single tree, such that the age of the rootnode is a point value, whereas the full Bayesian analyses include uncertainty in this parameter. Accordingly, the estimates from LSD were very close to those used to generate the data (within 5Â years of the true value), and those from the full Bayesian method always included the true value within their credible interval. The estimates for the demographic parameters, r and Î¦, had credible intervals that always included the true value for both methods, with mean values that often matched those used to generate the data (Fig.Â 5a). Interestingly, in 10 randomly selected simulation replicates, we found that the credible intervals for the demographic parameters were very similar for both methods, with the hybrid approach sometimes producing more precise estimates. We found no estimation biases in any of the methods (Fig.Â 5a).
We conducted a second set of simulations of data with no temporal structure. To do this, we generated similar sequence alignments as described above, but we assigned random sampling times in our analyses in LSD and in BEAST2. This means that the molecular clock calibration is effectively uninformative. The age of the rootnode was over estimated by both methods. In LSD this bias was of over three orders of magnitude, whereas in BEAST2 it ranged between half and three orders of magnitude. The value of Î¦ was similarly overestimated in both methods. The growth rate, r, was underestimated by several orders of magnitude with the hybrid approach, but it tended to be overestimated with the full Bayesian method (Fig.Â 5b). A key result about the simulations with no temporal structure is that Î¦ was always incorrectly estimated, and the true value of r was only contained within the 95% credible interval in about 14% of the analyses using the full Bayesian method. Moreover, the estimates with the hybrid approach often displayed larger discrepancies with the correct values.
Computational demands of the Bayesian and the hybrid methods
The hybrid approach was several times faster than the full Bayesian approach. For example, the computation time for each randomisation of the V. cholerae data set each was about 2Â h using BEAST2, where as those in LSD took 1.23Â s (sec). However, a key aspect of the daterandomisation test in LSD is that the tree topology and branch lengths are fixed for all randomisations, where as they are reestimated for each randomisation in BEAST2. For the V. cholerae data set, a complete analysis using the hybrid approach took: 10.06Â min (min) to infer a maximum likelihood tree in PhyML, 1.23Â s to estimate the evolutionary rate and timescale in LSD, and 5Â min to infer r and Î¦ in BEAST2 to obtain effective sample sizes (ESS) of over 200 for all parameters (drawing 1â€‰Ã—â€‰10^{7} steps, with 1Â min per 10^{6} steps), for a total of about 15Â min, and 1/12th of the time required in BEAST2. Analysis of the full S. dysenteriae dataset from [11], the largest data set in our study, took 10.6Â s to analyse in LSD and 1Â h infer r and Î¦ BEAST2 (drawing 5â€‰Ã—â€‰10^{7} steps, with 1.2Â min per 10^{6} steps), for the 329 sampled sequences.
Discussion
Our results demonstrate that, as long as a strong temporal signal is present, the hybrid and fully Bayesian methods can produce congruent estimates of evolutionary parameters, even in cases where the data display substantial rate variation among lineages. These methods also yielded similar estimates of demographic parameters in data sets with strong temporal signal, indicating the hybrid approach is a reliable alternative to full Bayesian analyses. However, r appears to be more robust than Î¦ to mild differences in estimates of the rate and timescale. This probably occurs because the age of the rootnode plays an important role in the population size under the coalescent. In particular, the effective population size, and therefore Î¦, are known to scale positively with the age of the rootnode [37].
Obtaining congruent estimates between the two methods depends on whether the data meet certain criteria. In practice, it is important to verify that the trees have high branch support and that the data have strong temporal structure. The trees inferred here were highly supported, but it is likely that the hybrid approach will produce misleadingly precise estimates (i.e. with narrow confidence intervals) if branch support is low, because the demographic parameters will still be conditioned on a single, and possibly incorrect, tree obtained in step 1 that does not capture uncertainty in the topology. In contrast, in such circumstances the Bayesian method will simply integrate over phylogenetic uncertainty and yield wider credible intervals. Our simulations illustrate ideal conditions, in which the data evolve under the correct model and have strong temporal structure. In this case, we find that both methods produce accurate estimates with similar precision.
Our simulations of data with no temporal structure demonstrate, not only that the hybrid and full Bayesian methods will produce different estimates, but that they both tend to be inaccurate. In the absence of temporal structure, LSD often produces rate estimates at the lower threshold of the program, which was 10^{âˆ’â€‰10} here. This means that the timescale of the chronogram is overestimated. The value of Î¦ is also overestimated, which occurs because this parameter scales positively with the age of the rootnode [37]. Although, we found that r was also overestimated, this parameter is determined by the distribution of branches in the tree, such that its error is less predictable. The full Bayesian method produced estimates with smaller bias. We used uniform priors for Î¦ and r, and the prior for the age of the root was determined by the coalescent prior. It is likely that these parameters, especially Î¦, will be affected by different choice of priors. For empirical data with low temporal structure, the hybrid approach will likely be misleading because it is conditioned on a single tree which is probably incorrect. In such cases, it may be necessary to use the full Bayesian method approach because it is possible to include sources of molecular clock calibration via prior parametric distributions, at the expense of much higher computational demands. For instance, a reasonable calibration on the age of the rootnode might be sufficient to overcome low temporal structure and to obtain reliable estimates for Î¦ and r. To investigate this, it is important to verify that there exists a difference between the prior and posterior for parameters of interest (see Boskova et al. [38] for an investigation of the prior and posterior in Bayesian phylodynamics).
Our results show that the daterandomisation test in LSD appears to be as effective as it is in BEAST2, with the advantage of being much less computationally demanding. As a result, it is possible to use a larger number of replicates, which can improve the power of the test. Moreover, the sampling times under a Bayesian analysis of sequentially sampled data are informative about the tree topology. That is, they impose a high prior probability on trees that cluster sequences with similar sampling times, which can render the daterandomisation test unreliable, with an increase in type I error [39]. Moreover, in some phylodynamic models, the estimate of the age of the rootnode and the evolutionary rate are determined by a combination of the sequence data and their sampling times [38], such that assessing temporal structure via the date randomisation test is not trivial. The daterandomisation test in LSD does not suffer from these problems because sequence data alone, not tip dates, are used to infer the tree topology in maximum likelihood.
Critically, the rates estimated using the daterandomisation in test in LSD are not necessarily unimodal in their distribution. This occurs because a lack of temporal structure usually leads to very low rate estimates, which affects randomisations in LSD and in BEAST2. In the case of LSD, very low values for the rate will correspond to the lower threshold set in the program [17], which we arbitrarily set at 10^{âˆ’â€‰10} subs/site/year, such that most randomisations will have this value. As such, a reasonable approach to interpret the daterandomisation test in LSD is to ensure that the rate estimate with the correct sampling times is higher than those from at least 95% of the randomisations, following the frequentist onetailed pvalue of Î±â€‰=â€‰0.05.
Conclusions
As shown here, hybrid methods offer an attractive alternative to full Bayesian approaches for genomescale data sets with very large numbers of samples. The accuracy and precision of both methods are comparable, but hybrid methods can perform an analysis in a about an eighth of the time required for full Bayesian analyses. Nevertheless, some steps of the hybrid method used here require oversimplifications of the evolutionary process. For example, LSD always assumes a strict molecular clock, such that it is impossible to assess amonglineage rate variation or to pinpoint potential biological causes for why lineages have different rates. The choice of whether to use a hybrid method should be made based on what parameters a user wishes to interrogate. In the context of molecular epidemiology, demographic parameters (r and Î¦) and divergence time information are of primary interest, all of which appear robust to some amonglineage rate variation.
In this study, we used a simple demographic model, the exponentialgrowth coalescent. This model appears to be well suited when outbreak data are sampled at an early stage, but it makes several assumptions, including that the population of susceptible hosts is constant and that there is no population structure [6]. A better understanding of the data used here requires more sophisticated phylodynamic models, such as those that include changes in diversification parameters over time [40], and migration [41]. To this end, our results suggest that harnessing the power of such models and largescale genome sequencing can be done through hybrid approaches.
Methods
Data collection
Our bacterial data sets consisted of publically available genome data. We obtained all of our genomewide SNP alignments from a previous studies [11, 25, 27, 35]. These data sets are freely available online (github.com/sebastianduchene/bacteria_genomic_rates_data). These data have had regions with evidence of recombination removed using Gubbins v2 [42].
Phylogenetic analyses under the fully Bayesian approach
We analysed the sequence alignments in BEAST v2.4 using the sampling times for calibration, the GTRâ€‰+â€‰Î“ substitution model, the exponentialgrowth coalescent tree prior, and two clock models; the strict and the UCLN. We used the default priors for all parameters. Our Markov chain Monte Carlo (MCMC) sampling scheme consisted of a chain length of 5â€‰Ã—â€‰10^{8} steps, sampling every 10^{4} steps. We verified that the ESS for all parameters was at least 200. To determine whether the data had temporal structure, we conducted a daterandomisation test by randomising the sampling dates 20 times and repeating the analyses [33].
Phylogenetic analyses using the hybrid approach
We inferred phylogenetic trees using maximum likelihood in PhyML v3.1. We used the GTRâ€‰+â€‰Î“ substitution model, and a search strategy that combines the nearestneighbour interchange and subtree prune and regraft algorithms. To assess branch support, we calculated the aLRT score for each branch. To visually assess temporal structure, we conducted a regression of the roottotip distances as a function of the sampling times using TempEst v1.5 [34]. To determine the optimal root in this program we selected the position that maximised R^{2}.
We analysed the maximum likelihood trees (i.e. phylograms) in LSD v0.3 to infer the evolutionary rate and timescale. We set the sampling times as calibrations and allowed the program to determine the optimal position of the root. We constrained the branching times of the estimated chronograms such that daughter nodes must be younger than their parent nodes. To obtain an uncertainty around estimates of times and rates, we conducted 100 parametric bootstrap replicates of the branch lengths, as implemented in the program. Therefore, the uncertainty corresponds to the 95% confidence interval of the parametric bootstrap values. We conducted a daterandomisation test 100 times by randomising the sampling times in the â€˜dateâ€™ file and running LSD each time. In this version of the test, the phylogenetic tree topology and branch lengths are fixed.
We used the chronograms estimated in LSD to infer demographic parameters in BEAST2. This consists in setting the input file to calculate the posterior as the likelihood of the tree given the model parameters multiplied by the priors on the parameters. In the exponential growth coalescent there are two parameters; Î¦ and r. We used an MCMC chain length of 1â€‰Ã—â€‰10^{7} sampling every 10^{4} steps, and we verified that all parameters had ESS values of at least 200.
Simulations
We simulated whole genome sequence alignments using the parameters from our S. dysenteriae data set. To do this, we took the highest clade credibility tree from this data set inferred in BEAST2 and simulated the evolutionary rate using NELSI [29], according to an UCLN clock model. We used a mean rate of 10^{âˆ’â€‰6} subs/site/year and a standard deviation of 10^{âˆ’â€‰7}. We used SeqGen v1.3 [43] to simulate genome sequence alignments of 3,750,125 nucleotides using the GTRâ€‰+â€‰Î“ substitution model with the mean parameter estimates for the empirical S. dysenteriae data. Finally, we extracted the SNPs from these alignments and analysed using the same method as for our empirical data. For our simulations with no temporal structure we set random sampling times for our analyses in LSD and BEAST2. In all cases, we conducted a daterandomisation test in LSD, as used in our empirical data analysis.
Availability of supporting data
The datasets generated and/or analysed during the current study are available in the github repository, github.com/sebastianduchene/bacteria_genomic_rates_data
Abbreviations
 ABC:

Approximate Bayesian Computation
 aLRT:

Local Likelihood ratio test for branch support
 LSD:

Leastsquares dating
 MCMC:

Markov chain Monte Carlo
 UCLN:

Uncorrelated lognormal clock
References
Gire SK, Goba A, Andersen KG, Sealfon RSG, Park DJ, Kanneh L, et al. Genomic surveillance elucidates Ebola virus origin and transmission during the 2014 outbreak. Science. 2014;345:1369â€“72.
Holmes EC, Dudas G, Rambaut A, Andersen KG. The evolution of Ebola virus: insights from the 2013â€“2016 epidemic. Nature. 2016;538:193â€“200.
Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969â€“73.
Bouckaert R, Heled J, KÃ¼hnert D, Vaughan T, Wu CH, Xie D, et al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:e1003537.
Ho SYW, DuchÃªne S. Molecularclock methods for estimating evolutionary rates and time scales. Mol Ecol. 2014;23:5947â€“75.
Volz EM, Koelle K, Bedford T. Viral phylodynamics. PLoS Comput Biol. 2013;9:e1002947.
du Plessis L, Stadler T. Getting to the root of epidemic spread with phylodynamic analysis of genomic data. Trends Microbiol. 2015;23:383â€“6.
Rieux A, Balloux F. Inferences from tipcalibrated phylogenies: a review and a practical guide. Mol Ecol. 2016;25:1911â€“24.
Stadler T, Kouyos R, von Wyl V, Yerly S, BÃ¶ni J, BÃ¼rgisser P, et al. Estimating the basic reproductive number from viral sequence data. Mol Biol Evol. 2012;29:347â€“57.
Wong VK, Baker S, Pickard DJ, Parkhill J, Page AJ, Feasey NA, et al. Phylogeographical analysis of the dominant multidrugresistant H58 clade of Salmonella Typhi identifies interand intracontinental transmission events. Nat Genet. 2015;47:632â€“9.
Njamkepo E, Fawal N, TranDien A, Hawkey J, Strockbine N, Jenkins C, et al. Global phylogeography and evolutionary history of Shigella dysenteriae type 1. Nat Microbiol. 2016;1:16027.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and postanalysis of large phylogenies. Bioinformatics. 2014;30:1312â€“3.
Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQTREE: a fast and effective stochastic algorithm for estimating maximumlikelihood phylogenies. Mol Biol Evol SMBE. 2015;32:268â€“74.
Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59:307â€“21.
Price MN, Dehal PS, Arkin AP. FastTree 2â€“approximately maximumlikelihood trees for large alignments. PLoS One. 2010;5:e9490.
To TH, Jung M, Lycett S, Gascuel O. Fast dating using leastsquares criteria and algorithms. Syst Biol. 2016;65:82â€“97.
DuchÃªne S, Geoghegan JL, Holmes EC, Ho SYW. Estimating evolutionary rates using timestructured data: a general comparison of phylogenetic methods. Bioinformatics. 2016;32:3375â€“9.
Kumar S, Hedges SB. Advances in time estimation methods for molecular data. Mol Biol Evol. 2016;33(4):863â€“9.
Volz EM, Frost SDW. Scalable relaxed clock phylogenetic dating. Virus Evol. 2017;3(2):vex025. https://doi.org/10.1093/ve/vex025.
Ho SY, DuchÃªne S. Molecularclock methods for estimating evolutionary rates and timescales. Mol ecol. 2014;23(24):5947â€“65.
Stadler T. Mammalian phylogeny reveals recent diversification rate shifts. Proc Natl Acad Sci. 2011;108:6187â€“92.
HÃ¶hna S, Landis MJ, Heath TA, Boussau B, Lartillot N, Moore BR, et al. RevBayes: Bayesian phylogenetic inference using graphical models and an interactive modelspecification language. Syst Biol. 2016;65:726â€“36.
Poon AFY. Phylodynamic inference with kernel ABC and its application to HIV epidemiology. Mol Biol Evol. 2015;32:2483â€“95.
Saulnier E, Alizon S, Gascuel O. Assessing the accuracy of approximate Bayesian computation approaches to infer epidemiological parameters from phylogenies. PLoS Comput Biol. 2017;13:e1005416.
Merker M, Blin C, Mona S, DuforetFrebourg N, Lecher S, Willery E, et al. Evolutionary history and global spread of the Mycobacterium tuberculosis Beijing lineage. Nat Genet. 2015;47:242â€“9.
Devault AM, Golding GB, Waglechner N, Enk JM, Kuch M, Tien JH, et al. Secondpandemic strain of Vibrio cholerae from the Philadelphia cholera outbreak of 1849. N Engl J Med. 2014;370:334â€“40.
Baines SL, Holt KE, Schultz MB, Seemann T, Howden BO, Jensen SO, et al. Convergent adaptation in the dominant global hospital clone ST239 of methicillinresistant Staphylococcus aureus. MBio. 2015;6:e00080â€“15.
Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006;4:699â€“710.
Ho SYW, DuchÃªne S, DuchÃªne D. Simulating and detecting autocorrelation of molecular evolutionary rates among lineages. Mol Ecol Resour. 2015;15:688.
DuchÃªne S, DuchÃªne DA, Di Giallonardo F, Eden JS, Geoghegan JL, Holt KE, et al. Crossvalidation to select Bayesian hierarchical models in phylogenetics. BMC Evol Biol. 2016;16:115.
Anisimova M, Gascuel O. Approximate likelihoodratio test for branches: a fast, accurate, and powerful alternative. Syst Biol. 2006;55:539â€“52.
Ramsden C, Holmes EC, Charleston MA. Hantavirus evolution in relation to its rodent and insectivore hosts: no evidence for codivergence. Mol Biol Evol. 2009;26:143â€“53.
DuchÃªne S, DuchÃªne DA, Holmes EC, Ho SYW. The performance of the daterandomization test in phylogenetic analyses of timestructured virus data. Mol Biol Evol. 2015;32:1895â€“906.
Rambaut A, Lam TT, Carvalho LM, Pybus OG. Exploring the temporal structure of heterochronous sequences using TempEst (formerly pathOgen). Virus Evol. 2016;2:vew007.
DuchÃªne S, Holt KE, Weill FX, Le Hello S, Hawkey J, Edwards DJ, et al. Genomescale rates of evolutionary change in bacteria. Microb Genomics. 2016;2:e000094.
Boskova V, Bonhoeffer S, Stadler T. Inference of epidemiological dynamics based on simulated phylogenies using birthdeath and coalescent models. PLoS Comput Biol. 2014;10:e1003913.
Rosenberg NA, Nordborg M. Genealogical trees, coalescent theory and the analysis of genetic polymorphisms. Nat Rev Genet. 2002;3:380.
Boskova V, Stadler T, Magnus C. The influence of phylodynamic model specifications on parameter estimates of the Zika virus epidemic. Virus Evol. 2018;4:vex044.
Murray GGR, Wang F, Harrison EM, Paterson GK, Mather AE, Harris SR, et al. The effect of genetic structure on molecular dating and tests for temporal signal. Methods Ecol Evol. 2015;7:80â€“9.
Stadler T, KÃ¼hnert D, Bonhoeffer S, Drummond AJ. Birthâ€“death skyline plot reveals temporal changes of epidemic spread in HIV and hepatitis C virus (HCV). Proc Natl Acad Sci. 2013;110:228â€“33.
KÃ¼hnert D, Stadler T, Vaughan TG, Drummond AJ. Phylodynamics with migration: a computational framework to quantify population structure from genomic data. Mol Biol Evol. 2016;33:2102â€“16.
Croucher NJ, Page AJ, Connor TR, Delaney AJ, Keane JA, Bentley SD, et al. Rapid phylogenetic analysis of large samples of recombinant bacterial whole genome sequences using Gubbins. Nucleic Acids Res. 2014;43:e15.
Rambaut A, Grass NC. SeqGen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees. Bioinformatics. 1997;13:235â€“8.
Funding
SD was supported by a McKenzie fellowship from the University of Melbourne. ZAD is funded by strategic award #106158 from the Wellcome Trust of Great Britain. KEH is supported by fellowship #1061409 from the NHMRC of Australia.
Availability of data and materials
All the software used in this study is freely available and open source. The data are all available online github.com/sebastianduchene/bacteria_genomic_rates_data
Author information
Authors and Affiliations
Contributions
SD, DD, JLG, and KEH conceived and designed the experiments. SD, JLG, JH and ZD analysed the data. SD wrote the manuscript with input from all the authors. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional file
Additional file 1:
Figure S1. Roottotip regression for all data sets. The blue points correspond to tips in the tree. The black line represents the linear regression of roottotip distance as a function of the sampling time. The roottotip distance is measured by fitting the root of the tree that maximises R^{2}. (PDF 81 kb)
Rights and permissions
Open Access This 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.
About this article
Cite this article
Duchene, S., Duchene, D.A., Geoghegan, J.L. et al. Inferring demographic parameters in bacterial genomic data using Bayesian and hybrid phylogenetic methods. BMC Evol Biol 18, 95 (2018). https://doi.org/10.1186/s1286201812105
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1286201812105