In this study, 87 individuals of Mexican ascendancy were included. Of them, 31 were diagnosed with refractory psychosis and treated with CLZ (CLZ-treated patients), and 56 had early onset psychiatric disorders. This last group of patients has not had any previous psychopharmacological treatment and will be called as drug-naive patients. Samples from both groups of patients were recruited from the outpatient clinics. The CLZ-treated patients with refractory psychosis were recruited from the Instituto Nacional de Neurología y Neurocirugía Manuel Velasco Suárez (INNNMVS), and drug-naive patients from the Hospital Psiquiátrico Infantil Juan N. Navarro and the Hospital General Dr. Gustavo Rivorosa Pérez. At least one specialized psychiatrist performed the diagnosis. Individuals with intellectual disabilities, heavy drinkers, smokers, or substance abusers were excluded. In addition, we excluded drug-naive patients who had begun psychopharmacological treatment. The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Research and Bioethics Committees of INNNMVS (protocol 38/19), 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. Furthermore, all individuals signed informed consents. For those patients under the age of 18 years, an informed assent letter was signed by both parents/guardians and the minor.

Blood DNA extraction, bisulfite-conversion (Zymo, Irvine, CA, United States), and microarray hybridization (Infinium MethylationEPIC BeadChip, Illumina, San Diego, CA, United States) were performed according to the manufacturer’s protocol. We used the Genome Reference Consortium Human Build 37 (GRCh37/hg19) as position reference for all the analyses. The fluorescence intensities were measured with the iScan instrument and transformed into idat files with the algorithm implemented in the GenomeStudio software (Illumina, San Diego, CA, United States). Specialized staff carried out microarray analyses in the INMEGEN. The ChAMP (24) package was used for quality control. Briefly, the probes using the following filters were removed: (i) detection p greater than 0.01, (ii) less than three beads in less than 5% of the samples, (iii) all non-CpG sites, (iv) single nucleotide polymorphism (SNP)-associated probes, (v) probes associated with sex chromosomes, and (vi) multi-hit probes. In addition, samples with a ratio greater than 0.1 were removed, resulting in one individual filtered out. After quality control, a matrix of beta values was constructed.

We normalized the previous beta value matrix by beta-mixture quantile normalization (25) and calculated epigenetic ages with the ENmix package (26). This package allowed us to calculate three epigenetic clocks: Hannum (7), Horvath (8), and DNAm PhenoAge (9). The epigenetic ages calculated with the Horvath, Hannum, and PhenoAge clocks and the chronological age in both groups were compared by correlation analysis. Then, we calculated the delta of age (Δage = epigenetic age − chronological age), subtracting the chronological age from the estimated epigenetic age in all samples and by each clock. Finally, we compared the Δage of each clock between CLZ-treated patients and drug-naive patients. First, stratified analysis by sex and age was carried out, comparing the CLZ-treated patients vs. drug-naive patients. The CLZ-treated patients’ group was later subdivided into subgroups according to the time’s mean of treatment with CLZ (>9 and ≤ 9 years of treatment; i.e., few and many years of treatment, respectively) and according to the mean chronological age (cut-off value = 37 years; i.e., young and older patients). Due to the small sample size, and to achieve subgroups with a similar number of patients (n), the mean values of the years of CLZ treatment and chronological age of patients were used. The student’s t-test was applied for the analysis of these four subgroups.

We removed the batch effect and adjusted blood cell proportion on the normalized beta values matrix with a single-value deconvolution (champ.SVD and champ.runCombat) and a reference base analysis (champ.refbase), respectively (24, 27). Once adjusted, we analyzed differentially methylated sites analysis by linear models implemented in champ.DMP function of the ChAMP package (28). First, we analyzed differentially methylated sites considering significance at a p lower than 1e-08 (a summary of statistics is available in Supplementary Table 1). Then, after identifying differentially methylated sites, a differentially methylated region analysis was carried out with the ProbeLasso algorithm (29) (a summary statistics of the differential methylation region analysis is available in Supplementary Table 2).

We performed pathway enrichment after extracting genes resulting from differentially methylated sites and regions through the WebGestalt platform (WEB-based GEne SeT AnaLysis Toolkit) with free access (30) (all pathways of this enrichment are available on Supplementary Table 3). First, we analyzed protein–protein interactions between gene products identified using the web application String (31) and plotted the interactions with Cytoscape 3.8.1 software (32). Then, a sub-enrichment analysis of the longevity regulating pathway was performed in String using the KEGG pathway database (33).

In this study, eighty-seven individuals were included and classified into two groups: the CLZ-treated patients with refractory psychosis (n = 31), and the drug-naive patients with psychiatric disorders (n = 56) (Table 1). The following comorbidities were observed in CLZ-treated patients: diabetes (22.6%), epilepsy (16.1%), hyperlipidemia (12.9%), hypercholesterolemia (6.4%), hypothyroidism (6.4%), hypertension (9.7%), and obsessive-compulsive disorder (3.2%). Details about the comorbidities found in the naive patients have been previously reported (34).

The variable sex did not show a significant correlation with epigenetic age; however, it was considered in the model as previous reports have demonstrated differences in the epigenetic age between men and women (23). The years of CLZ treatment correlated with epigenetic age in the PhenoAge clock (t = 2.6041, r² = 0.4415, p = 0.01458); thus, these two variables were included in the stratification of the patients. The chronological age positively correlated with the epigenetic age in all clocks, Hannum (r² = 0.679, p = 4.96e-13), Horvath (r² = 0.662, p = 2.85e-12), and PhenoAge (r² = 0.664, p-value = 2.26e-12) (Figure 1). After obtaining the Δage for all samples, we found that the CLZ-treated patients’ group had a smaller Δage compared with drug-naive patients (Figure 2). This difference was significant for the three analyzed clocks: Horvath (t = −3.23, p = 0.0018), Hannum (t = −4.83, p = 6.0740e-06), and PhenoAge (t = −2.72, p = 0.0079). Hannum’s clock showed the most significant value. The mean Δage in the CLZ-treated patients was 13.64, 9.65, and 1.54 years for Horvath, Hannum, and PhenoAge clocks, respectively. In drug-naive patients, the mean Δage was 20.77, 18.84, and 8.15 years in the same mentioned order for the three clocks.

