Skip to main content

Methylation-driven model for analysis of dinucleotide evolution in genomes

Abstract

Background

CpGs, the major methylation sites in vertebrate genomes, exhibit a high mutation rate from the methylated form of CpG to TpG/CpA and, therefore, influence the evolution of genome composition. However, the quantitative effects of CpG to TpG/CpA mutations on the evolution of genome composition in terms of the dinucleotide frequencies/proportions remain poorly understood.

Results

Based on the neutral theory of molecular evolution, we propose a methylation-driven model (MDM) that allows predicting the changes in frequencies/proportions of the 16 dinucleotides and in the GC content of a genome given the known number of CpG to TpG/CpA mutations. The application of MDM to the 10 published vertebrate genomes shows that, for most of the 16 dinucleotides and the GC content, a good consistency is achieved between the predicted and observed trends of changes in the frequencies and content relative to the assumed initial values, and that the model performs better on the mammalian genomes than it does on the lower-vertebrate genomes. The model’s performance depends on the genome composition characteristics, the assumed initial state of the genome, and the estimated parameters, one or more of which are responsible for the different application effects on the mammalian and lower-vertebrate genomes and for the large deviations of the predicted frequencies of a few dinucleotides from their observed frequencies.

Conclusions

Despite certain limitations of the current model, the successful application to the higher-vertebrate (mammalian) genomes witnesses its potential for facilitating studies aimed at understanding the role of methylation in driving the evolution of genome dinucleotide composition.

Background

The k-mer abundance analysis is widely used in genomics research [1,2,3,4,5,6,7,8,9,10]. The term k-mer refers to all possible substrings (in the 5′–3′ direction) of length k in a DNA sequence and, therefore, the k-mer frequency is a good variable for characterizing the composition of a genome’s DNA sequence. Generally, analysis of the 2-mer (i.e., dinucleotide) frequency will provide more abundant information on genome composition than a simple statistic of the 1-mer (i.e., single nucleotides A, C, G, and T) frequency.

It has been shown that in the vertebrate genomes, the CpG dinucleotide is present at a lower frequency than expected [11,12,13,14]. The reason for this is thought to be due to a high C-to-T mutation rate at the methylated CpG sites [14,15,16,17,18]. 5-Methylcytosine (5mc), a cytosine modified by addition of a methyl group on the fifth position of the cytidine ring, can spontaneously deaminate to form thymine (i.e., 5mc to T mutation). Unlike the cytosine (C) to uracil (U) mutation arising from the spontaneous deamination of cytosine, the 5mc to T mutation is rarely recognized and removed by DNA repair enzymes [13]. Therefore, the high level of methylation at the CpG sites explains why the mutation rate of C to T is 10–50 times higher than that of C to other nucleotides [11, 14].

Although several statistical analyses have revealed a negative correlation between the CpG and TpG/CpA levels [19,20,21], to the best of our knowledge, a rigorous method to predict the effects of CpG dinucleotide depletion on the changes in frequencies of all 16 dinucleotides (i.e., ApA, ApC, ApG, ApT, CpA, CpC, CpG, CpT, GpA, GpC, GpG, GpT, TpA, TpC, TpG, and TpT) and in the GC content (GC%) of vertebrate genomes is still lacking. Inspired by the substitution models [22, 23] commonly used in phylogenetic analyses and based on the neutral theory of molecular evolution, in this paper we propose a mathematical model, called the methylation-driven model (MDM), to investigate the effects of the methylation-induced CpG decay on the evolution of the genome dinucleotide composition and GC content.

Results

For the 10 vertebrate genomes, the statistical results of the frequencies/proportions (%) of the 16 dinucleotides and the GC content are listed in Table 1, the corresponding expected values obtained by application of MDM to the initial genome state with 50% GC content and 6.25% proportion of each dinucleotide are listed in Table 2, and the application results for the other two initial genome states, i.e., with 40 and 60% GC contents, are shown in Supplementary Tables 1 and 2 (Additional file 1), respectively.

Table 1 Observed frequencies/proportions of the 16 dinucleotides and GC contents
Table 2 Expected/calculated proportions/frequencies of the 16 dinucleotides and GC contents obtained by MDM (GCini %  = 50%)

Comparison between the expected and observed trends of frequency/proportion changes relative to the initial proportions reveals that, when the initial genomes have a GC content of 50% and proportion for each dinucleotide of 6.25% (see Supplementary Table 3, Additional file 1), most of the 16 dinucleotides in most studied genomes (with the exception of TpA in all 10 genomes, GpA/TpC in cattle and sheep genomes, and ApG/CpT in zebrafish genome; see Tables 1 and 2) have a good consistency between the expected and observed changing trends. This indicates that, on the one hand, 50% GC content could be a rational assumption for the initial state of vertebrate genomes, and on the other hand, our model can achieve a good performance in predicting the changing trends in frequencies of most dinucleotides caused by the methylation-induced CpG to TpG/CpA mutations. It should be noted that, for the dinucleotide TpA, although its proportion/frequency should not be affected by cytidine methylation, the observed frequency in the eight mammalian genomes is either slightly higher or lower than the assumed initial proportion of 6.25% (ranging between 6.01 and 6.59%; see Table 1). Interestingly, in the zebrafish genome, the observed TpA frequency (8.06%) is significantly higher than the assumed initial proportion (6.25%), and in the chicken genome, the observed frequency (5.98%) is lower than those of the eight mammalian genomes.

