A Worked 3 species example

The Ornstein-Uhlenbeck process with migration: evolution with interactions


The Ornstein–Uhlenbeck (OU) process plays a major role in the analysis of the evolution of phenotypic traits along phylogenies. The standard OU process includes drift and stabilizing selection and assumes that species evolve independently. However, especially in plants, there is ample evidence of hybridization and introgression during evolution. In this work we present a statistical approach with analytical solutions that allows for the inclusion of adaptation and migration in a common phylogenetic framework. We furthermore present a detailed simulation study that clearly indicates the adverse effects of ignoring migration. Similarity between species due to migration could be misinterpreted as very strong convergent evolution without proper correction for these additional dependencies. Our model can also be useful for studying local adaptation among populations within the same species. Finally, we show that our model can be interpreted in terms of ecological interactions between species, providing a general framework for the evolution of traits between “interacting” species or populations.

Keywords : Migration, Ornstein–Uhlenbeck process, phylogenetic comparative methods (PCMs), trait evolution

1 Introduction

Comparing phenotypic traits or gene expression patterns across species or populations is intrinsically difficult since species or populations are inherently non–independent entities: species are related through phylogenies and populations through gene genealogies (Felsenstein, 1985). This historical component of the variation within and between species will also need to be accounted for if our aim is to understand trait evolution and in particular test whether traits evolved neutrally or under natural selection.

The first proposed model of continuous trait evolution was the Brownian motion (BM) model (Edwards, 1970; Felsenstein, 1985). This is essentially a neutral model of evolution with the trait oscillating around the ancestral value with linearly increasing variance over time. The lack of any stabilization in the process was quickly noticed and hence the Ornstein –Uhlenbeck model was proposed as an alternative (Felsenstein, 1988; Hansen, 1997). Since the seminal papers of Edwards and Felsenstein (Edwards, 1970; Felsenstein, 1985), various forms of the Ornstein-Uhlenbeck process have been used to assess the relative importance of genetic drift and selection in phenotypic trait evolution across species. Initial versions incorporating only drift and stabilizing selection, but recent forms of the OU process now consider simultaneously genetic drift and both stabilizing and directional selection (Beaulieu et al., 2012; Held et al., 2014). Until recently, however, species or populations were assumed to evolve independently, neither interacting nor exchanging genetic material. Since migration between populations or hybridization between species will tend to homogenize phenotypes directly by sharing the underlying genetic information, the application of such models at the species level (i.e. between populations within species, see discussion by Stone et al., 2011) could pose a serious problem. More generally, hybridization between closely related species or divergence with gene flow during species formation is likely more frequent than initially thought (Abbott et al., 2013), leading to the same problem at the phylogenetic scale where BM and OU models where classically applied. Ecological interactions between species can also affect phenotype evolution if the strength and/or direction of selection on the phenotype of a focal species is affected by the phenotypes of other species. For instance, competition can lead to character displacement and diversification. On the contrary, mimicry, as migration, will tend to homogenize phenotypes. Two recent studies combined phylogenetic structure and interactions during trait evolution. Nuismer and Harmon (2015) described an interaction framework combined with an OU evolution process. However the interactions between the species occur only through “the covariance among population mean phenotypes for species and ”. Drury et al. (in press 2016) further related to this model but also included species interaction in the “deterministic part” of the stochastic dynamics in their framework. In our work we focus on the migration problem and consider a selection—migration model where different species or populations are allowed to exchange genes thus affecting each others’ primary optima (in the spirit of Hansen, 1997). We show that Drury et al. (in press 2016)’s model is a special subcase of ours, which can thus also be applied to the species interaction problem. We then used simulations to assess the consequences of ignoring migration when testing for the presence of selection.

2 Phylogenetic models of selection with migration

2.1 Selection model

The Ornstein-Uhlenbeck process is defined by the following stochastic differential equation (SDE) for the trait ,


