DNA Methylation Analysis in Monozygotic Twins Discordant for ALS in Blood Cells

Introduction

Studies investigating DNA methylation in blood cells in ALS have been accumulating in the last few years, indicating the importance of epigenetic studies in the field of motor neuron diseases (MND).^1-3 Despite the presence of numerous genomic sites that show static patterns of DNA methylation, there are also methylation sites in the genome that are regulated dynamically and influenced by lifestyle and environmental factors. Furthermore, comorbidities that are not associated with the neurodegenerative disease can have an impact on the DNA methylome, which makes it rather difficult to detect true changes in methylation associated with neurodegenerative diseases. Additionally, a complex interplay between genetic and epigenetic factors makes it even more difficult to detect the subtle changes in DNA methylation in blood cells associated with neurodegeneration. Our design, given the genetic and clinical nature of the study subjects, provides an ideal setup for matching both genotype and rearing environment within each pair of MZ twins, with the unaffected co-twins being the control group. This study design accounts for not only known potential confounding factors, such as genetics, age, and sex, but latent sources of unwanted variability, as well. ⁴ Overall, the approach we used in this work, also called the case co-twin design, is superior to the ordinary case-control design in terms of the power of a significance test, ⁵ providing a unique opportunity to detect subtle DNA-methylation changes associated with the motor neuron disease.

Methods

Results

Demographics

The current study is composed of 7 pairs of MZ twins (n = 14). Though the initial sample size at the time of study design was 16 (ie, 8 twin pairs), one of the pairs turned out to be dizygotic (DZ) during the verification stage of the zygosity information, and was thereby excluded from the downstream analyses. Of these 7 pairs (Table 1), 3 were female and 4 were male with chronological ages ranging from 38 to 74 years ( x ¯ = 57.43, s = 10.64, x ¯ = 58). The age at onset ranged from 36 to 73 years ( x ¯ = 55.84, s = 11.46, x ¯ = 56.6). The C9orf72 hexanucleotide repeat expansion (HRE) was identified in 3 twin pairs (1 male and 2 female twin pairs), 1 twin pair (female twin pair) showed a variant in the TBK1 (TBK1:c.217A>C, p.I73L) gene, 1 twin pair (male twin pair) was identified with a variant in the NEK1 (NEK1:c.1634T>G, p.M545T) gene and 1 twin pair (male twin pair) was identified with a variant in the NEFH gene (NEFH: c.1054C>A, p.R352S). In total there were 3 spinal ALS cases, 3 bulbar ALS cases, and one ALS-FTD case.

OPEN IN VIEWER

Table 1.

Demographics for the MZ twin set.

Twin pair ID	Disease status	Gender	Age	Mutation	Age of onset	BMI
Twin1	Unaffected	female	62	C9Orf72 HRE	0	22.9
	Affected	female	62	C9Orf72 HRE	60	19.2
Twin2	Unaffected	male	53	NEK1: p.M545T	0	23.9
	Affected	male	53	NEK1: p.M545T	53	20.2
Twin3	Unaffected	male	38	wt	0	23.3
	Affected	male	38	wt	37	19.5
Twin4	Unaffected	male	54	NEFH: p.R352S	0	23.1
	Affected	male	54	NEFH: p.R352S	52	22.3
Twin5	Unaffected	female	58	TBK1: p.I73L	0	28.3
	Affected	female	58	TBK1: p.I73L	57	27.6
Twin6	Unaffected	female	63	C9Orf72 HRE	0	23.1
	Affected	female	63	C9Orf72 HRE	62	18.7
Twin7	Unaffected	male	74	C9Orf72 HRE	0	22.8
	Affected	male	74	C9Orf72 HRE	73	21

Abbreviations: HRE, hexanucleotide repeat expansion; wt, wild type for known ALS genes; BMI, body mass index.

Intra-pair methylation profiles of the mutated ALS genes in the corresponding twin pairs

To determine whether there is any discrepancy in the methylation status of an upstream CpG island spanning the promoter region and exon 1 of each of the mutated ALS genes between affected and unaffected co-twins in the corresponding twin pair, we used the associated probes in the EPIC data and observed a highly similar pattern in the methylation profile of each of these regions per intra-pair comparison (Figure 1). As expected, there is a lack of DNA methylation at CpG islands encompassing the first exons of the TBK1 and NEK1 genes (Figure 1b and c). Importantly, as the probes mapping to CpG islands associated with the C9orf72 gene were filtered out during preprocessing (Supplemental Figure 1), we confined our analysis to this subgroup of ALS genes mutated in our MZ twin set (Table 2). Overall, there appears to be a clear consistency in the methylation profiles of the mutated ALS genes between co-twins in the mutation-positive twin pair, in line with the previous literature. ³⁴

