Genetic constraints predict evolutionary divergence in Dalechampia blossoms

If genetic constraints are important, then rates and direction of evolution should be related to trait evolvability. Here we use recently developed measures of evolvability to test the genetic constraint hypothesis with quantitative genetic data on floral morphology from the Neotropical vine Dalechampia scandens (Euphorbiaceae). These measures were compared against rates of evolution and patterns of divergence among 24 populations in two species in the D. scandens species complex. We found clear evidence for genetic constraints, particularly among traits that were tightly phenotypically integrated. This relationship between evolvability and evolutionary divergence is puzzling, because the estimated evolvabilities seem too large to constitute real constraints. We suggest that this paradox can be explained by a combination of weak stabilizing selection around moving adaptive optima and small realized evolvabilities relative to the observed additive genetic variance.

. Name, abbreviation (Abbr.), geographical location given by the state in Mexico and GPS coordinates (WGS84), and meters above sea level (Elevation  Table S4. Posterior distribution of the Tovar G-matrix for the functional traits, provided in the external file "G_posterior_Tovar_functional_traits.csv" (comma delimited). The traits are on mean standardized scale. Each row is one sample from the joint posterior distribution. The columns are the elements of the G-matrix given row wise (or column vise, the matrix is symmetrical). The format of the column names are best described by an example. The column with name GAD:GAD gives the posterior distribution of the variance in GAD, while the column with name GAD:GSD gives the posterior distribution of the covariance between GAD and GSD. Note that GA and BA are simplified to GA and BA. The file can be read into R by the command: G = read.csv("filename.csv"). The posterior median of each matrix element can be computed, and printed in matrix format, by the command: matrix(apply(G, 2, median), ncol = sqrt(ncol(G))). Table S5. Posterior distribution of the Tulum G-matrix for the functional traits, provided in the external file "G_posterior_Tulum_functional_traits.csv". See table S4 for explanation. Table S6. Posterior distribution of the Tovar G-matrix for the bract traits, provided in the external file "G_posterior_Tovar_bract_traits.csv". See table S4 for explanation. Table S7. Posterior distribution of the Tulum G-matrix for the bract traits, provided in the external file "G_posterior_Tulum_bract_traits.csv". See table S4 for explanation. Figure S1 Scaling relationship between evolutionary rate (σ 2 rate ) and the evolvability estimates for a population of the other species for the functional traits. Grey circles represent 1000 uniformly distributed random directions (selection gradients). The solid line indicates the isometric relationship (a slope of 1) through the mean of the random directions. Crosses are the measured traits, squares are the directions with highest or lowest trait divergence, and black circles are the directions with highest and lowest evolvability out of the 1000 random directions. The circles, crosses, and squares are the modes from the posterior distributions and the grey lines give the 95% highest posterior density intervals. The vertical dotted lines are the posterior modes for the parameters named above each plot, and the vertical thick grey bars are their 95% highest posterior density intervals (see Table 2). The differences between minimum and maximum evolvability and the lowest and highest evolvability of the random directions are partly due to sampling error and partly due to bias in the estimates of e min and e max . Figure S2 Scaling relationship between evolutionary rate (σ 2 rate ) and the evolvability estimates for a population of the other species for the bract traits. See Figure S1 for explanation of symbols. Figure S3. Scaling relationship between respondability, computed along selection gradients, and evolutionary rates (σ 2 rate ), computed along corresponding response vectors, for the functional traits. Unit length response vectors are given by Gβ/||Gβ||, where β is the column vector of the corresponding unit length selection gradient and ||x|| denotes the norm of x. See Figure S1 for explanation of symbols. Note that the traits are not included in this figure as we do not have information about the direction of selection. Figure S4. Scaling relationship between respondability, computed along selection gradients, and evolutionary rates (σ 2 rate ), computed along corresponding response vectors, for the bract traits. See Figure S3 for further explanation.

Appendix A: Tracking optimum model
Here we consider a simplified model of a trait tracking an optimum that moves according to an Ornstein-Uhlenbeck process. We assume the selection gradient on the trait is = −2 ( − ), where is (log) trait mean and is the optimum. The parameter is the (mean-scaled) curvature of an assumed quadratic fitness function. The model is where , , and are differentials in the trait mean, optimum and time, respectively, is the evolvability, is white noise (i.e. the stochastic differential of a unit variance Brownian motion), and and are parameters describing the Ornstein-Uhlenbeck process. By use of Ito's formula for stochastic differentiation we can then form the stochastic differentials: By integrating these over the joint distribution of and , we can form differential equations in the first and second moments of the distribution where E[ ] denotes expectation of . Because the model is a system of linear stochastic differential equations, the distribution of the variables is multivariate normal, and the time development of these moments fully specify the time development of the distribution. These equations can be solved to find the time development of the system, but for the present purpose it suffices with the stationary distribution, which can be found by setting all the differentials to zero and solve for the moments. This gives where = 2 ⁄ is the stationary variance of the optimum, . Note that if the dynamics of the optimum converges on a white noise, i.e. that → ∞ and stays constant, the variance in the trait mean goes to zero. This is because the optimum moves so fast that the trait can not track it. More generally, if the characteristic time of the optimum, 1 ⁄ , is much shorter (i.e. faster) than the characteristic time of the adaptation, 1 2 ⁄ , then the variance of the trait mean goes to zero, because the trait mean can not track the movements of the optimum. If, on the other hand, the adaptation is much faster than the movements of the optimum, the trait mean will track the optimum perfectly, and the variance of the trait mean converges on the variance of the optimum, .
If the Ornstein-Uhlenbeck process for the optimum converges on a Brownian motion, that is, → 0, and is finite, the variance of the optimum is no longer stationary and evolves as . Under these conditions the variance of the trait mean evolves as where Exp is the exponential function and the approximation is for large t. This means that the variance of the trait mean eventually settles on the same rate of increase as the variance of the optimum, but with a constant lag that is inversely proportional to the evolvability.

Appendix B: Construction of the Phylogeny
We constructed the phylogenetic relatedness as follows. Genomic DNA was extracted from leaf tissue using E-Z 96 Pland DNA Kit (Omega Bio-Tek, Norcross, GA, USA). Multiplex PCR reactions were prepared in a total volume of 10 μl containing 5 μl of 2x type-it Microsatellite PCR Kit (QIAGEN, Hilden, Germany), 1 μl of extracted DNA (10 ng) and 1 μl primer (vary in concentration 0.1-0.4 μM). The thermal condition was set to 95 ° C for 5 min for initial denaturation; 10 initial cycles as touchdown at 94 ° C for 30 s, 60-50 ° C for 45 s, 72 ° C for 45 s; 25 cycles at 94 ° C for 30 s, 50 ° C for 45 s, 72 ° C for 45 s; and final extension at 72 ° C for 10 min in an ABI 9600 thermal cycler (Applied Biosystems, Foster City, California, USA). The 5' end of forward primers for 70 microsatellite loci was labeled with four different fluorophores (6-FAM, HEX, NED, and PET); and fragment analysis was conducted with an ABI 3130xl Genetic Analyzer (Applied Biosystems). Microsatellite alleles were scored using GeneMapper software version 4 (Applied Biosystems). Seventy microsatellite markers (Table S8) (Table S9). A neighbour-joining phylogenetic tree (Saitou & Nei, 1987) was constructed at the population level based on the genetic distance D A described in Nei et al.(1983) using Populations 1.2.31 program freely available at http://bioinformatics.org/~tryphon/populations/. The branch lengths of the phylogenetic tree were nearly of equal length when the tree was rooted on the centre, and was further ultrametrized by the function "chronopl" in the R-package "ape" (Paradis et al. 2004) with the parameter lambda set to unit, which gives a clock-like model.   Appendix C: The R-package "evolvability" The package evolvability contains functions for calculating evolvability parameters described in Houle (2008, 2009). The package is distributed through http://cran.rproject.org/. Updates and additions of this package will be published on the CRAN package repositories. The following functions are currently implemented: evolvabilityMeans calculates average evolvability, conditional evolvability, respondability, autonomy, integration, maximum evolvability and minimum evolvability from a given G-matrix using the approximation equations in Hansen and Houle (2008) evolvabilityMeansMCMC calculates the same as evolvabilityMeans, only over a posterior distribution of a G-matrix. evolvabilityBeta calculates evolvability, conditional evolvability, respondability, autonomy, and integration along each given selection gradient. The function also calculates averages, minimum, and maximum values over a set of given selection gradients.
evolvabilityBetaMCMC calculates the same as evolvabilityBeta, only over a posterior distribution of a G-matrix.
evolvabilityBetaMCMC2 calculates evolvability, conditional evolvability, respondability, autonomy, and integration along one selection gradient when both the G-matrix and the selection gradient is given by a posterior distribution.
randomBeta generates a set of selection gradients uniformly distributed in a hypersphere of dimension k. meanStdG mean standardizes a G-matrix given the corresponding means. Note that mean standardization is only meaningful for traits on a ratio or log-interval scale. Traits on these scales have a meaningful (unique and non-arbitrary) zero value. meanStdGMCMC mean standardizes the posterior distribution of a G-matrix by the corresponding posterior distribution of the means.