DNA methylation changes at infertility genes in newborn twins conceived by in vitro fertilisation

The association of in vitro fertilisation (IVF) and DNA methylation has been studied predominantly at regulatory regions of imprinted genes and at just thousands of the ~28 million CpG sites in the human genome. We investigated the links between IVF and DNA methylation patterns in whole cord blood cells (n = 98) and cord blood mononuclear cells (n = 82) from newborn twins using genome-wide methylated DNA immunoprecipitation coupled with deep sequencing. At a false discovery rate (FDR) of 5%, we identified one significant whole blood DNA methylation change linked to conception via IVF, which was located ~3 kb upstream of TNP1, a gene previously linked to male infertility. The 46 most strongly associated signals (FDR of 25%) included a second region in a gene also previously linked to infertility, C9orf3, suggesting that our findings may in part capture the effect of parental subfertility. Using twin modelling, we observed that individual-specific environmental factors appear to be the main overall contributors of methylation variability at the FDR 25% IVF-associated differentially methylated regions, although evidence for methylation heritability was also obtained at several of these regions. We replicated previous findings of differential methylation associated with IVF at the H19/IGF2 region in cord blood mononuclear cells, and we validated the signal at C9orf3 in monozygotic twins. We also explored the impact of intracytoplasmic sperm injection on the FDR 25% signals for potential effects specific to male or female infertility factors. To our knowledge, this is the most comprehensive study of DNA methylation profiles at birth and IVF conception to date, and our results show evidence for epigenetic modifications that may in part reflect parental subfertility.


Background
As the frequency of in vitro fertilisation (IVF) treatment increases worldwide, much research effort has focused on exploring both short-and long-term health outcomes associated with conception via IVF, with contradictory results. A number of studies have observed associations with adverse perinatal and obstetric outcomes, including low birth weight, preterm birth, perinatal mortality, congenital malformations, placental complications, and increased frequency of imprinting disorders such as Angelman syndrome and Beckwith-Wiedemann syndrome [1][2][3][4]. On the other hand, parallel efforts have reported that these associations are not attributed to IVF treatment itself, but rather to multiple pregnancy or parental subfertility, both common factors in IVF births [5,6]. Further research is required to identify potential factors associated with conception via IVF, including not only health outcomes but also biological consequences such as epigenetic modifications.
Given that birth weight and imprinting disorders are controlled at least in part by epigenetic factors [7,8], IVF may have an influence on epigenetic profiles, potentially resulting in changes that persist well after birth and over the life course. Epigenetic mechanisms are considered possible mediators of the developmental origins of health and disease [9]; therefore, an assessment of the influence of IVF on DNA methylation profiles may give some insights into mechanisms underlying potential related health outcomes. Establishment of DNA methylation profiles in the germ line and embryo takes place early in development [10]. Theoretically, this epigenetic reprogramming could therefore be influenced by IVFrelated interventions that occur very early, prior to blastocyst implantation. Indeed, induction of ovulation, embryo culturing, and cryopreservation, among others, have all been linked to specific alterations in DNA methylation in mice, although results are somewhat inconsistent [11][12][13].
Most studies in humans comparing naturally and IVFconceived newborns have interrogated DNA methylation alterations targeting almost exclusively imprinted differentially methylated regions (DMRs). These studies have reported increased epigenetic variability at the KvDMR1, PEG1, and H19 DMRs in umbilical cord blood [14], hypomethylation of the H19 and MEST DMRs in placenta [15], and hypomethylation of the H19 DMR in buccal epithelium [16] in individuals conceived by IVF. Highthroughput approaches using bead array technology have also interrogated DNA methylation in IVF in a genome-wide manner. Katari et al. [17] reported differential methylation at 78 genes in cord blood and 40 in placenta with at least two differentially methylated CpG sites (P ≤ 0.08) when looking across the promoters of 736 genes (GoldenGate Array, Illumina) in ten cases and 13 controls. A more extensive study using the promoterenriched Illumina Infinium HumanMethylation27 bead array in cord blood samples from ten IVF cases and eight controls reported a total of 24 genes with at least two differentially methylated CpG sites (P < 0.05) [18]. More recently, a study used the genome-wide Illumina Infinium HumanMethylation450 bead array in samples from 38 IVF-conceived newborns followed by fresh embryo transfer, 38 IVF-conceived followed by cryopreserved embryo transfer, 18 born to subfertile parents after conception by intrauterine insemination, and 43 controls born to fertile parents [19]. This platform interrogates CpG sites across the whole genome, although with a limited coverage since it targets gene-centric annotations [20]. The authors identified differential methylation at multiple sites, including metastable epialleles.
Here, we interrogated evidence for differential methylation between IVF and non-IVF newborn twins in a more comprehensive manner by conducting epigenomewide association scans (EWAS) [21] using methylated DNA immunoprecipitation followed by deep sequencing (MeDIP-seq) [22] genome-wide in samples from cord blood, and its mononuclear fraction, collected at birth from IVF and non-IVF twins. The use of twins in this study allowed the partition of the observed variance in DNA methylation into genetic and environmental factors. The approach also avoids potential spurious associations due to an imbalanced number of multiple and single pregnancies between conception method groups.

