Placental DNA methylation signatures of maternal smoking during pregnancy and potential impacts on fetal growth

Maternal smoking during pregnancy (MSDP) contributes to poor birth outcomes, in part through disrupted placental functions, which may be reflected in the placental epigenome. Here we present a meta-analysis of the associations between MSDP and placental DNA methylation (DNAm) and between DNAm and birth outcomes within the Pregnancy And Childhood Epigenetics (PACE) consortium (N = 1700, 344 with MSDP). We identify 443 CpGs that are associated with MSDP, of which 142 associated with birth outcomes, 40 associated with gene expression, and 13 CpGs are associated with all three. Only two CpGs have consistent associations from a prior meta-analysis of cord blood DNAm, demonstrating substantial tissue-specific responses to MSDP. The placental MSDP-associated CpGs are enriched for environmental response genes, growth-factor signaling, and inflammation, which play important roles in placental function. We demonstrate links between placental DNAm, MSDP and poor birth outcomes, which may better inform the mechanisms through which MSDP impacts placental function and fetal growth.

gestational age 126 127 Almost 1 in 10 pregnancies are impacted by the effects of maternal smoking during 128 pregnancy (MSDP), with state-specific prevalence ranging from as low as 1.8% to as 129 high as 27.1% in the USA 1 , while in Europe, the prevalence of MSDP is estimated to 130 range between 4.2% and 18.9% 2 . Consequently, the numerous health effects of MSDP 131 remain a significant public health concern. The impact of this exposure on fetal 132 development has been the source of significant investigation, resulting in MSDP being 133 recognized as a cause of multiple negative pregnancy and birth outcomes 3 . 134 The mechanisms that underlie this reproductive and developmental toxicity are 135 partially understood and include molecular and anatomical changes of the placenta 4,5 . These studies have begun to characterize the impact that MSDP has on the 153 human placental epigenome but have been limited by small sample sizes and have not 154 adjusted for placental tissue heterogeneity. We aimed to address these gaps by 155 performing a fixed effects meta-analysis examining the relationships between MSDP 156 and variations in the placental methylome across seven independent studies that are 157 members of the PACE consortium. We also aimed to gain insights into the potential 158 biological and health-related impacts of these associations by performing additional 159 analyses with nearby mRNA expression, as well as functional, regulatory, and 160 phenotypic enrichment, and a secondary meta-analysis of the associations between 161 placental DNAm and birth outcomes. This study represents the largest and most Glucose regulation in Gestation and Growth (Gen3G) 18 Table S1). 183