For each of the 10 tested genomes, the comparison between the observed and expected (obtained with the initial GC content of 50%) frequencies/proportions of the 16 dinucleotides shows an acceptable conformance for most of the 16 dinucleotides (Supplementary Fig. 1, Additional file 1). In order to check whether our model could replicate the observed frequencies/proportions of the 16 dinucleotides, for each of the 10 vertebrate genomes, we have performed the paired t-test of the null hypothesis regarding the differences between the predicted (Table 2) and observed frequencies (Table 1) of the 16 dinucleotides. The results (see Supplementary Fig. 1, Additional file 1) show that all the 10 P-values are close to 1 (> 0.05), indicating that the null hypothesis cannot be rejected at the 95% confidence level and, hence, the differences between the predicted and observed frequencies are statistically acceptable for each tested genome.

In order to evaluate the relative difference between the expected (Proportionexp) and observed (Proportionobs) dinucleotide frequencies/proportions, the value of ΔProportion, defined as the ratio of |(Proportionexp - Proportionobs)| to Proportionobs, was calculated. Figure 1 shows the ΔProportion values of all 16 nucleotides (calculated using Proportionexp obtained from the assumed initial genome state with GC content of 50%) for the 10 vertebrate genomes investigated. For the eight mammalian genomes, most of the 16 dinucleotides have the ΔProportion values either lower than 10% (ApG/CpT, ApT, CpC/GpG, CpG, GpA/TpC, and TpA) or clustering around 20% (ApC/GpT, CpA/TpG, and GpC), indicating a high similarity between the expected and observed frequencies/proportions of them; only two dinucleotides (i.e., ApA/TpT) in few mammal genomes have the ΔProportion values greater than 30%, which indicate relatively large deviations of the expected frequencies from the observed ones. In contrast, the results for the lower-vertebrate (zebrafish) genome do not seem to be satisfactory because several dinucleotides (e.g., CpC/GpG and ApA/TpT) have the expected proportions/frequencies that radically deviate from the observed ones (ΔProportion > 35%). Also worth noting is that the ΔProportion values of some dinucleotides in the lower-vertebrate genomes (e.g., ApA/TpT, ApG/CpT, ApT, CpC/GpG, and TpA in the zebrafish genome; ApT, CpC/GpG, and GpC in the chicken genome) do not cluster around those of the mammal genomes. Overall, the above results indicate that, despite the different application effects on the higher-vertebrate (mammalian) and lower-vertebrate genomes, our model displays a better performance when applied to the mammalian genomes.

Fig. 1
figure 1

Relative differences between the expected and observed proportions of the 16 dinucleotides

Changes in the GC content due to cytidine methylation are only related to the number of CpG to TpG/CpA mutations. The relative difference between the expected (GC%exp) and observed (GC%obs) GC contents was evaluated by calculating the ΔGC% value: ΔGC% = |(GC%exp - GC%obs)| / GC%obs. Figure 2 shows ΔGC% values (calculated using GC%exp obtained from the assumed initial state of a genome with 50% GC content) for the 10 genomes, out of which nine have ΔGC% values less than 10%, indicating a high agreement between the calculated and observed GC contents. The only exception is the zebrafish genome, for which the high ΔGC% value (24.4%) indicates a large deviation of the predicted GC content from the observed one.

Fig. 2
figure 2

Relative differences between the expected and observed GC contents

Discussion

In this paper, we have proposed a mathematical model based on the neutral theory of molecular evolution to analyze the effects of the methylation-induced CpG to TpG/CpA mutations on the evolution of the genome dinucleotide composition and GC content. The model hypothesizes that the neutral mutations (i.e., non-CpG-to-TpG/CpA mutations) would have no effect on the evolution of genome composition. What needs to be highlighted is that the model established here mainly focuses on the methylation-driven evolution of the dinucleotide composition. The previously proposed substitution models, such as the JC69 [22], K80 [23], and TN93 [2] models, despite being widely used in molecular phylogenetic analyses and genome evolution studies, use the rates of the single-nucleotide (1-mer) substitutions as the main parameters. Moreover, these models cannot answer the question of why the observed CpG frequency in the vertebrate genomes is much lower than that expected from the GC content. Although the high rate of the methylation-induced CpG to TpG/CpA mutations can explain the globally reduced frequency of the CpG dinucleotide compared with its expected frequency, there has been a lack of theoretical models with which to predict the mutation effect on the genome composition in terms of 2-mers. Therefore, we propose the MDM through which the methylation-induced changes in frequencies/proportions of the 16 dinucleotides and the GC content can be predicted based on an assumed initial state of a genome.

