BACKGROUND
the evolution of the coding exome is a major driving force of functional divergence.
in the past, a coding gene was considered as a basic unit for biological regulations and molecular functions.
as such, in the majority of evolutionary studies, the “functional unit” targeted by natural selection is presumed to be a gene.
however, with the advances in molecular biology and high-throughput sequencing technologies, it has gradually become clear that alternative transcript isoforms of the same gene  can be spatio-temporally regulated, and convey fairly divergent biological functions  <cit> .
in other words, in many cases, a “transcript” rather than a “gene” is the biologically functional unit.
the importance of transcript isoforms is particularly significant in complex organisms because they have highly developed networks of transcript/protein isoforms  <cit> .
transcript isoforms of the same gene differ from each other by alternatively spliced exonic regions.
in cases where transcript isoforms convey distinct biological functions, the alternatively spliced exonic regions are crucial for the between-isoform functional divergences.
these exonic regions should be accordingly targeted by natural selection.
therefore, the biological functions of alternative  exonic sequences and the selection pressure thereon can be revealed by examining the evolutionary patterns of these sequences  <cit> .
we previously examined the determinants of exonic evolutionary rates in mammals and arabidopsis species.
the biological factors that affect exonic evolutionary rates were found to differ between these two lineages  <cit> .
in addition, we discovered that in mammals, the position of an exon  in the transcript is significantly associated with the evolution of the exonic sequence in accordance with the level of dna methylation at cpg dinucleotides   <cit> .
this is probably because the position of an exon is related to its biological function , thus making the exon selectively constrained for the function mediated by cpg methylation, or prone to the mutagenesis effect of cpg methylation  <cit> .
however, whether this proposition is also true for arabidopsis remains unexplored.
plant coding exons differ from their mammalian counterparts in several aspects.
firstly, alternative rna splicing is less well developed, and plays a less important role in exon evolution in plants than in mammals  <cit> .
secondly, on average, a plant gene includes fewer but longer exons than a mammalian gene  <cit> .
thirdly, the effective population sizes of plants  are considerably larger than those of mammals   <cit> , giving rise to a higher efficiency of natural selection on plant exonic sequences.
given these differences, we expect the evolutionary patterns of arabidopsis exons at different positions to diverge from those of their mammalian counterparts.
in this study, we systematically examined the mutational effects and of cpg methylation and its correlations with exonic evolutionary rates for a. thaliana coding exons at different positions.
our results indicate that first, last, and internal coding exons of a. thaliana have fairly different evolutionary patterns in this regard.
the three exon groups diverge significantly in their liability to cpg methylation-related mutagenesis.
furthermore, the cpg methylation-evolutionary rate correlations differ significantly among the three exon groups.
these correlations also differ significantly between arabidopsis species and mammals.
in addition, we found exon length to be a cpg methylation-independent indicator of exonic evolutionary rates in arabidopsis species.
our results suggest that intrinsic exonic features  may be important determinants for the evolution of a. thaliana coding exons, and that cpg methylation may play different biological roles in the coding exons of mammals and arabidopsis species.
RESULTS
the mutagenic effect of cpg methylation for exons at different positions
to examine the mutagenic effect of cpg methylation on a. thaliana coding exons, we calculated the pearson’s coefficient of correlation between the level of cpg methylation  and the cpg o/e ratio .
of note, here we do not include methylation at chg or chh sites  because these two types of methylation account for a minority of the plant methylome  <cit> , and may have minor effects on the evolution of arabidopsis exons.
cpg methylation can significantly increase the rate of cytosine-to-thymine  transitions, leading to a decreased number of cpg dinucleotides.
therefore, mcg density is expected to be negatively correlated with the cpg o/e ratio.
furthermore, a larger absolute value of the coefficient of correlation  indicates a stronger mutagenic effect of dna methylation  <cit> .
here the pearson’s  correlation is employed to show this quantitative relationship.
of note, when dealing with the methylome data, we applied a set of filtering criteria to ensure data quality and to reduce variations in the estimation of mcg density .
these filtering criteria lead to differences in the number of analyzable exons among the four methylome datasets .
interestingly, as shown in figure 1a, despite the variations in the number of analyzable exons and the level of cpg methylation, mcg density is consistently negatively correlated with the cpg o/e ratio across datasets.
and the r value falls within a relatively narrow range .
this observation confirms the mutagenic effect of cpg methylation on a. thaliana coding exons.
note that the background dataset was not filtered for methylation data, but was filtered for the definitions of first, last, and internal exons and for the length threshold for calculations of evolutionary rates .
na: not applicable.
note that the “bisulfite seq.
read depth” is defined as the total length of the bisulfite sequencing reads divided by the size of the a. thaliana genome.
the next question to ask is whether the mutagenic effect of cpg methylation differs for first, last, and internal exons.
to address this issue, we evaluated the pearson’s correlations as described above separately for each of the three exon groups.
unexpectedly, as shown in figure 1b, although the mcg-cpg o/e correlations remain negative across the three exon groups, the strongest and the weakest mutagenic effect occur, respectively, in internal and first exons.
this is to the contrary of what was previously observed for mammals, where the strongest mutagenic effect of cpg methylation occurs in first exons, while the weakest in internal exons  <cit> .
we then ask how the variation in mcg-related mutagenic effect may have affected the evolutionary rates of the three exon groups.
to this end, we first examined the evolutionary rates  separately for the three exon groups between a. thaliana and a. lyrata based on datasets s1 ~ s <dig> .
the four sperm methylome datasets used here are appropriate for this evolutionary analysis because only the mutations that occur in germ line cells can be propagated over generations, thus leaving observable changes in the inter-species comparison.
figure  <dig> shows that among the three exon groups, first exons have the largest median dn and dn/ds ratio, followed by last exons, and finally by internal exons .
for ds, internal exons have the lowest median value, and first and last exons have similar values.
the evolutionary rate profiles apparently are inconsistent with the mcg-related mutagenic effect profiles in first, last, and internal exons.
specifically, in internal exons, we observe the co-occurrence of a high median mcg density , a strong mutagenic effect of cpg methylation , and low evolutionary rates as compared with the other two exon groups .
one possible explanation is that the strong selection pressure imposed on internal exons has significantly constrained the mcg-related mutations from occurring in this exon group.
this appears to be true judging from the low dn/ds ratio in internal exons as compared with the other two exon groups .
interestingly, the median ds is also the lowest in internal exons, suggesting that synonymous substitutions are subject to strong purifying selection in this exon group.
the correlations between cpg methylation and exonic evolutionary rates
to investigate in more details the correlations between cpg methylation and exonic evolutionary rates, we evaluated spearman’s correlations between mcg density and each of the three evolutionary measurements  separately for the four sperm methylome datasets.
figure  <dig> shows that for first exons, mcg density has weak positive correlations with dn, ds and the dn/ds ratio.
this observation suggests that the mutagenic effect of cpg methylation influences both synonymous and nonsynonymous sites in first exons, but the effects may be relatively small.
for last exons, the spearman’s coefficients of correlation  are consistently negative between mcg density and each of dn, ds, and the dn/ds ratio across datasets.
this observation is somewhat surprising because a negative mcg density-dn  correlation indicates that a higher level of cpg methylation is accompanied with a reduced rate of sequence evolution.
in mammals, by comparison, the mcg density-ds correlation is positive, while the mcg density-dn and mcg density-dn/ds correlations are negative in last exons  <cit> .
meanwhile, for internal exons, the mcg density-ds correlations are consistently negative across datasets.
however, the mcg density-dn correlations are weakly negative in this exon group.
unexpectedly, the mcg density-dn/ds ratio correlations are positive in internal exons.
these observations appear to suggest that the synonymous sites in internal exons are subject to certain mcg-associated selective constraints, thus heavily methylated exons tend to have lower ds values.
however, such selective constraints may be weaker for the nonsynonymous sites in the same exonic regions.
the decrease in ds and the relative stasis in dn lead to the increase in the dn/ds ratio as mcg density increases in internal exons .
these results also differ from what were reported for mammals, where the mcg density-ds correlation is positive, and the mcg- dn and mcg density-dn/ds ratio correlations are both negative in internal exons  <cit> .
the differences in the mcg density-evolutionary rate correlations between a. thaliana and mammals suggest that the biological roles of cpg methylation may have diverged between the two lineages.
we also conducted the same analyses while controlling for four potential confounding factors  or alternative exon  type)  <cit> .
the results remain virtually the same .
the associations between exon length and exonic sequence evolution
we have shown that the mcg density-cpg o/e correlation is the most marked in internal exons.
one important question is whether this observation has actually resulted from certain sampling biases.
we noticed that our filtering criteria for dealing with the methylome data  tend to retain longer exons in the datasets.
the average lengths of the selected exons in s1 ~ s <dig> range from  <dig>  to  <dig>  base pairs , which are considerably longer than the previously reported exome-wide average of ~250 bp in a. thaliana <cit>  and the background value  in this study .
to investigate whether exon length has affected our results, we divided each of the methylome datasets  further into five subgroups of approximately equal sizes according to exon length, and evaluated the mcg density-cpg o/e correlations separately for first, last, and internal exons for each length subgroup.
surprisingly, as shown in figure  <dig>  first exons show a clear length-dependent decrease in the r value between mcg density and the cpg o/e ratio, indicating stronger mutagenic effects of cpg methylation on longer first exons.
by comparison, internal and last exons do not show similarly clear trends.
since exon length is related to the mutagenic effect of cpg methylation, we then ask whether exon length is correlated with the evolutionary rates of first, last, and internal exons.
accordingly, we evaluated spearman’s correlations between exon length and dn, ds, and the dn/ds ratio separately for the three exon groups.
of note, here we take the background dataset  instead of the methylome datasets to avoid potential sampling biases.
interestingly, table  <dig>  shows that for first exons, exon length is significantly correlated with both dn and ds, although the ρ value is smaller in the exon length-dn correlation.
meanwhile, the exon length-dn/ds ratio correlation is statistically insignificant for this exon group.
this observation seems to suggest that the increased mutagenic effect of cpg methylation  in longer exons has led to increases in both dn and ds in this exon group.
this conjecture, however, is only partly true, as will be shown later.
meanwhile, for last and internal exons, the correlations between exon length and the three evolutionary measurements are all statistically significant.
yet a noteworthy difference between last and internal exons is that the ρ value of the exon length-dn correlation is about four times larger than that of the exon length-ds correlation for internal exons.
this difference in ρ value is smaller in last exons.
d
n
/
d
s
ratio,
d
n
, and
d
s
before  and after  controlling for four potential confounding factors
note that this table is based on the background dataset.
the above observations may be confounded by other biological factors.
for example, ases are known to have increased dn and the dn/ds ratios as compared with cses  <cit> .
therefore, the increase in dn and the dn/ds ratio in longer exons might have resulted from an increase in the proportion of ases.
meanwhile, the proportion of repetitive elements  is also correlated with evolutionary rates because these elements are subject to relaxed selective constraints  <cit> .
similar comments also apply to intrinsically disordered protein regions  <cit> .
the next factor to consider is exonic expression level, which has been shown to be an important determinant of dn and the dn/ds ratio  <cit> .
we thus conducted partial spearman’s correlation analyses while simultaneously controlling all of these four factors .
as shown in table  <dig> , the results remain virtually the same.
the last but a critical factor to control is the level of cpg methylation.
to evaluate the influence of cpg methylation, we have to employ the sperm methylome datasets , which include considerably fewer but longer exons as compared with the background dataset .
using the four sperm methylome datasets, we again conducted partial spearman’s correlation analysis between exon length and evolutionary rates while simultaneously controlling for mcg density, ase/cse exon type, proportion of repetitive elements/disordered region, and exonic expression level.
as shown in additional file  <dig>  for first exons, the results are similar but the ρ values are decreased.
this observation indicates that the length dependence of mcg-related mutagenic effect in first exons  accounts for part but not all of the length dependence of dn and ds in this exon group.
meanwhile, for last exons, the exon length-ds correlation becomes statistically insignificant in s <dig> and s <dig>  whereas the exon length-dn and the exon length-dn/ds correlations remain statistically significant with decreased ρ values.
these results imply that mcg-related mutations may account for part of the length dependence of ds, dn, and the dn/ds ratio in last exons.
however, the decreases in ρ value and the level of statistical significance may also be ascribable to the decrease in sample size and the bias in exon length.
by comparison, for internal exons, all of the correlations remain statistically significant with two notable changes as compared with the results in table 2:  the ρ values of the exon length-dn and the exon length-dn/ds correlations are increased; and  the ρ values of the exon length-ds correlations turn negative.
therefore, for internal exons, mcg-related mutations appear to be an important factor affecting ds.
nevertheless, mcg-related mutations cannot explain the length dependence of dn and the dn/ds ratio in this exon groups.
taken together, our results indicate that the correlations between exon length and the evolutionary measurements  are unaffected by the ase/cse exon type, proportion of repetitive elements/disordered region, and exonic expression level in any of the three exon groups.
however, the level of cpg methylation may account for part of the exon length-evolutionary rate correlations differentially for first, last, and internal exons.
in summary, exon length appears to be a cpg methylation-independent indicator for dn in all of the three exon groups, and for the dn/ds ratio in last and internal exons of a. thaliana.
discussion
we have shown that for the coding sequences of a. thaliana, the mutagenic effects of cpg methylation differ between exons at different relative positions.
among the three compared exon groups , the highest cpg methylation level and the strongest mutagenic effect of cpg methylation both occur in internal coding exons  despite the most stringent selective constraint  on this exon group .
first coding exons, quite to the opposite, have the lowest level of cpg methylation and suffer the weakest mutagenic effect of cpg methylation, yet evolve the most rapidly.
interestingly, we show that mcg density is  positively correlated with ds, dn, and dn/ds ratio in first exons, yet the same correlations are significantly negative for last exons.
for internal exons, the correlations are negative, weakly negative, and positive for ds, dn, and dn/ds ratio, respectively .
the mutagenic effect of cpg methylation cannot fully explain these observations.
apparently, selection pressure has played a major role here.
we have previously reported that in mammals, cpg methylation may have different biological roles in first, last, and internal coding exons  <cit> .
similar comments may also apply to arabidopsis species – that first exons are more liable to the mutagenic effects, yet the other two exon groups are more affected by the regulatory functions of cpg methylation.
noticeably, however, the correlations between mcg density and evolutionary rates actually diverge significantly between arabidopsis species and mammals  <cit> .
one riveting difference is that for internal exons, the mcg density-ds and mcg density-dn/ds correlations are quite to the opposite between the two lineages.
such divergences appear to suggest that the biological roles of cpg methylation in coding exons have diverged significantly between the two lineages.
we also report here that exon length is an indicator of evolutionary rates of coding exons in arabidopsis species.
and this is not confounded by the ase/cse exon type, the proportion of repetitive elements, the proportion of intrinsically disordered regions, or exonic expression level.
one may suspect that this observation has resulted from alignment errors, leading to increased dn and dn/ds ratios in longer exons.
however, this is unlikely to be the case for two reasons.
firstly, the compared species - a. thaliana and a. lyrata - are very closely related.
the median dn value of first exons  is smaller than  <dig>  .
alignment errors may be a minor issue for sequence pairs with such a high level of similarity.
secondly, the length dependence of dn and dn/ds ratio is unlikely to result from the alignment between paralogous exonic sequences.
this is because to observe such length dependence, we should have systematically aligned orthologous sequences for shorter exons but paralogous sequences for longer exons.
we perceive no possible reasons why this may happen.
another possible explanation for the length dependence of dn and dn/ds ratio is annotation error.
however, this may not be a major problem judging from the small evolutionary rates as shown in figure  <dig>
the coding exons of animal and plant genes differ from each other in a number of biological features.
one example is microrna  targeting sites.
previous studies have reported that genes targeted by more mirnas tend to be under stronger selective constraints  <cit> .
a recent study indicated that in mammals, approximately 2% of the synonymous sites were selectively constrained for such regulatory sequences as splicing motifs, enhancers, and mirna target sites  <cit> .
for a. thaliana, it was predicted that ~75% of mirna target sites were located in cds  <cit> .
in comparison, only  <dig> % and  <dig> % of mirna targets were predicted to reside in cds in human and mouse, respectively  <cit> .
one important question is whether differential mirna targeting is the true reason for the differences in the mcg density-evolutionary rate correlations between arabidopsis species and mammals .
recall that the differences between the two lineages lie mainly in the mcg density-ds correlations in internal and last exons.
these correlations are significantly positive in mammals but negative in arabidopsis.
this divergence implies that for internal and last exons in mammals, the principal biological role of mcg is mutagenesis.
in arabidopsis, however, mcg density may be associated with other selection-constrained biological functions.
if the divergence in mcg density-ds correlations is to be ascribed to the higher proportion of mirna target sites in the cds of arabidopsis, three prerequisites should be fulfilled:  in the internal and last exons of arabidopsis, mcg density must be positively correlated with the probability of mirna targeting;  mirna targeting must be significantly constrained by selection in the two exon groups of arabidopsis; and  this mirna targeting-related selection affects only synonymous sites in internal and last exons of arabidopsis.
an example of mirna-mediated dna methylation has been reported for rice  <cit> .
the authors discovered that a specific group of 24-nucleotide  mirnas could mediate dna methylation within a ~80-nt region around the target sites.
however, only five such targets were identified.
and most of the methylation occurred in the chh or chg context  <cit> .
a follow-up study published lately showed that  <dig> of 24-nt mirnas exhibited elevated chh methylation  around their target sites  <cit> .
these studies imply that mirna targeting may lead to an increased level of dna methylation in the gene body of plants .
of note, nevertheless, each mirna was predicted to have only one target site in the target gene.
furthermore, only  <dig> of the  <dig> target sites were located in cds  <cit> .
meanwhile, a recent study suggested that the mirna target sites in cds were subject to negative selection  <cit> .
these observations seem to suggest a connection between mirna targeting and the mcg density-ds correlations in plants.
however, we speculate that the influences of mirna targeting might be insubstantial for three reasons.
first, only a relatively small number  of mirnas have been reported to cause dna methylation at the target sites.
and most of them occur outside of cds.
mirna-mediated methylation in cds thus may be uncommon.
second, the sequences that are subject to mirna-mediated methylation account for a minority  in light of the average cds length of ~1300 bp in the a. thaliana genome  <cit> .
certainly, we cannot exclude the possibility that a methylation-inducing mirna has multiple target sites in one gene, or that a gene is targeted by multiple methylation-inducing mirnas.
in such cases, the effects of mirna targeting will undoubtedly be non-negligible.
nevertheless, these scenarios were not observed in the recent studies  <cit> .
the overall influences of mirna targeting on cds methylation thus might be immaterial.
third, the identified mirna-mediated dna methylation occurred mostly in the chh or chg contexts  <cit> .
since we focus on methylation at cpg dinucleotides, the influences of mirna-mediated methylation on our analysis should be fairly limited.
another potential confounding factor in the mcg density-evolutionary rate analysis is the level of protein phosphorylation.
phosphorylated amino acid residues have been known to evolve more slowly than those unphosphorylated  <cit> .
since the motifs for phosphorylation differ between arabidopsis and mammals  <cit> , the evolutionary rates of coding exons in the two lineages may be differentially affected by phosphorylation-related constraints.
however, phosphorylation occurs at amino acid residues.
the selective constraints at the amino acid level influence dn but not ds.
note that the mcg density-dn correlations are generally similar between mammals and arabidopsis .
therefore, phosphorylation appears to have no significant effects on the differences in the mcg density-dn correlations between the two lineages.
one may suspect that the correlations between exon length and dn and dn/ds ratio have resulted from functional biases.
this is because exons of different lengths may belong to genes of different functional categories.
to examine this possibility, we divided the background dataset  into five length subgroups and conducted an all-to-all pairwise comparison of gene ontology functional categories between the five subgroups of internal exons using fatigo  <cit> .
as shown in additional file  <dig>  although the five length subgroups of internal exons differ from one another in view of gene ontology annotations, we do not observe any particular trend that may cause the length dependence of dn and dn/ds ratio.
we also examined whether the correlations between mcg density and evolutionary rate could differ between different functional categories.
we classified the analyzed genes according to the third level of “molecular function” of gene ontology, and calculated the correlations for nine functional groups that included ≥  <dig> genes.
note that one gene can be assigned to multiple functional groups.
the sum of genes in all of the functional groups thus outnumbers the analyzed genes.
the mcg density-evolutionary rate correlations in individual functional groups are similar to what we observed in figure  <dig> .
therefore, functional bias may not be a major concern in our analysis.
the correlations between exon length and evolutionary rates in arabidopsis species have been previously observed  <cit> .
however, the underlying mechanism remains unclear.
here we show that first, last, and internal coding exons diverge from each other in terms of the exon length-dn/ds ratio correlation – the correlation is stronger in internal exons than in last exons, and is statistically insignificant in first exons.
the length dependence of dn/ds ratio in last and internal exons remains statistically significant after controlling for potential confounding factors .
of note, for last and internal exons, this length dependence occurs because longer exons have a larger increase in dn than in ds when compared with shorter exons.
this increase in dn is probably unrelated to structural-functional reasons, for the proportion of disordered protein region  does not significantly affect the exon length-dn/ds ratio correlations.
it will be interesting to test the evolutionary neutrality of exons of different lengths when adequate polymorphism data become available.
meanwhile, it has been recently reported that in human, transcription factor binding sites  frequently reside in coding exons, and may significantly affect the evolution of these exonic sequences  <cit> .
the same comment may also apply to a. thaliana.
however, currently no base-resolution tfbs datasets are available for a. thaliana.
we may revisit this issue and investigate whether the density of tfbs is associated with the observed length dependence of dn and dn/ds ratio when such datasets are accessible.
one important issue is that we analyzed only one plant species in this study.
whether the observations in a. thaliana can be applied to other plant species remains unknown.
to address this issue, we retrieved three genome-scale methylome datasets of rice .
two of the datasets were derived from young panicles  <cit> , and the other was derived from leaves  <cit>  .
our analysis confirmed the mutagenic effect of mcg on coding exons and the stronger mutagenic effect on non-first exons than on first exons in rice .
the evolutionary rates of first, last, and internal exons were similar to what we observed in a. thaliana .
intriguingly, however, the correlations between mcg density and evolutionary rates were fairly different between rice  and a. thaliana .
particularly, in view of the mcg density-ds correlations in last and internal exons, rice was similar to mammals  <cit>  but not to a. thaliana.
of note, the rice methylome data were derived from panicles and leaves but not gamete cells.
whether the identified mcgs and the associated substitutions are heritable is therefore questionable.
to be sure, we cannot rule out the possibility that the differences in mcg density-ds correlations between a. thaliana and rice represent genuine divergences in the biological roles of mcg.
adding to the complexity of this issue is that the domesticated rice  has been artificially selected.
it will be interesting to re-examine this topic when the gamete methylome datasets of both cultivated and wild rice are available.
CONCLUSIONS
the mammal-arabidopsis divergence in the association between dna methylation and coding exon evolution is unexpected.
dna methylation is a major source of genomic sequence mutation on one hand, and an important transcriptional/splicing regulator on the other hand.
our results imply that this balance between biological roles of dna methylation in coding exons may have differed significant between arabidopsis and mammals in a length- and position-dependent manner.
the detailed evolutionary mechanisms and functional outcomes are worth further explorations.

 the evolution of the coding exome is a major driving force of functional divergence. in the past, a coding gene was considered as a basic unit for biological regulations and molecular functions. as such, in the majority of evolutionary studies, the “functional unit” targeted by natural selection is presumed to be a gene. however, with the advances in molecular biology and high-throughput sequencing technologies, it has gradually become clear that alternative transcript isoforms of the same gene can be spatio-temporally regulated, and convey fairly divergent biological functions. in other words, in many cases, a “transcript” rather than a “gene” is the biologically functional unit. the importance of transcript isoforms is particularly significant in complex organisms because they have highly developed networks of transcript/protein isoforms. transcript isoforms of the same gene differ from each other by alternatively spliced exonic regions. in cases where transcript isoforms convey distinct biological functions, the alternatively spliced exonic regions are crucial for the between-isoform functional divergences. these exonic regions should be accordingly targeted by natural selection. therefore, the biological functions of alternative exonic sequences and the selection pressure thereon can be revealed by examining the evolutionary patterns of these sequences. we previously examined the determinants of exonic evolutionary rates in mammals and arabidopsis species. the biological factors that affect exonic evolutionary rates were found to differ between these two lineages. in addition, we discovered that in mammals, the position of an exon in the transcript is significantly associated with the evolution of the exonic sequence in accordance with the level of dna methylation at cpg dinucleotides. this is probably because the position of an exon is related to its biological function, thus making the exon selectively constrained for the function mediated by cpg methylation, or prone to the mutagenesis effect of cpg methylation. however, whether this proposition is also true for arabidopsis remains unexplored. plant coding exons differ from their mammalian counterparts in several aspects. firstly, alternative rna splicing is less well developed, and plays a less important role in exon evolution in plants than in mammals. secondly, on average, a plant gene includes fewer but longer exons than a mammalian gene. thirdly, the effective population sizes of plants are considerably larger than those of mammals, giving rise to a higher efficiency of natural selection on plant exonic sequences. given these differences, we expect the evolutionary patterns of arabidopsis exons at different positions to diverge from those of their mammalian counterparts. in this study, we systematically examined the mutational effects and of cpg methylation and its correlations with exonic evolutionary rates for a. thaliana coding exons at different positions. our results indicate that first, last, and internal coding exons of a. thaliana have fairly different evolutionary patterns in this regard. the three exon groups diverge significantly in their liability to cpg methylation-related mutagenesis. furthermore, the cpg methylation-evolutionary rate correlations differ significantly among the three exon groups. these correlations also differ significantly between arabidopsis species and mammals. in addition, we found exon length to be a cpg methylation-independent indicator of exonic evolutionary rates in arabidopsis species. our results suggest that intrinsic exonic features may be important determinants for the evolution of a. thaliana coding exons, and that cpg methylation may play different biological roles in the coding exons of mammals and arabidopsis species.  the mutagenic effect of cpg methylation for exons at different positions to examine the mutagenic effect of cpg methylation on a. thaliana coding exons, we calculated the pearson’s coefficient of correlation between the level of cpg methylation and the cpg o/e ratio. of note, here we do not include methylation at chg or chh sites because these two types of methylation account for a minority of the plant methylome, and may have minor effects on the evolution of arabidopsis exons. cpg methylation can significantly increase the rate of cytosine-to-thymine transitions, leading to a decreased number of cpg dinucleotides. therefore, mcg density is expected to be negatively correlated with the cpg o/e ratio. furthermore, a larger absolute value of the coefficient of correlation indicates a stronger mutagenic effect of dna methylation. here the pearson’s correlation is employed to show this quantitative relationship. of note, when dealing with the methylome data, we applied a set of filtering criteria to ensure data quality and to reduce variations in the estimation of mcg density. these filtering criteria lead to differences in the number of analyzable exons among the four methylome datasets. interestingly, as shown in figure 1a, despite the variations in the number of analyzable exons and the level of cpg methylation, mcg density is consistently negatively correlated with the cpg o/e ratio across datasets. and the r value falls within a relatively narrow range. this observation confirms the mutagenic effect of cpg methylation on a. thaliana coding exons. note that the background dataset was not filtered for methylation data, but was filtered for the definitions of first, last, and internal exons and for the length threshold for calculations of evolutionary rates. na: not applicable. note that the “bisulfite seq. read depth” is defined as the total length of the bisulfite sequencing reads divided by the size of the a. thaliana genome. the next question to ask is whether the mutagenic effect of cpg methylation differs for first, last, and internal exons. to address this issue, we evaluated the pearson’s correlations as described above separately for each of the three exon groups. unexpectedly, as shown in figure 1b, although the mcg-cpg o/e correlations remain negative across the three exon groups, the strongest and the weakest mutagenic effect occur, respectively, in internal and first exons. this is to the contrary of what was previously observed for mammals, where the strongest mutagenic effect of cpg methylation occurs in first exons, while the weakest in internal exons. we then ask how the variation in mcg-related mutagenic effect may have affected the evolutionary rates of the three exon groups. to this end, we first examined the evolutionary rates separately for the three exon groups between a. thaliana and a. lyrata based on datasets s1 ~ s. the four sperm methylome datasets used here are appropriate for this evolutionary analysis because only the mutations that occur in germ line cells can be propagated over generations, thus leaving observable changes in the inter-species comparison. figure  shows that among the three exon groups, first exons have the largest median dn and dn/ds ratio, followed by last exons, and finally by internal exons. for ds, internal exons have the lowest median value, and first and last exons have similar values. the evolutionary rate profiles apparently are inconsistent with the mcg-related mutagenic effect profiles in first, last, and internal exons. specifically, in internal exons, we observe the co-occurrence of a high median mcg density, a strong mutagenic effect of cpg methylation, and low evolutionary rates as compared with the other two exon groups. one possible explanation is that the strong selection pressure imposed on internal exons has significantly constrained the mcg-related mutations from occurring in this exon group. this appears to be true judging from the low dn/ds ratio in internal exons as compared with the other two exon groups. interestingly, the median ds is also the lowest in internal exons, suggesting that synonymous substitutions are subject to strong purifying selection in this exon group. the correlations between cpg methylation and exonic evolutionary rates to investigate in more details the correlations between cpg methylation and exonic evolutionary rates, we evaluated spearman’s correlations between mcg density and each of the three evolutionary measurements separately for the four sperm methylome datasets. figure  shows that for first exons, mcg density has weak positive correlations with dn, ds and the dn/ds ratio. this observation suggests that the mutagenic effect of cpg methylation influences both synonymous and nonsynonymous sites in first exons, but the effects may be relatively small. for last exons, the spearman’s coefficients of correlation are consistently negative between mcg density and each of dn, ds, and the dn/ds ratio across datasets. this observation is somewhat surprising because a negative mcg density-dn correlation indicates that a higher level of cpg methylation is accompanied with a reduced rate of sequence evolution. in mammals, by comparison, the mcg density-ds correlation is positive, while the mcg density-dn and mcg density-dn/ds correlations are negative in last exons. meanwhile, for internal exons, the mcg density-ds correlations are consistently negative across datasets. however, the mcg density-dn correlations are weakly negative in this exon group. unexpectedly, the mcg density-dn/ds ratio correlations are positive in internal exons. these observations appear to suggest that the synonymous sites in internal exons are subject to certain mcg-associated selective constraints, thus heavily methylated exons tend to have lower ds values. however, such selective constraints may be weaker for the nonsynonymous sites in the same exonic regions. the decrease in ds and the relative stasis in dn lead to the increase in the dn/ds ratio as mcg density increases in internal exons. these results also differ from what were reported for mammals, where the mcg density-ds correlation is positive, and the mcg- dn and mcg density-dn/ds ratio correlations are both negative in internal exons. the differences in the mcg density-evolutionary rate correlations between a. thaliana and mammals suggest that the biological roles of cpg methylation may have diverged between the two lineages. we also conducted the same analyses while controlling for four potential confounding factors or alternative exon type). the results remain virtually the same. the associations between exon length and exonic sequence evolution we have shown that the mcg density-cpg o/e correlation is the most marked in internal exons. one important question is whether this observation has actually resulted from certain sampling biases. we noticed that our filtering criteria for dealing with the methylome data tend to retain longer exons in the datasets. the average lengths of the selected exons in s1 ~ s range from to base pairs, which are considerably longer than the previously reported exome-wide average of ~250 bp in a. thaliana and the background value in this study. to investigate whether exon length has affected our results, we divided each of the methylome datasets further into five subgroups of approximately equal sizes according to exon length, and evaluated the mcg density-cpg o/e correlations separately for first, last, and internal exons for each length subgroup. surprisingly, as shown in figure  first exons show a clear length-dependent decrease in the r value between mcg density and the cpg o/e ratio, indicating stronger mutagenic effects of cpg methylation on longer first exons. by comparison, internal and last exons do not show similarly clear trends. since exon length is related to the mutagenic effect of cpg methylation, we then ask whether exon length is correlated with the evolutionary rates of first, last, and internal exons. accordingly, we evaluated spearman’s correlations between exon length and dn, ds, and the dn/ds ratio separately for the three exon groups. of note, here we take the background dataset instead of the methylome datasets to avoid potential sampling biases. interestingly, table  shows that for first exons, exon length is significantly correlated with both dn and ds, although the ρ value is smaller in the exon length-dn correlation. meanwhile, the exon length-dn/ds ratio correlation is statistically insignificant for this exon group. this observation seems to suggest that the increased mutagenic effect of cpg methylation in longer exons has led to increases in both dn and ds in this exon group. this conjecture, however, is only partly true, as will be shown later. meanwhile, for last and internal exons, the correlations between exon length and the three evolutionary measurements are all statistically significant. yet a noteworthy difference between last and internal exons is that the ρ value of the exon length-dn correlation is about four times larger than that of the exon length-ds correlation for internal exons. this difference in ρ value is smaller in last exons. d n / d s ratio, d n, and d s before and after controlling for four potential confounding factors note that this table is based on the background dataset. the above observations may be confounded by other biological factors. for example, ases are known to have increased dn and the dn/ds ratios as compared with cses. therefore, the increase in dn and the dn/ds ratio in longer exons might have resulted from an increase in the proportion of ases. meanwhile, the proportion of repetitive elements is also correlated with evolutionary rates because these elements are subject to relaxed selective constraints. similar comments also apply to intrinsically disordered protein regions. the next factor to consider is exonic expression level, which has been shown to be an important determinant of dn and the dn/ds ratio. we thus conducted partial spearman’s correlation analyses while simultaneously controlling all of these four factors. as shown in table , the results remain virtually the same. the last but a critical factor to control is the level of cpg methylation. to evaluate the influence of cpg methylation, we have to employ the sperm methylome datasets, which include considerably fewer but longer exons as compared with the background dataset. using the four sperm methylome datasets, we again conducted partial spearman’s correlation analysis between exon length and evolutionary rates while simultaneously controlling for mcg density, ase/cse exon type, proportion of repetitive elements/disordered region, and exonic expression level. as shown in additional file for first exons, the results are similar but the ρ values are decreased. this observation indicates that the length dependence of mcg-related mutagenic effect in first exons accounts for part but not all of the length dependence of dn and ds in this exon group. meanwhile, for last exons, the exon length-ds correlation becomes statistically insignificant in s and s whereas the exon length-dn and the exon length-dn/ds correlations remain statistically significant with decreased ρ values. these results imply that mcg-related mutations may account for part of the length dependence of ds, dn, and the dn/ds ratio in last exons. however, the decreases in ρ value and the level of statistical significance may also be ascribable to the decrease in sample size and the bias in exon length. by comparison, for internal exons, all of the correlations remain statistically significant with two notable changes as compared with the results in table 2: the ρ values of the exon length-dn and the exon length-dn/ds correlations are increased; and the ρ values of the exon length-ds correlations turn negative. therefore, for internal exons, mcg-related mutations appear to be an important factor affecting ds. nevertheless, mcg-related mutations cannot explain the length dependence of dn and the dn/ds ratio in this exon groups. taken together, our results indicate that the correlations between exon length and the evolutionary measurements are unaffected by the ase/cse exon type, proportion of repetitive elements/disordered region, and exonic expression level in any of the three exon groups. however, the level of cpg methylation may account for part of the exon length-evolutionary rate correlations differentially for first, last, and internal exons. in summary, exon length appears to be a cpg methylation-independent indicator for dn in all of the three exon groups, and for the dn/ds ratio in last and internal exons of a. thaliana. discussion we have shown that for the coding sequences of a. thaliana, the mutagenic effects of cpg methylation differ between exons at different relative positions. among the three compared exon groups, the highest cpg methylation level and the strongest mutagenic effect of cpg methylation both occur in internal coding exons despite the most stringent selective constraint on this exon group. first coding exons, quite to the opposite, have the lowest level of cpg methylation and suffer the weakest mutagenic effect of cpg methylation, yet evolve the most rapidly. interestingly, we show that mcg density is positively correlated with ds, dn, and dn/ds ratio in first exons, yet the same correlations are significantly negative for last exons. for internal exons, the correlations are negative, weakly negative, and positive for ds, dn, and dn/ds ratio, respectively. the mutagenic effect of cpg methylation cannot fully explain these observations. apparently, selection pressure has played a major role here. we have previously reported that in mammals, cpg methylation may have different biological roles in first, last, and internal coding exons. similar comments may also apply to arabidopsis species – that first exons are more liable to the mutagenic effects, yet the other two exon groups are more affected by the regulatory functions of cpg methylation. noticeably, however, the correlations between mcg density and evolutionary rates actually diverge significantly between arabidopsis species and mammals. one riveting difference is that for internal exons, the mcg density-ds and mcg density-dn/ds correlations are quite to the opposite between the two lineages. such divergences appear to suggest that the biological roles of cpg methylation in coding exons have diverged significantly between the two lineages. we also report here that exon length is an indicator of evolutionary rates of coding exons in arabidopsis species. and this is not confounded by the ase/cse exon type, the proportion of repetitive elements, the proportion of intrinsically disordered regions, or exonic expression level. one may suspect that this observation has resulted from alignment errors, leading to increased dn and dn/ds ratios in longer exons. however, this is unlikely to be the case for two reasons. firstly, the compared species - a. thaliana and a. lyrata - are very closely related. the median dn value of first exons is smaller than. alignment errors may be a minor issue for sequence pairs with such a high level of similarity. secondly, the length dependence of dn and dn/ds ratio is unlikely to result from the alignment between paralogous exonic sequences. this is because to observe such length dependence, we should have systematically aligned orthologous sequences for shorter exons but paralogous sequences for longer exons. we perceive no possible reasons why this may happen. another possible explanation for the length dependence of dn and dn/ds ratio is annotation error. however, this may not be a major problem judging from the small evolutionary rates as shown in figure  the coding exons of animal and plant genes differ from each other in a number of biological features. one example is microrna targeting sites. previous studies have reported that genes targeted by more mirnas tend to be under stronger selective constraints. a recent study indicated that in mammals, approximately 2% of the synonymous sites were selectively constrained for such regulatory sequences as splicing motifs, enhancers, and mirna target sites. for a. thaliana, it was predicted that ~75% of mirna target sites were located in cds. in comparison, only % and % of mirna targets were predicted to reside in cds in human and mouse, respectively. one important question is whether differential mirna targeting is the true reason for the differences in the mcg density-evolutionary rate correlations between arabidopsis species and mammals. recall that the differences between the two lineages lie mainly in the mcg density-ds correlations in internal and last exons. these correlations are significantly positive in mammals but negative in arabidopsis. this divergence implies that for internal and last exons in mammals, the principal biological role of mcg is mutagenesis. in arabidopsis, however, mcg density may be associated with other selection-constrained biological functions. if the divergence in mcg density-ds correlations is to be ascribed to the higher proportion of mirna target sites in the cds of arabidopsis, three prerequisites should be fulfilled: in the internal and last exons of arabidopsis, mcg density must be positively correlated with the probability of mirna targeting; mirna targeting must be significantly constrained by selection in the two exon groups of arabidopsis; and this mirna targeting-related selection affects only synonymous sites in internal and last exons of arabidopsis. an example of mirna-mediated dna methylation has been reported for rice. the authors discovered that a specific group of 24-nucleotide mirnas could mediate dna methylation within a ~80-nt region around the target sites. however, only five such targets were identified. and most of the methylation occurred in the chh or chg context. a follow-up study published lately showed that of 24-nt mirnas exhibited elevated chh methylation around their target sites. these studies imply that mirna targeting may lead to an increased level of dna methylation in the gene body of plants. of note, nevertheless, each mirna was predicted to have only one target site in the target gene. furthermore, only of the target sites were located in cds. meanwhile, a recent study suggested that the mirna target sites in cds were subject to negative selection. these observations seem to suggest a connection between mirna targeting and the mcg density-ds correlations in plants. however, we speculate that the influences of mirna targeting might be insubstantial for three reasons. first, only a relatively small number of mirnas have been reported to cause dna methylation at the target sites. and most of them occur outside of cds. mirna-mediated methylation in cds thus may be uncommon. second, the sequences that are subject to mirna-mediated methylation account for a minority in light of the average cds length of ~1300 bp in the a. thaliana genome. certainly, we cannot exclude the possibility that a methylation-inducing mirna has multiple target sites in one gene, or that a gene is targeted by multiple methylation-inducing mirnas. in such cases, the effects of mirna targeting will undoubtedly be non-negligible. nevertheless, these scenarios were not observed in the recent studies. the overall influences of mirna targeting on cds methylation thus might be immaterial. third, the identified mirna-mediated dna methylation occurred mostly in the chh or chg contexts. since we focus on methylation at cpg dinucleotides, the influences of mirna-mediated methylation on our analysis should be fairly limited. another potential confounding factor in the mcg density-evolutionary rate analysis is the level of protein phosphorylation. phosphorylated amino acid residues have been known to evolve more slowly than those unphosphorylated. since the motifs for phosphorylation differ between arabidopsis and mammals, the evolutionary rates of coding exons in the two lineages may be differentially affected by phosphorylation-related constraints. however, phosphorylation occurs at amino acid residues. the selective constraints at the amino acid level influence dn but not ds. note that the mcg density-dn correlations are generally similar between mammals and arabidopsis. therefore, phosphorylation appears to have no significant effects on the differences in the mcg density-dn correlations between the two lineages. one may suspect that the correlations between exon length and dn and dn/ds ratio have resulted from functional biases. this is because exons of different lengths may belong to genes of different functional categories. to examine this possibility, we divided the background dataset into five length subgroups and conducted an all-to-all pairwise comparison of gene ontology functional categories between the five subgroups of internal exons using fatigo. as shown in additional file although the five length subgroups of internal exons differ from one another in view of gene ontology annotations, we do not observe any particular trend that may cause the length dependence of dn and dn/ds ratio. we also examined whether the correlations between mcg density and evolutionary rate could differ between different functional categories. we classified the analyzed genes according to the third level of “molecular function” of gene ontology, and calculated the correlations for nine functional groups that included ≥  genes. note that one gene can be assigned to multiple functional groups. the sum of genes in all of the functional groups thus outnumbers the analyzed genes. the mcg density-evolutionary rate correlations in individual functional groups are similar to what we observed in figure . therefore, functional bias may not be a major concern in our analysis. the correlations between exon length and evolutionary rates in arabidopsis species have been previously observed. however, the underlying mechanism remains unclear. here we show that first, last, and internal coding exons diverge from each other in terms of the exon length-dn/ds ratio correlation – the correlation is stronger in internal exons than in last exons, and is statistically insignificant in first exons. the length dependence of dn/ds ratio in last and internal exons remains statistically significant after controlling for potential confounding factors. of note, for last and internal exons, this length dependence occurs because longer exons have a larger increase in dn than in ds when compared with shorter exons. this increase in dn is probably unrelated to structural-functional reasons, for the proportion of disordered protein region does not significantly affect the exon length-dn/ds ratio correlations. it will be interesting to test the evolutionary neutrality of exons of different lengths when adequate polymorphism data become available. meanwhile, it has been recently reported that in human, transcription factor binding sites frequently reside in coding exons, and may significantly affect the evolution of these exonic sequences. the same comment may also apply to a. thaliana. however, currently no base-resolution tfbs datasets are available for a. thaliana. we may revisit this issue and investigate whether the density of tfbs is associated with the observed length dependence of dn and dn/ds ratio when such datasets are accessible. one important issue is that we analyzed only one plant species in this study. whether the observations in a. thaliana can be applied to other plant species remains unknown. to address this issue, we retrieved three genome-scale methylome datasets of rice. two of the datasets were derived from young panicles, and the other was derived from leaves. our analysis confirmed the mutagenic effect of mcg on coding exons and the stronger mutagenic effect on non-first exons than on first exons in rice. the evolutionary rates of first, last, and internal exons were similar to what we observed in a. thaliana. intriguingly, however, the correlations between mcg density and evolutionary rates were fairly different between rice and a. thaliana. particularly, in view of the mcg density-ds correlations in last and internal exons, rice was similar to mammals but not to a. thaliana. of note, the rice methylome data were derived from panicles and leaves but not gamete cells. whether the identified mcgs and the associated substitutions are heritable is therefore questionable. to be sure, we cannot rule out the possibility that the differences in mcg density-ds correlations between a. thaliana and rice represent genuine divergences in the biological roles of mcg. adding to the complexity of this issue is that the domesticated rice has been artificially selected. it will be interesting to re-examine this topic when the gamete methylome datasets of both cultivated and wild rice are available.  the mammal-arabidopsis divergence in the association between dna methylation and coding exon evolution is unexpected. dna methylation is a major source of genomic sequence mutation on one hand, and an important transcriptional/splicing regulator on the other hand. our results imply that this balance between biological roles of dna methylation in coding exons may have differed significant between arabidopsis and mammals in a length- and position-dependent manner. the detailed evolutionary mechanisms and functional outcomes are worth further explorations. 