Genome-wide DNAm meta-analyses 184
We produced four statistical models for each CpG site, modeling the 185 associations between DNAm with both any and sustained MSDP, with and without 186 adjustment for putative cellular heterogeneity which was estimated with a reference free 187 deconvolution algorithm 22 . Models were adjusted for maternal age, parity, and maternal 188 education. Genomic inflation factors from the meta-analyses ranged from λ =2.0 to 2.8 189 (Excel Table S2, Supplemental Figure S1). Any MSDP was associated with DNAm at 190 532 CpG sites after Bonferroni adjustment (Excel Table S3), while sustained MSDP 191 was associated with DNAm at 568 CpG sites (Excel Table S4). After adjusting for 192 putative cellular heterogeneity, any MSDP was associated with 878 CpGs (Excel Table  193 S5), while sustained MSDP was associated with 894 CpGs (Excel Table S6). Among 194 these, 10.7% of models for any MSDP and 3.9% of models for sustained MSDP, not 195 adjusted for cell type, yielded heterogeneity p-values < 0.01, while only 5.5% of the 196 models for any MSDP and 3.5% for sustained MSDP, adjusted for putative cellular 197 heterogeneity, produced heterogeneity p-values < 0.01. Thus, heterogeneity in the 198 associations across cohorts was lowest in the analyses that were adjusted for putative 199 cellular heterogeneity and the CpGs from these models that yielded Bonferroni-200 significant associations (1224 unique CpGs) were carried forward for secondary 201 analyses. 202 Overall, the Bonferroni significant CpGs were distributed throughout the 203 genome, and the majority (68%) of CpGs exhibited lower DNAm in association with 204 any and sustained MSDP ( Figure 1). Of the 548 CpGs that yielded Bonferroni 205 significant associations with both any and sustained MSDP, the absolute values of the 206 parameter estimates were greater for the models of sustained MSDP compared to the 207 estimates for any MSDP (Excel Table S7). In a secondary analysis that restricted the 208 any MSDP models to the same three cohorts that contributed to the sustained MSDP 209 models (EDEN, GENEIDA, and INMA), we again found that almost all CpGs (544 out 210 of 548) exhibited increased magnitudes of association with increased duration of 211 exposure (Excel Table S7, Supplemental Figure S2). 212 The most notable association was observed at cg27402634, located upstream of 213 the LEKR1 gene and the non-coding RNA LINC00886, which showed the largest 214 differential DNAm and smallest p-values in all meta-analyses. Placentas that were 215 exposed to any MSDP had 22.95% lower DNAm (95% CI: 21.01-24.99% lower 216 DNAm; p-value = 5.35E-119) and those exposed to sustained MSDP had 25.78% lower 217 DNAm (95% CI: 24.02-27.54% lower DNAm; p-value = 7.22E-181) when compared to 218 mothers that did not smoke at all during pregnancy. Though all cohorts observed 219 substantial hypomethylation with MSDP at this CpG, the actual estimates of the 220 associations were highly variable between cohorts for models of any MSDP 221 (heterogeneity p-value = 2.66E-15), but relatively consistent for models of sustained 222 MSDP (heterogeneity p-value = 1.55E-01) ( Figure 2A). Overall, we consistency in the 223 associations across cohorts for the vast majority of the 1224 CpGs that yielded 224 Bonferroni-significant associations via meta-analysis: 92% and 97% of these CpGs 225 yielded heterogeneity p-values > 0.01, for their associations with any and sustained 226 MSDP respectively. In addition to cg27402634, we highlight those relationships that 227 yielded that largest magnitudes of association: (|β Any MSDP |) > 0.05 for cg26843110 228 (EDC3), cg20340720 (WBP1L), and cg17823829 (KDM5B) ( Figure 2B-2D). We 229 identified numerous other noteworthy relationships but due to the large number of 230 Bonferroni-significant associations from our meta-analyses, we highlight the 231 relationship among the 20 most statistically significant CpGs from the primary meta-232 analysis of any MSDP going forward (Table 2). 233

Expression Quantitative trait methylation (eQTM) analyses 234
We tested whether the DNAm levels at MSDP-associated CpGs were associated 235 with the expression of nearby (± 250kb of the CpG) mRNA from 194 placental samples 236 in the RICHS cohort, and corrected for multiple testing with the less conservative false 237 discovery rate (FDR) due to the smaller sample size. We identified 170 associations 238 between DNAm and gene expression within a 5% FDR (Excel Table S8), including 141 239 unique CpGs and 140 unique mRNAs; 72.4% of the eQTMs exhibited inverse 240 associations, while inverse associations were most common and statistical significance 241 was strongest for CpGs that were closer to the transcription start site (Supplemental 242 Figure S3). Among the top 20 CpGs from our meta-analysis, eight were identified as 243 eQTMs (5% FDR) (Table 3). Higher placental DNAm at cg27402634, was associated 244 with lower expression (β = -2.32, and FDR < 5%) of LEKR1. Additionally, DNAm at 245 cg17823829 (annotated to KDM5B) was most strongly associated lower expression of 246 PPFIA4 (FDR < 5%), and DNAm at cg26843110 (annotated to EDC3) was nominally 247 associated with higher expression of CSK, while cg20340720 (annotated to WBP1L) 248 was not associated with any nearby mRNA expression. 249

Functional and regulatory enrichment analyses 250
Enrichment analyses were performed to gain insights into which biological 251 processes may be impacted by these MSDP-associated CpGs. We performed gene-set 252 enrichment analyses, which included 565 genes annotated to CpGs associated with any 253 MSDP and 581 genes for sustained MSDP. 143 and 25 pathways were significantly (q-254 value < 0.05) enriched among CpGs associated with any MSDP (Excel Table S9) and 255 sustained MSDP (Excel Table S10), respectively. In both gene-sets, "signaling by nerve 256 growth factor (NGF)" was the most significant pathway (q-value = 2E-04 for any 257 MSDP and q-value = 9.6E-03 for sustained MSDP). Other significantly enriched 258 pathways were related to growth factors (VEGF, EGF, PDGF, IGF1R), hormones 259 (aldosterone, insulin, TSH, GnRH), interleukins (IL2, IL4, IL7), myometrial 260 contraction, vascular smooth muscle contraction, thrombin and platelet activation, 261 signaling and aggregation. We also tested whether the genes annotated to MSDP-262 associated CpGs were enriched for regulatory regions of specific transcription factors 263 (TFs). Most notably, our MSDP-associated CpGs were enriched for genes regulated by 264 GATA1 and GATA2, as well as by RUNX1, TP63, SMAD4, AR, TP53 or PPARG 265 (Excel Tables S11-S12). 266 We then examined whether the MSDP-associated CpG sites were enriched for 267 CpG island locations, allele-specific germline differentially methylated regions 268 (gDMR) 23 , regulatory features from the placenta specific 15-chromatin state annotation 269 from ROADMAP 24 , or placenta specific partially methylated domains (PMD) 25

