Associations between infant sex and DNA methylation across umbilical cord blood, artery, and placenta samples

ABSTRACT

DNA methylation (DNAm) is vulnerable to dysregulation by environmental exposures during epigenetic reprogramming that occurs in embryogenesis. Sexual dimorphism in environmentally induced DNAm dysregulation has been identified and therefore it is important to understand sex-specific DNAm patterns. DNAm at several autosomal sites has been consistently associated with sex in cord blood and placental foetal tissues. However, there is limited research comparing sex-specific DNAm across tissues, particularly differentially methylated regions (DMRs). This study leverages DNAm data measured using the Illumina HumanMethylation450 BeadChip in cord blood (N = 179), placenta (N = 229), and umbilical artery samples (N = 229) in the PRogramming of Intergenerational Stress Mechanisms (PRISM) cohort to identify autosomal DMRs and differentially methylated positions (DMPs). A replication analyses was conducted in an independent cohort (GEO Accession GSE129841). We identified 183, 257, and 419 DMRs and 2119, 2281, and 3405 DMPs (p_Bonferroni < 0.05) in cord blood, placenta, and artery samples, respectively. Thirty-nine DMRs overlapped in all three tissues, overlapping with genes involved in spermatogenesis (NKAPL, PIWIL2 and AURKC) and X–inactivation (LRIF1). In replication analysis, 85% of DMRs overlapped with those identified in PRISM. Overall, DMRs and DMPs had higher methylation levels among females in cord blood and artery samples, but higher methylation levels among males in placenta samples. Further research is necessary to understand biological mechanisms that contribute to differences in sex-specific DNAm signatures across tissues, as well as to determine if sexual dimorphism in the epigenome impacts response to environmental stressors.

KEYWORDS:

• Epigenetics
• DNA methylation
• EWAS
• sex
• birth cohort

Introduction

DNA methylation (DNAm) represents the most studied epigenetic marker reflecting prior and current environmental exposures and future disease risks. The epigenome is particularly vulnerable to dysregulation by environmental exposures during embryogenesis, a time of rapid cellular differentiation and extensive epigenetic reprogramming [1]. Demethylation across functional elements occurs after fertilization, followed by remethylation after implantation [2]. As conceptualized in the developmental origins of health and disease hypothesis, conditions during development, including nutrition, stress, and environmental chemicals, can alter epigenetic marks resulting in health consequences throughout the lifespan [3]. For example, in epidemiological studies, in utero dietary factors and exposure to environmental pollutants have been associated with changes in DNAm of foetal tissues [4] which have been shown to persist into adulthood [5,6]. Methylation patterns differ between foetal tissues; embryos have an overall bimodal distribution of DNAm and the placenta exhibits a higher proportion of loci with intermediate DNAm levels [2,7], particularly within gene bodies [8].

Epidemiological studies document sexual dimorphism in environmentally induced disorders [9–12], DNAm dysregulation [13–15], and gestational epigenetic age acceleration measured in placenta [16]. Epigenetics represents a potential mechanism through which sexually dimorphic in utero programming effects manifest [17]. Foetal sex is increasingly recognized as a key modulator of placental development and adaptation to environmental conditions throughout pregnancy [14,17–20]. Sex differences in the methylome are well established; however, it is less understood how these differences translate across foetal tissues. DNAm is involved in X inactivation in females, a mechanism for dosage compensation [21]. In addition, epigenome-wide association studies (EWAS) have found sex differences in autosomal DNAm patterns of foetal tissues [22–24]. For example, an EWAS of cord blood (N = 111) measuring DNAm using the Illumina Infinium HumanMethylation450 BeadChip (450K) identified 3,031 differentially methylated positions (DMPs) (false discovery rate (FDR)-adjusted p < 0.05), with predominantly (83%) higher methylation among females [22]. Similarly, a study of cord blood (N = 71) identified 2,332 autosomal DMPs, 66% of which had higher methylation among females [25]. Conversely, an EWAS of placenta samples (N = 84) identified 91 DMPs (p < 1 × 10⁻⁷) including 82% with lower methylation in females compared to males [24]. This finding has been corroborated in a study of four placentas using whole-genome oxidative bisulphite sequencing in which Gong et al. observed global hypermethylation among males [8].

Although sex differences appear to occur throughout the autosomal methylome and sex-specific associations between environmental exposures and differentially methylated regions (DMRs) have been reported [26], fewer studies have reported sex-associated DMRs. Yousefi et al. identified 2,471 sex-associated DMRs annotated to autosomes in cord blood, which predominantly (83%) had higher mean methylation levels among females [22]. Maschietto et al. assessed DMRs by preidentified and tiling regions, identifying 171, 131, 127 and 333 DMRs in CpG islands, promotors, genes, and tiling regions, respectively, all of which had a larger proportion of higher methylation among females [25]. Gong et al. conducted DMR analysis of placenta DNAm by 5kb tiles, identifying a cluster of DMRs located on chromosome 8 near the gene CSMD1 with higher mean methylation levels among males (72% in males vs. 36% in females) [8].

Investigating sex differences in the methylome, particularly across tissues, is important in understanding the role of epigenetic regulation in determining phenotype, as well as in understanding sex-specific effects of environmental exposures on epigenetic dysregulation, particularly during development. Sex-specific regional methylation patterns may provide more biological insight than individual DMPs because DNAm measured at individual loci is subject to measurement error, DNAm is correlated over short distances [27], and CpG islands commonly occur in functionally relevant chromosomal regions including promoters [28]. Considering this motivation, our primary objective was to identify and compare DMRs associated with infant sex in cord blood, placenta, and umbilical artery samples by leveraging DNAm data measured using the Illumina Infinium HumanMethylation450 BeadChip as part of the PRogramming of Intergenerational Stress Mechanisms (PRISM) pregnancy cohort. These tissues have distinct roles in foetal development which may influence epigenetic signatures; the umbilical cord connects the foetus to the placenta and the umbilical artery transports blood from the foetus to the placenta. To determine if results replicated in an independent cohort, we analysed DMRs in a publicly available cord blood dataset from the Gene Expression Omnibus (GEO) repository (GEO Accession GSE129841). We also conducted an EWAS to identify DMPs in each tissue.

Methods

Study population:

The PRISM pregnancy cohort has been previously described [29]. Briefly, PRISM is a cohort of mother-child pairs designed to investigate associations between perinatal stress, other environmental exposures, and child development. Women ≥18 years old and in the first or second trimester of pregnancy were recruited from the Beth Israel Deaconess Medical Center in Boston, MA from March 2011 to August 2012 and from the Icahn School of Medicine at Mount Sinai in New York City from November 2012 to August 2014. Additional inclusions criteria were having a singleton birth and English or Spanish speaking. Exclusion criteria were HIV+, reported drinking ≥7 alcoholic drinks/week before pregnancy recognition, or any alcohol during pregnancy. Maternal sociodemographics was collected at the first visit.

Ethics:

The PRISM study protocol was approved by the Institutional Review Boards at Brigham and Women’s Hospital (BWH) and the Icahn School of Medicine at Mount Sinai. All participants provided written informed consent in their primary language.

Sample collection:

Sample collection and DNAm measurement of cord blood and placenta tissues has previously been described [30,31]. Immediately following birth, cord blood, umbilical cord, and placenta samples were collected. Cord blood was separated into plasma and buffy coat using centrifugation and stored at −20°C. Blood was drained from cord samples, and samples were stored in saline solution at 4°C. Within 24 hours, samples were dissected by separating the arteries and vein, and stored at −80°C until DNA isolation. Immediately after birth, the placenta was collected [32]. Four ~1 cm³ samples were extracted from the foetal side ~4 cm from the cord insertion site and 1–1.5 cm below the foetal membrane, and the decidua and foetal membranes were removed. Samples were rinsed in a cold phosphate-buffered saline bath, cut into ~0.1 cm³ pieces, and placed in RNAlater RNA Stabilization Reagent (Quiagen, Germantown, MD) for ≤ 24 hours; excess RNAlater was then removed and samples were stored at −80°C until DNA isolation.

DNA was isolated using the Gentra Puregene kit (Quiagen, Germantown, MD) according to the manufacturer’s instructions and quantified using an Implen Nanophotometer Pearl (Westlake Village, CA).

Epigenome-wide DNAm analysis:

DNAm measurement:

Sufficient DNA was available from 181, 234, 231 cord blood, placenta, and artery samples, respectively (Supplemental Figure 1). DNA was bisulphite converted and DNAm was measured using the 450K array (Illumina, San Diego, CA) according to the manufacturer’s instructions by the Partners HealthCare Translational Genomics Core Labs. The 450K array measures single nucleotide-resolution DNAm at >485,000 methylation sites [33]. To minimize batch effects, stratified randomization was used to randomize samples to chips and plates; the position of samples within plates was also determined by randomization.

Quality control and DNAm data processing:

Raw methylation .idat files were processed using the R package ewastools [34,35]. The following quality control steps were implemented: samples that were outliers in principal components analysis (PCA) and potentially contaminated samples were removed, genetic fingerprinting using 65 Illumina SNP probes was performed and compared across tissues from the same donor, sex was confirmed using X and Y chromosome probe intensities, and dye bias was corrected for using RELIC [36]. One cord blood sample, one placenta sample, and two artery samples were removed in quality control checks. Observations with detection p-values < 0.05 were removed [35]. Density distributions of Beta-values were examined and samples that did not follow the expected bimodal distribution were removed (cord blood n = 1; placenta n = 4) (Supplemental Figure 1). Technical replicates from 6, 5 and 4 for cord blood, placenta, and artery samples, respectively, were included across chips and plates. DNAm levels of replicated samples were highly correlated, indicating good precision (cord blood r_Pearson range: 0.98–0.99; placenta r_Pearson range: 0.98–0.99; placenta r_Pearson range: 0.98–1.00).

Beta-values were transformed into M-values using the formula log₂(Beta/(1-Beta)). PCA was performed and the top PCs were plotted to identify batch effects of plate, row, well, and column. Missing values were imputed as the probe median M-value, data were adjusted for batch effects of plate using the ComBat [37] implemented in the sva package [38], and imputed values were replaced as missing.

Our analyses focused on autosomal chromosomes; probes annotated to X and Y chromosomes (11,648) were removed. In addition, probes previously identified as cross-reactive (29,233) by Chen et al. [39], probes located within 10 base pairs of a known SNP with minor allele frequency ≥ 0.05 in an American admixture population (34,167), and probes with missing values in ≥ 10% of samples were dropped (cord blood: 559; placenta: 426; artery: 401). The number of samples and probes available for analysis were 179 and 414,940 for cord blood, 229 and 415,073 placenta, and 229 and 415,092 for artery samples, respectively (Supplemental Figure S1).

Four ancestry-related PCs were estimated using EPISTRUCTURE [40]. Cell type heterogeneity was adjusted for using ReFactor, a reference-free approach that uses PCA implemented on a sparse set of methylation sites to estimate components that can serve as a proxy for cell-type proportion estimates [41]. Before implementing ReFactor, missing values were imputed as the probe median M-value and M-values were transformed to Beta-values. The number of ReFactor components (k) was selected based on the ratio between adjacent eigenvalues, which provides information on axes of variance in data, as implemented by the estimate_k function (plots of the eigenvalue ratios are provided in Supplemental Figures S6–S8). We used six components for cord blood, eight components for placenta, and six components for artery. We used 500 sites (i.e., number of informative features: t) to calculate refactor components. The selected value of t was confirmed using the estimate_t function. For sensitivity analyses, cell-type proportions of cord blood samples was estimated using the Houseman projection method [42] implemented through ewastools using the Gervin et al. [43] and Bakulski et al. [44] cord blood reference datasets. The Gervin et al. reference set provides estimates for the proportions of CD8 + T cells, CD4 + T cells, natural killer (NK) cells, B cells, monocytes, and granulocytes, whereas the Bakulski et al. reference set also includes nucleated red blood cell (nRBC) estimates. As described in Gervin et al. [43] nRBC DNAm profiles are not fully characterized and may have substantial inter-individual variability, which could contribute to decreased validity of estimated nRBC proportions. For placenta samples, the proportions of trophoblasts, stromal cells, Hofbauer cells, endothelial cells, nRBC, and syncytiotrophoblasts were estimated using the minfi projectCellType function and the third trimester (term) Yuan et al. placenta reference dataset [45]. The reference dataset includes 600 CpGs with methylation level measured using the Illumina EPIC methylation array and therefore 574 probes common to our placenta data and the reference dataset were used as input. Reference datasets are not available for artery.

Data analysis:

Descriptive statistics:

We calculated descriptive statistics (medians and IQRs for continuous variables, frequencies for categorical variables) for all covariates. Differences between tissues in the characteristics of mother-infant pairs were assessed using the Mann–Whitney test for continuous variables and χ² test for categorical variables. Pearson correlations for DNAm levels across tissues were calculated for each individual with DNAm available for multiple tissues.

Differentially methylated regions:

Differentially methylated positions were identified using the R package DMRcate [46]. This method fits a linear model to all probes and applies a Gaussian kernel smoother to test statistics within λ base pairs. Smoothed test statistics and p-values are calculated using the Satterthwaite method and significant probes are determined after an FDR adjustment [47]. Significant probes within λ base pairs are grouped to determine DMRs. DMRcate is agnostic to the direction of association at individual probes between DNAm and the predictor. We implemented DMRcate on M-values with λ = 1,000, the FDR cut-off = 0.01, and the minimum number of probes in a DMR = 10. We selected a more stringent FDR cut-off and greater number of probes to define a DMR than the default parameters to allow us to identify DMRs more likely to have biological relevance. In our primary analysis, models were adjusted for ReFactor components.

Differentially methylated positions:

We identified DMPs associated with sex using the R package limma which implements linear models with empirical Bayes smoothing of standard errors. Limma was applied to M-values adjusting for ReFactor components. We applied a Bonferroni correction to p-values using the p.adjust function in R with calculates adjusted p-values as p × the number of tests (i.e., p_i,Bonferroni = p_i × N).

Gene ontology analysis:

GO analysis was conducted using the GOregion [48] function for DMRs or the GOmeth [49] function for DMPs in the missMethyl package [50]. GOregion and GOmeth identify GO terms with an overrepresentation of genes containing DMPs accounting for differences in the a priori probabilities of genes to include differential methylation based on representation on the 450K array. For each tissue, all probes tested and DMRs or DMPs (p_Bonferroni < 0.05) were used as arguments.

Sensitivity analyses:

Although we did not observe associations between infant sex and maternal age, race/ethnicity, or gestational age (Table 1), sensitivity analyses were performed adding these variables as precision covariates to the DMRcate and limma models. We also performed sensitivity analyses adjusting for cell-type proportions estimated using the Gervin et al. [43] and Bakulski et al. [44]. reference datasets for cord blood and the Yuan et al. reference dataset for placenta [45]. To determine if differences among infants with data of each tissue type may affect results, we performed analyses restricted to infants with data available for all three tissues.

Table 1. PRISM participant characteristics by tissue and associations with infant sex.

	Cord blood (N = 179)			Placenta (N = 229)			Artery (N = 229)
	N	%	p, association with sex ^a	N	%	p, association with sex ^a	N	%	p, association with sex ^a	p, difference between tissues ^b,c
Maternal characteristics
Age at delivery (yrs), median (IQR)	32.0	(28.6, 34.8)	0.37	31.3	(26.6, 34.6)	0.66	31.5	(26.7, 34.7)	0.77	0.42
Education ^d
≤ High school diploma/GED	34	19%	0.33	61	27%	0.14	61	27%	0.20	0.12
> High school diploma/GED	143	81%		164	73%		164	73%
Race/ethnicity ^e
White	68	38%	0.67	75	33%	0.38	78	34%	0.85	0.67
Black	69	39%		91	40%		91	40%
Hispanic	23	13%		43	19%		39	17%
Other	19	11%		18	8%		19	9%
Infant characteristics
Sex
Male	96	54%	-	126	55%	-	124	54%	-	0.96
Female	83	46%		103	45%		105	46%
Gestational age (wks), median (IQR)	39.1	(38.2, 40.0)	0.91	39.1	(38.1, 40.0)	0.35	39.1	(38.1, 40.0)	0.64	0.81
Preterm ^f	13	7%	0.76	21	9%	0.37	19	8%	0.56	0.79
Birth weight (g), median, (IQR)	3,380	(3,073, 3,620)	0.87	3,315	(3,010, 3,580)	0.66	3,330	(3,050, 3,610)	0.39	0.67
Low birth weight ^g	11	6%	0.32	15	7%	0.50	14	6%	0.61	0.98

a. P-value for associations with infant sex; Mann–Whitney test for continuous variables; χ² test for categorical variables. b. Kruskal–Wallis Test for continuous variables; χ² test for categorical variables. c. Participants with cord blood, placenta, and artery samples; N = 159; participants with cord blood and placenta samples: N = 165; participants with cord blood and artery samples: N = 169; participants with placenta and artery samples: N = 218. d. Cord blood: N = 177; placenta: N = 225; artery: N = 225. e. Placenta: N = 227; artery: N = 227. f. < 37 weeks gestation. g. < 2,500 g.

Replication analysis:

We tested for DMRs associated with sex in an independent cord blood dataset (GEO Accession GSE129841) including DNAm from 60 males and 54 female infants. Data were imported to R using GEOquery and minfi, and Beta-values were extracted from the GenomicRatioSet. Data were processed as previously described by Kadakia et al. [51]. Briefly, samples with > 5% of probes with detection p-values > 0.01 were dropped. Probes with detection p-values > 0.01 in > 5% of samples, probes measuring SNPs or non-CpG sites, probes located near SNPs, and probes located on sex chromosomes were removed. Data were corrected for background noise and channel colour bias and quantile normalized. Batch effects were adjusted for using ComBat. We additionally removed cross-reactive probes identified by Chen et al. [39]. and probes with minor allele frequencies > 0.05 located < 10 base pairs from the target site, resulting in 390,729 probes for analysis. ReFactor was applied to estimate six cell type-related components. DMRs were identified using DMRcate adjusting for ReFactor components with the FDR cut-off = 0.01, λ = 1,000, and the minimum number of CpGs = 10. DMPs were identified using limma adjusting for ReFactor components. We evaluated overlap between DMRs and DMPs (p_Bonferroni < 0.05) identified in the GEO and the PRISM cord blood analyses. In addition, we used a lookup approach to assess replication of PRISM cord blood DMPs with an FDR adjustment for multiple comparisons among probes included in both analyses.

All analysis were conducted in R 4.0.0 [52].

Results

Participant characteristics:

For these analyses, DNAm was available from N = 179 cord blood, N = 229 placenta, and N = 229 umbilical artery samples (Supplemental Figure S1). Overlap between tissues was 165 mother-infant pairs for cord blood and placenta, 169 for cord blood and artery, and 218 for placenta and artery; 159 mother-infant pairs had data available for all tissues. The median age of mothers at delivery was 31.3–32.0 years for each tissue (Table 1). Approximately 40% of mothers were Black, one-third were White, and one-quarter were Hispanic or other. Slightly over half of infants were male (cord blood: 54%; placenta: 55%; artery: 54%), and less than 10% were preterm or had low birth weight (cord blood: 6%; placenta: 7%; artery: 6%). Maternal age at delivery, education, race, gestational age, and birth weight were not significantly different between sexes or tissues (p < 0.05).

We assessed associations between infant sex and cell-type related components (i.e., ReFactor components) (Supplemental Excel Table S1). Infant sex was associated with the second component in cord blood (female vs. male: B = −3.15, p = 0.046), the seventh component in placenta samples (B = 0.35, p = 0.037), and the fifth component in artery samples (B = 0.42, p = 0.021). Using cellular proportions estimated using the Houseman regression calibration method, infant sex was associated with the proportion of CD8 + T cells (Bakulski et al. reference set 44 B = −0.02, p = 0.002) in cord blood, and with the proportion of trophoblasts (Yuan et al. reference dataset 45 B = 0.02, p = 0.037), Hofbauer cells (B = −0.003, p = 0.007), and syncytiotrophoblasts (B = −0.02, p = 0.034).

DNAm levels within individuals were highly correlated across tissues. For cord blood and placenta samples, the range of r_Pearson was 0.82–0.88 (N = 165; p < 0.001); for cord blood and artery samples, the range of r_Pearson was 0.88–0.93 (N = 169; p < 0.001); and for placenta and artery samples, the range of r_Pearson was 0.85–0.90 (N = 218; p < 0.001).

Differentially methylated regions:

In models adjusted for cellular heterogeneity, we identified 183, 257, and 419 DMRs in cord blood, placenta, and artery samples, respectively (summarized in Table 2) The top-10 significant DMRs for each tissue are listed in Table 3 and all DMRs are included in Supplemental Excel Tables S2-S4. DMRs were located on all chromosomes except chromosome 9. Within all three tissues, the highest proportion of DMRs was located in chromosome 6 (cord blood: 32, 17%; placenta: 52, 20%; artery: 57, 14%). DMRs also included multiple imprinted genes (cord blood: PLAGL1, INS); placenta: TP73, HOXA4, DLX5, GNAS, LRRTM1, OSBPL5, CDKN1C; artery: PLAGL1, DLK1, DLX5, KCNQ1, CMTM1, ATP10A, NNAT) [53]. Although the direction of association may not be consistent across a DMR, DMRs identified in cord blood and artery samples showed a trend towards higher methylation levels in females compared to males (cord blood mean difference in methylation for females vs. males > 0 in 163 DMRs (89%); artery mean difference in methylation > 0 in 391 DMRs (93%)) (Figure 1). However, DMRs in placenta samples tended to have higher methylation among males compared to females (mean difference in methylation for females vs. males > 0 in 100 DMRs (39%)).

