Novel Differentially Methylated Regions Identified by Genome-Wide DNA Methylation Analyses Contribute to Racial Disparities in Childhood Obesity

The magnitude of the childhood obesity epidemic and its effects on public health has accelerated the pursuit of practical preventative measures. Epigenetics is one subject that holds a lot of promise, despite being relatively new. The study of potentially heritable variations in gene expression that do not require modifications to the underlying DNA sequence is known as epigenetics. Here, we used Illumina MethylationEPIC BeadChip Array to identify differentially methylated regions in DNA isolated from saliva between normal weight (NW) and overweight/obese (OW/OB) children and between European American (EA) and African American (AA) children. A total of 3133 target IDs (associated with 2313 genes) were differentially methylated (p < 0.05) between NW and OW/OB children. In OW/OB children, 792 target IDs were hypermethylated and 2341 were hypomethylated compared to NW. Similarly, in the racial groups EA and AA, a total of 1239 target IDs corresponding to 739 genes were significantly differentially methylated in which 643 target IDs were hypermethylated and 596 were hypomethylated in the AA compared to EA participants. Along with this, the study identified novel genes that could contribute to the epigenetic regulation of childhood obesity.


Introduction
In 2016, approximately 340 million children and adolescents between the ages of 5 and 19 years, as well as 40 million children under the age of 5 years, were overweight or obese [1]. From childhood and throughout adulthood, this global epidemic of overweight and obesity poses substantial health risks. Overweight and obese children are a direct consequence of a long-term positive energy balance, which has its proximal determinants in the complex interactions among a person's genetic makeup, lifestyle choices, an environment that promotes obesity, and social factors [2][3][4]. Even if the frequency of obesity among youth has stabilized, it has been significantly higher since the late 1980s for some population subgroups. During the period 1988-1994, only 8.6% of non-Hispanic White girls and 9.7% of White boys aged 2-19 years were obese, whereas obesity among non-Hispanic Black girls (14.5%) and Black boys (14.8%) was much higher. This trend is still continuing, with obesity prevalence being higher among non-Hispanic Black children (24.8%) compared to non-Hispanic White children (16.6%) in 2017-2020 [5]. These statistics signify the fact that obesity has disproportionately affected various ethnic groups, and it is important to understand the factors causing this disparity.