Phenotype enrichment analyses 281
The genes annotated to our MSDP-associated CpGs were tested for phenotype 282 enrichment using data from the database of Genotypes and Phenotypes (dbGAP). Our 283 MSDP-associated genes were enriched for numerous phenotypes in dbGAP, including 284 several adiposity phenotypes (body mass index (BMI), waist-hip ratio, and obesity), 285 blood pressure, cardiovascular diseases, type 2 diabetes, asthma and respiratory function 286 (Excel Tables S13-S14). 287

Proximity to genetic variants for birth outcomes 288
We aimed to understand whether the MSDP-associated CpGs that we identified 289 co-localize within the same genomic regions as genetic variants that have been 290 associated with birth outcomes via genome-wide association studies (GWAS). Thus we 291 investigated whether MSDP-associated CpGs were within ± 0.5 Mb (1 Mb window) of 292 single nucleotide polymorphisms (SNPs) that have been associated with birth weight 293 (BW), birth length (BL), head circumference (HC) and gestational age (GA) 26-31 , which 294 have been added as annotations to the meta-analysis results files (Excel Tables S3-S6). 295 Of the 324 birth outcome SNPs in autosomal chromosomes, 94 SNPs (83 loci) and 108 296 (97 loci) were within 0.5 Mb of CpGs that were associated with any or sustained 297 MSDP, respectively (Excel Tables S15). Overall ~16% of the 1224 MSDP-associated 298 CpG sites were within 0.5 Mb of birth outcome SNPs, including cg27402634 (LEKR1), 299 cg26843110 (EDC3) and cg20340720 (WBP1L). We also explored whether our MSDP-300 associated CpGs may be biased by methylation quantitative trait loci (mQTLs), in 301 which SNPs influence the methylation levels at nearby CpGs. Two studies have 302 examined this question in human placenta, identifying 866 32 and 4,342 33 placental 303 mQTLs. Our findings did not appear to be biased by genetic variation as only 9 of the 304 1,224 MSDP-associated CpGs are previously characterized placental mQTLs. 305

Association of DNAM at smoking associated loci and smoking related birth outcomes 306
We also performed a second meta-analysis to examine the relationships between 307 DNAm with gestational age at birth, preterm birth, BW, BL, and HC z-scores, outcomes 308 that are known to be related to maternal smoking. Of the 1224 CpGs tested, 341 309 (27.9%) were related to at least one birth outcome after Bonferroni-adjustment 310 (0.05/1224). The majority of birth outcome associations were related to gestational age 311 at birth (298 CpGs) (Excel Table S16). Preterm delivery, for which only two cohorts 312 could contribute data (EDEN and NHBCS), produced similar associations found for 313 GA, though fewer CpGs were statistically significant (Excel Table S17). We also found 314 that numerous loci were associated birth size z-scores, with the majority of these being 315 associated with BW (43 CpGs) (Excel Table S18), followed by BL (20 CpGs) (Excel 316  Table S19) and HC (4 CpGs) (Excel Table S20). Some of the CpGs associated with GA 317 were also associated with birth size measurements, even though they were standardized 318 for GA, suggesting independent associations with both gestational duration and fetal 319 growth (Supplemental Figure S8A), including 6 CpGs (annotated to SYNJ2, PXN, 320 PTPRE, IGF2BP2, 4q21.1) shared between GA and BW, and 2 (annotated to POLR3E 321 and LOC441869) shared between GA and BL. Among the CpGs that yielded at least 322 one Bonferroni-significant association with birth outcomes, CpGs that tended to have 323 positive associations with birth outcomes clustered together and were typically 324 hypomethylated with MSDP, while CpGs that exhibited inverse associations with birth 325 outcomes tended to be hypermethylated with exposure to MSDP (Supplemental Figure  326 S8B). 327 Among our top 20 CpGs that were associated with any MSDP, 5 were associated 328 with GA at birth and 4 were associated with BW z-scores after Bonferroni adjustment 329 (Table 3). DNAm at cg27402634 (LEKR1) and cg20340720 (WBP1L), both located 330 close to BW-SNPs and for which MSDP associated with lower DNAm, were associated 331 with larger BW (p-value = 6.71E-07, and p-value = 2.42E-07). On the other hand, 332 DNAm at cg26843110 (EDC3; hypomethylated in response to MDSP and also close to 333 BW-SNPs) and at cg17823829 (KDM5B; hypermethylated) were associated with longer 334 and shorter gestational ages at birth, respectively (p-value = 5.09E-12, and p-value = 335 9.12E-06). Forest plots of BW z-scores and GA for these four CpGs are shown in 336 We summarize the results for all secondary analyses with a circos plot, for those 338 548 CpGs that yielded Bonferroni significant associations with both any and sustained 339 (FDR < 5%) and at least one birth outcome (Bonferroni-significant), which have been 341 annotated with the gene symbols from their respective eQTM models. 342

