Epigenetic And Transcriptomic Alterations in Offspring Born To Women With Type 1 Diabetes (The EPICOM Study)

doi:10.21203/rs.3.rs-1046258/v1

PDF

Research Article

Epigenetic And Transcriptomic Alterations in Offspring Born To Women With Type 1 Diabetes (The EPICOM Study)

https://doi.org/10.21203/rs.3.rs-1046258/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 18 Nov, 2021

You are reading this latest preprint version

Background: Offspring born to women with pregestational type 1 diabetes (T1DM) are exposed to an intrauterine hyperglycemic milieu and has an increased risk of metabolic disease in later in life. In this present study we hypothesize that in utero exposure to T1DM alters offspring DNA methylation and gene expression, thereby altering their risk of future disease.

Design: Follow-up study using data from the Epigenetic, Genetic and Environmental Effects on Growth, Metabolism and Cognitive Functions in Offspring of Women with Type 1 Diabetes (EPICOM) collected between 2012-2013.

Setting: Exploratory sub study using data from the nationwide EPICOM study.

Participants: Adolescent offspring born to women with T1DM (n=20) and controls (n=20) matched on age, sex and postal code.

Main outcome measures: This study investigates DNA methylation using the 450K-Illumina Infinium assay® and RNA expression (RNA sequencing) of leucocytes from peripheral blood samples.

Results: We identified 9 hypermethylated and 5 hypomethylated positions (p < 0.005, |DM-value| > 1). RNA expression profiling identified 38 up- and 1 down-regulated genes (p < 0.005, log2FC ≥ 0.3.). Functional enrichment analysis revealed enrichment in ontologies related to diabetes, carbohydrate and glucose metabolism, pathways including MAPK1/MAPK3 and MAPK family signaling and genes related to T1DM, obesity and atherosclerosis. Lastly, by integrating the DNA methylation and RNA expression data, we identified six genes where relevant methylation changes corresponded with RNA expression (CIITA, TPM1, PXN, ST8SIA1, LIPA, DAXX).

Conclusions: Findings suggest the possibility for intrauterine hyperglycemia to impact later life methylation and gene expression, a profile that may be linked to the increased risk of metabolic disease.

Endocrinology & Metabolism

Epigenetics & Genomics

Type 1 Diabetes

Pregnancy

Epigenetics

Methylation

Transcriptome

The “Developmental Origins of Health And Disease” hypothesis describes how an adverse early life environment may lead to increased risk of adult disease ¹. Offspring born to women with pregestational type 1 diabetes (T1DM) are exposed to a hyperglycemic intrauterine milieu and recent studies have shown an increased risk of developing type 2 diabetes, cardiovascular disease and obesity ^2–4. The pathophysiological mechanisms behind this increased risk in the offspring are unclear, and whether a link to intrauterine T1DM exposure exists remains to be elucidated. Epigenetic modifications have been proposed as a possible mechanism by which maternal diabetes induces long-term effects ⁵. However, our knowledge on epigenetic programming of offspring born to women with T1DM is still limited.

Methylation of cytosine is the best-studied epigenetic mark, with the addition of a methyl group to the C5 position of the cytosine-ring, and generally occurs on cytosines followed by a guanine (CpG). Genomic regions with an overrepresentation of CpG sites are termed CpG islands and are found at 60-70% of human gene promoters ⁶. Hypermethylation of CpG islands is involved in transcriptional repression, thus, changes in DNA methylation adds level of regulation to the DNA molecule without changing the DNA sequence and via this mechanism DNA transcription can be altered ⁷. Even though CpG sites are more sparse in gene bodies than in promotors, studies have described CpG hypermethylation in gene bodies to positively correlate with RNA expression ⁸. Changes in DNA methylation is well known to be involved in cell differentiation and embryogenesis and is considered to be part of developmental biology ⁹. Besides CpG methylation, epigenetic alterations also include histone modifications and non-coding RNA ⁵.

Recent studies in cord blood and peripheral blood have linked especially maternal overweight but also, to a lesser degree, maternal gestational diabetes to offspring epigenetic alterations ^10–12. Few studies have evaluated the impact of maternal T1DM on offspring epigenetic changes, and included studies exploring an adult cohort including both offspring born to women with gestational diabetes and T1DM, using muscle, adipose tissue or preadipocytes, described increased expression of microRNA-15a and microRNA-15b in skeletal muscle tissue and decreased leptin promotor methylation in subcutaneous adipocytes in offspring of T1DM women ^12–16. Only one study by Gaultier et al. explored the more genome-wide DNA methylation profile using Illuminas Human Methylation 27 Beadchip in describing an association between intrauterine exposure to maternal T1DM and kidney function in the adult offspring ¹⁷.

In this study we hypothesized that in utero exposure to maternal T1DM alters global DNA methylation and RNA expression and thereby alters the offspring’s risk of later disease development.

Through medical records and the original Diabetes Association registry we retrieved information regarding pregnancy and birth for the majority of both Index and Controls (Table 1). The Index and Controls included in the exploratory cohort were similar on most parameters including BMI, oral glucose tolerance test (OGTT) fasting plasma glucose, OGTT 120 min plasma glucose, fasting insulin, HOMA-IR, HbA_1c, systolic and diastolic blood pressure and birth weight, although significant differences were present in some of these parameters in the larger EPICOM group ¹⁸. Only gestational age at birth was significantly different for the intrauterine T1DM exposed children compared to the controls (at that time it was common practice to induce labor ~2-3 weeks prior to the expected delivery date in women with diabetes) (Gestational age: 259 (range=217–278) days vs. 278 (258–301) days, p <0.001). Also, maternal age at birth, maternal BMI at birth and parity were similar between the women with T1DM and the mothers of the Controls. For most of the women with T1DM we also had knowledge of their glycaemic regulation pregestationally and during the pregnancy. For the RNAseq data one sample from the control cohort was missing. This however, did not change the overall phenotypic results.

Table 1

