Detecting Individual Sites Subject to Episodic Diversifying Selection

The imprint of natural selection on protein coding genes is often difficult to identify because selection is frequently transient or episodic, i.e. it affects only a subset of lineages. Existing computational techniques, which are designed to identify sites subject to pervasive selection, may fail to recognize sites where selection is episodic: a large proportion of positively selected sites. We present a mixed effects model of evolution (MEME) that is capable of identifying instances of both episodic and pervasive positive selection at the level of an individual site. Using empirical and simulated data, we demonstrate the superior performance of MEME over older models under a broad range of scenarios. We find that episodic selection is widespread and conclude that the number of sites experiencing positive selection may have been vastly underestimated.


Introduction
Following the introduction of computationally tractable codonsubstitution models [1,2] nearly two decades ago, there has been sustained interest in using these models to study the past action of natural selection on protein coding genes.Positive selection can be inferred whenever the estimated ratio (v) of non-synonymous (b) to synonymous (a) substitution rates significantly exceeds one (reviewed in [3] and [4]).In the original models, the v ratio was shared by all sites in an alignment, providing little power to detect the signature of positive selection.Indeed, even among classical examples of positively selected genes [5,6,7], most substitutions are expected to be neutral or deleterious [8].Consequently, relatively few genes in which mean v estimates are significantly greater than one are expected to exist, e.g.only 35=8079 were found in a human -chimpanzee genome-wide comparison [9].
Random effects codon-substitution models [10] permitted v to vary from site to site, which made it possible to identify instances when positive selection had acted only upon a small proportion of sites.Such site-level models can detect which positions in a sequence alignment may have been influenced by diversifying positive selection, e.g.[11,12].However, these models posit that diversifying selective pressure at each site remains constant throughout time, i.e. affects most lineages in the phylogenetic tree, (Figure 1A), and there are very few cases where this assumption is biologically justified (see [13,14,15,16] for examples of models that allow selection to vary throughout the tree).When a site evolves under purifying selection on most lineages, site methods which assume v is constant over time may be unable to identify any episodic positive selection, since they will likely infer vv1 [17].It has been noted that positive selection is more readily identified in smaller alignments: counterintuitively, including additional sequences may cause sites to no longer be detected [18,19].This phenomenon could be readily explained by purifying selection on some lineages masking the signal of positive selection on others.
We present a mixed effects model of evolution (MEME), based on the broad class of branch-site random effects phylogenetic methods recently developed by our group [20].MEME allows the distribution of v to vary from site to site (the fixed effect) and also from branch to branch at a site (the random effect, Figure 1B).Our approach provides a qualitative methodological advance over existing approaches which integrate site-to-site and lineage-tolineage rate variation, e.g. the branch-site methods [17] or codonbased covarion models [13].MEME can reliably capture the molecular footprints of both episodic and pervasive positive selection, a task for which current models are not well suited.Using empirical sequence data sets spanning diverse taxonomic categories and gene functions, along with comprehensive simulations, we demonstrate that MEME matches the performance of traditional site methods when natural selection is pervasive, and that MEME reliably identifies episodes of diversifying evolution affecting a small subset of branches at individual sites, where site methods often report purifying selection at the same site.For most empirical data sets analyzed here, episodic selection appears to be the dominant form of adaptive evolution.The biological implications of this type of selection are discussed for each specific data set.We conclude by providing practical guidelines for applying MEME to biological data, and argue that while it is possible to reliably identify sites or branches subject to episodic diversifying selection, statistical power to detect individual branchsite pairs evolving adaptively is inherently limited by a small sample size available for such inference.

Methods
At its core, our approach uses phylogenetic models to describe the evolution of codon characters along a branch in a phylogeny by a continuous-time stationary Markov process.Given a phylogenetic tree t, with B branches and a vector of relative branch length parameters T~(t i ,i~1 . . .B), the probability of changing from codon i to j at a site along branch b in time t b , is recorded in the (i,j) element of the transition matrix M b (t b )~e Qt b , where Q is the rate matrix.The elements Q~q ij È É parameterize the instantaneous rate of substitution of codon i with codon j: Figure 1.The standard random effects approach and samples.A) The standard random effects approach, in which the rates vary randomly over sites but are constant over branches.Different values of v are showed in different colors.B) Samples from our new random effects approach [20], used by MEME, in which the rate on each branch is drawn independently of the rate on any other branch.All possible assignments of rates to sites are considered.doi:10.1371/journal.pgen.1002764.g001

