Epigenome-wide analysis in newborn blood spots from monozygotic twins discordant for cerebral palsy reveals consistent regional differences in DNA methylation

Background Cerebral palsy (CP) is a clinical description for a group of motor disorders that are heterogeneous with respect to causes, symptoms and severity. A diagnosis of CP cannot usually be made at birth and in some cases may be delayed until 2–3 years of age. This limits opportunities for early intervention that could otherwise improve long-term outcomes. CP has been recorded in monozygotic twins discordant for the disorder, indicating a potential role of non-genetic factors such as intrauterine infection, hypoxia-ischaemia, haemorrhage and thrombosis. The aim of this exploratory study was to utilise the discordant monozygotic twin model to understand and measure epigenetic changes associated with the development of CP. Methods We performed a genome-wide analysis of DNA methylation using the Illumina Infinium Human Methylation 450 BeadChip array with DNA from newborn blood spots of 15 monozygotic twin pairs who later became discordant for CP. Quality control and data preprocessing were undertaken using the minfi R package. Differential methylation analysis was performed using the remove unwanted variation (RUVm) method, taking twin pairing into account in order to identify CP-specific differentially methylated probes (DMPs), and bumphunter was performed to identify differentially methylated regions (DMRs). Results We identified 33 top-ranked DMPs based on a nominal p value cut-off of p < 1 × 10−4 and two DMRs (p < 1 × 10−3) associated with CP. The top-ranked probes related to 25 genes including HNRNPL, RASSF5, CD3D and KALRN involved in immune signalling pathways, in addition to TBC1D24, FBXO9 and VIPR2 previously linked to epileptic encephalopathy. Gene ontology and pathway analysis of top-ranked DMP-associated genes revealed enrichment of inflammatory signalling pathways, regulation of cytokine secretion and regulation of leukocyte-mediated immunity. We also identified two top-ranked DMRs including one on chromosome 6 within the promoter region of LTA gene encoding tumour necrosis factor-beta (TNF-β), an important regulator of inflammation and brain development. The second was within the transcription start site of the LIME1 gene, which plays a key role in inflammatory pathways such as MAPK signalling. CP-specific differential DNA methylation within one of our two top DMRs was validated using an independent platform, MassArray EpiTyper. Conclusions Ours is the first epigenome-wide association study of CP in disease-discordant monozygotic twin pairs and suggests a potential role for immune dysfunction in this condition. Electronic supplementary material The online version of this article (10.1186/s13148-018-0457-4) contains supplementary material, which is available to authorized users.


