Genome-Wide DNA Methylation Analysis in Male Methamphetamine Users With Different Addiction Qualities

This paper aimed to explore the genome-wide DNA methylation status of methamphetamine (MA) abusers with different qualities to addiction and to identify differentially methylated candidate genes. A total of 207 male MA abusers with an MA abuse frequency of ≥10 times and an MA abuse duration of ≥1 year were assigned to the high MA addiction quality group (HMAQ group; 168 subjects who met the diagnostic criteria for MA dependence according to the DSM-IV) or to the low MA addictive quality group (LMAQ group; 39 subjects who did not meet the criteria for MA dependence). In addition 105 healthy controls were recruited. Eight HMAQ subjects, eight LMAQ subjects, and eight healthy controls underwent genome-wide DNA methylation scans with an Infinium Human Methylation 450 array (Illumina). The differentially methylated region (DMR) data were entered into pathway analysis, and the differentially methylated position (DMP) data were screened for candidate genes and verified by MethyLight qPCR with all samples. Seven specific pathways with an abnormal methylation status were identified, including the circadian entrainment, cholinergic synapse, glutamatergic synapse, retrograde endocannabinoid signaling, GABAergic synapse, morphine addiction and PI3K-Akt signaling pathways. SLC1A6, BHLHB9, LYNX1, CAV2, and PCSK9 showed differences in their methylation levels in the three groups. Only the number of methylated copies of CAV2 was significantly higher in the LMAQ group than in the HMAQ group. Our findings suggest that the circadian entrainment pathway and the caveolin-2 gene may play key roles in MA addiction quality. Further studies on their functions and mechanisms will help us to better understand the pathogenesis of MA addiction and to explore new targets for drug intervention.


INTRODUCTION
Methamphetamine (MA) is an amphetamine-type stimulant drug with a strong propensity for abuse, and it promotes the release of monoamine neurotransmitters in the central and peripheral systems (1). A low dose of MA can induce a series of short-term effects, including euphoria, an increase in energy, positive emotions, increases in heart rate and blood pressure, mydriasis, peripheral high fever, a loss of appetite and weight, improvements in mental acuity, and social and sexual disinhibition (2). A high dose or repeated abuse of MA can cause mental disorders with typical symptoms, including paranoid delusions, auditory hallucinations, increased activity, and stereotypic behavior (3).
To date, MA abuse remains a serious social problem worldwide. The report of the United Nations Office on Drugs and Crime (UNODC) shows that amphetamine-type stimulants, especially MA, rank as the second most widely used illegal drugs worldwide, and 0.7% of the world's population used amphetamine-type stimulants in 2016 (4). According to the 2017 Report on Drug Control in China, by the end of 2017, more than 1.54 million individuals abused synthetic drugs (especially MA), accounting for 60.2% of drug abusers, which did not include double hidden drug users (5). Numerous clinical observations have indicated that different individuals experiencing the same level of repeated MA exposure show different qualities to addiction. A large sample epidemiological survey also exhibited the same situation. In 2012, the Substance Abuse and Mental Health Services Administration (SAMHSA) estimated that 12 million individuals (∼4.7% of the population aged 12 years or older) had used MA once in their lifetime, 1.2 million individuals (∼0.4% of the population aged 12 years or older) had used MA in the past year, and 0.44 million individuals (∼0.2% of the population aged 12 years or older) had used MA in the past month (6). Moreover, 0.54 million individuals (also ∼0.2% of the population aged 12 years or older, roughly close to the proportion of the individuals using MA in the past month) in the United States in 2012 met the diagnostic criteria for stimulant drug abuse or dependence according to the DSM-IV (6).
However, the mechanisms for differences in quality to addiction remain unclear. Over the years, classic DNA polymorphism studies have not clearly resolved this issue. Most of the studies have focused on the potential susceptible genetic variations for addiction and have found that genetic polymorphisms of multiple candidate genes are significantly associated with susceptibility to addiction (7). It should be noted, however, that the currently identified risk locus associated with addiction susceptibility could account for only a moderate degree of the genetic variance (<5%) (8). Moreover, the results of molecular genetic studies are non-repeatable because of the small sample size and significant heterogeneity. In general, differences in quality to addiction among different individuals are currently thought to arise from genetic predisposition and the environment (e.g., acute drug exposure). At present, epigenetic changes are generally considered to be manifestations of geneticenvironment interactions. Epigenetic refers to heritable but reversible regulation of various genetic functions, including gene expression mediated by modifications of DNA and histones (9). Various aspects of epigenetic changes are very stable, making them ideal vectors for addiction susceptibility research (10).
Disturbances in epigenetics, particularly DNA methylation, have been implicated in the pathophysiology of several psychiatric disorders, such as eating disorders (11), depression (12), schizophrenia (13), cocaine addiction (14), and alcohol dependence (15). Epigenetic research on MA addiction is at a relatively early stage, far behind that of cancer or other complex disorders, and far behind that of cocaine, heroin, and other substance abuse. A review discussed recent advances in epigenetic mechanisms underlying amphetamine-or MAinduced behavioral, transcriptional, and synaptic plasticity, and found drug exposure induces major epigenetic modifications in a very complex manner, such as histone acetylation and methylation, DNA methylation (16). DNA methylation is an important gene regulatory mechanism underlying MAinduced learning and memory alterations (17). The genomewide epigenetic approaches will accelerate the speed of discovery in the field of addiction. However, there has only a few genome-wide DNA methylation study on MA. Desplats et al. (18) reported genome-wide profiling of DNA methylation in postmortem frontal cortex tissue of with HIV-1 infection and MA dependence patients showed differential methylation on genes related to neurodegeneration; dopamine metabolism and transport; and oxidative phosphorylation. Cadet et al. (19) reported genome-wide DNA hydroxymethylation analysis of rats with MA exposure indicated that changes in differentially hydroxymethylated regions and increased expression of specific potassium channels in the NAc may promote abstinence from drug-taking behaviors. González et al. (20) reported that the global DNA cytosine methylation (5-mC) levels in genomic DNA were significantly increased in the medial prefrontal cortex (mPFC) of mice treated with a 7-days repeated MA administration.
This study aimed to explore the genome-wide DNA methylation status of individuals with different qualities to MA abuse through a two-step experiment, including screening for candidate genes with genome-wide DNA methylation microarrays and verifying abnormally methylated genes from a large number of samples using MethyLight quantitative polymerase chain reaction, which were divided into samples from high MA addictive quality and low MA addictive quality groups.