Subjects and sample collection
The study included 47 IVF and 60 non-IVF newborn twins (from 54 twin pairs) from the Peri/postnatal Epigenetic Twins Study (PETS), Melbourne, Australia. Recruitment and full study procedure have been described previously [23,24]. Cord blood was collected at birth and used to process mononuclear cells by Ficoll gradient centrifugation as described previously [25]. Whole blood cells (WBCs) from cord blood were available for a total of 98 twins (40 IVF and 58 non-IVF) and cord blood mononuclear cells (CBMCs) for a total of 82 twins (35 IVF and 47 non-IVF). Maternal age and method of conception were determined via questionnaire at recruitment (18-20 weeks gestation). Twins of mothers who said yes to IVF or intracytoplasmic sperm injection (ICSI) treatment were classified as IVF regardless of the use of ovulation induction medication or other fertility treatments. Maternal smoking status was collected via questionnaire on recruitment and at 24 and 36 weeks of pregnancy. Birth weight was collected during the immediate neonatal period. Zygosity and chorionicity were determined by physical examination of the inter-placental membranes at birth, and by genetic test when required, as described previously [23,24]. Pregnancy complications were recorded and are shown in Additional file 1: Table S1.

DNA methylation profiling
MeDIP-seq was performed at BGI-Shenzhen, Shenzhen, China. Extracted DNA was fragmented using a Covaris sonication system and sequencing libraries were prepared from 5 μg fragmented genomic DNA. End repair, <A > base addition and adaptor ligation steps were performed using Illumina's Single-End DNA Sample Prep kit. Adaptor-ligated DNA was immunoprecipitated by anti-5mC using a commercial antibody (Diagenode) and MeDIP products were validated by quantitative PCR. MeDIP DNA was purified with ZYMO DNA Clean & Concentrator-5 columns and amplified using adaptormediated PCR. DNA fragments between 200 and 500 bp in size were gel-excised, and the amplification quality and quantity were evaluated by Agilent BioAnalyzer analysis. The libraries were subjected to highly parallel 50-bp single-end sequencing on the Illumina GAII platform. All sequencing data passed initial quality checks for base composition (no exclusions) using FASTQC v0.10.0. For each individual,~30 million reads were generated and mapped onto hg19 using BWA. After removing duplicates, we filtered data using quality score Q10. We quantified methylation levels using MEDIPS [26], producing the mean relative methylation score (RPM) in 500-bp bins (overlap of 250 bp) across the genome. Altogether, there were 11,524,145 windows and these were used for the analyses. Bins with RPM values of zero in more than 50% of the samples were excluded, resulting in 9,592,803 (WBC) and 9,285,089 (CBMC) bins used in downstream analyses.

Epigenome-wide IVF-DMR analyses
Normalised (N(0,1)) methylation scores in each genomic bin were regressed using a linear mixed-effects model to account for twin structure (lme4 package [27] in R [28]). Tissue type, birth weight, sex, maternal smoking, 260/ 280 ratio, DNA concentration, and the loadings of the first five principal components were used as covariates and included as fixed effects in the model. Family and zygosity were included as random effects in the linear mixed model. The principal components were included to account for unknown sources of variation, such as cell heterogeneity. Correction for multiple testing was performed by a Benjamini-Hochberg false discovery rate (FDR) calculation.