Comparison with smoking-sensitive CpGs in cord blood 343
We assessed whether the DNAm signatures of MSDP in the placenta were 344 consistent with MSDP associations in cord blood previously reported by the PACE 345 consortium 10 . Only nine CpGs annotated to seven unique genes (AHRR, CYP1A1, 346 GNG12, PXN, RNF122, SLC23A2, and ZBTB4) yielded Bonferroni-significant 347 associations in both tissues, out of 1224 CpGs from our study and 568 CpGs from the 348 cord blood study (Table 4). Of note, the CpGs within CYP1A1 and RNF122 showed 349 opposite directions of association with MSDP in cord blood and placenta. We also 350 compared the parameter estimates from our study that yielded associations with MSDP 351 within 5% FDR to the parameter estimates of those CpGs described in cord blood also 352 within a 5% FDR 10 . There was no overall correlation (r 2 <0.1) of the regression 353 coefficients across these two tissues (Supplemental Figure S9). 354 355 Discussion 356 We identified 1224 CpG sites with placental methylation levels that were 357 associated with any or sustained MSDP. Differential DNAm was greater with increased 358 duration of exposure at all loci that were associated with both any and sustained MSDP, 359 and a large proportion of the MSDP-associated CpGs were related to birth outcomes. 360 Those CpGs that were observed to have higher DNAm associated with MSDP, tended 361 to be inversely associated with gestational age and birth size, while CpGs exhibiting 362 lower DNAm with MSDP tended to be positively associated with gestational age and 363 birth size. 364 The MSDP-associated loci with the most statistically significant association 365 Additionally, three CpG sites within CYP1A1 and RNF122 were identified in both meta-454 analyses but with different directions of association in cord blood versus placenta. 455 Interestingly, we observed CYP1A1 to be hypomethylated in placenta with exposure to 456 MSDP, which is consistent with studies of adipose, skin, and lung tissues 54 , but this 457 CpG was hypermethylated in the cord blood meta-analysis 10 . Additionally, the most 458 statistically significant association with MSDP in placenta (cg27402634, LEKR1) was 459 not associated with MSDP in the cord blood meta-analysis. These observations, and the 460 lack of overall correlation in regression coefficients when comparing placental and cord 461 blood responses to MSDP, suggest that there are unique tissue-specific molecular 462 responses to this exposure. 463 The above findings should be interpreted within the context this study's 464 limitations. MSDP was self-reported and subject to misclassification, though differential 465 misclassification likely would have biased our findings towards the null. We modeled 466 two different definitions for MSDP and found that the models with sustained MSDP 467 produced larger magnitudes of association, as has been previously found in studies of 468 blood DNAm 55 . Although these findings suggest that increased duration of MSDP is 469 associated with greater differential DNAm, we did not assess dose-response patterns (ie. 470 number of cigarettes or cotinine concentrations), which should be the focus of future 471 investigations. Our study predominantly consisted of samples from mother-infant pairs 472 of European ancestry, and thus additional studies involving diverse racial and ethnic 473 backgrounds are needed in order to improve the generalizability of these findings. 474 While it is unlikely that placental DNAm would influence maternal smoking, the 475 observed associations between DNAm levels and reproductive outcomes could be due 476 to reverse causation. Placenta is a heterogeneous tissue with multiple different cell 477 types 56 that serve different functions and thus have different epigenetic states 57 . To 478 correct for this, we estimated and adjusted for variability in placental DNAm that may 479 be due to tissue heterogeneity. We utilized a data driven approach that was not based on 480 a methylome reference, as no references for placental cell-type methylomes are 481 currently available. Adjustments for these estimates of putative cellular admixtures did 482 reduce heterogeneity in the meta-analyses, thus improving the consistency of observed 483 associations between MSDP and DNAm across these independent cohorts. However, it 484 is possible that residual confounding may have influenced some of our results. 485 Despite these limitations, our study had numerous strengths, including a large 486 overall sample-size and seven independent studies to identify these relationships. We 487 used harmonized definitions of exposure variables and covariates, standardized 488 protocols for quality control and pre-processing of DNAm data, and standardized 489 methods for estimating/adjusting for tissue heterogeneity and for statistical analyses. 490 We also performed secondary analyses involving mRNA expression, functional and 491 phenotype enrichment, overlap with GWAS hits for reproductive outcomes, and meta-492 analyses of DNAm variation with birth outcomes to provide biological and health-493 related interpretations of our findings. 494 We identified a DNAm signature of MSDP in the placenta that shows substantial 495 differences from that observed in cord blood, most notably the CpG in close proximity 496 to LEKR1. Many of the identified MSDP-associated loci are involved in biological 497 process that are known to play critical roles in placental development, including 498 vascularization, angiogenesis, and inflammation. Additionally, a large proportion of the 499 MSDP-associated CpGs were also associated with GA at birth, or birth size z-scores, 500 suggesting these placental epigenetic variations may be intermediate molecular markers 501 between MSDP and these outcomes. Further study is required to determine whether 502 these epigenetic variations are causal mediators within these relationships or reflecting 503 other processes. Beta-values were normalized via functional normalization 58 and beta-mixture quantile 534 normalization (BMIQ) 59 was applied to correct for the probe type bias. Cohorts applied 535 ComBat to remove batch effects when applicable. Probes that hybridize to the X/Y 536 chromosomes, cross-hybridizing probes and probes with SNPs at the CpG site, 537 extension site, or within 10 bp of the extension site with an average minor allele 538 frequency > 0.01 were filtered out 60 . Overall, 418,639 probes and 415,396 were 539 available for modelling any MSDP and sustained MSDP, respectively. Finally, DNAm 540 extreme outliers (<25 th percentile -3*IQR or >75 th percentile + 3*IQR across all the 541 samples) were trimmed. 542

