Linkage disequilibrium analysis of allelic heterogeneity in DNA methylation

Abstract

Heterogeneity of DNA methylation status among alleles is observed in various cell types and is involved in epigenetic gene regulation and cancer biology. However, the individual methylation profile within each allele has not yet been examined at the whole-genome level. In the present study, we applied linkage disequilibrium analysis to the DNA methylation data obtained from whole-genome bisulfite sequencing studies in mouse germline and other types of cells. We found that the methylation status of 2 consecutive CpG sites showed deviation from equilibrium frequency toward concordant linkage (both methylated or both unmethylated) in germline cells. In the imprinting loci where methylation of constituent alleles is known, our analysis detected the deviation toward the concordant linkage as expected. In addition, we applied this analysis to the transitional zone between methylated and unmethylated regions and to the cells undergoing epigenetic reprogramming. In both cases, deviation to the concordant-linked alleles was conspicuous, indicating that the methylation pattern is not random but rather concordant within each allele. These results will provide the key to understanding the mechanism underlying allelic heterogeneity.

Keywords:

• adjacent CpG sites
• allelic heterogeneity
• DNA methylation
• DNA methyltransferase
• germ cells

Introduction

DNA methylation is the most important epigenetic modification and is involved in cellular differentiation, fate determination, developmental plasticity, and environmental responses. The genome-wide pattern of methylation status is maintained through cell division and is thus inherited by the descendant cell population. While the methylation status of individual CpG sites is generally consistent among alleles, resulting in methylated or unmethylated CpG sites in all alleles, there are certain regions that show heterogeneity in methylation status among alleles. The methylation heterogeneity (intermediate methylation¹ or partially methylated domain²) is reported in primordial germ cells under epigenetic reprogramming,³ sperms,⁴ embryonic stem cells,⁵ induced pluripotent stem cells,⁶ placental cells,² and other somatic cells under normal and disease conditions.^7–9 These regions are highly gene enriched and often interact with other epigenetic markers such as histone modifications.¹ In embryonic stem cells, allelic heterogeneity is involved in maintaining the bimodal states of gene expression.⁵ Furthermore, allelic heterogeneity is frequently found in various types of cancer cells, and the possible relationship with aberrant gene expression in cancer cells is of particular interest.⁹ While previous studies have defined allelic heterogeneity based mainly on the methylation levels of individual CpG sites, the underlying methylation pattern within each allele has not yet been examined. Accordingly, it is important to determine the linear methylation pattern of consecutive CpG sites within each allele.

In the present study, we introduce a novel approach for the analysis of the methylation pattern within each allele using bisulfite sequencing reads. Characterization of the methylation pattern underlying allelic heterogeneity will provide the basis for understanding the dynamics of heterogeneity, its biological significance, and its roles in epigenetic gene regulation.

Results and Discussion

Coefficient of linkage disequilibrium D

We used the published data of whole-genome bisulfite sequencing experiments on mouse embryonic germ cells, developing and mature oocytes, sperms, ES cells, and astrocytes.^3,10,11 To characterize the methylation pattern within each allele underlying heterogeneity, we focused on the methylation status of the 2 adjacent CpG sites (CpG pair) residing in a single sequencing read, which is derived from a single allele. We examined the linkage of methylation states in a CpG pair using a modified version of the coefficient of linkage disequilibrium, hereafter designated as D. The coefficient D corresponds to the deviation of the observed frequency of alleles from their equilibrium frequency. The range of D depends on the allele frequency, i.e., D has a maximum range when the allele frequency of each locus is 50%, whereas it becomes 0 when the frequency of one allele is 100%. With respect to the linkage of methylation status, D has maximum range when the methylation level is intermediate, whereas it becomes 0 when all alleles are fully methylated or unmethylated (Fig. 1A). Because the region of allelic heterogeneity shows the intermediate methylation level, D is suitable for our analysis. Furthermore, the sign of D can be used to determine the type of linkage, i.e., both methylated or both unmethylated (concordant) or alternately methylated and unmethylated (discordant).

Figure 1. Coefficient of linkage disequilibrium D indicates the linkage between 2 adjacent CpG sites. (A) The methylation level and the range of D. The D has a maximum range at the intermediate methylation level. The sign of D corresponds to the concordant (red bars) or discordant (blue bars) linkage. (B) An example plot of methylation level and D in mouse embryonic stem cells at the genomic region chr1:134,200,000–134,600,000. Upper track: methylation levels of individual CpG sites (black). Lower track: D values of individual CpG pairs (red for positive and blue for negative values). Bottom track shows the gene annotation.