Participants
Two hundred and seven male MA abusers, with an MA abuse frequency of ≥10 times and an MA abuse duration of ≥1 year from 577 volunteer participants of Xinkaipu Compulsory Rehabilitation Institute (Changsha, Hunan) and White Mud Lake Compulsory Rehabilitation Institute (Xiangyin, Hunan), underwent a Structured Clinical Interview for DSM-IV Axis I Disorders (SCID-I/P) from March 2013 to January 2014. Compulsory rehabilitation is the primary form of treatment for illegal drug dependence (61.6% drug users received compulsory treatment by 2012) in China (21). The residential MA users were admitted in the treatment facilities after confirmed MA use by urine or hair test. Among them, 168 individuals were classified into the high MA addictive quality group (HMAQ group; MA-dependent patients according to the DSM-IV), and 39 individuals were classified into the low MA addictive quality group (LMAQ group; did not meet the diagnostic criteria for MA dependence). In addition, 105 healthy controls who received physical examinations at the Second Xiangya Hospital were recruited.
A complete medical history, physical examination and laboratory tests were obtained from all subjects. All of them were male, Han Chinese, aged 18-50 years, and without severe acute or chronic medical illnesses. The demographic data and substance abuse characteristics of the three groups are shown in Table 1. Polydrug abuse in methamphetamine users is very common (22), and this paper wants to focus on methamphetamine users. So, the subjects with other illicit drugs use more than three times, including but not limited to heroin, marijuana, and Ketamine, were excluded from this study. There were only 207 male MA abusers met the requirements in 577 volunteer participants from two compulsory rehabilitation institutes. The co-use of alcohol and tobacco were recorded in Table 1.
This study was conducted in accordance with the principles expressed in the Declaration of Helsinki and was approved by the Ethics Committee of the Second Xiangya Hospital. Prior to the study, the procedure was fully explained, and written informed consent was obtained from each subject. All the subjects were free to participate in or abstain from this study, and free to withdraw from this study anytime without threat of punishment.

Sample Processing, DNA Extraction, and Bisulfite Conversion
Fasting blood samples were collected and stored at −80 • C. Genomic DNA was isolated from blood samples using DNeasy blood & tissue kits (Qiagen, Valencia, CA, USA), and bisulfite modification was performed using the EZ DNA Methylation-Gold kit (Zymo, Irvine, CA) after all blood samples were collected. The bisulfite-converted DNA was resuspended in TE buffer and stored at −80 • C until analysis.