OPEN IN VIEWER

Figure 1.

The promoter proximal sites of neither NEFH nor NEK1 or TBK1 are differentially methylated between mutation-positive ALS-discordant twins. The methylation β-values of both the relative location of the promoter proximal regions and exon 1 for NEFH ((a), upper left), NEK1 ((b), upper right), and TBK1 ((c), bottom) plotted as a function of corresponding probe IDs. Methylation status of each gene was identified using EPIC arrays. Methylation profile of the promoter regions and exon 1 of NEFH, NEK1, and TBK1 does not indicate any pattern of differential methylation between affected and unaffected co-twins in the corresponding twin pairs, concordant for NEFH p.R352S, NEK1 p.M545T, and TBK1 p.I73L, respectively.

OPEN IN VIEWER

Table 2.

Pairwise correlations between individual samples and the gold standard.

Twin ID	Pearson’s R2
Twin1 (Unaff.)	0.905
Twin1 (ALS)	0.916
Twin2 (Unaff.)	0.920
Twin2 (ALS)	0.919
Twin3 (Unaff.)	0.918
Twin3 (ALS)	0.919
Twin4 (Unaff.)	0.914
Twin4 (ALS)	0.914
Twin5 (Unaff.)	0.920
Twin5 (ALS)	0.920
Twin6 (Unaff.)	0.910
Twin6 (ALS)	0.915
Twin7 (Unaff.)	0.916
Twin7 (ALS)	0.915

The gold standard is defined by averaging the beta values across the samples from the largest blood data set. The accepted threshold of R² = 0.800.

Global changes in CpG methylation across the twin pairs

After confirming the twin correlation on global CpG methylation levels using the methylation data, as well as sample sex and twin IDs (Supplemental Figure 2), we explored any relationship between genome-wide DNA methylation and chronological age within- or between MZ twin pairs (Figure 2a). With the mean levels of global DNA methylation being almost the same across all 14 samples, there appeared to be no particular trend associating genomic methylation on a broad scale with age, which also seems to be independent of disease status. Then, we sorted the data by sex and repeated the experiment, only to see the same outcome (Figure 2b): there is likely no relationship between sex and CpG methylation at the genome level, with or without disease status involved. Taken together, our observations on global alterations in CpG methylation across the twin pairs are mostly in agreement with the current literature.

OPEN IN VIEWER

Figure 2.

Global changes in CpG methylation across the twin pairs. (a) An attempt to explore any relationship between genome-wide DNA methylation and chronological age within- or between MZ twin pairs revealed that there is no particular trend associating genomic methylation on a broad scale with aging regardless of disease status. All twin pairs are ranked by chronological age. N(un)affected = 7. (b) It is also clear that there is no relationship between gender and CpG methylation at the genome level, with or without disease status involved. N(un)affected, female = 3, N(un)affected, male = 4. (c) Intra-pair variability analysis performed to identify probes with a sufficient (ie, Δβ-value ⩾ 0.10) intra-pair variability shared by all twin pairs in our data set revealed altered DNA methylation associated with the GRIK1 gene. Also, the twin pairs 1, 6, and 7 accommodate the highest number of pair-specific probes with sufficient intra-pair variability, with 6 being the most extreme. (d) Comparisons of blood cell type composition estimates using neither basic univariate (ie, wilcoxon test) nor more advanced multivariate statistics (ie, the Robust Partial Correlations (RPCs)) approaches reported differentially methylated cell types by disease status. Yet, an increase in B cell and granulocyte levels is apparent in ALS-affected co-twins. ns: not significant. N(un)affected = 7.

Probe-level differential methylation analysis

