DNA methylation differences at birth after conception through ART.

STUDY QUESTION: Is there a relation between ART and DNA methylation (DNAm) patterns in cord blood, including any differences between IVF and ICSI? SUMMARY ANSWER: DNAm at 19 CpGs was associated with conception via ART, with no difference found between IVF and ICSI. WHAT IS KNOWN ALREADY: Prior studies on either IVF or ICSI show conflicting outcomes, as both widespread effects on DNAm and highly localized associations have been reported. No study on both IVF and ICSI and genome-wide neonatal DNAm has been performed. STUDY DESIGN, SIZE, DURATION: This was a cross-sectional study comprising 87 infants conceived with IVF or ICSI and 70 conceived following medically unassisted conception. The requirement for inclusion in the study was an understanding of the Swedish language and exclusion was the use of donor gametes. PARTICIPANTS/MATERIALS, SETTING,
METHODS: Participants were from the UppstART study, which was recruited from fertility and reproductive health clinics, and the Born into Life cohort, which is recruited from the larger LifeGene study. We measured DNAm from DNA extracted from cord blood collected at birth using a micro-array (450k array). Group differences in DNAm at individual CpG dinucleotides (CpGs) were determined using robust linear models and post-hoc Tukey's tests. MAIN RESULTS AND THE ROLE OF CHANCE: We found no association of ART conception with global methylation levels, imprinted loci and meta-stable epialleles. In contrast, we identify 19 CpGs at which DNAm was associated with being conceived via ART (effect estimates: 0.5-4.9%, PFDR < 0.05), but no difference was found between IVF and ICSI. The associated CpGs map to genes related to brain function/development or genes connected to the plethora of conditions linked to subfertility, but functional annotation did not point to any likely functional consequences. LIMITATIONS, REASONS FOR CAUTION: We measured DNAm in cord blood and not at later ages or in other tissues. Given the number of tests performed, our study power is limited and the findings need to be replicated in an independent study. WIDER IMPLICATIONS OF THE
FINDINGS: We find that ART is associated with DNAm differences in cord blood when compared to non-ART samples, but these differences are limited in number and effect size and have unknown functional consequences in adult blood. We did not find indications of differences between IVF and ICSI. STUDY FUNDING/COMPETING INTEREST(S): E.W.T. was supported by a VENI grant from the Netherlands Organization for Scientific Research (91617128) and JPI-H2020 Joint Programming Initiative a Healthy Diet for a Healthy Life (JPI HDHL) under proposal number 655 (PREcisE Project) through ZonMw (529051023). Financial support was provided from the European Union's Seventh Framework Program IDEAL (259679), the Swedish Research Council (K2011-69X-21871-01-6, 2011-3060, 2015-02434 and 2018-02640) and the Strategic Research Program in Epidemiology Young Scholar Awards, Karolinska Institute (to A.N.I.) and through the Swedish Initiative for Research on Microdata in the Social And Medical Sciences (SIMSAM) framework grant no 340-2013-5867, grants provided by the Stockholm County Council (ALF-projects), the Strategic Research Program in Epidemiology at Karolinska Institutet and the Swedish Heart-Lung Foundation and Danderyd University Hospital (Stockholm, Sweden). The funders had no role in study design, data collection, analysis, decision to publish or preparation of the manuscript. The authors declare no competing interests. TRIAL REGISTRATION NUMBER: N/A.

Introduction

The usage of ART is increasing worldwide and its possible health consequences are a topic of intense study. Evidence is building for an association of ART with long-term health outcomes for the offspring, including autism spectrum disorders (Liu ) and cardiovascular health (higher blood pressure, suboptimal cardiac diastolic function and vessel thickness) (Guo ). A body of evidence links ART to short-term health effects, including low birthweight, placenta-associated anomalies (Vermey ), pregnancy complications (Qin ) and congenital malformations and imprinting disorders (Turkgeldi ). However, it is still debated if these associations stem from the application of ART or the underlying infertility leading couples to ART (Luke ). Animal studies have provided potential molecular mechanisms for these observations by highlighting that the ART procedure may induce changes to epigenetic marks (Morgan ; Wang ). Epigenetic marks, such as DNA methylation (DNAm), influence the transcription potential of genomic regions (Jaenisch and Bird, 2003). Human studies likewise point toward a link between early development, DNAm and (late-life) phenotypes, but the causality of these associations remains unknown (Tobi ).

To date, several studies on DNAm and ART have been performed. Multiple studies have focused on candidate gene regions and/or global methylation levels (Lazaraviciute ; Canovas ). Genome-wide efforts have focused on samples taken from extra-embryonic lineages (Xu ; Choufani ) and when focused on material from the infant itself suffered from a small sample size (Melamed ) or batch effects correlating with ART status (Estill ). Larger studies have focused on specific ART techniques, namely ICSI (El Hajj ) or IVF (Castillo-Fernandez ) or IVF and the less invasive gamete intra-fallopian transfer (GIFT) coupled with IUI (Novakovic ). It has been hypothesized that the technique used may matter for the possible consequences on DNAm patterns (Loke and Craig, 2016). Overall, these studies report conflicting outcomes, with two reporting widespread associations between ART and (cord) blood DNAm at birth (El Hajj ), which appear to fade with age (Novakovic ), while another study reported a DNAm difference at a single genomic locus only (Castillo-Fernandez ).