Variance decomposition of WBC IVF-DMRs
The contribution of additive genetic (A), common environmental (C), and unique environmental (E) factors to DNA methylation was estimated using the ACE model based on the classic twin design [29]. The model was fitted using the OpenMX statistical package [30]. RPM values without adjustment for covariates were used to estimate the ACE proportions.

Statistical analysis
Pairwise correlations and principal components analysis were performed using RPM values across all bins with values > 0 in at least 50% of the samples. Hierarchical clustering was performed using Euclidean distance as a measure of dissimilarity and average linkage clustering.

Validation analysis
Genomic DNA (500 ng) was bisulphite converted using the MethylEasy Exceed Rapid Bisulphite Modification Kit (Human Genetic Signatures, North Ryde, NSW, Australia). Primers to target the regions in TNP1 and C9orf3 were designed using the EpiDesigner tool (Sequenom Inc., Herston, QLD, Australia). The H19 CTCF6 region was the same used in a previous study [25]. Primers, genomic coordinates, and PCR conditions are shown in Additional file 1: Table S2. Methylation levels were determined by EpiTYPER on the MassARRAY System (Sequenom Inc., Herston, QLD, Australia). Statistical analysis considered the average of two to three technical replicates and were performed using data on single CpG sites.

Genome-wide methylation profiles in twins
We profiled DNA methylation levels from a total of 107 newborn twins (47 conceived via IVF and 60 conceived in vivo) in WBCs and CBMCs. Details of any fertility treatment used and demographic characteristics that represent potential confounders of DNA methylation levels at birth, such as sex, birth weight, maternal age, and maternal smoking status, are shown in Table 1. We first explored the genome-wide patterns of DNA methylation variability in the dataset. Principal component analysis was used to identify factors that were significantly associated with genome-wide variability in DNA methylation profiles. The first five principal components in the dataset, which explained~13% of the total variance in DNA methylation, were at least nominally associated (P < 0.05) with sample type (WBCs versus CMBCs), birth weight, maternal smoking, and conception method (Fig. 1a).
We next estimated the within twin-pair correlation patterns in methylation profiles of twin pairs available in both datasets using Pearson's correlation. In concordance with previous studies [7], we observed higher median correlation within monozygotic (MZ) twin pairs compared to dizygotic (DZ) twin pairs (Fig. 1b). Previous studies have shown that twin chorionicity can have an effect on within-pair DNA methylation differences, but not with consistent direction of effect across tissues [7,25,31]. In our study, we did not observe significant chorionicity-related methylation differences (Fig. 1b), but the number of MZ twins within chorionicity categories was relatively low (n = 8 monochorionic and n = 5 dichorionic pairs).
Interestingly, the method of conception showed methylation profile differences within MZ twin pairs. MZ IVF twins had higher median correlation compared to MZ non-IVF twins in WBCs, but the opposite trend was observed in CBMCs, and in both cases the MZ IVF sample was very small (n = 3).

