Examining the Vanishing Twin Hypothesis of Neural Tube Defects: Application of an Epigenetic Predictor for Monozygotic Twinning

Abstract

Strong associations between neural tube defects (NTDs) and monozygotic (MZ) twinning have long been noted, and it has been suggested that NTD cases who do not present as MZ twins may be the survivors of MZ twinning events. We have recently shown that MZ twins carry a strong, distinctive DNA methylation signature and have developed an algorithm based on genomewide DNA methylation array data that distinguishes MZ twins from dizygotic twins and other relatives at well above chance level. We have applied this algorithm to published methylation data from five fetal tissues (placental chorionic villi, kidney, spinal cord, brain and muscle) collected from spina bifida cases (n = 22), anencephalic cases (n = 15) and controls (n = 19). We see no difference in signature between cases and controls, providing no support for a common etiological role of MZ twinning in NTDs. The strong associations therefore continue to await elucidation.

Neural tube defects (NTDs) are a group of congenital malformations affecting around 1 per 1000 human pregnancies worldwide, although large variation exists across countries and time (Zaganjor et al., 2016). The most common forms are anencephaly and spina bifida. A strong association between monozygotic (MZ) twinning and NTDs—in particular anencephaly—has long been noted (Clarke et al., 1975; Holmes, 2012; Källén et al., 1994; Windham & Sever, 1982). The frequency of MZ twins in the general population is 0.3%−0.4% (Bulmer, 1970; Smits & Monden, 2011) but has been reported to be increased in NTD patient groups (Holmes, 2012). A summary of multiple international databases indicated that the rate of twinning was highest among infants with anencephaly (4.3%), followed by infants with encephalocele (3.3%), lumbosacral spina bifida (2.2%) and thoracic spina bifida (1.2%), and that twins were mainly MZ when zygosity had been determined (Källén et al., 1994). Though reports of exact percentages of discordant and concordant MZ pairs are lacking, it has been noted that concordant monoamniotic twins are rare, suggesting a high frequency of discordant MZ pairs (Holmes, 2012). In addition to the comprehensive multinational study by Källén et al., 1994, multiple additional papers reporting NTD and twinning rates have been published, including more recent articles, with mixed findings with respect to the association between NTDs and MZ twinning (Dawson et al., 2016; Holmes, 2012; Parker et al., 2010; Yu et al., 2019). Important factors that may cause inconsistencies include small sample size, whether rates were obtained before or after the introduction of folate supplementation during pregnancy, inclusion of children conceived with/without assisted reproductive technologies/fertility treatment, inclusion of live births versus fetal deaths and elective terminations, inclusion of individuals with chromosomal abnormalities, inclusion of higher order multiples in addition to twins, consideration as to whether the NTD is part of a syndrome or occurs in isolation (isolated vs. associated NTD cases) and distinction of subtypes (e.g., thoracic versus lumbosacral).

Multiple hypotheses have been proposed regarding mechanisms linking MZ twinning to NTDs. Clarke et al. (1975) suggested that unfavorable interaction between twin fetuses or between residual trophoblastic material from a previous miscarriage and another fetus may cause anencephaly and/or spina bifida in the other fetus. James (1976, 1981) proposed that the comorbidity could be explained by delays at critical stages of zygote development as a result of the MZ twinning event. He argued that separation of the embryo into two halves gives rise to two embryos with a slight developmental delay, which could cause the embryos to be deficient in oxygen or nutrients, leading to developmental arrest or NTDs, an earlier delay causing anencephaly and a later delay to spina bifida. Similarly, Myrianthopoulos (1978) hypothesized that the MZ twinning process disrupts the developmental clock of the two resulting embryos, creating disadvantages that make them more susceptible to teratogens. Noteworthy, recent studies of artificially produced human MZ twins by splitting of blastomeres indeed suggest that the resulting twin zygotes are developmentally delayed (Noli et al., 2015), but it is unknown whether this also applies to naturally occurring MZ twins. Schinzel et al. (1979) suggested that MZ twinning and early malformations share a common etiology, that is, an early insult may cause the separation of the embryos and lead to additional morphologic problems. Such a shared causal factor could be genetic or environmental; proposed environmental contributors include a disturbance in folate or other nutritional factors (Rivas et al., 1995). Cloacal exstrophy, which is often comorbid with spina bifida and associated with a high frequency of recorded vanishing twins (Dawson et al., 2016; Siebert et al., 2005), has been attributed to conjoined twinning, producing one affected individual and one vanishing twin (Casale et al., 2004). The following hypothesized mechanisms explicitly suggest that MZ twins will be discordant for NTDs or that affected MZ twins will have a vanishing co-twin: the unfavorable interaction between twin fetuses hypothesis (Clarke et al., 1975) and the conjoined twinning hypothesis (Casale et al., 2004). The developmental delay model by James (1976, 1981) and Myrianthopoulos (1978), supported by artificially produced MZ twins (Noli et al., 2015), and the common etiology model (Schinzel et al., 1979) are expected to produce NTD concordant-affected MZ twin embryos (who may or may not both survive), unless the MZ twinning event is somehow associated with a developmental asymmetry, causing only one of the embryos to experience a developmental delay or insult (thus producing discordant MZ twins).

