Association Study among Comethylation Modules, Genetic Polymorphisms and Clinical Features in Mexican Teenagers with Eating Disorders: Preliminary Results

Eating disorders are psychiatric disorders characterized by disturbed eating behaviors. They have a complex etiology in which genetic and environmental factors interact. Analyzing gene-environment interactions could help us to identify the mechanisms involved in the etiology of such conditions. For example, comethylation module analysis could detect the small effects of epigenetic interactions, reflecting the influence of environmental factors. We used MethylationEPIC and Psycharray microarrays to determine DNA methylation levels and genotype from 63 teenagers with eating disorders. We identified 11 comethylation modules in WGCNA (Weighted Gene Correlation Network Analysis) and correlated them with single nucleotide polymorphisms (SNP) and clinical features in our subjects. Two comethylation modules correlated with clinical features (BMI and height) in our sample and with SNPs associated with these phenotypes. One of these comethylation modules (yellow) correlated with BMI and rs10494217 polymorphism (associated with waist-hip ratio). Another module (black) was correlated with height, rs9349206, rs11761528, and rs17726787 SNPs; these polymorphisms were associated with height in previous GWAS. Our data suggest that genetic variations could alter epigenetics, and that these perturbations could be reflected as variations in clinical features.


Introduction
Eating disorders (EDs) are severe psychiatric disorders characterized by disturbances of eating behavior, affecting the health and quality of life of individuals. These disorders have an early teenage onset and a hereditary component [1]. EDs have a complex etiology in which genetic and environmental factors interact [2]. Genome-wide association studies (GWAS) and other genetic studies have revealed loci and single nucleotide polymorphisms (SNP) associated with ED [1,3,4]. The clinical characteristics of EDs have been associated in genetic studies. In this sense, significant genetic correlations have been reported in anorexia nervosa with psychiatric disorders, physical activity and metabolic, lipid and Nutrients 2021, 13, 3210 2 of 13 anthropometric traits [5]. In addition, genetic associations between ED and substance use have been described [6]. Additionally, the genetic-environment relationship in ED has been studied through DNA methylation, reporting perturbations in the methylation levels of some genes (DRD2, SLC6A3, POMC, OXTR, among others) [2,7,8]. However, genes do not function alone: on average, each gene is estimated to interact with another four or eight genes, and to be involved with 10 biological functions. Furthermore, recent studies suggest that gene networks provide the potential to identify hundreds of disease-related genes [9]. Analyzing gene-environment interactions in EDs could help us to identify the mechanisms involved in their etiology. Nowadays, new technologies evaluating thousands of genes apply statistic approaches that integrate different information sources from gene interactions (e.g., comethylation module construction) [10]. Comethylation modules are clusters of highly interconnected CpG sites. These modules are detected through the construction of a correlation network. Correlation networks are used to analyze large, high-dimensional data sets. These correlation networks are constructed on the basis of correlations among quantitative measurements (e.g., gene expression profiles, methylation levels) [11]. Comethylation modules are formed by using methylation data as quantitative measurements of gene-environment interactions [10]. Additionally, comethylation modules alleviate various testing problems which are inherent to microarray data analyses, and have been found to be useful for describing pairwise relationships among methylated genes [9][10][11]. In brief, comethylation modules (1) consider all genes as interconnected, (2) identify groups of CpG sites with similar methylation levels, (3) increase statistical power, and (4) detect small effects of epigenetic interactions [9][10][11]. Thus, evaluating correlations among genetic factors, comethylation modules, and clinical features in EDs could be a means by which to identify biological markers in such disorders. The objective of the present study was to detect comethylation modules from DNA methylation samples from children and teenagers with an ED, and to correlate these modules with clinical features and genetic variability.

Study Population
We included 63 subjects diagnosed with anorexia nervosa (AN), bulimia nervosa, (BN) or binge eating disorder (BED) using DSM 5 criteria [12]. Individuals were recruited in the outpatient center of the Children's Psychiatric Hospital "Dr. Juan N. Navarro" from May 2014 to August 2016. Inclusion criteria were subjects with at least three generations of Mexican lineage, 12-18 years of age, and individuals not using psychotropic or psychoactive drugs. The clinical features of the sample are descripted in Table 1. This study followed the principles of the Declaration of Helsinki. Sample recollection and processing were approved by the Ethics Committee of the Children's Psychiatric Hospital "Dr. Juan N. Navarro" with approval No. II3/01/0913 (11 October 2017), and by the Ethics Committee of the National Institute of Genomic Medicine (INMEGEN) with approval No. 06/2018/I.