Background
Cerebral palsy (CP) describes a group of motor impairment syndromes caused by lesions or anomalies of the developing brain [1]. It is non-progressive, but the severity of symptoms may change over time [2]. CP is the most common childhood physical disability [3] with a worldwide prevalence of 2.11 per 1000 live births [4]. In preterm infants (< 37 weeks gestation), the prevalence is higher, ranging from 5 to 92 per 1000 depending on gestational age [5]. The prevalence of CP in multiple births is almost four times that of singletons [6]. There are many factors that may be responsible for this increased risk in multiple birth pregnancies, with the most likely being low birth weight and preterm birth, both known risk factors for CP. The brain insult or anomaly resulting in CP may occur during the prenatal, perinatal or early postnatal period [1,7], and in many cases, the timing is unknown. Although newborns may be recognised as being at risk of CP, less than half of all children who are ultimately diagnosed are identified before 1 year of age, and only three quarters are identified before age 2 [8].
CP has a multifactorial pathogenesis and risk factors including intrauterine growth restriction or infection, placental abnormalities, inflammation, signs of fetal distress and genetic variation [9,10]. Although the latter explains a proportion of CP cases, particularly cerebral maldevelopments [9,11,12], many non-genetic factors likely play a role [7] though their respective contributions have not been comprehensively addressed [13][14][15].
Genetically identical monozygotic (MZ) twin pairs discordant for CP highlight the role of non-shared factors in the pathogenesis of CP [16]. Non-shared factors are often described as the differences in the intrauterine environment that influences the development of individual members of a twin pair [17,18]. Such within-pair variation can arise from differences in the length and morphology of the umbilical cord or placenta. This can affect growth rate and development of the individual twins leading to discordance in infection or inflammation, leading to disease discordance [13,19,20]. More generally, the study of phenotypically discordant MZ twins, matched for genetic variation, sex, gestational age and maternal factors, provides a great opportunity to examine the role of epigenetics in disease aetiology by allowing us to isolate the effect of such non-shared environmental factors [21,22].
Epigenetics refers to a range of modifications and processes that regulate the activity of DNA, including gene expression. Epigenetic variation has emerged as a candidate mediator of a range of health outcomes beginning in early life as part of the 'Developmental Origins of Health and Disease' (DOHaD) phenomenon. [23][24][25]. In DOHaD, environmental conditions in utero and during infancy alter the developmental trajectory of an individual, which manifests as specific chronic health phenotypes later in life.
DNA methylation is the most widely studied epigenetic process and is one of the mechanisms that are involved in tissue differentiation during early development. Despite this, previous studies have investigated the concordance in DNA methylation state between the brain and peripheral tissues, revealing many similarities [26][27][28]. Several epigenome-wide association studies (EWAS) have identified DNA methylation variation in the cord blood in association with later neurocognitive function and behaviour [29][30][31]. Similarly, the whole blood has been used to detect differential methylation patterns between affected and unaffected individuals in brain disorders such as schizophrenia [27], bipolar disorder [32,33] and Alzheimer's disease [34]. Animal studies have also reported that environmental factors affecting brain processes leave biomarker signatures in the blood with consistent methylation status across the brain and peripheral tissues [35]. We have previously shown that DNA methylation varies within pairs of MZ twins from birth [36,37]. In this study, we hypothesised that early life non-shared factors that play a role in the aetiology of CP in discordant MZ twins may be reflected in differences in DNA methylation across tissues including neonatal blood. Furthermore, we hypothesised that a subset of differential methylation will be located in genes previously implicated in CP aetiology [38][39][40][41][42], in particular, pathways involved in inflammation, hypoxiaischaemia and thrombosis.

Samples and DNA extraction
Participant CP-discordant twin pairs, who were suspected to be MZ on the basis of same sex and questionnaire data that measured concordance for hair and eye colour, were recruited through the Victorian Cerebral Palsy Register, a population-based registry of individuals born or receiving medical services in the Australian state of Victoria. Participants were excluded if either twin presented with another known neurological disorder. Written informed consent from the families was obtained. A single 1-cm-diameter dried blood spot was acquired from all participants from a neonatal newborn screening card collected 2 to 4 days after birth and then stored within the Victorian Clinical Genetics Service.
Genomic DNA (gDNA) was extracted from the dried blood spot samples using the ZR DNA Card Extraction Kit (Zymo Research, Irvine, CA, USA) with some modifications to the manufacturer's protocol. Briefly, eight 3-mm punches were from each 1-cm-diameter blood spot and were transferred to a 2-mL Eppendorf Safe-Lock microcentrifuge tube (Merck, Darmstadt, Germany) containing ZR BashingBeads. Four hundred microliters of PBS containing 40 μL of 20 mg/mL proteinase K (Sigma-Aldrich, St. Louis, Missouri, USA) was added, and samples were vortexed and centrifuged briefly, followed by incubation overnight at 37°C. Following incubation, 400 μL of ZR lysis solution was added to each tube. Punches were homogenised for 30 s at 4 m/s 2 using Thermo Savant FastPrep 120 Cell Disrupter System (Global Medical Instrumentation (GMI) Incorporation, Minnesota, USA). Tubes were centrifuged for 1 min at 10,000 rpm, and 390 μL of 2× digestion buffer and 10 μL of 20 mg/mL proteinase K were added. Tubes were mixed by inversion and incubated for 30 min at 55°C, then left to cool at ambient temperature for 3-4 min before centrifuging for 1 min at 8000 rpm. Six hundred fifty microliters of supernatant was added to 1.3 mL of DNA isolation buffer contained in a 5-mL Falcon tube (Thermo Fisher Scientific, MA, USA). This mixture was passed through the Zymo-Spin IC column by centrifuging for 1 min at 14,000 rpm, followed by the discard of flow-through liquid. The spin column was then washed twice by adding 200 μL of DNA wash buffer and centrifuged for 1 min at 14,000 rpm. Finally, 20 μL of DNA elution buffer (prewarmed at 55°C) was added to the column and incubated at ambient temperature for 15 min before final centrifugation for 2 min at 14,000 rpm. This was repeated, resulting in a final elution volume of 40 μL containing genomic DNA. DNA concentration was measured by spectrophotometry (Nanodrop, Wilmington, DE, USA) to allow calculation of the required volume of each sample for array analysis. The quality of the extracted gDNA samples was visualised using agarose gel electrophoresis.

