text:
Alternative Substrate Kinetics of Escherichia coli Ribonuclease P

Background: RNase P acts on multiple tRNA precursors, but little is known about substrate competition. Results: Alternative substrate kinetics reveal a simple competitive mechanism allowing relative rate constants to be quantified. Conclusion: The relative rates of processing of different substrates are restricted to a narrow range and reflect association. Significance: RNase P processing is tuned for uniform specificity and coupling to tRNA precursor biosynthesis. A single enzyme, ribonuclease P (RNase P), processes the 5′ ends of tRNA precursors (ptRNA) in cells and organelles that carry out tRNA biosynthesis. This substrate population includes over 80 different competing ptRNAs in Escherichia coli. Although the reaction kinetics and molecular recognition of a few individual model substrates of bacterial RNase P have been well described, the competitive substrate kinetics of the enzyme are comparatively unexplored. To understand the factors that determine how different ptRNA substrates compete for processing by E. coli RNase P, we compared the steady state reaction kinetics of two ptRNAs that differ at sequences that are contacted by the enzyme. For both ptRNAs, substrate cleavage is fast relative to dissociation. As a consequence, V/K, the rate constant for the reaction at limiting substrate concentrations, reflects the substrate association step for both ptRNAs. Reactions containing two or more ptRNAs follow simple competitive alternative substrate kinetics in which the relative rates of processing are determined by ptRNA concentration and their V/K. The relative V/K values for eight different ptRNAs, which were selected to represent the range of structure variation at sites contacted by RNase P, were determined by internal competition in reactions in which all eight substrates were present simultaneously. The results reveal a relatively narrow range of V/K values, suggesting that rates of ptRNA processing by RNase P are tuned for uniform specificity and consequently optimal coupling to precursor biosynthesis.

Ribonuclease P (RNase P) 2 is an essential ribonucleoprotein enzyme that is responsible for catalyzing the maturation of the 5Ј end of transfer RNAs (tRNAs) through site-specific hydrolysis of a phosphodiester bond in precursor tRNAs (ptRNAs) (1)(2)(3). The RNA subunit, termed P RNA, contains the active site (4,5), whereas the smaller protein subunit (C5 in Escherichia coli) is required for optimal molecular recognition and catalysis in vitro and is essential in vivo (6 -13). Although P RNA is a ribozyme, its mode of molecular recognition differs from other catalytic RNAs in two important ways. First, its biological role in ptRNA processing requires that it act in trans as a multiple turnover enzyme (14), whereas, other ribozymes, with the exceptions of the ribosome and spliceosome, undergo single turnover self-splicing or self-cleavage reactions (15)(16)(17). Second, RNase P processes multiple RNA substrates, including all ptRNAs in the cell, whereas other ribozymes, again with the exceptions of the ribosome and spliceosome, have one specific substrate (18,19). These characteristics are essential to RNase P function, as they are to the ribosome and spliceosome, and are common to many enzymes involved in RNA metabolism (20 -23). Therefore, understanding the multiple substrate recognition properties of RNase P can shed light on general principles of molecular recognition by other ribonucleoproteins and multisubstrate enzymes.
The ptRNA nucleotides contacted by RNase P have been determined by chemical interference and structure-function studies ( Fig. 1) (14,18,24,25). The recognition elements near the cleavage site include the 3Ј-RCCA sequence, a G(ϩ1)/ C(ϩ72), as the first base pair in the acceptor stem and the 2Ј-OH and nucleobase of a U(Ϫ1) residue 5Ј to the cleavage site. The substrate binding domain of P RNA also contacts 2Ј-OH groups in the T stem loop (26,27). The spacing of these contacts in the T stem loop in relation to the cleavage site results in an overall shape recognition of the substrate (28 -31).
Comparative analysis of E. coli tRNA gene sequences shows significant variation among the nucleotides identified as contact points with the enzyme. As shown in Fig. 1, alignment of the 87 ptRNA genes of E. coli K12 (19) reveals that the leader sequences (region I) and the acceptor stem (regions II and III) show only minimal sequence conservation. An exception is the 3Ј-CCA sequence that is recognized by the ribosome (32), aminoacyl-tRNA synthetases (33,34), and elongation factor Tu (35). Only two-thirds of E. coli ptRNAs (66/87) contain a G(ϩ1)/C(ϩ72), and a similar fraction (63/87) have an optimal U at the N(Ϫ1) position (19). The population of ptRNAs that contain all of the recognition elements is significantly smaller (ϳ50%; 42/87). These ptRNAs make up a canonical sequence, whereas a noncanonical ptRNA is missing one or more of these recognition elements. The adjacent base pair to these recognition elements is often a G(ϩ2)/C(ϩ71); however, this position is not known to contact RNase P. The 5Ј leader sequence shows no conserved motif; however, both binding and cleavage of model substrates by E. coli RNase P are sensitive to changes in the sequence of the 5Ј leader (36,37). Indeed, recent studies identified a protein-RNA interaction between the leader sequence and the Bacillus subtilis RNase P (38). A structure from Mondragón and colleagues (31) of the Thermotoga maritima RNase P bound to tRNA and leader products is consistent with the experimentally defined interface between enzyme and substrate drawn from biochemical studies. Although specific leader contacts are not resolved in the crystal structure, it generally corresponds with the perspective from cross-linking and structure-function studies (10,36).
Central to achieving a complete understanding of multiple substrate recognition by RNase P is the observation that catalysis by P RNA alone is sensitive to natural structural variation among ptRNAs that results in loss of RNA-RNA contacts between P RNA and ptRNA (36, 39 -41). Catalysis by the ribonucleoprotein holoenzyme, which forms additional leader sequence interactions, is less sensitive to sequence and structure variation among endogenous ptRNAs (36). A conformational change during substrate recognition has been documented for B. subtilis RNase P by Fierke and colleagues (9,25) in which the protein subunit facilitates via leader sequence contacts. This two-step mechanism for substrate binding may give rise to threshold effects resulting in similar rate constants for catalysis for substrates lacking optimal contacts with the enzyme (36,39,40,42).
Thus, detailed in vitro structure-function studies measuring binding and catalysis for model substrates have revealed basic principles of molecular recognition by RNase P. Nonetheless, information on the competition between different alternative substrates is needed to understand function in vivo. Here, we test a simple competitive model to describe the relative rates of ptRNA processing by RNase P and apply this model to evaluate the effect of natural, genomic variation in ptRNA substrate structure on relative processing rates. The results provide insight into the features of the kinetic mechanism of RNase P that govern its function in vivo that are relevant for other multiple substrate enzymes.