When modeling the genome evolution, it is inevitable to make assumptions regarding some parameters that cannot be directly obtained. As a matter of fact, the effectiveness of our model largely depends on the validity of the assumed parameters, which in turn depends on the composition characteristics of genomes of interest. For example, for a majority of the 16 dinucleotides and the GC content, their expected/calculated frequencies/content in the zebrafish genome are distinctly different from those in the mammalian genomes (Figs. 1 and 2), and this could be attributed to the large difference in the genome composition between zebrafish and mammals [24]. Specifically, the assumed initial state of the genome can have a large influence on the performance of our model. For example, when the dinucleotide proportions in the initial genome state were assigned based on the assumed GC contents of 40 and 60%, the expected/predicted values for more than half of the 16 dinucleotides and GC contents radically deviate from the observed ones in all the 10 tested genomes (Table 1 and Supplementary Tables 1 and 2, Additional file 1). Even in the case of the “reasonable” initial genome state with the assumed GC content of 50%, the expected proportions/frequencies for ApA/TpT differ from the observed ones by more than 25% in most of the 10 genomes (Fig. 1). We consider that it is the assumed initial proportions of ApA/TpT in the genome that leads to the large deviations from the observed values, while the model itself does not account for these.

It should be noted that there are still shortcomings in MDM regarding its application. Since the model is built based on the neutral mutation theory of molecular evolution, it assumes that the numbers of substitutions between any two single nucleotides are equal (e.g., the numbers of mutations from C to A and from A to C are the same) and, hence, all mutations except for CpG to TpG/CpA will have no effect on the evolution of the genome dinucleotide composition. In fact, mutation bias is prevalent in nature. In Kimura’s two-parameter model [23], rates of transition (α) and transversion (β) substitutions are different and, furthermore, there is no evidence indicating that the CpG to TpG/CpA mutation is the sole factor influencing the rates of transitions and transversions. Ignoring the effects of non-CpG-to-TpG/CpA mutations on the evolution of the genome dinucleotide composition can be expected to produce deviations from the observed values. Moreover, there are limitations in the parameter estimation methods of the presented model, which assume that the rates of the methylation-induced CpG to TpG/CpA mutations are independent of the sequence context of CpG. In fact, the CpG to TpG/CpA mutations are non-neutral [25] and their substitution rates are sequence context dependent [26, 27].

Conclusions

In this work, we have proposed a mathematical model to investigate the effects of the methylation-induced CpG to TpG/CpA mutations on the evolution of genome composition in terms of the 2-mers and GC content. The application of our model to the 10 vertebrate genomes has achieved a good consistency between the predicted and observed trends of changes in the GC content and frequencies of most of the 16 dinucleotides with respect to their assumed initial values; moreover, for the 10 tested genomes, quantitative evaluations of the relative differences in the dinucleotide frequency and GC content between the expected and observed values show a better performance of our model when applied to the mammalian genomes than to the lower vertebrate genomes. Despite the capability of MDM to quantify the effects of the methylation-induced CpG decay on the evolution of the genome dinucleotide composition and GC content, there are still limitations to the current model because i) the rates of the methylation-induced CpG to TpG/CpA mutations are dependent (rather than independent, as assumed in the current model) on the sequence context of CpG sites and, ii) the proportions of the 16 dinucleotides in the initial state of different vertebrate genomes may not be simply identical but depend on the genome composition characteristics. As a result, future efforts in improving the model should be directed toward i) improving the parameter estimation method to make the estimated parameters reflect the context-dependent rates of methylation-induced CpG to TpG/CpA mutations and, ii) realizing the customization of the initial dinucleotide proportions. These two points could be addressed by assigning appropriate weighting coefficients to the parameters as estimated by the Trinucleotide-method and using the genome-composition-based bias factors to calibrate the proportions of the 16 dinucleotides in the initial state of a genome of interest, respectively.

Methods

According to the neutral mutation theory of molecular evolution [28], most evolutionary changes and most of the variability within species at the molecular level are not caused by natural selection, but by random genetic drift of selectively neutral mutations. Since neutral mutations are those that do not affect the survival or reproduction of an organism, they are expected to have a minor effect on the genome dinucleotide composition in the long-term evolutionary process. However, in vertebrate genomes, CpG hypermutability caused by the cytosine methylation clearly exceeds the random expectation. In vertebrate genomes, the observed frequency of CpG dinucleotides is much lower than that expected based on the GC content, implying that the methylation-induced CpG to TpG/CpA mutations exert a substantial effect on the evolution of the genome dinucleotide composition.

MDM proposed in this study focuses on the effects of the methylation-induced CpG to TpG/CpA mutations on the changes in frequencies/proportions of all 16 dinucleotides and in the GC content. The basic hypothesis of MDM is that the cytidine methylation is a key factor that causes the different rates of the transition (α) and transversion (β) substitutions, which under Kimura’s model [23] are considered to drive genome evolution.