where is the standard Wiener process. The interpretation of the individual model parameters has been discussed at length (e.g. Butler and King, 2004; Hansen, 1997; Hansen et al., 2008). Shortly, is the optimal value for the trait, determines the magnitude of genetic drift affecting the trait, and controls the speed of the approach to the optimum, i.e. the strength of selection. The neutral case in the sense of absence of selection, , corresponds to the Brownian motion model. The expected value (mean) and variance of the process are


This means that after a long time the trait will stabilize around the optimal value with stationary oscillations of variance . A key parameter of this model is the half–life (Hansen, 1997), , which quantifies how much time is needed to lose half of the ancestral effect.

Figure 1: Left: a simulated pure birth tree on tip species (speciation rate , TreeSim R package; Stadler, 2009, 2011; R Core Team, 2013) , centre: Brownian motion simulated on this tree , right: Ornstein–Uhlenbeck process simulated on this tree (BM: , ; OU: , , , ; mvSLOUCH R package Bartoszek et al., 2012).

The trait follows the speciation pattern described by the phylogeny. The trait changes along each branch according to the diffusion process in (1) and at speciation it splits into two processes. Commonly these two daughter lineages evolve independently. Nearly all current phylogenetic comparative methods make this “no exchange or interaction between species” assumption. However, one can easily argue that this is not realistic biologically: species interact with each other due to e.g. ecological interactions, gene flow or hybridization.

2.2 Including migration

Assume we have species or populations under selective pressure present in the time interval and that there is gene flow between some of them. For simplicity and parallelism with other OU models, we will use the term species in the rest of the manuscript, though gene exchange is more likely to occur between populations of the same species. The mathematics of the model is the same whether one studies evolving species or populations, whereas the parameter interpretations and model implications can be different.

If the state of the process at time is a vector consisting of the current trait-values, then the state after an additional time is

The matrix describes the migratory gene flow between the species while the other parameters are the standard OU process parameters. Taking , we obtain the ordinary differential equation (ODE)

Equivalently, letting denote the matrix

in which all interactions between the species are coded, we have


Combined with the random noise we obtain an OU with migration stochastic model


where is a vector of independent Wiener processes. In Eq. (4) we may now recognize as a multivariate, –dimensional, OU process with “deterministic part” the matrix and optimum value . This means that one can obtain the moments of the process analytically. The knowledge of the mean and variance–covariance function then characterizes its normal law, as the OU model is a normal process.

We represent the eigen decomposition of the matrix by

Under the assumption that is non–singular the mean vector and covariance matrix will be (Bartoszek et al., 2012)


The remaining case, when is singular, amounts to some of its eigenvalues being . Assuming that eigenvalue is , in Eq. (5) we need to consider the limit

Then the mean and variance formulae of Eq. (5) hold with at appropriate diagonal entries of the diagonal matrix . In fact we will discuss an important submodel where this is the case, namely the model of Drury et al. (in press 2016).

The model setting laid out above is quite general. Each species tends towards its own optimum (which in particular may be the same as the optimum of all/some other species) with controlling the speed of the approach. Along its route it exchanges migrants (hence shares trait values) or more generally interacts with all/some other species as controlled by , hence . As a result the model is high–dimensional: with species the matrix is parametrized by parameters. This number increases further if changes over the phylogeny. Usually in a phylogenetic comparative method (PCM) setting there are only measurements – all from contemporary species. Therefore one needs to add further structure on to obtain a feasible framework for analysis, and furthermore make the parameters interpretable.

Let us write out in detail the migration component for species (remembering that in Eq. (4) )

Now let us assume that migration is conservative, i.e. that the number of immigrants and emigrants are equal. This translates to the sum of entries in each row of equalling , i.e.

Let us assume furthermore that the exchange of information is symmetric between all species, i.e. and for . This implies that the elements of are given by

Such a matrix is eigen decomposable – it has the simple eigenvalue and an eigenvalue of multiplicity equalling . The eigenvectors matrix and its inverse are