exons at different positions in the transcript or in different transcript isoforms may mutate at different rates due to variations in dna methylation level; and serve distinct biological roles, and thus be differentially targeted by natural selection. our results suggest that intrinsic features, including relative exonic position in the transcript and exon length, play an important role in the evolution of a. thaliana coding exons. fourthly, exon length is a cpg methylation-independent indicator of evolutionary rates, particularly for dn and the dn/ds ratio in last and internal exons. finally, the evolutionary patterns of coding exons with regard to cpg methylation differ significantly between arabidopsis species and mammals. firstly, the mutagenic effect of cpg methylation is the strongest for internal exons and the weakest for first exons despite the stringent selective constraints on the former group. secondly, the mutagenic effect of cpg methylation increases significantly with length in first exons but not in the other two exon groups. the evolution of the coding exome is a major driving force of functional divergence both between species and between protein isoforms. thirdly, cpg methylation level is correlated with evolutionary rates with markedly different patterns among the three exon groups. furthermore, cpg methylation is correlated with exonic evolutionary rates differentially between a. thaliana and animals, and may have served different biological roles in the two lineages. furthermore, intrinsic exonic features, such as exon length, may also affect the evolution of individual exons. dna methylationexon evolutionevolutionary raterelative exon positionexon length the correlations are generally positive, negative, and mixed for first, last, and internal exons, respectively. we observed fairly different evolutionary patterns in a. thaliana as compared with those reported for animals. here we examine how dna methylation at cpg dinucleotides, in the context of intrinsic exonic features, influences the evolution of coding exons of arabidopsis thaliana.
 the evolution of the coding exome is a major driving force of functional divergence both between species and between protein isoforms. exons at different positions in the transcript or in different transcript isoforms may mutate at different rates due to variations in dna methylation level; and serve distinct biological roles, and thus be differentially targeted by natural selection. furthermore, intrinsic exonic features, such as exon length, may also affect the evolution of individual exons. importantly, the evolutionary effects of these intrinsic/extrinsic features may differ significantly between animals and plants. such inter-lineage differences, however, have not been systematically examined.  here we examine how dna methylation at cpg dinucleotides, in the context of intrinsic exonic features, influences the evolution of coding exons of arabidopsis thaliana. we observed fairly different evolutionary patterns in a. thaliana as compared with those reported for animals. firstly, the mutagenic effect of cpg methylation is the strongest for internal exons and the weakest for first exons despite the stringent selective constraints on the former group. secondly, the mutagenic effect of cpg methylation increases significantly with length in first exons but not in the other two exon groups. thirdly, cpg methylation level is correlated with evolutionary rates with markedly different patterns among the three exon groups. the correlations are generally positive, negative, and mixed for first, last, and internal exons, respectively. fourthly, exon length is a cpg methylation-independent indicator of evolutionary rates, particularly for dn and the dn/ds ratio in last and internal exons. finally, the evolutionary patterns of coding exons with regard to cpg methylation differ significantly between arabidopsis species and mammals.  our results suggest that intrinsic features, including relative exonic position in the transcript and exon length, play an important role in the evolution of a. thaliana coding exons. furthermore, cpg methylation is correlated with exonic evolutionary rates differentially between a. thaliana and animals, and may have served different biological roles in the two lineages. dna methylationexon evolutionevolutionary raterelative exon positionexon length 