To analyze the methylation data from ALS-discordant MZ twins at the probe level using the case co-twin design, we used a mixed-effects model with fixed effects of age, sex, C9orf72 status, and cell type estimates, with random effects of batch (Supplemental Figure 3). As indicated above for biological age, though discordant MZ twins are matched for biological factors specific to individual pairs, the effects of these factors on the within-pair difference in CpG methylation patterns should be accounted for,^4,5 which is addressed in the Discussion section. Considering that the intercept is the key term of the analysis as the mean level of intra-pair methylation differences per probe when all fixed effect variables are set to zero, we plotted the frequency distribution of the intercept estimates obtained from the regression analysis. In line with the literature, ³⁸ they followed a clear pattern of normal distribution (Figure 3a), just as the M-values. Ultimately, a total of 35 probes were found significant at an FDR ⩽ 0.05 with varying levels of significance (Supplemental Figure 4). In an attempt to visually assess their degree of variation in methylation levels by disease status, we observed a heterogeneous profile of CpG methylation for the DMPs (Figure 3b). Subsequently, each of the underlying assumptions of our linear model was evaluated using a selected subset of probes for the quality of the model fit, performance characteristics, and diagnostic accuracy. From a technical aspect, a remarkable robustness of LMMs has been confirmed in the literature,^39,40 except for small biases in estimated parameters stemming from missing predictors or random effect components, even if the majority of the assumptions are violated. Our results revealed, though using a subset of probes, that the majority of the assumptions were met (Supplemental Figures 5 and 6). Next, we checked to see whether these DMPs are capable of distinguishing between ALS and non-ALS samples. The results of different methodologies converged on a low power of discrimination between sample groups given just these significant probes (Figure 3c and d). When the DMP-annotated genes were clustered based on the direction of the change in methylation, 2 groups emerged: (a) ATG9A, BCAT1, BARX2, C1orf103, CYGB, FES, MYH9, PARVB, and RNF4 as hypomethylated genes, (b) ATP6V0B, BOC, BUD31, C3orf75, FEZ2, IL17D, KILLIN, MEGF6, OSBPL10, PDAP1, PTEN, RSPO4, SLX4, UBLCP1, and USP20 as hypermethylated. Though not significant (FET P-value: .0699), protein association network analysis of these 24 genes performed using the STRING platform revealed the enrichment of a part of the ubiquitin (Ub)-proteasome pathway (UPP) of protein degradation, which is known to be dysfunctional in ALS.^41,42 On the other hand, “autophagy” and “axon guidance” were identified as the top deregulated signaling pathways (q-value = 0.012 and 0.020, respectively), and “actin cytoskeleton reorganization” was identified as the top enriched biological process (q-value = 0.003) in these DMP-annotated genes which has been strongly associated with ALS pathogenesis.^43-45 To detect differences up to 20% and 50% in CpG methylation between affected and unaffected co-twins with a minimum number of 7 subjects in a sample group based on the power of sample calculations, the calculated power of making a statistically significant comparison at an FDR ⩽ 0.05 turned out to be low, as expected. Yet, it was also observed for target differences of 20% and 50% that the probability of detecting at least 1 true positive significant probe out of 35 DMPs identified between affected and unaffected co-twins is about 10% and 61%, respectively given the current sample size (Supplemental Figure 7). In other words, approximately 4 of the DMPs identified within the scope of this study are expected to be true positive results with ⩾20% methylation difference between sample groups. Interestingly, our ROC analysis revealed using a classical AUC threshold of 0.70 that there are exactly 4 probes (Figure 4a) out of 35 DMPs (Supplemental Figure 8) that appear to classify the sample groups well. While 3 of these candidate probes are located in the “OpenSea clusters” of the genome, one of them with an AUC value of 0.755 is located near the TSS of the IL17D gene, which is implicated in positive regulation of cytokine production during inflammatory response ⁴⁶ (Figure 4b). Overall, our probe-level differential methylation analysis appears to provide a small but biologically relevant group of genes that are differentially methylated within the context of ALS.

OPEN IN VIEWER

Figure 3.

Probe-level differential methylation analysis identified 35 DMPs with varying levels of methylation changes between conditions. (a) As a QC step for the statistical approach (ie, a linear mixed-effects model) used in the differentially methylated probe analysis, the frequency distribution of the intercept estimates obtained from the regression analysis follows a clear pattern of normal distribution. (b) The degree of variation by disease status in methylation levels of the 35 DMPs identified within the scope of this work is highly variable: while there are DMPs with clear differences in methylation levels between ALS and non-ALS samples, there are also DMP that have roughly similar levels of methylation by disease status when raw βvalues are used for visual assessment, implying the importance of adjusting for key covariates such as age and gender in estimating the relationships between CpG methylation patterns and disease status. N(un)affected = 7. (c and d) When we tested whether these DMPs are capable of distinguishing between ALS and non-ALS samples, the results of different methodologies (ie, a heatmap (c) or a PCA (d) of raw methylation levels) indicated a low power of discrimination between sample groups given just the significant probes. In other words, methylation was roughly similar across samples for most DMPs, and samples did not seem to cluster by disease status. PC1 in (d) explains 34% of the variation in the input data of DMPs.