EXPERIMENTAL PROCEDURES
RNA Synthesis and RNase P Holoenzyme Reconstitution-E. coli P RNA and ptRNAs were generated by in vitro transcription from cloned DNA or PCR DNA templates. The ptRNA fMet47 and ptRNA Met82 substrates were synthesized from plasmids ptRNA605 and ptRNA608, respectively, as described previously (39). The ptRNA His37 , ptRNA Gln85 , ptRNA Gly62 , ptRNA Ser80 , ptRNA Leu76 , and ptRNA Ile1 substrates were synthesized from PCR products generated from E. coli genomic DNA using primers that introduce a T7 promoter and nine nucleotides of 5Ј leader sequence to generate the RNA sequences shown in Fig. 2: ptRNA Gln85 , forward, 5Ј-TAATACGACTCACTATAGGCCGGTTATTGGGGTA-TCGCC-3Ј, and reverse, 5Ј-TGGCTGGGGTACCTGGATTC-G-3Ј; ptRNA Gly62 , forward, 5Ј-TAATACGACTCACTATAG-GATCTCGAAGCGGGCGTAGTTC-3Ј, and reverse, 5Ј-TGG-AGCGGGCGAAGGGAATCG-3Ј; ptRNA Ser80 , forward, 5Ј-TAATACGACTCACTATAGGGTCATTCCGGAAGTGTG-GCCG-3Ј, and reverse, 5Ј-TGGCGGAAGCGCAGAGATTCG-3Ј; ptRNA Ile1 , forward, 5Ј-TAATACGACTCACTATA-GGACCTCTACAGGCTTGTAGCTC-3Ј, and reverse, 5Ј-TGGTAGGCCTGAGTGGACTTG-3Ј. PCR reactions contained 100 mM KCl, Tris-HCl, pH 8, 5 mM MgCl 2 , 100 M dNTP and were performed at 95°C, 30 s, 60°C, 60 s, and 72°C, 90 s for 30 cycles. The resulting DNA was purified by phenol and chlo-FIGURE 1. Secondary structure and sequence conservation of E. coli (K12) ptRNAs at regions contacted by RNase P. A, the conserved secondary structure of tRNA is shown. The RNase P ribonucleoprotein interacts with the acceptor stem and TC stem and loop of tRNA (black circles). The enzyme also contacts in the 5Ј leader (gray circles). The sequences identified as forming RNA-RNA contacts with the RNase P enzyme are shown as letters and include a U at N(Ϫ1), a G(1)-C(73) base pair at the top of the acceptor stem, and the 3Ј-RCCA sequence as described in the Introduction. B, the linear sequences of tRNA involved in substrate binding are separated into three regions for presentation of sequence variation. Genomic tRNA sequences and alignments were obtained from the genomic tRNA database (19). Sequence logos for regions I, II, and III were created using WebLogo (88). Region I (N(Ϫ10 to N(Ϫ1)) includes the protein binding site in the 5Ј leader and the nucleotide at N(Ϫ1) that contacts P RNA. Region II (N(1) to N(7)) is the 5Ј side of the acceptor step including the RNase P cleavage site 5Ј to N(1). Region III (N (66) to N(76)) includes the 3Ј side of the acceptor stem and the conserved RCCA sequence that interacts with P RNA. roform extraction and precipitated with ethanol as described previously (40).
In vitro transcription reactions contained ϳ10 g of plasmid or PCR template, 5 mM each NTP, 50 mM Tris-HCl, pH 8, 17.5 mM MgCl 2 , 10 mM DTT, 2 mM spermidine, and 1 unit/l of T7 RNA polymerase (Ambion) and were incubated overnight at 37°C. The resulting RNA products were recovered by phenol/ chloroform extraction and ethanol precipitation. The fulllength transcription products were purified by PAGE as described (40).
Kinetic Analyses of Single ptRNA Substrate Reactions-The ptRNA substrates were 5Ј-end-labeled with [␥-32 P]ATP and T4 polynucleotide kinase after dephosphorylation by alkaline phosphatase. RNase P holoenzyme reaction kinetics were measured under the following conditions: 50 mM Tris-HCl, pH 8, 100 mM NaCl, 17.5 mM MgCl 2 , and 0.005% Triton X-100. The P RNA and 5Ј-32 P-end-labeled ptRNA were renatured separately in the above buffer (omitting the Mg 2ϩ ) by incubation at 95°C for 4 min followed immediately by incubation at 37°C for 10 min. MgCl 2 was added to a concentration of 17.5 mM, and the incubation was continued for an additional 10 min. C5 protein was purified, and activity was tested by titrations into multiple turnover reactions containing a constant concentration of P RNA as described previously (43). C5 protein was added to a concentration equal to that of P RNA, and the 37°C incubation was continued for an additional 10 min. Reactions were initiated by mixing equal volumes of enzyme and substrate. Aliquots were removed at specific time intervals and quenched with 50 mM EDTA. The residual ptRNA substrate and 5Ј leader cleavage product were resolved by denaturing PAGE (15%). The conversion of substrate to product was quantified by Phos-phorImager analysis using a GE Healthcare Storm system and ImageQuant software. The fraction of reaction was determined as F ϭ For multiple turnover reactions, initial rates were measured and plotted versus ptRNA concentration and fit to the Michaelis-Menten equation as described under "Results." Transient kinetic experiments were performed in a similar manner, with the substrate and enzyme concentrations indicated in the figure legends. Single turnover kinetic analyses were performed as described previously (44) using the following reaction conditions: 50 mM MES, pH 6, 100 mM NaCl, 17.5 mM MgCl 2 , and 0.005% Triton X-100. F was plotted versus time and fit to a single exponential function, where F o is the maximal fraction of reacted substrate. Experiments were measured using the single turnover conditions described above at pH 6 with either 1 nM ptRNA fMet47 or 1 nM ptRNA Met82 and 100 nM RNase P. At 30 s after initiating, the reaction the volume was divided in half. One aliquot was added to a tube containing a substrate quench giving a final concentration of 5 M B. subtilis ptRNA Asp . Both reactions were continued, and the kinetic data for both pre-chase and post-chase were fit to Equation 1 (45).

Determination of Relative Rate Constants by Internal
Competition-To detect the formation of products from two substrates independently in the same reaction, ptRNA fMet47 and ptRNA Met82 were modified by the addition of two extra G nucleotides to the 5Ј end of the leader sequence, giving rise to ptRNA fMet47 (ϩ2) and ptRNA Met82 (ϩ2). The two additional residues allow their products to be separated on 15% PAGE and quantified independently from the products of the other ptRNAs used in this study.
Internal competition reactions contained two, three, or eight ptRNAs as indicated in the text and in the legends for Figs. 6 -9. Reactions containing both ptRNA fMet47 and ptRNA Met82 at a range of concentrations from 10 to 100 nM were performed with one or the other substrates containing the two extra G residues in the 5Ј leader. Additional reactions containing ptRNA fMet47 and ptRNA Met82 (ϩ2) at 100 nM in the presence of increasing concentrations of nonradiolabeled ptRNA Leu76 from 10 to 3000 nM were also analyzed by monitoring the formation of the 5Ј end-labeled ptRNA fMet47 and ptRNA Met82 (ϩ2). Similarly, reactions containing all eight of the substrates shown in Fig. 2 were conducted with each ptRNA present at 100 nM (800 nM total). For these reactions, trace concentrations of radiolabeled reference substrate ptRNA Met82 (ϩ2) and one of the remaining seven substrates were included to follow product formation. Relative rate constants ( r k ϭ (V/K)/(V/K) reference ; see below) were determined using analytical methods based on Equation 2 and Scheme 1, as described in the following sections.

RESULTS
Application of Competitive Alternative Substrate Kinetics to ptRNA Processing by RNase P-As illustrated in Scheme 1, a simple competitive multiple turnover mechanism allows the competition between different ptRNA substrates for processing by RNase P to be quantified (46 -48). A single population of RNase P (E) combines with multiple ptRNA substrates (S 1 , S 2 , S 3 . . . S N ) to form individual RNase P-ptRNA complexes (ES 1 , ES 2 , ES 3 . . . ES N ) that react with rate constants V 1 , V 2 , V 3 . . . V N to form tRNA and leader products that together are represented by P 1 , P 2 , P 3 . . . P N . We apply the convention of using V and V/K as the fundamental multiple turnover kinetic parameters. The parameter V is the rate constant for reaction of ES to form products and regenerate free enzyme and is equivalent to k cat . The V/K is the second order rate constant at limiting substrate concentrations (i.e. k cat /K m ) (49, 50). Importantly, both S 1 and S 2 must compete with the remaining population of sub- SCHEME 1 strates, which act as competitive inhibitors (46 -48, 51). As a result, the expression for the ratio of the rates for conversion of S 1 and S 2 to products simplifies to Thus, the ratio of the observed rates of product formation for the two substrates depends on the ratio of their V/K values and their concentrations. The designation r k is used below to refer to the ratio of the V/K values for an experimental or unknown substrate relative to a reference substrate ( r k ϭ (V/K)/ (V/K) reference ) (51). As indicated under "Experimental Procedures," the ptRNA Met82 (ϩ2) substrate is used as the primary reference in this study.
There are two key consequences of Scheme 1 and consequently Equation 2 that are important in considering the in vivo function of RNase P (46 -48, 51). First, the relative V/K values and consequently the observed rates of any two substrates will be independent of the presence or concentration of alternative substrates. The reason for this is that the additional substrates act as competitive inhibitors decreasing the concentration of free enzyme available for all substrates equally. Second, the relative processing rates will depend on the V/K values of the two substrates regardless of the enzyme concentration or whether either substrate concentration is saturating. These considerations highlight that the second order rate constant at limiting substrate is an essential parameter in understanding the biological function of RNase P, as it is with other enzymes.
Accordingly, we set out to test whether this simple competitive model describes the relative rates of ptRNA processing by RNase P and to evaluate the effect of natural, genomic variation in ptRNA substrate structure on the kinetics of competi-tion reactions containing multiple substrates. As described in the following sections, we first measured the V and V/K values for two well characterized canonical and noncanonical ptRNAs using standard steady state reactions of uniform RNA populations. We used pre-steady state and single turnover kinetic analyses to determine the reaction steps that limit V and V/K. Reactions containing mixtures of both substrates were analyzed, and the simple competitive model described above was validated. Using an internal competition approach based on this model, we determined the relative rate constants for eight different ptRNAs representing the range of ptRNA structural variation at sites of RNase P contact occurring in the E. coli genome.
Comparison of the Multiple Turnover Kinetics of ptRNA Met82 and ptRNA fMet47 Processing by E. coli RNase P-The substrates ptRNA Met82 and ptRNA fMet47 (Fig. 2) were selected as representative examples of canonical and noncanonical ptRNAs, respectively. Both ptRNAs have similar sequence length and base composition; however, they differ significantly at the nucleotides contacted by the P RNA subunit of RNase P. The ptRNA fMet47 , an initiator tRNA, has an A instead of an optimal U at N(Ϫ1) and a C(ϩ1)-A(ϩ72) pair at the cleavage site that results in a Ͼ900-fold decrease in the rate of catalysis by the P RNA subunit alone (36). In contrast, the RNase P holoenzyme binds both ptRNA Met82 and ptRNA fMet47 with equivalent equilibrium binding constants and processes them with similar single turnover rate constants (36). The metal ion and pH dependence of the single turnover reactions of both substrates are also comparable (37). To isolate the effects of tRNA sequence and structure that contact RNase P from secondary effects due to flanking sequences that are idiosyncratic to individual ptRNAs, we used a standard substrate structure containing the tRNA and 10 additional leader sequences (Fig. 2).
To evaluate the V/K for processing of ptRNA Met82 and ptRNA fMet47 by RNase P, the observed initial rates for both substrates were plotted against their concentrations and fit to the Michaelis-Menten equation (Fig. 3).
The steady state kinetic parameters V and K for both substrates are highly similar (V Met82 ϭ 0.11 Ϯ 0.01 s Ϫ1 ; V fMet47 ϭ 0.14 Ϯ 0.01 s Ϫ1 and K Met82 ϭ 310 Ϯ 60 nM; K fMet47 ϭ 280 Ϯ 40 nM), resulting in an r k ratio near unity (ϳ0.9, where r k ϭ (V fMet47 / K fMet47 )/(V Met82 /K Met82 )) ( Table 1). Fitting complete time courses of the multiple turnover reactions of ptRNA Met82 and ptRNA fMet47 to the integrated Michaelis-Menten equation shows evidence of product inhibition (Fig. 4). An additional approach to measure V/K from multiple turnover reactions is to analyze progress curve data using the integrated Michaelis-Menten equation.
Although the multiple turnover time courses for RNase P cleavage of ptRNA Met82 and ptRNA fMet47 fit well to the above equation, the values of V/K determined using the initial rate data do not predict the observed time courses for either substrate (Fig. 4, dotted lines). It is observed that the kinetics are significantly slower and that a much larger K value is obtained from fitting. These features are hallmarks of product inhibition, and thus we fit the progress curve data to the integrated equation including product inhibition.
Equilibrium binding studies as well as competitive single turnover inhibition experiments indicate that the K d for the   tRNA Met82 is 150 nM and that the K d for the tRNA fMet47 is 100 nM (36,39). Using these values for K i in the above equations provides a much improved fit of the data (Fig. 4, solid line). The values of V and K for the two substrates obtained by this method are ϳ2-fold lower than those obtained from analysis of the initial rate data; however, the values of V/K are highly similar.
The multiple turnover kinetic parameters, V and V/K, estimated by both approaches are highly similar for both substrates despite their significant difference in structure (Table  1). Next, we asked whether the similar V and V/K values for the two substrates reflect the same or different rate-limiting steps.
Pre-steady State Kinetic Analyses to Evaluate the Reaction Step That Limits V-The kinetics of ptRNA cleavage at increasing concentrations of RNase P were examined to determine the predominant form of the enzyme that is populated at steady state (ES or EP, in Scheme 1). As shown in Fig. 3B, for reactions in which RNase P (5 and 10 nM) and the ptRNA substrate (500 nM) are both present at concentrations in excess of K d (Ͼ1 nM), there is a linear increase in product concentration that extrapolates back to the origin. Reactions with either 5 nM or 10 nM RNase P result in product formation that increases linearly with no evidence for a pre-steady state burst. A simple interpretation of this result is that the net rate constant for dissociation of products and regeneration of free enzyme (k 3 in Scheme 2) is faster than substrate cleavage (k 2 ). Therefore, ES is the predominate form of the enzyme that accumulates at steady state. This result contrasts with the kinetic mechanism of B. subtilis RNase P, which is limited by product release for a canonical ptRNA Asp substrate (52).
Single Turnover Kinetics to Evaluate the Reaction Step That Is Rate-limiting for V/K-An important observation relevant to the reaction mechanism of E. coli RNase P is that the observed K (where K ϭ (k off ϩ V)/k on in Scheme 1; i.e. K m ) (53) from multiple turnover kinetic analyses is greater than the independently measured equilibrium dissociation constant, K d (310 nM versus 0.5 nM for ptRNA Met82 and 280 nM versus 0.3 nM for ptRNA fMet47 ) (Tables 1 and 2). This result implies that the net rate constant for cleavage to regenerate free enzyme (V ϭ k 2 k 3 /k 2 ϩ k 3 ) in Scheme 2) is fast relative to substrate dissociation (k Ϫ1 ) (53)(54)(55). It follows that at limiting substrate concentration, the rate of multiple turnover could therefore be limited by substrate association (50).
To test these predictions, we determined the relative magnitudes of the rate constant for catalysis, k 2 , and the rate constant for substrate dissociation, k Ϫ1 , using a "sequential mixing" or "isotope trapping" experiment (45). The RNase P-ptRNA complex was formed by mixing limiting substrate (1-2 nM) with a saturating concentration of enzyme (100 nM). At an intermediate time, an excess of nonradiolabeled substrate was added. If k 2 is fast relative to k Ϫ1 , then there will be little dissociation of the remaining RNase P-ptRNA complexes over the remaining time course of the reaction, and correspondingly no effect on the accumulation of product. Alternatively, if substrate dissociation is fast relative to catalysis (k Ϫ1 Ͼ Ͼ k 2 ), then the addition of nonradiolabeled substrate will quench the formation of radiolabeled product.
The dependence of the observed pseudo-first order rate constant on enzyme concentration showed saturable behavior as predicted based on Scheme 1 (data not shown). These data allowed reaction conditions to be determined under which all of the radiolabeled ptRNA is present in the ES complex. As shown in Fig. 5, A and B, the addition of a cold substrate chase after formation of ES did not result in quenching or a change in reaction kinetics. In contrast, the addition of the excess nonradiolabeled substrate at the start of the reaction resulted in the expected slow, multiple turnover kinetics. Therefore, we concluded that substrate dissociation is negligible over the remaining time course of the reaction (k 2 Ͼ Ͼ k Ϫ1 ).
An important implication of this result is that the substrate association rate constant, k 1 , can be measured from the concentration dependence of the single turnover reaction (k obs versus [E]) (54). Fitting the dependence of k obs to [E] at concentrations below K 1/2 permits k 1 and k Ϫ1 to be estimated as the slope and intercept (Fig. 5C). The kinetic parameters for both substrates are similar (1.9 Ϯ 0.1 ϫ 10 6 and 1.5 Ϯ 0.4 ϫ 10 6 M Ϫ1 s Ϫ1 for ptRNA Met82 and ptRNA fMet47 , respectively) ( Table 2). The estimates for k Ϫ1 from second order analyses are less than the V for both substrates. In this case, the observed K for the multiple turnover reaction will be approximated by V/k 1 (53). The experimentally measured values of these kinetic parameters result in calculated K values of 57-172 nM for ptRNA Met82 and 71-117 nM for ptRNA Met47 . These calculated values are within 2-fold of the experimentally observed K determined from analysis of initial rate data (Tables 1 and 2). It is possible that differences in the reaction pH or errors in the determination of concentrations of substrate and enzyme account for this difference.
From the definitions for V and K, above, it follows that V/K ϭ k 2 k 1 /(k Ϫ1 ϩ k 2 ) (50). Thus, when k 2 Ͼ Ͼ k Ϫ1 , then V/K Ϸ k 1 . Therefore, the most simple interpretation of the pre-steady state and single turnover results is that the cleavage step (k 2 ) is rate-limiting for V (i.e. at saturating substrate concentrations) and that V/K reflects the association step (k 1 ) for both ptRNAs at limiting substrate concentrations.   MARCH 22, 2013 • VOLUME 288 • NUMBER 12