Estimates of putative cellular heterogeneity 543
Placental putative cellular heterogeneity was estimated from DNAm data using a 544 reference-free cell-mixture decomposition method 61 . The number of surrogate variables 545 ranged from 2 to 5 depending on the cohort. Models for differential DNAm were 546 corrected for the number of surrogate variables minus one to reduce multi-collinearity. 547

Genome-wide differential DNAM analyses 548
Within each cohort, robust linear regression from the "MASS" package 62 in R 549 were used to account for potential heteroskedasticity while testing the associations 550 between normalized DNAm beta values at each CpG with any MSDP and sustained 551 MSDP. Models were adjusted for maternal age, parity, maternal education and cohort-552 specific variables first unadjusted for putative cellular heterogeneity then adjusted for 553 cellular heterogeneity. We performed inverse variance-weighted fixed-effects meta-554 analyses using METAL 63 . The meta-analysis was performed independently by two 555 groups to ensure consistent results. CpGs not retained in at least 2 cohorts were filtered 556 out. We used the Bonferroni adjustment to control for multiple testing. To examine 557 whether increased duration of exposure (sustained smoking versus any smoking) 558 yielded increased magnitudes of association, we calculated the percent change in the 559 coefficients between the two models (|β Sustained | -|β Any |)/|β Any | * 100. Secondary analyses 560 were only performed on CpGs that yielded Bonferroni significant associations with any 561 or sustained MSDP in models that were adjusted for putative cellular heterogeneity. 562