Illumina Infinium HumanMethylation450 arrays
Following bisulphite conversion of genomic DNA, genomewide analysis of DNA methylation was assessed using HM450 (Illumina, San Diego, CA, USA), at the Department of Pathology, University of Melbourne. Hybridisation and scanning were performed following the manufacturer's instructions. Statistical analysis was performed using the R statistical programming language (http://www.Rproject.org) in conjunction with Bioconductor packages developed for the analysis of methylation arrays.

Preprocessing of Illumina Infinium 450K array data
The raw intensity data (IDAT files) were imported into R (3.3.1; http://cran.r-project.org/). Data quality was assessed using the minfi (v1.20.2) Bioconductor package [43]. From 485,512 HM450 probes, 67,120 were removed based on either (1) poor performance (mean detection p value of > 0.01, n = 14,056); (2) probes containing either a single nucleotide polymorphism (SNP) at the target CpG site or at the single nucleotide extension site (n = 16,307); (3) probes that map to multiple locations in the genome (n = 27,120), [44]; and (4) or to sex chromosomes (n = 9637). Samples were also evaluated using a modified version of the Houseman method [45,46] implemented in minfi, to estimate the cell type composition. The Wilcoxon signed-rank statistical test was used to compare the difference in cell type proportion between CP cases and controls. The analysis was completed before a cord blood reference panel was widely available, so cohorts used an adult whole blood reference [43] to estimate the proportion of B cells, CD8+ T-cells, CD4+ T-cells, granulocytes, NK cells and monocytes in each sample. The data was normalised using subset-quantile within array normalisation (SWAN) [47]. Covariates such as birth weight, birth order and postnatal age (in days) at which newborn screening cards were created were assessed as potential confounders.

Differential methylation analysis
Beta (β) values (proportion of the methylated signal over the total signal) were converted to M-values, the log 2 ratio of the intensities of the methylated signal versus the unmethylated signal. Differential methylation analysis was performed using remove unwanted variation (RUVm) [48], implemented in the missMethyl R package [49] taking into account the twin relationships. RUVm is a data-driven method of controlling for unwanted technical and biological variation in regression analyses using an empirically determined set of negative control probes assumed not to be associated with the biological factor of interest. p values were adjusted to control for the false discovery rate (FDR) using the Benjamini-Hochberg method [50]. Differentially methylated probes (DMPs) were considered significant if they fell within the FDR threshold of 0.1. We also investigated the top-ranked DMPs with an unadjusted p value less than 1 × 10 −4 [51].

Identification of differentially methylated regions
Differentially methylated regions (DMRs) were identified using the bumphunter package [43,52]. The cut-off value, which is a user-defined numeric value that determines the upper and lower bounds of the genomic profiles that will be used as candidate regions, was set to 0.02 and the number of permutations set to 1000.

Functional annotation and pathway analysis
Gene ontology and pathway analysis were performed using the gometh function from the missMethyl package [49], which appropriately takes into account the variable number of HM450 probes associated with each gene. Gene ontology enrichment was performed for the 1000 topranked DMP-associated genes. The KEGG option of the gometh function in missMethyl was used to provide further insights into relevant biological processes associated with the top-ranked DMPs.

Validation of differentially methylated regions
Site-specific validation was performed using the Sequenom MassArray EpiTYPER (Agena Biosciences). T7-tagged primers were designed for two regions (Additional file 1) using the Sequenom EpiDesigner package [53]. Forward primer sequences contained a 10 base 5′ tag (AGGAAGA-GAG) and reverse primers a 31 base 5′ tag (CAGTAA-TACGACTCACTATAGGGAGAAGGCT). In silico, assay prediction was performed using the BiocLite MassArray package. DNA used for validation was the same as that used for the HM450 analysis. Bisulphite treatment of genomic DNA was accomplished using the MethylEasy Xceed Kit (Human Genetic Signatures, North Ryde, Australia). One microliter of bisulphite-converted product was amplified in triplicate for each sample using the FastStart kit (Roche, Mannheim, Germany) in 15 μL of reagents with thermal cycling conditions as follows: 95°C for 10 min; 5 cycles of 95°C for 10 s, 60 for 30 s and 72°C for 2 min; 40 cycles of 95°C for 10 s, 62 for 30 s and 72°C for 90 s; and final extension at 72°C for 7 min. Raw data generated from the MassArray EpiTYPER was cleaned using a Microsoft Excel macro developed in-house [54,55]. The median value of triplicates was determined, and any replicates > 10% from the median were removed as previously described [54,55].

Within-twin pair analysis
To explore the top-ranked CpGs within each twin pair and compare them across twin pair groups, probes were ranked according to delta beta (Δβ, the difference in DNA methylation of the CP minus non-CP twin) within pairs, and the top-ranked 100 probes were compared across all 15 twin pairs. The genes corresponding to all probes with a Δβ value > 0.5 were then compared across the 15 twin pairs. Gene ontology analysis was performed on the top 1000 probes from each twin pair, and common ontologies between twin pairs were identified.

Subject characteristics
The study cohort consisted of 16 CP-discordant twin pairs (ten male and six female) for which pre-screening suggested a high probability of monozygosity (Table 1). All were tested for genetic zygosity using data from 65 SNPs from the Infinium arrays. The variability of SNPs for one twin pair (pair no. 9003) was substantially larger than the remaining samples and was therefore assigned as dizygotic (DZ). This pair was excluded from further analysis. Five subtypes of CP were reported: spastic diplegia (6), spastic quadriplegia (3), spastic hemiplegia (3), dyskinesia (2) and ataxia (1). The severity of CP ranged from mild (independently ambulant) to severe (wheelchair dependent), and the underlying neuropathology included white matter (11), grey matter brain injury (2) and both white and grey matter mixed injury (2). Three twin pairs were born at term (37-41 weeks), while all other twin pairs were born preterm (< 37 weeks).

Global DNA methylation profiles in CP-discordant monozygotic twins
Global DNA methylation (average β value across all probes) within twin pairs was compared using a pairwise Pearson correlation for the 418,392 probes remaining after filtering and quality control (see the 'Methods' section) for all 15 twin pairs. Within-pair methylation

Top-ranked CP-associated DMPs
Cleaned data was explored by principal component (PC) analysis which revealed few (6/54) significant correlations (p < 0.05, r < 0.6 shaded in Additional file 3) between the top six principal components of DNA methylation and nine technical (e.g. age at which Guthrie card was created) and biological (e.g. sex, subtype of CP) covariates. This suggested that none of the covariates tested was consistently associated with DNA methylation. In addition, none were associated with the largest principal component of variation within the dataset, PC1. Multi dimensional scaling (MDS) plots of the first three dimensions of the processed methylation data also showed that chip location and position on the 450K array were not found to affect methylation data (Additional file 4). Apart from the above covariates, it is known that celltype heterogeneity within the whole blood can confound epigenome-wide analyses. Therefore cell-type composition within CP-discordant pairs was evaluated. The levels of CD8 + T cells (CD8T) and CD4 + T cells (CD4T), B cells, natural-killer (NK) cells, monocytes and granulocytes were compared between the two groups (CP cases and normal co-twins). There was no statistically significant difference in the estimated cell-type proportions of CD8 + T cells, NK cells, B cells and monocytes (p > 0.05). However, the proportion of CD4 + T cells was lower (p = 0.002), and the proportion of granulocytes was found higher (p = 0.021) in CP cases relative to normal co-twins (Additional file 5).
To take into account potential sources of unwanted variation (such as cell-type composition), the genomewide analysis was performed using RUVm, which adjusts for biological and technical variation using a set of datadriven negative control probes [46,48,56,57]. This analysis did not identify any significant CP-associated DMPs after adjusting for multiple testing. Nevertheless, as this is an exploratory epigenome-wide study of CP, we focused on the characteristics of the top-ranked DMPs based on a nominal p-value cut-off of p < 1 × 10 −4 as used by others [51]. This resulted in a list of 33 topranked DMPs, corresponding to 25 genes ( Table 2), most of which showed a consistent direction in the majority of twin pairs (> 12/15). The average difference in methylation (Δβ = CP twin minus unaffected twin) ranged from + 0.6 to + 11.9% and from − 2.5 to − 12.4% (Table 2). Figure 1 shows the within-pair differences in methylation for the top ten DMPs.
The top-ranked probe cg00376816 (average Δβ = 11.6%, p = 4.57 × 10 −6 ) was located on chromosome 19, within the gene body of the HNRNPL gene encoding the heterogeneous nuclear ribonucleoprotein L. Others included cg04242728 (in the 5′ end of the TBC1D24 gene, ranked 3) and cg19607845 (in the gene body of FBXO9 gene, ranked 13). Probes located near the gene body of immune and inflammatory genes, such as Ras association domain family member 5 (RASSF5), major histocompatibility complex DM alpha-chain (HLA-DMA), CD3D and kalirin (KALRN) genes, were also among the top-ranked 33 DMPs.
To identify enriched biological processes or molecular functions, we performed gene ontology analysis on genes associated with the top-ranked 1000 DMPs. The top 20 gene ontologies ranked by nominal p-value were 'regulation of immune response' , 'lymphocyte activation' , 'differentiation and aggregation and T cell activation' with the top two ontologies related to cell-cell adhesion processes (Table 3; Additional file 6). Enriched disease pathways, as reported by KEGG analysis, included MAPK signalling (19 associated genes from 245 genes in the KEGG pathway list; p value 3.6 × 10 −10 ), cytokine-cytokine receptor interaction (13 associated genes from 240 genes; p value 1.3 × 10 −08 ) and Ras signalling (15 associated genes from 218 genes; p value: 9.8 × 10 −08 ).
We also identified DMRs associated with CP [52] ( Table 4). The top-ranked DMR, spanning 434 bp and with a p value of 5.6 × 10 −4 , was located on chromosome 6 ( Fig. 2). This DMR spans 12 probes (average Δβ = 3.7%) within the coding region of the LTA gene, approximatelỹ 800 bp downstream of the transcription start site (TSS). LTA codes for the lymphotoxin-alpha protein otherwise known as tumour necrosis factor beta (TNF-β). Other top DMRs include those within LTBP1, CD300, CHST11 and LIME1. Gene ontology analysis of the top-ranked DMR-associated genes revealed similar findings to topranked DMPs (Table 5; Additional file 7). We found an over-representation of inflammatory signalling pathways, also similar to that of the top-ranked DMPs. The top pathways included TNF and TGF-beta signalling and cytokine-cytokine receptor interaction. The nuclear factor kappa-light-chain-enhancer of activated B cell (NF-κB) signalling pathway was also enriched.

Validation of DMRs
LIME1 and LTA DMRs were selected for validation as they were highly ranked, had large, consistent effect sizes across all pairs and were biologically relevant to CP. Three CpGs from the HM450 platform contained within three CpG units on the MassArray Epityper (consisting of seven CpG sites in total) were tested for the LIME1 DMR, and four CpGs contained within three units on the MassArray Epityper (four CpG sites in total) were tested for the LTA DMR, both in regions being approximately 200 base pairs upstream of the transcriptional start site (TSS200) and likely to be in gene promoters. Scatter plots were generated to assess the validity of the data for the three LIME1 and LTA probes within a 250base-pair region across both the HM450 and the EpiTYPER platforms (Additional file 8). Pearson's correlation coefficients were determined, and the significance of the correlation was assessed for each probe-CpG unit comparison. Five out of six probes (cg21201401, cg06653796, cg14597739, cg11586857, cg21999229) had a positive correlation between the two platforms. Among these, one probe out of three for LIME1 and two of three for LTA had moderate correlations (r > 0.5). All moderate correlations were also significant with p < 0.05 (r = 0.88, p = 1.

Differential methylation analysis within individual twin pairs
Since CP is a highly heterogeneous disorder [58], it is possible that a subset of disease-associated DNA methylation patterns may be specific to each proband. To test this [59], we determined the top CP-associated CpGs for each twin pair ranked by absolute differences in DNA methylation (Δβ > 0.5) and looked at significant gene ontologies common to multiple pairs (Additional file 9). Gene ontologies corresponding to cell adhesion were found in 5/15 twin pairs (Additional file 10). Similarly, common CpGs with a within-pair methylation difference > 50% were found in multiple twin pairs (Table 6) in genes such as BICD2, HLA-DPB2, RPTOR and PIK3CG (Additional file 11), involved in neuronal cell migration, muscular atrophy or muscle contraction pathways and immune response and inflammatory pathways, respectively [60][61][62][63]. Notably, two different CpG sites within twin pairs 4 and 8 corresponding to the WWTR1 gene had an absolute methylation difference of greater than 50% with the same direction of effect.
Affected CpG sites were located near the 5′ end of the gene in both pairs (Table 6).

Discussion
This exploratory study represents an initial step towards investigating potential CP-associated epigenetic differences, with the longer-term aim of identifying predictive biomarkers with clinical utility. We identified DNA methylation differences in dried blood spots from 15 CP-discordant MZ twin pairs and found differential methylation at several gene loci associated with hypoxia signalling, inflammation and cell adhesion. These pathways had been previously linked to CP, consistent with part of our hypothesis. Pairwise global DNA methylation difference between CP and non-CP members of each pair measured for Fig. 1 DNA methylation differences (CP minus normal twin) for the top ten CP-specific DMPs. The direction of methylation difference for each probe was consistently in the same direction (hyper-or hypomethylated) across > 12/15 twin pairs. The dotted line represents the average methylation difference, and methylation difference is represented in absolute value (effect size of 0.1 is 10% methylation difference) each twin group and comparison of the size of DNA methylation difference made between groups allowed for an assessment of how variable the differences in methylation may be for different cases of CP. Our results are consistent with previous studies (e.g. [59]) that have indicated that neurodevelopmental disorders such as autism spectrum disorder are not associated with systemic within-twin pair differences in global DNA methylation.
We tested site-specific DNA methylation patterns across the genome for their association with CP, taking CP status within discordant MZ pairs into account. Although no probe reached an adjusted statistical significance at FDR < 0.1, the top-ranked DMPs, at a nominal cut-off of p < 1 × 10 −4 , were enriched for the cellular processes of inflammation, cell adhesion and immune response. These showed the direction of effect across most or all discordant twin pairs. This was in accordance with EWAS of other neurodevelopmental disorders including twins discordant for autism spectrum disorders [59], asthma [51], depression [64] and aggressive behaviour [65].  Average within-pair DNA methylation differences of up to 12.5% were observed, comparable to previous findings in neurodevelopmental disorders, ranging from 1.5 to 12%. Furthermore, we validated five of the six CpG probes, with positive correlations across platforms. Three of these five probes (one from LIME1 and two from LTA) showed a strong and significant cross-platform correlation, indicating the validity of methylation values between platforms. Genes associated with top-ranked DMPs and DMRs were enriched for similar ontologies and pathways, namely immune response, lymphocyte-mediated immunity, interferongamma production and regulation of immune response.
Our study revealed top-ranked DMRs associated with genes that play a role in inflammation, such as LTA/TNFβ and LIME1, supporting part of our hypothesis that inflammation plays a key role in CP aetiology. Genetic variants of LTA have been implicated in multiple studies as being associated with risk for CP [66][67][68]. LTA plays an important role in inflammation and brain development, mediating preterm birth and white matter brain injury [69]. It has been implicated that inflammation [70] and increased levels of its isoform TNF-α were found in children with CP compared to healthy controls [71]. LIME1 gene links T and B cell signalling to the activation of tyrosine and MAP kinases [72]. Other top DMR-associated genes such as LTBP1 and CD300 are also known mediators of inflammatory pathways such as ERK signalling pathway, interleukin-3,-5 pathway, B cell receptor signalling pathway and other chemokine signalling [72][73][74]. Genes associated with top-ranked DMPs also showed involvement in other key inflammatory pathways such as Ras signalling (WDFY4), MAP kinase signalling (CD3D) and interleukin-3,-5 signalling (KALRN). Analysis within individual twin pairs also revealed associations to the WWTR1 gene, which is involved in the activation of TGF-β signalling pathway, an inflammatory pathway that regulates neural survival and death [75,76]. Gene ontology analysis of top-ranked DMPs showed enrichment of genes involved in regulation of immune response pathways such as those involved in signalling or T-cell activation. It is noteworthy that intrauterine infection is a known risk factor for CP and that many inflammatory cytokines have been shown to be critical to the risks associated with CP [38,41,60]. Perinatal brain injury can be induced by a range of insults such as hypoxic-ischaemic injury or infection [77]. An in utero infection such as chorioamnionitis may trigger an innate immune response, resulting in elevated cytokine levels. Cytokines in the fetal blood may contribute to a systemic fetal inflammatory response with eventual penetration across the blood-brain barrier resulting in an inflammatory cascade in the fetal brain [78]. Brain injury induced by neonatal hypoxia-ischaemia also involves key components of inflammation such as immune cells, chemokines, cytokines and cell adhesion molecules [79]. Therefore, we suggest that inflammation may play a role in perinatal brain damage associated with CP.
We also observed an enrichment of CP-associated DMPs and DMRs in gene ontologies associated with cell adhesion, and this was also observed in individual twin pairs. Aberrant expression of cell adhesion molecules has been reported in muscle biopsies of both children and adults with CP [39]. Previous whole exome and whole genome sequencing studies have also illustrated the potential role of cell adhesion in CP by identifying genetic variants in novel candidate genes which function as neural adhesion molecules essential for neurite outgrowth and axon guidance [9,12]. The NF-κB transcription factor signalling pathway was common in both DMPs and DMRs.
Our results are consistent with a previous gene expression study in newborn blood spot samples from children with CP [80], which identified up-regulation of inflammatory pathways in preterm children who later developed the disorder. Other similarities between the two studies include variation in genes involved in T-cell and B-cell receptor signalling pathways and cytokine-cytokine receptor interaction, all of which were shown to have a dysregulation in CP cases [80]. However, we found no evidence for an association with increased thyroid function in preterm-born CP cases as hypothesised and reported previously [80].
Another top-ranked DMP was HNRNPL, which likely plays a role in response to hypoxia via regulation of the vascular endothelial growth factor (VEGF) gene [81]. Hypoxia is known to hinder normal development and maturation of the brain and can cause white matter injury in preterm born infants [82] resulting in CP [83]. This finding supports our hypothesis that epigenetic alterations in genes involved in hypoxic pathways play a role in the aetiology of CP. Two high-ranking DMPs lie within the TBC1D24 and FBOX9 genes respectively, and both have previously been associated with epilepsy [84,85]. In 29% of CP cases in Victoria, Australia, epilepsy is comorbid with CP [5]. These results may suggest a potentially shared aetiology between epilepsy and CP [86,87].
Our findings agree with those of others showing a link between early life DNA methylation state and neurodevelopmental and cognitive outcomes [29,88], which would allow for early diagnosis and facilitate timely intervention. Currently, MRI scans, assessment tests such as the General Movements Assessment and interventions such as environmental enrichment, early developmental, early motor and physiotherapy interventions are used to inform strategies for early intervention in high-risk groups, such as preterm born children, [89][90][91][92][93].
Given that CP is a highly heterogeneous condition, this study highlights the importance of using epigenetic biomarkers to distinguish and detect underlying pathways across the disorder. For individuals, such an epigenetic state at birth could be used to estimate risk for subsequent development of overt CP.
The strength of using MZ twins is that they are matched for parental age, age, sex, season of birth and genetic factors. Although twins have a higher risk of CP than singletons, the causative mechanisms, such as thrombosis, and infection in the mother, the placenta or the umbilical cords are likely to be similar [7]. Studying twins discordant for CP allows genetic and environmental components to be partitioned from each other and provides a unique opportunity to evaluate the importance of non-shared environmental factors such as umbilical cord or placental function during early development in isolation. It is possible that only twin of a pair may develop an infection or inflammation of the umbilical cord or placenta ( [13,19,20,94]). Such non-shared environmental factors are known to influence the development of individual members of a twin pair. This pilot study also highlights the importance of analysing DNA methylation in dried blood spots, which are collected at birth and stored by many countries and the potential for developing future predictive diagnostic tests [55,95,96].
Although each tissue has a subset of CpGs whose DNA methylation patterns are tissue-specific, DNA methylation changes that are concordant between the blood and brain have been detected in previous studies [32,97,98]. One example is where they identified parallel changes in DNA methylation between the brain and blood in 30% of the genes implicated in Parkinson's disease [98]. It was also shown that a DNA methylation module exists in key ageing-related regulatory genes both in the brain and blood [99]. In addition, animal studies have reported that an early environment resulting in a brain disorder can alter DNA methylation in the same gene across the brain and peripheral tissues [35].
There are some limitations to this study. While similar in sample size to many comparable twin studies of brain-related disorders [58,59,64,100], we acknowledge that larger sample sizes of 25 twin pairs or more are preferable to detect a mean effect size of at least 8% methylation (FDR = 0.05) [101]. As CP is a heterogeneous condition, the small sample size of our cohort, with a lack of CP concordant and healthy twin pairs for comparison, also limits our capability to understand the biological mechanisms of brain injury that may be specifically associated with CP subtype. The use of peripheral tissue and blood also limits our capability to pinpoint the mechanism of CP. Despite the fact that the brain and blood arise from separate cell lineages, and are thought to be epigenetically distinct, many epigenetic studies are often conducted in the blood due to ease of availability [102]. Previous investigations of methylomic variation across the blood and brain tissue from different regions of the brain have found distinct differences in gene expression and DNA methylation patterns [26,[103][104][105]. Studies have also shown the inconsistencies in DNA methylation markers from the blood in predicting brain DNA methylation status [27,106]. However, evidence from animal studies have also shown that blood DNA methylation patterns may in fact reflect patterns in the brain in a subset of genes [35,107], suggesting that peripheral epigenetic marks may reflect disease mechanisms in some cases. Examples where methylation levels correlate between blood and brain have been reported in Parkinson's disease, depression, schizophrenia, bipolar disorder and autism [107,108]. The blood is also particularly useful in investigating disease biomarkers [98] and is an important peripheral tissue to consider for neurological disorders, as it is easily accessible to assist in diagnosis. Another limitation is that our data apply to twins only, and we cannot yet generalise our findings more broadly, as there is evidence that risk factors and associated mechanisms leading to CP may be different in twins compared to singletons [109]. To overcome this, our analysis will be repeated in further sets of twins and singletons. Only with this information can we then start to put together risk models for predicting CP at the time of birth. This approach will provide a unique opportunity to identify a biomarker to predict neurodevelopmental outcomes such as CP.

Conclusion
This study provides the first evidence that environmentmediated differential methylation in genes involved in known processes such as hypoxia and inflammation, and perhaps processes such as cell adhesion, may contribute to the development of CP. Our data also pave the way for larger studies to use DNA methylation data in risk models to help predict CP before the onset of overt symptoms and therefore provide a chance for timely ameliorative interventions.