Table 2. Summary of differentially methylation regions (DMRs) and positions (DMPs) associated with infant sex. DMRs identified using DMRcate with the FDR cut-off = 0.01, λ = 1,000, and the minimum number of CpGs = 10. DMPs identified using limma applied to M-values. All models adjusted for cellular heterogeneity using ReFactor components.

	# of DMRs	# (%) of DMRs with higher mean methylation among females	# of DMPs	# (%) of DMPs with higher methylation among females
Cord blood ^a	183	163 (89%)	2,119	1,948 (92%)
Placenta ^b	257	100 (39%)	2,281	291 (13%)
Artery ^c	419	391 (93%)	3,405	3,179 (93%)

DMR: differentially methylated region. a. N = 179. b. N = 229. c. N = 229.

Table 3. Ranking of top-10 DMRs associated with infant sex. DMRs identified using DMRcate modelling females vs. males adjusting for ReFactor components with the FDR cut-off = 0.01, λ = 1,000, and the minimum number of CpGs = 10. DMRs ordered by Fisher combined probability transform of the individual CpG FDRs.

Chr	Start position ^a	End position	# of CpGs	Fisher	Mean difference in % methylation ^b	Overlapping genes
Cord blood (N = 179)
6	32,116,841	32,122,436	101	5.36E-77	1.47	PPT2, PPT2-EGFL8, PRRT1
20	30,071,726	30,074,596	12	9.20E-68	5.75	REM1, LINC00028
15	75,018,150	75,019,376	27	6.09E-62	2.39
11	2,292,000	2,293,665	37	3.97E-61	2.67	ASCL2
1	2,344,133	2,345,894	22	1.29E-60	1.79	PEX10
3	195,942,717	195,943,713	10	1.83E-60	3.38	SLC51A, PCYT1A
6	28,583,655	28,584,464	18	1.85E-57	3.05	SCAND3
17	48,911,410	48,913,424	21	3.15E-57	2.63	WFIKKN2
1	78,511,140	78,512,129	12	2.74E-50	1.78	GIPC2
6	28,601,269	28,601,519	13	1.08E-45	6.99	RP11-373N24.2
Placenta (N = 229)
1	3,619,574	3,624,292	18	5.58E-77	−3.79	TP73
20	61,146,941	61,149,399	22	4.71E-70	−3.01	C20orf166, C20orf166-AS1
12	130,501,516	130,502,586	12	9.48E-48	−4.55
1	53,904,065	53,905,934	15	8.66E-46	2.46	SLC25A3P1
21	27,106,440	27,108,257	13	1.61E-45	−0.67	GABPA, ATP5J
11	45,713,046	45,719,826	31	2.15E-43	−2.30
8	146,227,573	146,229,440	20	3.99E-43	−4.53	ZNF252P-AS1, ZNF252P
12	53,241,452	53,243,489	10	4.33E-40	−2.49	KRT78
1	174,967,663	174,968,980	18	1.91E-39	1.26	CACYBP
6	33,864,689	33,869,169	23	2.20E-38	−2.81	RP3-468B3.2, LINC01016
Artery (N = 229)
17	46,668,715	46,680,461	52	7.78E-231	3.86	HOXB-AS3, HOXB3, HOXB5, HOXB6
2	119,602,212	119,617,128	70	1.28E-191	2.75	EN1
1	2,344,089	2,345,894	24	3.15E-166	2.75	PEX10
15	83,951,663	83,954,849	31	2.44E-143	2.20	RP11-382A20.4, BNC1
2	177,012,117	177,017,797	35	2.09E-127	4.95	HOXD3, MIR10B, HOXD4
13	49,792,685	49,796,489	14	1.58E-111	3.02	MLNR
17	46,799,479	46,805,563	30	1.35E-87	1.90	PRAC2, PRAC1, MIR3185, HOXB13
17	46,655,796	46,660,940	23	3.83E-87	4.37	HOXB-AS3, HOXB3, HOXB4, MIR10A
6	144,385,261	144,386,848	24	8.61E-80	1.77	PLAGL1
15	75,018,150	75,019,376	27	2.93E-77	1.88

Chr = chromosome. a. GRCh37/hg19 assembly. b. Mean difference in % methylation for females vs. males.

Figure 1. Overlap between DMRs in cord blood, placenta, and artery samples associated with infant sex. DMRs identified using DMRcate modelling females vs. males and adjusting for ReFactor components with the FDR cut-off = 0.01, λ = 1,000, and the minimum number of CpGs = 10. DMRs with consistently higher mean methylation among females are represented in Orange, DMRs with consistently higher mean methylation among males are represented in purple, and DMRs with different mean directions of association between tissues are represented in grey. Mean differences in % methylation for each DMR are displayed in the scatter plot. Cord blood: N = 179; placenta: N = 229; artery: N = 229.

The number of DMRs that overlapped between tissues (i.e., share at least one CpG) are shown in Figure 1. We identified 39 DMRs overlapping across all cord blood, placenta, and artery samples (Table 4). The direction of association between mean methylation levels and infant sex was consistent across tissues for 32 overlapping DMRs, with 26 having consistently higher methylation among females and 6 having higher methylation among males (Figure 1). Eight (21%) the DMRs overlapping across tissues were located on chromosome 6, including 6 DMRs within the human leukocyte antigen (HLA) region (chr6:28,477,797–33,448,354). To illustrate this trend in association, Beta-values of the DMR chr6:28,601,269–28,601,519 are plotted in Figure 2.

Table 4. Overlapping differentially methylated regions (DMRs) associated with infant sex in cord blood, placenta, and artery samples. DMRs identified using DMRcate modelling females vs. males adjusting for ReFactor components with the FDR cut-off = 0.01, λ = 1,000, and the minimum number of CpGs = 10. DMRs ordered by chromosomal coordinate.

			Cord blood	Placenta	Artery