Model construction

We assume that in a DNA sequence, each CpG to TpG/CpA substitution is independent from its context. In the case of NpCpG (N = A, C, G, or T), the outcomes of the methylation-induced CpG to TpG mutation on the dinucleotide components can be dissected as follows: the number of CpG is reduced by 1, the number of TpG is increased by 1, and the changes in the numbers of other dinucleotides depend on the nucleotide type of N; for example, if N is A, the number of the dinucleotide ApC will be reduced by 1, while that of ApT will be increased by 1. Let PA, PC, PG, and PT represent the probabilities of N beating A, C, G and T, respectively, then PA + PC + PG + PT = 1. The changes in the numbers of all 16 dinucleotides in the context of NpCpG upon mutation can be described by the matrix D:

$$ \mathrm{D}=\left({d}_{ij}\right)=\left(\begin{array}{cccc}{d}_{AA}& {d}_{AC}& {d}_{AG}& {d}_{AT}\\ {}{d}_{CA}& {d}_{CC}& {d}_{CG}& {d}_{CT}\\ {}{d}_{GA}& {d}_{GC}& {d}_{GG}& {d}_{GT}\\ {}{d}_{TA}& {d}_{TC}& {d}_{TG}& {d}_{TT}\end{array}\right)=\left(\begin{array}{cccc}0& -{P}_A& 0& {P}_A\\ {}0& -{P}_C& -1& {P}_C\\ {}0& -{P}_G& 0& {P}_G\\ {}0& -{P}_T& 1& {P}_T\end{array}\right) $$
(1)

where dij represents the change in the probability of a dinucleotide composed of nucleotides i and j upon NpCpG to NpTpG mutation (i, j = A, C, G, or T), and a negative element, e.g., dAC = −PA, indicates that the number of the corresponding dinucleotide ApC is reduced, and vice versa. Similarly, changes in the numbers of 16 dinucleotides upon CpGpM to CpApM (M = A, C, G, or T) mutation can be represented by the matrix D; where P’A, P’C, P’G, and P’T denote the probabilities of M beating A, C, G, and T, respectively, and P’A + P’C + P’G + P’T = 1.

$$ {\mathrm{D}}^{\prime }=\left({d}_{ij}^{\prime}\right)=\left(\begin{array}{cccc}{d}_{AA}^{\prime }& {d}_{AC}^{\prime }& {d}_{AG}^{\prime }& {d}_{AT}^{\prime}\\ {}{d}_{CA}^{\prime }& {d}_{CC}^{\prime }& {d}_{CG}^{\prime }& {d}_{CT}^{\prime}\\ {}{d}_{GA}^{\prime }& {d}_{GC}^{\prime }& {d}_{GG}^{\prime }& {d}_{GT}^{\prime}\\ {}{d}_{TA}^{\prime }& {d}_{TC}^{\prime }& {d}_{TG}^{\prime }& {d}_{TT}^{\prime}\end{array}\right)=\left(\begin{array}{cccc}P{\prime}_A& P{\prime}_C& P{\prime}_G& P{\prime}_T\\ {}1& 0& -1& 0\\ {}-P{\prime}_A& -P{\prime}_C& -P{\prime}_G& P{\prime}_T\\ {}0& 0& 0& 0\end{array}\right) $$
(2)

The principle of complementary base pairing dictates that the number of NpCpG to NpTpG mutations on the forward strand is the same as that of CpGpM to CpApM mutations on the reverse strand; furthermore, this principle dictates that the CpG to TpG mutations caused by methylation on the reverse strand corresponds to the CpG to CpA mutations on the forward strand, and vice versa. As a result, both the CpG to TpG and CpG to CpA mutations observed on any one of the two strands are the consequence of cytidine methylation of the CpG dinucleotides, although the summation of their numbers doubles the number of the actually depleted CpG dinucleotides. If the number of CpG to TpG/CpA mutations is set to H, the changes in the numbers of the 16 dinucleotides can be calculated by the matrix Q:

$$ \mathrm{Q}=\frac{\mathrm{H}}{2}\left(\mathrm{D}+{\mathrm{D}}^{\prime}\right)=\frac{\mathrm{H}}{2}\left(\begin{array}{cccc}{P}_A^{\prime }& {P}_C^{\prime }-{P}_A& {P}_G^{\prime }& {P}_T^{\prime }+{P}_A\\ {}1& -{P}_C& -2& {P}_C\\ {}-{P}_A^{\prime }& -{P}_C^{\prime }-{P}_G& -{P}_G^{\prime }& {P}_G-{P}_T^{\prime}\\ {}0& -{P}_T& 1& {P}_T\end{array}\right) $$
(3)

The values of the eight parameters PA, PC, PG, PT, P’A, P’C, P’G, and P’T vary between 0 and 1. Before estimating the parameters of the model, we cannot determine the positive or negative values of \( \Big({P}_C^{\prime }-{P}_A \)) and (\( {P}_G-{P}_T^{\prime } \)). However, it is clear from the matrix Q that for 14 of the 16 dinucleotides (with the exception of ApC (P’C - PA) and GpT (PG - P’T)), their changes in number can be directly inferred without requiring parameter estimation. Table 3 lists the changing trends in the numbers of the 14 dinucleotides upon CpG to TpG/CpA mutations inferred directly from the matrix Q.

