Epigenome-850K-wide profiling reveals peripheral blood differential methylation in term low birth weight

Abstract

Aim: We investigate the genome-wide DNA methylation (DNAm) patterns of term low birth weight (TLBW) neonates. Methods: In the discovery phase, we assayed 32 samples (TLBW/control:16/16) using the EPIC 850k BeadChip Array. Targeted pyrosequencing of in 60 samples (TLBW/control:28/32) using targeted pyrosequencing during the replication phase. Results: The 850K array identified TLBW-associated 144 differentially methylated positions (DMPs) and 149 DMRs. Nearly 77% DMPs exhibited hypomethylation, located in the opensea and gene body regions. The most significantly enriched pathway in KEGG is sphingolipid metabolism (hsa00600), and the genes GALC and SGMS1 related to this pathway both show hypomethylation. Conclusion: Our analysis provides evidence of genome-wide DNAm alterations in TLBW. Further investigations are needed to elucidate the functional significance of these DNAm changes.

Plain language summary

This study looked at the DNA of babies born after 37 weeks of pregnancy but weighing less than 2500 grams. We found that these babies had lower levels of DNA methylation, which might change how their bodies handle fats.

TWEETABLE ABSTRACT

This study offers a comprehensive exploration of the DNA methylation landscape in TLBW. We identified TLBW-associated 144 DMPs and 149 DMRs. Several DMPs related genes involved in crucial processes such as immune response, transcriptional regulation and sphingolipid metabolism.

Article highlights

Epigenetic regulation

• Epigenetic regulation, particularly DNA methylation may play an important role in the development of term low birth weight (TLBW).

Genome-wide methylation patterns

• Our study presents a landmark investigation, marking the first-ever comprehensive comparison of genome-wide methylation patterns in the peripheral blood of TLBW and term normal birth weight (TNBW) neonates.

Differential methylation analysis

• For the first time, we identified 144 differentially methylated positions and 149 differentially methylated regions in the peripheral blood of TLBW neonates.

Implicated genes

• This novel revelation highlights the genes implicated in term low birth weight, notably KY, CEP164, ATG7, CHSY1, GADD45G, PRKAR1B and SH3BP5.

Distinctive methylation profile

• A striking observation was the predominance of hypomethylation, encompassing almost 77% of differentially methylated positions, suggesting that reduced DNA methylation levels constitute a hallmark feature of TLBW neonates.

Immune & inflammatory responses

• Our results resonate with the perturbed methylation of genes governing immune and inflammatory responses in TLBW neonates.

Sphingolipid metabolism

• Our findings accentuate the intriguing possibility of sphingolipid metabolism's involvement in the intricate mechanism underpinning TLBW.

Interdisciplinary investigations

• Future investigations are poised to combine diverse omics disciplines, including genomics, metabolomics and proteomics, to provide a more comprehensive perspective on TLBW's epigenetic mechanisms.

Keywords: :

• 850K
• DNA methylation
• epigenetics
• pyrosequencing
• term low birth weight
• TLBW

1. Background

Term low birth weight (TLBW), defined as neonates delivered at full term weighing less than 2500g, is still prevalent and remains an intractable global public health problem. In 2020, 19.8 million newborns, accounting for approximately 14.7% of all births globally that year, suffered from low birth weight (LBW) [1]. Notably, among neonates who are small for gestational age (SGA) in low- and middle-income countries, a striking 45% are classified as TLBW, while merely 6% are classified as preterm [2]. These newborns face a significantly elevated risk of morbidity and mortality within the first month of life and also remain vulnerable to the onset of other health conditions later in life. Accumulating evidence supports that TLBW is a risk factor for the development of adult kidney and cardiovascular disease [3–5]. Understanding the underlying mechanisms of TLBW and identifying associated risk factors are crucial steps toward reducing morbidity and improving their long-term health outcomes.

The major determinants for the significant variation in fetal birth weight are still not fully understood. TLBW is commonly associated with factors contributing to intrauterine growth restriction (IUGR). These factors can be classified into several general categories: socio-demographic factors; medical risks before pregnancy; risks of the current pregnancy (such as gestational hypertension and diabetes; maternal nutrition, placental causes); prenatal care; environmental and behavioral risks (involving exposure to pollutants, smoking, alcohol consumption or illicit drug use during pregnancy) [6].