Chr	Start position ^a	End position	Mean difference in % methylation ^b	Mean difference in % methylation ^b	Mean difference in % methylation ^b	Overlapping genes
1	26,232,706	26,234,015	1.34	0.54	1.26	STMN1, MIR3917
1	111,506,021	111,507,160	−0.04	0.04	−0.19	RP11-96K19.5, LRIF1
1	154,244,759	154,245,414	0.61	0.49	0.57	HAX1
1	174,967,663	174,968,980	1.29	1.26	1.12	CACYBP
2	21,265,912	21,267,334	3.32	−2.76	2.28	APOB
2	27,850,964	27,852,231	0.64	1.11	0.46	ZNF512, RP11-158I13.2, GPN1, CCDC121
2	73,461,068	73,461,773	−0.08	−0.13	−0.10	CCT7
3	46,758,813	46,759,698	−3.80	−1.88	−3.05	PRSS50
6	27,258,466	27,259,237	1.31	1.30	1.17	POM121L2
6	28,226,885	28,227,482	2.09	2.41	2.17	NKAPL, ZKSCAN4
6	28,601,269	28,601,519	6.99	1.80	6.58	RP11-373N24.2
6	28,972,348	28,973,828	0.96	2.77	1.21	ZNF311
6	30,094,947	30,095,802	1.67	−1.91	0.78
6	31,631,638	31,633,028	0.87	0.31	0.75	CSNK2B, GPANK1
6	32,116,841	32,122,436	1.47	0.78	0.95	PPT2, PPT2-EGFL8, PRRT1
6	33,244,976	33,246,488	0.68	−0.83	0.74	B3GALT4
7	24,796,981	24,798,175	1.94	0.99	1.09	DFNA5
7	27,182,493	27,185,136	−2.20	0.41	1.04	HOXA-AS3, HOXA3, RP1-170O19.22, HOXA5, HOXA6
7	50,517,714	50,518,921	1.78	1.53	1.36	FIGNL1
7	95,025,194	95,027,158	1.28	2.17	1.09	PON1, PON3
7	150,019,955	150,021,553	2.52	1.41	1.62	LRRC61, ACTR3C
8	22,132,563	22,133,076	2.90	1.27	3.85	PIWIL2, CTD-2530N21.4
10	131,263,962	131,265,796	0.77	0.03	0.96	MGMT
11	2,291,652	2,293,665	2.67	1.41	1.00	ASCL2
12	1,099,823	1,101,044	0.53	0.32	0.50	ERC1, RP11-359B12.2
13	36,871,646	36,872,346	1.28	1.36	1.18	SOHLH2, CCDC169-SOHLH2, CCDC169
13	114,291,977	114,292,869	−0.79	−0.70	−0.71	TFDP1
15	75,018,150	75,019,376	2.39	2.36	1.88
16	89,557,017	89,558,207	−0.12	−0.66	−0.01	SPG7
17	7,486,551	7,487,249	1.08	0.92	1.47	MPDU1, AC113189.5
17	33,913,953	33,914,828	0.72	0.51	0.66	AP2B1
17	37,123,638	37,124,558	4.56	0.59	4.31	FBXO47
17	48,911,410	48,913,424	2.63	−1.47	1.43	WFIKKN2
19	19,625,264	19,627,257	1.56	1.05	2.06	NDUFA13, CTC-260F20.3, YJEFN3, TSSK6
19	55,972,504	55,973,778	−1.36	−0.98	−2.16	ISOC2
19	57,630,034	57,631,478	1.66	1.52	1.78	USP29
19	57,741,325	57,742,444	3.28	3.56	4.74	AURKC
20	25,128,211	25,129,562	1.71	−1.94	1.78	RP4-738P15.1
21	27,106,440	27,108,257	−0.78	−0.67	−0.85	GABPA, ATP5J

Chr = chromosome. a. GRCh37/hg19 assembly. b. Mean difference in % methylation for females vs. males.

Figure 2. Scatter plot of DMR chr6: 28,601,269–28,601,519.

Results of GO analyses are summarized in Table 5; a complete list of significant GO terms is included in Supplemental Excel Table S5. In cord blood, we identified one GO term with on overrepresentation of genes including DMRs (homophilic cell adhesion via plasma membrane adhesion molecules), and in placenta samples, we identified nine closely related GO terms involving organ and skeletal system morphogenesis (p_FDR < 0.05). A total of 39 GO terms were identified for artery samples, including embryo development; organ, skeletal system, and eye morphogenesis; male meiotic nuclear division; DNA binding and regulation of gene expression by RNA polymerase II; and chromosome and chromatin (p_FDR < 0.05).

Table 5. Summary of Gene Ontology (GO) analyses of differentially methylated regions (DMRs) and positions (DMPs). One, 3, and 39 GO terms were identified with an overrepresentation of DMRs associated with infant sex in cord blood, placenta, and artery samples, respectively (p_FDR < 0.05). Twenty-eight and 9 GO terms were identified with an overrepresentation of DMPs associated with infant sex in cord blood in placenta and artery samples, respectively (p_FDR < 0.05). No GO terms were associated to infant sex in cord blood (p_FDR > 0.05). Closely related GO terms have been combined; for a complete list of sex-associated GO terms, see Supplemental Tables S2 and S4.

	Cord blood	Placenta		Artery
	DMRs	DMRs	DMPs	DMRs	DMPs
Biological process	Homophilic cell adhesion via plasma membrane adhesion molecules	Organ morphogenesis, including skeletal system morphogenesis	Tissue development including epidermis development and keratinization; System development; Peptide cross-linking; Monocyte, lymphocyte, and eosinophil chemotaxis; Response to chemokine, including chemokine-mediated signalling pathway	Tissue development; Embryo development; Pattern-specific processes including regionalization; Organ morphogenesis, including skeletal system and eye morphogenesis; Cranial nerve development; Male meiotic nuclear division
Molecular function			Structural constituent of skin epidermis; Serine-type endopeptidase inhibitor activity; Chemokine receptor binding; Voltage-gated cation channel activity	Regulation of gene expression including regulation of transcription by RNA polymerase II; DNA binding transcription factor activity including RNA polymerase II-specific TAP1 binding	Transcription regulatory region and sequence-specific DNA binding; DNA binding transcription factor activity including RNA polymerase II-specific
Cellular component			Keratin filament and cornified envelope	Chromosome and chromatin

Table 6. Top-10 DMPs associated with infant sex ranked by p-value. DMPs identified using limma applied to M-values modelling females vs. males and adjusting for ReFactor components. Effect size estimates from limma applied to Beta-values. Bolded CpGs were identified in more than one tissue.