Isolation of Salivary DNA
The Oragene Geno-Tek saliva collection kit (Catalog #OGR-500; Ottawa, ON, Canada) was used to collect saliva. The saliva samples were incubated in a water bath at 50 • C for three hours, as directed by the manufacturer. The PrepIT.L2P DNA isolation kit (Catalog #PT-L2P-5; DNAgenotek, Ottawa, ON, Canada) was used to isolate DNA from a 500 µL aliquot. Each sample was labeled and kept at −20 • C. The quantity of isolated DNA was measured using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Inc., Wilmington, DE, USA).
For further evaluation, the samples were sent to the University of Minnesota Genomics Centre for the Illumina Infinium MethylationEPIC BeadChip method analyses.

Sodium Bisulfite Conversion and Infinium Arrays
Using the EZ DNA methylation kit, sodium bisulfite was used to treat DNA (250-750 ng) (Zymo Research, Irvine, CA, USA). The manufacturer's recommended methodology was followed for measuring DNA methylation using the Illumina Human MethylationEPIC BeadChip (Illumina, CA, USA). The bisulfite-treated samples were amplified, fragmented, purified, and hybridized onto the EPIC BeadChip. Using the Illumina HiScan System, the arrays were cleaned and scanned.
With Illumina's GenomeStudio software V2011.1, raw IDAT data were processed, and background normalization using negative control probes produced methylation β-values that were used for all subsequent analyses. For processing the EPIC data, we used Methyla-tionEPIC v-1-0 B2 manifest.

Data Analysis
From GenomeStudio by Illumina, signal intensities and raw methylation levels were extracted. Probes that had data from just two beads or fewer or those that had signal detection p-values higher than 0.01 were eliminated. To obtain methylation β levels, the signal intensities were standardized, and the noise was eliminated using negative control probes. The methylation β values were derived as the ratio of methylation probe intensity to the overall intensity. IBM SPSS version 24.0 and the R macro-package were used for the statistical analysis. From the 862,927 CpG sites analyzed, a total of 823,645 target IDs were used for further analysis after a quality check. We used a t-test analysis, as shown in Table 1, to measure the difference between weight and category of race. Multinomial hierarchical regression analysis was performed to measure the relationship between the BMI z-score and target IDs, adjusting for covariates (Table S4). To obtain the differences in the methylation status of significant target IDs between NW and OW/OB subjects, an independent sample t-test was conducted. A similar calculation was performed to determine the target IDs that were significantly different between EA and AA subjects. To account for the false discovery rate, a Q-value for each sample was calculated, where a Q-value of 0.05 indicates that 5% of significant tests will result in false positives. Differentially Methylated Regions (DMR) were identified based on the Q-value. The percent difference was calculated by subtracting the average β methylation value of the target ID of OW/OB children and AA children from the average methylation of NW and EA children, respectively, and multiplying it by a hundred. Furthermore, using the list provided by the DisGeNET (https://www.disgenet.org, accessed on 12 August 2022) platform, one of the largest databases for genes and genomes, we chose the target IDs of the 2822 obesity-related genes for further analysis.

Study Participants
For the Infinium MethylationEpic BeadChip Array, we used DNA samples from 31 children: 11 NW and 20 OW/OB. The participants were also divided into different races: 16 EA and 15 AA. Table 1 shows the general characteristics of the study population. The mean age and height were not statistically significant among both groups and, as expected, OW/OB children had higher anthropometric measurements compared to NW children. The body weight (kg) of OW/OB (40.33 ± 12.72) was statistically (p < 0.001) higher compared to NW (27.16 ± 1.89) children. Similarly, the BMI (kg/m 2 ) of OW/OB (2.03 ± 0.10) was significantly greater than that of NW (0.09 ± 0.31). However, between the racial groups (i.e., EA and AA), we did not see any statistical differences in the anthropometric measurements.

Identification of Differentially Methylation Regions
To demonstrate the DMR between NW and OW/OB children and between EA and AA children, we conducted an independent sample t-test between the groups. Figure 1 shows the workflow of the differentially methylated regions (based on the Q-value). A total of 3133 target IDs (associated with 2313 genes) were differentially methylated (p < 0.05) between NW and OW/OB children (Table S1). Among them, 792 target IDs were hypermethylated, and 2341 were hypomethylated in OW/OB children compared to NW participants. From those target IDs, we identified 366 target IDs corresponding to 326 unique genes that are associated with obesity that were differentially methylated between NW and OW/OB children (Table S2). The target IDs that are non obesity related and are differentially methylated between NW and OW/OB children are listed separately in Table  S3. To determine the role of maternal education and family income on the significant target IDs between NW and OW/OB children, we conducted a multinominal hierarchical regression analysis between the obesity-related target IDs and the BMI z-score of the children. Out of the 366 significant obesity-related target IDs, 40 target IDs were not significant after adjusting for maternal education, family income, and race (Table S4). Similarly, for the racial groups EA and AA, a total of 1239 target IDs corresponding to 739 genes (Table S5) were significantly differentially methylated of which 643 target IDs were hypermethylated and 596 were hypomethylated in AA in comparison to EA participants. A total of 116 target IDs were obesity-related corresponding to 91 genes (Table S6), and 1123 target IDs were nonobesity related (Table S7).
To determine the role of maternal education and family income on the significant target IDs between NW and OW/OB children, we conducted a multinominal hierarchical regression analysis between the obesity-related target IDs and the BMI z-score of the children. Out of the 366 significant obesity-related target IDs, 40 target IDs were not significant after adjusting for maternal education, family income, and race (Table S4). Similarly, for the racial groups EA and AA, a total of 1239 target IDs corresponding to 739 genes (Table S5) were significantly differentially methylated of which 643 target IDs were hypermethylated and 596 were hypomethylated in AA in comparison to EA participants. A total of 116 target IDs were obesity-related corresponding to 91 genes (Table S6), and 1123 target IDs were nonobesity related (Table S7).

Methylation Distribution and Classification Analysis
Methylated cytosines can be found in gene bodies, 3' untranslated regions (UTRs), other/open sea areas derived from genome-wide association studies, CpG islands, shores, shelves, open seas, and sites surrounding the transcription sites (200 to 1500 bp, 5' UTR, and exons 1) for coding genes. The regions 0-2 kb from CpG islands is referred to as shores, the regions 2-4 kb from CpG islands are referred to as shelves, and the other/open sea regions are isolated CpG sites in the genome that do not have a specific name [22]. Figure 2 shows the distribution of the genomic location of the significant differentially methylated regions between NW and OW/OB children and between EA and AA children, where most of the target IDs in both groups were islands, 1506 and 292 for BMI and racial group, respectively. Tables 2 and 3 show some clusters of target IDs of the top obesityrelated significant genes with a methylation difference of more than 10%. Some of the top hits for genes differentially methylated between NW and OW/OB children were PCDHGA4, BRD2, PTPRN2, and DIP2C, while those between the racial groups EA and AA were PTPRN2, TNXB, MDLL1, and PRDMI6.

Methylation Distribution and Classification Analysis
Methylated cytosines can be found in gene bodies, 3' untranslated regions (UTRs), other/open sea areas derived from genome-wide association studies, CpG islands, shores, shelves, open seas, and sites surrounding the transcription sites (200 to 1500 bp, 5' UTR, and exons 1) for coding genes. The regions 0-2 kb from CpG islands is referred to as shores, the regions 2-4 kb from CpG islands are referred to as shelves, and the other/open sea regions are isolated CpG sites in the genome that do not have a specific name [22]. Figure 2 shows the distribution of the genomic location of the significant differentially methylated regions between NW and OW/OB children and between EA and AA children, where most of the target IDs in both groups were islands, 1506 and 292 for BMI and racial group, respectively. Tables 2 and 3 show some clusters of target IDs of the top obesity-related significant genes with a methylation difference of more than 10%. Some of the top hits for genes differentially methylated between NW and OW/OB children were PCDHGA4, BRD2, PTPRN2, and DIP2C, while those between the racial groups EA and AA were PTPRN2, TNXB, MDLL1, and PRDMI6.
At the genomic scale, the Manhattan plots show the p-values for the complete epigenomewide association studies (EWAS). The p-values for each target ID are shown in Figure 3a,b by chromosome and position on the chromosome in genomic order (x-axis). The number on the y-axis represents the p-values log10 value, which is equal to one more zero after the decimal point. The Manhattan plots were created using the qqman software in the R-package to illustrate the chromosomal distribution of the significant associations, and volcano plots are shown in Figure 4a,b to visualize the overall size and direction of the associations. From the volcano plots, some of the top hits for the BMI group included cg20242624, cg2230998, and cg19714851, while for the racial groups, the top target IDs were cg03111560 and cg09125754.   Table 3. Cluster of target IDs that were significantly differentially methylated between European American (EA) and African American (AA) children.

