CpG Substitution Rate Exhibits Substantial Regional Variation

We analyzed approximately 38 million orthologous sites from human, chimpanzee, and baboon obtained by aligning genomic DNA segments from these species. We used only high quality Bacterial Artificial Chromosome (BAC) based sequences in all of our analyses. Sequences that resulted from low-coverage whole genome shotgun sequencing projects, which may harbor errors in sequencing and assembly, were not used. We then extracted only non-coding sequences (see Methods) to analyze patterns of substitution rate variation free from the effect of natural selection. Our final non-coding data set included approximately ∼14.7 million aligned sites from 17 chromosomes.

We used a parsimony method [19]–[23] to identify sites that have been part of a CpG dinucleotide in the recent past (“CpG sites”). In addition, using the same method, we identified sites that have been part of a GpC dinucleotide (“GpC sites”, see the Methods section for details). Because GpC dinucleotides consist of the same bases (C and G) as CpG dinucleotides, while not involved in DNA methylation [24], they are often used as a dinucleotide control for CpG sites [13],[17],[25]. We also analyzed all sites that have not been a part of a CpG dinucleotide during the given evolutionary timescale (“non-CpG sites”, [22]). Note that GpC sites are a subset of non-CpG sites. Results from non-CpG sites were similar to that from GpC sites (see below).

It is important to note that a great majority of CpG sites in certain regions, called CpG islands, of mammalian genomes are typically hypo-methylated [26] and hence do not undergo methylation-origin mutation process. Therefore, it is crucial to exclude CpG islands from our analyses. We used similar but slightly more stringent criteria to those proposed by Takai and Jones [27], a widely used method, to identify and exclude CpG islands in our data (see Methods). We chose to use a shorter length cutoff than the conventional length as proposed in an earlier study [28], so that our method will be conservative in terms of masking true CpG islands.

Table 1 describes the numbers of CG->TA transition substitutions in CpG and GpC sites in our dataset. Note that even though there are over an order more GpC sites than CpG sites in our data, the total numbers of CpG and GpC transition substitutions are similar. This observation confirms that CpG substitutions occur much more frequently than other types of substitutions in the human genome [2],[29],[30]. For simplicity, we refer to the proportion of the number of transition substitutions to the total number of (CpG or GpC) sites as the “rate of substitution” in the rest of our paper, which implies that the unit of timescale is since the divergence of the genomes of humans and chimpanzees. The rate of CpG substitution in intergenic regions and introns are 14.29±0.4%, and 13.88±0.5%, respectively. The rate of GpC substitution in intergenic regions and introns are 1.12±0.037% and 0.98±0.043%, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Description of the Dataset.

https://doi.org/10.1371/journal.pcbi.1000015.t001

When all single nucleotide substitutions are considered, their rates vary substantially among different regions, more than expected solely from stochastic effects [2]. To examine regional variation of CpG substitution rates, we plotted the rate of CpG substitution in 50 kb segments of non-coding regions with at least 10 kb aligned sites (Figure 1A). The mean CpG substitution rate in these segments is 15.3%. The observed standard deviation is 6.3%, which is significantly greater than the standard deviation expected under a model that assumes uniform CpG substitution rate in all the segments (95% confidence interval [CI] 4.2%–5.8%; see Methods). As expected, GpC sites and non-CpG sites also exhibited substantial variation (Figure 1B). These results remained the same when we changed the size of the windows examined.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Histogram of the rate of CG->TA substitutions in non-coding regions.

The rate of CG->TA transitions in CpG sites (A), and GpC sites (B) in 50 kb segments of non-coding regions having at least 10,000 aligned sites. Variation of CpG substitution rate among non-coding regions is significantly greater than that expected under a uniform substitution rate model. A similar result was obtained for GpC sites (see text).

https://doi.org/10.1371/journal.pcbi.1000015.g001

The rates of CpG substitution and those of non-CpG substitutions are significantly correlated in our sample (Pearson correlation coefficient ρ = 0.32; P<0.001). Because CpG substitution rate did not follow a normal distribution, we also analyzed log-transformed data and obtained a similar result (ρ_tr = 0.31; P<0.001).

CpG Substitution Rate in Non-Coding Regions of Primate Genomes Is Negatively Correlated with G+C Content