JOURNAL OF BIOLOGICAL CHEMISTRY 8347
Competitive Alternative Substrate Kinetics of ptRNA Met82 and ptRNA fMet47 Processing by RNase P-As introduced above, in competitive multiple turnover reactions, the relative rates for two competing ptRNAs are expected to be determined by their relative V/K values ( r k ϭ (V/K)/(V/K) reference ) and their concentrations (46 -48, 51). Also, it follows that the presence of additional substrates will decrease the observed rates for all substrates in the reaction due to competition for free enzyme, but should not affect the r k value for any two substrates (46,56,57). We tested the competitive alternative substrate model for RNase P by analyzing the competitive kinetics of reactions containing both ptRNA Met82 and ptRNA fMet47 .
To simultaneously measure the reaction kinetics of two ptRNAs in the same reaction, we used a reference substrate in which two additional G residues are added to the 5Ј end of the leader sequence. This modification allows the products of the ptRNAs to be distinguished by their mobility on denaturing PAGE and quantified individually. An example of the primary data from this approach for ptRNA Met82 (ϩ2) and ptRNA fMet47 is shown in Fig. 6A. The precursor band contains both substrates as these species are not resolved under these gel conditions. However, the products from the two substrates are readily distinguished and quantified, allowing relative rates of product formation to be measured. To address the effect of the additional nucleotides on ptRNA processing, the relative rate constants for comparison of ptRNA Met82 with ptRNA Met82 (ϩ2) and comparison of ptRNA fMet47 with ptRNA fMet47 (ϩ2) were also measured (data not shown) and were observed to be 0.9 -1.3. Thus, the presence of the additional nucleotides required to distinguish the products from two substrates has essentially no effect on the rate of RNase P processing.
To compare the competitive kinetics of ptRNA fMet47 and ptRNA Met82 , the observed rates of product formation were determined for reactions containing substrate concentration ratios (in nM) of 10:100, 100:10 10:10 and 100:100 (ptRNA fMet47 : ptRNA Met82 ). The product ratios from at least three time points taken under steady state conditions were averaged and then corrected for the relative substrate concentrations. Additionally, the collection of data for the observed rates as a function of the relative rates of the two substrates were fit to the logarithmic form of Equation 2, log(v 2 /v 1 ) ϭ log r k ϩ log͑S 2 /S 1 ͒ (Eq. 6) where v 2 /v 1 is the ratio of the observed initial rates for ptRNA fMet47 relative to ptRNA Met82 (ϩ2). Analysis of the data in this manner allows determination of r k from the combined data set. As shown in Fig. 6B, the data for both the ptRNA fMet47 / ptRNA Met82 (ϩ2) and the ptRNA fMet47 (ϩ2)/ptRNA Met82 combination of substrates fit this relationship as predicted. Fitting to Equation 4 yields an r k value of 0.5 ((V/K) fMet47 /(V/K) Met82 ) for the reaction in which the ptRNA Met82 was modified to contain the additional two leader nucleotides. As a control, the r k was measured in competitive reactions in which ptRNA fMet47 instead of ptRNA Met82 was lengthened to distinguish the products from the two substrates. A similar value of 0.6 was observed consistent with the value measured in which the ptRNA Met82 (ϩ2) was used for the reference substrate.
An additional prediction of the internal competition model is that the addition of a third substrate will not affect the r k value for these two substrates. Accordingly, we tested the effect of increasing concentrations of a third substrate on the observed rates of ptRNA Met82 (ϩ2) and ptRNA fMet47 product formation. In Scheme 1, additional substrates act as competitive inhibitors that decrease the observed rate of processing of both substrates by competing for free enzyme. In Fig. 6C, ptRNA Leu76 is added as a competitive alternative third substrate in reactions containing ptRNA fMet47 and ptRNA Met82 (ϩ2) as the reference substrate. Increasing concentrations of nonradiolabeled ptRNA Leu76 , which binds to RNase P with similar affinity as the FIGURE 5. Single turnover kinetics of ptRNA Met82 and ptRNA fMet47 processing by RNase P. A, single turnover sequential mixing experiment with initial concentrations of 1 nM ptRNA Met82 and 100 nM RNase P. At the time indicated by the vertical dotted line, the reaction was divided, and one fraction was combined with a high concentration (5 M) of nonradiolabeled ptRNA (circles). Time points were continuously collected from the remaining fraction (squares). As a control, an identical reaction was combined with nonradiolabeled substrate before the addition of enzyme (triangles). B, single turnover sequential mixing experiment using ptRNA fMet47 performed as described in panel A. C, second order analysis of RNase P binding of ptRNA Met82 and ptRNA fMet47 to increasing concentrations of RNase P. The pseudo-first order rate constants (k obs ) determined for a single turnover reaction containing 2.5-10 nM RNase P concentrations are plotted versus [E]. These data are fit to a linear function k obs ϭ k 1 [E] ϩ k Ϫ1 to determine the rate constants reported in Table 2. Error bars indicate S.D.
other two ptRNAs in the reaction, decrease the observed rates of ptRNA fMet47 and ptRNA Met82 (ϩ2) processing as expected. The data fit to a simple competitive inhibition model derived from Scheme 1, where S 1 is the concentration of the labeled substrate and S 2 and S 3 are the concentrations of the competitive alternative, unlabeled substrates in Fig. 6C. Analysis of the observed rates data for ptRNA Met82 and ptRNA fMet47 in the presence of 10 -3000 nM ptRNA Leu76 allows the K value for ptRNA Leu76 to be estimated. A value of ϳ300 nM is obtained, which is similar to the values measured by analysis of initial rate data for reactions containing ptRNA Met82 or ptRNA fMet47 alone (Table 1). Nonetheless, as demonstrated in Fig. 6C, inset, the ratio of the observed rates, the r k for ptRNA fMet47 referenced to ptRNA Met82 (ϩ2), is independent of the presence and concentration of a competing substrate. Because the r k value for two competing substrates is insensitive to a third competing substrate, the internal competition approach could be used to determine the r k values for substrates in reactions containing more complex populations.