In recent years, there has been a growing interest in investigating the role of genetics in the regulation of birth weight. Studies have identified several genetic variants associated with birth weight. Genome-wide association studies (GWAS) have played a crucial role in this area, revealing hundreds of genomic loci that are associated with fetal growth [7,8]. Maternal TCF7L2 Type 2 diabetes–risk genotypes at rs7903146 are associated with increased offspring birth weight [9]. Elevated placental expression of the imprinted PHLDA2 gene has been reported to be associated with LBW [10]. Maternal MTHFR C677T polymorphism was also associated with LBW [11].

Epigenetic regulation, encompassing processes such as DNA methylation (DNAm), histone post-translational modifications and noncoding RNAs, has emerged as a key mechanism in understanding the intricate interplay between genetics and the environment in disease pathogenesis [12–14]. DNAm, the most typical process of epigenetic regulation, occurs throughout the life cycle from embryo to adult and can directly inhibit transcription or indirectly suppress gene expression through transcriptional silencing [15,16]. Notably, epigenetic changes, particularly changes in DNAm profiles, have been repeatedly proposed as a potential molecular mechanism linking various prenatal risk factors to adverse maternal outcomes, such as LBW and preterm birth.

Several studies have tried to explore changes in the DNAm profile of TLBW neonates. Hayashi et al. analyzed umbilical cord blood from 4 TLBW cases and found 483 hyper-differentially methylated genes (DMGs) and 35 hypo-DMGs in the TLBW promoter region [17]. Similarly, MAO et al. detected differentially methylated regions (DMRs) in cord blood and placental DNA of TLBW neonates [18]. Tekola-Ayele et al. found that placenta methylation at 15 CpGs was associated with birthweight [19]. Moreover, specific gene methylation patterns have been reported to be associated with TLBW. For instance, Tian et al. found that per one standard deviation increase in the DNAm of CpG 11 at FGFR2, the risk of TLBW increased by 1.64-times [20]. Additionally, cord blood and maternal blood AHRR DNAm may be independently associated with TLBW in different ways [21].

Examination of cord blood or placental tissue may serve as indirect evidence of intrauterine conditions. Changes in DNA methylation patterns in peripheral blood provide more direct evidence. However, limited by the difficulty of collecting peripheral blood samples from newborns after birth, these past studies mainly studied the methylation levels of cord blood or placental tissue. Many studies mainly focus on target genomic regions, are relatively small in scale, and do not fully consider the influence of gestational age factors. There is currently a lack of large-scale methylation profiling analysis of postnatal peripheral blood of TLBW. In our study, we conducted a two-stage case-control study. For the first time, we performed genome-wide DNA methylation analysis using peripheral blood samples from TLBW neonates. We identified differentially methylated positions (DMPs) and performed bioinformatics analysis to explore the potential signaling pathways.

2. Materials & methods

2.1. Study population

We conducted a two-stage case-control study comprising a discovery stage and a replication phase. This study was approved by the China-Japan Friendship Hospital ethics committee board (number: 2023-KY-057). All methods were carried out in accordance with relevant guidelines and regulations. Parental informed consent was obtained in writing from each participant. The study design flowchart is depicted in Figure 1. The TLBW group included neonates (birth weight <2500 g, gestational age ≥37 weeks) born at a tertiary general hospital in Beijing, China between July 2019 and April 2023. We excluded mothers who smoked or consumed alcohol during pregnancy, had pregnancies resulting from in vitro fertilization, experienced multiple births or had neonates with birth defects. The healthy control comprised singleton full-term neonates with normal birth weight (term normal birth weight, TNBW group, 2500 g< birth weight <4000 g, gestational age ≥37 weeks) born at the same hospital during the same period. We meticulously matched the controls to the case group based on gestational age, sex and the mother's delivery age. The discovery phase included 32 neonates (16 in the TLBW group and 16 in the TNBW group), while the verification stage consisted of 60 neonates (28 in the TLBW group and 32 in the TNBW group).

Figure 1. Flowchart of study design.

We reviewed the hospital's electronic medical record system to obtain clinical information, including gender, gestational age, birth weight, parity, length and head circumference. We also gathered information regarding the mother's age at delivery, pre-pregnancy weight and height, weight gain during pregnancy, pre-pregnancy BMI, mode of delivery and the newborn father's age.

2.2. Sample collection & DNA extraction