Evaluation Instruments
BED was screened with the QEWP-R (Questionnaire on Eating and Weight Pattern-Revised) [13]. AN was screened with EAT-26 (Eating Attitudes Test) [14]. We evaluated the presence of psychiatric comorbidity with the Spanish version of MINI Kid (Mini International Neuropsychiatric Interview for Children and Adolescent) [15]. A pedopsychiatrist performed all ED diagnoses.

DNA Extraction and Microarray Analysis
After diagnostic evaluation of each individual, a blood sample was collected using an EDTA tube; DNA was subsequently extracted from this sample. We used the salting-out method from the Gentra Puregene Blood (Qiagen, Germantown, MD, USA) commercial kit. DNA extraction quality and integrity were evaluated by analysis with a NanoDrop spectrophotometer (Thermofisher, Waltham, MA, USA) and 2% agarose gel. DNA samples met the following quality criteria: visible genomic DNA band, 230/260 and 260/230 ratios >1.8, concentration >50 ng/µL, and no signs of DNA degradation. For genotypification, we hybridized DNA with commercial microarray Infinium Psycharray Beadchip (Illumina, San Diego, CA, USA). For methylation analysis, DNA was bisulfite converted using an EZ DNA Methylation Kit (Zymo Research, Irvine, CA, USA). Converted DNA was hybridized with the Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA, USA). Each microarray was processed in the Microarray and Expression Unit of the National Institute of Genomic Medicine.

Quality Control of Genotypification Data
We transformed fluorescence intensities from the Psycharray into genotypes using the GenomeStudio (v. 2.0) software, and quality control was done with the PLINK (v. 1.9) toolset [16]. We eliminated: (1) SNPs with less than 95% genotype calls, (2) individuals with less than 95% genotype calls, (3) individuals with sex discrepancy, (4) SNPs located in chromosomes X and Y, (5) SNPs with less than 0.05 minor allele frequency (MAF), (6) SNPs deviating from Hardy-Weinberger equilibrium (p < 1 × 10 −6 ), and (7) SNPs with A/T and C/G alleles. Subsequently, filtrated data were exported to the R (v. 4.0) software [17], and we removed SNPs with missing data and SNPs without homozygous individuals to minor alleles. Only 193,314 SNPs passed the quality control.

Quality Control of DNA Methylation Data
The fluorescence intensities of the MethylationEPIC microarray were transformed into idat files, which were filtered with ChAMP pipeline (v.2.18) [18] for R (v. 4.0) software. Quality control removed: (1) probes with detection p-value > 0.01, (2) probes with <3 beads in at least 5% of samples per probe, (3) non-CpG probes, (4) multihit probes, (5) probes located in chromosome X and Y, and (6) individuals with sex discrepancy in their genotypification data. We converted filtered methylation data into β-values, which were normalized using the BMIQ (Beta-Mixture Quantile Normalization) method [19]. Afterwards, we evaluated the presence of the batch effect with a singular value decomposition (SVD) method. We preserved CpG sites with standard deviation (SD) > 0.05, keeping 105,393 sites. Likewise, we made many cut points in SD (0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.15, 0.20) to find an optimal point for the construction of comethylation modules.

Comethylation Modules Construction
In order to identify comethylation modules, we processed methylation values with the WGCNA (Weighted correlation network analysis) package [11,20] (R software). Later, we applied the means method to achieve hierarchical clustering, and eliminated individuals with atypical samples. For the final analysis, we only considered 50 subjects. Furthermore, an analysis of network topology was used to determine a soft-thresholding power less than 20 with suitable independence (>0.8) and mean connectivity (<1000). The CpG sites with SD > 0.06 in their methylation values had the best network topology. We applied blockwiseModules function (WGCNA package [11,20]) to detect comethylation modules, using a minimum module size of 175 and a threshold of 20. In a subsequent analysis, we discarded the CpG sites clustered in the grey module, and identified a new set of modules. The 11 constructed comethylation modules included 11,418 CpG sites. Each module was automatically assigned a color by WGCNA, indicating its size. The grey module was discarded from further analysis as it groups unassigned CpG sites to other modules; thus, these sites are unrelated.