OPEN IN VIEWER

Figure 4.

In silico validation of DMP analysis results indicate differential methylation of IL17D in ALS blood. (a) The ROC analysis results of 35 DMPs using a classical cutoff value at AUC ⩾ 0.70 identified 4 probes that appear to classify the sample groups well, just as expected after the results of the power of sample calculation. The corresponding probe IDs are shown in red at the top of each plot. N(un)affected = 7. (b) Graphical representation of the DMP with the ID number “cg13993853” overlapping the IL17D promoter region, as indicated by the adjacent activating epigenetic marks (eg, CpG islands, H3K27Ac, and DNaseI hypersensitivity), as well as the neighboring TF binding sites. Unlike the other 3 candidate probes that are located in the “OpenSea clusters” of the genome, this probe with an AUC value of 0.755 is the only candidate probe that is located near the TSS of a gene.

Discussion

In this study, our goal was to provide further biological insight into ALS pathogenesis at the DNA methylation level using whole blood samples of ALS-discordant MZ twins. We have shown that a few CpG locations in the genome carry differential methylation marks associated with ALS. In addition, we observed differential methylation of neither genomic bins of different sizes nor blood cell types associated with ALS in our data set. Nevertheless, global patterns of genomic methylation across the twin pairs, including an increasing intra-pair discrepancy in methylation age and a decreasing intra-pair correlation on CpG methylation with chronological aging, a low degree of variation in within-pair CpG methylation, and elevated levels of B cells and neutrophils in ALS-affected co-twins, were all in agreement with the current literature. Our results imply that the methylation landscape of whole blood in ALS is a dynamic and highly heterogeneous network at the genome level, with only some commonly found perturbations and affected pathways contributing to the disease.

As rarely as these study subjects can be found, our twin set presumably harbors invaluable information on the molecular etiology of ALS, independent of the genetic mutations in the ALS gene each twin pair carries. More specifically, the ultimate reason why we chose a study design as such, that is, with identical twins that have remarkably different phenotypes, was because co-twins in a pair are matched for familial factors, including but not limited to baseline genetic sequence, age, sex, and rearing environment, thereby reducing the variability between pairs of subjects, but are still discordant for the disease. This leaves room for epigenetic components that might be involved in the disease pathogenesis. We observed rather consistent methylation profiles of the mutated ALS genes between co-twins in the mutation-positive twin pair, implying causative mutations detected in each twin pair are independent of differential CpG methylations identified within the scope of this study.

Another limitation of the current work relates to the LMM that was found to produce heavily inflated P-values (likely a high number of false positives ⁴⁷ ) within the context of the probe-level differential methylation analysis. As statistically precarious as this situation may be, such inconsistencies can be explained in the context of our experimental design, as well as our analysis workflow, independent of modeling errors. However, it should also be noted that as we preferred the case co-twin design to the classical twin design for DMP analysis chiefly due to our initial experimental design, the whole procedure of data preparation and processing mentioned in the “Differentially methylated probe (DMP) analyses” part of the Methods section, composed mainly of transformations, value adjustments, selection of predictor variables (both fixed and random effect covariates) in the LMM, model fitting, and extraction of key metrics for downstream analysis, has been planned and performed based on the developers’ guidelines. ⁵ This is also evident by the total number of DMPs identified in this work, which is in line with many other discordant twin studies in the literature.^34,48,49 Back to the potential reasons for inflated P-values, outlier probes in a data set are widely known to cause instability in statistical tests, which can result in exceptionally high standard errors, and thereby, inflated P-values. The preprocessQuantile() function used for data normalizations during preprocessing stage of the analysis workflow also fixes outliers, but only those that have low (ie, close to zero) signal intensities, not high. In this regard, it is likely that the significant results are, to some extent, driven by technical artifacts or outliers. In fact, only 4 of the 35 DMPs identified in this study could be validated using ROC analysis, just as predicted using the power of sample calculations as the total number of true positive findings from DMP analysis with ⩾20% methylation difference between sample groups. This is why we decided to prioritize the genomic locations associated with these 4 DMPs over those associated with all DMPs while making biological interpretations of our results. Another point is that we already performed a comprehensive search trying most of the combinations of predictor variables to find the subset with the best evaluation criterion using not only simple linear model but also mixed-effects model, with the MDS/PCA results in mind, before constructing the final version of our multivariate regression model. In this respect, singularity issues and multicollinearity were also taken into account as selection criteria. Lastly, our power calculations demonstrated that the small sample size we had for the EWAS presented here made it rather difficult to obtain any significant results (Supplemental Figure 7). So, there is also a chance, though not the most likely scenario, that small sample sizes lead to such statistical biases as an inflated type-1 error, as reported in many different clinical trial studies before. ⁵⁰ Taken together, the statistical approach used in this work, though not minimizing the risk of false positive findings being reported together with true positive results, allowed us to identify a small subset of DMPs that were later validated bioinformatically.