Peripheral blood samples were collected from all participants prior to feeding and treatment within 4 h of birth. Each sample (2 ml) was collected in EDTA tubes and then stored at -80℃ prior to analysis. DNA was extracted from the peripheral blood samples using the DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer's protocol. DNA concentration and purity were examined by agarose gel electrophoresis and NanoDrop 2000.

2.3. 850K BeadChip DNA methylation profiling

We performed bisulfite conversion on the purified DNA samples using the DNA bisulfite conversion kit (Zymo Research) following the manufacturer's guidelines and protocol. Bisulfite-converted DNA was amplified, fragmented and then hybridized to the Illumina Infinium Human MethylationEPIC BeadChip(850K), which permits the simultaneous analysis of methylation at over 850,000 CpG sites throughout the genome. The BeadChip arrays were scanned on the Illumina iScan System, and raw methylation data were extracted using Illumina's Genome Studio Methylation Module. The software calculates the degree of DNAm from the raw signal value and expresses it as a β value ranging from 0 (unmethylated) to 1 (fully methylated). The β value represents the proportion of methylated cytosines at a specific CpG site. Subsequent processing and analysis of the methylation data were performed using the R programming language.

The Bioconductor R package minfi was utilized for data quality control. We conducted a series of filtering steps to exclude low-quality and irrelevant probe sites, as shown in Figure 2A. We filtered out probes with p-values greater than 0.01. We excluded non-CpG probes, sex-chromosome probes, probes with a single nucleotide polymorphism (SNPs) at the CpG site (minor allele frequency (MAF) ≥ 5%, and probes with the number of NBeads <3 in more than 5% of samples.

Figure 2. Quality control of the methylation array data. (A) Flowchart of the 850K chip probes filtration process. A total of 812215 finalized probes were generated. (B) Density plot of raw methylation β values distributions. Each line represents the mean of each individual in that group combination. (C) Principal component analysis (PCA) of all probes from the two groups of samples. Different groups are represented by different colors. All samples are in the range of 4 standard deviations of principal component 1 and principal component 2, respectively.

We applied several methods to correct for technical biases and batch effects. First, we utilized the BMIQ algorithm from the R package ChAMP. This algorithm performs Infinium type I and type II probe bias correction, ensuring that the methylation measurements from different probe types are normalized and comparable. Next, we conducted principal component analysis (PCA) of the raw beta values for the 850K BeadChips. It allowed us to understand the natural clustering of samples based on their methylation profiles. To correct for batch effects, we employed the ‘champ.runCombat’ function from the ChAMP package. This function utilizes a singular value decomposition (SVD) method to remove any unwanted technical variations that might arise due to different experimental batches or processing conditions. Furthermore, we implemented the SmartSVA algorithm for cell heterogeneity correction. We also used sex and gestational age as covariates for differential methylation analysis correction.

DMPs analysis was performed using the champ.DMP function of the ChAMP package. We set the criteria for statistical significance as follows: a p-value < 0.001 and |Δβ| >0.05 [22–24]. Specifically, a Δβ greater than 0.05 indicated hypermethylation, whereas a Δβ less than -0.05 indicated hypomethylation. In addition to single-site DMP analysis, we also utilized the DMRcate package to identify DMRs. A DMR was considered significant when it included at least 7 probes within a window of 200–10000 bp and the corrected p value after Benjamini & Hochberg correction was <0.05.

Gene Ontology (GO) analysis is a widely used method for functional enrichment studies, including three main categories: biological process (BP), molecular function (MF) and cellular component (CC). Similarly, the Kyoto Encyclopedia of Genes and Genomes (KEGG) is a valuable database that provides information on biological pathways, genetics, metabolism and environmental interactions. We conducted GO and KEGG pathway enrichment analyses of the DMPs-related genes using the ClusterProfiler package. The significance of enrichment was assessed using the hypergeometric test, and a p value less than 0.05 was considered as the threshold for determining the significant enrichment of differential genes in specific GO terms and KEGG pathways.

2.4. Validation of DNA methylation by pyrosequencing

Four significant CpG sites (cg09530217, cg16955800, cg24312537, cg27041096) covering four genes were validated with pyrosequencing. These CpGs were chosen for their significance, Δβ and proximity to a gene promoter. Genomic DNA was extracted from the peripheral blood of the replication sample (28 TLBW: 32 matched controls) using the DNeasy Blood & Tissue Kit (Qiagen). The sodium bisulfite conversion was performed with EZ DNA Methylation kits (Zymo Research) following the manufacturer's instructions. Primers were designed using the PyroMark Assay Design software 2.0 (Additional file 4). The target fragment was amplified with PCR. Pyrosequencing reactions were performed using a PyroMark Q96 system (QIAGEN).

2.5. Statistical analysis

Descriptive statistics for demographic information and clinical variables are reported. Continuous variables were reported as mean and standard deviation (Mean ± SD), while categorical variables were expressed as n (%) in frequency and percentage. The Shapiro-Wilk test was used to assess the normality of quantitative variables. For normally distributed continuous variables, we performed parametric t-tests to compare differences between groups. For skewed continuous variables, nonparametric tests were used. The chi-square test was used to measure differences in categorical variables. All statistical analyses were carried out using R software (version 4.1.1).

3. Results

3.1. Description of the study population

The baseline demographic and clinical characteristics of the study participants are summarized in Table 1 and Additional file 1. In the discovery phase, there were seven boys and nine girls in both the TLBW and TNBW groups. The gestational ages of the two groups were comparable (263.6 ± 4.2 days vs. 265.9 ± 3.9 days, p = 0.099). Although parity, mother's age at delivery, mother's height and weight before delivery, maternal pregnancy complications(such as pre-eclampsia and gestational diabetes mellitus), gestational weight gain, prenatal BMI and delivery mode did not show significant differences between the two groups (all p > 0.05), strong differences were found for birth weight (2287.2 ± 154.3 g vs 3026.9 ± 273.6 g), length (46.7 ± 1.1 cm vs 49.8 ± 0.9 cm) and head circumference (31.7 ± 1.2 cm vs 33.3 ± 0.9 cm).

Table 1. Clinical characteristics of the study subjects.

Variables	Methylation 850K BeadChip analysis		p-value	Bisulfite pyrosequencing validation		p-value
	TLBW (n = 16)	TNBW(n = 16)		TLBW (n = 28)	TNBW (n = 32)
Infant characteristics
Gestational age (days), mean ± SD	263.6 ± 4.2	265.9 ± 3.9	0.099	264.1 ± 6.1	265.5 ± 4.5	0.071
Sex, n (%)
Male	7 (43.8)	7 (43.8)	1	11 (39.3)	14 (43.8)	0.726
Female	9 (56.2)	9 (56.2)		17 (60.7)	18 (56.2)
Birth weight (g), mean ± SD	2287.2 ± 154.3	3026.9 ± 273.6	<0.0001	2272.0 ± 170.8	3003.3 ± 264.3	<0.0001
Gravidity, mean ± SD	1.1 ± 0.5	1.6 ± 0.8	0.011	1.3 ± 0.7	1.6 ± 0.8	0.023
Parity, mean ± SD	1.1 ± 0.3	1.3 ± 0.5	0.151	1.1 ± 0.4	1.3 ± 0.6	0.172
Length (cm), mean ± SD	46.7 ± 1.1	49.8 ± 0.9	<0.0001	46.6 ± 1.3	49.7 ± 1.1	<0.0001
Head circumference (cm), mean ± SD	31.7 ± 1.2	33.3 ± 0.9	<0.0001	31.6 ± 1.0	33.5 ± 0.8	<0.0001
Maternal characteristics
Age at delivery(years), mean ± SD	32.5 ± 3.3	34.0 ± 3.8	0.281	32.5 ± 3.3	33.2 ± 4.2	0.494
Maternal height, mean ± SD	161.7 ± 5.9	164.5 ± 3.5	0.074	162.8 ± 5.2	164.5 ± 3.4	0.151
Maternal weight before delivery, mean ± SD	66.4 ± 8.8	69.8 ± 8.0	0.148	67.6 ± 9.7	71.2 ± 10.7	0.138
Gestational weight gain (kg), mean ± SD	10.1 ± 2.7	11.1 ± 4.8	0.637	10.0 ± 3.7	10.8 ± 4.0	0.468
Pre-pregnancy BMI, mean ± SD	25.4 ± 2.9	25.8 ± 3.0	0.874	25.5 ± 3.2	26.3 ± 4.0	0.681
Delivery mode, n (%)
Forceps delivery	1 (6.3)		0.741	2 (7.2)	2 (6.3)	0.925
Spontaneous vaginal delivery	4 (25)			9 (32.1)	9 (28.1)
Caesarean section	11 (68.7)			17 (60.7)	21 (65.6)
Paternal age	36.1 ± 4.5		0.044	33.9 ± 3.8	35.3 ± 4.1	0.128

In the validation stage population, there were no significant differences in gender (p = 0.726) and gestational age (264.1 ± 6.1 days vs 265.5 ± 4.5 days, p = 0.071) between the two groups. Parity, maternal and paternal clinical characteristics were also matched. However, birth weight (2272.0 ± 170.8 g vs 3003.3 ± 264.3 g), body length and head circumference were also significantly lower in the TLBW group compared with the TNBW group.

4. 850k DNA methylation landscapes

After filtering and quality control, we obtained a final set of 812,215 sites that were suitable for subsequent analysis (Figure 2A). The density plot of β values demonstrated a relatively similar overall DNAm distribution for all samples in both the TLBW group and healthy controls (Figure 2B). PCA based on all CpG sites did not reveal any significant separation between the two groups (Figure 2C).

By conducting a comprehensive comparison of genome-wide peripheral blood DNAm patterns between TLBW and TNBW, we identified 144 significant DMPs (Additional file 2). The mean difference in DNAm levels between the two groups ranged from -34.6% to 12%.

The Manhattan plot (Figure 3A) displays the distribution of the identified DMPs across different chromosomes. The points above the red dashed line represent the 144 DMPs that passed our screening threshold. These DMPs are predominantly found on chromosomes 1, 6, 2 and 10. Further annotation of the most significant 10 DMPs revealed seven genes. On chromosome 3, the most differentially methylated genes are KY, ATG7 and SH3BP5, while the other four genes are CEP164 on chromosome 11, CHSY1 on chromosome 15, GADD45G on chromosome 9 and PRKAR1B on chromosome 7. The volcano map (Figure 3B) visually represents these DMPs, with red dots representing 111 hypomethylated sites and green dots representing 33 hypermethylated sites, while gray points represent metabolites that were detected but did not show significant differences. Notably, the 10 most prominent DMPs are cg19280785, cg26155439, cg24594555, cg12120750, cg09530217, cg26452434, cg00772932, cg23605961, cg14554423 and cg19282986.

Figure 3. Differentially methylated position analysis in peripheral blood DNA between TLBW and TNBW group. (A) Manhattan plot. Arrows indicate the genes where the top 10 most significant DMPs were located. (B) Volvano plot. Red spots reveal hypermethylated DMPs (N = 111) with △β≥ 0.05 and p-value < 0.001, green spots represent hypomethylated DMPs(N = 33) with △β≤ -0.05 and p-value < 0.001. (C) Distribution of DMPs according to CpG island contextual distribution. (D) Distribution of DMPs on functional genomic regions.

DMP: Differentially methylated position.

In terms of genomic context, the majority of the DMPs were located in the opensea (69%), followed by shore regions (18%) (Figure 3C). Approximately 23% of the CpG sites were found in the promoter region (TSS200, TSS1500, 5’UTR and 1stExon), 42% in the gene body, and more than 35% in the 3’UTR and intergenic regions (IGR). When examining the distribution of statistically distinct hypomethylated and hypermethylated CpG sites, it was observed that hypermethylated positions were more concentrated in the intergenic regions (46%), while hypomethylated positions were more abundant in the gene body (42%) region, with a smaller proportion in TSS1500 (17%) (Figure 3D). Overall, in the peripheral blood of TLBW neonates, differences in DNAm appear to be predominantly characterized by hypomethylation.

In addition to identifying DMPs, our study also revealed 149 DMRs, as shown in Additional file 3. These DMRs are associated with 147 genes. The top 10 significant DMRs, involving 11 genes, are presented in Table 2. It is noteworthy that 4 of these DMRs are located on chromosome 6.

4.1 Methylation validation

The congruence between the average methylation levels of the four differential DNA methylation sites, as gathered through the Illumina Infinium Human Methylation 850 k Bead Chip and Pyrosequencing, is evident from Figure 4. The pyrosequencing outcomes for these four CpG sites are detailed in Table 3. At the cg23371584 site, the mean methylation rate for the TLBW and TNBW groups stood at 53.725% ± 10.770% and 58.443% ± 7.717%, respectively. The cg16955800 site exhibited methylation rates of 69.594% ± 5.236% and 66.812% ± 5.326% for the TLBW and TNBW groups, respectively. Moving to the cg24312537 site, methylation rates were observed at 32.010% ± 10.655% for TLBW and 36.046% ± 10.344% for TNBW. Finally, the cg27041096 site displayed methylation rates of 48.965% ± 3.941% for TLBW and 51.192% ± 4.534% for TNBW

Table 2. TOP 10 differentially methylated regions.

DMR	CHR	Start	End	Width	Genes	CpGs	Adj. p-value
DMR_2	6	29648161	29649084	923	–	cg22494932, cg08022281, cg08041448, cg10648573, cg07134666, cg11383134, cg03198009, cg13835168, cg19636627, cg25978138, cg15570656, cg02157626, cg03449857, cg12644888, cg04071440, cg25699073, cg24100841, cg11747594, cg15708526, cg20228636, cg16885113, cg00588198	0.0010
DMR_3	6	31868724	31869588	864	ZBTB12	cg27576679, cg08173520, cg20651916, cg22306579, cg11055649, cg11349594, cg27574828, cg11801452, cg03438644, cg06849349, cg24717854, cg22128918, cg11299683, cg24826636, cg11049439, cg20000994, cg09505722, cg13296093, cg08583411, cg26687082	0.0014
DMR_4	18	44561718	44562525	807	KATNAL2, TCEB3B	cg16643261, cg14276258, cg20879768, cg07426472, cg08211738, cg16053920, cg24883732, cg13993676, cg22940153, cg07614190, cg17887685, cg20960383	0.0017
DMR_5	6	31867698	31868207	509	ZBTB12	cg00058449, cg00805874, cg10917426, cg08975528, cg16463880, cg14562426, cg05680710, cg07910050, cg14734916, cg09788778, cg06636203, cg06169557, cg00459243, cg12484688, cg25110523, cg11645762, cg00889295, cg07249939, cg17243044	0.0016
DMR_6	5	78281684	78282099	415	ARSB	cg26802063, cg11541897, cg17845680, cg24893837, cg25318491, cg18856359, cg18828365, cg22224068, cg02582925	0.0023
DMR_7	1	110254603	110254919	316	GSTM5	cg20803293, cg25210835, cg24467349, cg12858902, cg14377951, cg25520068, cg05376982, cg25593510, cg02797771, cg23719124	0.0027
DMR_8	11	111250093	111250744	651	POU2AF1	cg24049888, cg18316498, cg11362935, cg17390301, cg12012195, cg07534554, cg19148410, cg01212326, cg20140034, cg10694914, cg19602649, cg05905611	0.0027
DMR_9	6	33048444	33048919	475	HLA-DPB1	cg19990651, cg01132696, cg06437840, cg25511667, cg12893780, cg02692313, cg10850215, cg15829535, cg17588455, cg25491704, cg14870156, cg02197634, cg05798641, cg00570906, cg14801692, cg13581859	0.0027
DMR_11	12	10182937	10183364	427	CLEC9A	cg26262442, cg17469479, cg02930518, cg20098659, cg19223467, cg11059944, cg13830870, cg22032092	0.0043
DMR_12	2	3642400	3642967	567	COLEC11	cg17872886, cg12135269, cg19867917, cg18105675, cg00835279, cg26615126, cg10945539, cg24755996	0.0044

p-values were calculated using the Benjamini-Hochberg multiple testing correction.

CHR: Chromosome; DMR: Differentially methylated region.

Table 3. The 4 CpGs mean methylation level(%) performed by DNA bisulfite pyrosequencing.

CpG site	TLBW group	TLBW group	Difference	p-value
cg09530217	53.725 ± 10.770	58.443 ± 7.717	-4.719	0.054
cg16955800	69.594± 5.236	66.812 ± 5.326	2.782	0.046
cg24312537	32.010 ± 10.655	36.046 ± 10.344	-4.036	0.142
CG27041096	48.965 ± 3.941	51.192 ± 4.534	-2.227	0.048

5. Enrichment analysis

The statistical difference in DMP-related gene enrichment pathways is not significant and results can be found in Additional file 5. A total of 73 genes associated with DMPs were annotated to 145 GO terms with statistically significant differences. These genes were classified into three major functional categories: 83 related to BP, 40 related to MF and 22 related to CC. The top ten entries in each of the three GO categories are shown in Figure 5A. For cellular components, the genes were significantly enriched in the cell membrane and nuclear chromatin. In terms of molecular function, the genes were significantly enriched in transcriptional regulation and protein binding. The biological process analysis revealed significant enrichment in cytokine production and signaling pathway regulation, among other processes. A bubble plot (Figure 5B) was used to visualize the top 30 enriched pathways analyzed by GO. The most significant enrichment category in BP was “GO:0032743 positive regulation of interleukin-2 production,” while the most significant pathway in MF was “GO:0000980 RNA polymerase II distal enhancer sequence-specific DNA binding.” In the CC category, the most significant enrichment was related to “GO:0005886 plasma membrane.”

Figure 4. The correlation between the methylation levels of 850 k array data and the validation data by pyrosequencing.

Figure 5. Functional annotation analysis of DMPs. (A) Gene ontology (GO) functional annotation. (B) The top 30 GO items pathway enrichment analysis bubble chart. the horizontal axis indicates pathway enrichment score, the vertical axis represents the pathway name and the dot size indicates the number of differential methylation-related genes involved in the pathway. (C) The top 30 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The dot size represents the number of enriched genes, and the dot color represents the p-value.

GO: Gene ontology.

Additionally, 42 genes associated with DMPs were annotated to KEGG pathways. The top 30 significantly enriched KEGG pathways are shown in Figure 5C. Among these pathways, the most important one was the ‘hsa00600 sphingolipid’. Additional file 6 depicts the hsa00600 pathway. Within this pathway, the genes associated with the differential methylation sites are GALC and SGMS1, with the corresponding methylation sites being cg17192520 and cg02638724, respectively. Remarkably, both sites exhibit hypomethylation.

6. Discussion

Epigenetic modifications have emerged as pivotal players in unraveling the intricate interplay between environmental factors and disease onset [25]. It is widely recognized that genetic changes, particularly alterations in DNAm levels, can profoundly impact the expression of key genes by influencing the binding affinity of transcription factors to their target sites, thus contributing to the onset and progression of various diseases [26,27]. Our study presents a landmark investigation, marking the first-ever comprehensive comparison of genome-wide methylation patterns in the peripheral blood of TLBW and TNBW neonates. Notably, our sample selection was carefully balanced with regard to maternal factors that might influence DNAm, ensuring the robustness and reliability of our findings.

For the first time, we identified 144 DMPs and 149 DMRs in the peripheral blood of TLBW neonates. This novel revelation highlights the genes implicated in term low birth weight, notably KY, CEP164, ATG7, CHSY1, GADD45G, PRKAR1B and SH3BP5. These genes encompass a spectrum of cellular processes, ranging from autophagy and cell survival to immune response, potentially substantiating their profound involvement in the development of TLBW. ATG7, a pivotal autophagy effector enzyme, collaborates with other ATG proteins to regulate immunity, cell death, protein secretion and even the cell cycle and apoptosis pathways [28,29]. CEP164 is a key player in the signaling cascade activated by DNA damage [30]. CHSY1 is secreted by the patient's fibroblasts and is required for the synthesis of the chondroitin sulfate moiety. In vitro, elimination of CHSY1 by RNAi knockdown resulted in enhanced osteogenesis of fetal osteoblasts [31]. GADD45G‘s epigenetic silencing in tumors is well-established [32,33]. Evidence suggests that PRKAR1B p.R115K variant affects cAMP signaling and may favorably regulate the lipoprotein profile of African-American youth, protecting them from some adverse metabolic outcomes [34]. The role of these genes in TLBW neonates has not been reported yet. The significant up-regulation or down-regulation of target genes at different methylation sites may eventually lead to the occurrence and development of TLBW, and further experimental verification is needed.

A striking observation was the predominance of hypomethylation, encompassing almost 77% of DMPs, suggesting that reduced DNA methylation levels constitute a hallmark feature of TLBW neonates. In general hypomethylation induces activation whereas hypermethylation induces repression [35]. Global DNA hypomethylation is also a common epigenetic alteration in cancer [36]. Besides that, we found that hypomethylation predominantly affected the gene body and promoter regions, indicating potential ramifications on gene expression modulation during early development.

The GO analysis unveiled that DMPs-related genes were significantly enriched in biological processes involving cytokine production and signaling pathway regulation, notably the positive regulation of IL-2 production. Given that adverse maternal outcomes like preterm birth and low birth weight often results from compromised placental function, such as placental inflammation [37,38], our results resonate with the perturbed methylation of genes governing immune and inflammatory responses in TLBW neonates. Furthermore, molecular functional analysis highlighted enrichment for transcriptional regulation and protein binding, which may contribute to dysregulated gene expression in TLBW neonates.

Our exploration of cellular components highlighted significant enrichment in the cell membrane, while the KEGG pathway analysis spotlighted the sphingolipid metabolism pathway as the most significantly enriched. Lipid mediators of inflammation, such as prostaglandins and leukotrienes, have been extensively investigated and characterized for their roles in immune responses. However, the focus has also extended to another critical class of lipids - sphingolipids. Sphingolipids, integral components of cell membranes and lipoproteins, have garnered significant attention due to their versatile functions. Accumulating evidence underscores the pivotal role of sphingolipids as mediators in the intricate orchestration of apoptosis, cellular stress responses and inflammation regulation [39,40]. Dysregulation of sphingolipid metabolism has been found to be associated with various psychiatric and neurological diseases as well as cognitive impairment [41,42]. The altered metabolism of sphingolipids, marked by augmented synthesis of sphingosine-1-phosphate (S1P), has been associated with the development of the asthmatic phenotype in house dust mite allergy [43]. We also found that both the cell membrane pathway and the sphingolipid metabolism pathway feature the gene SGMS1. The product encoded by SGMS1 serves as the primary sphingomyelin synthase located on the Golgi apparatus, catalyzing the reversible transfer of phosphatidylcholine head groups from phosphatidylcholine to ceramide in sphingomyelin biosynthesis. Our findings accentuate the intriguing possibility of sphingolipid metabolism's involvement in the intricate mechanism underpinning TLBW. In our subsequent studies, we analyzed changes in lipid levels in TLBW, and excitingly, we found that sphingomyelin, Hex3Cer, Hex2Cer and CerAS all showed specific upregulation. These results have not yet been published, and further experimental design is needed to validate them.

While our study offers significant insights into the intricate molecular mechanisms driving the epigenetic regulation of TLBW, it is crucial to acknowledge certain limitations inherent to our research. Foremost among these is that the size of our cohort has increased compared with the single-digit size in the past, and it has been carefully balanced in terms of maternal factors, but because it is a single-center study, it may affect the generalizability of our findings. Furthermore, the cross-sectional design of our study, while valuable in capturing the epigenetic landscape at birth, presents a challenge in establishing a definitive causal relationship between differential DNAm and TLBW.

7. Conclusion

In summary, this study offers a comprehensive exploration of the DNAm landscape in TLBW neonates. Marking a significant milestone, we have identified 144 DMPs between the TLBW and TNBW groups using the human methylation 850k beadchip. Among these differential methylation sites, several have associations with genes involved in crucial processes such as immune response, transcriptional regulation and sphingolipid metabolism. These findings underscore the potential pivotal role of DNAm in the intricate pathogenesis of TLBW. Future investigations are poised to combine diverse omics disciplines, including genomics, metabolomics and proteomics, to provide a more comprehensive perspective on TLBW's epigenetic mechanisms.

Author contributions

Conception and design: J Liu, Q Zhang; administrative support: Q Zhang; provision of study materials or patients: Q Sun, D Liu, HX Liang, F Ye; collection and assembly of data: D Liu, YM Chen; data analysis and interpretation: J Liu, Q Sun; manuscript writing: J Liu, Q Zhang. All authors read and approved the final manuscript.

Financial disclosure

This work was funded by National High Level Hospital Clinical Research Funding (2022-NHLHCRF-LX-01-0301), and the Clinical research project of Beijing research ward construction (2022-YJXBF-04-01-03). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Writing disclosure

No writing assistance was utilized in the production of this manuscript.

Ethical conduct of research

This study was approved by the China-Japan Friendship Hospital ethics committee board (number: 2023-KY-057). All methods were carried out in accordance with relevant guidelines and regulations. Parental informed consent was obtained in writing from each participant.

Competing interests disclosure

The authors have no competing interests or relevant affiliations with any organization or entity with the subject matter or materials discussed in the manuscript. This includes employment, consultancies, honoraria, stock ownership or options, expert testimony, grants or patents received or pending, or royalties.

Data availability statement

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Additional information

Funding

This work was funded by National High Level Hospital Clinical Research Funding (2022-NHLHCRF-LX-01-0301), and the Clinical research project of Beijing research ward construction (2022-YJXBF-04-01-03).