From these representations one easily obtains the moments of the process, Eq. (5).

The above assumptions on the process may seem restrictive. They, however, perfectly describe the “mean field” corresponding to the Wright island model in population genetics or to a situation where a number of exchangeable species constantly interacts with each other. The magnitude of migration is measured by . At the species level, this corresponds to a model of local adaptation as analyzed by Hendry et al. (2001). However, they only considered two populations.

2.3 Including the phylogeny

So far we have considered the evolution of traits without any branching structure – a fixed number of species evolve and interact within a given time frame. Ideally we would like to model the phylogenetic or populations history setup – where the trait evolves with time and species not only interact but also speciate (or populations split) and new lineages, that also may interact, arise.

In order to trace the evolution of traits in the species tree and to derive the law of the contemporary sample we introduce some additional notation. By we denote the time (the origin of the tree is ) of the th speciation event, i.e. at time the number of species changes from to . Hence is the duration of time when there were lineages present. In particular if the tree has no root branch but starts at the first bifurcation event we will have . We assume that in each of these intervals all parameters are constant. A parameter with subscript , e.g. , indicates the value of the parameter at the time when there are lineages present in the phylogeny. After each speciation event the two daughter lineages inherit their ancestral trait value and the traits in both lineages begin to evolve over time following the assumed OU dynamics with possible interactions. The subscripted trait vector will be used to underline that at time there are lineages present. Just after speciation the two daughter lineages are identical, hence they have to share the same mean, variance and covariance with other lineages. In principle this causes a momentary singularity of the process, which however immediately disappears. To formalize the mean and covariance of the phylogeny traits process at this instance, as changes dimension to , we introduce a branching operator , which increases the size of a vector by one element and the size of a square matrix by adding a row and a column. This operator copies the appropriate entry of the traits’ mean vector and the appropriate row and column of the variance–covariance matrix when speciation happens, resulting in a recursive (in terms of speciation events) formula for the mean and variance

In principle the operator should be indexed by the phylogeny and speciation number but we suppress this notation to avoid clutter.

Figure 2: An evolving trait on a phylogenetic tree with migrations indicated, by horizontal lines.

2.4 Including measurement error and missing observations

Our setting is targeted to maximum likelihood optimizations, which allow for analyzing measurement errors and missing observations.

Trait measurements at the tips of the phylogeny are typically species’ means obtained as averages over many individuals. The resulting intra–species variation adds a sample variance attached to each observation, which is straightforward to include but could be highly misleading to ignore (e.g. Hansen and Bartoszek, 2012). From the vector of observations, for species we can estimate the variance for the –th species as

where is the average for the –th species. These variances are then added to the diagonal of the variance–covariance matrix of the contemporary species, , given by Eq. (2.3). This is a common approach used in phylogenetic comparative software that allows for measurement error (see e.g. Bartoszek et al., 2012).

To correctly handle species with missing measurements, the rows of the design matrix, expectation vector, Eq. (2.3) and rows and columns of the variance–covariance matrix, Eq. (2.3) that correspond to measurements with missing values should be removed. In the case of each species being described by a single trait this is the same as just removing the “missing value” species from the data set. However, in the case of a multivariate analysis, when observations are just missing from some traits, there is no need to remove a whole species from the analysis. One just removes the missing components by removing the appropriate rows and columns as described above (Bartoszek et al., 2012).

3 Simulations — effects of migration and selection

To assess the effects of migration we simulated populations of individuals with migration between them. We considered a fixed star phylogeny of species with height equal to and tips. We do not have a phylogenetically structured population but rather a population descending from a common ancestor and all individuals are constantly interacting.

We draw trait values from a normal distribution with mean vector and variance–covariance matrix defined by Eqs. (2.3) and (2.3) for different parameter values. Time was always assumed as , the initial value , . The optimum was either constant at or , or half the population had or while the other half had . The adaptation parameter varied as and the migration parameter . Here, is an arbitrary value on the trait scale. For the special case this is consistent with the common interpretation of as a proportion of individuals coming from abroad.