No. Gene
No. of Significant Target IDs BAT5 5   At the genomic scale, the Manhattan plots show the p-values for the complete epigenome-wide association studies (EWAS). The p-values for each target ID are shown in Figure 3a,b by chromosome and position on the chromosome in genomic order (x-axis). The number on the y-axis represents the p-values log10 value, which is equal to one more zero after the decimal point. The Manhattan plots were created using the qqman software in the R-package to illustrate the chromosomal distribution of the significant associations, and volcano plots are shown in Figure 4a,b to visualize the overall size and direction of the associations. From the volcano plots, some of the top hits for the BMI group included cg20242624, cg2230998, and cg19714851, while for the racial groups, the top target IDs were cg03111560 and cg09125754.

Pathway Analysis
The extended list of genes with p < 0.05 for NW vs. OW/OB were significantly enriched for Gene Ontology (GO) processes involved in carbohydrate kinase activity (p = 0.000059) and sequence-specific double-stranded DNA binding (p = 0.000388). Table 4 represents the genes involved in each of these processes. While for the genes that are significantly different between EA and AA, GO processes such as cholesterol binding (p = 0.009942), and sterol binding (p = 0.006447) were significantly enriched (Table 5). We also used the same list of genes to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. For the genes that were differentially methylated between the NW and OW/OB groups, pathways such as Hippo signaling pathway (p = 0.000429), folate biosynthesis (p = 0.001785), type II diabetes mellitus (p = 0.004951), and MAPK signaling (p = 0.008373) were enriched, as shown in Table 6. For the gene that was differentially methylated between EA and AA, pathways such as the Ras signaling pathway (p = 0.0126) and biosynthesis of unsaturated fatty acids (p = 0.016499) were enriched (Table 7).