CpG	Chr	Position ^a	Gene	Feature category	P	PBonferroni	Female	Male	Mean difference in % methylation (95% CI)^b
Cord blood (N = 179)
cg12691488	1	243,053,673			2.11 ×10^−120	8.75 ×10^−115	7.83	32.29	−24.09	(−24.98, −23.20)
cg11643285	3	16,411,667	RFTN1	Body	1.51 ×10^−114	6.27 ×10^−109	91.84	75.59	15.92	(15.37, 16.48)
cg03691818	12	53,085,038	KRT77	Body	1.09 ×10^−74	4.51 ×10^−69	12.97	2.58	9.50	(8.73, 10.26)
cg17238319	3	16,428,391	RFTN1	Body	1.31 ×10^−74	5.42 ×10^−69	85.14	70.45	13.65	(12.73, 14.56)
cg03618918	1	160,865,097			8.10 ×10^−72	3.36 ×10^−66	69.55	79.99	−10.19	(−10.82, −9.55)
cg25294185	11	65,487,814	RNASEH2C	Body	6.92 ×10^−71	2.87 ×10^−65	4.46	14.57	−10.06	(−10.92, −9.20)
cg26921482	16	2,570,283	AMDHD2	TSS200	1.45 ×10^−65	6.01 ×10^−60	46.69	17.98	27.43	(25.50, 29.36)
cg26355737	13	114,292,172	TFDP1	Body	2.90 ×10^−60	1.20 ×10^−54	81.42	89.04	−7.82	(−8.44, −7.20)
cg04946709	16	59,789,030	LOC644649	Body	3.45 ×10^−58	1.43 ×10^−52	70.86	82.86	−11.99	(−13.03, −10.94)
cg06710937	13	23,489,940			5.27 ×10^−58	2.19 ×10^−52	9.40	2.74	6.61	(5.88, 7.35)
Placenta (N = 229)
cg02325951	14	89,878,619	FOXN3	Body	1.41 ×10^−75	5.85 ×10^−70	60.36	70.44	−10.30	(−11.00, −9.61)
cg12691488	1	243,053,673			1.57 ×10^−67	6.50 ×10^−62	7.40	19.41	−12.57	(−13.79, −11.35)
cg17612569	21	27,107,221	GABPA;ATP5J;ATP5J;ATP5J	TSS200;5'UTR;1stExon;Body	6.47 ×10^−65	2.69 ×10^−59	2.97	9.82	−8.33	(−9.12, −7.54)
cg13192180	17	65,717,569	NOL11	Body	1.82 ×10^−63	7.57 ×10^−58	72.86	86.14	−12.88	(−13.96, −11.80)
cg07355069	5	43,289,541	HMGCS1	3ʹUTR	3.49 ×10^−62	1.45 ×10^−56	85.23	93.19	−7.49	(−8.12, −6.86)
cg03691818	12	53,085,038	KRT77	Body	3.81 ×10^−53	1.58 ×10^−47	12.11	4.16	6.67	(5.93, 7.40)
cg02606018	12	10,658,281			8.07 ×10^−53	3.35 ×10^−47	23.38	29.68	−6.65	(−7.29, −6.01)
cg12780322	4	57,334,163	SRP72	Body	1.08 ×10^−48	4.48 ×10^−43	8.80	5.34	3.27	(2.89, 3.65)
cg15817705	1	209,406,063			1.06 ×10^−43	4.41 ×10^−38	60.52	68.49	−7.50	(−8.33, −6.68)
cg03986557	1	202,430,474	PPP1R12B;PPP1R12B	Body;TSS1500	1.24 ×10^−43	5.16 ×10^−38	90.96	84.28	7.01	(6.21, 7.81)
Artery (N = 229)
cg12691488	1	243,053,673			7.34 ×10^−143	3.05 ×10^−137	5.74	30.84	−25.02	(−26.03, −24.00)
cg11643285	3	16,411,667	RFTN1	Body	1.32 ×10^−120	5.48 ×10^−115	85.31	70.11	15.05	(14.46, 15.65)
cg03691818	12	53,085,038	KRT77	Body	2.59 ×10^−101	1.08 ×10^−95	13.31	2.49	10.01	(9.26, 10.76)
cg03618918	1	160,865,097			1.39 ×10^−95	5.76 ×10^−90	68.63	78.35	−9.81	(−10.32, −9.30)
cg26921482	16	2,570,283	AMDHD2	TSS200	7.88 ×10^−92	3.27 ×10^−86	43.94	13.27	29.01	(27.31, 30.72)
cg02325951	14	89,878,619	FOXN3	Body	1.78 ×10^−87	7.37 ×10^−82	61.70	72.73	−11.08	(−11.74, −10.41)
cg17612569	21	27,107,221	GABPA;ATP5J;ATP5J;ATP5J	TSS200;5'UTR;1stExon;Body	3.97 ×10^−79	1.65 ×10^−73	3.10	12.16	−10.35	(−11.23, −9.48)
cg25294185	11	65,487,814	RNASEH2C	Body	6.43 ×10^−79	2.67 ×10^−73	5.13	14.88	−10.15	(−10.92, −9.39)
cg04946709	16	59,789,030	LOC644649	Body	4.55 ×10^−75	1.89 ×10^−69	69.59	80.20	−10.19	(−10.95, −9.42)
cg25568337	6	157,098,338	ARID1B	TSS1500	5.17 ×10^−73	2.15 ×10^−67	19.62	11.01	8.11	(7.51, 8.70)

a. GRCh37/hg19 assembly. b. Mean difference in % methylation for females vs. males.

Sensitivity analyses:

Results of sensitivity analyses are summarized in Supplemental Excel Table S6, including the total number of DMRs identified and DMRs overlapping with our primary analyses. Overall, the number of DMRs identified and directions of association were consistent with our primary analyses. Including the covariates of ancestry-related PCs, maternal age, and gestational age, we identified 179, 277, and 421 DMRs in cord blood, placenta, and artery samples, respectively, of which 96%, 88%, and 95% overlapped with DMRs in our primary analyses. In analyses limited to infants with data available for all three tissues (N = 159), we identified fewer DMRs than analyses including all samples but with high overlap (cord blood: 169 DMRs with 90% overlapping with primary analyses; placenta: 183 DMRs with 89% overlapping; artery: 226 DMRs with 98% overlapping).

We also conducted sensitivity analyses for cellular heterogeneity corrections in cord blood and placenta by estimating the cell proportions using the Houseman regression calibration method[42]. In cord blood, adjusting for CD8 + T cells, CD4 + T cells, NK cells, B cells, monocytes, and granulocytes, we identified 159 DMRs (90% overlapping with primary analyses); including nRBCs, we identified 202 DMRs (77% overlapping). In placenta samples, adjusting for trophoblasts, stromal cells, Hofbauer cells, endothelial cells, nRBCs, and syncytiotrophoblasts, we identified 162 DMRs (81% overlapping).

Replication analysis:

We conducted analyses to identify DMRs in a cord blood DNAm dataset from an independent cohort (GEO Accession GSE129841)[51]. In analyses adjusted for ReFactor components, we identified 47 DMRs associated with infant sex (Supplemental Excel Table S7). As observed in our primary analyses, 44 DMRs (94%) had higher mean methylation levels among females, and the highest proportion (19%) were located in chromosome 6. There was substantial overlap between DMRs identified in the GEO analyses and analyses of associations in PRISM cord blood samples; 40 (85%) of GEO DMRs overlapped with DMRs identified in PRISM, 38 of which had consistent median mean methylation differences between sexes and 16 of which had complete overlap between cohorts.

Differentially methylated positions:

In our primary analyses adjusted for cellular heterogeneity using ReFactor components, we identified 2,119 DMPs associated with infant sex in cord blood, 2,281 DMPs in placenta, and 3,405 DMPs in artery samples (p_Bonferroni < 0.05, Table 2). Quality checks and results visualization for each tissue are in Supplemental Figures S2-S4 and Figure 3; DMPs (p_Bonferroni < 0.05) are in Supplemental Excel Tables S8-S10. In cord blood and artery samples, DMPs predominantly had higher methylation levels among females compared to males (92% and 93%, respectively); however, in placenta samples, only 13% of DMPs had higher methylation among females compared to males. The highest proportion of DMPs were located near transcription start sites (TSS200 or TSS1500, cord blood: 34%; placenta: 23%; artery: 34%) or in gene bodies (cord blood: 29%; placenta 34%; artery: 30%), followed by untranslated regions (5ʹ or 3ʹ UTR, cord blood: 14%; placenta: 16%; artery: 14%) and exons (cord blood: 7%; placenta: 6%; artery: 6%). Across tissues, approximate 30% of DMPs were located in intergenic regions (cord blood: 29%; placenta: 33%; artery: 28%).