We undertook an epigenome-wide association study (EWAS) on ART, comparing DNAm patterns in cord blood of children conceived via ART (N = 87) (Iliadou ) with that of children from medically unassisted conceptions (MUC, N = 70) (Almqvist ). Since we had detailed information on the ART technique used for 77 of the children conceived via ART, we investigated possible different outcomes for IVF and ICSI on DNAm, which has been hypothesized to be important (Loke and Craig, 2016), but is yet to be tested. In addition, we explored the possible functional consequences of methylation differences using external datasets.

Materials and methods

Study subjects

ART controls: Born into Life

Study participants in the prospective longitudinal birth cohort study Born into Life were recruited from the larger LifeGene study (Almqvist ). LifeGene is a prospective cohort study with the aim to combine advances in modern biotechnology with information on individuals’ health and lifestyle collected through web-based questionnaires and biosamples. Recruitment of LifeGene participants aged 18–45 years was based on random sampling in the general population and they were invited to include their adult household members. Between the years 2010 and 2012, pregnant women who were already participating in the LifeGene study and living in Stockholm County were recruited to Born into Life (Smew ). The inclusion criteria for Born into Life were that the women had responded to baseline questionnaires from the LifeGene study, were pregnant and gave written informed consent. They were recruited both before and after 10–14 gestational weeks, but no later than 26–28 weeks. Originally, 107 pregnant women were included in Born into Life. We only included the 77 women for which cord blood had been successfully collected in this study.

Women in the Born into Life were asked to complete questionnaires regarding pregnancy, lifestyle and health at 10–14 and 26–28 gestational weeks. Data on maternal self-reported smoking during pregnancy and BMI (kg/m2) obtained at the first antenatal care visit in gestational weeks of 5–12 were collected from birth records. Data regarding highest attained educational level, ranging from mandatory secondary school to high school, university or other, were retrieved from the baseline LifeGene questionnaires. Maternal age at delivery, in years, was calculated from mothers’ date of birth. Parity and data regarding the infants’ sex, gestational age in weeks and mode of delivery, defined as vaginal delivery or cesarean section, were also collected from the birth records. Cord blood and placenta samples were obtained from delivery. Cord blood was aspirated 2 min after birth by the assisting midwife into a test tube (EDTA) and kept at room temperature while awaiting transportation. The birth records regarding both mother and child at delivery were collected from Danderyd Hospital, Stockholm.

ART cases: UppstART

The UppstART study has been described in detail elsewhere (Iliadou ). Participants of the UppstART study were recruited from three of the four fertility and reproductive health clinics in Stockholm (one public, two private) and one private clinic in Uppsala county, which also serves a large volume of patients from Stockholm. Recruitment took place from September 2011 to December 2013. IVF treatment(s) of the participants was followed until December 2014 or drop-out/consent withdrawn (n = 4), whichever came first. The participants were asked to answer a web-based baseline questionnaire within a few days of their clinic visit and prior to their IVF treatment start, which included an extensive list of questions on sociodemographic, anthropometric and lifestyle factors. Once the IVF treatment began and the participants reached the stage of oocyte retrieval, they were asked to respond to a second online follow-up questionnaire, with a shorter version of the baseline questionnaire, to identify any changes in lifestyle factors since the initiation of their treatment. Staff at seven delivery units in Stockholm and Uppsala were recruited to assist in collecting samples from UppStART participants during delivery. Delivery clinics were provided with a sample collection kit including tubes for collection of cord blood for DNA extraction (EDTA). Samples were stored in −20°C freezers at the delivery units until they were collected by UppStART study staff and deposited into the KI Biobank.

DNAm measurements