Pathway Analysis
The extended list of genes with p < 0.05 for NW vs. OW/OB were significantly enriched for Gene Ontology (GO) processes involved in carbohydrate kinase activity (p = 0.000059) and sequence-specific double-stranded DNA binding (p = 0.000388). Table 4 represents the genes involved in each of these processes. While for the genes that are significantly different between EA and AA, GO processes such as cholesterol binding (p = 0.009942), and sterol binding (p = 0.006447) were significantly enriched (Table 5). We also used the same list of genes to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. For the genes that were differentially methylated between the NW and OW/OB groups, pathways such as Hippo signaling pathway (p = 0.000429), folate biosynthesis (p = 0.001785), type II diabetes mellitus (p = 0.004951), and MAPK signaling (p = 0.008373) were enriched, as shown in Table 6. For the gene that was differentially methylated between EA and AA, pathways such as the Ras signaling pathway (p = 0.0126) and biosynthesis of unsaturated fatty acids (p = 0.016499) were enriched (Table 7).   Table 4. GO (Gene Ontology) processes that were enriched due to DMRs between NW and OW/OB children.    Lastly, we found 17 obesity-related genes ( Figure 5) that were differentially methylated between NW and OW/OB, as well as between EA and AA, children. Some of these genes were SLC24A3, PLOD2, RPTOR, GH2, RELN, LRRN1, VKORC1L1, and YAP1.   Table 7. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways that were enriched due to DMRs between EA and AA children.