Baseline characteristics of diabetes exposed children and controls
	Exposed to T1DM	Control children	P-value
Number of children	20	20
Male sex	10 (50%)	10 (50%)
Age (years)	17.8 (0.8)	18.1 (0.8)	0.26
Body mass index (kg/m²)	23.2 (2.6)	23.4 (2.2)	0.81
Oral glucose tolerance test (mmol/L)
time = 0 min	5.4 (0.4)	5.3 (0.3)	0.60
time = 120 min	6.0 (0.9)	6.4 (0.8)	0.20
S-Insulin (fasting) (pmol/L)	55.8 (17.3)	57.9 (19.8)	0.72
HOMA-IR	2.2 (0.7)	2.3 (0.8)	0.77
Waist circumference (cm)	75.8 (5.6)	76.4 (6.5)	0.77
Hip circumference (cm)	98.5 (7.4)	98.7 (5.3)	0.93
Systolic blood pressure (mmHg)	120.8 (9.6)	120.8 (11.6)	0.99
Diastolic blood pressure (mmHg)	63.4 (8.2)	63.7 (6.8)	0.89
Total cholesterol (mmol/L)	3.9 (0.7)	4.1 (1.1)	0.55
LDL cholesterol (mmol/L)	2.2 (0.6)	2.3 (1.0)	0.55
HDL cholesterol (mmol/L)	1.4 (0.4)	1.3 (0.4)	0.53
Triglycerides (mmol/L)^‡	0.8 (0.3-2.5)	0.9 (0.5-1.9)	0.24
HbA_1c (%) [mmol/mol]	5.2 (0.2) [33 (2.5)]	5.1 (0.3) [32 (2.8)]	0.15
Gestational age at birth (days)	259 (217–278) (n=18)	278 (258 – 301) (n=17)	< 0.001*
Birth weight (g)	3,482 (711) (n=19)	3,571 (542) (n=18)	0.68
Maternal Characteristics
Maternal age at birth (years)	28.6 (3.8) (n=20)	29.8 (3.7) (n=20)	0.30
Maternal BMI at birth (kg/m²)	22.9 (2.3) (n=16)	24.2 (4.4) (n=17)	0.47*
Parity ^#	1.5 (1–3) (n=20)	1.8 (1–4) (n=18)	0.55*
Pre-eclampsia	5 (26%) (n=19)	2 (12%) (n=17)	0.41^†
Maternal HbA_1c (%) [mmol/mol]
Pregestational	7.4 (1.3) [57] (n=18)	(–)
1st trimester	7.2 (1.2) [55] (n=17)	(–)
2nd trimester	6.3 (1.0) [45] (n=19)	(–)
3rd trimester	6.5 (0.9) [48] (n=18)	(–)

Differentially methylated positions and regions

Performing whole-blood genome-wide DNA methylation analysis using the 450K-Illumina Infinium assay we identified 14 differentially methylated positions (DMPs) with a p < 0.005 and |ΔM-value| > 1 (9 hypomethylated and 5 hypermethylated) between Index and Controls (Figure 1 and Supplemental Table S1 online) However, after Bonferroni correction, no single DMP reached the threshold of a family wise error rate (FWER) (< 0.05). Differentially methylated regions (DMRs) defined as methylation in groups of nearby positions have been proposed to be involved in transcriptional regulation. Therefore, we extended our analysis to the regional level. However, no DMRs were identified applying a Benjamini-Hochberg false detection rate (FDR) < 0.05.

RNA expression