These early hypotheses lead to a major new hypothesis: if some biomarker is found to distinguish MZ twins from other individuals, then at least a proportion of NTD cases (those where the cause of the disorder is connected to MZ twinning) would bear this biomarker. That is, individuals who are not recognized as the sole survivor of an MZ twin pregnancy would be identified as MZ twins. Such a biomarker has just been found in the form of a DNA methylation signature that can be derived from Illumina 450k and Illumina EPIC array data that predicts monozygosity well above chance, although not perfectly (van Dongen et al., in press). The signature was discovered in blood samples from adult twins (mean age = 34.9, SD = 11.3) and showed strong replication in blood samples from independent adult twin cohorts, and in buccal samples from twin children (mean age = 9.6, SD = 1.8). A prediction based on early hypotheses is therefore that NTD cases should have a methylation signature closer to that of MZ twins than of unaffected controls.

We have tested this hypothesis in anencephalic (n = 15), spina bifida (n = 22) and control (n = 19) fetuses sampled for five different tissues (four in the case of anencephalics), namely brain, kidney, muscle, placental chorionic villi and spinal cord (Price et al., 2016). None of the pregnancies were known to be a multiple pregnancy or vanishing twin pregnancy, but it should be noted that for vanishing twins to be recorded and/or noticed, the vanishing twin must at least survive up until the first Doppler or ultrasound screening, while twins who vanish earlier will go unnoticed. All samples were typed using the Illumina 450k methylation arrays, as were the samples from MZ, dizygotic (DZ) twins and their relatives on which the monozygosity predictor was trained. Two versions of the predictor were evaluated, one optimally distinguishing MZ from DZ twins (predictor A) and the other distinguishing MZs from all others (predictor B).

Methods

We previously detected a robust DNA methylation signature of MZ twinning based on whole-blood Illumina 450k array data from adult twins, which showed remarkable replication in four independent twin cohorts with whole-blood methylation data (Illumina 450k array) from adult twins and in buccal DNA methylation data (EPIC array) from twin children (full details in van Dongen et al., in press). In a subsequent meta-analysis of the whole-blood datasets (N = 5723), 834 differentially methylated positions (DMPs) associated with MZ twinning were detected. The robustness of these results suggested that a DNA methylation-based score for MZ twinning could open up new avenues to investigate the link between vanishing MZ twins and various congenital disorders (including NTDs) through an epigenetic biomarker. To this end, we previously trained a DNA methylation-based classifier of MZ twinning using penalized regression models (elastic net) on blood Illumina 450k array DNA methylation data from the Netherlands Twin Register (NTR). We compared models based on two input sets (genomewide methylation sites vs. meta-analysis DMPs) and trained on two phenotypes: MZ versus DZ twins, and MZ twins versus everyone else (including DZ twins and family members of twins). Regressions returned predictors based on 232–1867 methylation sites. In NTR test data from blood (which were left out of the training data set), the area under the curve (AUC) ranged from 0.69 to 0.77, with up to 84% of MZ twins correctly classified, up to 57% of DZ twins correctly classified and up to 63% of family members correctly classified as non-MZ. We tested prediction in two independent data sets: Brisbane Systems Genetic Study (BSGS) adults (blood, 450k array) and NTR children (buccal, EPIC array). AUCs ranged from 0.67 to 0.80 in BSGS, and from 0.63 to 0.76 in buccal data from NTR children. In all data sets, the predictors trained on the smaller set of CpGs from the meta-analysis performed best.