Table 3 CpG to TpG/CpA mutation-caused changing trends in the numbers of dinucleotides

Parameter estimation

In the matrix Q, the values for the parameters PA, PC, PG, PT, P’A, P’C, P’G, and P’T are unknown. To determine the change in the number of each dinucleotide as a function of H (i.e., the number of CpG to TpG/CpA mutations), each parameter value should be estimated. Based on the assumption that the rate of the methylation-induced CpG-to-TpG/CpA mutations is independent of the sequence context of CpG sites, we propose two methods for parameter estimation.

Trinucleotide-method: Using a simple ratio approach, the parameters PA, PC, PG, and PT can be calculated as the proportion of ApCpG, CpCpG, GpCpG, and TpCpG among all NpCpG trinucleotides, respectively. Accordingly, P’A, P’C, P’G, and P’T can be obtained through calculating the proportion of CpGpA, CpGpC, CpGpG, and CpGpT among all CpGpM trinucleotides, respectively. Since most CpG islands remain unmethylated in normal cells [29, 30], they are excluded from parameter estimation. In addition, the CpG sites located within the coding regions could also be excluded due to the high selection pressure on these regions. Nevertheless, since the proportion of CpGs in the coding regions out of the total CpGs is small, we expect that the inclusion of the coding-region CpGs would have a negligible effect on the estimation results. Through counting all the eight trinucleotides in the context of NpCpG and CpGpM, the parameters can be calculated by eqs. (4) and (5):

$$ \left.\begin{array}{c}{P}_A=\frac{S_{\mathrm{ACG}}}{S_{\mathrm{ACG}}+{S}_{\mathrm{TCG}}+{S}_{\mathrm{CCG}}+{S}_{\mathrm{GCG}}}\\ {}{P}_C=\frac{S_{\mathrm{CCG}}}{S_{\mathrm{ACG}}+{S}_{\mathrm{TCG}}+{S}_{\mathrm{CCG}}+{S}_{\mathrm{GCG}}}\\ {}{P}_G=\frac{S_{\mathrm{GCG}}}{S_{\mathrm{ACG}}+{S}_{\mathrm{TCG}}+{S}_{\mathrm{CCG}}+{S}_{\mathrm{GCG}}}\\ {}{P}_T=\frac{S_{\mathrm{TCG}}}{S_{\mathrm{ACG}}+{S}_{\mathrm{TCG}}+{S}_{\mathrm{CCG}}+{S}_{\mathrm{GCG}}}\end{array}\right\} $$
(4)
$$ \left.\begin{array}{c}P{\prime}_A=\frac{S_{\mathrm{CGA}}}{S_{\mathrm{CGA}}+{S}_{\mathrm{CGT}}+{S}_{\mathrm{CGC}}+{S}_{\mathrm{CGG}}}\\ {}P{\prime}_C=\frac{S_{\mathrm{CGC}}}{S_{\mathrm{CGA}}+{S}_{\mathrm{CGT}}+{S}_{\mathrm{CGC}}+{S}_{\mathrm{CGG}}}\\ {}P{\prime}_G=\frac{S_{\mathrm{CGG}}}{S_{\mathrm{CGA}}+{S}_{\mathrm{CGT}}+{S}_{\mathrm{CGC}}+{S}_{\mathrm{CGG}}}\\ {}P{\prime}_T=\frac{S_{\mathrm{CGT}}}{S_{\mathrm{CGA}}+{S}_{\mathrm{CGT}}+{S}_{\mathrm{CGC}}+{S}_{\mathrm{CGG}}}\end{array}\right\} $$
(5)

where SACG, SCCG, SGCG, and STCG are the numbers of ApCpG, CpCpG, GpCpG, and TpCpG, respectively, and SCGA, SCGC, SCGG, and SCGT are the numbers of CpGpA, CpGpC, CpGpG, and CpGpT, respectively, in a genome sequence.

GC-method: Because of the sequence complementarity between the forward and reverse strands of DNA [31], the parameters have the following relationships: PA ≈ P’T, PC ≈ P’G, PG ≈ P’C, and PT ≈ P’A. Therefore, the eight parameters can be estimated based on the GC content in a genome. If the GC content is p, then the AT content is (1 – p), and the following relations can be obtained:

$$ \left.\begin{array}{c}{P}_A=\frac{1-\mathrm{p}}{2}\\ {}{P}_C=\frac{\mathrm{p}}{2}\ \\ {}{P}_G=\frac{\mathrm{p}}{2}\ \\ {}{P}_T=\frac{1-\mathrm{p}}{2}\end{array}\right\} $$
(6)
$$ \left.\begin{array}{c}P{\prime}_A\approx \frac{1-\mathrm{p}}{2}\\ {}P{\prime}_C\approx \frac{\mathrm{p}}{2}\ \\ {}P{\prime}_G\approx \frac{\mathrm{p}}{2}\ \\ {}P{\prime}_T\approx \frac{1-\mathrm{p}}{2}\end{array}\right\} $$
(7)