Earlier studies reported conflicting results on the relationship between the G+C content and the CpG substitution rate of a genomic region. Some studies have proposed a negative relationship [13],[17],[18] while others observed no correlation [19]. We analyzed this relationship by dividing the non-coding regions into 6 equal size bins based on their G+C content and plotted the rate of CG->TA substitution at CpG sites in each bin against its average G+C content (Figure 2A). To avoid the effect of variation in methylation efficiency, sites in transposable elements were excluded from analyses in this and the following section.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Negative correlation between the rate of CG->TA substitution and G+C content of non-coding regions.

(A) Non-coding regions (intergenic and introns) were partitioned into six equal-sized bins based on their G+C contents. The rates of CG->TA substitutions in CpG sites of these bins are negatively correlated with their G+C contents. The negative relationship holds when introns were analyzed separately. In case of intergenic regions the relationship was not significant. Nevertheless, we found a negative trend. (B) GpC substitution rates in non-coding regions exhibited a negative relationship with G+C contents. When divided into intergenic regions and introns, however, the relationships were not significant, although there was a clear negative trend. Refer text for r² values and P- values.

https://doi.org/10.1371/journal.pcbi.1000015.g002

We observed a significant negative correlation between the G+C content of the non-coding regions and the CpG substitution rate (Figure 2A; r² = 0.642, P = 0.032). When the CpG data was partitioned into introns and intergenic regions, the negative relationship with G+C content was significant in introns (r² = 0.71, P = 0.021), but not in intergenic regions although there was a clear negative trend (r² = 0.29, P = 0.15). The average length of the intron segments is 5.5 kb and of the intergenic segments is 119 kb.

We also investigated rates of CG->TA substitutions at GpC sites, as a control for sites that are not affected by DNA methylation [24]. We observed a negative correlation with the G+C content of the non-coding regions (Figure 2B; r² = 0.86, P = 0.004). When the GpC data was partitioned into introns and intergenic regions, the relationships were not significant [r² = 0.313 (P = 0.14), r² = 0.341 (P = 0.10) for introns and intergenic regions, respectively]. Nevertheless, we observed negative trends in both introns and intergenic regions (Figure 2B).

Distance-Decaying Relationship between CpG Substitution Rate and G+C Content

The negative relationship between CpG substitution rate and G+C content (see above) is consistent with the thermodynamic requirement during deamination process (see Introduction). However, given that transition rates at other sites (such as GpC sites) also exhibit negative relationship with G+C content (see above) and that thermodynamic hypothesis does not necessarily explain long-range effect of G+C content on CpG substitution rates, we proposed the “distance-decaying” hypothesis.

We now present our main results on the distance-decaying influence of G+C content on rates of CpG substitution. We performed a sliding window analysis using a window size of 200 bps and a step size of 25 bps (partially overlapping windows) to analyze the relationship between G+C content at varying distances from the CpG site and the rate of CpG substitution. Because the average length of the introns in our dataset was only ∼5.5 kb [small as compared to the average length of intergenic regions (119 kb)], a large proportion of CpG sites in introns will have a large portion of their sliding windows lie in exons (the sliding windows extended up to 5 kb around each CpG site, see below). Therefore, we used only intergenic CpG sites in this analysis, to minimize the effects of constraints related to coding for amino acids and natural selection. When we included a subset of intronic CpG sites that are at least 3 kb away from an exon in our analyses, the results were the same (Figure S5). Thus, the observed pattern is a common feature of non-coding regions in primate genomes.

At each distance, we binned CpG sites with similar G+C content (as measured from the 200 bps sliding window) at that distance from the site. The cutoffs used for binning were <38%, 38%–45%, 45%–52%, >52%. These cut-offs divided the data into similar bin size and also roughly corresponded to the traditional definition of isochores [31]. For each distance i, and each G+C content bin, we then calculated the rate of CpG substitution (as the proportion of CpG sites that are mutated). More formally, for each G+C content bin B, we considered:

The results of this analysis are shown in Figure 3A. We observed a clear effect of G+C content of windows close (less than 2,000 bps in each direction) to the target CpG site. Higher G+C content in the window lowers the CpG substitution rate (as expected from the negative relationship between G+C content and CpG rate, reported in the previous section). This effect is the most pronounced at distances very close to the target CpG site. For example, ∼17% of the CpG sites exhibiting low G+C content (GC<38%; red color curve in Figure 3A) at distance 100 bps are mutated, while the same measure for CpG sites with high G+C content (GC>52%; black color curve in Figure 3A) at distance 100 bps is ∼11%. As we move farther away from the CpG site (Figure 3A; left to right along the X axis), the rate of CpG substitution in the low G+C bin (red color curve) and the high G+C bin (black color curve) progressively become closer to each other, displaying the distance-decaying effect of G+C content on the rate of CpG substitution. This effect appears to vanish at around 2,000 bps from the CpG site (i.e., 4,000 bps around the CpG site).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Sliding window analysis of the relationship between CpG substitution rate and G+C content of windows.

(A) At each distance (along the X-axis), CpG sites were divided into four bins based on G+C content at that distance from the site (as measured from the G+C content of the 200 bps window centered at that distance; G+C<38: red curve, 38< = G+C<45%: green curve, 45< = G+C<52%: blue curve, G+C>52%: black curve). The proportion of CpG sites mutated in each of these bins is plotted as a function of distance from the site. At distances closer to the CpG site, the rate of substitution in high local-GC bins (black curve) is clearly lower compared to that of low local-GC bins (red curve). This relationship progressively declines as we move farther away from the site, suggesting a distance-decaying relationship between G+C content and CpG substitution rate. In case of GpC sites, we do not observe a distance-decaying effect (see inset). (B) Results of the chi-square test for the independence of the rate of CpG substitution and the G+C content of the windows at each distance, in log scale. The blue line indicates the P-value cutoff of 0.05 [or log₁₀ (P-value) = −1.30]. The P-values are very low at distances close to the CpG site, and progressively become larger as the distance from the CpG site increases (distance-decaying effect). The rate of CpG substitution becomes independent of the G+C content [log₁₀ (P-value)>−1.30] after ∼2,000 bps from the CpG. (C): Results of the chi-square test for the independence of the rate of GpC substitution and the G+C content of the windows at each distance, in log scale. Again, the blue line indicates the P-value cutoff of 0.05 (or log₁₀ (P-value) = −1.30). The rate of GpC substitution becomes independent of the G+C content [log₁₀ (P-value)>−1.30] at a distance very close to the GpC site, and no distance-decaying effect was observed.

https://doi.org/10.1371/journal.pcbi.1000015.g003

We used a chi-square test to determine whether the rate of CpG substitution is dependent on G+C content, at each distance point (Figure 3B). The distance-decaying relationship is clear when we observe the P-value obtained at each distance. We observe very low P-values (P<<0.05 or log₁₀ (P-value)<<−1.30; which means that the rate of CpG substitution is highly dependent on G+C content) at distances close to the CpG site. As we move farther from the CpG site, the P-value increases progressively, becoming insignificant at ∼2,000 bps, consistent with the results observed in Figure 3A.

We obtained similar results with different overlapping window sizes [for example, window size of 25 bps and step size of 5 bps (Figure S1)], and when 6 bins were used instead of 4 bins (however, using 6 bins decreased the sample sizes and thus increased the fluctuations in the figures). Similar results were obtained with non-overlapping windows (results not shown).

Note that even though there is a distance-decaying relationship between CpG substitution rate and G+C content in immediate neighboring sites, the CpG substitution rate in each GC category differ, although not significantly, from each other even after 2,000 bps. In particular, the CpG substitution rate in low GC bins (red curve in Figure 3A) did not converge with the CpG substitution rate in high GC bins (black curve in Figure 3A) even after 2,000 bps. In the following section, we show that this may be a consequence of long-range (global) nucleotide composition.

As a control for other dinucleotide sites, we performed the sliding window analysis using GpC sites. In contrast to the case of CpG sites, we do not observe a distance-decaying relationship between G+C content and the rate of GpC substitution (GpC curves in Figure 3A). The P-values for testing the independence of the substitution rate at GpC sites and the G+C content was significant at distances >2,000 bps, suggesting a negative relationship between GpC substitution rate and long-range G+C content (see the following section for more details).

The Distance-Decaying Relationship Persists after Correcting for Variation in Global G+C Content

Although the distance-decaying negative relationship between CpG substitution rate and G+C content subsides at ∼2,000 bps (i.e., 4,000 bps around the CpG site; Fig. 3B), the rate of CpG substitution in low G+C content bin continued to remain higher than that of high G+C content bin (Figure 3A). GpC rates also differed between high and low G+C contents (Figure 3A and 3C).