An overly underestimated concept in discordant MZ twin studies is the intra-pair dynamics of epigenetic drift. Epigenetic drift can be defined as the divergence of the epigenome with biological age due to stochastic alterations in methylation. ⁵¹ This paradigm has often been overlooked especially when the methylation changes are investigated, with chronological age being part of the equation, in identical twins who do not share the disease phenotype between co-twins. It is critical to note that the effect of the age-associated intra-pair difference in DNA methylation status must not be disregarded in such studies though the biological factors specific to individual pairs are matched out in MZ twins. ⁵² Empirical evidence confirming this effect came from the observation that twin correlation on genome-wide DNA methylation declines with age (Supplemental Figure 2). If this increasing divergence of the epigenome with age actually roots from stochastic alterations, then co-twins in identical twin pairs should be affected differently along the natural course of epigenetic drift, leading to diminishing within-pair twin correlations on global DNA methylation. From a technical perspective, this approach might as well serve as a QC check in MZ twin studies and be used to monitor sample compliance.

Considering our analysis approach for probe-level differential methylation, one could easily claim that age and sex are matched out in the case co-twin design and that including them as predictor variables in the final model seems excessive. Yet, such an argument does not apply to the EWAS as: (a) intra-pair differences in the epigenome are likely to be higher in older twin pairs, just as mentioned above, (b) it has been known that intra-pair epigenetic differences in particular regions of the genome vary by sex. ⁵ As a result, fixed effect covariates as such can improve the regression analysis considerably when included in the model. From a practical standpoint, the regression models we tested in the absence of such predictors did not seem to yield any better results, either. Further examples can be found in other studies with similar experimental designs in the literature.^4,34,48

The only probe with sufficient (ie, Δβ-value ⩾ 0.10) intra-pair variability shared by all twin pairs in our data set, cg26023019, was found to be hypermethylated in affected co-twins and map to GRIK1 which is involved in the excitatory neurotransmitter receptors in the mammalian brain. The glutamate signaling pathway is known to be intimately related to both acute and chronic neurodegenerative disorders, including ALS. ⁵³ In this regard, several lines of evidence also suggest that excessive activation of the glutamate receptors results in excitotoxicity.^53-55 With the GRIK1 gene being hypermethylated (or presumably transcriptionally suppressed) in ALS samples, our findings implicate the activation of a compensatory mechanism responsible for restoring the overly active glutamate signaling back to normal levels in ALS. Further research is needed to confirm this hypothesis and to understand the potential role of GRIK1 in ALS biology.

Finally, 1 of 4 DMPs identified within the scope of this study and later verified by ROC was found to be hypermethylated in affected co-twins and located near the TSS of the IL17D gene, which is implicated in positive regulation of cytokine production during an inflammatory response. Interestingly, recent evidence has shown that IL17 directly affects human motor neuron survival and directly contributes to motor neuron degeneration in ALS.^56,57 Considering that IL17D is vastly expressed in various brain regions ⁵⁸ and that blood is a reasonable surrogate for brain tissue in DNA methylation studies, ⁵⁹ IL17D has the potential to serve as a blood-based biomarker.

Conclusion

Despite the small sample size of the current study, we identified a group of DMPs that might serve as potential epigenetic loci for future ALS studies. Especially for the studies investigating different phenotypes (eg, fast vs slow ALS progressors) irrespective of the genetic cause underlying the disease, testing for these DMPs might provide meaningful biological outcomes and could further validate the findings of the current study.

Supplemental Material