For example, if the GC content is 40%, then PA = PT = P’T = P’A = 0.3 and PC = PG = P’C = P’G = 0.2.

It should be noted that in the practical application of MDM, the trinucleotide-method should be preferred over the GC-method. In the application example below, the parameters are estimated based on the current state of a genome and are assumed to be constant during the genome evolution. In fact, the cumulative mutations of CpG to TpG/CpA will inevitably lead to a reduction in the GC content. Therefore, if the GC-method is used, the four estimated parameters (i.e., PC, PG, P’C, and P’G) will be prone to errors due to the large changes in the GC content between the assumed initial state and current state of a genome.

Application example

The constructed MDM was applied to predict the effects of CpG to TpG/CpA mutations on the evolution of the dinucleotide composition in vertebrate genomes starting from three assumed initial states. The source code used is publicly available at https://github.com/sparkhonghe/MDM.

Ten vertebrate genomes, of which eight are from mammals (higher vertebrates) and two from non-mammals (lower vertebrates, i.e., Gallus gallus (chicken) and Danio rerio (zebrafish)) (see Supplementary Table 4, Additional file 1), were obtained from the genome database of the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov). To avoid artifacts arising from sex-specific effects, only the sequence from autosomes was included in the statistics of the frequencies/proportions (%) of the 16 dinucleotides and the GC content in these 10 genomes.

There are three basic assumptions in application of MDM: i) the loss of CpG dinucleotides in a vertebrate genome is exclusively caused by cytidine methylation; ii) the rate of the methylation-induced CpG to TpG/CpA mutations is independent of the sequence context of CpG sites; and iii) the total number of dinucleotides in the genome remains constant during evolution. In addition, the model’s application also needs to make assumptions about the proportions/frequencies of the 16 dinucleotides in the initial state of a genome, which can be obtained based on the assumed initial GC content (GCini%). In this study, three assumed initial genome states, i.e., in which the proportions of the 16 dinucleotides were obtained (see eq. (8) below) according to the GCini% of 40% (see Supplementary Table 5, Additional file 1), 50% (Supplementary Table 3, Additional file 1), and 60% (Supplementary Table 6, Additional file 1), respectively, were tested. It should be noted that, under the same GC content, all the 10 genomes possess the same initial state in terms of the frequency/proportion of each dinucleotide.

Given the GC content, the frequencies/proportions of the 16 dinucleotides (NpMini%) in the assumed initial state of a genome can be calculated using the following equation:

$$ {\mathrm{N}\mathrm{pM}}_{\mathrm{ini}}\%=\left({\mathrm{N}}_{\mathrm{num}}\ast {\mathrm{M}}_{\mathrm{num}}\right)/{\mathrm{L}}^2 $$
(8)

where M and N = A, C, G, or T, Nnum and Mnum represent the numbers of nucleotides N and M in the initial genome, respectively, and L is the genome length. As a result, Nnum/L and Mnum/L designate the frequencies of N and M, respectively. For example, assuming that the GC content is 50% and genome length is 100, then in the initial state of this genome, the number of each of the four nucleotides (A, C, G and T) is 25, and the frequency of each of the 16 dinucleotides (NpMini%) is 25 * 25/1002 = 6.25%. Note that the eq. (8) can also be used for estimating the expected frequency of CpG in the current state of a genome if we take Nnum/L and Mnum / L as the observed frequencies of C and G (i.e., one-half of the observed GC content), respectively, in the genome.

As rationalized above, the GC content can vary during the evolution of a genome due to continuous CpG depletion. Therefore, we adopt the “Trinucleotide-method” to estimate the parameters PA, PC, PG, PT, P’A, P’C, P’G, and P’T in the matrix Q. The total number of the depleted CpG dinucleotides (H) in a genome can be calculated using the following equations:

$$ \mathrm{H}=\left(\mathrm{CpG}{\%}_{\mathrm{ini}}-\mathrm{CpG}{\%}_{\mathrm{obs}}\right)\ast {\mathrm{N}}_{\mathrm{N}\mathrm{pM}} $$
(9)
$$ {\mathrm{N}}_{\mathrm{N}\mathrm{pM}}=\mathrm{L}-1 $$
(10)

where CpG%ini and CpG%obs is the assumed initial proportion (see Supplementary Tables 3, 5, and 6, Additional file 1) and the observed frequency of CpG (Table 1) in the genome, respectively, NNpM is the total number of dinucleotides in the genome, and L is the genome length after removing gaps (Supplementary Table 4, Additional file 1).