For the mouse reference genome (mm10), the total number of CpG pairs with the distance of <100 bp is 14,137,359, of which the median distance is 28 bp. We first examined the relationship between D and the methylation level of the CpG pairs (Fig. 1B). In the examined germline cells, D generally showed positive values at intermediate methylation levels (0.25–0.75) and approached zero at higher or lower methylation levels (Fig. 2A). The positive D values at the intermediate methylation levels indicate that the methylation status of CpG pairs at the regions of heterogeneity tend to show concordant linkage. This tendency was conspicuous in developing germline cells but not as evident in embryonic stem cells and adult astrocyte cells (Fig. 2B), whereas the means were significantly different between all samples examined (${\displaystyle P\le {10}^{-10}}$. t-test). The embryonic germ cells undergo epigenetic reprogramming during which methylation levels change dramatically.³ Thus, it is possible that positive D values observed in germline cells may be the result of global changes in methylation status. In contrast, methylation patterns are relatively stable in adult somatic cells. If the stochastic changes continuously accumulate on stable methylation patterns in somatic cells, it might be expected that D would approach equilibrium. Thus, the value of D close to zero observed in adult somatic cells indicates that epigenetic information carried by DNA methylation is continuously degraded along with cell division and differentiation.¹² Importantly, heterogeneity among alleles does not indicate stochastic methylation within each allele, i.e., the methylation pattern within each allele is not necessarily an interspersed mixture of methylated and unmethylated sites. In fact, the CpG pairs showing negative D values are rare and distributed sporadically, and there was no continuous region with highly discordant methylation pattern. The fact that positive D prevailed in the regions showing allelic heterogeneity indicates that the methylation pattern within each allele is generally concordant rather than stochastic.

Figure 2. Germline cells showed positive D values. (A) Boxplots of D at various methylation levels in male primordial germ cells at embryonic day 10.5. (B) Boxplot for D at intermediate methylation level (0.25–0.75) in various cell types (female and male primordial germ cells at embryonic days 10.5, 13.5, and 16.5, non-growing oocytes, germinal vesicle oocytes, sperm cells, embryonic stem cells, and astrocytes).

Heterogeneity in imprinting loci

Since the sign of D corresponds to concordant or discordant linkage, we consider that D might be suitable for analyzing the linear pattern of methylation within each allele. To test this, we focused on the imprinting loci where maternal and paternal alleles show a distinct methylation status. In the region involved in the regulation of paternal expression, the methylation level is low in sperm-derived alleles and high in oocyte-derived alleles.^13–15 Accordingly, the sexually differential methylation pattern in gametes results in allelic heterogeneity in the zygote, which might be inherited by the descendant cells (Fig. 3A). As expected, D was close to zero in both sperms and oocytes, reflecting their uniform methylation status in imprinting loci (Fig. 3B). In contrast, the deviation to concordant linkage was conspicuous in descendant cells, such as ES cells and astrocytes, as well as in 10.5 embryonic day germ cells where methylation erasure had not yet completed (Fig. 3B). The same tendency was additionally observed in several known imprinting loci (see Fig. S1). The result of D reflecting the linkage of methylation status indicates that D is suitable to infer the spatial pattern of methylation within alleles constituting heterogeneity.

Figure 3. D in known imprinting loci. (A) Schematic representation of allele constitutions at paternally imprinted loci for sperms, oocytes, and zygotes. Alleles in sperms and in oocytes are fully unmethylated and methylated, respectively. Alleles in zygotes are a mixture of unmethylated and methylated. (B) Plots of methylation levels (black) and D (red for positive and blue for negative values) at paternally imprinted loci (Peg10 and Kcnq1ot1) in sperms, oocytes, male and female primordial germ cells (embryonic day 10.5), embryonic stem cells, and astrocytes. While D is close to zero in sperms and oocytes, D becomes positive in primordial germ cells, embryonic stem cells, and astrocytes at imprinting loci. Located at the bottom is the track for gene annotation.

Heterogeneity in transitional zones