We next sought a data set with methylation data on NTDs and identified the study reported by Price et al., 2016. Second trimester (14–26 weeks) human placental chorionic villi, kidney, spinal cord, brain and muscle were collected from 19 control, 22 spina bifida and 15 anencephalic fetuses in British Columbia, Canada. DNA was extracted and typed using the Illumina Infinium Human Methylation450 (450k) array, making it comparable with the twin samples analyzed in the zygosity studies (van Dongen et al., in press). The filtered and normalized DNA methylation beta-values (see Price et al., 2016, for detail) were downloaded from NCBI’s gene expression omnibus (GEO; accession number GSE69502). This data set consisted of 179 samples: 52 chorionic villus samples, 44 kidney samples, 32 spinal cord samples, 20 brain samples and 31 muscle samples and contained 442,091 methylation sites. We discarded methylation sites that had missing values in any of the 179 samples from NTD cases or controls (remaining n methylation sites = 354,030). Next, we fine-tuned the previously described predictors of MZ twinning by training them on the subset of CpGs detected in our previous meta-analysis of MZ twinning that were also present in the NTD data set (n = 692 methylation sites); resulting prediction models included 213 methylation sites for predictor A (trained on data from MZ and DZ twins) and 242 methylation sites for predictor B (trained on data from MZ twins, DZ twins and relatives) and had AUCs of 0.75 (Supplemental Figure 1) and 0.76 (Supplemental Figure 2), respectively, in the test data from the NTR (whole blood).

As previously described, DNA methylation beta values of individual CpGs were standardized (z scores) in training and test data sets. For the NTD data set, standardization was applied to each CpG within each tissue group separately before applying the predictors (similar to the training data set, where standardization was also applied within tissue, i.e., whole blood, to express within-tissue individual differences). Next, the epigenetic predictors for MZ twinning were applied to the three groups (controls, spina bifida and anencephaly) and five tissues (four for anencephaly; note that brain is absent for the anencephalics) to derive a continuous DNA methylation score for each sample. Note that confirmed MZ twins exhibit a higher MZ-DNA methylation score, see Supplementary Figures S1 and S2 based on data from the NTR. We tested whether the MZ-methylation score differed between NTD cases and controls, while correcting for tissue using generalized estimation equation (gee) models, with DNA methylation score as outcome variable, and NTD diagnosis (case vs. control) and tissue (dummy-coding) as predictors. A second model additionally included the interaction between diagnosis and tissue, to test whether a difference in methylation score was evident within (a) particular tissue(s). Gee models were fitted with the R package gee, with the following specifications: Gaussian link function (for continuous data), 100 iterations, and the ‘exchangeable’ option to account for the correlation structure within fetus (different tissues from the same fetus). Four gee models were run (one for predictor A, one for predictor B, one model comparing spina bifida cases versus controls and one model comparing anencephaly cases versus controls). Based on previous reports of higher rates of MZ twins among NTD cases (especially for anencephaly), our hypothesis was that NTD cases would exhibit a higher average MZ-methylation score compared to controls, with a more pronounced difference in anencephaly. Because none of the NTD samples were derived from known twin pregnancies, the hypothesized pattern would be in line with NTDs being associated with a vanishing MZ twin.

Results

Boxplots are shown for each tissue (brain, kidney, muscle, placental chorionic villi, spinal cord) in each of the three patient groups (anencephalics, controls, spina bifida) for predictor A in Figure 1 and for predictor B in Figure 2.

Fig. 1. Results Predictor A (trained to distinguish monozygotic [MZ] from dizygotic [DZ] twins). Predictor A was trained on whole-blood DNA methylation data from the Netherlands Twin Register to distinguish MZ twins from DZ twins and is based on 213 methylation sites. The predictor was applied to second trimester samples from NTD fetuses and controls. The boxplots show the distribution of continuous MZ-methylation scores based on Predictor A in five tissues from NTD and control fetuses.