The CpG island, which is defined as a stretch of DNA sequence with length > 200 bp, GC content > 50%, and observed-to-expected CpG ratio > 0.6 [13, 32], was identified using the CpG Island Searcher [33]. After removing CpG islands, NpCpG and CpGpM trinucleotides in each of the 10 vertebrate genomes were counted using an in-house Java program (for results, see Supplementary Table 7, Additional file 1), and the eight parameters were then obtained with eqs. (4) and (5). Table 4 lists the estimated parameters for all the 10 vertebrate genomes.

Table 4 Parameters estimated by statistics of the trinucleotides NpCpG and CpGpM

Once the eight parameters (Table 4) were obtained and the number of CpG to TpG/CpA mutations (H) was determined (eq. (9)), the changes in the numbers of the 16 dinucleotides in a genome from the assumed initial state to the current state can be derived from the matrix Q (eq. (3)). Finally, for each of the 16 dinucleotides, its predicted/expected number in the current genome state was obtained by adding the predicted number of changes to the assumed initial number, followed by converting to the proportion of the total number of all 16 dinucleotides to facilitate comparison.

The expected GC content (GC%exp) in a genome was obtained using the following equation:

$$ \mathrm{GC}{\%}_{\mathrm{exp}}={\mathrm{GC}}_{\mathrm{ini}}\%-\mathrm{H}/\mathrm{L} $$
(11)

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Abbreviations

MDM:

Methylation-driven model

References

  1. Tamura K. The rate and pattern of nucleotide substitution in Drosophila mitochondrial DNA. Mol Biol Evol. 1992;9:814–25.

    CAS  PubMed  Google Scholar 

  2. Tamura K, Nei M. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol Biol Evol. 1993;10:512–26.

    CAS  Google Scholar 

  3. Compeau PE, Pevzner PA, Tesler G. How to apply de Bruijn graphs to genome assembly. Nat Biotechnol. 2011;29:987–91.

    Article  CAS  Google Scholar 

  4. Dubinkina VB, Ischenko DS, Ulyantsev VI, Tyakht AV, Alexeev DG. Assessment of k-mer spectrum applicability for metagenomic dissimilarity analysis. BMC Bioinformatics. 2016;17:38.

    Article  Google Scholar 

  5. Fiannaca A, La Rosa M, Rizzo R, Urso A. A k-mer-based barcode DNA classification methodology based on spectral representation and a neural gas network. Artif Intell Med. 2015;64:173–84.

    Article  Google Scholar 

  6. Mohamed Hashim EK, Abdullah R. Rare k-mer DNA: identification of sequence motifs and prediction of CpG Island and promoter. J Theor Biol. 2015;387:88–100.

    Article  CAS  Google Scholar 

  7. Lee D, Karchin R, Beer MA. Discriminative prediction of mammalian enhancers from DNA sequence. Genome Res. 2011;21:2167–80.

    Article  CAS  Google Scholar 

  8. Meher PK, Sahu TK, Rao AR. Identification of species based on DNA barcode using k-mer feature vector and random forest classifier. Gene. 2016;592:316–24.

    Article  CAS  Google Scholar 

  9. Ounit R, Wanamaker S, Close TJ, Lonardi S. CLARK: fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers. BMC Genomics. 2015;16:236.

    Article  Google Scholar 

  10. Wang R, Xu Y, Liu B. Corrigendum: recombination spot identification based on gapped k-mers. Sci Rep. 2016;6:35331.

    Article  CAS  Google Scholar 

  11. Antequera F, Bird A. Number of CpG islands and genes in human and mouse. Proc Natl Acad Sci U S A. 1993;90:11995–9.

    Article  CAS  Google Scholar 

  12. Furano AV, Walser JC. Mutation rate of non-CpG DNA. In: eLS. Chichester: Wiley; 2009. https://doi.org/10.1002/9780470015902.a0021740.

    Book  Google Scholar 

  13. Gardiner-Garden M, Frommer M. CpG Islands in vertebrate genomes. J Mol Biol. 1987;196:261–82.

    Article  CAS  Google Scholar 

  14. Ioshikhes IP, Zhang MQ. Large-scale human promoter mapping using CpG islands. Nat Genet. 2000;26:61–3.

    Article  CAS  Google Scholar 

  15. Bird A. DNA methylation de novo. Science. 1999;286:2287–8.

    Article  CAS  Google Scholar 

  16. Duret L, Galtier N. The covariation between TpA deficiency, CpG deficiency, and G+C content of human isochores is due to a mathematical artifact. Mol Biol Evol. 2000;17:1620–5.

    Article  CAS  Google Scholar 

  17. Saxonov S, Berg P, Brutlag DL. A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters. Proc Natl Acad Sci U S A. 2006;103:1412–7.

    Article  CAS  Google Scholar 

  18. Scarano E, Iaccarino M, Grippo P, Parisi E. The heterogeneity of thymine methyl group origin in DNA pyrimidine isostichs of developing sea urchin embryos. Proc Natl Acad Sci U S A. 1967;57:1394–400.

    Article  CAS  Google Scholar 

  19. Jabbari K, Bernardi G. Cytosine methylation and CpG, TpG (CpA) and TpA frequencies. Gene. 2004;333:143–9.

    Article  CAS  Google Scholar 

  20. Upadhyay M, Samal J, Kandpal M, Vasaikar S, Biswas B, Gomes J, et al. CpG dinucleotide frequencies reveal the role of host methylation capabilities in parvovirus evolution. J Virol. 2013;87:13816–24.

    Article  CAS  Google Scholar 

  21. Xiang S, Liu Z, Zhang B, Zhou J, Zhu BD, Ji J, et al. Methylation status of individual CpG sites within Alu elements in the human genome and Alu hypomethylation in gastric carcinomas. BMC Cancer. 2010;10:44.

    Article  Google Scholar 

  22. Jukes TH, Cantor CR. Evolution of protein molecules. In: Munro HN, editor. Mammalian protein metabolism. New York: Academic; 1969. p. 21–132.

    Chapter  Google Scholar 

  23. Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111–20.

    Article  CAS  Google Scholar 

  24. Howe K, Clark MD, Torroja CF, Torrance J, Berthelot C, Muffato M, et al. The zebrafish reference genome sequence and its relationship to the human genome. Nature. 2013;496:498–503.

    Article  CAS  Google Scholar 

  25. Schmidt S, Gerasimova A, Kondrashov FA, Adzhubei IA, Kondrashov AS, Sunyaev S. Hypermutable non-synonymous sites are under stronger negative selection. PLoS Genet. 2008;4:e1000281.

    Article  Google Scholar 

  26. Mugal CF, Ellegren H. Substitution rate variation at human CpG sites correlates with non-CpG divergence, methylation level and GC content. Genome Biol. 2011;12:R58.

    Article  CAS  Google Scholar 

  27. Aggarwala V, Voight BF. An expanded sequence context model broadly explains variability in polymorphism levels across the human genome. Nat Genet. 2016;48:349–55.

    Article  CAS  Google Scholar 

  28. Kimura M. The neutral theory of molecular evolution: a review of recent evidence. Jpn J Genet. 1991;66:367–86.

    Article  CAS  Google Scholar 

  29. Bird AP. CpG-rich islands and the function of DNA methylation. Nature. 1986;321:209–13.

    Article  CAS  Google Scholar 

  30. Dunham I, Shimizu N, Roe BA, Chissoe S, Hunt AR, Collins JE, et al. The DNA sequence of human chromosome 22. Nature. 1999;402:489–95.

    Article  CAS  Google Scholar 

  31. Sueoka N. Two aspects of DNA base composition: G+C content and translation-coupled deviation from intra-strand rule of a = T and G = C. J Mol Evol. 1999;49:49–62.

    Article  CAS  Google Scholar 

  32. Takai D, Jones PA. Comprehensive analysis of CpG islands in human chromosomes 21 and 22. Proc Natl Acad Sci U S A. 2002;99:3740–5.

    Article  CAS  Google Scholar 

  33. Takai D, Jones PA. The CpG island searcher: a new WWW resource. In Silico Biol. 2003;3:235–40.

    CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable

Funding

This research has been funded by the National Natural Sciences Foundation of China (Nos. 31370715 and 31160181) and Programs for Excellent Young Talents (XT412003) and Donglu Scholar in Yunnan University.

Author information

Authors and Affiliations

Authors

Contributions

JHS and SQL conceived and designed the study. JHS wrote the paper. JHS, SMA and SQL reviewed and edited the manuscript. All authors read and approved the manuscript.

Corresponding author

Correspondence to Shu-Qun Liu.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Supplementary Table 1

. Expected/calculated proportions/frequencies of the 16 dinucleotides and GC contents obtained by MDM (GCini% = 40%). Supplementary Table 2. Expected/calculated proportions/frequencies of the 16 dinucleotides and GC contents obtained by MDM (GCini% = 60%). Supplementary Table 3. Proportions/frequencies of the 16 dinucleotides in the assumed initial state of genomes with GCini% = 50%. Supplementary Table 4. Information of the 10 vertebrate genomes. Supplementary Table 5. Proportions/frequencies of the 16 dinucleotides in the assumed initial state of genomes with GCini% = 40%. Supplementary Table 6. Proportions/frequencies of the 16 dinucleotides in the assumed initial state of genomes with GCini% = 60%. Supplementary Table 7. Numbers of the trinucleotides NpCpG and CpGpM in the 10 vertebrate genomes. Supplementary Fig. 1. Comparison between the observed and expected frequencies/proportions of the 16 dinucleotides. Note that the expected frequencies were obtained using GCini% = 50%. P-value shown in the inserted box was obtained by performing the paired t-test on the observed and expected frequencies of the 16 dinucleotides for each genome.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, JH., Ai, SM. & Liu, SQ. Methylation-driven model for analysis of dinucleotide evolution in genomes. Theor Biol Med Model 17, 3 (2020). https://doi.org/10.1186/s12976-020-00122-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12976-020-00122-x

Keywords