**1. 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

**Citation:** Nolasco-Rosales, G.A.; Martínez-Magaña, J.J.; Juárez-Rojop, I.E.; González-Castro, T.B.; Tovilla-Zarate, C.A.; García, A.R.; Sarmiento, E.; Ruiz-Ramos, D.; Genis-Mendoza, A.D.; Nicolini, H. Association Study among Comethylation Modules, Genetic Polymorphisms and Clinical Features in Mexican Teenagers with Eating Disorders: Preliminary Results. *Nutrients* **2021**, *13*, 3210. https:// doi.org/10.3390/nu13093210

4

**\***

Academic Editors: Daniel-Antonio de Luis Roman and Ana B. Crujeiras

Received: 9 July 2021 Accepted: 7 September 2021 Published: 15 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 sugges<sup>t</sup> 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–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–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.

### **2. Materials and Methods**

### *2.1. 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.


**Table 1.** Clinical features of study population.


**Table 1.** *Cont.*

Features of subjects who satisfied inclusion criteria. Continuous variables are expressed as mean ± standard deviation. Categorical variables are expressed as *n* (%).

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.

### *2.2. 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.

### *2.3. 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.

### *2.4. 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 × <sup>10</sup>−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.

### *2.5. 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.

### *2.6. 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.

### *2.7. 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.

#### *2.8. 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<sup>2</sup> and *p* values with *cor* and *corPvalueStudent* functions and set the significant value *p* < 5 × 10−<sup>3</sup> for clinical data and *p* < 5 × 10−<sup>8</sup> 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.

### **3. Results**

### *3.1. 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

(56.34%) in the purple module to 2430 sites (69.27%) in the turquoise module. Table 2 shows the details of functional annotation of CpG sites from comethylation modules regarding gene location.


**Table 2.** Comethylation module CpG sites classification.

CpG sites were annotated with the IlluminaHumanMethylationEPICanno.ilm10b4.hg19 [21] package. Data expressed as *n* of sites, (%) by rows. TSS: Transcription Start Site. UTR: Untranslated Region. ExonBnd: Exon Boundaries.

> 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.

### *3.2. Enriched Pathways on Each Module*

We found significant enriched pathways of genes annotated on the CpG sites from turquoise and blue modules (Table S3). Genes in the turquoise module were enriched for the longevity regulating pathway (adjusted *p* value = 0.0047), GnRH (Gonadotropin-releasing hormone) signaling pathway (adjusted *p* value = 0.0042), glioma (adjusted *p* value = 0.0126), cholinergic synapse (adjusted *p* value = 0.0091), human cytomegalovirus infection (adjusted *p* value = 0.0126), and endocytosis (adjusted *p* value = 0.0126). Genes from blue comethylation module were enriched in pathways for Th1 and Th2 cell differentiation (adjusted *p* value = 6.8672 × <sup>10</sup>−7), allograft rejection (adjusted *p* value = 0.0185), endometrial cancer (adjusted *p* value = 0.0111), and TNF signaling pathway (adjusted *p* value = 0.0007). Another enrichment pathways within the same module included the AGE-RAGE signaling

pathway in diabetic complications (adjusted *p* value = 0.0033), phosphatidylinositol signaling system (adjusted *p* value = 0.0074), glioma (adjusted *p* value = 0.0365), longevity regulating pathway (adjusted *p* value = 0.0325), human cytomegalovirus infection (adjusted *p* value = 0.0008), and focal adhesion (adjusted *p* value = 0.0039).

#### *3.3. Correlations of Modules with Clinical Features in Our Population*

In our results, seven clinical features and comorbidities correlated with different comethylation modules (Figure 1). The yellow comethylation module correlated with body mass index (BMI) zscore (R<sup>2</sup> = 0.47, *p* = 0.0006), conduct disorder (R<sup>2</sup> = −0.41, *p* = 0.0030), and psychotic disorder (R<sup>2</sup> = −0.45, *p* = 0.0010). Meanwhile, the purple comethylation module correlated with gender (R<sup>2</sup> = −1, *p* < 1 × <sup>10</sup>−50), suicidal behavior (R<sup>2</sup> = 0.41, *p* = 0.0030), and attention-deficit/hyperactivity disorder (R<sup>2</sup> = −0.59, *p* = 6 ×10−6). Finally, the black module correlated with height (R<sup>2</sup> = 0.4, *p* = 0.0040). Notably, clinical features did not correlate with more than one module at a time.

**Figure 1.** Correlations between modules and clinical features. MDD: major depressive disorder. SB: suicidal behavior. ADHD: attention deficit hyperactivity disorder. Red borders indicate significant correlations. Significance indicates *p*-values < 5 × 10−3.

### *3.4. 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

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).


#### **Table 4.** Annotations of SNPs correlated with modules.

dbSNP codes were annotated with InfiniumPsychArray-24v1-3\_A1\_b150\_rsids file. Coding regions were annotated with Ensembl Variant Effect Predictor. Data expressed as *n* of SNPs, (%) by columns.

### *3.5. 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).

**Figure 2.** Phenotypes associated with correlated single nucleotide polymorphisms (SNP) in PheWAS.

### **4. 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–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–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 sugges<sup>t</sup> 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–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–61]. Another enriched pathway was associated with adiponectin (adiponectin/CaMKK/AMPK) [62,63]. Many authors sugges<sup>t</sup> 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.

### **5. 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 sugges<sup>t</sup> 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/2072-664 3/13/9/3210/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.

**Author Contributions:** Conceptualization, H.N. and A.D.G.-M.; methodology, G.A.N.-R. and I.E.J.- R.; software, G.A.N.-R. and J.J.M.-M.; validation, I.E.J.-R. and C.A.T.-Z.; formal analysis, G.A.N.-R. and J.J.M.-M.; investigation, A.R.G., E.S., G.A.N.-R. and T.B.G.-C.; resources, H.N.; data curation, G.A.N.-R.; writing—original draft preparation, G.A.N.-R.; writing—review and editing, G.A.N.-R., J.J.M.-M., I.E.J.-R., T.B.G.-C., C.A.T.-Z., A.R.G., E.S., D.R.-R., A.D.G.-M. and H.N.; visualization, G.A.N.-R.; supervision, J.J.M.-M. and I.E.J.-R.; project administration, A.D.G.-M.; funding acquisition, H.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Fundación Gonzalo Río Arronte, gran<sup>t</sup> number S591, and Instituto Nacional de Medicina Genómica (INMEGEN), gran<sup>t</sup> number 06/2018/I.

**Institutional Review Board Statement:** This study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committees of the Psychiatric Hospital "Dr. Juan N. Navarro" (protocol code No. II3/01/0913; 11 October 2017), and National Institute of Genomic Medicine (INMEGEN) (protocol code No. 06/2018/I; June 2018).

**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.

**Acknowledgments:** We want to acknowledge Microarray and Expression Unit of National Institute of Genomic Medicine for their technical support. Germán Alberto Nolasco-Rosales and David Ruiz-Ramos are students of Master of Biomedical Sciences degree in Juárez Autonomous University of Tabasco and were supported by CONACYT.

**Conflicts of Interest:** The authors declare no conflict of interest.