We now present the most illustrative histograms of the effects of migration, while all the simulation results can be found in the supplementary material. In the simulations we can immediately see that, as expected, the migration parameter mixes the individuals. When all have the same optimum value then increasing brings all individuals closer around the mean, by decreasing the observed trait variance. This effect is thus similar to increasing selection so that neglecting migration could affect selection estimates. In fact this increase of the selection effect is clearly visible in the simulated trajectories presented in Fig. 3. In the contrasting situation of different optima then an increase of the migration parameter pushes individuals away from their optimum and towards an average of these optima. This is again visible in the simulated trajectories of Fig. 3. When migration is weak the trait values still manage to tend towards their optima. However in the case of strong migration (, ) the traits move away from their optima.

Figure 3: Simulated trajectories on the phylogeny depicted at the top. First column: , second column: , third column , first row: all species have , second row: half of species have , other half . The solid line branches have as the optimum while the dashed have . In all presented simulations , , . We can clearly observe that migration has an averaging out effect on the “optimum” value.
Figure 4: Histograms of estimates of values. First column: , second column: , third column , first row: all individuals have , second row: half of individuals have , other half . In all presented simulations , , .

The simulations presented in Figs. 3 and 4 suggest that not taking into consideration migration can seriously bias parameter estimation. We did a simulation study to test this and the resulting histograms are presented in Figs. 7, 8, 9, 6, 5 and Tab. 1. For each parameter combination we simulated –tip Yule trees with speciation rate (actually with speciation rate and rescaled all branch lengths by ) using TreeSim. On each tree we simulated a trait process under the law defined by the parameter combination. If (i.e. usual OU process) we used the mvSLOUCH package and when R code that drew trait values from the law described in this paper. Then, using mvSLOUCH we estimated parameters of the OU process based on the simulated data. The mvSLOUCH package does not take migration into account, hence we can observe the consequence of ignoring migration effects. When there were two values, half the tip species had one value and half the other. We tried to cluster this as far as possible (the trees are random). Nodes with indices (in phylo format) had one value and those with indices the other. We drew (with equal probability) which half had which value of . Afterwards the regimes were painted back on the tree using the Fitch algorithm (Fitch, 1971; Swofford and Maddison, 1987) implemented in the fitch.mvsl() function in mvSLOUCH.

The results of the simulation study show that the consequences of ignoring migration in a phylogenetic comparative study are as those indicated in Fig. 4 (populations evolving under an OU with migration). Estimates of (Fig. 5), the selection parameter are very heavily inflated. This is as the positive parameter brings the observations closer together (compare Figs. 4, 8) — indicating a stronger selection effect, but to a different optimum when there are multiple optima. We can see that with a single optimum it is difficult to estimate . This is not surprising, Cressler et al. (2015) observed the same issue. However, we can see that is much easier to estimate when there are multiple optima (Fig. 5, bottom left). The intuition is that when there is a switching of optima, especially near the tips, one can observe how quickly the trait moves from one optimum to the other — giving information on . A similar situation is for estimation of (Fig. 6). The presence of migration makes it more difficult to estimate and multiple optima make it easier.

In OU phylogenetic models it is often easier to estimate the stationary variance of the process, . As with strong selection the variance of the contemporary sample should be essentially given by this formula (cf. Bartoszek and Sagitov, 2015). With no migration (Fig. 7, left column) the parameter is well estimated. However when we increase migration the estimate of this value can differ significantly. The maximum likelihood estimation procedure tries to fit such values so that will correspond to the variance of the contemporary sample. As migration significantly decreases the variability around the mean (see Fig. 8) the estimation procedure (unaware of migration) will fixate in a parameter region where is small.