Genome-wide DNAm data were generated using the Illumina Infinium Human Methylation 450K BeadChip (450k array). A total of 500 ng of genomic DNA isolated from cord blood was bisulfite treated using the EZ-96 DNA methylation kit (Zymo Research, Orange County, CA, USA). We used the D-optimum criterion to assign samples over two 96-well plates and individual 450k arrays, ensuring even distributions of ART cases and controls, ART methods (IVF or ICSI), sex, gestational age and birth month across the two 96-well bisulfite plates and each 450k array. The 450K arrays were measured at ServiceXS (Leiden, The Netherlands). The quality of the generated 450K array data was assessed using both sample dependent and sample independent quality metrics using the Bioconductor package MethylAid (van Iterson ) with default settings. We used the Bioconductor package omicsPrint (van Iterson ) to check for sample duplications and mixtures, the absence of family relations and to identify the probes where single nucleotide polymorphisms (SNPs) influence the measurement of the methylated or unmethylated signal (Zhou ). These probes were kept as, although not our main interest, any difference in those probes might hint at interesting genetic differences. Sample sexes were checked using the X-chromosomal CpGs. One sample was not the correct sex and was deleted. We used principal component analysis (PCA) and hierarchical clustering on the raw autosomal beta values to search for outliers and suspect patterns in the dataset, finding none. The cell proportions of cord blood were imputed using the “Identifying Optimal Libraries” (IDOL) algorithm (Koestler ) on the ‘FlowSorted.CordBloodCombined.450k’ reference set using the estimateCellCounts2() function in the FlowSorted.Blood.EPIC R package (Salas ; Gervin ). Normalization of the dataset was performed by NOOB background and color correction in combination with Functional Normalization (Fortin ) using four principal components (eigenvalue >1). All measurements with <3 beads (0.05% of probes), <1 intensity value (0.012% of probes) and a detection P-value >0.01 (0.11% of CpGs) were set as missing. The measurement success rate per sample was >99%. CpGs with probes that did not map to unique genomic locations (Chen ) or with a <95% measurement success rate (0.23% of CpGs) were then removed. We used custom scripts that add on the functions from the minfi package (Aryee ) and implements parallelization where possible (https://github.com/molepi/DNAmArray).

Outlier detection was performed by PCA and hierarchical clustering on the autosomal beta values. Hierarchical clustering identified one possible outlier. This individual had a very low gestational age (30 weeks), the lowest in the dataset, and this clustering is therefore likely based on biological reasons. High-quality DNAm data were obtained for 157 individuals, including 70 MUC from Born into Life, six ART from Born into Life (specific type of ART unknown) and 81 ART from UppStART (specific type of ART unknown for N = 4), for a total of 441 836 autosomal CpGs. The datasets generated and analyzed during the current study are not publicly available due to Swedish privacy and data safety laws but are available from A.N.I. and C.A. on reasonable request and after meeting legal requirements.

Measures of global methylation

Genome-wide average DNA methylation (GWAM) (Li ) was calculated by averaging all beta-value measurements across the autosomes for each individual. The R package REMP (Zheng ) was used to infer DNAm at either ALU or LINES-1 sequences and then the average across all ALU or LINES-1 elements was calculated for each neonate.

Transcription factor binding site enrichment

We calculated transcription factor binding site (TFBS) enrichments using the R package PWMEnrich (Stojnic and Diez, 2018) and using binding motifs from motifDB. We calculated position weight matrices using the DNA base background frequencies calculated for the CpGs tested with 25 bp flanking sequences. Enrichment was tested relative to this background.

Statistics

We compared MUC and ART cohort descriptives via Kruskal–Wallis (maternal BMI, maternal socio-economic status (SES), years until index pregnancy, CD4T, natural killer (NK) and nucleated red blood cells: nRBC), Student’s t-test (maternal and gestational age, birthweight, CD8T, monocytes (mono), granulocytes (gran) and B cell (Bcell) proportions) assuming equal or unequal variance, where appropriate, and chi-square tests (child sex, parity, preterm births and maternal smoking). DNAm was always analyzed as a beta value, which is reported as a percentage (e.g. ×100). GWAM, ALU and LINES-1 methylation levels were compared via t-test (MUC versus ART) or ANOVA (MUC versus IVF versus ICSI).

We used the R package cate (Wang and Zhao, 2015) to identify latent variables (e.g. ‘hidden variables’) influencing the relation between DNAm and ART (yes/no) or ART method (IVF/ICSI/MUC). Two latent variables were identified (P < 0.001) for ART and four latent variables for ART method (P < 0.001). In both cases, the first two latent variables correlated with the proportion of nRBC in blood (rho > 0.94, P < 0.001) and the other two latent variables correlated with CD4T or B cell proportions in blood and several batch effects (rho > 0.3, P < 0.001) arguing, together with directed acyclic graphs (DAG) (Supplementary Data), for a basic model adjusting for batch effects and imputed cell proportions.

We performed EWAS on being conceived by ART or not using robust linear regression (rlm) from the R MASS package with White’s estimator for robust standard errors, as implemented in the R package sandwich (Zeileis, 2006) (which leads to a model robust for outlying beta values and heteroscedasticity). To test for differences in neonates conceived via IVF or ICSI or without medical assistance (MUC), we used type II ANOVA with White’s estimator for robust standard errors to account for any heteroscedasticity and unequal variances between groups using the R car package. We employed heteroscedasticity robust Tukey’s tests to test for DNAm differences between IVF, ICSI and MUC post-hoc. In all instances, we used the beta value of an individual CpG dinucleotide as the outcome and ART status as the dependent. We adjusted for the sex of the individual, height of the sample on the 450k array glass slide (coded as a continuous variable from 1 to 6), scan batch and bisulfite plate and the imputed cell proportions (CD4T, CD8T, Mono, NK, Bcell, Gran, nRBC). We ran additional sensitivity analyses by omitting the cell proportions or expanding the adjustment with minimal adjustment sets coming from DAG analysis (Supplementary Data) with the following models: years until the index pregnancy, maternal age, pre-pregnancy BMI, educational attainment as proxy for SES, smoking history, and gestational length and parity. In addition, we performed look-ups in other published meta-EWAS of prenatal exposures.

Regions were tested with the same adjustments as defined above in a linear mixed model with a random effect for individual and a factorial covariate denoting the specific CpG dinucleotide of each DNAm measurement in the R lme4 package using compound symmetry as correlation structure and the R lmerTest package for the Satterthwaite's approximation of the degree of freedom of the fixed effects within each linear mixed model.

All analyses were performed in R (version 3.6.1) (R Core Team, 2019). All P-values reported are two-sided and multiple testing correction is done using false discovery rate (FDR).

Ethical approval

Ethical approval for Born into Life was granted by the Regional Ethics Review Board in Stockholm, Sweden. Written informed consent from both parents was obtained for all study participants. The UppstART study has been approved by the regional ethics review board at the Karolinska Institutet (Dnr 2011/230-31/1, Dnr 2011/1427-32, Dnr 2012/131-32, Dnr 2012/792-32, 2013/1700-32, 2014/1956-32, Dnr 2015/1604-32). Women and their partners were approached by the clinic nurse and asked to participate in the study. To facilitate the process of informed consent, the couple was provided with information approved by the regional ethical board, both verbally and in written format, about the purpose of the study, methods, possible risks, and that participation was voluntary. Additionally, participants were informed that they could withdraw from the study at any time with no impact to their medical care. The requirement for inclusion in the study was an understanding of the Swedish language and exclusion was the use of donor gametes.

Results

Pre-processing and normalization resulted in DNAm data of 157 infants from cord blood, 87 of which were conceived with the help of ART (Table I). We had detailed information on the ART technique used for 77 of the 87 newborns conceived by ART and with DNAm data, of which 44 were conceived via IVF and 33 via ICSI. The mothers who conceived with the help of ART were not different in terms of BMI (P = 0.25) and age (P = 0.25), but had a lower SES (P < 0.001) and more were former smokers (P < 0.001) than mothers with a MUC. It took on average 2 years longer for the mothers who conceived through ART to become pregnant than for mothers with MUC (P < 0.001). There was no difference between ART and MUC in length of gestation (P = 0.16) and parity (P = 0.20). The percentage of male newborns in the MUC group was higher (61.4% versus 50.6%), although this difference was not significant (P = 0.23).

Table I

Characteristics of the cohorts in a study of DNAm and ART.

Variable	MUC	ART	P-value	IVF	ICSI
N	70	87		44	33
N frozen1	–	10		7	3
N blastocyst (N frozen)2	–	8		5 (3 frozen)	3 (2 frozen)
BMI in kg/m2 (SD)	22.8 (3.5)	23.4 (3.4)	0.25	23.5 (3.3)	23.4 (3.0)
Age in years (SD)	32.4 (3.6)	33.1 (3.6)	0.25	33.0 (3.7)	33.0 (3.7)
SES3 (SD)	3.0 (0.3)	2.6 (0.8)	<0.001*	2.7 (0.7)	2.3 (0.9)
Smoking (%)			<0.001*
Never	65 (92.8)	56 (64.4)		27 (61.4)	23 (69.7)
Ever	4 (5.8)	26 (29.9)		17 (38.6)	9 (27.3)
NA	1 (1.4)	5 (5.7)		0 (0.0)	1 (3.0)
Years to index pregnancy			<0.001*
Mean (SD)	0.2 (0.7)	2.6 (1.3)		2.7 (1.2)	2.4 (1.5)
Min.–max.	0–4	0–7		1–7	0–7
Parity (%)			0.20
First	47 (67.1)	62 (71.4)		34 (77.3)	22 (66.7)
Second	21 (30.0)	17 (19.5)		10 (22.7)	7 (21.2)
Third	2 (2.9)	1 (1.1)		0 (0.0)	1 (3.0)
NA	0 (0.0)	7 (8.0)		0 (0.0)	3 (9.1)
N male children (%)	43 (61.4)	44 (50.6)	0.23	22 (50.0)	16 (48.5)
Gestational age in weeks	39.1 (1.6)	39.5 (2.0)	0.16	39.2 (2.3)	39.7 (1.6)
N preterm births4	3	3	0.99	3	0

The number of embryos that were frozen and thawed before placement in utero.

The number of embryos that were cultured in vitro to the blastocyst stage before placement in utero. Five of these blastocysts were frozen and thawed.

SES graded on a four-level scale based on the highest attained educational level.

Number of infants born preterm (gestational age of <36 weeks).

P < 0.001 in a comparison between ART and MUC using Kruskal–Wallis (SES, years until index pregnancy) and chi-square (maternal smoking) tests.

DNAm, DNA methylation; MUC, medically unassisted conception; SES, socio-economic status.

DNAm is a key determinant of cell identity (Jaenisch and Bird, 2003) and cord blood consists of multiple cell types that influence DNAm variation. Seven major cell types, including nRBCs, were imputed from the genome-wide DNAm data of the newborns (Gervin ). In concordance with a prior report on unassisted and ICSI newborns (El Hajj ), there were no differences in these seven cell proportions between ART and MUC cord blood samples (Pnominal > 0.13, Supplementary Fig. S1).

Global DNAm comparison between ART and unassisted conceptions

First, we investigated GWAM (Li ) in the cord blood of neonates conceived via IVF, ICSI or MUC. There was no difference between these three groups (P = 0.89, Supplementary Fig. S2). Next, we investigated the average methylation of ALU and LINES-1 elements, an often used proxy for global DNAm (Zheng ), finding no difference between groups (P > 0.2, Supplementary Fig. S3).

DNAm comparison between ART and MUC

We compared DNAm at 441 836 autosomal CpGs between 87 ART and 70 MUC neonates (Supplementary Fig. S4). Nineteen CpGs were associated with being conceived via ART in a model adjusting for sex, batches and cell heterogeneity (PFDR < 0.05, Table II, Supplementary Fig. S5). The mean absolute differences between MUC and ART are small from a molecular point of view (β = 0.48–4.89%), and medium to large when looked at in terms of effect size (0.44–1.00 SD), which is similar to that found for other prenatal exposures such as prenatal smoking (Joubert ), folate use (Joubert ), famine (Tobi ) and hypertension and pre-eclampsia (Kazmi ). The associations extended to neighboring probes (Fig. 1) for the proximal promoter of AK054845, a gene of unknown function (five CpGs Pnominal < 0.05), and the proximal promoter of a smaller variant of GABRB3 (NM_001191321.2), which is expressed in the brain (eight CpGs Pnominal < 0.05) and has a neuro-developmental role (Tanaka ).

Table II

Results of the ART epigenome-wide association study.

CpG ID	Location hg19	Methylation (SD)1	Nearest gene2	Distance2	Estimate (SE)3	Effect size (SD)	P	P FDR
cg27266479	chr1:9294882	32.8 (1.9)	H6PD	0	−1.94 (0.27)	−1.00	3.18e−13	8.28E−08
cg04811592	chr3:69834386	89.6 (1.6)	MITF	0	0.99 (0.21)	0.62	1.59e−06	0.039
cg24959663	chr5:10578618	71.2 (5.2)	ANKRD33B	0	3.91 (0.54)	0.74	3.75e−13	8.28E−08
cg22916646	chr5:162672583	69.0 (2.6)	–	–	2.18 (0.36)	0.83	1.22e−09	1.08E−04
cg01500567	chr6:44355777	4.2 (0.6)	CDC5L	0	0.48 (0.1)	0.74	5.46e−07	0.020
cg00478390	chr7:150703765	89.4 (2.1)	NOS3	0	−0.92 (0.19)	−0.44	7.37e−07	0.024
cg03207674	chr7:1523569	92.7 (1.1)	INTS1	0	0.72 (0.14)	0.64	1.92e−07	9.43E−03
cg17123384	chr7:83379152	82.0 (4.0)	–	–	2.78 (0.51)	0.70	6.39e−08	3.53E−03
cg19347588	chr10:3868336	92.9 (1.2)	KLF6	40862	0.81 (0.16)	0.66	4.49e−07	0.018
cg07569385	chr13:20766226	10.8 (2.3)	GJB2	0	1.53 (0.31)	0.68	1.24e−06	0.036
cg06485032	chr13:22615064	71.4 (5.3)	AK054845	0	−3.47 (0.67)	−0.65	2.59e−07	0.011
cg13051607	chr15:22956714	86.1 (2.1)	CYFIP1	0	1.39 (0.24)	0.65	1.15e−08	7.24 E−04
cg01251603	chr15:26874098	76.7 (5.3)	GABRB3	0	−4.89 (0.8)	−0.92	1.18e−09	1.08 E−04
cg15066197	chr15:26874202	85.8 (4.8)	GABRB3	0	−4.71 (0.79)	−0.99	2.48e−09	1.83 E−04
cg14859324	chr15:26874363	92.3 (2.9)	GABRB3	0	−2.26 (0.46)	−0.78	7.68e−07	0.024
cg06450634	chr16:30430044	67.0 (3.3)	ZNF771	0	3.26 (0.46)	0.98	2.15e−12	3.17E−07
cg08783253	chr17:40996565	73.7 (5.3)	AOC2	42	−2.93 (0.61)	−0.55	1.48e−06	0.038
cg16771467	chr18:55315872	98.5 (0.2)	ATP8B1	0	0.11 (0.02)	0.55	1.34e−06	0.037
cg14560133	chr19:51199453	28.8 (2.0)	SHANK1	0	−1.42 (0.3)	−0.71	2.07e−06	0.048

The results of the ART versus MUC epigenome-wide association study are ordered on genomic location (hg19). A negative estimate means a lower DNAm in the ART group.

Mean methylation (beta value*100) and SD.

Nearest gene within 100k nucleotides and the distance to the transcription start site. Human Genome Nomenclature consortium approved gene names: H6PD, hexose-6-phosphate dehydrogenase/glucose 1-dehydrogenase; MITF, melanocyte inducing transcription factor; ANKRD33B, ankyrin repeat domain 33B; CDC5L, cell division cycle 5 like; NOS3, nitric oxide synthase 3; INTS1, integrator complex subunit 1, KLF6, Kruppel-like factor 6; GJB2, gap junction protein beta 2; AK054845, transcript for long intergenic non-protein coding RNA 540; CYF1P1, cytoplasmic FMR1-interacting protein 1; GABRB3, Gamma-Aminobutyric Acid Type A Receptor Subunit Beta3; ZNF771, Zinc Finger Protein 771; AOC2, Amine Oxidase Copper Containing 2; ATP8B1, ATPase Phospholipid Transporting 8B1; SHANK1, SH3 and Multiple Ankyrin Repeat Domains 1.

Estimate and SE followed by columns with the effect size (in SDs) and P-value, with and without false discovery rate (FDR) correction of the estimate of ART versus MUC in the model: Beta ∼ ART (yes/no) + height on micro-array slide + scan batch + bisulfite plate + sex + CD4T + CD8T + Mono + Bcell + NK + Gran + nRBC. NK: natural killer cells; nRBC, nucleated red blood cells, Mono: monocytes, Gran: granulocytes.

Figure 1.

Manhattan plot of the epigenome-wide association study for ART status. A Manhattan plot showing the –log10 P-values (y-axis) for the association of DNAm at individual CpG dinucleotides across the 22 autosomal chromosomes (x-axis). CpG dinucleotides 5kb up- and downstream of the lead association have been colored in the same color as the lead association. H6PD, hexose-6-phosphate dehydrogenase/glucose 1-dehydrogenase; MITF, melanocyte inducing transcription factor; ANKRD33B, ankyrin repeat domain 33B; CDC5L, cell division cycle 5 like; NOS3, nitric oxide synthase 3; INTS1, integrator complex subunit 1; KLF6, Kruppel-like factor 6; GJB2, gap junction protein beta 2; AK054845, transcript for long intergenic non-protein coding RNA 540; CYF1P1, cytoplasmic FMR1-interacting protein 1; GABRB3, Gamma-Aminobutyric Acid Type A Receptor Subunit Beta3; ZNF771, Zinc Finger Protein 771; AOC2, Amine Oxidase Copper Containing 2; ATP8B1, ATPase Phospholipid Transporting 8B1; SHANK1, SH3 and Multiple Ankyrin Repeat Domains 1; DNAm, DNA methylation.

Omitting the adjustment for cell heterogeneity had little to no influence on the effect estimates for all 19 CpGs (Supplementary Table SI) just like additional adjustment for maternal and pregnancy characteristics (Supplementary Table SII) or exclusion of the neonates of which the embryo was frozen and/or transferred as a blastocyst (Supplementary Table SIII). The 19 CpGs were not among those previously reported in other EWASs (FDR corrected P-value of <0.05) with prenatal smoking (Joubert ), folate use (Joubert ), maternal hypertension and pre-eclampsia (Kazmi ), maternal BMI (Sharp ) and birthweight (Küpers ). Moreover, a look-up of these previously reported CpGs in our study did not yield more nominally associated CpGs with ART than may be expected by chance (P > 0.07).

DNAm comparison between IVF, ICSI and MUC

To investigate any differences between IVF and ICSI, we performed post-hoc Tukey’s tests for DNAm differences between neonates conceived via IVF and ICSI for these 19 CpGs. No differences were found (PFDR-19 tests > 0.15). Next, we extended our analysis to all CpGs by performing an ANOVA analysis to test for differences in DNAm across the 441 836 autosomal CpGs between neonates conceived with IVF (N = 44), ICSI (N = 33) or MUC (N = 70). This yielded five CpGs already identified in the EWAS for ART status and in all cases, a post-hoc Tukey’s test showed that there was a DNAm difference for all five CpGs for both the MUC versus IVF and MUC versus ICSI comparison, thus finding no evidence for IVF or ICSI specific associations.

Functional annotation

We performed KEGG and Gene Ontology tests (Phipson ), finding no enrichment for the genes linked to the 19 CpGs associated with ART, nor when we relaxed the significance threshold to suggestive associations (P < 10−5, N = 37). We did not find significant enrichments of single transcription factor (TF) binding motifs (Stojnic and Diez, 2018), although multiple CpGs overlapped TF binding sites for TFs with a role in early development, such as TFAP2C (Sharma ) overlapping cg27266479 mapping to H6PD, and RAX (Bennett ) overlapping cg16771467 and mapping to ATP8B1. Look-up of the 19 CpGs in reference data (Bonder ) uncovered no link between DNAm and gene expression in adult blood.

We measured DNAm in cord blood and average DNAm levels may vary between tissues, but still reflect variation in other tissues (Slieker ). There was little to no correlation of methylation at the 19 CpGs between adult whole blood and tissues from 16 cadavers (Slieker ; Relton ), with the notable exception of cg14560133 (SHANK1, β = −1.42% (SE = 0.3%), P = 2.1 × 10−6), which showed moderate correlation with 5/8 tissues (rho > 0.53, P < 0.04, Fig. 2). The three CpG dinucleotides at gamma-aminobutyric acid receptor subunit beta-3 (GABRB3) (cg01251603, cg15066197, cg14859324), a gene active in the brain, showed weak to moderate correlation between adult blood and four brain regions (rho = 0.39–0.62, P < 0.001, N = 71–74) in another reference dataset (Hannon ).

Figure 2.

Correlation between adult blood DNAm and internal tissues. The spearman correlation plotted for eight tissues from the body (N = 16, P < 0.1) and four from the brain (N > 71, P < 0.01) in colors ranging from bright red (rho = −1.0) to dark blue (rho = 1.0). The light gray lines along the rows identify the three CpG dinucleotides from the gamma-aminobutyric acid receptor subunit beta-3 (GABRB3) region associated with ART. The light gray line running along one column denotes the separation of the two tissue reference datasets used. ‘Fat Sub’: subcutaneous fat, ‘Left Myocard’: left myocardium, ‘Muscle Skel’: skeleton muscle, ‘prefr.cortex’: prefrontal cortex, ‘ent.cortex’: entorhinal cortex, ‘s.t.gyrus’: superior temporal gyrus.

Imprinted regions and meta-stable epialleles

Multiple candidate gene studies have looked at DNAm or gene expression differences of imprinted genes between ART and MUC with mixed results. Therefore, we looked at 374 CpGs on the 450k array known to overlap imprinted differentially methylated regions (Yuen ). There was no overlap with the CpGs showing a (suggestive) association with ART in our study (P < 10−5, N = 37). These 374 CpGs were distributed over 59 regions. There was no difference in DNAm between MUC and ART infants when we tested each of these 59 regions (PFDR > 0.86). Similarly, so-called meta-stable epialleles (MEs) are hypothesized to be especially sensitive to the early prenatal environment (Kessler ). There were 187 CpGs overlapping MEs, none of which were associated with ART (PFDR > 0.06).

Prior studies

We cross-referenced our results with those from prior studies on one form of ART. The sole region associated with IVF in the methylated DNA immunoprecipitation sequencing study from Castillo-Fernandez was not covered by our genome-scale screen with the 450k array. Next, El Hajj found 4730 CpG dinucleotides associated with ICSI after FDR adjustment (which were not available from their article) of which two regions, out of five regions selected for validation with pyrosequencing, were validated in an additional sample. The seven CpGs in these two regions, which cover the proximal promoters of ATG4C and SNORD114-9, were not associated with ART in our study (P > 0.23) or with ICSI only (post-hoc Tukey’s test: P > 0.40).

Using the 850k array, Novakovic identified 2340 CpG dinucleotides at which DNAm in Guthrie card blood was associated with IVF and GIFT coupled with IUI, 1228 of which were covered in our 450k study. Of these, 38 were nominally associated (P < 0.05) and with the same direction of effect in our study. This decreased to 12 CpGs when we tested for DNAm differences between MUC and IVF only. Novakovic also identified three differentially methylated regions where the difference in methylation associated with IVF remained in an adult sample. CpG dinucleotides at CHRNE (chr17:4803506–4805392), PRSS16 (chr6:27185896–27186199) and TMEM18 (chr2:731073–732037) showed a similar direction of effect in our study as that of Novakovic and the association was almost nominally significant for PRSS16 (cg10279314: β = 2.3% (SE = 1.2%), P = 0.057, cg09395805: β  = 3.1 (SE = 1.6), P = 0.049, cg07555084: β  = 2.2 (SE =1.2), P = 0.074, Supplementary Table SIV).

Discussion

We identified 19 CpGs at which DNAm was associated with ART and the association extended to neighboring CpG dinucleotides at AK054845 and GABRB3. We found no indication for a disparity between IVF or ICSI on DNAm patterns. The associations were robust to adjustment for cellular heterogeneity and maternal characteristics, and did not overlap with those loci associated with other common prenatal exposures. We found no evidence that imprinted loci and MEs are especially sensitive to ART. The 19 CpGs can be annotated to genes with a role in the brain or the plethora of conditions that can be linked to subfertility, but we found no clear functional implications of the variation in DNAm at these 19 CpGs.

The 19 CpG dinucleotides could be mapped to 15 genes. Four of the 19 were located in a proximal promoter and here the neighboring CpG dinucleotides were also nominally associated with ART. In other cases, the associations were either limited to one CpG dinucleotide, which may stem from the sparse coverage of the 450K array across the methylome and in other cases is consistent with the fact that DNAm acts through the altering of the binding potential of a specific TFBS (Bonder ). Indeed, none of the CpGs were located in a CpG island, which is consistent with data from other exposures showing that associations between DNAm and prenatal and/or environmental conditions are enriched at regulatory regions such as enhancers and CpG poor promoters that show intermediate levels of methylation (Tobi ).

The mean absolute differences between MUC and ART are small from a molecular perspective (β = 0.48–4.89%), while effect sizes were medium to large effect from an epidemiological perspective (0.44–1.00 SD). This is in line with results from (meta-)analyses on other prenatal exposures such as prenatal smoking (Joubert ), folate use (Joubert ), famine (Tobi ) and hypertension and pre-eclampsia (Kazmi ). Such differences are thought to exert an effect through the modulation of gene networks (Stoger, 2008; Jiao ) or to mark other larger molecular differences in a regulatory region, such as in histone modifications, and animal experimental data suggest that small absolute differences in methylation of medium effect size may still exert a measurable effect on gene expression (Lillycrop ).

Several of the CpGs can be linked to genes with a role in brain function and/or development. ZNF771 (cg06450634) has been identified as a driving factor behind a large gene network in the brain (Maulik ) and CYFIP1 (cg13051607) (Hsiao ), GABRB3 (cg01251603, cg15066197, cg14859324) (Tanaka ) and SHANK1 (cg14560133) (Sungur ) are all autism spectrum disorder candidate genes, a phenotype tentatively associated with ART (Liu ). Although these observations are of potential interest, it should be noted that only for CpGs near GABRB3 the DNAm level correlated between blood and brain regions (Hannon ) and for none of the CpGs was DNAm level associated with expression of the nearest gene in whole blood data. In addition, several CpG dinucleotides can be linked to genes connected to the plethora of conditions that can be linked to subfertility. H6PD (cg27266479) is a candidate gene for polycystic ovary syndrome (Martínez-García ). INTS1 (cg03207674) has a crucial role in the developing blastocyst, as inhibition of INTS1 function causes growth arrest (Hata and Nakayama, 2007). NOS3 (cg00478390) knockout mice are used as an in vivo model of (recurrent) embryo loss, as nitric oxide metabolism plays an important role in implantation (Pallares and Gonzalez-Bulnes, 2010). Although the function of the genes is arguably plausible in the context of ART, we could not detect a correlation between DNAm level in blood and other tissues at these CpGs, although the available reference data does not have tissues relevant to subfertility (e.g. ovary, placenta and uterus).

The latter may hint that the various medical reasons for ART may underlie the associations, rather than the ART process itself, which is a key question in the study of ART (Luke ). We employed DAG (Krieger and Davey Smith, 2016) and DAG indicated that ascertainment of a direct effect is possible with a small minimal adjustment set (Supplementary Figs S6 and S7). There was little to no effect on the effect estimates from adjustment for various maternal characteristics, gestational age and the number of years until pregnancy (as a proxy for in-/subfertility). In addition, we did not find any evidence for different or stronger effects of ICSI on DNAm as has been hypothesized (Loke and Craig, 2016) and the fertility clinics from which we recruited all used the same culture media for ICSI and IVF. These last two items might argue that not the (reasons for) infertility but the ART process itself may explain the associations.

Our analysis did not find widespread genome-wide differences like other genome-scale studies to date (El Hajj ; Novakovic ), but is consistent with a genome-scale study on only IVF using immunoprecipitation of methylated DNA with next-generation sequencing (Castillo-Fernandez ). We extensively compared individual results, as far as possible, and found little overlap. One explanation is that genome-scale DNAm studies of ART to date had a relatively low statistical power and hence a substantial false-negative rate. This is also true for our own study, despite being one of the larger studies to date. The field would benefit from a future meta-analysis using ART case-control and regular (birth) cohort studies. An alternative explanation is that the differences may also stem from the different scope of the measurement techniques used (Castillo-Fernandez ) and the different sources of genomic DNA (Guthrie card versus whole cord blood) (Novakovic ). In addition, it is possible that different media were used between countries or that our results are different because only a small proportion of our IVF and ICSI groups consisted of embryos that were frozen and/or cultured extensively in comparison to other studies. We were underpowered to test the influence of this latter aspect as most ART procedures in our study entailed fresh embryos that were not cultured extensively in vitro.

Our study is one of the larger genome-scale DNAm studies on ART to date and we are the first to study both IVF and ICSI on this scale. Nonetheless, our study has important limitations to consider. First, DNAm is one of the drivers of cell identity (Jaenisch and Bird, 2003) and we measured cord blood, which may be a tissue less relevant in relation to study outcomes although it may still mark processes in relevant tissues due to mitotic inheritance (Heijmans and Mill, 2012). We used the latest methods to impute the seven major cell types in cord blood (Gervin ) and in line with an earlier report (El Hajj ), these imputed cell types did not differ between infants conceived with or without ART and adjustment did not alter the effect estimates. However, most of these major cell types consist of subsets of more specialized cells (subsequently showing unique DNAm profiles at increasingly select loci) and their influence could not be investigated. However, none of the 19 CpG dinucleotides is linked to a gene with a role in (auto)immune function or hematopoiesis. Another major driver of DNAm variation is genetic variation (Bonder ) and we did not control for genetic variation in our analyses. DNAm of the 19 CpG dinucleotides was not measured by probes where genetic variation influences the actual measurement of DNAm (Zhou ). A large meta-analysis has identified SNPs influencing the DNAm levels (Bonder ) and 10 out of the 19 CpG dinucleotides have one or more SNPs shown to influence DNAm levels. Although the variance explained by these SNPs was small, we cannot completely exclude the possibility that (part of) the difference in DNAm is explained by genetic variation as we have not measured this. However, prior work on the influence of the prenatal environment on DNAm variation shows that this effect of the prenatal environment can be completely independent of, and additive to, the influence of local genetic variation influencing DNAm at the same genomic regions (Tobi ). Also important is to consider the fact that an ART population is inherently different from the general population (Luke ). Despite controlling for possible confounding factors, and cross-checking our results with other EWAS on prenatal complications/exposures, unresolved confounding remains a possibility. Our study differs from some of the earlier ART studies, which is useful for triangulation between studies (Lawlor ), as our control group has a higher rather than a lower SES, as is normal in most of studies on ART.

The data presented here showcases modest and specific DNAm differences that are associated with ART, of which the functional relevance in adult tissues is unknown. We did not find any difference in DNAm patterns between IVF and ICSI. Our study found little evidence for the hypothesis that ART, be it IVF or ICSI, leads to widespread disturbances of DNAm patterns or for the hypothesis that ICSI has a different/larger relation with DNAm patterns. Our findings warrant cautious interpretation given the sample size and the subsequent low power due to the number of tests performed, the tissue studied and the unknown functional consequences of the identified DNAm differences.

Supplementary data

Supplementary data are available at Human Reproduction online.

Data availability

The datasets generated and analyzed during the current study are not publicly available due to Swedish privacy and data safety laws but are available from A.N.I. and C.A. on reasonable request and after meeting legal requirements.

Click here for additional data file.