Discussion
In this study, we looked at the DMRs between NW and OW/OB children and also between EA and AA races. Overall, 17 obesity-related overlapping genes were found to be differentially methylated within both groups. Furthermore, the target IDs that were differentially methylated between NW and OW/OB children showed no effects on maternal education and family income.
Along with this, 366 CpG islands were identified for 326 genes that were related to obesity, among which the gene PCDHGA4 (protocadherin γ subfamily A, 4) demonstrated the highest number of target IDs-52 target IDs that were differentially methylated between NW and OW/OB children. Previously, this gene has been associated with dysglycemia [23] but, according to our knowledge, this is the first study that shows its association with childhood obesity. The Brd2 (bromodomain containing 2) deficient model for "metabolically healthy" obesity has garnered attention recently and provides a novel interpretive tool. This model has an epigenetic basis for obesity. Targeted disruption of the Brd2 gene, which is present in the class II major histocompatibility complex (MHC) nearer to Tnf, results in extreme obesity with hyperinsulinemia, as well as hypoglycemia, hyperadiponectinemia, and improved glucose tolerance, which is quite different from other animal models of obesity [24,25]. BRD2 had 18 target IDs that were hypomethylated in the OW/OB children compared to NW and based on previous literature, hypomethylation can cause a decrease in gene expression; therefore, the OW/OB children were hypomethylated. A significant islet autoantigen for type 1 diabetes is encoded by PTPRN2 (protein tyrosine phosphatase, receptor type N2). Previous genetic research has demonstrated a strong correlation between it and obesity. With 12 hypomethylated and 3 hypermethylated sites in OW/OB children, this gene was one of the top hits in this study; unfortunately, the only two CpG islands that previously [26] showed an association between childhood obesity and BMI were not observed in this study, but the results showed novel sites that could play a role in gene expression. We also assessed the DMRs between the two racial groups EA and AA; the gene PTPRN2 had 10 CpG sites hypomethylated in AA children compared to EA. Since for both groups in our study BMI and the racial groups showed clusters of CpGs hypomethylated in the later groups, the results provide us with insight into the potential role of PTPRN2 hypomethylation linked with childhood obesity and racial disparities.
To understand the processes behind these associations, gene-and network-specific functional experiments are needed to determine the significance of pre-pregnancy obesityrelated methylation levels identified in the TAPBP (TAP binding protein) gene. Recently, it was discovered that among five-year-old children, the hypermethylation of numerous CpG sites in TAPBP at birth was linked to an elevated risk of cardiometabolic outcomes (such as insulin sensitivity). Another study demonstrates that hypermethylation of TAPBP regulatory sequences may play a significant role in the association between maternal obesity and cardiometabolic dysfunction in children [27,28]. Our study showed 13 hypomethylated and 1 hypermethylated region in the OW/OB children compared to NW, implying further studies are needed to understand the epigenetic changes related to TAPBP and obesity. According to one study [29], DMR contains the gene FDFT1, which codes for the essential enzyme squalene synthase, which is involved in sterol biosynthesis. Individuals with abdominal obesity have been discovered to have higher visceral fat levels of FDFT1 expression. It has been suggested to treat hyperlipidemia with FDFT1 inhibitors. Furthermore, homozygous knockout mice have demonstrated that FDFT1 is critical for liver function, central nervous system development, and cholesterol production. Children under the age of 10 who breastfed for more than three months have a large DMR across the gene FDFT1. Therefore, breastfeeding may, through methylation, have favorable effects on later life occurrences such as obesity, heart disease, and neurological function. In our study, FDFT1 showed 11 DMRs that were hypomethylated in the OW/OB children. Further studies in obesity and breastfeeding through epigenetic changes in the FDFT1 gene can provide new insight [30][31][32]. Looking at the DMR between EA and AA, TNXB was one of the top hits with nine different target IDs differentially methylated. It has been shown that with various child anthropometric parameters, EWAS have found hundreds of sites that were differentially methylated, with similarities between several studies for the four genes HDAC4, PRLHR, TNXB, and PRDM16 [33][34][35][36][37]. Given that children's obesity outcomes have been linked to the genes PRLHR and TNXB, these genes need further study.
Along with this, the top hits for the DMR between the NW and OW/OB children also included some novel genes that were not identified in previous studies to have been linked to childhood obesity and epigenetics. Those genes were TGIF1, PITX2, COL9A3, PAX6, NKK6-2, SMC4, PCHDB16, TRIM27, BAHCC1, and NFIX, and for the DMR between the EA and AA, genes such as MDLL1, PRDMI6, CSMD1, and TRIM31 were among the top hits, but little to no information was found in the literature.
One of the findings of this study are the enrichment genes associated with various pathways and processes. Divers' pathology and physiological processes, such as tissue repair, wound healing, tissue size, and tissue regeneration, are influenced by the Hippo signaling system, which has been one of the top pathways to be enriched in this study, as well as in previous research [38]. Histone acetylation and deacetylation, dysfunctional miRNAs, abnormal DNA methylation, and aberrant activation of inflammatory genes in the Hippo signaling pathway are all indications of epigenetic alterations [39]. An intervention study [40] observing the effect of a hypocaloric diet on the DNA methylation patterns showed that the MAPK signaling pathway was enriched and hypomethylated in women with a hypocaloric diet, the result of our study falls in line with this, where a cluster of genes with DMRs between NW and OW children were enriched for the MAPK pathway.
Emerging literature shows that the biosynthesis of unsaturated fatty acids is affected by the epigenome. The liver of mice fed a diet high in linolenic acid during pregnancy showed an increase in the DNA methylation status of the Fads2 promoter. In rats, increasing the amount of total maternal fat consumed during pregnancy and breastfeeding caused the Fads2 promoter to remain hypermethylated in the offspring's liver and aorta. However, adult rats that consumed more fish oil experienced momentary, reversible hypermethylation of Fads2. The rodents on high-fat diets also had their placentas and adipose tissue's histone methylation levels changed. Docosahexaenoic acid supplementation during pregnancy caused minor alterations in the overall DNA methylation of leukocytes from cord blood. A high-fat diet changed the DNA methylation state of several genes in young men's skeletal muscle [41][42][43]. The genes that were differentially methylated between EA and AA showed enrichment in the biosynthesis of the unsaturated fatty acids pathway.
Our study has a few limitations as well, such as the relatively small number of participants. Nonetheless, we detected significant epigenetic differences. We did not consider the variations in the genome background of the subjects, which may explain a proportion of the epigenetic variation. Additionally, we want to establish that salivary samples can be used to determine DNA methylation. Previous research shows DNA methylation patterns, both at the gene-specific and global (hydroxyl) methylation levels, are comparable between blood and saliva. The genome-wide DNA methylation profiles of saliva in adults and adolescents are more than 90% similar to those in blood, according to a number of independent studies [44].
To conclude, we found novel genes that were differentially methylated between NW and OW/OB children and between EA and AA ethnicity, which further opens the window to explore the top hits and the overlapping genes in more detail. To our knowledge, this is the first study that explored the epigenomic differences in children who are NW and OW/OB, along with differences in the racial groups EA and AA. Our study also highlights the importance of the consideration of health disparities in epigenetics and obesity, as not many studies have been conducted.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/genes14051098/s1, Table S1: DMRs between NW and OW/OB children; Table S2: Obesity related DMRs between NW and OW/OB children; Table S3: Non obesity-related target IDs between NW and OW/OB children; Table S4: Multinominal hierarchical regression analysis between obesity-related target IDs and BMI z-scores of children adjusted for race, maternal education, and family income; Table S5: DMRs between EA and AA children; Table S6: Obesityrelated DMRs between EA and AA children; Table S7: Non-obesity-related DMRs between NW and OW/OB children.
Author Contributions: Conceptualization, J.R.B., X.W. and T.G.; formal analysis, P.P. and V.S.; funding acquisition, J.R.B., X.W. and T.G.; methodology, P.P. and V.S.; project administration, T.G.; supervision, T.G.; visualization, J.R.B. and T.G.; writing-original draft preparation P.P.; writing-review and editing, all. All authors have read and agreed to the published version of the manuscript. Informed Consent Statement: Informed consent was obtained from all participants involved in the study.

Data Availability Statement:
The datasets used in the current manuscript are available from the corresponding author upon request.