Estimating the value of the optimum parameter, , was one the original motivations behind phylogenetic comparative methods. This parameter (or vector of parameters if there are multiple optima) is usually the easiest to estimate (Bartoszek et al., 2012; Bartoszek and Sagitov, 2015; Cressler et al., 2015). In fact when there is a single optimum, migration just makes it easier to estimate it: migration decreases sample variance and increases the speed of convergence of the mean. However, when there are multiple optima (Fig. 9, bottom row) ignoring migration makes it difficult to identify the different optima. The averages of the subsamples (i.e. those under the different regimes) are much closer together, as migration mixes the individuals of the different populations. With very strong migration (Fig.9, bottom right) we can actually observe a very interesting effect: the ancestral and selective regime effects interact with each other. Each dot on the scatter plot is the estimated pair. When the ancestor was at then the value of , was estimated around while in we may observe a lagging ancestral effect. If there were no ancestral effects we would expect (Fig. 4) as migration is symmetric between all species. However, when there are ancestral effects they will show up in both estimates. One will have an interplay between the force towards the average of , ancestral effects and weak (compared to migration) selection effects.

In summary, ignoring migration can have a large impact on the identification of the strength of selection and the magnitude of disruptive effects. If populations are mixing and are otherwise in separate selective regimes it is impossible to estimate the optima of these regimes correctly. However, if migration is only between individuals inside the same selective regime, then identification of its optimum value is significantly easier.

Setup Parameter Mean Variance MSE rbias rMSE
Table 1: Summary of simulation results. For each setup and parameter we report the mean, variance, mean square error (MSE), relative bias (rbias) and relative MSE (rMSE) of the estimates. The means are visualized on histograms. The parameter is the sample variance, refers to the sample variance of individuals, to individuals. Migration is set to (setup 1 and 4), (setup 2 and 5), or (setup 3 and 6).
Figure 5: Histograms of estimates of . First column: , second column: , third column , first row: all individuals have , second row: half of individuals have , other half . In all presented histograms , , . The vertical line represents the true value. Note that the histograms are not all on the same scale. The left and right–most bars are values equal or more extreme than the appropriate –value.
Figure 6: Histograms of estimates of . First column: , second column: , third column , first row: all individuals have , second row: half of individuals have , other half . In all presented histograms , , . The vertical line represents the true value. Note that the histograms are not all on the same scale. The left and right–most bars are values equal or more extreme than the appropriate –value.
Figure 7: Histograms of estimates of the stationary variance of the OU model. First column: , second column: , third column , first row: all individuals have , second row: half of individuals have , other half . In all presented histograms , , . Hence the stationary variance is . The vertical line represents the true value. Note that the histograms are not all on the same scale. The left and right–most bars are values equal or more extreme than the appropriate –value.
Figure 8: Histograms of sample variances of data simulated under the migration model. First column: , second column: , third column , first row: all individuals have , second row: half of individuals have , other half . In all presented histograms , , . In the second row we plot separately the sample variances of individuals (black) and individuals (gray). The vertical line is the OU process stationary variance . Note that the histograms are not all on the same scale. The left and right–most bars are values equal or more extreme than the appropriate –value.
Figure 9: Histograms of estimates of . First column: , second column: , third column , first row: all individuals have , second row: half of individuals have (black estimates), other half (gray estimates). In all presented histograms , , . In the second row we plot separately the individuals (black) and individuals (gray). The vertical lines represent the true values. In the case of strong migration (bottom right figure) we do not present the histogram, but plot the actually estimated values in each of the runs. Note that the histograms in the two rows are not on the same scale. The left and right–most bars are values equal or more extreme than the appropriate –value.

It is important to notice that in the modelling setup considered here the layout of the regimes is known. Only the values of the optima are considered unknown. An interesting model would of course be where the layout is not known but is estimated from the data. However this is a very complex problem, with possible identifiability issues (see e.g. Ingram and Mahler, 2013b; Khabbazian et al., 2016). Allowing for modelling the regime layout is a further development step.

4 Comparison with ecological interaction between species