Author Summary
Identifying regions of protein coding genes that have undergone adaptive evolution is important to answering many questions in evolutionary biology and genetics.In order to tease out genetic evidence for natural selection, genes from a diverse array of taxa must be analyzed, only a subset of which may have undergone adaptive evolution; the same gene region may be under stabilizing or relaxed selection in lineages leading to other taxa.Most current computational methods designed to detect the imprint of natural selection at a site in a protein coding gene assume the strength and direction of natural selection is constant across all lineages.Here, we present a method to detect adaptive evolution, even when the selective forces are not constant across taxa.Using a variety of well-characterized genes, we find evidence suggesting that natural selection is generally episodic and that modeling it as such reveals that many more sites are subject to episodic positive selection than previously appreciated.
q ij (a,b,P,H)~a d(i,j) counts the number of nucleotide differences between codons i and j. a and b parameterize the rates of synonymous and nonsynonymous substitutions, respectively.h ij (comprising H) are the nucleotide mutational biases, which we model using the 5parameter general time reversible nucleotide model.p ij (comprising P) denote the equilibrium frequency parameters.Our estimate (denoted throughout as P P) uses nine position-specific frequency parameters for the target nucleotides [1], corrected for the absence of stop codons using the CF3|4 estimator [21].The likelihood of observing the site is calculated using the pruning algorithm [22] given the data, the tree (t), the instantaneous rate matrix (Q), and the branch lengths (T).
To model the evolution of a site in an alignment in a manner that treats the non-synonymous rate (b) at each branch b as a random draw from one of K selective categories, we introduce a variable, c b , which can take values from 1 . . .K.An assignment of categories to all B branches, is described by the configuration vector C~(c 1 , . . .,c B ) of branch categories.We assume that the category on each branch is independent of that on all other branches, and that each category has an associated probability, p(c b ), for each branch.Next, we seek to marginalize the likelihood of each site D over all branch configuration vectors: Since this sum is over possible configurations, it has B K terms, and would appear infeasible, unless B is small.However, if we assume that branch categories are independent, p(C)~P B b~1 p(c b ), then the sum can be computed directly using the pruning algorithm by replacing the transition matrices with mixtures of transition matrices (see [20] for the derivation).If M b is the transition matrix on branch b, and we denote Felsenstein's algorithm, which computes the probability of observing D given a transition probability matrix for every branch, as F (M 1 , . . .,M B ), then: where M c b b associates a transition matrix at each branch with a category.We have thus constructed a tractable model where the process at every branch is a random draw from a set of K categories.
In [20], we used this result to develop a model where each branch had a set of v values and proportion parameters common to all sites.The goal was to identify lineages with a proportion of sites evolving with vw1.Here, we let each site have a set of free parameters governing the strength of selection for two discrete categories, and weights for each category, and these parameters are shared for all branches at that site.The goal is to detect sites where a proportion of lineages are evolving with vw1.