Illumina Infinium ® HumanMethylation450 BeadChip Microarray Analysis
We randomly selected eight subjects (dosage: 675.22 ± 93.26 mg/time) in the HMAQ group, eight subjects (dosage: 384.32 ± 87.75 mg/time) in the LMAQ group, and eight subjects in the healthy control group for genome-wide DNA methylation analysis. The bisulfite-converted DNA was then hybridized to the Infinium HumanMethylation450 BeadChip microarray (Illumina, San Diego, USA) in Shanghai Biostar Genechip, Inc. (Shanghai, China) using Illumina-supplied reagents and protocols. This method simultaneously profiles the methylation status for > 485,000 CpG sites at a single-nucleotide resolution. The methylation status of each individual CpG locus was calculated as the ratio of the fluorescent signals (β = Max (M,0)/[Max (M,0) + Max (U,0) + 100]), which ranged from 0 (no methylation) to 1 (complete methylation), using the average probe intensity of the methylated (M) and unmethylated (U) alleles. The quality control information is provided in the supplement as Supplementary Figures 1-3. The original data of the chip were preprocessed using the ChAMP package of R software with the BMIQ method for type 2 probe calibration, the ComBat method of the SVA package for the correction of multiple batch effects, the RefbaseEWAS for celltype heterogeneity calibration, the minfi package for functional normalization function analysis, and the Lasso method for differentially methylated region (DMR) analysis. The screening of differentially methylated loci screening used the pooled t-test with p < 0.05 and |beta.difference|>0.14 as the threshold for significant differences.