So far we built our model by considering migration between species or populations. However, the structure of the model allows for a wider interpretation as that of interaction between species mediated through their phenotypes, especially ecological interactions. Our model applies as far as these interactions are described by linear combinations of the phenotype of the different species, which then correspond to a particular matrix .

For instance, Drury et al. (in press 2016) used a special submodel of Eq. (4) to describe competition between evolving species and study Anoles morphology. Their model is the special case , with their . However, as they focused on the effect of competition, in their analysis. This essentially is a Brownian motion with migration model (without selection). We can see that the process does not have a stationary distribution as one of the eigenvalues of . In fact its variance tends to infinity. A slight difference between the model of Eq. (5) and Drury et al. (in press 2016)’s is that they did not assume that all species interact with each other. Instead they considered subsets of species interacting with each other. This can be seen in their Fig. , where different Anoles species interact if they live on the same island. However this setup means that and are block–diagonal with different “” in different blocks, but same and . Hence the eigen decomposition is the same as of “our” except that it is done independently for each block.

If we assume on the contrary that , the model in Drury et al. (in press 2016) can be interpreted as a mimicry model where selection favour a common phenotype among different species, used for example as a warning message against predators (Mallet and Joron, 1999). For each species , the equivalent ODE can be written as:


which corresponds to our equation (3) with and an matrix with elements and . We remark that this can be generalized to a weighted mean of phenotypes, hence species-dependent values of , if the different species do not contribute equally to the definition of the optimum (for example as a function of their abundance).

As for competition, this mimicry process does not have a stationary distribution as one of the eigenvalues is zero. However, contrary to the situation in Drury et al. (in press 2016), we should note that as the remaining eigenvalues are positive in this case, there is a stationary component in the mimicry process. If we condition on one species then the conditional distribution of the remaining ones will tend to stationarity (compare with the mvSLOUCH model Bartoszek et al., 2012). This can be biologically interpreted that if one arbitrarily chooses one species then the remaining ones will “shadow” it with stable fluctuations. This one species defines a randomly changing primary optimum (Hansen et al., 2008) for the others. It is important to notice that this does not imply that the remaining species will be in any way more similar to each other. Rather it implies that the deviations of these species from the arbitrarily chosen focal species will be exhibiting the same stochastic behaviour. The mimicry trait can randomly diverge but the variation among species for this trait will reach a stationary distribution. Biologically what is selected for is not the absolute value of the trait but the fact that all species are sufficiently similar.

Our model generalizes these special cases (where ) to more realistic conditions where there is stabilizing selection around fixed optima in addition to selection (stabilizing or disruptive) due to species interactions. Moreover, both migration and ecological interaction effects could be encoded in the matrix, although the two components would not be necessarily identifiable without additional biological information. Our model thus provides a general framework to study the evolution of traits under selection, drift and species interactions sensu lato.

Finally, we also note that Drury et al. (in press 2016) numerically solve a system of ODEs, their Eqs. (a) and (b), in order to obtain the mean and covariance structure of the system and in effect are able to calculate the likelihood. However the solution of their ODE system is exactly our Eq. (5). Hence an analytical solution is possible instead of a numerical one. As long as has a nice enough structure (e.g. is block diagonal) the solution can be written in closed form. For a general the solution will be in closed form up to the eigen decomposition of (for which there are numerous numerical libraries).

5 Discussion

The Ornstein–Uhlenbeck process has been the model of choice to assess the relative importance of genetic drift and selection in phenotypic trait evolution across species. The availability of large phenotypic datasets – in particular, genomewide transcriptomic datasets, that inherently are equivalent to a myriad of phenotypic traits – has led to a renewed interest in the OU process and its application to phenotypic traits evolution (e.g. Rohlfs et al., 2014; Rohlfs and Nielsen, 2015; Nourmohammad et al., 2015).