Enrichment Analysis of Modules
The CpG sites inside modules were annotated using IlluminaHumanMethylationEPI-Canno.ilm10b4.hg19 package [21]. Genes of CpG sites were extracted and enriched using the WebGestalt online tool [22]. We accomplished the enrichments with Over-Representation Analysis of KEGG Database (Kyoto Encyclopedia of Genes and Genomes) [23]; this was considered to be significant with an adjusted p-value by FDR ≤ 0.05.

Correlation of Comethylation Modules with Clinical Features and SNPs
The eigengene of each comethylation module was correlated with clinical data and SNPs using Pearson's correlations. We calculated the R 2 and p values with cor and corPval-ueStudent functions and set the significant value p < 5 × 10 −3 for clinical data and p < 5 × 10 −8 for SNPs. In order to find associations between SNPs and phenotypes, we used the PheWAS tool on the GWAS Atlas website [24]. Associations were considered to be significant for any phenotype with a p < 1 × 10 −10 .

Description of Comethylation Modules
There were 11 modules in our study. Modules were turquoise (5073 sties), blue (2928 sites), brown (193 sites), yellow (166 sites), green (151 sites), red (150 sites), black (148 sites), pink (145 sites), magenta (135 sites), purple (111 sites), and grey (2218 sites). The CpG sites of these modules were located in 4005 genes. According to relative position to gene, the gene body was the most common annotated location, ranging from 40 sites  Table 2 shows the details of functional annotation of CpG sites from comethylation modules regarding gene location. Concerning the distribution of CpG sites with respect to CpG islands, a majority of comethylation modules corresponded to Open Sea (71.11-84.56%, 85-4207 sites); on the other hand, the purple comethylation module had a high percentage of CpG sites annotated on islands (15.32%, 17 sites) ( Table 3). CpG sites in comethylation modules were heterogeneously distributed among chromosomes (Table S1). Additionally, we observed two groups of comethylation modules given the methylation levels from beta values. One group had partially methylated values (0.2 < β value < 0.8) (turquoise, blue and purple), while the other group had hypermethylated values (β value ≥ 0.8) (brown, yellow, green, red, black, pink and magenta) (Table S2). CpG sites were annotated with the IlluminaHumanMethylationEPICanno.ilm10b4.hg19 [21] package. Data expressed in n of sites, (%) by rows.

Correlations of SNPs with Modules
Seven comethylation modules had correlations with any SNP (brown, green, yellow, magenta, red, black and pink) (Table S4). SNPs were located mostly in intronic regions, ranging from 33.96% in the red module (18 SNPs) to 55.56% in the yellow module (15 SNPs). Another frequent location was intergenic regions, ranging from 14.81% in the yellow module (4 SNPs) to 31.71% in the black module (13 SNPs). The most frequent location in Nutrients 2021, 13, 3210 7 of 13 the green module was intergenic regions (28.13%, 9 SNPs), followed by intronic regions (21.88%, 7 SNPs).
The most correlated SNPs (89.95%, 206 SNPs) were in nonprotein-coding transcript regions, while 10.05% (23 SNPs) were in protein-coding regions (synonymous and missense). Seventeen SNPs (7.42%) were annotated as missense variants; the red and black modules had four missense SNPs each. Meanwhile, six correlated SNPs (2.62%) were annotated as synonymous variants, with two SNPs per module (brown, yellow, and red). We observed 10 correlated SNPs (4.37%) in regulatory regions, although no SNPs were found in the yellow comethylation module. The magenta comethylation module had no SNPs annotated in the upstream and downstream regions. Finally, correlated SNPs annotated in 3 untranslated regions (UTR) were the least frequent (2 SNPs, 0.87%), found within the magenta and red modules ( Table 4).