Figure 3. Volcano and Manhattan plots for DMPs associated with infant sex.

The top-10 most significant DMPs in each tissue are listed in Table 6; estimates from models of Beta-values are included to ease results interpretation. Notably, in cord blood and artery samples, a smaller proportion of the top-10 most significant DMPs had higher methylation among females (5 DMPs in cord blood, 3 DMPs in placenta, and 4 in artery samples). cg12691488, located in an intergenic region of chromosome 1, was the most significant DMP in cord blood and artery samples and the second most significant DMP in placenta samples (cord blood: p_Bonferroni = 8.75 × 10⁻¹¹⁵; placenta: p_Bonferroni = 6.50 × 10⁻⁶²; artery: p_Bonferroni = 3.05 × 10⁻¹³⁷). Additional highly significant DMPs identified across tissues were annotated to FOXN3 (cg02325951, placenta: p_Bonferroni = 5.85 × 10⁻⁷⁰; artery: p_Bonferroni = 7.37 × 10⁻⁸²), RFTN1 (cg11643285, cord blood: p_Bonferroni = 6.27 × 10⁻¹⁰⁹; artery: p_Bonferroni = 5.48 × 10⁻¹¹⁵), KRT77 (cg03691818, cord blood: p_Bonferroni = 4.51 × 10⁻⁶⁹; placenta: p_Bonferroni = 1.58 × 10⁻⁴⁷; artery: p_Bonferroni = 1.08 × 10⁻⁹⁵), RNASEH2C (cg25294185, cord blood: p_Bonferroni = 2.87 × 10⁻⁶⁵; artery: p_Bonferroni = 2.67 × 10⁻⁷³), AMDHD2 (cg26921482, cord blood: p_Bonferroni = 6.01 × 10⁻⁶⁰; artery: p_Bonferroni = 3.27 × 10⁻⁸⁶), and GABPA; ATP5J (cg17612569, placenta: p_Bonferroni = 2.69 × 10⁻⁵⁹; artery: p_Bonferroni = 1.65 × 10⁻⁷³). For each of these DMPs, the direction of association was consistent across tissues.

We observed overlap of DMPs across tissues, with the largest overlap (805 DMPs) in cord blood and artery samples (Supplemental Figure S5). A total of 127 DMPs were common across all three tissues. Substantially higher overlap was observed among the most significant DMPs. When comparing the top 100 DMPs ranked by p-value for each tissue, 29 DMPs were identified in all three tissues, and an addition 37 DMPs were identified in both cord blood and artery samples.

In GO analysis, no GO term was associated with infant sex in cord blood (p_FDR > 0.05). However, 28 and 9 GO terms had an overrepresentation of DMPs in placenta and artery samples, respectively (p_FDR < 0.05) (summarized in Table 5; complete results listed in Supplemental Excel Table S11). GO terms associated with sex in placenta samples were related to development, particularly of the epithelium, and immune function (e.g., chemokine-mediated signalling pathways), whereas significant GO terms in artery samples were related to transcription regulation, including DNA binding by RNA polymerase II.

Sensitivity analyses:

Sensitivity analyses of DMPs are summarized in Supplemental Table S6. In models including ancestry-related PCs, maternal age at delivery, and gestational age, we identified 2,086 DMPs in cord blood (p_Bonferroni < 0.05; 93% identified in our primary analyses); 2,805 DMPs in placenta samples (79% identified in our primary analyses); and 3,210 DMPs in artery samples (98% identified in our primary analyses). Restricting analyses to infants with data available for all three tissues, we identified fewer DMPs but observed substantial overlap with DMPs in our primary analyses: of the 2,012, 1,416, and 1,586 DMPs (p_Bonferroni < 0.05) in cord blood, placenta, and artery, respectively, 95%, 92%, and 97% were identified in our primary analyses. Similar to our primary analyses, among the top-100 DMPs, 29 were common to cord blood, placenta, and artery analyses; however, among all DMPs, only 69 were common across tissues.

Overall results were consistent when adjusting for Houseman-estimated cell-type proportions. For analyses of cord blood including the proportions of CD8 + T cells, CD4 + T cells, NK cells, B cells, monocytes, and granulocytes we identified 2,192 DMPs (p_Bonferroni < 0.05), of which 83% were also significant in ReFactor-adjusted analyses (p_Bonferroni < 0.05). Using estimates from the Bakulski et al. dataset [44], which includes nRBC proportions, we found 2,114 DMPs (p_Bonferroni < 0.05), including 84% overlapping with DMPs identified in ReFactor-adjusted analyses (p_Bonferroni < 0.05). In placenta analyses, including the proportions of trophoblasts, stromal cells, Hofbauer cells, endothelial cells, nRBCs, and syncytiotrophoblasts, we identified fewer DMPs than in ReFactor-adjusted analyses (1,325 vs. 2,281, p_Bonferroni < 0.05); however, 96% of these DMPs were also identified in our ReFactor-adjusted analyses.

Replication analysis:

We identified 514 DMPs in the GEO cord blood dataset using linear models adjusted for ReFactor components (p_Bonferroni < 0.05; Supplemental Excel Table S12), of which 432 (84%) were significant in cord blood in PRISM (p_Bonferroni < 0.05). As we observed in PRISM cord blood analyses, DMPs primarily had higher methylation among females (441; 86%). Additionally, we used a lookup approach of the 2,029 PRISM cord blood DMPs included in the replication analysis. Among these probes, 1,941 (96%) were significant at p_FDR < 0.05 (corrected for 2,029 comparisons) and 2,027 had a consistent direction of association between cohorts.

Consistency with previous studies:

We compared our results to DMPs reported previously reported in 450K-based EWAS of cord blood and placenta samples reported by Yousefi et al., Maschietto et al., and Martin et al. [22,24,25]. A total of 4,562 unique DMPs were identified in cord blood and placenta tissues in these three studies. Among 4,473 DMPs identified in cord blood by Yousefi et al. and Maschietto et al., 4,715 were included in our analyses of cord blood in PRISM, 2,679 (50%) of which were significant at p_Bonferroni < 0.05 and 4,870 (91%) were nominally significant at p < 0.05. Among the 91 DMPs identified in placenta samples by Martin et al., all were significant at p_Bonferroni < 0.05 in our analyses of placenta samples in PRISM.

Discussion

This study of differential methylation levels by infant sex is unique in assessing epigenetic patterns across three birth tissues. We identified 183, 257, and 419 DMRs in cord blood, placenta, and artery samples, respectively. We replicated analyses in cord blood DNAm data from an independent cohort (GEO Accession GSE129841); 85% of DMRs overlapped with those identified in PRISM.