IVF-DMRs in CBMCs and WBCs
In order to identify tissue-independent and tissue-specific IVF-associated DMRs, we compared DNA methylation profiles in WBCs and CBMCs in relation to method of conception adjusting for birth weight, sex, maternal smoking, and the first five principal components, which partly capture cell heterogeneity. Epigenome-wide analyses of DNA methylation in relation to method of conception did not identify genome-wide significant signals in the CBMCs subset or in the combined CBMC and WBC datasets, after correction for multiple testing. In WBCs alone, one significant DMR was observed at a FDR of 5% (Fig. 2). This was located~3 kb upstream of TNP1 (chr2:217,726,751-217,727,250), which encodes a transition nuclear protein that replaces histones and is subsequently replaced by protamines during spermiogenesis. A deletion in the promoter region of this gene, which reduces its expression, has been reported in infertile men [32]. Methylation upstream of TNP1 might have an impact on its expression. In mice, methylation changes during spermatogenesis have been observed at TNP1, which suggests a role of methylation in the regulation of this gene [33]. To explore the biological characteristics of the top-ranked results in the IVF epigenome-wide analyses we selected a more liberal threshold of FDR 25%, at which 46 IVF-DMRs were identified (Table 2). Interestingly, the third-ranked DMR genome-wide (Additional file 1: Figure S1) was located in the first intron of C9orf3 (chr9:97,504,001-97,504,500), which has been associated with polycystic ovary syndrome in women [34] and development of erectile dysfunction after radiotherapy for prostate cancer in men [35]. Another signal within this list was located in intron 1 of STOX2 (chr4:184,814,001-184,814,500), whose reduced expression has been implicated in pre-eclampsia [36].
Since adverse perinatal outcomes may be associated with maternal age, we further adjusted for this covariate and observed that the 46 FDR 25% WBC IVF-DMRs remained significant ( Table 2). The non-IVF group included a small number (n = 4) of newborns conceived with other types of fertility treatments not equivalent to IVF, such as gamete intrafallopian transfer (GIFT) and ovarian stimulation. We re-analysed the 46 FDR 25% WBC IVF-DMRs excluding GIFT (n = 2) and non-IVF ovarian stimulation (n = 2) controls and observed that conclusions remained unchanged (Additional file 1: Table S3).
Hierarchical clustering using DNA methylation levels at these 46 FDR 25% DMRs alone grouped twins by method of conception, assigning 38 out of 40 IVF twins and 57 out of 58 non-IVF twins to the correct group (Fig. 3). We also explored these signals with respect to functional annotations. A total of ten FDR 25% WBC IVF-DMRs overlapped CpG sites previously shown to be dynamic during development [37], 20 overlapped DNase I hypersensitivity sites (wgEncodeRegDnaseClusteredV3) [38], one overlapped a CpG island (cpgIslandEx) [39],  and none overlapped with candidate metastable epialleles [40] (Fig. 3). Cell type-specific DNA methylation can impact the profiles observed in a population of cells, such as in a whole blood sample, and we therefore accounted for blood cell type heterogeneity using a twofold approach. First, we performed principal component analysis on the methylation levels of the entire set of WBC samples, and our main EWAS analyses above are corrected for the first five principal components, which likely capture variation attributed to technical and biological factors, potentially including cell heterogeneity. To assess whether the first five principal components capture cell heterogeneity, blood cell subtype counts were obtained through automatic differential counting for a subset of the WBC samples (n = 54 twins, 22 IVF, and 32 non-IVF) and these were compared against the distributions of the first five principal components. The proportion of neutrophils, eosinophils, and lymphocytes were associated (P < 0.05) with the loadings of the second, third, and fourth principal  components, respectively (Fig. 1a). Therefore, since the EWAS model used in this study took into account the loadings of the first five principal components, these analyses already take into account the influence of cell heterogeneity to a certain extent. Second, we re-analysed the 46 FDR 25% WBC IVF-DMRs in the subset of 54 WBC samples with available cell counts, adjusting for the proportion of neutrophils, eosinophils, monocytes, and lymphocytes. We also performed analyses adjusting for the loadings of the first five principal components within this dataset alone. Most results were concordant when comparing across all models (Additional file 1: Table S4) and only five out the 46 FDR 25% WBC IVF-DMRs were not significant (P > 0.05) after adjusting for cell proportions (chr8:120,972,0 01-120,972,500, chr7:2,487,251-2,487,750, chr18:23,69 5,001-23,695,500, chr12:4,310,251-4,310,750, and chr1 4:104,067,251-104,067,750).