The MEME test for episodic diversifying selection
The fitting of MEME to an alignment of coding sequences proceeds in three stages: First, the MG94|REV codon model with an alignment-wide v~b=a is fitted to the data using parameter estimates under a GTR nucleotide model as initial values.Although in some cases nucleotide branch lengths may be a good approximation to codon branch lengths [23,24], recent results indicate that in other instances, nucleotide models can significantly underestimate branch lengths and possibly bias downstream inference [25].The resulting maximum likelihood estimates, Ĥ H and t t b , for each branch b[1 . . .B, are used in the site-by-site analyses in the next two steps.Thus we are assuming that the relative branch length and mutational bias parameters are shared across sites and are well approximated by those estimated under a simpler codon model.However, the absolute branch lengths also depend on the site-and model-specific rate parameters below.
Second, at each site, we first fit the alternative random effects model of lineage-specific selective pressure with two categories of b: b { ƒa and b z (unrestricted).The probability (p(c b ) in equation 1) that branch b[1 . . .B is evolving with b b ~b{ , is 0ƒq { ƒ1, and the complementary probability that it is evolving with b b ~bz is q z ~1{q { .By equation 1, the phylogenetic likelihood at a site, marginalized over all 2 B possible joint assignments of b b , is equivalent to computing the standard likelihood function with the following mixture transition matrix for each branch b: Consequently, the alternative substitution model includes four parameters for each site, inferred jointly from all branches of the tree: b { ,b z ,q { and a.These form the fixed effects component of the model.Estimating a separately for each site accounts for the site-to-site variability in synonymous substitution rates [26].
Lastly, at every site, we fit the model from the previous step, but with b z ƒa: our null model.Using simulated data, we determined that an appropriate asymptotic test statistic for testing most worstcase null of of b z ~b{ ~a is a 0:33 : 0:3 : 0:37 mixture of x 2 0 ,x 2 1 and x 2 2 (see Text S1).Mixture statistics of this form often arise in hypothesis testing where model parameters take values on the boundaries of the parameter space, and closed-form expressions for mixing coefficients are difficult to obtain [27].
Throughout the manuscript, we compare MEME to the fixed effects likelihood approach, introduced in [24] (see Text S1 for motivation).The procedure used by FEL differs from MEME in that a single pair of a,b rates are fitted at each site (no variation over branches) in Step 2, and the test in Step 3 is to determine if a=b.Positive selection is inferred by FEL when b bwâ a and the pvalue derived from the LRT is significant, based on the x 2 1 asymptotic distribution.

Detecting individual branches subject to diversifying selection at a given site
If the LRT indicates that a particular site (s) is subject to episodic diversifying selection, it may be of interest to explore which branches at that site have undergone diversification.The empirical Bayes (EB) procedure originally used to identify individual sites subject to diversifying selection in random effects models [28], can be readily adapted here.To compute the empirical posterior probability at branch b that b~b z wa, we apply Bayes' theorem, using D s to denote the data at site s and Ĥ H to denote all the maximum likelihood parameter estimates from the alternative MEME model fitted to site s: To compute the two likelihood terms P(D s Db b ~bz ) and P(D s Db b ~b{ ), we use q { ~0 and q { ~1, respectively, for the model assigned to branch b in equation 2. The rest of the branches employ the matrices fitted under the alternative model of MEME.
for each branch b, we evaluate the empirical Bayes factor for the event of observing positive selection at each branch: When EBwKw1, sequence data increase the prior odds of observing selection at the branch.We do not recommend using this type of inference other than for the purposes of data exploration, even for large values of K (e.g.100).Intuitively, all the information contributing to the estimate of EB is derived from observing the evolution along a single branch at a single site (i.e. from a sample with size &1).To quantify this supposition, we simulated sequence data using the vertebrate rhodopsin phylogeny and branch lengths, applied positive selection of varying strength to five branches in the tree selected a priori (see Text S1), and applied the EB procedure to infer the identity of selected branches.

Model assessment
To assess the performance of MEME on both simulated and empirical data, we selected the fixed effects likelihood method (FEL [24]) as the most appropriate reference test for pervasive diversifying selection, because FEL most closely matches the assumptions made by MEME (see Text S1).We simulated data sets under a number of scenarios: refer to Text S1 for details of simulation strategies.
Assessing the rates of false positives.Under the scenario where each site was evolved under the worst-case null hypothesis of constant v~1, MEME had well controlled rates of false positives at test p-value of 0:05 (Figure S1, also see Text S1 for the empirical derivation of the asymptotic distribution of the test statistic for this hypothesis).MEME appears to be conservative for smaller sample sizes (numbers of sequences, N), but not for larger samples.The rates of false positives were v0:01 (N~8), 0:01 (N~16), 0:03 (N~32), 0:04 (N~64), and 0:05 (N~128 and 256).We also analyzed simulations based on seven large (N~517{640) phylogenies downloaded from TreeBase (http:// www.treebase.org).The rate of false positives remained well controlled (0:047{0:053) at a nominal p-value of 0.05, suggesting that further increasing the number of taxa does not lead to a degradation of Type I error rates.
A further analysis using 36 trees from a variety of published studies downloaded from TreeBase, to simulate 10 replicates from each tree (see Text S1 and Tables S1 and S2 for details), revealed that MEME is generally conservative for alignments of with low pairwise divergence (e.g.v0:1 nucleotide substitutions per site), nominal for those with medium to high pairwise divergence (0:1{0:4 nucleotide substitutions per site), and nominal to slightly anti-conservative for higher pairwise divergence (w0:4 nucleotide substitutions per site), although this relationship is influenced by other factors.Overall, we conclude that false positive rates of MEME, are well controlled in the setting of the most pessimistic (strict neutral) null.
Constant selection pressure at individual sites.At nominal p~0:05 MEME consistently tracked FEL on sequence alignments simulated under the lineage-constant model assumed by FEL (Table S3), losing several percentage points of power because of its more conservative test statistic.Because each simulated alignment contained a subset of sites generated under the null (neutral model), we could derive empirical estimates of the size of the test and set the nominal p-value to achieve a Type I error rate of 5%.When calibrated to deliver a 5% Type I error rate, MEME held a small edge in power.This finding is not surprising, because at a fixed Type I rate, MEME should find every site found by FEL, and resolve FEL borderline cases affected by stochastic variation in v throughout the tree.
Variable selection pressure at individual sites.The difference in power between MEME and FEL became stark when selection at individual sites varied among lineages, with each branch evolving under positive selection (v z ) with probability q z , and negative selection (v { ) with complimentary probability 1{q z .For every combination of independent simulation parameters (v { ,v z ,q z ), MEME had more power to detect sites under episodic diversifying selection (Table 1).Both methods gained power with an increasing proportion of positively selected lineages and/or a greater degree of diversification.The largest differences between MEME and FEL were observed when a small proportion of lineages (q z ~0:1) were subjected to diversifying selection.Regardless of the strength of background purifying selection, FEL was effectively powerless (power 0{10%) to detect episodes of positive selection under any of the three phylogenetic simulation scenarios, whereas MEME achieved low (4{53% when v z ~4), modest (15{95% when v z ~12), and excellent (37{100% when v z ~36) power.Under these conditions, the power of MEME increased with the alignment size, whereas the power of FEL remained very low.Although FEL gained appreciable power when 25% (or 50%) of the lineages were subject to diversification, its power was on average only &24% (&67%) of that realized by MEME.
Taken together, the constant and variable selection pressure simulations demonstrate the uniform superiority of MEME over a standard test for diversifying positive selection.MEME has well controlled rates of false positives, has power comparable to FEL when selective forces are uniform at individual sites, and gains a large power advantage when these forces are variable, as is undoubtedly the case in most biological data sets.
Power and accuracy of the empirical Bayes procedure to identify branches subject to diversifying selection at a single site.Our exploratory simulations (see Figure S2) suggest that it is difficult to accurately identify individual positively selected branches at an individual site.We restricted the analysis to only those sites, which were found to be under episodic diversifying selection by MEME (pƒ0:05) and set the threshold of 20 for the empirical Bayes factor to call an individual branch selected.The best results are achieved when selected branches are placed in the background of strongly conserved lineages (v~0:1) -an individual branch is correctly detected in approximately 25% of cases, while at least one selected branch is found in 89:8% of cases (see Figure S3).However, while none of the negatively selected background branches are reported in more than 5% of cases, in 55% of cases at least one background branch was falsely detected as positively selected.In a more difficult case of neutrally evolving background, the EB procedure performs considerably worse: at least one select branch is found in 55:6% of cases, whereas at least one background branch is detected in 86:5% instances.18 background neutral branches are reported as selected at over 5% frequency, while the 5 positively selected branches are identified at 3:4{26% of selected sites.

Empirical data
To gauge the comparative performance of MEME and FEL when identifying sites subject to pervasive diversifying selection, we used a collection of 16 protein-coding alignments, representing a diverse array of taxa, genes subject to differing levels of conservation, and a range of data set sizes (Table 2).In 12=16 alignments analyzed, MEME identified all the sites inferred by FEL to be under diversifying positive selection and found between Table 1.Comparative performance of FEL and MEME on simulated data where v varies along phylogenetic lineages.

Japanese encephalitis virus env
Vertebrate rhodopsin Camelid VHH  1 (e.g.West Nile virus NS3) and 48 (Diatom SIT) additional sites that were subject to episodic diversifying selection (Table 2).In four data sets, 1{7 sites identified by FEL with p-values close to 0:05 were missed by MEME.Note that MEME p-values for these sites remained in the 0:05{0:07 range (Table 2), i.e. marginally significant.
Sites identified by both methods tended to have a greater average proportion of lineages under selection (0:59, measured by the mean of MLE estimates of q z ); sites found only by MEME experienced more episodic selection (0:10).In 9 data sets (Table 2), sites that FEL inferred to be under purifying selection are instead identified by MEME as likely to have been subjected to episodic diversifying selection.Almost universally (Tables S4, S5, S6, S7, S8, S9, S10, S11, S12, S13, S14, S15, S16, S17, S18, S19), such sites had a smaller estimated proportion of positively selected lineages (v10%).This behavior is consistent with the relative performance of the two tests on simulated data and corroborates the expectation that MEME has greater power to identify sites when only a proportion of lineages evolved under positive selection.Vertebrate rhodopsin, Japanese encephalitis virus env, and Camelid VHH are investigated in detail below; for a discussion other genes, see Text S1.

Vertebrate rhodopsin
The vertebrate rhodopsin (a low-light vision protein) data set was previously experimentally investigated for the substitutions that modulate the wavelength of the light absorbed by the molecule (l max , [18]).The authors asserted that, because none of the 12 sites that they had determined as affecting l max by sitedirected mutagenesis were detected by site-level computational methods, ''statistical tests of positive selection can be misleading without experimental support.''Other authors reanalyzed the same data set more comprehensively and went even further, questioning the utility of v-based methods for detecting experimentally validated sites, because ''most of the current statistical methods are designed to identify codon sites with high v values, which may not have anything to do with functional changes.The codon sites showing functional changes generally do not show a high v value'' [29].The validity of this generalization has been correctly questioned with a simple counter-argument that the sites detected by computational methods may also be functionally important, because the change in l max is unlikely to be the sole determinant of adaptation [17].
The MEME analysis of this gene suggests another obvious alternative, also expounded by previous studies [17]: the failure of the original computational analysis [18] to identify functionally important sites results from the fact that these sites have been subjected to episodic selection, which is masked by predominantly purifying selection elsewhere in the tree.Indeed, among three sites that alter l max found by MEME (96, 183 and 195, versus none found by FEL), no more than 13% of the branches exhibited vw1 (Table S17); at these sites, the average v is less than 1.We note that, because adaptive evolution will not always adhere to a single, simple scenario of episodic diversifying selection, we do not expect MEME to find all 12 sites experimentally confirmed to alter l max .For example, three of the nine missed sites (83,194,292) appear to have been subjected to partial selective sweeps and have been detected using a specialized model of directional evolution [29].
Three sites from this alignment can be used to illustrate how the inclusion of lineage variability modifies inference of selection (Figure 2).Site 54 was inferred to have experienced pervasive nonsynonymous substitutions throughout its evolutionary history.Both FEL and MEME detect this site as positively selected (p~0:02).Sixty three percent of the lineages at this site evolved with b z wa, whereas the remainder were conserved (a~b { ~0), according to MEME.The log-likelihood of the site is only marginally higher for MEME, which suggests that MEME behaves like FEL at sites with Figure 2. Individual sites of the vertebrate rhodopsin alignment used to illustrate similarities and differences between FEL and MEME.Branches that have experienced substitutions, based on most likely joint maximum likelihood ancestral reconstructions at a given site, are labeled as count of synonymous substitutions:count of non-synonymous substitutions.The thickness of each branch is proportional to the minimal number of single nucleotide substitutions mapped to the branch.Branches are colored according to the magnitude of the empirical Bayes factor (EBF) for the event of positive selection: red -evidence for positive selection, teal -evidence for neutral evolution or negative selection, black -E ˆno information.See Methods for more detail.All three sites were identified as experiencing positive diversifying selection by MEME.FEL reported site 54 as positively selected, site 273 as neutral, and site 210 as negatively selected.doi:10.1371/journal.pgen.1002764.g002 ''canonical'' patterns of diversifying selection, corroborating the simulation results.
At codon 273, FEL obtained a maximum likelihood estimate of bwa, but failed to infer positive selection, as the signal was not statistically significant (p~0:70).MEME, on the other hand, allocated 0:04 (0.013-0.10: 95% confidence interval obtained by latin hypercube sampling importance resampling [30]) of branches to a rate class with a~0:0,b z ~9:49 (2.94-6726) and inferred positive selection (p~0:03).The difference in log-likelihoods between MEME and FEL is 4:9 points: MEME fits significantly better, based on a 2-degrees of freedom likelihood ratio test (p~0:007).The maximum likelihood estimates of individual model parameters have large associated errors (although in all posterior samples we obtained b z wa), as is expected for inference based on a single site.This has also been noted by Yang and dos Reis [17].The point estimates themselves, however, are immaterial for inferring whether or not a site is positively selected, since the likelihood ratio test is used for that purpose.Perhaps the most dramatic example of the added power of MEME is illustrated by site 210.At this site, the evolutionary history is replete with non-synonymous substitutions along deep lineages followed by extensive synonymous evolution, indicative of purifying selection.There is also a small clade with repeated synonymous and nonsynonymous substitutions.Averaging over all branches, FEL determined that the site, overall, is under negative selection (p~0:01).MEME reported that 89% of the branches were under a very strong selective constraint (a~2:13,b { ~0:0), but that the remaining 11% 5:5{18:6% ð Þwere under strong diversifying selection (b z ~26:5 10:1{6519 ð Þ ).The log-likelihood improvement is now 13:4 at the cost of two parameters, which is highly significant (pv0:001).Site 210 is the ideal illustration of why it is undesirable to average v over all lineages: bursts of diversification followed by conservation will most likely be missed by traditional site methods.

Japanese encephalitis virus env
No evidence for selection was found in this envelope gene in previous analyses [28], and FEL found only one site under positive selection.Despite the low levels of divergence among a relatively small number of taxa (23 isolates), MEME found episodic selection at sites called negatively selected by FEL (Table S12).Two of these sites fall within a beta-barrel epitope known to be involved in escape from neutralizing antibodies [31].Sites 33 and 242 showed evidence of repeated toggling at terminal lineages.Remarkably, site 33 -likely a part of a neutralizing antibody epitope [32] changed from isoleucine to leucine on 6 terminal lineages; site 242 changed from phenylalanine to serine on 5 terminal lineages.These substitutions co-occur on three terminal lineages.Evidence of recombination was detected in this alignment, and corrected for using a partitioning approach (details on how MEME can correct for recombination are in Text S1).

Camelid VHH
The camelid VHH data set comprises partial variable domain sequences (germline alleles) of llama and dromedary heavy chain only antibodies (Table S3).11 of 16 sites in the variable complementarity determining regions (CDR) 1 (sites 26-33) and 2 (sites 51-58) were found to be under diversifying selection by MEME (2/16 were detected by FEL and 2 more were marginally significant).Because CDR regions are driven to diversify in order to provide a broad basis of antigen recognition, positive selection is expected to be commonplace in the CDRs [33].MEME was able to uncover selective signatures at a majority of those sites.Of the remaining 19 sites classified by MEME as positively selected, six were associated with VHH family differentiation [34].Unlike standard antibodies, which must maintain relatively conserved framework regions (FR) involved in binding heavy and light chains to form functional tetramers, VHH antibodies are free of such functional constraints.A previous analysis of camelid VHH for evidence of positive selection using counting methods [35] reported evidence for positive selection at a single site (14) in FR1 (sites 1-25 in Table S3), but this analysis could find no clear evidence of positive or negative selection at 49 FR sites.In contrast, MEME inferred episodic selection at six sites in FR1, six sites in FR2 (sites 34-50), and 7 sites in FR3 (sites 59{96).The well-known lack of power of counting methods to detect even pervasive selection [17] likely hampered the previous study.

Effect of sequence sampling
Although a previous analysis of 38 vertebrate rhodopsin sequences found no sites under selection at posterior probability §95% [18], the same authors found 7 selected sites in the subset of 11 squirrelfish sequences, and 2 selected sites when the subset of 28 fish sequences was analyzed.These results run counter to the expectation that more data should provide greater power to detect selection.MEME, on the other hand, detects more selected sites when more sequences are included.One site is identified in the squirrelfish alignment, 9 in the fish alignment, and 19 in the complete rhodopsin alignment.All but 5 sites detected in the subset alignments are also identified in the full alignment (Table S20).Allowing v to vary over branches at least partially mitigates the pathology of constant-v models which effectively rely on an average v for inferring selection.A similar pattern is seen in the analysis of the influenza A virus H3N2 hemagglutinin sequences, where site-level methods also appear to be sensitive to sequence sampling ( [19], see Text S1 and Table 23).

Discussion
We have presented a mixed effects model of evolution, MEME, and a statistical test for detecting the signal of past episodic positive selection from molecular sequence data.Our model corrects the biologically unrealistic assumption that selective pressure, as measured by the v ratio, remains constant over lineages.Based on comprehensive simulations and empirical analysis of an array of taxonomically diverse genes, MEME can be recommended as a replacement for existing site models.MEME matches the performance of older approaches when natural selection is pervasive, but possesses greater power to identify sites where episodes of positive selection are confined to a small subset of branches in a phylogenetic tree.
Our results suggest that it may be necessary to revise previous estimates of the proportion of sites under positive selection in many genes.Using the FEL method, which assumes constant selective pressure at a site, we are able to detect 63 sites across all 16 empirical alignments.MEME identifies 51 of these sites (the remaining 12 are borderline significant) and 186 additional sitesnearly 4 times as many as FEL.For individual data sets (e.g.Drosophila adh and Diatom SIT, Table 2), the differences may be even more dramatic.The greater power of MEME indicates that selection acting at individual sites is considerably more widespread than constant v models would suggest.It also suggests that natural selection is predominantly episodic, with transient periods of adaptive evolution masked by the prevalence of purifying or neutral selection on other branches.We emphasize that MEME is not just a quantitative improvement over existing models: for 56 sites in our empirical analyses, we obtain qualitatively different conclusions.FEL asserts that these sites evolved under significant purifying selection, but MEME is able to identify the signature of positive selection on some branches.Furthermore, MEME is less sensitive to sampling effects that plague existing positive selection detection tools [18,19].Variable levels of purifying selection pressure across different lineages prevented these older methods from detecting instances of episodic positive selection; MEME is able to peer through the fog of purifying selection.
It is important to bear in mind that the mixture x 2 statistic used to calculate the p-values reported here is based on a null model under which all sites are evolving neutrally.This, however, is not biologically realistic: the null hypothesis against which sites ideally ought to be screened is one under which sites are evolving either neutrally or under purifying selection.But the proportion of sites evolving under negative selection and the strength of this selection are unknown and vary from case to case, which means that such a null hypothesis would be very sensitive to modeling assumptions that cannot be justified in general.Instead, the neutral null hypothesis represents a worst case scenario for our inference, so that the p-values we obtain are upper bounds of the true p-values.This ensures that our inference is conservative.Even in the worst (and biologically unrealistic) case for MEME, namely when selective pressures are constant throughout the phylogeny, the loss of power compared to FEL is minimal: a site with FEL pvalues between 0:0346 and 0:05 will be missed by MEME, since its p-value will be w0:05 for the same ranges of the likelihood ratio test statistic (LRT).In our simulation scenarios under the assumption of constant v, this translates to no more a 5% loss in power (Table S3).
Our inference is performed in a site-wise rather than an alignment-wide manner, and we therefore control the site-wise rather than the family-wise error rate.We do not recommend combining the results of multiple site-wise inferences to perform alignment-wide inference.To aid interpretation of the results while taking account of multiple testing, we calculate the false discovery rate [36]; the resulting q-value upper bounds are reported alongside their corresponding p-value upper bounds in Tables S4, S5, S6, S7, S8, S9, S10, S11, S12, S13, S14, S15, S16, S17, S18, S19.This gives an upper bound on how many of the reported sites can be expected to be false discoveries: for instance, of the 30 sites reported in Table S5 we expect no more than 5 (14%) to be false, and probably far fewer because of the conservative choice of null model.We emphasize that q-values are usually much larger than their corresponding p-values and caution that p-values (regardless of whether they have been corrected for multiple testing) cannot be used to estimate an expected number of false discoveries in the same way.
MEME is a conceptual advance over the first generation of random effects models designed to detect episodic selection (called ''branch-site models'' in the literature [17]).MEME does not require a priori designation of, or an exhaustive search for, the branches under selection, and it allows each site to have its own selective history.Whereas branch-site models make restrictive a priori assumptions about how v values are distributed across the tree -sometimes leading to very poor statistical performance [20] -MEME treats the selective class on each branch as a random effect that is marginalized over in the likelihood calculation.
For computational tractability, MEME assumes that the value taken by v on each branch is independent of that on any other branch, i.e. selective pressures between branches are uncorrelated.This assumption could potentially be violated: for example, if v changes very slowly across the phylogeny, then v values on neighboring branches will be correlated.Further research is needed to understand how inference of selection would be affected if these correlations were directly accounted for, and whether the additional model and computational complexity would be justified.In practice, MEME could be combined with models of directional selection to improve power, e.g.[15,16].Unlike covarion models [37,13], MEME does not allow v to change in the middle of a tree branch.The effect of this restriction is unclear, but it could be tested by implementing a mixed effects covarion model, where switching rates and proportion of time spent under vw1 are estimated at an individual site.
The ability of MEME, or similar substitution model-based methods, to accurately infer the identity of individual branches subject to diversifying selection at a given site seems unavoidably limited.Most of the information that such inference might be based on is limited to character substitutions along a single branch at a single site, i.e. one realization of the Markov substitution process.Selection along terminal branches in the context of negatively selected background can be detected more reliably than selection along interior branches among neutrally evolving background lineages.However, we caution that despite obvious interest in identifying specific branch-site combinations subject to diversifying selection, such inference is based on very limited data (the evolution of one codon along one branch), and cannot be recommended for purposes other than data exploration and result visualization.This observation could be codified as the ''selection inference uncertainty principle'' -one cannot simultaneously infer both the site and the branch subject to diversifying selection.In this manuscript, we describe how to infer the location of sites, pooling information over branches; previously [20] we have outlined a complementary approach to find selected branches by pooling information over sites.
Finally, although MEME is considerably more powerful than existing methods at detecting bursts of selection, it still requires that a measurable proportion of lineages (5{10%) experience non-synonymous evolution at a site.When a single substitution modifies an adaptive trait and is subsequently fixed, we expect v based methods to have very little power.Specialized methods which make use of change in allele frequencies [15,16], or between and within-population diversification patterns [38], will be required in such cases.

(PDF)
Table S20 Test p-values for positively selected sites found by MEME in a set of 38 vertebrate rhodopsin sequences analyzed with REL methods in Yokoyama2008fk.Sites with pƒ0:05 are shown in bold.The partial ordering of subsets is as follows: Squirrelfish 5 Fish 5 All, Coelacanth and tetrapods 5 All.Sites found to be under positive selection with posterior probability of w95% (M8 model) in Yokoyama2008fk in at least one of the subsets are marked with ? .(PDF) Table S21 Test p-values for positively selected sites found by MEME in a set of 86 influenza A virus hemagglutinin sequences (Set 3) and its various subsets, analyzed with REL methods in Chen2011fk.Sites with pƒ0:05 are shown in bold.The partial ordering of subsets is as follows: Set Text S1 Supplementary methods, results, and discussion.(PDF)

N
(S) reports the number of sequences (codons) in the alignment.M z (M { ) refers sites found by MEME to be positively (negatively) selected (pƒ0:05).F z (F { ) denote sites found by FEL to be positively (negatively) selected (pƒ0:05).F 0 references sites that are classified as neutrally evolving by FEL.Values in parentheses for the M { F z column show the mean p-values for FEL and MEME on this set of sites, respectively.Values reported in the rightmost column count the number of sites where MEME fits significantly better than FEL, based on a 2-degrees of freedom LRT (pƒ0:05).Abbreviations: IAV = Influenza A virus, JEV = Japanese encephalitis virus.doi:10.1371/journal.pgen.1002764.t002

4 5 1 5 5 2 5
Set Set 3, Set 5 Set Set 3, Set 6 5 Set 3, Set 7 5 Set 3. Sites found to be under positive selection with posterior probability of w95% (M3 model) in Chen2011fk in at least one of the subsets are marked with ? .(PDF)

Table 2 .
Power to detect sites under selection (p~0:05) are reported for FEL and MEME (in boldface) for each unique combination of negative selection (v { ), positive selection (v z ), and proportion of branches under positive selection (q z ) parameters.doi:10.1371/journal.pgen.1002764.t001Comparative performance of MEME and FEL on 16 empirical alignments (see Results and Text S1 for an extended discussion of each individual case).