Expression quantitative trait methylation (eQTM) loci 563
We performed expression quantitative trait methylation (eQTM) 64  ChIP-X database, and dbGaP database, respectively. EnrichR results were ranked using 588 the combined score (P-value computed using Fisher exact test combined with the z-589 score of the deviation from the expected rank). 69 Enrichment for regulatory features was 590 assessed with the hypergeometric test, and P-values were Bonferroni-corrected for 15 591 (placental chromatin 15-states) and 6 (relation to CpG island) tests, respectively. 592

Overlap of MSDP-sensitive CpG sites and birth outcome SNPs 593
Co-localization between MSDP-associated CpGs in placenta with previously 594 identified BW, BL, HC and GA SNPs from the largest genome-wide association studies 595 (GWAS) to date 26-31 was assessed using the GenomicRanges package in R 70 . We 596 identified which CpGs were located within 1 Mb windows (± 0.5 Mb) surrounding each 597 of the 324 autosomal SNPs, which correspond to 280 potential unique loci. Unique loci 598 were defined based on the criteria in Warrington et al. 2019 31 , and linkage 599 disequilibrium in Europeans (r 2 > 0.1 in < 2Mb). 600

Association between DNAm and birth outcomes 601
Within each cohort, robust linear regression models were utilized to test the 602 association between normalized DNAm beta values at each CpG as the independent 603 variable and gestational age at birth (inverse normal transformation of sex residuals), 604 BW z-scores, BL z-scores, and HC z-scores as the dependent variables. Logistic 605 regression was used to examine the relationships between DNAm and pre-term birth 606 (defined as <37 weeks of gestation). Birth size z-scores were calculated using 607 international references from the INTERGROWTH-21 st Project 71 and standardized by 608 both gestational age and newborn sex. Models were adjusted for maternal age, parity, 609 maternal education, cohort-specific variables (see Supplemental Methods) and putative 610 cellular heterogeneity. Inverse variance-weighted fixed-effects meta-analyses 63 were 611 again used to estimate pooled associations. Multiple testing was controlled with the 612 Bonferroni adjustment (α = 0.05/1224). 613

Comparison of MSDP-sensitive CpGs in placenta with cord blood 614
We examined the consistency between MSDP-sensitive CpGs in placenta and in 615 cord blood 10 . First we compared the coefficients from the models for sustained MSDP 616 in cord blood, unadjusted for cellular heterogeneity, to results for both any and 617 sustained MSDP in placenta, adjusted for cellular heterogeneity, using Pearson 618 correlation coefficients. 619 All DNAm data processing and analyses were conducted in R, with the 620 exception of the meta-analyses which were performed with METAL. 621 622 623  Table 2. Meta-analysis results from models of any and sustained MSDP, for the 20 CpGs with the smallest 830 p-values for the model of any MSDP, and the percent increase in effect size between sustained and any 831 MSDP (% Change); all models adjusted for maternal age, parity, education and putative cellular 832 heterogeneity; CpGs that were not annotated with a gene name in the Illumina 450K annotation file have  833 been annotated with their genomic region (ie. 4q12  Table 3. Results from eQTM models, DNAm versus gestational age at birth, and DNAm versus BW Z-score 839 models that yielded at least nominally significant associations, among the 20 CpGs that yielded the most 840 statistically significant associations with any MSDP in the primary meta-analysis (sorted in the same order 841 as   4: Differentially methylated CpGs with MSDP in placenta (from our meta-analysis) and cord blood 848 (from a published PACE meta-analysis 10 ) that yielded Bonferroni-significant associations in both tissues. Volcano and Manhattan plots of the association between any MSDP (total N=1700, exposed = 855 344) (A) and sustained MSDP (total N=795, exposed = 163) (B) with placental DNAm adjusted for maternal 856 age, parity, maternal education and putative cellular heterogeneity. In the volcano plots, the x-axis shows the 857 difference in DNAm (effect) with a possible range between 0 and 1, while the x-axis in the Manhattan plot 858 represents genomic location; both plots share the same y-axis with -log 10 (P). Bonferroni thresholds for 859 statistical significance are shown as blue dots and a blue horizontal line, for volcano and Manhattan plots, 860 respectively. The y-axes were truncated to a minimum p-value of 1*10 -40 (or maximum -log 10 (P) of 40), to 861 allow for better visualization of the majority of our results. The CpGs that were impacted by y-axis 862 truncation are indicated with arrows and annotated with their actual p-values. 863 864 865 866