Determination of Relative Rate Constants for ptRNAs in Complex Substrate Populations by Internal Competition-It follows
from Scheme 1 and the observations documented above that the presence of additional substrates, regardless of their number or concentration, should also have no effect on the relative rates of processing of any two substrates in the population. To test this concept, we generated five ptRNA substrates in addition to ptRNA Met82 , ptRNA fMet47 , and ptRNA Leu76 . Substrates were selected to span the range of ptRNA structure variation encountered by E. coli RNase P in vivo, and their secondary structures are shown in Fig. 2. Among these, similar ptRNA His and ptRNA Ser substrates have served as substrates for analyzing the determinants of specificity adjacent to the site of 5Ј processing (58).
We used the same approach, described in the preceding section, of distinguishing between the products of two substrates by analyzing the relative rate constants of ptRNA Met82 (ϩ2) and ptRNA fMet47 . Because the two substrates of interest are the only species that are radiolabeled, their products alone are detected. As shown in Fig. 9, the r k determined by this method (0.3) is within error of the value of 0.5 determined by analysis of the two substrates alone. Thus, the presence of additional competing substrates in the reaction does not have an appreciable effect on the magnitude of the relative V/K for ptRNA fMet47 and ptRNA Met82 (ϩ2).
Next, we determined the r k values for the remaining seven substrates using the ptRNA Met82 (ϩ2) as the reference substrate. As shown in Fig. 7, the r k value for the ptRNA Leu76 substrate is readily determined by this approach. This substrate has an r k of 3.5, indicating faster processing of ptRNA Leu76 over the reference ptRNA Met82 when they compete for RNase P processing. For this particular substrate, the ptRNA can be resolved from the unreacted reference ptRNA Met82 (ϩ2). This allows the change in the relative concentrations of the residual substrates to be quantified as well. As shown in Fig. 7B, we took advantage FIGURE 6. Competitive multiple turnover reactions containing both ptRNA fMet47 and ptRNA Met82 (؉2). A, PAGE analysis of the products of a reaction containing 5Ј 32 P end-labeled ptRNA fMet47 and ptRNA Met82 (ϩ2). The two precursors run as a single band indicated by a bracket denoting the presence of both substrates. The two leader sequence products are indicated by lines with or without the additional guanosines that identify the product from ptRNA Met82 (ϩ2). B, plot of the observed multiple turnover rate constants (v obs ) for ptRNA fMet47 (open symbols) and ptRNA Met82 (filled symbols) as a function of the relative concentrations of the two substrates. The data are fit to the log form of Equation 2 (Equation 5). C, plot of the observed multiple turnover rate constants (v obs ) for ptRNA fMet47 (open symbols) and ptRNA Met82 (filled symbols) as a function of the concentration of the third substrate ptRNA Leu76 . The data are fit to a mechanism in which ptRNA Leu76 acts as a competitive inhibitor (Equation 6). The inset shows the individual r k values determined from dividing the observed rate for ptRNA fMet47 by the observed rate for ptRNA Met82 at each of the different ptRNA Leu76 concentration. The solid line and dashed lines represent the average and standard deviation, respectively, calculated from this data set.
of internal competition analyses typically used to measure the relative rate constants for isotope effect measurements (51). The slower reacting ptRNA will become progressively enriched in the residual substrate population, and the relative rate constant can be determined by analyzing the change in substrate ratio as a function of the fraction of reaction. Using the ratio of residual precursor concentrations derived from the ratio of radiolabeled precursor bands, the r k for ptRNA Leu76 was determined by fitting to where R 0 is the ratio at the start of the reaction and R s is the ratio at the fraction of reaction (f) of the reference substrate (Fig. 7B)  (51, 59). The fraction of reaction for ptRNA Met82 (ϩ2) is determined from the intensity of its precursor and product bands. As expected, the faster rate constant for the ptRNA Leu76 substrate results in faster depletion of this substrate from the residual precursor population relative to the slower reacting ptRNA Met82 (ϩ2). As a result, the ptRNA Leu76 /ptRNA Met82 (ϩ2) ratio becomes progressively smaller as the reaction progresses. An essentially identical r k value of 3.4 is obtained from the fitting of the data shown in Fig. 7B.
Interestingly, in the course of experiments to determine the r k for ptRNA Ser80 , we detected two cleavage products in addition to correct RNase P cleavage at the mature tRNA 5Ј end. As shown in Fig. 8A, the reaction of ptRNA Met82 (ϩ2) yields a single product as expected, whereas the ptRNA Ser80 substrate   Fig. 2. The large variable arm of ptRNA Ser80 results in a substrate that is 15 nucleotides longer than ptRNA Met82 (ϩ2), which can be resolved under these gel conditions (15% PAGE). The ptRNA Ser80 substrate is cleaved by RNase P to give three products: the correct cleavage product resulting from cleavage 5Ј to N(ϩ1) (SER P1), miscleavage one nucleotide 5Ј into the leader sequence (SER P2), and miscleavage four nucleotides 5Ј to the correct cleavage site (SER P3). All of these products are resolved from the single cleavage product resulting from processing at the authentic 5Ј end of ptRNA Met82 (MET82). B, plot of product accumulation versus time for the products indicated in panel A showing the initial rates of product formation for the Ser P1 (open circles), Ser P2 (filled circles), and Ser P3 (open triangles) products relative to the accumulation of the product from ptRNA Met82 processing (squares). C, secondary structure of ptRNA Ser80 with arrows indicating the location of the three cleavage sites in ptRNA Ser80 by RNase P.
gives rise to three products (Fig. 8, labeled P1, P2, and P3). The P1 product maps to the expected site for RNase P processing between N(Ϫ1) and N(1). The P2 product is derived from miscleavage one nucleotide 5Ј to the authentic site, yielding a product one nucleotide smaller. Cleavage to give the P3 product occurs five nucleotides upstream of the correct site. RNase P cleavage in the leader sequence is not expected, although several studies have demonstrated the ability of the RNase P holoenzyme to cleave unstructured RNA, but with sequence or structure specificity that is not yet well defined (29,60). Alternatively, cleavage may result from alternative RNA folding (61,62). The r k for the miscleavage at P2 occurs at essentially the same rate as P1 (both have an r k of ϳ0. 6). Surprisingly, the r k value for miscleavage of the ptRNA Ser80 substrate at P3 occurs with an r k that is significantly higher (2.2). Although the precursors of both substrates can be resolved, the fact that ptRNA Ser80 reacts to form multiple products precludes determination of its relative rate constant by analysis of precursor ratios by Equation 6.
Nonetheless, as demonstrated in Fig. 8B, the relative rates of accumulation of the three products of ptRNA Ser80 are readily distinguished. For the ptRNA Gly62 , ptRNA Ile1 , ptRNA His37 , and ptRNA Gln85 substrates, the r k values were determined relative to ptRNA Met82 (ϩ2) from analyzing the initial rates of formation of the products. The r k values for all eight substrates, shown as the natural log to provide a linear scale, are compared in Fig. 9 together with the values for the alternative products for ptRNA Ser80 .

DISCUSSION
Understanding the competition between multiple, different cognate ptRNAs for processing by RNase P is important for understanding its function in vivo. This point is similarly true for the broad class of RNA-processing enzymes that recognize multiple cognate substrates. In this study, steady state kinetic analyses of two representative canonical and noncanonical ptRNA substrates were used to provide a framework for establishing a basic alternative substrate kinetic model for E. coli RNase P. For both substrates, the rate-limiting step for V is the substrate cleavage step; however, the rate-limiting step for V/K is substrate association. Because substrate binding is essentially irreversible relative to cleavage, both substrates compete for RNase P processing based on their relative rates of association. Using internal competition reactions, in which two or more substrates are present in the same reaction, we validated the basic features of the alternative substrate kinetic model. The results demonstrated that the relative rates of processing of two substrates directly reflect their relative V/K values (i.e. the specificity constant) and their concentrations. Also, the presence of additional substrates reduced the observed rate of processing of individual ptRNAs, but did not affect the observed r k for the two substrates being compared. Thus, for the set of substrates examined here, the E. coli RNase P enzyme follows simple alternative substrate kinetics; assuming that the results for ptRNA fMet47 and ptRNA Met82 are representative, the competition is governed by similar association rate constants. The results represent an important advance because they shed light on the biological function of RNase P in tRNA biosynthesis and provide a framework for quantifying relative processing rates in complex populations of competing tRNA precursors.
Biosynthesis of the translational machinery, including tRNA, consumes most of the resources in rapidly dividing cells (63). Importantly, the distribution of tRNA species is not uniform. The tRNAs that are present at higher concentrations are those that recognize the preferred codons of the genes encoding the highly expressed proteins of rapidly growing bacteria (64,65). This correspondence of codon usage and tRNA abundance is believed to maximize translation efficiency and therefore growth rates (66,67). The turnover of mature tRNA does not appear to be a major mechanism for the modulation of tRNA abundance (68). Thus, the steady state levels of tRNAs are largely set at the transcriptional level. In E. coli, tRNA biosynthesis and rRNA biosynthesis are tightly regulated by the stringent response, which senses the accumulation of uncharged tRNA and negatively regulates the initiation of transcription of tRNA and rRNA operons (69). The ptRNA substrates for RNase P can be transcribed individually or as part of polycistronic RNAs containing additional ptRNAs, rRNA, or mRNA (18,19). Precursors to individual tRNAs that are part of an rRNA primary transcript are thought to be released during the course of rRNA maturation by endonucleolytic cleavage by RNase III or RNase E (70). Separation of individual tRNA precursors from transcripts containing mRNAs or other tRNAs is accomplished primarily by RNase E (71,72).
The narrow range of V/K values for the model ptRNA substrates described here suggests that the relative rates of processing different ptRNAs in the cell will be proportional to the abundance of each precursor, assuming minimal influence of additional flanking RNA sequences. From this perspective, the alternative substrate kinetics of RNase P are tuned to be directly and uniformly proportional to rates of precursor bio- synthesis. It follows that the enzyme has a negligible role in influencing the steady state abundance, but rather functions to maintain the distribution set by precursor biosynthesis despite significant differences in the structure and context of each individual substrate.
This model is specific to ptRNA substrates, and indeed, there may be circumstances under which the rates of RNase P processing are modulated to contribute to regulation of gene expression. There are non-ptRNA substrates of RNase P including riboswitches (73)(74)(75) that are likely to have different kinetic properties related to their unique function and so are not accounted for in the scenario described above for ptRNA substrates. Additionally, the context of different RNase P substrates within larger polycistronic transcripts may influence relative processing rates (19) and clearly deserve further attention. For example, a suppressor tRNA substrate the length of the leader sequence has an effect on cell growth potentially due to effects on in vivo processing efficiency (76). Recently, Mohanty and Kushner (77,78) have identified two tRNA polycistrons in which RNase P is the primary processing event separating the individual tRNA units. RNase P cleaves 4 -7 nucleotides downstream of the CCA determinant, generating a substrate for additional processing nucleases in the substrate for tRNA Leu5 (79). The structures and binding modes that underlie RNase P cleavage at alternative sites in non-ptRNA substrates are poorly understood. Both bacterial (29,60) and eukaryotic (80,81) RNase P have been observed to catalyze cleavage at multiple sites distinct from authentic ptRNA-processing sites, and these forms of alternative specificity clearly deserve further attention. For example, we note that the cleavages of the ptRNA Ser80 substrate demonstrated here occur at a V/K that is within the range of processing at the correct site.
The results not only provide insight into issues specific to RNase P function, but also draw attention to parallels with the recognition capabilities of the ribosome. Studies by Uhlenbeck and colleagues (82,83) demonstrated uniform affinities, dominated by similar association constants, of elongation factor Tu for aminoacyl-tRNAs despite differences in tRNA and amino acid structure and chemical properties. Detailed analysis of misacylated tRNAs revealed thermodynamic compensation between the contributions of the amino acids and the tRNA moieties of the aminoacyl-tRNA to binding. The ribosomal A-site shows specificity for both the amino acid and the tRNA portions of their aminoacyl-tRNA substrates (84,85). Structure-function analysis of chimeric tRNAs presented that each tRNA sequence has coevolved with its anticodon to tune ribosome affinity to a value that is the same for all tRNAs (86). Thus, the observation of similarly uniform rates of ptRNA processing by RNase P provides an additional example of the tuning of RNA recognition to accommodate the structural variation in tRNA necessary for its function in aminoacylation and translation.
Fierke and colleagues (9, 87) demonstrated a conformational change during substrate binding by B. subtilis RNase P. Given the potential for the induced fit model to contribute to specificity, it is important to consider how the occurrence of a conformational change could impact alternative substrate competition. We assume that for cognate ptRNAs, any conformational change upon binding is favorable, whereas for noncognate RNAs, this step is unfavorable. As described previously, induced fit decreases the V/K for both a cognate and a noncognate substrate when the chemical step is rate-limiting and therefore does not intrinsically provide specificity (48). On the other hand, a conformational change can provide specificity when a binding step or product release step is rate-limiting for the cognate substrate, whereas the chemical step is rate-limiting for a noncognate RNA (48). As suggested by our results, association is likely to be broadly rate-limiting for ptRNA processing by E. coli RNase P. Thus, the presence of a conformational change could clearly enhance specificity over noncognate substrates.
In this study, we have taken a new approach to analyzing substrate recognition by RNase P by applying the perspective of alternative substrate kinetics. In addition to providing insight into the enzymatic behavior that underlies its biological function, the framework described here is useful for extracting relative rates by internal competition. With the simple competitive substrate kinetics of RNase P established, more broad application of competitive kinetics may be considered. In principle, internal competition methods are applicable to very large populations of substrates so long as reaction progress and the ratios of their precursors or products can be quantified.