Until recently species were assumed to evolve independently of each other in OU models. However, this is not always realistic and new models have started to appear where interactions between species are permitted (Nuismer and Harmon, 2015; Drury et al., in press 2016). These studies were motivated by the presence of ecological interactions between species, while our motivation in the present work stemmed from the occurrence of demographic interactions, mainly the exchange of migrants between species. Hybridization and introgression are now well documented in both plants and animals (Abbott et al., 2013) and have played an important role in speciation and in the resulting species trees. As pointed out recently by Solís-Lemus et al. (2016) gene trees discordances have often been solely attributed to incomplete lineage sorting but gene flow is another possibility that has received less attention. However, ignoring gene flow can lead to the reconstruction of the wrong evolutionary history. This, of course, will also apply when reconstruction phenotypic trait evolution and it is therefore important to integrate gene flow in addition to ecological interactions in OU models.

An interesting finding of our study is that the linking of species through ecological or demographic interactions lead to the same modified OU process compared to the OU model with species interaction. In addition to modeling migration, our model also extend previous models of species interaction. Moreover, it also has implications for our ability to estimate the importance of stabilizing selection. Our simulations showed that ignoring migration when it is actually present can lead to severely overestimate selection. Below we discuss some consequences and possible extensions of our results.

The analytical derivations and the simulation study presented here show the (intuitively obvious) effect of gene flow. The gene flow between species clearly leads them to become more similar. The means of the different species tend to a weighted average of the optima of individual species and the variance of the species decays very quickly. Furthermore it is intrinsically difficult to distinguish between the effect of the pull towards the environmental optimum due to selection for a given species from the push effect of the interactions with different species. Hence, analyzing data where gene flow does occur between species with a OU process without migration estimation procedure could naturally lead to the conclusion that adaptation is rapidly moving to an optimum that is common for all species and the inference of evidence of convergent evolution. For instance, the program SURFACE that explicitly aims at identifying cases of convergent evolution by fitting OU models does not incorporate gene flow in the model (Ingram and Mahler, 2013a). Even other approaches to assess convergence in phenotypic evolution do not explicitly incorporate gene flow or refer to it as a possible contributing factor (Moen et al., 2016). Although, the observation that convergence, for instance in Anolis is spectacular at a regional scale but is not observed at a larger scale would fit well with an effect of gene flow on convergence pattern (references in Moen et al., 2016). In many cases the assumption that gene flow did not occur between species will be reasonable, but the evidence in favor of the existence of gene exchange with congeners in may species, especially in plants, suggest that care is required when drawing inferences based on a OU model that does not consider gene flow. It should also be noted that the trees used in the phylogenetic reconstruction are often based on the assumption that gene flow is absent and on a limited number of markers, often on solely cytoplasmic markers, and that they may therefore simply miss it. The availability of genomic data and of multiple individuals within species together with an increasing awareness of gene exchange among species and of the limits of the OU model (Cooper et al., 2016) should, however, alleviate this issue.

One of the original motivations for introducing phylogenetic comparative methods was to “correct” for the correlations between observations induced by the evolutionary tree. It is important to keep in mind that these correlations are in fact more than a “nuisance” factor since the dependency structure allows one to identify the phylogenetic half–lives (Hansen et al., 2008), something which would have been impossible from an i.i.d. sample. However, if the correlations are stronger, due to interaction effects, then the “tree correction” is not sufficient. Identity through similarity is mistaken for identity through convergence, resulting in over–estimation of selection. The rule of thumb is that if very fast adaptation is observed to very few regimes then one should consider if this might not be the result of past/present migration or interactions effects.