Fig. 2. Results Predictor B (trained to distinguish monozygotic [MZ] twins from all others including parents/sibs). Predictor B was trained on whole-blood DNA methylation data from the Netherlands Twin Register to distinguish MZ twins from all others (dizygotic twins + nontwins) and is based on 242 methylation sites. The predictor was applied to second trimester samples from NTD fetuses and controls. The boxplots show the distribution of continuous MZ-methylation scores based on Predictor B in five tissues from NTD and control fetuses.

Visual inspection of these plots suggests a remarkable consistency in mean methylation scores for a given tissue between the three patient groups, despite the small sample sizes, suggesting no strong methylation signature for MZ twinning in the two NTD groups. This is confirmed in statistical analysis. There was a slight effect for predictor A in the expected direction (higher MZ-methylation score in NTD cases), though not significant (Table 1). For predictor B, the effect went in the wrong direction (ns, Table 2). Results from interaction models (tissue × diagnosis) are presented in Supplementary Tables S1 (predictor A) and S2 (Predictor B). Visual inspection of Figures 1 and 2 further suggests no difference in the number of outliers (samples with an extreme MZ epigenetic score) between NTDs and controls, providing no support for the presence of a strong MZ-methylation signature in a subset of NTD cases.

Table 1. Results from gee model with main effects of diagnosis and tissue for Predictor A (trained to distinguish MZ from DZ twins)

Table 2. Results from gee model with main effects of diagnosis and tissue for Predictor B (trained to distinguish MZ twins from all others including parents/sibs)

Discussion

This work does not provide support for the vanishing MZ twin hypothesis of NTDs, but we note many limitations that may temper this negative conclusion—most importantly, the relatively small sample size, the imperfect predictive performance of our MZ twinning algorithm and unknown performance in tissues other than blood and buccal, and that both MZ twinning and NTDs may have multiple, heterogeneous causes so that MZ twinning and NTDs may be connected in some cases but not all (i.e., only a proportion of the NTD cases included in our study may have resulted from a vanishing MZ twin pregnancy, while others may be true singletons). Most if not all twin cases with NTDs recorded in the literature are the result of a pregnancy where both co-twins survived at least up till a certain stage (some reports include only live births, others also include selectively terminated pregnancies), and it has been suggested that the frequency of discordant MZ pairs is high (Holmes, 2012). Since the patients included in our study are not from a known twin pregnancy, they are either not twins or they resulted from twin pregnancies in which only one embryo survived the early stages of pregnancy. We do not know if such cases may represent more severe or even atypical twin cases.

Another possibility is that the DNA methylation signature we have developed is a post facto product of successful MZ twinning and does not describe or characterize the state of fission, nor cause it. Further, the DNA methylation signature was discovered in healthy MZ twins who resulted from a successful splitting event. Possibly, NTDs, and perhaps in particular cases with a vanishing twin, may be linked to an imperfect splitting event and a weaker MZ epigenetic signature.

Larger samples, dissection of heterogeneity, including fetal deaths, and more precise methylation algorithms trained on multiple tissues may allow us to tease these possibilities apart. Future epigenetic studies may also address the hypothesis of a delayed biological clock as a mechanism connecting NTDs and MZ twinning. To this end, epigenetic age or epigenetic mitotic clock algorithms could be applied to tissues from NTDs and MZ twins to test the hypothesis that their clocks are delayed. Based on data from the same NTD samples included in the current study, a previous study on placental gestational age epigenetic clocks showed that NTD placentas show no difference in epigenetic aging (Lee et al., 2019).

We note that although we found no support for a strongly increased MZ-methylation score in NTD fetuses, our previously reported epigenetic findings in MZ twins did point to possible connections with NTDs (van Dongen et al., in press). First, MZ-DMPs were enriched near genes involved in cell adhesion including (proto)cadherins, which have also been implicated in closure of the neural tube (Copp & Greene, 2010). Second, we found enrichment of MZ-DMPs in transcription factor motifs involved, among other processes, in ‘tube development’. Third, MZ-DMPs were enriched for CpGs that had been previously associated with maternal folate, which is the most characterized environmental risk factor for NTDs. Finally, one of the top differentially methylated loci in MZ twins encodes Psk9, a protein that was identified as a putative serum biomarker for NTDs (An et al., 2015).