Correlated SNP PheWAS
Regarding clinical features, BMI, body fat, and height were the most frequent phenotypes associated with SNPs ( Figure 2).
Concerning psychiatric disorders, we found several SNPs to be associated with three comethylation modules. The brown comethylation module had a SNP associated with depressive symptoms and neuroticism (rs4598994). Meanwhile, the pink comethylation module was associated with depressive affect (rs4800995); two SNPs were associated with schizophrenia (rs3129012 and rs356971) in the black comethylation module. Moreover, seven correlated SNPs with four comethylation modules were associated with autoimmune diseases. The magenta comethylation module was associated with rheumatoid arthritis and Crohn's disease (rs1893217). Likewise, the red (rs3095345) and pink (rs9267546 and rs9267547) comethylation modules were associated with rheumatoid arthritis and type 1 diabetes. The black module was associated with primary sclerosing cholangitis (rs3129012 and rs356971), autoimmune vitiligo, and systemic lupus erythematosus (rs356971). Finally, the yellow and black comethylation modules were correlated with clinical features in our population (BMI zscore, conduct disorder, psychotic disorder, and height); likewise, these comethylation modules were correlated with SNPs associated with similar phenotypes. Meanwhile, the yellow module correlated with one SNP (rs10494217) associated with the waist-hip ratio in PheWAS, and the black comethylation module correlated with three SNPs associated with height (rs9349206, rs11761528 and rs17726787). Concerning psychiatric disorders, we found several SNPs to be associated with three comethylation modules. The brown comethylation module had a SNP associated with depressive symptoms and neuroticism (rs4598994). Meanwhile, the pink comethylation module was associated with depressive affect (rs4800995); two SNPs were associated with schizophrenia (rs3129012 and rs356971) in the black comethylation module. Moreover, seven correlated SNPs with four comethylation modules were associated with autoimmune diseases. The magenta comethylation module was associated with rheumatoid arthritis and Crohn's disease (rs1893217). Likewise, the red (rs3095345) and pink (rs9267546 and rs9267547) comethylation modules were associated with rheumatoid arthritis and type 1 diabetes. The black module was associated with primary sclerosing cholangitis (rs3129012 and rs356971), autoimmune vitiligo, and systemic lupus erythematosus (rs356971). Finally, the yellow and black comethylation modules were correlated with clinical features in our population (BMI zscore, conduct disorder, psychotic disorder, and height); likewise, these comethylation modules were correlated with SNPs associated with similar phenotypes. Meanwhile, the yellow module correlated with one SNP (rs10494217) associated with the waist-hip ratio in PheWAS, and the black comethylation module correlated with three SNPs associated with height (rs9349206, rs11761528 and rs17726787).

Discussion
Several studies have evaluated the clinical features, genetic variants, and single DNA methylation sites involved in ED, but none has considered these factors together to date [1,5,7]. As such, little information is available about the integration of these different levels of biological information. As far as we know, this is the first study to integrate clinical features, genetic variants, and DNA methylation using a comethylation network analysis in teenagers with EDs.