Adaptive migration models should be widely applicable in evolutionary biology. One immediate application is to include population structure in phylogenetic comparative studies. The tips of the tree are usually species but inside them there may be distinct subpopulations that can have diverged significantly if local adaptation is strong (Savolainen et al., 2013). Neither treating them as one (including intra–species variance Hansen and Bartoszek, 2012) or separate entities is fully satisfactory. In the first case one ignores the sub–population structure and assumes common adaptive regimes for all. In the second approach the problem is taken care of, but one ignores the fact that one is looking at a single species. The different sub–populations may mix, reproduce and generally interact. The migration model allows for correct treatment of such situations, with the size of controlling how much flow there is between the different sub–populations. More generally, the migration matrix could be parameterized using additional information on gene flow among diverging populations, typically using molecular markers. The OU model with migration is also a generalization of the two-population model proposed by Hendry et al. (2001) to analyze local adaptation. First, it generalizes it to any number of populations. Second, it allows to include a branching structure and not only a single stationary migration matrix structure. It thus provides an interesting basis for further local adaptation models.

Beyond the effect of migration on trait evolution, we have showed the equivalence with species interaction models and how it can be applied to competition or mimicry. We can also consider the use of OU with migration model to study e.g. host–parasite or plant–pollinator systems (combining interactions with phylogenies in such systems has been recently considered, e.g. Minoarivelo et al., 2014; Poisot and Stouffer, 2016). With the current phylogenetic comparative methods one would need to treat the two clades separately, e.g. modelling the evolution of the host and parasite trait on either of the trees. In a migration setup the two phylogenies can be included separately, with a very distant time of coalescence. Then the appropriate host–parasite pairs interact through the matrix.

Currently evolutionary models assume that speciation is an instantaneous event. This might be the case on long–time scales but with recent clades the time during which species split could be comparable with the tree height. Furthermore, phylogenetic tree estimation software, can result in non–binary splitting. Such multiple radiations are commonly interpreted not as true multi–furcations but rather as resulting from a lack of information. In both cases the migration framework allows one to use the tree in a phylogenetic comparative analysis. Around the point of multi–furcation/species splitting then one would model that all involved lineages interact with each other.

All of the above setups require the availability of software to estimate parameters of migration models. Such are most certainly under development. However, their inference power remains to be assessed. With a contemporary sample only it is already difficult to make inference on the value of in the absence of migration. One would not expect the addition of migration to make estimation any easier. It is a question whether without fossil measurements it will be possible to separately estimate the effects of and or only a combination of them. This immediately leads to the question whether different regimes optima are estimable or whether we can only estimate their linear combinations (species averages). There might be some hope however, as our simulation study indicated that (in the situation of ) with multiple regimes is well identifiable. Also the linear transformation of the optima vector depends only on and  — hence if one is able to identify them, or rather appropriate ratios of them (see Appendix) then adaptive regimes should be identifiable. All of this, however, requires a careful study alongside inference software development.


KB was supported by the Knut and Alice Wallenbergs Foundation. SG was supported by the French CNRS and the Marie Curie IEF Grant ‘SELFADAPT’ 623486. ML acknowledges support from the Swedish Research Council: VR grants 2012–04999 and 2015–03797 to which this work is related.

Appendix A Worked species example

Let us consider the simple model as in Fig. 10.

Figure 10: A model with constant migration between the populations. It is modelled as a three–dimensional stochastic process on the time interval . On the time interval all three coordinates evolve equally, on the time interval a split occurs two coordinates evolve equally while the third evolves in a correlated fashion. Then at the start of interval the two identical processes split into two populations and gene flow is allowed between all three populations.

We will go through all the calculations for this toy species model in detail. Let us write that the length of time interval is . Assume that the gene flow rates are the same everywhere and also let us assume that on each branch we have a constant (see notation in Fig. 10). We consider the process time interval by time interval.

  • On time interval we have only one species so no migration can take place. Therefore the system is described by the SDE

    and then at time

  • On the time interval we have two lineages between which there is gene flow. Therefore the system has to be described by a bivariate process

    To solve it we recall the following eigen decomposition of the “deterministic part” matrix

    Then we have at the end of the interval

  • On the time interval all three lineages are interacting hence the system has to be described by a trivariate process

    To solve it we recall the following eigen decomposition of the “deterministic part” matrix

    and the solution at the end of the interval is