Indeed, both CpG and GpC rates were significantly negatively correlated with G+C contents of 5 kb, 20 kb, and 100 kb blocks around each dinucleotide (results not shown). Also, the G+C contents of 5 kb, 20 kb, and 100 kb blocks are all positively correlated (results not shown), presumably because of the isochore structure of primate genomes. In the remainder of this paper, we refer to the G+C content of 100 kb segments around each CpG site as GC_global, and the negative relationship between GC_global and CpG substitution rate as a “global effect”. We tested if the distance-decaying effect of local (<4,000 bps) G+C content exists even after controlling for global GC effect (by removing variation in GC_global) by the following two analyses.

In the first analysis, we performed the aforementioned sliding window experiment with the G+C content of the windows normalized by the G+C content of the 100 kb segment flanking the CpG site. Precisely, we define the normalized G+C content of a window, denoted GC_norm, as:where GCi denotes the G+C content of the window at distance i from the CpG site and GC_global denotes the G+C content of the 100 kb region surrounding the CpG site. The value GC_norm represents the “relative” G+C content of the i-th window with respect to the G+C content of 100 kb segment within which it belongs.

Figure 4 shows the results of the sliding window analysis, where the G+C contents of the 200 bps windows (GCi) are normalized. The GC_norm cutoffs we used were <0.9, 0.9–1.1, 1.1–1.25, >1.25. These cutoffs were chosen because they divided the data into approximately equal size bins. Even when variation in GC_global was removed, we found the same decaying effect of local G+C content on the rate of CpG substitution (CpG curves in Figure 4A, and 4B). The G+C content of windows closer to the CpG site affected the rate of substitution more compared to that of windows farther away from the CpG sites (CpG curves in Figure 4A and 4B). In contrast to the CpG curves in Figure 3A, the CpG curves in Figure 4A converged, suggesting that the non-convergence in Figure 3A is in fact caused by the global effect. The curves converged at a distance of ∼1,500 bps, which is slightly lower than that observed in Figure 3A. In the case of GpC sites, as expected, there was no distance-decaying relationship and the curves converged at a distance very close to the site (GpC curves in Figure 4A and 4C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Sliding window analysis of the relationship between CpG substitution rate and normalized G+C content.

The same experiment as in Figure 3 with the G+C content of each window normalized with respect to GC_global (removing global effect). (A) The distance-decaying effect of G+C content on the rate of CpG substitution persists even after removing the global effect. In case of GpC substitutions, there was no distance-decaying effect. (B) Results of the chi-square test for the independence of the rate of CpG substitution and the G+C content of the windows. The blue line indicates log₁₀ (P-value) = −1.30. The distance-decaying effect subsided after ∼1,500 bps. (C) Results of the same experiment as in (B), but for GpC sites. There is no distance-decaying effect, as expected.

https://doi.org/10.1371/journal.pcbi.1000015.g004

In our second analysis, we first focused on the distribution of GC_global. GC_global exhibits a bimodal distribution with means of ∼39% and ∼48%, respectively (Figure S2A). A similar distribution of G+C content was observed in case of 100 kb segments surrounding GpC sites (Figure S2B). These bimodal distributions occur because the number of CpG sites and GpC sites are expected to increase with the G+C content of the non-coding region. To reduce variation in GC_global, we analyzed these two distributions (which we call as low- GC_global and high-GC_global regions) separately.

We performed the sliding window analysis to test for the distance-decaying effect of local G+C content in low-GC_global and high-GC_global regions (without normalizing). In the case of CpG sites in low-GC_global regions (Figure S3A and S3B), we obtained similar results as in Figure 4A. In the case of CpG sites in high-GC_global regions, there were large fluctuations in the trend, especially in the low local-GC bins (Figure S4A; red CpG curve). These fluctuations are caused by the reduced sample size in this bin (data not shown). Nevertheless, the distance-decaying effect is clear when we consider the curves of high and low local-G+C content bins (black and red curves CpG, respectively, in Figure S4A). Consistent with the results from Figure 4, the distance from the CpG site at which the curves converged was ∼1,500 bps (Figure S3A). Again, GpC sites did not exhibit any distance-decaying effect (Figures S3A, S3C, S4A, S4C).

Thus, even after removing the global effect, the distance-decaying effect of local G+C content on CpG substitution rate persists, and it subsides after approximately 1,500 bps in each direction from each CpG site. These results were robust to using the flanking 5 kb or 20 kb window for the normalization, different window sizes, and non-overlapping windows (results not shown).

Other Factors That Cause Variation of CpG Substitution Rates

The previous sections analyzed the effect of the efficiency of deamination step on CpG substitution rates. In this section we briefly discuss other factors that cause regional heterogeneity of CpG substitution rates.

For this purpose we consider that a CpG substitution occurs in a three-step process (this is a simplification, but for our purposes this view suffices). First, DNA methylation in mammals (and other warm-blooded vertebrates) occurs specifically at the cytosine bases of CpG dinucleotides. Second, the methylated cytosine undergoes spontaneous deamination. Deaminated methyl-cytosine is identical to thymine, creating a C to T mismatch at the CpG dinucleotide. Third, if this mismatch is left un-repaired, it causes a C to T transition (or a G to A transition in the complementary strand) in the next replication cycle. Based upon such molecular mechanism of the origin of CpG substitution, their rates may vary due to differences in (i) levels of germline methylation, (ii) efficiency of deamination of a methylated cytosine, and (iii) repair efficiency of a transition that occurred via deamination of a methylated cytosine.

A strong support for the first cause, differential methylation, is the observation that transposable elements have generally higher rates of CpG substitutions [19]. Proliferation of transposable elements can potentially impose several deleterious effects on a genome, because insertion of transposable elements can disrupt gene regulation, or cause deleterious recombination events, translocations and other rearrangements [32]. Hence, some organisms may use DNA methylation as a defense mechanism against the proliferation of transposable elements [19],[33].

In our data, we found that transposable elements as a whole exhibit ∼14% higher CpG substitution rate as compared to that of non-repetitive regions (see Text S1). Overall, there is a significant positive correlation between CpG rates and the proportion of transposable elements in a non-coding segment (ρ = 0.142; P<0.001 and ρ_tr = 0.199; P<0.001), consistent with the idea that transposable elements have higher CpG substitution rates as a consequence of increased methylation.

Interestingly, when transposable elements were divided into different classes of elements (LTR, LINE, SINE and DNA elements), SINEs showed similar level of CpG rates to non-repetitive sequences. This was also observed by Meunier et al. [19], who proposed that it might be a consequence of the complex methylation pattern of SINEs in germline cells [34]. However, another possibility is that SINEs show lower CpG rate because their G+C content is in general higher compared to other classes of transposable elements, in light of our finding that there is a negative relationship between CpG substitution rates and G+C contents (see above). SINEs have the highest G+C content among the classes of transposable elements considered in our analyses (50.7%, while 36.6%, 44.0%, 36.4%, for LINEs, LTRs and DNA transposons). This mechanism can at least partially explain why SINEs have lower CpG substitution rates. Because a major proportion of CpG sites in the transposable elements of introns were contributed by SINEs, the rate of CpG substitution in intronic transposable elements is comparable to that of intronic non-repetitive regions (see Text S1).

If the efficiency of the third step, namely mismatch repair, affects regional variation of CpG substitution rates, we expect to see a significant negative correlation between recombination rates and CpG substitution rates. Indeed, we did find a significant negative correlation between recombination rates and the rates of CpG transition substitution (Spearman’s ρ = −0.165, P<0.001). This can be due to the direct effect of biased gene conversion acting on mutations at CpG sites, repairing the TG mismatches (caused by deamination of methyl-cytosine) to correct CG base pairs, the efficiency of which depends on the rate of recombination. Alternatively, biased gene conversion may increase G+C content, which in turn reduces the intensity of DNA melting, leading to a lower rate of substitution at CpG sites. In the first scenario, a significant correlation between CpG substitution rate and recombination rate is expected independent of the effect of G+C content. In the second scenario, the effect of recombination may disappear after correcting for G+C content. In our case, there was no significant correlation after removing the effect of G+C content using partial correlation analysis (ρ_{recombination, CpG rate | GC} = −0.03, P>0.05). However, this question needs to be revisited with more data, most certainly when a higher-quality whole genome alignment of human, chimpanzee, and rhesus macaque becomes available.

General Approach

We can use a parsimony method [19]–[23] to infer CpG sites and determine the substitution rate in humans and chimpanzees using baboon as an outgroup. However, we cannot distinguish substitutions caused by methylation followed by deamination versus replication errors. In primates, CpG substitutions are markedly more frequent than other single nucleotide substitutions [2],[21],[29]. Therefore we assumed methylation followed by deamination is the primary cause of all CpG substitutions in these genomes.

Human, Chimpanzee and Baboon Sequences

We analyzed ∼38 Mbps of sequence data obtained from two different sources; 22 Mbps from sequences orthologous to the human chromosome 7 sequenced by the NISC comparative sequencing group, and ∼15 Mbps obtained from other chromosomes from database mining.

Sequence Annotation and Alignment

Sequences orthologous to human chromosome 7 were aligned using the Threaded Blockset Aligner program [45]. The “best chains” (Text S1) from other chromosomes were aligned using the Multiz program [45]. Non-coding regions (introns and intergenic regions) were identified using gene annotations included in the Known Genes and Ensembl Genes tables of the hg17 assembly of the human genome at UCSC Genome Browser [46]. Intergenic and intronic sequences likely to be selectively constrained [the 5′ and 3′ untranslated regions, first introns and small (<250 bps) introns or intergenic intervals] were excluded. In addition, we excluded alternatively spliced introns based on gene annotations from the UCSC genome browser. In particular, if the span of an intron is different in different transcripts of the same gene, it was removed. The above methods yielded ∼14.7 Mbps of aligned sites in the non-coding regions of the human genome (Table 1). Recombination rates were obtained from [47].

Identification of CpG Islands

CpG islands are regions of the genome where a majority of CpG sites are not methylated. Because methylation of the cytosine in CpG dinucleotides is a prerequisite for the CG->TA mutation to occur (see Introduction), it is crucial for our analyses to exclude CpG sites in CpG islands. Takai and Jones [27] proposed that a good definition of CpG islands is- a region of the genome with (a) G+C content >55%, (b) length >500 bps and (c) observed/expected proportion of CpG dinucleotides (OE) >0.65. The condition length >500 bps was used by Takai and Jones to eliminate the possibility of falsely calling a CpG rich regions generally associated with Alu elements as CpG islands. However, there may be some CpG islands less that 500 bps in length.

In this study, we changed the length constraint to >200 bps to err on the side of caution to safely eliminate most of the CpG islands. We masked out CpG islands identified using the algorithm by Takai and Jones [27] with parameters G+C content >55%, OE>0.65, length >200.

Identification of CpG Sites

To identify CpG sites, we used a parsimony method. Specifically, CpG sites are the middle base of the sites having the following human/chimpanzee/baboon patterns: XNG/XCG/XCG or XCG/XNG/XCG, where X is any nucleotide except G. Given the evolutionary distances considered here, in spite of the hypermutability of CpG sites, such a definition is shown to be quite accurate by a simulation study in Meunier and Duret [22]. GpC sites are sites having the following human/chimpanzee/baboon pattern GNY/GCY/GCY or GCY/GNY/GCY, where Y is any nucleotide except G. The restrictions on X and Y were imposed to avoid overlapping CpG and GpC sites. Non-CpG sites are defined as sites not preceded by a C and not followed by a G. Sites following the complementary patterns of the above definitions were also considered as CpG and GpC and non-CpG sites.

Substitution Rate Estimates and Statistical Tests

To obtain better estimates of substitution rates, we pooled together substitutions in human and chimpanzee lineages. The rate of substitutions for a particular class (CpG or GpC) of site was estimated by dividing the number of substitutions by the total number of sites in that class.

Simulation of Uniform Rate Model

To test if the observed variation in CpG substitution rate in the 50 kb segments is greater than that expected under a uniform substitution rate model, we performed the following simulation. For each segment (s), we kept the number of CpG sites (n_s) the same as that observed in our data and sampled the number of substitutions from the binomial distribution b(n_s,p), where p is the probability of observing a CpG substitution. p was kept constant across segments (0.153 in the case of CpG sites). We simulated 1,000 replicates and found the standard deviation of CpG substitution rate among the segments in each replicate. The range between 2.5% quantile and 97.5% quantile was taken as the 95% confidence interval of the standard deviation under uniform substitution rate model. A similar simulation was also performed for GpC and non-CpG sites.