Variance decomposition of WBC IVF-DMRs
Given that epigenetic changes were potentially affecting infertility genes, we wanted to investigate if the findings may capture a genetic signature affecting DNA methylation that could be transmitted to offspring. We applied twin variance decomposition analyses to partition the total epigenetic variance into additive genetic (A) and common (C) and unique (E) environmental components (ACE) [29]. The ACE model was used to determine the contribution of genetics, shared intrauterine environment due to shared maternal influences, and non-shared (twin-specific) or stochastic factors to epigenetic variation. The mean contribution of additive genetic effects (narrow-sense heritability) to DNA methylation across the genome in different tissues from newborns has been previously estimated to be between 0.05 and 0.12 [7]. Here we estimated the average genome-wide narrow-sense heritability for DNA methylation in WBCs at 0.06. At the 46 FDR 25% WBC IVF-DMRs, the major contributors to DNA methylation variation were non-shared or stochastic events (Fig. 3). However, several FDR 25% IVF-DMRs had evidence for heritability (A > 0.4), suggestive of genetic effects underlying specific IVF-associated DNA methylation changes. These included an intronic region in DEC1 (chr9:118,148,751-118,149,500), a region 33 kb away from XYLT1 (chr16:17,161,751-17,162,250), and an intergenic region in chromosome 12 (chr12:4,310,251-4,310,750). When looking at the two DMRs associated with infertility genes, DNA methylation variation showed no evidence for genetic effects (A = 0) near TNP1, while heritability at the DMR in C9orf3 was estimated at 0.25.

Effects of IVF on imprinting
Previous studies have explored DNA methylation patterns in IVF births specifically at imprinting control regions (ICRs). We therefore assessed whether there was an enrichment of differential methylation effects at 34 known ICRs [41] in our genome-wide results, but no enrichment was observed (P > 0.05). However, when we explored individual signals at candidate IVF-DMRs we were able to replicate one previously reported ICR IVFassociated DMR. Concordantly with previous IVF methylation studies in placental tissue [15] and buccal epithelium [16], we observed hypomethylation in IVF twins at the sixth CTCF binding site within the H19/IGF2 (H19 CTCF6) DMR (Additional file 1: Figure S2). This association was observed in CBMCs (P = 0.01), but not in WBCs.

Effects of ICSI
ICSI is a technique in IVF used to treat couples with male-factor infertility [42]. In contrast to conventional IVF where fertilisation occurs by placing spermatozoa near an egg, ICSI consists of the direct injection of a selected single sperm cell into the egg. This manipulation may introduce additional risk factors [43]. To assess the effect of ICSI on the 46 FDR 25% WBC IVF-DMRs we adjusted for the use of this technology and also compared the ICSI and the conventional IVF groups separately against the non-IVF group. After adjustment for ICSI, the association weakened at several FDR 25% IVF-DMRs ( Table 3), suggesting that ICSI or paternal infertility might have a role in these methylation changes. One FDR 25% IVF-DMR signal (chr1:85,522,251-85,522,750) appeared stronger after adjustment, suggesting either a female infertility effect or that ICSI prevents or corrects a methylation change that occurs in conventional IVF. This DMR was located upstream of WDR63, a gene mainly expressed in testis, fallopian tube, and adrenal gland [44].

Validation of IVF-DMRs
We pursued validation of the differential methylation signals at the top associated DMR (located~3 kb upstream of TNP1) and at the third-ranked DMR (located in C9orf3), both in or near genes previously linked to infertility. Altogether, four CpG sites were targeted for validation using Sequenom's EpiTYPER technology.
For the DMR in C9orf3, we were able to target two CpG sites within the most-associated 500-bp bin in this locus (Additional file 1: Figure S3). We assayed methylation levels in 36 MZ twins included in the discovery EWAS and observed significantly higher methylation in the IVF group, concordant with the MeDIP-seq analysis, at both tested CpG sites in the C9orf3 locus (P = 0.02 and 0.03, respectively), therefore validating this signal using a different methylation profiling approach (Additional file 1: Figure S4).
For the TNP1 DMR we were unable to target CpGs within the most associated 500-bp bin, and we therefore selected two of the closest CpG sites contained within the second most associated DMR in that locus (Additional file 1: Figure S3). Within the sample of 36 MZ twins we also observed higher methylation in the IVF group, consistent with the MeDIP-seq signal, with effects close to nominal significance (P = 0.08; Additional file 1: Figure S4). However, correlation between the MeDIP-seq signal at the most-associated DMR in TNP1 and the EpiTYPER methylation values was, as expected, relatively low as we were unable to target CpG sites within this most-associated DMR (correlation of 0.18 and 0 at the two tested CpG sites). We profiled additional samples from DZ twin pairs but did not obtain validation of the signal. We also considered the effect of ICSI compared to conventional IVF in MZ twins in the validation dataset. We observed significantly higher methylation in the ICSI group at the first CpG of the targeted region near TNP1 and at the first CpG site of C9orf3 (Additional file 1: Figure S5).
Lastly, we also compared methylation in relation to conception method at the H19 CTCF6 DMR in a reduced subset of CBMCs samples (n = 42 twins) using EpiTYPER. When comparing IVF to non-IVF twins (Additional file 1: Figure S6) we observed a difference with the same direction of effect as in the MeDIP-seq analysis, although not significant (P = 0.19). Interestingly, when comparing naturally conceived twins to twins that were conceived with any type of medical help (Additional file 1: Figure S6), i.e. not exclusively IVF, the difference reached nominal significance (P = 0.04), suggesting that differential methylation at this region is associated with parental subfertility rather than IVF conception.