Subsequently, we analysed gene expression by RNAseq to test if any differentially expressed genes (DEGs), coding or non-coding RNA, were present when comparing our two cohort groups (Index vs. Controls). We identified 39 coding DEGs (38 up regulated and 1 down regulated) (Table 2, Figure 2a+b and Supplemental Figure S1 online) (p < 0.005 and a log2FC ≥ 0.3). In addition, we also found 20 non-coding DEGs ((7 up regulated tags and 13 downregulated tags) (p < 0.005 and a log2FC ≥ 0.3) (Supplemental Table S2 online). No DEGs reached an FDR < 0.05.

Table 2

List of differentially expressed autosomal coding genes between type 1 diabetes exposed offspring and controls.
Gene-ID	Log2Fold Change	Log CPM	P-value	Chromosome Name	Start Position	End Position	External Gene Name
ENSG00000103148.11	1.402	4.682	0.0003637	16	134273	188859	NPRL3
ENSG00000036448.5	1.238	3.216	0.001291	8	1993155	2113475	MYOM2
ENSG00000163735.6	0.778	3.821	0.0003347	4	74861359	74864496	CXCL5
ENSG00000120885.15	0.754	3.931	7.9319E-05	8	27454434	27472548	CLU
ENSG00000259207.3	0.729	5.136	0.001563	17	45331212	45421658	ITGB3
ENSG00000005961.13	0.696	3.951	0.003207	17	42449548	42466873	ITGA2B
ENSG00000184922.9	0.687	5.329	0.004019	17	43298811	43324687	FMNL1
ENSG00000161911.7	0.682	1.937	0.001432	6	41117080	41122075	TREML1
ENSG00000113140.6	0.670	4.702	0.0003625	5	151040657	151066726	SPARC
ENSG00000175746.4	0.670	2.299	0.003372	15	39542885	39547046	C15orf54
ENSG00000138798.7	0.650	3.438	0.0008035	4	110834040	110933422	EGF
ENSG00000122786.15	0.642	2.518	0.004715	7	134429003	134655479	CALD1
ENSG00000100385.9	0.610	3.173	0.001113	22	37521878	37571094	IL2RB
ENSG00000082781.7	0.566	3.489	0.002852	3	124480795	124620265	ITGB5
ENSG00000198478.6	0.565	6.232	0.0013672	6	80341000	80413372	SH3BGRL2
ENSG00000168497.4	0.564	6.971	0.001083	2	192699028	192711981	SDPR
ENSG00000151693.5	0.564	2.812	0.0004620	2	9346894	9545812	ASAP2
ENSG00000185909.10	0.549	2.620	0.001593	3	49209044	49213917	KLHDC8B
ENSG00000140022.5	0.547	4.040	4.5483E-05	14	81727000	81902809	STON2
ENSG00000019582.10	0.545	6.485	0.003180	5	149781200	149792492	CD74
ENSG00000101162.3	0.516	7.425	0.001578	20	57594309	57601709	TUBB1
ENSG00000090975.8	0.506	2.549	0.002654	12	123468027	123634562	PITPNM2
ENSG00000095303.10	0.505	5.532	0.0005005	9	125132824	125157982	PTGS1
ENSG00000111644.3	0.497	2.010	0.001323	12	6747241	6756626	ACRBP
ENSG00000169313.9	0.483	3.076	0.003989	3	151055168	151102600	P2RY12
ENSG00000163737.3	0.471	4.875	0.004818	4	74846794	74847841	PF4
ENSG00000141522.7	0.459	3.114	0.003189	17	79825597	79829282	ARHGDIA
ENSG00000054793.9	0.441	3.264	0.002419	20	50213053	50385173	ATP9A
ENSG00000107438.4	0.440	3.137	0.004217	10	96997329	97050781	PDLIM1
ENSG00000143368.9	0.429	2.080	0.004332	1	149895209	149900236	SF3B4
ENSG00000153071.10	0.419	4.220	0.0003941	5	39371780	39462402	DAB2
ENSG00000067225.13	0.401	6.571	0.0002313	15	72491370	72524164	PKM
ENSG00000143537.9	0.399	2.580	0.001539	1	155023042	155035252	ADAM15
ENSG00000146192.10	0.387	5.061	0.003886	6	36973422	36996846	FGD2
ENSG00000177830.13	0.361	2.254	0.002612	11	867357	915058	CHID1
ENSG00000157823.12	0.339	3.560	0.004107	15	90373831	90437574	AP3S2
ENSG00000075624.9	0.335	9.389	0.003498	7	5566782	5603415	ACTB
ENSG00000172757.8	0.309	6.236	0.000952	11	65590493	65629497	CFL1
ENSG00000185261.9	-0.383	5.944	0.002067	5	93488671	93954309	KIAA0825

Functional enrichment analysis

To gain mechanistic insight into functional importance of the coding DEGs found above, we used enrichment analysis to associate these changes to disease phenotypes (DO), Biological Processes (GOBP) and biological pathways (REACTOME). We observed enrichment in ontologies relating to acquired metabolic disease, glucose metabolism disease, diabetes mellitus, carbohydrate metabolism disease, disease of metabolism and clustering around genes including CXCL5, EGF and SPARC (DO) (Figure 3). For biological processes and pathways, the analysis revealed enrichment relating to striated muscle tissue development, muscle tissue development, growth and developmental growth (GOBP) and biological pathways (REACTOME) including ECM proteoglycans, RAF/MAP kinase cascade, MAPK1/MAPK3 signaling and MAPK family signaling (Figure 3 and Supplemental Figure S2 online (cnet plot of GOBP)).

Functional enrichment analysis was also carried out with identified DMPs as input. Again, enrichment in in ontologies and pathways relating to carbohydrate metabolism and glucose transportation and pathways including Glucagon-like Peptide-1 (GLP1) regulated insulin secretion were in the top-enriched terms (Supplemental Table S3 and S4 online).

Correlation between DNA methylation and gene expression

DNA methylation has been shown to play a regulatory role in mediating gene expression, and if located within a gene promotor, higher levels of methylation usually repress transcriptional expression ⁸. To further explore the consequences of altered DNA methylation patterns, we therefore compared the DNA methylation data with our gene expression data (Figure 4 and Supplemental Table S5 online). First, we visualized correlations between changes at the group level. That is, the mean change in methylation level plotted against the mean change in gene expression of the gene annotated to the methylation site. DEGs overlapping with DMPs (annotated to closest gene to the methylation site), between T1DM exposed offspring and controls (p-value < 0.05) were plotted in a scatter plot based on log2FC and ΔM-value.

This analysis revealed that 93 DMPs annotated to 76 DEGs showed an inverse correlation where DNA hypomethylation correlated to an increased RNA expression (Figure 4, lower-right (blue)). On the contrary, DNA hypermethylation correlated to a decreased RNA expression, 14 DMPs annotated to 12 DEGs were found to be inversely correlated (Figure 4, upper-left (red)). Also, 201 DMPs annotated to 129 genes showed a positive correlation between DNA hypermethylation and increased RNA expression. As the second step of our correlation analysis, we performed a Spearman rank correlation for each DMP-DEG pair, using each subject’s methylation level at the involved DMPs and expression level (normalized counts) for the corresponding RNA expression data as input. Here, our analysis showed that 2 DMPs and their corresponding genes had significant inverse correlation between DNA hypomethylation and increased RNA expression (p < 0.1) (CIITA and TPM1, Table 3). DNA hypermethylation of CpG sites localized within the gene body has been described as being associated with increased RNA expression and in our analysis, we found 4 annotated DMPs localized within the gene body and their corresponding genes displaying positive correlation between DNA methylation and RNA expression (p < 0.1) (PXN, ST8SIA1, LIPA and DAXX, Table 3).

Table 3

Top candidate genes from the two-step correlation analysis.
	Probe name	Gene region	Relation to CPG island	Gene symbol	rho	P-value	FDR-P-value
DNA hypomethylation and positive RNA expression	cg06871213	Body	Nshore	CIITA	-0.3490	0.03002	0.06005
DNA hypomethylation and positive RNA expression	cg04194852	TSS1500	Island	TPM1	-0.3237	0.04493	0.42438
DNA hypermethylation and positive RNA expression (located in genebodies)	cg04256697	Body;TSS1500;Body	OpenSea	PXN	0.3933	0.01379	0.05514
	cg24112692	Body	OpenSea	ST8SIA1	0.3686	0.02152	0.26539
	cg13931663	Body;Body	OpenSea	LIPA	0.3020	0.06206	0.56834
	cg20618109	Body;Body;TSS1500;Body;Body	N_Shore	DAXX	0.2909	0.07271	0.581641

This exploratory study suggests the existence of epigenetic marks that possibly influence the long-term health of offspring born to women with pregestational T1DM and elucidate intuitive pathways to be involved. Combining these epigenetic marks with gene expression data in offspring exposed to T1DM in utero is novel and provides new insight into the mechanistic consequences of epigenetic changes. However, none of the findings were present after correction for multiple testing.

In a cohort of prospectively studied offspring of women with T1DM we identified 14 DMPs (p < 0.005 and |ΔM-value| > 1), 9 being hypomethylated and 5 hypermethylated, between T1DM exposed offspring and their matched unexposed controls. We also sought to identify any overall differences in RNA expression between index and controls, and were able to identify a set of predominantly up regulated genes (p < 0.005 and a log2FC ≥ 0.3).

The functional enrichment analysis showed clustering around the genes CXCL5, EGF and SPARC which is associated with diabetes, obesity and insulin secretion and enriched in disease ontology terms of metabolic disease, diabetes and glucose metabolism disease. These findings seem to fit the hypothesis of altered intrauterine programming of offspring born to women with T1DM ^19–21. The functional enrichment analysis points towards muscle tissue development and developmental growth to be associated with intrauterine exposure to T1DM and that the MAPK signaling pathway could play a role in the pathogenesis. The link between intrauterine exposure to T1DM and later life alterations in skeletal muscle has previously been described by Houshmand-Oeregaard et al. and even though no studies have described a link between intrauterine exposure to diabetes and later life alterations in the MAPK pathway, the central placement of MAPK in insulin signaling makes it a plausible object for pathogenetic changes ^13,22.

Our access to both DNA methylation and gene expression data made it possible to explore the impact of epigenetic changes induced by intrauterine T1DM exposure on gene expression. We chose to perform a two-step analysis where we first explored the overall correlation between DMPs and DEGs at group level. During step two, we studied the individual correlation between ΔM-value at the involved DMP and log2FC for the corresponding DEG. No studies within the area of fetal programming as a consequence of maternal T1DM have to our knowledge, performed a similar linkage study between DNA methylation and RNA expression data. However, in cancer studies similar methodology has been used ²³. Applying this method, provided us with six genes where changes in DNA methylation between offspring exposed to maternal T1DM and controls correlated with relevant changes in gene expression (Table 3). PXN, which is localized at chromosome 12 has been linked to progression of T1DM and LIPA at chromosome 10 to the composition of lipoproteins ^24,25. Mutations in TPM1 at chromosome 15 is associated with both the development of the metabolic syndrome and cardiac hypertrophy ^26,27. CIITA is associated with MHC class II (MHCII) expression and MHCII antigen presentation in adipocytes and is reported to trigger early adipose inflammation and insulin resistance as well as inducing changes to energy expenditure in skeletal muscle ^28,29.

In our study we chose to match our participants so that they did not overtly differ in phenotypic appearance. This approach resulted in a very homogenous group where the major difference was in T1DM exposure status. The matching was an attempt to exclude any phenotypic differences to be the cause of differences in methylation status. However, a consequence of the matching could, in theory, result in the observed changes in methylation and gene expression as being protective towards the development of later life metabolic disease, as none of the participants in the exploratory cohort had overt metabolic disease in contrast to the participants in the full EPICOM study ¹⁸. Also, it is well known that gestational age is associated with DNA-methylation and as gestational age differed between our two cohorts, this could have influenced our findings ³⁰.

Epigenetic changes in relation to intrauterine exposure to either nutrition or maternal diabetes has previously been described in animal studies but the results have been difficult to apply to human studies ^31,32. In recent years methylation analyses of cord blood have shown promising results, but the use of cord blood methylation status as a biomarker of long-term offspring risk of disease development has yet to be fully described ^11,33.

Here we used leucocytes from peripheral blood from both T1DM exposed and non-exposed offspring. The most likely target tissue for the pathogenesis of long-term consequences of being born to a mother with T1DM may, however, be pancreatic β-cells, and muscle- and fat-tissue. These three tissues are not as readily available and it is debated whether peripheral blood leucocytes reflects global epigenetic changes induced by maternal pregestational T1DM. Studies of tissue-specific differentially methylated regions have shown conservation of DNA methylation patterns across different tissues ^34,35. However, replication of our findings in target tissues would further strengthen our hypothesis. In our study we did not have access to specific white cell counts for each participant. Instead, we used the EstimateCellCounts algorithm from the minfi package, which enables us to minimize confounding related to cell composition, a strength compared to prior studies within the same area.

The samples used in our study were collected at a mean age of ~18 years. This provides time to develop phenotypic characteristics. However, it also provides the possibility that any epigenetic changes have occurred during the time period from birth to follow-up. In the EPICOM study we didn’t have access to cord blood samples and it is therefore not possible for us to state, that the epigenetic findings in our study have been present since birth. Lifestyle and diet in families where the mother has T1DM could be different compared to families without diabetes and this could affect the epigenetic findings ^36,37.

When interpreting our finding one must consider the fact that no single DMP reached the FWER adjusted value of 0.05. This is probably a consequence of our modest sample size in the exploratory cohort and an obvious limitation of our study. Instead we used the same pragmatic threshold for DMPs as West et al. providing an opportunity to compare findings in related studies ³⁸.

The strength of this study is that it was performed on a group of prospectively studied offspring born to women with T1DM. For both diabetes exposed offspring and their matched controls, we have extensive knowledge of intrauterine exposure, birth, neonatal period as well as current health status. Our women with T1DM are well characterized, and no women with T2DM or gestational diabetes were included in our cohort. This enables us to explore associations between both fetal and long-term consequences of being born to mother with T1DM and epigenetic and transcriptomic alterations.

In conclusion, our data indicate intrauterine hyperglycemia to impact later life methylation and gene expression, a profile that may be linked to the increased risk of metabolic disease. However comprehensive follow-up studies using a genome-wide approach and including relevant target tissue is needed.

The EPICOM cohort:

This study is a part of the EPICOM (Environmental Versus Genetic and Epigenetic Influences on Growth, Metabolism and Cognitive Functions in Offspring of Women with Type 1 Diabetes) study (ClinicalTrials.gov registration no. NCT01559181). From 1992-1999, pregnancies in women with pregestational T1DM were prospectively reported to a registry managed by the Danish Diabetes Association. Information regarding diabetes status and pregnancy outcome was reported to the registry by local obstetricians at eight hospitals in Denmark, who were responsible for antenatal care and delivery for pregnant women with T1DM. Coverage of cases from the reporting centers spanned from 75 to 93%, as described by Jensen et al. ³⁹. The Danish Diabetes Association Register consists of 1326 records of children born to women with T1DM. Only children born after 24 completed gestational weeks were registered. As part of the EPICOM study we invited one child per mother (the oldest) to participate in a clinical study concerning metabolic risk. Controls were recruited from the background population and were matched by sex, age and postal number. Details from this study have previously been described by Vlachova et al. ¹⁸.

Sample inclusion

Among the oldest participants of the EPICOM cohort we selected 20 offspring of women with pregestational T1DM (Index children) and their 20 matched controls (Controls) as our exploratory cohort of DNA methylation and RNA expression (Table 1). Besides fulfilling the matching criteria for inclusion in the EPICOM cohort, we made sure that our exploratory cohort did not differ in body mass index (BMI), glucose tolerance, HbA_1c, cholesterol or triglyceride levels and that the index child and their control were age matched. DNA methylation and RNA expression profiling were performed on leucocytes isolated from peripheral blood samples.

DNA isolation and the 450K-Illumina Infinium assay®

EDTA treated peripheral blood samples were collected from the participants when they participated in the EPICOM study and were stored as EDTA treated at -80°C until use. Genomic DNA was extracted from peripheral blood using QIAmp® Mini Kit (Qiagen, Germany). For each sample, 1 µg of genomic DNA was bisulfite-converted using Zymo EZ DNA Methylation Kit according to the manufacturer’s recommendations. DNA methylation level was measured using the 450K-Illumina Infinium assay (Illumina, Inc.) at Aros Applied Biotechnology A/S.

(Pre-)processing of the 450K-Illumina Infinium assay data

All analyses were performed in R statistics, version 3.6.1 and R package minfi (version 1.32.0) was used for normalization, analysis and visualization ⁴⁰. Detection p-values were calculated to identify failed positions with a p-value cut-off > 0.01. Probes were removed if they failed in more than 20% of the samples (n = 350). No samples were identified as failed, as the proportion of failed probes did not exceed 1% for any single sample. Density-bean-plots were used to identify outliers and minfi’s inbuilt function was used to evaluate data with respect to extreme methylation outliers (> 3 SD away from the median). We performed background normalization and control normalization implementing the pre-processing choices of Genome Studio®. Next, we applied subset-quantile-within-array-normalization correcting for technical differences between Infinium type I and II assay design allowing both within-array and between-sample normalization. Cross reactive probes (n = 29.541), probes with SNPs documented in C or G of the target (n = 18.284) probes on sex chromosomes (n = 12,312) were excluded, leaving 415,009 probes. Methylation values were calculated as M-values (logit [beta]) (Equation (I)) ⁴¹.

$$M-value=log2\left(\frac{Beta}{1-Beta}\right)$$

Multidimensional scaling plots were evaluated to identify clusters of samples behaving differently than expected. Finally, the probes were annotated to the human genome version 19 using the enhanced Illumina annotation method developed by Price et al. ⁴².

Estimated differential cell counts

We adjusted for differences in cell proportions using minfi’s estimateCellCounts that implements Houseman et al.’s regression calibration algorithm ⁴³. This algorithm returns relative proportions of CD4+ and CD8+ T-cells, natural killer cells, monocytes, granulocytes and B-cells in each sample.

Identifying differentially methylated positions

To identify positions where methylation is associated with intrauterine exposure to maternal T1DM we fitted a linear model, which utilizes a generalized least squares model (lmFit of R-package limma) allowing for missing values ⁴⁴. Significance was evaluated using F-test. The sample variances were estimated using an empirical Bayes approach with shrinkage towards the means. A Bonferroni adjusted FWER below 0.05 was considered significant. We applied the model both without and with adjustment for estimated relative cell proportions (CD4+ and CD8+ T-cells, natural killer cells, monocytes, granulocytes, and B-cells).

Identifying differentially methylated regions (DMRs)

We used DMRcate to identify any autosomal DMRs ⁴⁵. DMRcate identifies and ranks the most differentially methylated regions across the genome based on kernel smoothing of the differential methylation signal. The model performs well on small sample sizes and builds on the well-established Limma package, allowing us to incorporate estimated cell proportions as covariates. A FDR below 0.05 with a |ΔM-value| < 1 was considered significant.

RNA sequencing (RNA-seq) sample preparation

Blood samples were drawn using PAXgene Blood RNA Tubes® and placed 2 hours at room temperature, sequentially stored overnight at -20 degrees before being stored at -80 degrees until analysis.

RNA-Seq library construction and sequencing

Whole transcriptome, strand-specific RNA-Seq libraries were prepared from total-RNA using the Ribo-Zero Globin technology (Illumina, Inc.) for depletion of rRNA and globin mRNA followed by library preparation using the ScriptSeq technology (Illumina, Inc.). Depletion and library preparation were automated on a Sciclone NGS (Caliper, Perkin Elmer) liquid handling robot. The total-RNA (1.7 µg per sample) was subjected to Baseline-ZERO DNase prior to depletion. Total-RNA was purified using Agencourt RNAClean XP Beads before and after DNase treatment followed by on-chip electrophoresis on a LabChip GX (Caliper, Perkin Elmer) and by UV measurements on a NanoQuant (Tecan). Cytoplasmic and mitochondrial rRNA as well as globin mRNA were removed from 400 ng DNAse treated total RNA using the Ribo-Zero Globin Gold Kit (Human/Mouse/Rat, Illumina, Inc.) following the manufacturer’s instructions and quality of the depleted RNA was estimated on a LabChip GX (Caliper, Perkin Elmer). Synthesis of strand-specific RNA-Seq libraries were conducted using the ScriptSeq v2 kit (Illumina, Inc.) following the recommended procedure, and the qualities of the RNA-Seq libraries were estimated by on-chip electrophoresis (HS Chip, LabChip GX, Caliper, Perkin Elmer) of a 1 µL sample. The DNA concentrations of the libraries were estimated using the KAPA Library Quantification Kit (Kapa Biosystems). The RNA-Seq libraries were multiplexed paired-end sequenced on an NextSeq 500 (75 + 6 + 75 bp) using a high output flowcell.

RNA-Seq analysis

Paired de-multiplexed fastq files were generated using bclfastq (v.2.20 Illumina) and initial quality control was performed using FastQC. Adapter trimming was conducted using the GATK ReadAdaptorTrimmer tool followed by mapping to the human genome (hg19) in addition to transcripts from databases on non-coding RNAs (mirbase, mitranscriptome, rfam, snornabase, tjumirna and trnascanse) using Bowtie and then further analyzed using Tophat, and Cufflinks and HTSeq-count (union method) ^46,47. HTSeq-count (union method) was applied to produce raw counts which were then submitted for differential expression analysis in R using edgeR ⁴⁸. All non-informative features were filtered out by removing features with less than one counts per million (CPM) in 39 samples removing 98334 features, leaving 13842 for downstream analysis. A generalized linear model was fitted and p-values and log fold changes (Log2FC) were retrieved from the individual comparisons of Index vs. Controls. A Benjamini-Hochberg FDR below 0.05 was considered significant.

Functional enrichment analysis

To investigate possible biological functions of the RNA expression changes in our cohort, we performed gene set enrichment analysis (GSEA) using Clusterprofiler ⁴⁹. The input genes were ranked based on the magnitude of changes in the Index vs. Controls comparison (log2FC). Disease associations were identified by Disease Ontology (DO) enrichment analysis, while Gene Ontology Biological Processes (GOBP) and REACTOME were used for functional and pathway enrichment. The top-enriched terms for each enrichment analysis were sorted according to p-value and presented as barplots. Furthermore, gene-concept networks, that depicts linkages of the genes and biological terms (DO, GOBP and REACTOME terms), were made by evoking the cnetplot function of Clusterprofiler.

For methylation data, GOBP and REACTOME enrichment analysis were performed on all differentially methylated positions (p < 0.005) using the methylGSA package for R ⁵⁰. The probes were annotated to the human genome version 19 as previously described.

Correlation between DNA methylation and gene expression

DNA methylation has been shown to play an important role in modulating gene expression. Therefore, we correlated changes in methylation and gene expression of the gene closest to the methylation site. Shared differentially expressed genes (DEGs) and DMPs (closest gene to methylation site) between Index and Controls (p < 0.05) were plotted in a scatter plot based on log2FC and ΔM-value. In the second part of the analysis, we used Spearman’s rank correlation to analyze the correlation between the specific methylation level (M-value) and RNA expression level (normalized counts) for shared genes. We considered the correlations to be interesting (either negative or positive) if p < 0.1.

Analysis software

Statistics was done using R 3.6.1 (R Foundation for Statistical Computing, Vienna, Austria) with Bioconductor 3.9.⁵¹ DNA methylation data was analysed using the minfi ⁴⁰, DMRcate ⁴⁵ and Limma ⁵² package 1.32.0, and RNAseq data using edgeR ⁴⁸. Functional gene enrichment analysis was carried out using the R packages Clusterprofiler and MethylGSA ^49,50. Graphics were made using basic R functions ggbio, Gviz, DEseq, DEXSeq and ggplot2. The package knitr was used for data documentation.

Ethics declarations

All participants provided informed consent. The study was approved by The Danish Data Protection Agency and approved by The Danish National Committee on Health Research Ethics (M-20110239). All clinical investigations were conducted according to the Declaration of Helsinki. The research has been registered at ClinicalTrials.gov (no. NCT01559181).

Data availability

All sequence and methylation data from the EPICOM study are deposited at the EGA MOMA repository (https://ega-archive.org/dacs/EGAC00001000145) and will be available upon request within the General Data Protection Regulation.

Acknowledgements

The Danish Diabetes Association is acknowledged for originally assisting in the creation of a registry of pregnant women with T1DM. PD contributed to the establishment of the original registry. In addition, data collection in the original registry was performed by Lars Mølsted-Pedersen, Joachim Klebe, Niels Hahnemann, Margrethe Møller, Jes G. Westergaard, Hans Gjessing, Jens Kragh Mostrup, K. H. Frandsen, Edna Stage, Anders Thomsen, Thea Lousen, Kresten Rubeck Petersen, Bjarne Øvlisen, Jan Kvetny, Hedvig Poulsen. Apart from PD, the original registry working group included Henning Beck-Nielsen, Anders Frøland, Lars Mølsted-Pedersen, Joachim Klebe and Carl Erik Mogensen. Data collection during the EPICOM study was performed by Karen Mathiesen and Merethe Møller and Birgitte Trolle, Margaret Gellet and Pamela Celis performed bisulfite converting and DNA and RNA isolation.

JJ is acknowledged for generating Fig. 3 and 4 and Supplementary Fig. S1.

Author contributions

SK, ZL and BB contributed with data collection in the EPICOM cohort. SK, CH, ZL, BB, PD, KH, DMJ, KS, and CHG all contributed substantial to conception and design of the study. SK, AS, JJ, CT, SV and CHG analysed the data, and SK drafted the manuscript and designed the tables. All authors were involved in the interpretation of the data and critically revised the article. All authors had full access to the data in the study and take full responsibility for the integrity of the data and the accuracy of the data analysis. SK is the guarantor of the study and, as such, had full access to all data and affirms that the manuscript is honest, accurate, and transparent and no important aspects of the study have been omitted. The data has not been published previously and is not under consideration for publication elsewhere. The data has partially been presented at the 47^th Annual Diabetes and Pregnancy Study Group Meeting.

Additional Information

Grants and fellowships supporting the study: European Foundation for the Study of Diabetes, Danish Ministry of Science, Innovation and Higher Education, Lundbeck Foundation, Aarhus University, Danish Diabetes Academy, Beckett Foundation, the Novo Nordisk Foundation, Danielsen Foundation, Hede Nielsen Foundation and Central Denmark Region.

Competing interests statement: CHG has given talks for Novo Nordisk (NN). PD is participating in multicenter and multinational clinical studies on the use of insulin in pregnant women with pre-existing diabetes in collaboration with NN; no personal honorarium is involved. None of the study funders had any role in study design, collection, analysis, interpretation, preparation of the manuscript or in the decision to submit the article for publication.

1. Hochberg, Z. et al. Child health, developmental plasticity, and epigenetic programming. Endocr Rev 32, 159–224 (2011).

2. Clausen, T. D. et al. High Prevalence of Type 2 Diabetes and Pre-Diabetes in Adult Offspring of Women With Gestational Diabetes Mellitus or Type 1 Diabetes. Diabetes Care 31, 340–36 (2008).

3. Rijpert, M. et al. Risk Factors for Childhood Overweight in Offspring of Type 1 Diabetic Women With Adequate Glycemic Control During Pregnancy. Diabetes Care 32, 2099–2104 (2009).

4. Manderson, J. G. et al. Cardiovascular and metabolic abnormalities in the offspring of diabetic pregnancy. Diabetologia 45, 991–96 (2002).

5. Hjort, L. et al. Diabetes in pregnancy and epigenetic mechanisms—how the first 9 months from conception might affect the child’s epigenome and later risk of disease. Lancet Diabetes Endocrinol. 7, 796–806 (2019).

6. Illingworth, R. S. & Bird, A. P. CpG islands--’a rough guide’. FEBS Lett 583, 1713–1720 (2009).

7. Szyf, M. & Bick, J. DNA methylation: a mechanism for embedding early life experiences in the genome. Child Dev. 84, 49–57 (2013).

8. Jones, P. A. Functions of DNA methylation: Islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 13, 484–492 (2012).

9. Guo, H. et al. The DNA methylation landscape of human early embryos. Nature 511, 606–610 (2014).

10. Opsahl, J. O. et al. Epigenetic signatures associated with maternal body mass index or gestational weight gain: A systematic review. J. Dev. Orig. Health Dis. (2020). doi:10.1017/S2040174420000811

11. Jönsson, J. et al. Lifestyle intervention in pregnant women with obesity impacts cord blood DNA methylation which associates with body composition in the offspring. Diabetes 1–202 (2020).

12. Kelstrup, L., Hjort, L., Houshmand-Oeregaard, A., Clausen, T. D. & Ninna, S. Gene expression and DNA methylation of PPARGC1A in Muscle and Adipose Tissue from Adult Offspring of Women with Diabetes in Pregnancy. Diabetes 3, 1–41 (2016).

13. Houshmand-Oeregaard, A. et al. Increased expression of microRNA-15a and microRNA-15b in skeletal muscle from adult offspring of women with diabetes in pregnancy. Hum. Mol. Genet. 27, 1763–1771 (2018).

14. Hansen, N. S. et al. Fetal hyperglycemia changes human preadipocyte function in adult life. J. Clin. Endocrinol. Metab. 102, 1141–1150 (2017).

15. Houshmand-Oeregaard, A. et al. Differential adipokine DNA methylation and gene expression in subcutaneous adipose tissue from adult offspring of women with diabetes in pregnancy. Clin. Epigenetics 9, 1–12 (2017).

16. Houshmand-Oeregaard, A. et al. DNA methylation and gene expression of TXNIP in adult offspring of women with diabetes in pregnancy. PLoS One 12, 1–18 (2017).

17. Gautier, J. et al. Kidney Dysfunction in Adult Offspring Exposed In Utero to Type 1 Diabetes Is Associated with Alterations in Genome-Wide DNA Methylation. 1–14 (2015). doi:10.1371/journal.pone.0134654

18. Vlachová, Z. et al. Increased metabolic risk in adolescent offspring of mothers with type 1 diabetes: the EPICOM study. Diabetologia 58, 1454–1463 (2015).

19. Chavey, C. & Fajas, L. CXCL5 drives obesity to diabetes, and further. Aging (Albany. NY). 1, 674–677 (2009).

20. Maachi, H. et al. HB-EGF signaling is required for glucose-induced pancreatic β-cell proliferation in rats. Diabetes 69, 369–380 (2020).

21. Kos, K. & Wilding, J. P. H. SPARC: A key player in the pathologies associated with obesity and diabetes. Nat. Rev. Endocrinol. 6, 225–235 (2010).

22. King, G. L., Park, K. & Li, Q. Selective insulin resistance and the development of cardiovascular diseases in diabetes: The 2015 Edwin Bierman Award Lecture. Diabetes 65, 1462–1471 (2016).

23. Diéz-Villanueva, A. et al. DNA methylation events in transcription factors and gene expression changes in colon cancer. Epigenomics 12, 1593–1610 (2020).

24. Prashanth, G., Vastrad, B., Tengli, A., Vastrad, C. & Kotturshetti, I. Identification of hub genes related to the progression of type 1 diabetes by computational analysis. BMC Endocr. Disord. 21, 1–65 (2021).

25. Arnaboldi, L. et al. LIPA gene mutations affect the composition of lipoproteins: Enrichment in ACAT-derived cholesteryl esters. Atherosclerosis 297, 8–15 (2020).

26. Savill, S. A., Leitch, H. F., Harvey, J. N. & Thomas, T. H. Inflammatory adipokines decrease expression of two high molecular weight isoforms of tropomyosin similar to the change in type 2 diabetic patients. PLoS One 11, 1–13 (2016).

27. Mazzarotto, F. et al. Reevaluating the Genetic Contribution of Monogenic Dilated Cardiomyopathy. Circulation 387–398 (2020). doi:10.1161/CIRCULATIONAHA.119.037661

28. Zhang, S. Y. et al. Adrenomedullin 2 improves early obesity-induced adipose insulin resistance by inhibiting the class II MHC in adipocytes. Diabetes 65, 2342–2355 (2016).

29. Fang, M. et al. HDAC4 mediates IFN-γ induced disruption of energy expenditure-related gene expression by repressing SIRT1 transcription in skeletal muscle cells. Biochim. Biophys. Acta - Gene Regul. Mech. 1859, 294–305 (2016).

30. Merid, S. K. et al. Epigenome-wide meta-analysis of blood DNA methylation in newborns and children identifies numerous loci related to gestational age. Genome Med. 12, 1–17 (2020).

31. Waterland, R. & Jirtle, R. L. Transposable Elements: Targets for Early Nutritional Effects on Epigenetic Gene Regulation. Mol. Cell. Biol. 23, 5293–5300 (2003).

32. Li, C. C. Y. et al. Maternal obesity and diabetes induces latent metabolic defects and widespread epigenetic changes in isogenic mice. Epigenetics 8, 602–11 (2013).

33. Quilter, C. R. et al. Impact on offspring methylation patterns of maternal gestational diabetes mellitus and intrauterine growth restraint suggest common genes and pathways linked to subsequent type 2 diabetes risk. FASEB J. 28, 4868–79 (2014).

34. Slieker, R. C. et al. Identification and systematic annotation of tissue-specific differentially methylated regions using the Illumina 450k array. Epigenetics Chromatin 6, 26 (2013).

35. Byun, H. M. et al. Epigenetic profiling of somatic tissues from human autopsy specimens identifies tissue- and individual-specific DNA methylation patterns. Hum. Mol. Genet. 18, 4808–4817 (2009).

36. Fraga, M. F. et al. Epigenetic differences arise during the lifetime of monozygotic twins. Proc. Natl. Acad. Sci. U. S. A. 102, 10604–9 (2005).

37. Feil, R. & Fraga, M. F. Epigenetics and the environment: Emerging patterns and implications. Nat. Rev. Genet. 13, 97–109 (2012).

38. West, N. A., Kechris, K. & Dabelea, D. Exposure to Maternal Diabetes in Utero and DNA Methylation Patterns in the Offspring. Immunometabolism 3, 1–9 (2013).

39. Jensen, D. M. et al. Outcomes in type 1 diabetic pregnancies: a nationwide, population-based study. Diabetes Care 27, 2819–23 (2004).

40. Aryee, M. J. et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363–9 (2014).

41. Du, P. et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics 11, 587 (2010).

42. Price, M. E. et al. Additional annotation enhances potential for biologically-relevant analysis of the Illumina Infinium HumanMethylation450 BeadChip array. Epigenetics and Chromatin 6, 1–15 (2013).

43. Houseman, E. A. et al. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinformatics 13, (2012).

44. Smyth, G. K. Linear Models for Microarray Data. Bioinforma. Comput. Biol. Solut. using R Bioconductor, R. R. Gentlem, 397–420 (2005).

45. Peters, T. J. et al. De novo identification of differentially methylated regions in the human genome. 1–16 (2015).

46. Trapnell, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578 (2012).

47. Anders, S., Pyl, P. T. & Huber, W. HTSeq-A Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2015).

48. Robinson, M. D., Mccarthy, D. J. & Smyth, G. K. edgeR : a Bioconductor package for differential expression analysis of digital gene expression data. 26, 139–140 (2010).

49. Yu, G., Wang, L. G., Han, Y. & He, Q. Y. ClusterProfiler: An R package for comparing biological themes among gene clusters. Omi. A J. Integr. Biol. 16, 284–287 (2012).

50. Ren, X. & Kuan, P. F. methylGSA: a Bioconductor package and Shiny app for DNA methylation data length bias adjustment in gene set testing. Bioinformatics 35, 1958–1959 (2019).

51. Gentleman, R. C. et al. Bioconductor : open software development for computational biology and bioinformatics. Genome Biol. 5, (2004).

52. Gentleman, R., Carey, V., Huber, W., Irizarry, R. & Dudoit, S. Smyth, G.K. In. in Bioinformatics and Computational Biology Solutions using R and Biocunductor. 397–420 (2015).

Competing interest reported. CHG has given talks for Novo Nordisk (NN). PD is participating in multicenter and multinational clinical studies on the use of insulin in pregnant women with pre-existing diabetes in collaboration with NN; no personal honorarium is involved. None of the study funders had any role in study design, collection, analysis, interpretation, preparation of the manuscript or in the decision to submit the article for publication.

PDF

Version 1

posted 18 Nov, 2021

You are reading this latest preprint version

Epigenetic And Transcriptomic Alterations in Offspring Born To Women With Type 1 Diabetes (The EPICOM Study)

Status:

Version 1

Abstract

Figures

Introduction

Results

Differentially methylated positions and regions

RNA expression

Functional enrichment analysis

Correlation between DNA methylation and gene expression

Discussion

Methods

The EPICOM cohort:

Sample inclusion

DNA isolation and the 450K-Illumina Infinium assay®

(Pre-)processing of the 450K-Illumina Infinium assay data

Estimated differential cell counts

Identifying differentially methylated positions

Identifying differentially methylated regions (DMRs)

RNA sequencing (RNA-seq) sample preparation

RNA-Seq library construction and sequencing

RNA-Seq analysis

Functional enrichment analysis

Correlation between DNA methylation and gene expression

Analysis software

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1