The stratified analysis by sex and age, after comparing CLZ-treated patients vs. drug-naive patients, showed significant differences in the three clocks only for CLZ-treated male patients with mean Δage values of 13.15, 10.07, and 1.60 years for the Horvath, Hannum, and PhenoAge clocks, respectively, showing corresponding values of p of 0.0376, 0.0073, and 0.0370 (Supplementary Figure 1A). The mean Δage in CLZ-treated female patients was 14.30, 9.06, and 1.47 years for the Horvath, Hannum, and PhenoAge clocks, respectively. However, significant differences were only found in Hannum’s clock (t = −3.2253, p = 0.0035). Thus, we did not identify differences between age subgroups (Supplementary Figure 1B).

We explored Δage differences in CLZ-treated patients, subdivided according to the mean value of chronological age (37.74 years) and the mean treatment time with CLZ (9.22 years). Significant differences of Δage for the Horvath and Hanum clocks were observed between patients by years of CLZ treatment and chronological age (younger vs. older patients, and few years of treatment vs. many years of treatment, respectively). In addition, other significant differences were found for the Hannum clock in the following comparisons: older patients with more years of treatment (>9.22 years) vs. younger patients with more years of treatment (p = 0.0021); older patients with many years of treatment vs. young patients with few years of treatment (p = 0.0263); and older patients with few years of treatment vs. young patients with many years of treatment (p = 0.0236) (Figure 3A). Regarding Horvath’s clock, we only observed differences between the group of older patients with many years of treatment and younger patients with many years of treatment (p = 0.019) (Figure 3B).

We identified 44,716 differentially methylated sites between CLZ-treated patients and drug-naive patients. Of note, 39,204 CpG sites (87.70%) were hypomethylated in the CLZ-treated patients, and the remaining (5,512 CpG sites, 12.30%) were hyper-methylated.

According to the position of the differentially methylated sites in relation to the gene structure, the following distribution was observed: 9,573 (24.42%) hypomethylated sites were located at 200 nucleotides of the transcriptional start site (TSS200), 7,247 (18.49%) sites on TSS1500, 7,126 (18.18%) on gene body, 5,332 (13.6%) on 5′-UTR, 4,626 (11.8%) on first exon, 306 (0.78%) on 3′-UTR, and 22 (0.06%) on exon binding sites. The distribution of the hypermethylated sites was as follows: 3,111 (56.44%) sites were annotated to gene body, 266 (4.83%) on 5′-UTR, 265 (4.81%) on TSS1500, 262 (4.75%) on 3′- UTR, 138 (2.5%) were located on the TSS200, 75 (1.36%) on first exon, and 55 (0.1%) on exon binding sites.

The location of differentially methylated sites concerning CpG islands was as follows: 29,895 (66.86%) sites on CpG islands, 7,860 sites (17.56%) on shores, 5,808 sites (12.99%) on the open sea, and 1,153 sites (2.58%) on shelves. Regarding hypomethylated sites, 29,020 of them (74%) were found on CpG islands, 6,822 (17.4%) on shores, 2,892 sites (7.38%) on the open sea, and 470 (1.2%) sites on shelves. The distribution of the hypermethylated sites was: 2,916 (52.9%) sites located on the open sea, 1,038 (18.83%) sites on shores, 875 (15.87%) sites on CpG islands, and 683 (12.39%) sites on shelves.

In total, sixty-five pathways enriched at the hypomethylated sites were identified, and the main ones were involved in the biosynthesis of glycosaminoglycan, glioma, longevity regulatory pathway, VEGF signaling pathway, and circadian rhythm. Additionally, the hypermethylated sites are enriched on ABC transporters and endocytosis pathways (Table 2). Supplementary Table 2 contains the enrichment of all the observed pathways. Among the enriched pathways at the hypomethylated sites, the longevity pathway was identified, which, in turn, interacts with others, such as the AMPK and the insulin signaling pathways. The longevity regulating pathway occupied the 21st place of the ranked list generated by the pathway enrichment; it showed an overlap of 74 genes out of an expected 89 [p = 1.25e-7, false discovery rate (FDR) = 1.14e-6, and enrichment ratio = 1.4578]. We sub-analyzed the longevity regulating pathway because CLZ has been associated with a reduced mortality risk compared with other antipsychotics. For this purpose, a sub-enrichment analysis using protein–protein interaction on this pathway (Figure 4) was performed.

A total of 286 differentially methylated regions were identified, of which 247 regions were hypomethylated and 39 hypermethylated. Then, the FUMA platform [Functional Mapping and Annotation of Genome-Wide Association Studies (GWAS)] (35) was used to prioritize and interpret GWAS results of the 307 annotated genes in this region. This analysis revealed top genes, such as NOTCH4, MICA, TRIM27, PBX2, and FKBPL, as enriched in GWAS of schizophrenia, CLZ-induced agranulocytosis/granulocytopenia in treatment-resistant schizophrenia, and bipolar disorder among other psychiatric and non-psychiatric conditions (Supplementary Figure 2). In this enrichment analysis of the differentially methylated regions, we did not find any pathway with a statistically significant association; nevertheless, in the nominal association, twelve pathways were identified (Supplementary Figure 3).