extfieldsofstudy: Medicine Biology

text:
Linking Great Apes Genome Evolution across Time Scales Using Polymorphism-Aware Phylogenetic Models

The genomes of related species contain valuable information on the history of the considered taxa. Great apes in particular exhibit variation of evolutionary patterns along their genomes. However, the great ape data also bring new challenges, such as the presence of incomplete lineage sorting and ancestral shared polymorphisms. Previous methods for genome-scale analysis are restricted to very few individuals or cannot disentangle the contribution of mutation rates and fixation biases. This represents a limitation both for the understanding of these forces as well as for the detection of regions affected by selection. Here, we present a new model designed to estimate mutation rates and fixation biases from genetic variation within and between species. We relax the assumption of instantaneous substitutions, modeling substitutions as mutational events followed by a gradual fixation. Hence, we straightforwardly account for shared ancestral polymorphisms and incomplete lineage sorting. We analyze genome-wide synonymous site alignments of human, chimpanzee, and two orangutan species. From each taxon, we include data from several individuals. We estimate mutation rates and GC-biased gene conversion intensity. We find that both mutation rates and biased gene conversion vary with GC content. We also find lineage-specific differences, with weaker fixation biases in orangutan species, suggesting a reduced historical effective population size. Finally, our results are consistent with directional selection acting on coding sequences in relation to exonic splicing enhancers.

Introduction
The increased availability of sequenced genomes both from closely related species and from individuals of the same species, offers a great opportunity to study the speciation and evolutionary history of populations at different timescales, provided we can properly model the process of sequence evolution using inter-and intraspecific data together. The role of mutation and selection are of particular interest in this context. Mutation introduces genetic diversity, the raw material of evolution. Natural selection, along with neutral fixation biases and random genetic drift, can cause alleles newly introduced by mutations to increase in frequency and reach fixation. For comparative analysis that aim to detect selection and identify functional elements, disentangling the contribution of these forces is important.
In the past, phylogenetic methods focused on interspecies data, whereas population genetics was mainly concerned with intraspecies patterns. Classical population genetics methods can test the presence of selection, but do not include divergence data from multiple species (except as outgroups, see e.g., McDonald and Kreitman 1991;Schneider et al. 2011). Standard phylogenetic models instead infer substitution rates but not mutation rates and fixation biases.
There are a few exceptions, for example, the mutation-selection codon model of Yang and Nielsen (2008). This model assumes the same nucleotide mutational process for all codon positions, and estimates a fitness parameter for each codon, allowing to test the presence of selection on codon usage from interspecies data.
Some methods use both population genetics and phylogenetics models. For example, it is possible to estimate phylogenetic trees by reconstructing the genealogies of individuals from different species using the multi-species coalescent. Liu (2008) assumes no recombination within genes and free recombination among genes. RoyChoudhury et al. (2008) assume no new mutations, so that all the divergence among taxa originates from change in allele frequencies in standing variation. This method has been recently generalized to allow new mutations along the population tree, but is still limited to bi-allelic sites (Bryant et al. 2012). All these coalescentbased methods assume neutrality. Wilson et al. (2011) proposed a combined phylogeneticpopulation genetics approach that analyzes population data from different species and estimates a distribution of selective coefficients in coding regions. Recently, Gronau et al. (2013) developed a model similar to that of Wilson et al. (2011) and applied it to noncoding sequence data. Both the latter methods assume a standard substitution model along the phylogeny relating the species, and require all polymorphisms to be recent. In fact, the population genetics model allowing for intraspecific differences is only used at the tips of the tree, and not in inner phylogenetic branches and nodes.
Here, we introduce a POlymorphism-aware phylogenetic MOdel (PoMo), that, similarly to the model of Wilson et al. (2011), uses both polymorphism and divergence data simultaneously. However, for PoMo, we do not assume that polymorphisms originate from recent mutations. In our phylogenetic continuous-time Markov chain, polymorphisms are present both at terminal and ancestral nodes of the species tree. In this way, we can naturally account for ancestral shared polymorphisms (Clark 1997) and incomplete lineage sorting (Maddison and Knowles 2006;Pollard et al. 2006). Furthermore, by not assuming stationarity, reversibility, context-independence or mutational strand-symmetry, our model can describe complex mutational scenarios (Hwang and Green 2004;Polak and Arndt 2008). We show using simulations that, with our new model, we can accurately infer relative mutation rates and fixation biases, and that the inferences are robust to changes in demography.
One of the most intriguing aspects of the human genome is its exceptional heterogeneity. Base substitution rates differ among nucleotides, nucleotide contexts, genomic regions, and chromosomes (for a review see Hodgkinson and Eyre-Walker 2011). Knowing the intensity and variability of both mutation and fixation biases, whether due to selection or other forces such as biased gene conversion, is fundamental for interpreting evolutionary patterns (Ratnakumar et al. 2010). For example, coding sequence is a major determinant of fitness and adaptation (Eyre-Walker and , but undergoes peculiar evolutionary forces, with transcribed sequences evolving differently from the rest of the genome, showing, for example, strand-specific substitution rates (Hwang and Green 2004). It is therefore appealing to use synonymous sites as a neutral reference for coding sequence evolution, although selection can affect evolution of synonymous sites involved in the splicing process (Chamary et al. 2006;Parmley and Hurst 2007). Furthermore, mutation and fixation biases can have severe consequences on the fitness of individuals and populations Hodgkinson and Eyre-Walker 2011).
We performed a comprehensive study of evolutionary patterns of synonymous sites in great apes (humans [Homo sapiens], chimpanzees [Pan troglodytes], and orangutans [Pongo abelii and Pon. pygmaeus]). By using PoMo on polymorphism and divergence data simultaneously, we were able to overcome the limitations of previous studies, in particular disentangling the contributions of mutation and fixation biases to the evolution of synonymous sites. We first estimate global patterns of coding sequence evolution in great apes genomewide, including a comparison of lineage-specific trends. Then, we show evidence in favor of variation in mutation and fixation rates between genomic regions with different base composition, contributing to the long-standing debate regarding the origin and maintenance of GC-content variation (Eyre-Walker and Hurst 2001). Finally, we consider variation in evolutionary patterns within exons, examining evidence suggesting recent directional selection on synonymous sites.