Discussion
Since IVF procedures are carried out during an important period of epigenetic reprogramming in early development, we hypothesised that IVF may induce epigenetic differences that persist to birth. We were able to identify significant and suggestive DMRs related to IVF conception (IVF-DMRs) in WBCs, although our results suggest that at least some of these changes may be linked to parental subfertility, which is confounded with IVF treatment. The observation that IVF-DMRs were identified close to genes implicated in fertility and reproduction suggests that a genetic signature influencing DNA methylation could be transmitted from parent to offspring. To assess this further, we estimated the heritability of the IVF-DMRs. We observed that the IVF-DMR located in C9orf3, a gene associated with polycystic ovary syndrome, was estimated to have a heritability at 25% and eight other FDR 25% WBC IVF-DMRs showed heritability greater than this (Fig. 3).
Epigenetic states of metastable epialleles in mammals are mitotically inherited after establishment in early development, therefore shared across tissues, and can cause expression variability within isogenic individuals [45]. A study in humans looking for systematic interindividual variation in DNA methylation across tissues from two different lineages identified 109 candidate metastable epialleles [40]. Nutritional conditions during conception have been shown to be important to the establishment of epigenetic states at some of these metastable epialleles [46]. If an influence of IVF on the epigenetic marks of these alleles exists, it could potentially cause long lasting effects. A previous study, which included newborns from single and multiple pregnancies, identified DNA methylation differences in IVF conception at candidate metastable epialleles, although at different epialleles to those affected by maternal nutritional factors [19]. In our study, none of the 109 candidate metastable epialleles overlapped with the 46 FDR 25% WBC IVF-DMRs. This discrepancy could be attributed to differences between single and multiple pregnancies or to low power to detect such changes.
Our results also showed that IVF-DMRs, including hypomethylation of the regulatory region of H19, were generally not shared between WBCs and CBMCs. This observation suggests that the epigenetic differences reported here likely did not appear during early development or that these effects are not fixed and can revert in a cell type-specific manner. CBMCs, in contrast to WBCs, lack the granulocyte fraction, which is the predominant group of cells in the blood. Thus, the IVF-DMRs may be granulocyte-specific or at least in part influenced by this group of cells.
To date, there has been mixed evidence on the effect of IVF at imprinted genes and their regulatory regions. Some studies have reported DNA methylation changes or increased variability at these imprinted regions [14][15][16], while others have reported no associated changes [47,48]. We observed that there is not an overall destabilisation of methylation patterns in ICRs, but specific DMRs, such as the H19 DMR, can show a weak but nominally significant association with the method of conception. Previous studies have reported similar observations, that is, changes in methylation at some imprinted regions, but not in the majority [19,48]. It is unknown if these changes occur due to IVF since imprinting defects have been previously described in sperm of infertile men, including hypomethylation of the H19 CTCF6 DMR [49]. Loke et al. [16] reported that hypomethylation at this locus in buccal epithelium of newborns in the IVF group was driven by the subgroup conceived by ICSI. However, it is difficult to dissect whether the observed effect on DNA methylation of ICSI-conceived newborns is due to the technique itself or to male infertility. Whitelaw et al. [50] found higher levels of SNRPN methylation in buccal cells of ICSIconceived newborns and these were associated with longer duration of infertility in the parents. In our data, we observed that the difference at the H19 CTCF6 DMR was greater when considering any type of medical help during conception, supporting the idea that parental subfertility is the driver of methylation changes at this region. Information about the indication for assisted reproductive technology, the use of donor eggs or sperm, and the fertility status of parents in the control group would be required to further assess the effect of parental subfertility. Adverse perinatal outcomes and increased frequency of imprinting disorders have also been observed in offspring of couples with a history of subfertility that were able to conceive naturally [51][52][53]. However, studies that controlled for parental subfertility by comparing siblings in which one was conceived naturally and the other by IVF also observed an effect [54]. It is likely, therefore, that both parental subfertility and IVF may induce epigenetic changes, as observed in another genome-wide study that found DNA methylation differences between IVF-conceived newborns and a group conceived through intrauterine insemination (infertile controls), but also between the latter and naturally conceived newborns (fertile controls) [19]. In addition, a study looking at 37 candidate CpG sites identified seven that were differentially methylated when comparing an IVF-conceived group born to parents without male infertility that used donor oocytes to naturally conceived newborns [55].
Finally, two IVF-DMRs associated with infertility (TNP1 and C9orf3) were targeted for validation. Differential methylation was validated at the C9orf3 gene. However, validation of the TNP1 region was hampered by our inability to target CpG sites within the mostassociated DMR in this locus. We attempted validation at TNP1 by targeting CpG sites in the neighbouring 500bp bin and observed consistent direction of association close to nominal significance.
In this study the non-IVF group included a set of twins conceived after GIFT and another set conceived after ovarian stimulation not followed by IVF. GIFT and ovarian stimulation are fertility treatments not equivalent to IVF since fertilisation still occurs in the fallopian tubes. We showed that our results were not affected by the inclusion of these data, potentially because they were represented in small numbers, only four out of 58 samples.

