Biased estimates of phylogenetic branch lengths resulting from the discretised Gamma model of site rate heterogeneity

A standard method in phylogenetic reconstruction for representing variation in substitution rates between sites in the genome is the discrete Gamma model (DGM). Relative rates are assumed to be distributed according to a discretised Gamma distribution, where the probabilities that a site is included in each class are equal. Here, we identify a serious bias in the branch lengths of reconstructed phylogenies when the DGM is used, with the magnitude of the effect varying with the number of sequences in the alignment. We demonstrate the existence of the bias in both simulated dataset and using real HIV-1 sequences; in both cases branch lengths are overestimated. We show that the alternative "FreeRate" model, which assumes no parametric distribution and allows the class probabilities to vary, is not subject to the issue. We further establish that the reason for the behaviour is the equal class probabilities, not the discretisation itself. We explore the mathematics of for the phenomenon, showing how maximum likelihood branch lengths under the DGM may differ from the true ones used to generate the tree, and that branch lengths will be overestimated where there is a long-tail of fast-evolving sites in the true rate distribution, which is the usual situation in real datasets. We recommend that the DGM be retired from general use. While FreeRate is an immediately available replacement, it is known to be difficult to fit, and thus there is scope for innovation in rate heterogeneity models.