We then applied linkage disequilibrium analysis to various genomic regions showing heterogeneity. We first analyzed the transitional zone between fully methylated and fully unmethylated regions. In mouse oocytes, the transitional zone found at the boundary of a gene body often shows allelic heterogeneity.¹⁰ There are 2 possible explanations for the methylation pattern underlying heterogeneity in the transitional zone: 1) methylation activity is stochastic in the transitional zone, resulting in each allele consisting of an interspersed mixture of methylated and unmethylated CpG sites, or 2) the transition between methylated and unmethylated regions occurs at a distinct boundary in each allele, whereas the boundary position differs among alleles (Fig. 4A). In the former case, the methylation pattern in the transitional zone is expected to show no deviation (${\displaystyle D\cong 0}$), whereas in the latter case, there is a deviation to concordant linkage (${\displaystyle D>0}$). Accordingly, we examined the deviation at the transitional zones found in the methylome of mouse germinal vesicle oocytes.^10,11 In almost all regions examined, the CpG pair in the transitional zone showed a clear deviation to concordant linkage (Fig. 4B). These results indicate that methylation activity is not stochastic in the transitional zone, but the transitional boundaries differ among alleles. It is known that DNA methylation is found not only in promoter regions but additionally in the gene bodies. In addition, the boundary between methylated and unmethylated regions often corresponds to the boundary between the gene body and intergenic regions.¹⁶ Furthermore, gene body methylation is associated with intragenic chromatin modifications such as H3K9me3.¹⁷ Thus, the present finding that the transitional boundary differs among alleles may indicate that the boundary of gene body methylation differs among cells, possibly depending on the degree of histone modifications and/or activity of gene expression.

Figure 4. D in the transitional zone. (A) Two models explaining allelic heterogeneity in the transitional zone, i.e., boundary position is stochastic (upper panel) or methylation activity is stochastic (lower panel) at the transitional zone. (B) Plots showing the methylation level (upper panel, black) and D (lower panel, red for positive and blue for negative values) at the transitional zone (chr1:156,200,000–156,300,000), where methylation levels showed fluctuation around 0.5 and D showed conspicuous deviation to positive values. Bottom track shows the gene annotation.

Heterogeneity during reprogramming

We next examined the heterogeneity observed in the embryonic germ cells during the epigenetic reprogramming. Reprogramming of the DNA methylome consists of 2 phases, namely whole-genome demethylation and de novo methylation.^18–20 The allelic heterogeneity becomes most obvious at the period between these 2 phases.³ Knowledge of the methylation pattern underlying allelic heterogeneity may assist in understanding the mechanism of de novo methylation acting in the epigenetic reprogramming. We considered 2 hypotheses: 1) de novo methylation progresses stochastically at individual CpG sites, resulting in the spatially interspersed methylation patterns in all alleles, or 2) de novo methylation occurs stochastically among alleles; however, once present, it propagates rapidly, resulting in the mixture of fully methylated and unmethylated alleles (Fig. 5A). It is expected that deviation should not be observed in the former case (${\displaystyle D\cong 0}$), whereas deviation to the concordant pattern may be observed in the latter case (${\displaystyle D>0}$). We examined the deviation in male primordial germ cells (embryonic day 13.5 and 16.5) and in sperm cells, each corresponding to the progressive stages of reprogramming.³ On day 13.5, no deviation in male primordial germ cells is observed when methylation is fully erased and reprogramming has not yet started (Fig. 5B). In sperm cells where de novo methylation is completed and most of the CpG sites are fully methylated, only a small part showed heterogeneity with deviation to the concordant linkage (Fig. 5B). When reprogramming was in progress and almost all regions showed allelic heterogeneity, clear deviation to the concordant linkage was observed (Fig. 5B). These results support the latter hypothesis of allelic heterogeneity being caused by mixture of methylated and unmethylated alleles. It has been reported that in 16.5-day male primordial germ cells, DNA methyltransferases, Dnmt3a, Dnmt3b, and Dnmt3l show a burst of expression.²¹ Thus, it is possible that allelic heterogeneity may reflect the heterogeneity in timing of their expression burst among the cell population. It is known that mouse primordial germ cells show heterogeneity in the status of differentiation and proliferation,²² further supporting our hypothesis.