Limitations
There are several limitations in this study. First, it is known that cell composition may represent a confounding variable in EWAS [56]. Our results use principal component analysis anticipating that these will capture cell heterogeneity, and follow-up of our findings in a subset of twins with available cell counts showed that the majority of findings remained significant after adjustment for cell heterogeneity. Second, although MeDIP-seq has the strength of genomewide coverage, it lacks base-pair resolution, instead generating methylation scores across genomic regions. However, it has been reported that methylation of neighbouring CpG sites is correlated over distances up to 1000 bp [57], suggesting that the approach may be able to capture a good proportion of the methylation variance in a genomic region. Third, although this study includes a sample size larger than most previous studies exploring IVF, contemporary EWAS study designs generally require larger numbers of cases and controls to achieve sufficient power to detect small to moderate effect sizes [21,58]. Lastly, our approach cannot conclusively determine the cause of the observed IVF-associated methylation changes. Future studies of IVFassociated regions in animal models, where genetic differences and infertility diseases can be discarded, could help identify if these changes were caused by IVF itself.

Conclusions
We observed evidence for differences in DNA methylation between IVF and non-IVF twins on a genome-wide scale. A strength of this study design is that it allowed us to also estimate the contribution of genetic and environmental factors towards DNA methylation levels at the IVF-associated loci. The inclusion of only twin pregnancies also avoided biases present in studies that consider single and multiple pregnancies together. Multiple pregnancies are more common after IVF. Therefore, the differences observed when studying singleton and twin births together may be confounded with the higher risks of adverse perinatal outcomes in multiple pregnancy births, rather than IVF itself. Nevertheless, we were unable to dissect whether methylation changes were likely caused by IVF, or were due to the underlying parental subfertility, or other factors. These scenarios require further study exploring the stability of these DMRs over time, their relation with gene expression, and their potential role in health and disease.

Additional file
Additional file 1: Supplementary tables and figures. Table S1 Pregnancy complications and other outcomes. Table S2 PCR assays. Table S3 Re-analysis of FDR 25% WBC IVF-DMRs excluding other fertility treatments from non-IVF group. Table S4 Re-analysis of FDR 25% WBC IVF-DMRs in subset with cell counts. Figure S1 C9orf3 FDR 25% IVF-DMR. Figure S2 H19 CTCF6 IVF-DMR replication. Figure S3 CpG sites targeted for validation. Figure S4 IVF versus non-IVF DNA methylation differences using Sequenom's EpiTYPER technology. Figure S5 ICSI versus conventional IVF DNA methylation differences using Sequenom's EpiTYPER technology.