There was substantial overlap between DMRs across tissues, with 39 overlapping DMRs identified in all three tissues, 84 identified in cord blood and artery samples, 30 identified in placenta and artery samples, and 11 identified in cord blood and placenta samples. Among DMRs overlapping across tissues, several were annotated to autosomal genes related to characteristics of sexual dimorphism. The TSS200 and TSS1500 regions of genes involved in spermatogenesis, including NFkB activating protein-like (NKAPL)[54], fidgetin-like1 (FIGNL1) [55], piwi like RNA-mediated gene silencing 2 (PIWIL2) [56], and aurora kinase C (AURKC) [57] were annotated to DMRs common across all three tissues with consistently higher mean methylation among females. Similarly, male meiotic nuclear division was a significant GO term for DMRs in artery samples (p_FDR < 0.05). A DMR overlapping with ligand dependent nuclear receptor interacting factor 1 (LRIF1) [58], involved in X–inactivation, was also identified in all tissues with higher methylation in cord blood and artery and lower methylation in placenta among females. In addition, sex-associated GO terms for placenta and artery regional analyses were related to organ and skeletal system development.

We also identified 2,119, 2,281, and 3,405 sex-associated DMPs in cord blood, placenta, and artery samples, respectively (p_Bonferroni < 0.05). In replication analyses in the GEO cord blood dataset, 84% of DMPs were significant in PRISM (p_Bonferroni < 0.05). Our results were also consistent with DMPs previously identified in EWAS of foetal tissues [22,24,25]; 50% and 100% of previously-reported DMPs in cord blood and placenta samples, respectively, were significant in our study (p_Bonferroni < 0.05). We observed greater overlap in DMPs across tissues and with previous studies among the most significant sites. Of the 127 DMPs (p_Bonferroni < 0.05) common to cord blood, placenta, and artery samples, 29 were among the top-100 significant sites in all three tissues. In addition, among the top-10 significant cord blood DMPs, seven were among the top-10 found in artery samples and six were among the top-10 identified in cord blood by Yousefi et al. [22].This proportion of overlap among highly significant DMPs with large effect sizes may indicate a strong epigenetic signature of sex at specific CpGs. However, two of these DMPs (cg17238319 and cg25294185) may also be a technical artefact caused by cross-hybridization of sex chromosomes to array probes target autosomal regions, as previously suggested [59].

Although we observed overlap in DMRs and DMPs across tissues, the trend in directions of association was not consistent. Overall, higher mean methylation was observed among females compared to males in cord blood and artery samples. DNAm of placenta samples was highly and positively correlated with that of cord blood and artery samples; however, higher mean methylation of DMRs and DMPs in placenta samples was observed among males compared to females. This direction of association in placenta has previously been reported [24] however, the biological reason is not understood. The placenta is unique as a maternal-foetal interface that functions in nutrient and waste transfer, thermo-regulation, and gas exchange. As the organ involved in foetal growth and resource allocation with balance between the maternal and paternal genomes, the placenta has high expression of paternal imprinted genes [60]. Considering the invasive nature of foetal extravillous trophoblasts, HLA genes expressed by the placenta are involved in moderating maternal immune tolerance [61]. Higher expression of autosomal immune-regulating genes in female placentas has also been observed, including interleukin 2 receptor, beta (IL2RB) and C-X-C motif chemokine ligand 1 (CXCL1) [62]. Although we did not observe differential methylation of imprinted genes previously reported as differentially expressed in the placenta [63], we did identify placental sex-related DMRs annotated the imprinted genes TP73, HOXA4, DLX5, GNAS, LRRTM1, OSBPL5, CDKN1C. In addition, we identified a DMP annotated to IL2RB (cg21307484, p_Bonferrni = 0.012; mean difference in methylation for females vs. males B = −3.76), and GO analysis of DMPs indicated an overrepresentation of GO terms related to chemokine activity among placenta DMPs, including chemokine-mediated signalling pathway (GO:0070098) and chemokine receptor binding (GO:0042379).

This study had several limitations. We were not able to assess sex differences in gene expression, which may provide more information about the biological relevance of differential methylation. In addition, use of bisulphite conversion to measure DNAm does not allow for differentiation between 5-methylcytosine and 5-hydroxymethylcytosine, an intermediate nucleotide modification in active demethylation [64]. 5-hydroxymethylcytosine is substantially less abundant in the human genome [65], but may provide information about the dynamics of DNAm patterns, particularly during periods of epigenetic reprogramming. Our study was also limited by the use of the 450K array to measure DNAm levels. Although this array includes over 99% of RefSeq genes [33], probes are not evenly distributed throughout the genome. Biases in representation could impact analyses of DMRs, particularly when restricted to regions containing a higher number of CpGs. It should also be noted that our choice of DMRcate for analysing regional associations may affect results. There are currently several available methods for detecting DMRs but no established gold-standard. Although DMRcate may have lower power than other approaches, evaluations have shown good control of Type-I errors [66] and greater performance measured using the area under the precision–recall curve than other methods [46]. Our replication analysis was limited in being restricted to cord blood due to publicly available DNAm datasets with sufficient sample sizes for other foetal tissue types. In addition, different preprocessing methods were applied to our cord blood dataset and the GEO replication dataset (e.g., quantile normalization was applied to GEO Accession GSE129841). Although this could limit replication, consistency in results between our primary and replication analyses demonstrates robustness of our findings.

The study was novel in including DNAm measured in multiple foetal tissues from the same cohort, allowing us to compare sex-associated DMRs and DMPs identified in cord blood, placenta, and artery samples. In addition, this was the first study to our knowledge to report associations between sex and DNAm measured in artery samples. Our study was strengthened by having a larger sample size than previous EWAS of infant sex in cord blood [22,25] or placenta samples [24], which increased statistical power. Furthermore, results of sensitivity analyses adjusting for maternal age, ancestry, and gestational age or Houseman cell-type proportions, as well as results of replication analyses in an independent birth cohort, were consistent with our primary analyses, indicating robustness of our findings.

Further research is necessary to understand biological mechanisms that may contribute to differences in DNAm signatures across tissues, as well as to determine if sexual dimorphism in the epigenome is related to differences in gene expression and response to environmental stressors.

Author contributions

Anne Bozack: conceptualization, methodology, formal analysis, investigation, writing – original draft, visualization; Elena Colicino: methodology, writing – review & editing; Allan Just: data curation; Robert Wright: writing – review & editing; Andrea Baccarelli: writing – review & editing; Rosalind Wright: conceptualization, writing – review & editing, supervision, project administration, funding acquisition; Alison Lee: conceptualization, writing – review & editing, supervision.

Acknowledgments

We would like to thank Dr. Andres Cardenas from the UC Berkeley School of Public Health for his expertise in analytical methods for conducting epigenome-wide association studies.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data Availability

R script for analysis of replication analyses using GEO Accession GSE129841 is available at https://github.com/annebozack/GEO_infant_sex_rep_analysis. Additional data available upon request to the corresponding author.

Supplementary material

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the National Institute of Environmental Health Sciences [R01 ES030302]; National Institute on Minority Health and Health Disparities [R01 MD013310]; National Heart, Lung, and Blood Institute [R01 HL114396]; National Heart, Lung, and Blood Institute [K23 HL135349]; National Heart, Lung, and Blood Institute [R01 HL095606]; The Grant [P30 ES023515] from the National Institute of Environmental Health Sciences.