Figure 5. Changes in D during epigenetic reprogramming. (A) Schematic representation of changes in allelic methylation patterns during epigenetic reprogramming from erased (upper panel) to reprogrammed (lower panel). Middle panels show 2 models explaining allelic heterogeneity at the time of reprogramming, i.e., the mixture of fully methylated and unmethylated alleles or the spatially interspersed methylation patterns in all alleles. (B) Plots showing methylation level and D at the genomic region chr1:165,490,000 –166,000,000 in male primordial germ cells at embryonic days 13.5 and 16.5 and sperms. In 13.5 day germ cells, methylation is erased, while in sperms, reprogramming is completed. In 16.5 day germ cells, when methylation levels showed intermediate values, D tended to show positive values, thereby supporting the latter hypothesis. Bottom track shows the gene annotation.

Length of CpG pair and deviation D

The deviation to concordant linkage observed in various biological conditions suggests that DNA methyltransferases do not act in a site-specific manner but rather act by simultaneously targeting a range of CpG sites. It is thus expected that if the distance between neighboring CpG sites increases, the degree of deviation to concordant linkage might decrease. To test this, we examined the relationship between the length of CpG pair and D (Fig. 6). As expected, a significant negative correlation was identified in all samples examined (${\displaystyle P\le {10}^{-10}}$. t-test for Pearson correlation coefficient), indicating that closer neighboring CpG sites are more likely to be simultaneously methylated. It has been suggested that DNA methyltransferases are recruited to the target position by interaction with various chromatin modifications and remodeling factors, resulting in the target region of methylation being of at least a nucleosome size.²³ Our results support the nucleosomal but not the site-specific action of DNA methyltransferases. A recent model further suggests that DNA methylation is not a static process but rather reflects the dynamics between the methylation and de-methylation processes.²³ Heterogeneity among regionally methylated and unmethylated alleles may be the consequence of the dynamic methylation and demethylation processes. Considered collectively, our approach to analyze the linkage between the methylation status of adjacent CpG sites has proven to be a powerful tool to understand the mechanism of the methylation processes and is possibly useful for application to various cancer and disease epigenomics studies.

Figure 6. D depends on the length of the CpG pair. Box plots show the relationship between the length of the CpG pair and D at an intermediate methylation level (0.25–0.75) in male primordial germ cells (embryonic day 16.5).

Methods

Read data

We used the read data from a whole-genome bisulfite sequencing experiment in mouse primordial germ cells, oocytes, sperm cells, embryonic stem cells, and astrocytes.^3,10,11 The sequence reads were downloaded from the NCBI short read archives (accessions: DRA000484, DRA000495, DRA000570, and DRA000607).

Mapping of bisulfite reads

The reads were mapped onto the mouse reference genome (mm10) and analyzed using the Bismark software.²⁴ CpG pairs with the distance being less than 100 bp on the mm10 reference were listed and used for subsequent analyses. For each sample, the methylation status of the CpG pair on a single sequencing read was extracted and summarized. We included only the CpG pairs with depths of 10–200.

Calculation of deviation D

For each CpG pair, both the methylation level and the coefficient of linkage disequilibrium D were calculated from the frequency of the 4 allele types using the equation, as described below (note ${\displaystyle f_{MM}\text{+}{f}_{UU}+f_{MU}\text{+}{f}_{UM}=1}$). In the present study, we defined D as the deviation for concordant alleles (${\displaystyle D_{con}}$).

${\displaystyle f_{MM}}$: Frequency of methylated concordant alleles

${\displaystyle f_{UU}}$: Frequency of unmethylated concordant alleles

${\displaystyle f_{MU}}$, ${\displaystyle f_{UM}}$: Frequency of discordant alleles

${\displaystyle M_{pair}}$: Methylation level of a CpG pair

${\displaystyle E_{con}}$: Expected frequency of concordant alleles

${\displaystyle D_{con}}$: Deviation of observed from expected for concordant alleles${\displaystyle M_{pair}=0.5\left(f_{MU}+f_{UM}\right)+f_{MM}}$${\displaystyle E_{con}=\left(f_{MM}+f_{MU}\right)\left(f_{MM}+f_{UM}\right)+\left(f_{UU}+f_{UM}\right)\left(f_{UU}+f_{MU}\right)}$${\displaystyle D_{con}=f_{MM}\text{+}{f}_{UU}-E_{con}=2·\left(f_{MM}·f_{UU}-f_{MU}·f_{UM}\right)}$

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.