Functional Annotation and Pathway Analysis and Potential Candidate Gene Selection
After the genome DNA methylation chip data were analyzed, the differentially methylated positions (DMPs) and differentially methylated regions (DMRs) were screened. DMRs were included in the pathway analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.ad.jp/kegg/), and DMPs were used to screen for candidate genes to perform MethyLight qPCR verification with large samples. To screen for candidate genes from the complex DMP results and to eliminate the impact of confounding factors, the following screening strategies were performed: first, the abnormally methylated sites should correspond to the specific gene name or the UCSC_REFGENE_ACCESSION number; second, the abnormally methylated sites should be located in TSS200 and TSS1500 of the gene promoter region; third, Gene Ontology (GO) analysis (http://www.geneontology.org) should confirm the function of the abnormally methylated genes related to the nervous system; and fourth, the methylation statuses of the candidate gene should be different in all pairwise comparisons of the three groups.

Verification of DNA Methylation With MethyLight qPCR
MethyLight qPCR was used to verify the methylation status of the candidate genes, and it is based on the principle of fluorescence-based real-time PCR with two probes. We prepared standard gene products of SLC1A6, BHLHB9, LYNX1, CAV2, and PCSK9 by the amplification of target DNA fragments and the purification of agarose gel DNA fragments using a Gel Extraction Kit and DNA Methylation Kit (CWBiotech, Beijing, China). Their primers and probes were designed by Beacon Designer 7.0, as shown in Table 2. MethyLight qPCR was performed on a PikoReal 96 Real-Time PCR System (Thermo Fisher Scientific Inc., US) using bisulfite-treated DNA and a GoldStar TaqMan Mixture (CWBiotech, Beijing, China). The MethyLight qPCR

Statistical Analysis
The chip data and MethyLight qPCR data were processed by R version 3.2.4 (R Statistics), and raw-intensity files were preprocessed and normalized together in quantiles using the ChAMP package with the BMIQ method for type 2 probe calibration, the ComBat method of the SVA package for the correction of multiple batch effects, the RefbaseEWAS for celltype heterogeneity calibration, the minfi package for functional normalization function analysis, and the Lasso method for differentially methylated region (DMR) analysis. The screening for differentially methylated loci used the pooled t-test with p < 0.05 and |beta.difference|>0.14 as the threshold for significant differences. The remaining data were analyzed using SPSS 19.0 software with a t-test, chi-squared test, Pearson correlation analysis, and covariance analysis with age as the covariant. The discrete trend of CAV2 methylation level in LMAQ group was very serious after the outliers were winsorized at level 1 and 99% (The sample size of LMAQ group is 39, only the maximum value was assigned as 99%). In covariance analysis, the CAV2 methylation level was adjusted with age as covariance from 1.05E6 ± 9.83E5 to 1.05E6 ± 8.51E4, and its discrete trend was relieved.

Analysis of Genome-Wide DNA Methylation Microarray Data in Patients With Different Addictive Qualities
After genomic DNA methylation chip data were analyzed, the genome wide methylation level in human chromosomes was shown in Figure 1. The distribution of differential methylation sites on chromosomes was shown in Supplementary Figure 9.
The heatmaps and volcano plots of differential methylation sites in pairwise comparisons of three groups were shown in Supplementary Figures 10, 11. All the chip data had  Figure 4. After the non-specific pathways were removed, there were seven pathways associated with the nervous system or neuronal function in the KEGG pathway analysis results of the methylation chip data, as shown in Table 3 and Figure 5.
After the DMP results were screened following the strategy described in the Methods section, 78 significantly differentially methylated genes remained. There were 32 significantly differentially methylated genes (15 genes with downregulated methylation levels and 17 genes with upregulated methylation levels) between HMAQ patients and healthy controls, 21 significantly differentially methylated genes (12 genes with downregulated methylation levels and nine genes with upregulated methylation levels) between LMAQ patients and healthy controls, and 25 significantly differentially methylated genes (14 genes with downregulated methylation levels and 11 genes with upregulated methylation levels) between HMAQ patients and LMAQ patients, as shown in Figure 6.
After further literature review and gene functional analysis, SLC1A6, BHLHB9, LYNX1, CAV2, and PCSK9 were selected for MethyLight qPCR verification with a large sample, and they all had significant differences in the pairwise comparisons of methylation chip data among the HMAQ group, LMAQ group, and healthy control group.
Comparison of Methylation Levels of SLC1A6, BHLHB9, LYNX1, CAV2, and PCSK9 in the Three Groups MethyLight qPCR verification of DNA methylation results of the five target genes indicated significant differences in the number of methylated copies of CAV2, BHLHB9, LYNX1, PCSK9, and SLC1A6 among the three groups (Table 4 and Figure 7). Pairwise comparisons indicated that the number of methylated copies of CAV2 was significantly higher in the LMAQ group than in the HMAQ group, while it was significantly lower than that in the healthy control group. The number of methylated copies of BHLHB9, LYNX1, PCSK9, and SLC1A6 was not significantly different between the HMAQ group and the LMAQ group. However, the number of methylated copies of BHLHB9 and LYNX1 was significantly higher in both the HMAQ group and the LMAQ group than in healthy control group, while the number of methylated copies of PCSK9 and SLC1A6 was significantly lower in both the HMAQ group and the LMAQ group than in the healthy control group.

Correlation Analysis of Candidate Gene Methylation Levels With MA Abuse Characteristics in the Three Groups
The number of methylated copies of the CAV2 gene was significantly associated with age and the total duration of drug abuse in the LMAQ group. The number of methylated copies of the BHLHB9 gene was significantly correlated with the age of first drug abuse and the total duration of drug abuse in the LMAQ group. No other correlation was identified in the HMAQ group or the healthy control group ( Table 5).

DISCUSSION
To the best of our knowledge, this is the first study to examine genome-wide DNA methylation in patients who abuse MA. The main findings of this study included (1) the important role of the circadian entrainment pathway in the pathogenesis of MA addiction inferred by the analysis of DMR data and (2) the important role of CAV2 in different qualities mechanisms of MA addiction revealed by DMP data.

Analysis of DMR Data Inferred the Role of the Circadian Entrainment Pathway in the Pathogenesis of MA Addiction
In this study, all DMRs of the genomic DNA methylation chip data in the pairwise comparisons of the three groups were included in the KEGG pathway analysis, which showed that seven specific pathways were statistically significant. Moreover, all seven pathways, including the Circadian entrainment (23), Cholinergic synapse (24), Glutamatergic synapse (25), Retrograde endocannabinoid signaling (26), GABAergic synapse (27), Morphine addiction (28), and PI3K-Akt signaling pathways (29), were identified to be associated with MA addiction or reward circuits.
The most significant finding of our current study was that the circadian entrainment pathway, which plays an important role in helping organisms cope with dynamic environmental changes and maintain survival with adaptive behavior, was involved in susceptibility to MA addiction (30). In mammals, the circadian entrainment pathway consists of two parts. One is a light signal mechanism that converts light signals into glutamic acid and pituitary adenylate cyclase activating peptide (PACAP) signals and stimulates the suprachiasmatic nucleus (SCN) of the hypothalamus to form a circadian rhythm. Then, the SCN connects the auxiliary oscillators of the whole brain and the body through peptide and neurotransmitter signals (including dopamine, glutamic acid, gamma-aminobutyric acid, etc.) and regulates the physiology and behavior of organisms in a 24-h cycle with molecular clocks throughout various cells throughout the body (31). Second, in the non-optical signaling mechanism, melatonin inhibits light-induced circadian rhythms by binding to melatonin receptors (MT1, MT2) to inhibit adenylate cyclase (AC) (32).
The circadian rhythm transcription factors are also deeply involved in the mechanism of the reward pathway. The SCN directly sends gamma-aminobutyric acid (GABA)-ergic projections to connect the lateral habenular nucleus (LHb), dorsomedial hypothalamic nucleus (DmH), medial anterior optic area (mPOA) and paraventricular thalamic nucleus (PVT), which can affect reward circuits through indirect neural connections (33). The crucial brain structural regions in the reward circuits, such as the ventral tegmental area (VTA), nucleus accumbens (NAc), prefrontal cortex (PFC), hippocampus (Hipp) and amygdala (Amy), have been shown to express proteins related to clock regulation (circadian transcription factors) (34) and nonclock regulation (melatonin receptors) (32) to maintain circadian rhythm stability and control emotion and reward behavior through neural circuits and molecular mechanisms (35).
Further, circadian rhythm transcription factors play an important role in the pathogenesis of substance addiction disorder. First, genome-wide association studies (GWASs) have shown close relationships between polymorphisms and other mutations in the circadian rhythm core genes and the pathogenesis of addiction disorder (35). Second, circadian rhythm disorders increase susceptibility to addiction. Repeated light/dark conversion can increase alcohol intake (36) and the preference for MA (37) in pre-exposure-treated rats. In humans, epidemiological studies reported that shift work or jet lag was associated with increased smoking and alcohol consumption (38,39). Third, substance abuse can alter the functional state of the circadian entrainment pathway. Acute or chronic substance abuse in rodents has been found to alter the expression of circadian rhythm genes, including per1, per3, cry1, bmal1, npas2, and clock, in brain areas such as the NAc, dorsal striatum and hippocampus (40). Finally, manipulation of circadian rhythm genes can significantly alter dopaminergic signaling and rewardrelated behavior. Clock 19 mice showed a high dopaminergic state characterized by increased tyrosine hydroxylase expression in the VTA region and increased dopamine release in the NAc region, with strong sensitization to cocaine, increased preference   | The list of differential DNA methylation genes screened out from methylation chip data. After the DMP results were screened following the (Continued) FIGURE 6 | strategy described in the Methods section, 78 significantly differentially methylated genes remained. Totaling 32 significantly differentially methylated genes (15 genes with downregulated methylation levels and 17 genes with upregulated methylation levels) between HMAQ patients and healthy controls, 21 significantly differentially methylated genes (12 genes with downregulated methylation levels and nine genes with upregulated methylation levels) between LMAQ patients and healthy controls, and 25 significantly differentially methylated genes (14 genes with downregulated methylation levels and 11 genes with upregulated methylation levels) between HMAQ patients and LMAQ patients.
for cocaine, increased goal-directed behavior, and increased motivation for the self-administration of cocaine (41). Another interesting finding is that the results of KEGG pathway analysis inferred the role of the morphine addiction pathway, instead of amphetamine addiction pathway, in the pathogenesis of MA addiction. Georgiou et al. (28) reported chronic methamphetamine use and abstinence can induce brain-region specific neuroadaptations of the m-opioid receptor (MOR) in mice. Regional manipulation of MOR expression in NAc and VTA of mice may be a novel approach to differential modulate MA behavioral sensitization (42). Guo et al. (43) reported acute intragastrical administration of naltrexone, a nonselective opioid receptor antagonist, could significantly reduce the cue-induced drug-seeking behavior of rats induced by intraperitoneal injection of MA after extinction training. Also, we checked the pathway annotation and including genes, and found morphine addiction pathway related to GABAergic synapse pathway and dopaminergic synapse pathway, but amphetamine addiction pathway only related to dopaminergic synapse pathway. Our results also reported other related pathways included cholinergic synapse pathway, glutamatergic synapse pathway, and GABAergic synapse pathway. This situation implies that MA addiction is more complicated than only dopaminergic mechanism.
There are few studies on genome-wide methylation of addictive substances, especially methamphetamine, and they are scattered in brain tissue samples of animal models and human blood samples, so the results are difficult to be inconsistent. Only Desplats et al. (18) reported genome-wide profiling of DNA methylation in postmortem frontal cortex tissue of with HIV-1 infection and MA dependence patients showed differential methylation on genes related to neurodegeneration; dopamine metabolism and transport; and oxidative phosphorylation. Fonteneau et al. (14) reported a genome-wide DNA methylation study with cocaine self-administering rats treated with DNA methyltransferase inhibitor, and the Ingenuity pathway analysis (IPA) of differentially methylated genes from the PFCx showed Akt/PI3K pathway and Rho-GTPase family was involved in the plasticity of axonal growth, synaptogenesis, and spine remodeling, underlying the effect of DNA methyltransferase inhibitor on the observed behavioral changes. Akt/PI3K pathway also highlighted in our results.
Meanwhile, there are a few DNA methylation researches of candidate genes in people or animals addicted to MA. Yuka et al. (44) found that the rate of DNA methylation at six CpG islands of SHATI/NAT8L promoter sites is significantly higher in MA users when compared to healthy subjects. Jayanthi et al. (45) reported that a single prior injection of MA enhanced MA self-administration (SA) and blocked SAinduced changes in DNA methylation and mRNA expression of potassium channels in the rat nucleus accumbens. MA produces a variety of epigenetic effects in the brain, which are seminal to establish long-lasting alterations in neuronal activity. Exposure to high and/or prolonged doses of MA induced persistently demethylation within alpha-synuclein gene (SNCA) promoter, which matches the increase in alpha-synuclein protein.
Demethylation was remarkable (10-fold of controls) and steady, even at prolonged time intervals being tested so far (up to 21 days of MA withdrawal) (46). MA-induced changes in long interspersed element-1 (LINE-1) partial methylation patterns are associated with MA-induced paranoia and could explain in part the pathophysiology of this type of psychosis. It is argued that MA-induced neuro-oxidative pathways may have altered LINE-1 partial methylation patterns, which in turn may regulate neuro-oxidative and immune pathways, which may increase risk to develop MA-induced paranoia (47). Veerasakul et al. (48) found that the methylation levels in the parvalbumin (PVALB) promoter were increased in MA-dependent patients using pyrosequencing, which might contribute to the GABAergic deficits associated with MA dependence. Hao et al. (49) reported that the methylation levels of the CHN2 gene, which encodes chimeric protein-2, which regulates axonal pruning, were increased in MA addiction patients. Nohesara et al. (29) demonstrated that DNA hypomethylation of the promoter regions of DRD3, DRD4, MB-COMT, and AKT1 was associated with increased expression of the corresponding genes in MAdependent patients with psychosis and, to a lesser extent, in non-psychotic MA-dependent patients. Xu et al. (50) measured the methylation levels of five CpGs (CpG1-5) on the BDNF promoter using pyrosequencing, and only the CpG5 methylation level was found to be significantly lower in MA addicts, which played a key role in the regulation of BDNF gene expression. Moreover, significant associations were identified between CpG5 methylation and addiction phenotypes, including tension-anxiety, anger-hostility, fatigue-inertia, and depressiondejection.

Screening and the Validation of DMP Data Revealed the Role of CAV2 in Different Qualities Mechanisms of MA Addiction
After further literature review and gene functional analysis, SLC1A6, BHLHB9, LYNX1, CAV2, and PCSK9 were selected from the DMP results for MethyLight qPCR verification with a large sample. The results indicated that the number of methylated copies of CAV2 in the LMAQ group was significantly higher than that in the HMAQ group, while the number of methylated copies of CAV2 was significantly lower than that in the HMAQ group and the healthy control group. The number of methylated copies of BHLHB9, LYNX1, PCSK9, and SLC1A6 was not significantly different between the HMAQ group and the LMAQ group. However, the number of methylated copies of BHLHB9 and LYNX1 was significantly higher in both the HMAQ group and the LMAQ group than in the healthy control group, while the number of methylated copies of PCSK9 and SLC1A6 was lower in both the HMAQ group and the LMAQ group than in the healthy control group.
In the LMAQ group, the number of methylated copies of the CAV2 gene was significantly correlated with age and the total duration of drug abuse. The number of methylated copies of the BHLHB9 gene was significantly correlated with the age of first drug abuse and the total duration of drug abuse in the LMAQ group. No other correlations were identified in the HMAQ group or the healthy control group. These five genes with abnormal methylation levels in MA abusers were all associated with neurological function and neurocognition.
It has been reported that only SLC1A6 and PCSK9 were associated with MA and/or ethanol abuse (51)(52)(53)(54). SLC1A6 is FIGURE 7 | Comparison of methylation levels of SLC1A6, BHLHB9, LYNX1, CAV2, and PCSK9 in three groups. Covariance analysis with age as covariance and Bonferroni post-hoc test as pairwise comparisons method showed that the significant differences in the methylated copies of CAV2, BHLHB9, LYNX1, PCSK9, and (Continued) FIGURE 7 | SLC1A6 in the three groups. The methylation copies of CAV2 were significantly higher in LMAQ than those in HMAQ groups, while the latter were significantly lower than those in healthy control groups. Methylation copies of BHLHB9, LYNX1, PCSK9, and SLC1A6 were not significantly different between the HMAQ and LMAQ groups. However, the methylation copies of BHLHB9 and LYNX1 were significantly higher in both HAMQ and LAMQ than those in healthy control groups, while the methylation copies of PCSK9 and SLC1A6 were significantly lower in both HAMQ and LAMQ than those in the healthy control groups.
also referred to as excitatory amino acid transporter 4 (EAAT4), a high-affinity glutamate transporter. The maintenance of low levels of extracellular glutamate and the avoidance of dosedependent neurotoxicity can be achieved only by glutamate transporter-mediated cellular uptake, such as SLC1A6, because there is no extracellular glutamate metabolic system in humans (55). Szumlinski et al. (51) identified the hyperglutamatergic state within the nucleus accumbens as a biochemical trait that corresponds to both genetic and idiopathic vulnerability for high MA preferences and consumption. Alshehri et al. (52) found that sequential exposure to ethanol and MA produced an additive effect in the downregulation of glutamate transporter 1 (GLT-1, also EAAT2) in the striatum and hippocampus rather than the glutamate aspartate transporter (GLAST, EAAT1) and cysteine/glutamate antiporter (xCT). Pickering et al. (53) found significantly decreased correct alterations and increased time in the choice area at the Y-maze center, as well as decreased mRNA expression of SLC1A6 in the medial prefrontal cortex of male Wistar rats after 2.5 g/kg ethanol treatment and 5 days of intermittent treatment (53). BHLHB9, also referred to as p60TRP (transcription regulator protein), is a member of the G-protein-coupled receptor (GPCR)-associated sorting protein (GPRASP) family. BHLHB9 can regulate the survival and differentiation of NGF-dependent neurons, and its expression level is decreased in the brains of Alzheimer's disease (AD) patients (56). BHLHB9, similar to GPRASP1, downregulates the expression of delta opioid receptor 1 (Oprd1) in neural stem cells (57). Mishra et al. demonstrated that BHLHB9 improved memory and learning abilities in BHLHB9 transgenic mice, which was attributed to increased synaptic connections and plasticity (57). Zhao et al. (58) reported that GPCRs interacted with β-site APP cleaving enzyme 1 (BACE1), a key secretase in AD pathogenesis. Moreover, the allosteric modulators and biased ligands of GPCRs displayed potential for the pharmaceutical treatment of AD patients, suggesting the possibility that GPCRs may be a therapeutic target for AD.
LYNX1 is a member of the Ly6/neurotoxin family associated with nicotinic acetylcholine receptors (nAChRs), which affect a wide array of biological processes, including learning and memory, attention, and addiction (59). Lynx1 modulates nAChR function in vitro by altering agonist sensitivity and desensitization kinetics. Dopaminergic cells of the substantia nigra pars compacta (SNc) and nearly all of the parvalbumin interneurons (GABAergic neurons) express lynx1 mRNA transcripts, whose dysfunctions are characteristics of psychiatric disorders (60). The titration of the lynx1 dosage can maintain neuronal health, while the lynx1 gene may participate in the trade-off between neuroprotection and augmented learning (61). Lynx1 prevents Aβ1-42-induced cytotoxicity in cortical neurons by competing with binding to nAChR subunits, and cortical Lynx1 levels are decreased in a transgenic mouse model with concomitant β-amyloid and tau pathology, which might have functional and pathophysiological implications in Alzheimer's disease (62).
PCSK9 is a member of the Bacillus subtilis protease (proprotein invertase) family, which is involved in the systemic control of blood cholesterol levels, neuronal differentiation, apoptosis, migration and nervous system development (63). Monoclonal antibodies that target PCSK9, a novel lipid-lowering approach, have recently been found to be accompanied by an increase in neurocognitive events (64). Picard et al. (65) found that PCSK9 expression levels in autopsy-confirmed AD brains were elevated in the frontal cortices compared to levels found in controls, both at the mRNA and protein levels. PCSK9 regulation seems to be under tight genetic control in females only, with specific variants that could predispose them to increased risk of AD. Hyperlipidemia regulates neuronal apoptosis in the hippocampal CA3 of mice by increasing PCSK9 and BACE1 expression, which may elucidate the role of lipid metabolism disorder in AD pathogenesis (66). A cross-tissue and cross-phenotypic analysis of genome-wide methylomic variation performed in alcohol use disorder (AUD) using samples from three discovery, four replication, and two translational cohorts indicated that replication in the AUD datasets confirmed PCSK9 hypomethylation. A translational mouse model of AUD showed that alcohol exposure caused a downregulation of PCSK9 (54).
CAV2, a membrane microdomain or "lipid raft, " has emerged as an essential functional module of the cell that is critical for the regulation of nerve growth factor (NGF) signaling and subsequent cell differentiation (67). BDNF increased the level of CAV2 in the rafts of hippocampal neurons for synapse development (68). Single-cell gene expression profiles showed downregulation of CAV2, GNB4, and lipase A in AD Rac1b-positive/p75 (NTR)-labeled basal forebrain (CBF) nucleus basalis neurons (69). CAV2 is also involved in the binding of the dopamine D1 receptor and the positive regulation of the dopamine receptor signaling pathway. Only two studies have reported the relationship between CAV1 and addiction. CAV-1 plays an important role in morphine-induced changes in the structural plasticity of primary cultured mouse cortical neurons (70). The compulsive-like intake behavior of MA was associated with an enhanced expression of D1Rs, increased activity of ERK1/2, and a decrease in Cav-1 expression in the dorsal striatum (71).
There were several limitations in this study. First, the subjects who abused MA in this study were recruited from the rehabilitation facilities in one province of China. They had not used MA for at least 3 months or even 2 years; hence, they were not acute users. Thus, there was a broad range of past frequencies and dosages of MA abuse. Second, in this study, we reported only the genome-wide DNA methylation status of MA abusers with different qualities to addiction and did not further explore its underlying mechanisms. Third, this is a DNA methylation research of peripheral blood cells of MA abusers. There are affirmative differences in the environment and differentiation between peripheral blood cells and neurons, but it is the only convenient method to study the DNA methylation alteration of alive MA abusers. Also, Uno et al. (72) reported that the ratios of Shati/Nat8l CpG island methylation were significantly decreased in both the nucleus accumbens and the peripheral blood of MA-induced murine models of schizophrenia-like phenotype compared with those of control mice. To our knowledge, this is the first study to show similar alterations of DNA methylation status in both the brain and peripheral blood of MA treated mice, which suggests the feasibility of this kind of research. For the next step, gene knockout or transgenic animals should be used to investigate the roles of circadian entrainment pathways and the caveolin-2 gene in the pathogenesis of MA addiction.

CONCLUSIONS
This study investigated the genome-wide DNA methylation status of different addiction qualities in MA abusers and identified seven specific pathways with abnormal methylation status, which were all reported to be associated with MA addiction or reward circuits, especially the circadian entrainment pathway, a relatively new pathogenesis. Moreover, CAV2 was the only gene among the five candidate genes screened from the genome methylation chip data that showed significant differences in methylation levels in all pairwise comparisons of the three groups. Further studies on their functions and mechanisms will help us to better understand the pathogenesis of MA addiction and to explore new targets for drug intervention.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in GEO database (GSE154971).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of the Second Xiangya Hospital. The participants provided their written informed consent to participate in this study. Written informed consent was obtained from each individual for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
WH was responsible for the study design and supervision. LL, TLu, HD, CZ, and TLi had participated in the acquisition, analysis, and interpretation of the data. LL wrote a draft of the manuscript. XZ critically revised the manuscript for important intellectual content and provided the technical support. All authors contributed to the article and approved the submitted version.