Discussion
Several studies have evaluated the clinical features, genetic variants, and single DNA methylation sites involved in ED, but none has considered these factors together to date [1,5,7]. As such, little information is available about the integration of these different levels of biological information. As far as we know, this is the first study to integrate clinical features, genetic variants, and DNA methylation using a comethylation network analysis in teenagers with EDs.
BMI is an important clinical characteristic of the individuals diagnosed with ED, and it has a high impact on metabolic profiles [25]. Low BMI is a diagnosis criterion for anorexia nervosa [12]; on the other hand, bulimia nervosa and binge-eating disorder are related to risk of overweight/obesity [26]. In our analysis, the yellow module could be important in the changes in BMI found in individuals with EDs. The yellow module was correlated with BMI and rs10494217. This SNP is a missense genetic variant changing a histidine aminoacid for an asparagine in position 50 of TBX15 gene (p.His50Asn). Previously, GWAS associated rs10494217 with waist-hip index, a variable related to BMI [27]. The TBX15 gene is a member of the T-box gene family, i.e., transcriptional regulators which play an important role in the development of skeletal elements of limbs, vertebral column, and head, as well as other organs [28,29]. Likewise, this gene was reported as a regulator of metabolism of adipose tissue and muscle fibers, and shown to indirectly regulate body fat and BMI [30][31][32]. TBX15 is highly expressed in adipose tissue, and it binds to the promoter of PRDM16 gene. PRDM16 is essential for the browning of adipose tissue; reduced expression of its protein promotes obesity with high-fat diet and increases visceral fat [33]. As a missense variant, rs10494217 could reduce the binding of TBX15 protein to the promoter of PRDM16, and thus disturb adipose tissue function and alter BMI in individuals diagnosed with an ED. A possible alteration of PRDM16 expression could induce epigenetic reprogramming, as it is found in the yellow module in this study. CpG sites of the yellow module were enriched in alpha-linolenic acid metabolism (PLA2G4E and PLB1) and VEGF signaling pathway (AKT3, NFATC2, PLA2G4E and SHC2); both pathways are involved in adipose tissue function and BMI [34][35][36][37].
The black module could also be related to BMI in ED-diagnosed individuals. This module was correlated with rs11761528, a SNP associated with BMI and androsterone sulfate metabolism [27,38]. rs11761528 is an intronic polymorphism of the ZKSCAN5 gene.
There is little information about ZKSCAN5 (zinc finger with KRAB and SCAN domains 5) and its mechanism; however, animal models suggest that this gene is correlated with adipocyte volume, systolic blood pressure, and cardiac mass [39]. Similarly, rs17726787 was previously associated with height and trunk fat-free mass in GWAS [24,40]. This SNP is an intronic variant of the CELF1 gene. Disturbances in gene expression of CELF1 are related with cardiopathies [41][42][43]. The black module was enriched in mTOR signaling pathway (IGF1R, LPIN1 and RPS6KA2), suggesting that genetic variations like rs11761528 and rs17726787 could alter the epigenetics of this pathway. The mTOR signaling pathway is essential for cardiac development [44,45]. Heart complications are frequent in anorexia nervosa patients, reaching 80% in some studies. Severe anorexia nervosa can change cardiac structure, although most structural abnormalities are reversible [46,47]. Nevertheless, there is a lack of analyses which explore the relationships between altered genes in the module, genetic variation, and cardiopathies in ED-diagnosed individuals.
Although schizophrenia was not correlated in our sample, we found genes and SNPs associated with the disorder. Schizophrenia and other psychiatric disorders are associated with anorexia nervosa [3,5]. Also, there is reportedly a high prevalence of schizophrenia among individuals with eating disorders [48]. CpG sites conforming the black module were enriched in morphine addiction (GABBR2, GABRP and PDE4B), and polymorphisms of the PDE4B gene have been associated with susceptibility to schizophrenia [49]. rs356971 and rs3129012 SNPs were correlated with the black module; these SNPs are associated with schizophrenia [50] and waist-hip index [27]. Furthermore, rs356971 and rs3129012 are associated with phenotypes related with the immunological system, hemoglobin concentration, white blood cells, and platelet count [51]. These SNPs are also associated with primary sclerosing cholangitis [52], autoimmune vitiligo [53], IgA deficiency [54], and systemic lupus erythematosus [55]. Immune-mediated mechanisms have been suggested in the development of EDs; an increased risk for autoimmune diseases in EDs has been reported [56,57]. Likewise, a locus in chromosome 12 was associated with anorexia nervosa, diabetes type 1, and autoimmune diseases [3].
Some modules could have comethylation of CpG sites which was not altered by a genetic effect in individuals diagnosed with an ED, like the turquoise and blue modules. These modules were enriched in pathways associated with the immunological system (Th1 and Th2 cell differentiation, TNF signaling pathway, focal adhesion). The same modules were also enriched in pathways related with development status. The construction of these modules could be influenced by the developmental stage of individuals in the sample, i.e., mainly teenagers [58]. One pathway is the GnRH (Gonadotropin Release Hormone) signaling pathway, which is activated at the beginning of pubertal development, and it depends on neuroendocrine signaling [59][60][61]. Another enriched pathway was associated with adiponectin (adiponectin/CaMKK/AMPK) [62,63]. Many authors suggest that adiponectin levels change with pubertal development [64,65]. Also, partial methylation of these modules suggests transcriptional activation of these pathways. The detection of these modules is more likely to be an effect of background epigenetic alterations and the cell development stage in the tissue (white blood cells) used for the analysis.
Our study has some limitations. Firstly, we note the absence of a control group. However, in this exploratory study, our primary aim was to detect comethylation modules in ED patients and assess the relationship between such modules and clinical phenotypes in ED. Second, we had a small sample with many variables evaluated. Although these conditions could affect the statistic power of our analysis, comethylation modules aggregate covarying CpGs and evaluate grouped CpGs, reducing the number of tests needed. Besides, WGCNA requires at least 20 samples to construct biologically meaningful modules [66].
Third, our data was from a sample made up of Mexican teenagers; therefore, our results should not be applied to all populations with EDs.

Conclusions
This is the first study integrating clinical features, genetic variants, and DNA methylation using comethylation network analysis in teenagers with ED. Our findings showed that two comethylation modules correlated with physical features as well as with SNPs previously associated with metabolic and psychiatric phenotypes. These data suggest that genetic variations could alter epigenetics, and that these perturbations could be reflected as variations in clinical features.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/nu13093210/s1, Table S1: CpG sites per chromosome and normalized value per chromosomic size, Table S2: Descriptive statistics of comethylation modules' beta values, Table S3: Functional enrichments of comethylation modules, Table S4: Characteristics of SNPs correlated with comethylation modules.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author, which were omitted due to privacy and ethical issues.