Trait evolution with incomplete lineage sorting and gene flow: the Gaussian Coalescent model

Most phylogenetic comparative methods use a species-level phylogeny, ignoring the effect of incomplete lineage sorting (ILS) and hemiplasy on the traits of interest. We consider here a trait controlled additively by one or more unknown loci. Their gene trees may differ from the species phylogeny due to ILS. The species phylogeny can be any phylogenetic tree or network or admixture graph on any number of taxa or populations. Our process of trait evolution accounts for introgression represented in the phylogeny by hybrid (or admixture) edges, which capture gene flow, admixture or hybridization events. The process models the trait in each individual, for any number of individuals from each taxon. Our model allows for polymorphism in the ancestral population at the root of the species phylogeny, and we can derive the heritable within-population variation expected from our model. Even if each locus evolves according to a Brownian motion, the joint distribution of the trait across all measured individuals is not generally Gaussian due to ILS. We provide a Gaussian approximation, named the Gaussian Coalescent (GC), and show how to compute its variance matrix efficiently using a single traversal of the species phylogeny. In simulations, this model is much more accurate than the model ignoring ILS. In simulations and on a data set of tomato floral traits, it is favored over the standard Brownian motion model with extra within-population variance. The GC model opens new avenues for various phylogenetic methods, accounting for hemiplasy and gene flow simultaneously on any general network or admixture graph. It is implemented in phylolm v2.7.0 and in PhyloTraits v1.2.0.