New Approaches
We developed a new approach PoMo that uses polymorphism and divergence data simultaneously, and can estimate relative mutation rates, disentangling them from fixation biases in their contribution to substitution patterns. Similar to classical phylogenetic approaches (reviewed in Whelan et al. 2001), our model is a continuous-time Markov chain. We assume that a phylogenetic tree relates the species considered, and that nucleotide sequences evolve along it. Phylogenetic methods usually include only a reference genome for each species considered, and ignore intraspecies diversity. Here, in contrast, we use data from multiple withinpopulation individuals to infer allele frequencies for each taxon ( fig. 1A and see Materials and Methods). Similar to classical phylogenetic approaches, we model all sites as independent and discard haplotype structure. This way, we ignore information regarding recombination events, but we also bypass the problem of accounting for all possible coalescent trees, which is of elevated computational complexity when considering large samples (Dutheil et al. 2009).
We include polymorphisms as states of the Markov chain, in addition to the four nucleotide states of classical nucleotide models. In fact, in our Markov chain, two nucleotides can be present simultaneously at one site for one species/population. If a polymorphism is present at a tip, it means that the corresponding species has an observed polymorphism at the corresponding site and at the corresponding allele frequency ( fig. 1B and see Materials and Methods).
Although classical models assume instantaneous substitutions, we separate the mutation and fixation processes. In fact, we model sequence evolution as a gradual process made by small allele frequency changes ( fig. 1C and see Materials and Methods). As in classical phylogenetic models, the states in inner nodes/branches are usually unknown. This uncertainty is accounted for by considering the probability of each possible combinations of ancestral states via the Felsenstein pruning algorithm (Felsenstein 1981). In PoMo, we add the possibility of polymorphisms at various allele frequencies at inner nodes and branches ( fig. 1D and see Materials and Methods). This means that we account for ancestral polymorphisms and in particular for ancestral shared polymorphisms ) and incomplete lineage sorting (when two speciation events are separated by a lapse of time not sufficient for polymorphisms to reach fixation, see Maddison and Knowles 2006). The parameters in PoMo do not merely describe substitution rate, but are also informative of mutation rates, fixation biases, root nucleotide frequencies and branch lengths. All these parameters are estimated by maximum likelihood (ML) ( fig. 1E and see Materials and Methods).
Although many genomes (including the human genome) are not in base composition equilibrium, most phylogenetic models assume equilibrium and reversibility for convenience (Squartini and Arndt 2008). Here, we do not assume equilibrium or reversibility. Furthermore, because mutations in human coding sequences are thought to be strand-asymmetric and context-dependent (Hwang and Green 2004;

Simulations
To assess the precision of our methods in parameter estimation, we performed forward population genetics simulations with simuPOP (Peng and Kimmel 2005) on a phylogenetic tree. Our simulations closely mimic the features (divergence and diversity) of the great ape data set (see Materials and Methods). We reliably inferred the simulated parameter values when more than 10 5 sites were provided ( fig. 2 and supplementary figs. S10-S12, Supplementary Material online), far fewer than in the real data set (% 2Â10 6 , see Materials and Methods). We observed errors at levels of % 5% or below for branch lengths and ancestral and equilibrium nucleotide frequencies, and at most 10% for relative mutation rates. The intensity of selection was slightly  (Felsenstein 1981), we sum probabilities over all combinations of states at inner nodes. (E) The likelihood of all sites is combined, and the process is iterated with different parameter values until we find those that maximize the likelihood. These values (mutation rates, fixation biases, and root nucleotide frequencies) are our final estimates. underestimated, and required more data for acceptable inference ( fig. 2 and supplementary fig. S12, Supplementary Material online). This bias is probably due to the small number of polymorphic states used and to the fact that we ignore sampling variance.
We measured the running time of our method on simulated data (supplementary table S8, Supplementary Material online). ML estimations required less than 30 min on a single processor (2.66 GHz 6-Core Intel Xeon) for both the basic PoMo10 and the asy-CpG-PoMo10b models for a genome scale data set (up to 500 kB). This suggests that PoMo could be applied to genome-wide data from dozens of species simultaneously, that the state space could be expanded to include a larger virtual population (table 1), or to incorporate more model parameters to describe complex evolutionary scenarios. Performance for simulated data. Mutational and frequency parameters simulated were as estimated in the highest GC-content bin (see Materials and Methods). Intensity of selection for GC versus AT was set to 4N e s ¼ 1. On X axis is the number of sites in the data set used. Each box plot represents 10 simulations. The errors in the estimation, on the Y axis, were calculated as the Euclidean distance between the vector of estimated parameters and the true values, normalized by the Euclidean norm of the true vector. (A) Error in estimation of: fixation biases (6 entries vector, one for each substitution type, blue box plot), non-CpG mutation rates (6 entries, red), and CpG hypermutability (single-entry vector, orange). (B) Error in estimation of: branch lengths (green), ancestral nucleotide frequencies (pink), and equilibrium nucleotide frequencies (yellow). (C) Estimates of: GC versus AT fixation bias (blue), GC* in sites not preceded by C and not followed by G (yellow), GC* in sites not preceded by C and followed by G (green), and GC* in sites preceded by C and followed by G (red). The horizontal dashed lines represent the respective true values used for the simulations.
We also performed simulations including up to four different demographic events (bottleneck, expansion, migration, and reduction) in the same phylogeny (see Materials and Methods). Our approach proved robust in these cases (supplementary figs. S13-S15, Supplementary Material online). Demographic events usually bias estimation of selection in population genetics methods that use only polymorphism data (Haddrill et al. 2005;Keightley and Eyre-Walker 2007;Zeng and Charlesworth 2009). In contrast, we not only consider the site frequency spectrum, but also divergence patterns.

Analysis of Great Apes Whole-Exome Data
We extracted synonymous sites from coding sequence alignments of different species (human, chimp, and orangutan) and different individuals within species (see Materials and Methods). From these sites, using PoMo, we inferred relative mutation rates, fixation biases, and nucleotide frequencies at equilibrium and at the root. We first estimated global patterns from whole-exome data, then focused on variation between lineages, regions with different GC content, and sites in different positions within exons.

Global Estimates of Mutation Rates
Using different variants of PoMo, we estimated relative mutation rates from the global data set. Unsurprisingly, transitions (mutations between A and G and between C and T) had higher rates than transversions (the remaining mutations), in particular when CpG context was not accounted for ( fig. 3A). We compared our results with those of Lynch (2010), who measured the rate of new deleterious nonsense and missense mutations in humans. We find notable differences; in particular, Lynch (2010)  Regions near the transcription start site undergo peculiar substitutional patterns, with, in particular, reduced CpG context effects (Polak and Arndt 2008). Consistent with these observations, our hypermutability estimates greatly differ between first exons and other exons ( fig. 3B). After removing first exons, our relative mutation rate estimates are very similar to the phylogenetic estimates of (Duret and Arndt 2008; fig. 3B), despite the fact that they did not account for fixation biases and did not restrict their analysis to coding sequences.
Previous studies have suggested the presence of strandasymmetric substitution rates in human transcribed sequences (Hwang and Green 2004;Polak and Arndt 2008), and different asymmetries in different regions of the transcript (Polak and Arndt 2008). For this reason, we included strand-specific mutation and fixation biases in PoMo (see Materials and Methods, model asy-CpG-PoMo10b). Our analyses with this model support the idea that substitutional asymmetries are due to mutation rates, and not fixation biases ( fig. 3C). Furthermore, we detected the same mutational asymmetries as predicted by Polak and Arndt (2008), and, again as expected based on the latter study, first exons show different asymmetries from the other exons (supplementary fig. S3, Supplementary Material online).

Global Estimates of GC-Biased Gene Conversion and Base Composition
Although fixation biases can be caused by directional selection, the genome-wide fixation bias favoring GC over AT alleles in mammals is generally attributed to GC-biased gene conversion (gBGC; Duret and Galtier 2009). The effect of gBGC is similar to selection (Nagylaki 1983), and therefore the intensity of gBGC is usually expressed in terms of 4N e s. After accounting for context dependencies and mutational asymmetries (see Materials and Methods, model asy-CpG-PoMo10b), our estimate of gBGC is 4N e s ¼ 0:62, and is different from 0 according to both the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). On the other hand, accounting for differences in fitness between G and C, and between A and T does not improve the fit of the model according to AIC or BIC. We note, however, that the inferred fixation biases depend on the mutation model (supplementary fig. S4, Supplementary Material online; Hernandez et al. 2007).
Estimates of root, present and equilibrium nucleotide frequencies show that base composition is not at equilibrium in great apes, with GC content decreasing over time (

Lineage-Specific Estimates
Global patterns presented earlier should be interpreted as averages along the phylogenetic tree considered. However, substitution patterns vary considerably. For example, Polak et al. (2010) showed that the equilibrium GC content (GC*) in orangutan genes is lower than that in human and chimpanzee. GC* is determined by the GC/AT bias in substitution rates, and differences in substitution biases can be caused by changes in mutation rates or fixation biases. We investigated whether orangutan differs in mutation rates and/or fixation biases from human and chimpanzee. We allowed parameters to differ between the human-chimp lineage (comprising the human and chimpanzee branches, and the branch of their ancestor) and the orangutan lineage (comprising the two orangutan branches, and the branch of their ancestor). Allowing lineage-specific fixation biases and mutation rates resulted in improvements in both AIC (25.31) and BIC (12.83) scores (see table 2 for details). Our estimate of gBGC intensity in human-chimp was almost double that in orangutan (% 0:7 vs. % 0:35; table 2), even when we allowed different

2253
Great Apes Genome Evolution using PoMo . doi:10.1093/molbev/mst131 MBE mutation rates in different lineages. We conclude that a lower gBGC, possibly due to reduced effective population size, contributed to the lower GC* estimated in orangutan.

Variation among Exons
One well-studied aspect of genomic variation in mammals is GC content (Bernardi 2000 (Bernardi 2000), have been suggested to cause variation in GC content in mammals, but until now, no analytical framework was available to infer the relative importance of these processes. Our new approach represents a great opportunity for disentangling, and quantifying, mutational, and fixation biases variation. We binned exons according to their GC content at synonymous sites (GC4) from lowest to highest, so that all bins have roughly the same number of sites (% 3:25Â10 5 ). GC4 strongly correlates with regional GC content (Clay et al. 1996

Mutation type
Relative rate Mutation rates from AC Mutation rates from GT C B A Fig. 3. Estimates of mutation rates in great apes. (A) Estimates of relative mutation rates by Lynch (2010) in humans (blue) and PoMo10 on great ape data (red). mAC stands for mutation rate from A to C, etc. Values on the Y axis represent mutation rates normalized by ð AC + AG + AT + CA + CG Þ. (B) Estimates of relative mutation rates by Duret and Arndt (2008) in human-chimp (blue), CpG-PoMo10 in great apes without first exons (red), and CpG-PoMo10 on first exons only (green). hmCT represents the hypermutability from C to T and from G to A in CpG context. Error bars in (A and B) show the profile likelihood 95% confidence intervals. (C) Mutation rates from A and C nucleotides (red) compared with mutation rates from G and T (blue). In both cases, we refer to the nucleotide on the sense strand. We paired reverse-complement mutation types to remark strand-asymmetries. All rates are estimated with asy-CpG-PoMo10b on the whole data (see Materials and Methods). estimated mutation rates and fixation biases for each bin separately. Although most mutation rates vary only slightly, CpG hypermutability shows very large differences, being strongest in GC-poor exons and weakest in GC-rich exons ( fig. 4A). Among the other mutation rates, most noticeably, AC increases with GC content. These results are robust to the number of bins used, and to the exclusion of short exons (supplementary figs. S6 and S7, Supplementary Material online).
gBGC also increases with GC content, ranging from % 0:2 to % 1:2 ( fig. 4B). Even after removing the potential biases coming from the first and second exons of each gene ( fig. 3B), mutation rates and gBGC still vary between the extreme GC bins according to both AIC and BIC scores (table 3; supplementary table S2, Supplementary Material online). Although we accounted for many mutational biases, modeling site variation in total mutation rates still resulted in a model improvement according to AIC and BIC (supplementary table S7, Supplementary Material online). This suggests that more context-dependent or cryptic factors in mutation rate variation exist (Hodgkinson et al. 2009).
It has been inferred that GC content of GC-rich regions is decreasing in mammals (Duret et al. 2002;Belle et al. 2004;Gu and Li 2006). Alvarez-Valin et al. (2004) claimed that this result can be explained with a bias in the method of inference (parsimony), and the presence of context-dependent mutations, regional variation, indels, and alignment errors. Here, we account for these problems by using an ML context-dependent model, and by analyzing synonymous sites of closely related species (which are expected to contain negligibly few alignment errors and indels). Although we observe that GC* is highest in GC-rich exons, the difference in GC* among bins is considerably smaller than the difference in present or root GC content, meaning that base composition is becoming homogenous across the genome (fig. 5). Furthermore, except for the GC-poorest bin, GC* is always lower than present and root GC content. Number of bins used and exclusion of short exons did not affect these results (supplementary fig. S8, Supplementary Material online).

Variation within Exons
Finally, we addressed the issue of variation in evolutionary patterns between exonic positions. Different regions within an exon can vary in substitutional and compositional trends. 5 0 -and 3 0 -ends of nonterminal exons are GC-poorer than exon centers, and have stronger codon usage bias (Willie and Majewski 2004). This has been interpreted as the effect of selection on splicing motifs (reviewed in Chamary et al. 2006). In agreement with this hypothesis, codon usage bias in exon boundaries fits splicing motifs more than in exon centers (although not for every amino acid, Parmley and Hurst 2007). We analyzed variation in mutation and fixation biases within exons. We excluded terminal exons and divided sites in 3 bins. The 5 0 -bin contained the first 5 synonymous sites in each exon; the 3 0 -bin contained the last 5 sites; the central bin contained all the remaining sites. On each bin, we then estimated mutation rates and fixation biases as before.
We observe a higher GC content in exon centers, but extremely similar equilibrium frequencies in all bins ( 6B). Models that allow for these differences are preferable according to AIC, but not BIC (table 4; supplementary table  S3, Supplementary Material online). Nevertheless, estimated differences are larger than expected just by error according to simulations (cf. supplementary figs. S9C and D, Supplementary Material online).

Discussion
Understanding intensity and variation of mutation and fixation biases is fundamental for the interpretation of evolutionary patterns. With our new model, PoMo, and with genome-scale data of within and between-species diversity, we disentangled and estimated mutational and fixation biases at synonymous sites of great apes.
Our estimates of mutation rates considerably differ from those of Lynch (2010), the only previous study, to our knowledge, that inferred relative mutation rates in humans while accounting for fixation biases ( fig. 3A). Part of the difference might be due to the timescale considered (recent mutations in Lynch 2010, polymorphism and divergence here), and to the fact that we also include data from other great apes. However, we think that most of the discrepancy derives from transversions being over-represented in the missense  fig. S1, Supplementary Material online). Our analysis shows that estimates of relative mutation rates using phylogenetic data, as in Duret and Arndt (2008), are preferable, but for further comparative studies we suggest to use estimates from models accounting for fixation biases such as PoMo.
The strongest fixation bias that we detected favors GC over AT. In mammals, this phenomenon is generally attributed to gBGC. We inferred slightly lower estimates of gBGC than previous studies. Spencer et al. (2006) estimated the intensity of gBGC in the human genome from the allele frequency spectrum, while Lynch (2010) contrasted mutational patterns with nucleotide frequencies. The first approach considers the recent evolutionary past, while the second is informative of the fixation bias in the long term, probably on the scale of hundreds of millions of years (Duret et al. 2002). Our approach is intermediate, as it considers polymorphisms and divergence, but not base composition. Lynch (2010) estimated a gBGC intensity of 4N e s % 0:99, which was within the range 0:5 < 4N e s < 1:3 from Spencer et al. (2006). We estimate 4N e s to be % 0:62. Although these values are not necessarily comparable, we recognized some additional causes for the small discrepancy. First, our method tends to slightly underestimate gBGC ( fig. 2; supplementary fig. S12, Supplementary Material online). Second, we studied humanchimpanzee-orangutan data, and not only human. We observed a lower intensity of gBGC on the orangutan lineage (% 0:35 vs. % 0:7 of the human-chimp lineage; table 2). A lower fixation bias in favor of GC nucleotides causes a shift in substitution rates towards AT, and therefore a reduction in GC*. A reduction in GC* in orangutan was previously observed by Polak et al. (2010). We suggest that the most likely explanation for reduced gBGC fixation bias 4N e s in orangutan is a difference in historical effective population size (N e ), and not in the molecular repair bias itself (s). This is consistent with studies suggesting that N e in orangutan was smaller than on human-chimp lineage: 65,000 AE 30,000 for the human-chimp ancestor and 45,000 AE 10,000 for the human-chimp-gorilla ancestor (Hobolth et al. 2007), in contrast to 26,800 AE 6,700 for the Bornean and Sumatran orangutan ancestor (Mailund et al. 2011 (B) Estimation of fixation biases with the strand-specific asy-CpG-PoMo10b. GC-sAT represents the apparent selective advantage of GC versus AT, sC-sA between C and A, sG-sA between G and A, and sT-sA between T and A. Furthermore, we investigated variation of mutation and fixation biases along the exome with respect to present GC content. Regional variation in GC content is one of the most fascinating aspects of mammalian genomes, but its causes and consequences are still not well understood. Although some authors suggested that selection is the cause of GC content variation origin and maintenance (Bernardi 2000), most studies proposed neutral explanations, such as variation in mutation rates (Fryxell and Zuckerkandl 2000;Fryxell and Moon 2005) or gBGC ). Here, we jointly estimated variation in mutation and fixation biases, therefore accounting for possible confounding effects of one on the other. We conclude that mutation rates (and in particular CpG hypermutability) and gBGC vary with base composition ( fig. 4). Nevertheless, GC content decreases, and over time becomes more homogeneous across the genome ( fig. 5), as concluded by previous studies as well (Duret et al. 2002(Duret et al. , 2006Meunier and Duret 2004). One of the possible explanations for the homogenization of GC content are changes in the recombination map, and therefore in gBGC intensity (Auton et al. 2012). Otherwise, a reduction in effective population size may have led to a decrease in the intensity and variation of gBGC effects. Additional studies on different mammalian clades are necessary to determine the most likely scenarios.
Finally, we measured differences in evolutionary patterns between exon centers and boundaries. Previous studies of base composition suggested that synonymous sites in different exonic positions are subject to different selective pressures due to splicing motifs (reviewed in Chamary et al. 2006). We confirmed these trends; in fact, we measured a fixation bias favoring A over T in boundaries of nonterminal exons, and T over A in exon centers ( fig. 6B). This observation cannot be explained with gBGC, and is consistent with expectations of the hypothesis of selection on exonic splicing enhancers, since those are A-rich and T-poor (Parmley and Hurst 2007). However, our findings are only marginally significant, and   Great Apes Genome Evolution using PoMo . doi:10.1093/molbev/mst131 MBE they might be improved by including more individuals from more species in the future.
In conclusion, we presented a new phylogenetic model, PoMo, that can infer mutation and fixation biases from patterns of polymorphisms and divergence in populations/ species related by any arbitrary history. We provide the software to replicate our analyses that assumes a sample of 10 sequences per species. This is a limitation of our implementation, and not of the model. In fact, the number of haplotypes considered from each population as well as the number of sites do not represent a considerable computational burden for our methods. Furthermore, as PoMo10 has fewer states than a codon model, it can be applied to phylogenies spanning several dozens of taxa (Seo and Kishino 2009;Gil et al. 2013).
PoMo can also be applied to the estimation of phylogenetic trees from population data. In fact, it has the potential to improve the resolution of short branches (relative to N e ), where classical methods for phylogenetic inference often fail due to incomplete lineage sorting and shared ancestral polymorphisms. These issues can already be accounted for by either making strong assumptions regarding recombination events (independent loci and no recombination within loci, see e.g., Heled and Drummond 2010) or by computationally demanding hidden Markov model approaches (Mailund et al. 2011). But, unlike PoMo, neither method is applicable to large numbers of individuals within-species, due to many possible coalescent trees.
Until now there are only few clades with genome-wide population data available from multiple species, such as primates and model organisms. This number is expected to grow in the near future, providing great opportunities to understand changes in evolutionary patterns.

Polymorphism-Aware Phylogenetic Models
Model Background Phylogenetic substitution models represent DNA evolution as a continuous-time Markov process along a phylogenetic tree t (for a review see Whelan et al. 2001). Different sites are generally assumed to evolve independently. Each point of t represents a taxon at an instant, and tree bifurcations correspond to speciation events. Another common assumption is homogeneity: the evolutionary process does not change through time and among species. Nucleotide substitution models associate the points of t to elements of the state space fA, C, G, Tg with certain probabilities. If a point of t is in state C, it means that the corresponding taxon, at the corresponding time, had nucleotide C at the considered site of the genome. The phylogeny tips correspond to present, observable states.
States assigned to taxa can change in time, and the continuous-time Markov process modeling these changes is defined by an instantaneous rate matrix Q. Each entry Q IJ of Q, with I 6 ¼ J, is the rate at which nucleotide I is replaced by nucleotide J. Given Q, for any branch b of length t in t it is possible to calculate the transition probability matrix PðtÞ ¼ e Qt . Entry P IJ ðtÞ of PðtÞ is the probability that the end of branch b is in state J, conditioned on the start being in state I. Matrix P(t) is used to calculate the likelihood LðÞ of any parameter values y via the Felsenstein pruning algorithm (Felsenstein 1981). The likelihood is the conditional probability of the data D given y. Data D consist of DNA sequence alignments. Parameter estimates can be obtained by ML, that is, by determining the parameter values y that maximize the likelihood function.

State Space
Our new models are similar in most aspects to standard phylogenetic nucleotide substitution models described earlier. The most important difference is that we expand the state space. We do not only include four states associated to nucleotides, but also further states representing polymorphisms. Nucleotide states in the new models represent sites with a fixed allele. Assuming at most two alleles per site per time per taxon, there exist six types of polymorphisms determined by the alleles simultaneously present in a taxon: fA,Cg, fA,Gg, fA,Tg, fC,Gg, fC,Tg, and fG,Tg. We define PoMo N (POlymorphism-aware MOdel with virtual population size N) as a phylogenetic model with N -1 polymorphic states for each type of polymorphism. PoMo N state space has therefore 4 + 6ðN À 1Þ elements, which means that any taxon at any considered time point can be assigned to any of the 4 + 6ðN À 1Þ states. The N -1 states associated to the same polymorphism type represent different allele frequencies within a virtual population of N haploid individuals. For example, polymorphic state i (1 i N À 1) of type fA,Cg represents a frequency of i/N for allele A and ðN À iÞ=N for allele C in a virtual population of N + 1 haploid individuals.

Definition of Instantaneous Rates
It is a common practice in population genetics to approximate the dynamics of a large real population with those of a small virtual population Kaiser and Charlesworth 2009;Zeng and Charlesworth 2009). To our knowledge, we present the first application of this approach to phylogenetics. We assume that the real population has effective population size N e , mutation rate per generation IJ from allele I to J, fitness parameter s I for allele I, and number of generations t. We define our virtual population as evolving according to a Moran model (Moran 1958) with population size N, mutation rate per generation IJ from allele I to J, fitness parameter s I for allele I, and number of generations t. The dynamics of a virtual population with properly scaled mutation rate (4N e % 4N) and selection coefficient (2N e s % 2Ns), are a good approximation of the real population, if time is scaled by the population size in both cases (t/N for the virtual population and t=N e for the real one) and if N is sufficiently large (Zeng and Charlesworth 2009, observed that even with N as small as 10 reasonable results could be achieved).
We make the assumption that the scaled mutation rate (4N e ) is low, and allow mutations only at monomorphic sites in our virtual population (Vogl and Clemente 2012). Therefore, while our model is four-allelic, only two alleles can be present simultaneously in one population at a site. We represent the polymorphic state with i virtual individuals carrying allele I, and N -i carrying allele J, as i I N À i J .
The probability that the virtual population in a polymor- Similarly, the probability to evolve from The probability with which a new allele is introduced in the virtual population, that is, of evolving from monomorphic state I to polymorphic state N À 1 I 1 J in one virtual generation is Within a single generation no other changes are allowed, in fact, virtual allele counts can only increase or decrease by one per generation. We call M N the matrix of probabilities of allele frequency changes for one generation. Matrix M N has dimension equal to the number of states, 4 + 6ðN À 1Þ. The last step in defining our model is transforming the Markov chain from discrete-time (in number of generations) into continuous-time. A continuous-time Markov chain is defined by its instantaneous rate matrix Q. We set the instantaneous rate matrix of our continuous-time process as Q N :¼ NðM N À IÞ, where I is the identity matrix. Then, the probabilities of state changes in coalescent time t/N will be given by P t where t as before represents the number of virtual generations, but can now take noninteger values. We list the entries of the rate matrix Q N in supplementary table S9, Supplementary Material online.
After fitting the matrix Q N to real data by ML, we estimated the scaled fitness parameters in the real population (2N e s I ) as 2Ns I . The four fitness parameters (s A , s C , s G , and s T ) are defined up to an additive constant, and therefore correspond to three free parameters. When only gBGC is expected to drive fixation biases, we set s A ¼ s T and s C ¼ s G , reducing the number of free parameters describing fitness differences to one. Likewise, we estimated the scaled mutation rates (4N e IJ ) as 4N IJ .

Root Frequencies
Stationarity and reversibility are common and mathematically convenient assumptions for phylogenetic models, but are often not realistic (Galtier and Gouy 1995;Yang and Roberts 1995;Akashi et al. 2006;Gu and Li 2006). Here, we do not assume them. Because of nonstationarity of our model, state frequencies might change along t, and root state frequencies p might differ from the observed frequencies. To define the 4 + 6ðN À 1Þ entries of p, we use three additional free parameters ( A , C , and G , where T :¼ 1 À A À C À G ) representing relative frequencies of fixed nucleotides at the root.
The root frequency of the polymorphic state with i virtual individuals carrying allele I, and N -i carrying allele J is: where K norm is a normalization factor, so that all root frequencies of polymorphic states sum up to pol . Under neutrality and rare mutations, the expected proportion of polymorphic sites with derived allele count i in a sample of N individuals is proportional to 1/i (e.g., see eq. 4.20 in Wakeley 2009). A root polymorphism can be derived from both alleles present in the population, we therefore take into account both possibilities in equation (4). The proportion of polymorphic states at the root, pol , is not a free parameter, but is set equal to the observed proportion of polymorphic states. In fact, pol could not be reliably estimated via ML, and its value did not affect the estimation of other parameters noticeably (supplementary information and table S1, Supplementary Material online). The root frequency of a fixed state I is I ð1 À pol Þ.

PoMo10
All results in this study are based on PoMo10 (PoMo N with N = 10). PoMo10 has 58 states: 4 fixed states and 54 polymorphic states. In fact, for each of the 6 pairs of alleles (fA,Cg, fA,Gg, fA,Tg, fC,Gg, fC,Tg, and fG,Tg), there are 9 polymorphic states for the possible allele counts (f9,1g, f8,2g, . . . ,f1,9g). Therefore, PoMo10 has lower computational cost than a standard codon model (61 states), allowing genome-wide analysis of phylogenies with considerable numbers of species. PoMo10 approximates real population dynamics with those of a virtual population of 10 individuals. Although this is a rough approximation, it is expected to be sufficient for parameter estimation (Zeng and Charlesworth 2009, also confirmed by our simulations results). Smaller values of N generally resulted in considerable For each species and site considered, we randomly extracted a sample of haploid size 10 (see Description of Data), and trivially associated the observed allele frequencies to the corresponding virtual frequency state, ignoring sampling variance. A limitation of PoMo10 is that it requires 10 sampled sequences for each species. A larger N is expected to result in improved estimates, but at computational costs (table 1).
PoMo10 Extensions: CpG Hypermutability and Strand-Asymmetry When we assume no context-dependency or strand-asymmetry, we define six free parameters to describe mutation biases, one for each unordered pair of nucleotides ( AC ; AG , etc.). Yet, mutation rates in mammals show strong dependency on the neighboring bases (Hwang and Green 2004). We extended PoMo10 to include the strongest context dependency, the hypermutability of CpG (nucleotide C followed by G) toward TpG or CpA. This is accounted for by an extra parameter h CT that describes the mutation rate from C to T and from G to A in a CpG context. We call this model CpG-PoMo10 (for detailed description of rates see supplementary information, Supplementary Material online). To estimate parameters of CpG-PoMo10, we only used synonymous sites whose preceding and following bases are constant among the considered species. We generally have three free parameters describing nucleotide frequencies at the root, but with context dependency this number increases to 12 to account for different frequencies in different CpG contexts.
We further extended CpG-PoMo10 to account for hypermutability of transversions in CpG context (from CpG to ApG, GpG, CpC, and CpT). The resulting model, CpG-PoMo10b, has two additional free mutational parameters (for details see supplementary information, Supplementary Material online).
Finally, we accounted for strand-specificity of mutation rates (Hwang and Green 2004;Polak and Arndt 2008). In the resulting model, asy-CpG-PoMo10b, the constraints for strand-symmetry (e.g., AG ¼ TC ) are relaxed, and therefore 9 extra free mutational parameters are necessary with respect to CpG-PoMo10b, for a total of 18 (for details see supplementary information, Supplementary Material online).

Model Implementation
In this study, we assume that the tree topology is known and fixed (supplementary fig. S17, Supplementary Material online), and we estimate all branch lengths in the 4-species rooted tree. Lists of number of free parameters for different models are included in tables 2-4.
Parameter estimation was performed via ML with the conjugate gradient algorithm implemented in HyPhy (Pond et al. 2005). For this purpose, we produced custom scripts in HyPhy Batch Language (supplementary file S1 [Supplementary Material online] describes asy-CpG-PoMo10b and PoMo10 for the case s G ¼ s C and s A ¼ s T ). We always estimated all free parameters simultaneously. The scripts that we provide require 10 sequences sampled from each species. This is a limitation of the present state of our software, but not of the model (table 1 and discussion). We used different starting points for the conjugate gradient iterations on the whole real data set, and observed consistency of different optimization runs (supplementary information and supplementary fig. S2, Supplementary Material online).

Description of Data
Great Apes Data Set We constructed an exome-wide, inter-and intraspecies data set of alignments of 4-fold degenerate (synonymous) sites from H. sapiens, P. troglodytes, Pon. abelii, and Pon. pygmaeus (respectively human, chimpanzee, and Sumatran and Bornean orangutan).
First, CCDS (Pruitt et al. 2009) alignments of H. sapiens, P. troglodytes, and Pon. abelii (references hg18, panTro2, and ponAbe2) were downloaded from the UCSC genome browser (http://genome.ucsc.edu, last accessed August 8, 2013). Only CCDS alignments satisfying the following requirements were retained for the subsequent analyses: divergence from human reference below 10%, no gene duplication in any species, start and stop codons conserved, no frame-shifting gaps, no gap longer than 30 bases, no nonsense codon, no gene shorter than 21 bases, no gene with different number of exons in different species, or genes in different chromosomes in different species (chromosomes 2a and 2b in nonhumans were identified with human chromosome 2). From the remaining CCDSs (9,695 genes and 79,677 exons), we extracted synonymous sites. We only considered third codon positions where the first two nucleotides of the same codon were conserved in the alignment, as well as the first position of the next codon.
Then, population data were added to the species alignments. Human single nucleotide polymorphisms (SNPs) from 59 Yoruban (Nigerian) individuals (haploid sample size 118) sequenced from the 1,000 genomes pilot project (1000 Genomes Project Consortium 2010) were downloaded (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/, last accessed August 8, 2013) and included the alignments. Similarly, we added SNP data of 10 western chimpanzee individuals (haploid sample size 20) sequenced by the PanMap project (Auton et al. 2012) and downloaded from ftp://birch.well. ox.ac.uk/haplotypes/ (last accessed August 8, 2013). Orangutan SNP data for the two species considered, each with five sequenced individuals (haploid sample size 10, Locke et al. 2011), were kindly provided to us by X. Ma (in preparation) and are now available online (http://www.ncbi. nlm.nih.gov/projects/SNP/snp_viewTable.cgi?type=contact& handle=WUGSC_SNP&batch_id=1054968, last accessed August 8, 2013). We sub-sampled 10 alleles without replacement for each species and site. The final total number of synonymous sites included was 1,950,006 (for more details on real data sets see supplementary tables S4-S6, Supplementary Material online).
The collection of all synonymous site alignments in the great apes data set in valid format for PoMo10 is provided as supplementary file S2, Supplementary Material online. Custom scripts to convert SNP data from VCF v4.0 format into PoMo10 states, and to convert multi-species alignments into HyPhy input files, are also provided as supplementary file S3, Supplementary Material online.

Simulations
We simulated a population of 50 diploid individuals evolving according to a phylogenetic tree (supplementary fig. S17, Supplementary Material online), where a branch bifurcation represents the duplication and split of a population. Evolution was simulated according to a Wright-Fisher model with sexual reproduction using simuPOP (Peng and Kimmel 2005) and custom Python scripts. Phylogeny and mutation rates were set so to have similar divergence and diversity levels to those in real data (supplementary table S4, Supplementary Material online).
First, we simulated five scenarios, in which we progressively added demographic events. The first scenario consisted of a constant-size population phylogeny. In the second scenario, we added a bottleneck on the human branch. In the third, we made a population expansion follow the bottleneck. In the fourth, we further added migration between the two orangutan species. Finally, we reduced population size in the second half of the chimpanzee branch (for further details about simulated demographic events see supplementary information, Supplementary Material online).
In a second set of simulations, we included CpG context. We used root frequencies and mutation rates as estimated from GC-rich and GC-poor bins (the extreme bins in figs. 4 and 5; for a detailed description of mutational parameters in simulations, see supplementary information, Supplementary Material online). For both the GC contents considered, we simulated three selective regimes with a GC versus AT fitness difference of 4N e s 2 f1,0, À 1g, for a total of six scenarios. For each scenario, we simulated 10 6 independent sites, and then sub-sampled data sets of varying sizes (ranging from 10 4 to 5 Â 10 5 sites), with 10 replicates for each size.