Expression Quantitative Trait Locus rs6356 Is Associated with Susceptibility to Heroin Addiction by Potentially Influencing TH Gene Expression in the Hippocampus and Nucleus Accumbens

Opioid addiction is a complicated and highly heritable brain disease. Dysfunction in dopaminergic signaling is involved in the pathogenesis of addictive disorders. Encoding a dopamine synthetase, the tyrosine hydroxylase (TH) gene has long been an interesting candidate in genetic association studies for opioid addiction. However, the mechanisms underlying associations of risk gene variants and opioid addiction remain unknown. In the present study, we first analyzed the association between TH gene variants and susceptibility and traits of heroin addiction in 801 patients with heroin addiction and 930 healthy controls. Methylation levels in the promoter region of the TH gene were detected and compared between the heroin addiction and healthy control groups. To reveal the potential mechanism of the association of TH gene variants and heroin addiction, correlations between the risk TH single nucleotide polymorphism (SNPs) for heroin addiction and the methylation and expression levels of the TH gene were examined. Our results demonstrated that SNP rs6356 was associated with susceptibility to heroin addiction. CpG TH_15 was hypermethylated in the heroin addiction group compared with the healthy control group. Notably, SNP rs6356 was correlated in an allele-specific manner with expression of the TH gene in the hippocampus and nucleus accumbens but not with methylation levels of CpG TH_15. Our findings suggest that the eQTL rs6356 was associated with susceptibility to heroin addiction by potentially affecting the expression of the TH gene in brain regions in the mesocorticolimbic dopamine system, including the hippocampus and nucleus accumbens.


Introduction
Substance addiction is a complex brain disease characterized by compulsive substance seeking, abstinence, and repeated relapse to substance use (Chisholm et al. 2021;Huang et al. 2021). Opioids, one of the most abused prescription medicines for pain management and illicit addictive drugs, cause the most harm in terms of substance misuse-related deaths, killing tens of thousands of people per year in America alone (Browne et al. 2020). In China, more than 730 thousand people were exposed to illicit opioids, mainly heroin, in 2020. Mounting evidence has suggested that genetic underpinnings and environmental factors contribute to the development of substance addiction (Gerra et al. 2021). Twin studies have revealed that the heritability of opioid addiction is approximately 70% (Goldman et al. 2005). Genome-wide association studies have proposed a large number of risk gene variants for opioid addiction (Song et al. 2020;Zhou et al. 2020). However, the mechanisms underlying of associations of risk gene variants and opioid addiction remain unknown.
Prolonged use of addictive substances could lead to adaptive alterations in neural plasticity in discrete brain regions, including the mesocorticolimbic dopamine system, which processes reward-and motivation-related effects (Hyman et al. 2006), resulting in the transition to addicted states (Stewart et al. 2021). DNA methylation, a type of epigenetic modification, represents a mechanism that affects the plasticity process and then the development of substance addiction that involves environmental stimuli and the regulation of gene expression patterns (Stewart et al. 2021). DNA methylation, especially in the promoter region of the gene, which is rich in CpG sites and transcription factor binding sites, is an important mechanism that regulates gene expression (Pathak et al. 2018). Genetic variants can be associated with the methylation and expression levels of a gene. These variants are known as methylation quantitative trait loci (mQTLs) (Villicaña and Bell 2021) or expression quantitative trait loci (eQTLs) (Degtyareva et al. 2021). Previous studies have revealed that mQTLs and eQTLs are enriched in risk genetic variants in many neuropsychiatric diseases, such as schizophrenia (Perzel Mandell et al. 2021;Zhao et al. 2018), depression (Ciuculete et al. 2020;Wang et al. 2020a), and substance addiction (Lin et al. 2020;Zhang et al. 2014). Thus, risk genetic variants that correlate with the methylation or expression of a gene may be an important mechanism underlying the association with susceptibility and traits of substance addiction.
Tyrosine hydroxylase (TH) is a key rate-limiting enzyme in dopamine and norepinephrine biosynthesis (Boundy et al. 1998;Walters et al. 2003). By regulating the levels of these neurotransmitters that have essential roles in neuronal activity, TH has been shown to be involved in many neuropsychiatric disorders, including opioid addiction (Jalali Mashayekhi et al. 2018;Vaillancourt et al. 2021). In rodents, repeated morphine administration increased TH expression in the hippocampus (Fang et al. 2017), an important brain region in the mesocorticolimbic dopamine system. An aldehyde dehydrogenase 2 (ALDH2) inhibitor could prevent the reinstatement of cocaine and alcohol intake by suppressing activated tyrosine hydroxylase signaling (Diamond and Yao 2015). Hence, the TH gene has long been an interesting candidate for genetic association studies of opioid addiction. However, few studies have elucidated the underlying mechanisms.
In the present study, we first determined the association between TH gene variants and susceptibility to heroin addiction. Then, the association of TH gene variants and traits of heroin addiction were analyzed. To reveal the potential mechanisms underlying associations between risk single nucleotide polymorphism (SNPs) in the TH gene and heroin addiction, we detected and compared the methylation level in the promoter region of the TH gene. Association analysis between the risk SNPs and methylation levels of CpG sites that showed differential methylation between patients with heroin addiction and healthy controls was performed. Finally, the risk SNPs were screened in the genotype-tissue expression (GTEx) database to demonstrate whether they were eQTLs in human brain regions.

Study Subjects
A total of 1731 unrelated Northwestern Han Chinese individuals participated in the present study. In the heroin addiction group, 801 patients in a methadone maintenance treatment (MMT) program were recruited from Xi'an Mental Health Center. All patients were diagnosed by senior psychiatrists according to the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-V). Patients with other substance addictions, including tobacco and alcohol, or other psychiatric disorders were excluded. Demographic information, including age, gender, marital status, education, occupation, income, and family relation (Tables S1 and S2), as well as traits of heroin, addiction including the age of onset for heroin use, heroin dose per time, heroin intake times per day, duration from first heroin exposure to addiction and effective methadone dose (Table S3), were collected by a questionnaire survey of the patients. Regarding the healthy control group, 930 healthy subjects were enrolled from the Health Examination Centers of the First Affiliated Hospital of Xi'an Jiaotong University and Tangdu Hospital. They had no substance (including tobacco and alcohol) addictions, no personal or family history of psychiatric disorders, and no severe organic diseases. There were no age or gender distribution differences between groups (Table S1).
All study subjects were well informed and provided written informed consent for participating in the study. The protocol of our study was approved by the Medical Ethics Committees of Xi'an Jiaotong University and the Xi'an Mental Health Center, in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki).

SNP Selection and DNA Extraction
For the analysis of associations between TH gene variants and the risk of heroin addiction, tag SNPs in the functional regions (promoter region, untranslated regions, exons, and intron-exon boundaries) of the TH gene with a minor allele frequency (MAF) > 0.05 and a minor genotype frequency > 0.01 in the Han Chinese population of the 1000 Genomes database were selected. For the association analysis of TH gene variants and traits of heroin addiction, all functional SNPs with frequencies meeting the above criteria were genotyped in the heroin addiction group. In addition, all SNPs met Hardy-Weinberg equilibrium (HWE) in the healthy control group.
Genomic DNA was extracted from the peripheral blood with the EZNA™ Blood DNA Midi Kit (Omega Bio-TeK, Norcross, GA, USA) in accordance with the manufacturer's instructions. Following quality control, the DNA was subjected to SNP genotyping and detection of DNA methylation.

SNP Genotyping and Detection of DNA Methylation
SNP genotyping was performed with the standard SNaPshot assay by Genesky Biotechnologies (Shanghai, China). The corresponding primers were first designed according to the sequences around the SNPs. Then, a fragment of DNA surrounding the SNPs was amplified by polymerase chain reaction (PCR). The final product was separated by capillary electrophoresis with the ABI 3130XL Genetic Analyzer (Applied Biosystems Co. Ltd., CA, USA). The GeneMapper 4.1 software was used to obtain the genotypes of the SNPs. To ensure the accuracy of genotyping, five percent of the samples were randomly selected and replicated.
For methylation detection, DNA was first bisulfite converted with the EZ DNA Methylation Gold™ Kit (Zymo Research Corporation, CA, USA). Then, the target region of the TH gene promoter, containing the CpG sites, was amplified with the corresponding primers (Table S4). Finally, the amplified products of the target fragments were analyzed using the Illumina HiSeq X Ten platform. The methylation levels of the CpG sites were defined as the percentage of methylated sites in all sequences containing this site. DNA methylation detection was also conducted by Genesky Biotechnologies.

Statistical Analysis
PLINK (version 1.09) was used for the association analysis of TH gene variants and risk of heroin addiction, with age and gender as covariates. The associations between TH gene variants and traits of heroin addiction were assessed using Matrix eQTL (website: http:// www. bios. unc. edu/ resea rch/ genom ic_ softw are/ Matrix_ eQTL) with the demographic information (age, gender, marital status, education, occupation, income per month, and family relation) as covariates. For the comparison of DNA methylation between groups, Student's t tests were used. Linear regression analysis with SPSS 22.0 was conducted for the association of genotypes of TH SNPs and levels of DNA methylation. The significance threshold was set at an adjusted P value < 0.05 after Bonferroni correction.

Associations Between TH Gene Variants and Susceptibility to Heroin Addiction
To determine associations between TH gene variants and the risk of heroin addiction, two tag SNPs (rs10770140 and rs6356) were genotyped in 930 healthy controls and 801 patients with heroin addiction. Our results revealed that SNP rs6356 was associated with susceptibility to heroin addiction. The distribution of genotypes of SNP rs6356 was different between groups (χ 2 = 8.636, P = 0.013). In the dominant inheritance model, the GG and GA genotype frequencies were lower in the heroin addiction group (26.8%) than in the healthy control group (33.0%) (χ 2 = 7.763, P = 0.005). Regarding the minor allele frequency, the frequency of the G allele was lower in the patients with heroin addiction (14.8%) than in the healthy controls (17.7%) (χ 2 = 5.660, P = 0.017). All the above significance tests survived Bonferroni correction (α = 0.05/2 = 0.025) ( Fig. 1 and Table 1).

Associations Between TH Gene Variants and Traits of Heroin Addiction
To reveal the associations between the genotypes of TH SNPs and traits of heroin addiction, five functional SNPs were genotyped in the heroin addiction group, and linear regression analysis was conducted between the SNP genotypes and five traits of heroin addiction. According to our findings, only the rs10770140 polymorphism was nominally associated with the age of onset for heroin use ( Fig. 2; beta = −1.347, P = 0.028). Unfortunately, this association was not significant after Bonferroni correction (α = 0.05/5 = 0.010). The detailed results of the associations between TH SNPs and traits of heroin addiction are shown in Table S5.

Differential Methylation in the Promoter Region of the TH Gene Between the Healthy Control and Heroin Addiction Groups
DNA methylation plays vital roles in substance addiction (Browne et al. 2020). To explore the potential mechanisms underlying the associations between TH gene variants and heroin addiction, we first detected and compared methylation levels in the promoter region of the TH gene between groups. Student's t test revealed that the methylation level of CpG TH_15 was significantly higher in the heroin addiction group than in the healthy control group (Fig. 3; t = −3.465, P < 0.001), even after Bonferroni correction (Table 2; α = 0.05/15 = 0.003). However, methylation levels of other CpG sites and the mean methylation level of the promoter region of the TH gene did not show differences between groups (all P > 0.05).

Association of SNP rs6356 and Methylation Level of CpG TH_15
Previous studies suggested that mQTLs, SNPs that have an association with the level of gene methylation, convey more vulnerability in human diseases, including substance addiction than randomly selected SNPs (Lin et al. 2020;Min et al. 2021). Based on the differential methylation results above, we then analyzed the association between rs6356 genotypes and methylation levels of CpG TH_15. Linear regression analysis demonstrated that there was no association between them (P > 0.05).

Expression Quantitative Trait Locus of rs6356 in Human Brain Regions
To further reveal the potential function of rs6356 in heroin addiction, we screened for it in the GTEx database (website: www. gtexp ortal. org). SNP rs6356 was found to be associated in an allele-specific manner with the expression level of the TH gene in human brain regions, including hippocampus (P = 3.83E-04), caudate (P = 6.50E-03), nucleus accumbens (NAc, P = 8.22E-03), hypothalamus (P = 1.59E-02), and putamen (P = 1.91E-02). The subjects with the AA genotype had the highest expression levels of TH gene in the above  Fig. 3 Differential methylation in CpG TH_15 between patients with heroin addiction and healthy controls brain regions, followed by the GA and GG genotype subjects. The detailed outcome of the eQTL results of rs6356 is shown in Table 3.

Discussion
TH is involved in the rewarding effects of addictive substances by indirectly regulating dopaminergic signaling (Fang et al. 2017;Vaillancourt et al. 2021). The TH gene has long been an important candidate in genetic association studies in substance addiction. Further investigations of the potential mechanisms underlying the correlations between TH gene variants and susceptibility to and traits of opioid addiction are needed. In the present study, we first determined the association between TH gene variants and susceptibility to and traits of heroin addiction. Methylation levels in the promoter region of the TH gene were detected and compared between the heroin addiction and healthy control groups. Finally, mQTL and eQTL analyses were performed for the risk SNPs in heroin addiction to explore their potential mechanism or function in opioid addiction. Long-term exposure to addictive substances induces neuronal adaptations, including structural, functional, and molecular changes in the brain (Evangelou et al. 2021;Fernàndez-Castillo et al. 2021;Salery et al. 2020;Snyder and Silberman 2021). Epigenetic and gene expression alterations are essential molecular changes responsible for the neuronal adaptations that contribute to the development of substance addiction (Fernàndez-Castillo et al. 2021). Previous studies have indicated that hypermethylation at several CpG sites in the TH gene was found in postmortem NAc samples in patients with cocaine addiction compared with healthy controls (Vaillancourt et al. 2021). In accordance with the above findings, we observed that methylation levels of CpG TH_15 were higher in the heroin addiction group than in the healthy control group. DNA methylation is one of the mechanisms that could inhibit gene transcription and expression (Li et al. 2020). With a potential role in regulating DNA methylation, mQTLs are believed to constitute more liability than random SNPs in diseases including heroin addiction . Unfortunately, we failed to detect any association between the risk SNP rs6356 for heroin addiction and the methylation level of TH_15. Hence, mechanisms other than methylation patterns may underlie the association of rs6356 and the risk of heroin addiction.
Our results demonstrated that SNP rs6356 is an eQTL that is allele-specific correlated with expression levels of the TH gene in the hippocampus and NAc. The A allele of rs6356 was associated with higher expression of the TH gene in the  above two brain regions. In addition, the A allele frequency in the patients with heroin addiction was found to be significantly higher than that in the healthy controls. This was consistent with previous findings in rodents that repeated morphine exposure caused upregulation of TH gene expression in the hippocampus (Fang et al. 2017) and locus coeruleus (Mashayekhi et al. 2018), two important brain regions in substance addiction. The hippocampus is involved in the formation of pathological contextual memories in substance addiction (Wang et al. 2020b). The NAc is responsible for the addictive substance-induced rewarding and reinforcing effects (Al-Hasani et al. 2021;Amaral et al. 2020). Increased dopamine release in brain regions in the mesocorticolimbic dopamine system, including the hippocampus and NAc, has been shown in all kinds of addictive substance use (Morales and Pickel 2012;Pinheiro et al. 2015). Increases in dopamine synthesis may be associated with the upregulation of TH expression in these brain regions. Thus, we speculate that the eQTL rs6356 is associated with heroin addiction by influencing the TH gene expression in brain regions, including the hippocampus and NAc. Several limitations should be considered when interpreting the results of our study. First, we used a candidate gene method to elucidate the role and potential function of TH gene variants in heroin addiction. Genome-wide association studies with large numbers of participants have the highest resolution in identifying risk gene variants for addictive diseases. Second, the eQTL results were obtained from the GTEx database, and we did not examine the expression level of the TH gene in the relevant brain regions due to a lack of postmortem brain samples. Third, functional studies should be conducted to verify the precise role of TH gene variants in the development of heroin addiction.
In summary, we found that the eQTL rs6356 was associated with susceptibility to heroin addiction. The CpG site TH_15 in the promoter region of the TH gene was hypermethylated in the heroin addiction group compared with the healthy control group. SNP rs6356 was associated in an allele-specific manner with expression levels of the TH gene in brain regions in the mesocorticolimbic dopamine system, including the hippocampus and NAc, but not the methylation level of CpG TH_15 in peripheral blood. Our findings suggest that the eQTL rs6356 is associated with opioid addiction, potentially by influencing TH gene expression in the hippocampus and NAc. These findings further extend the mechanisms underlying associations of genetic variants and addictive disorders.
Author Contribution Jianbo Zhang and Yongsheng Zhu conceptualized this study. Kena Wang and Hongbo Zhang analyzed the data and drafted the manuscript. Jinshan Ji and Rui Zhang contributed to the methodology. Wei Dang and Qiaoli Xie reviewed and edited the manuscript. All the authors reviewed and approved the final version of publication.
Funding This work was supported by the Natural Science Foundation of Shaanxi Province (No. 2021JQ-065, 2021JM-028, and 2020JQ-802) and the Young Talent Promotion Program of Shaanxi Science and Technology Association (No. 20200301).

Availability of Data and Materials
All data and materials will be available upon reasonable request. The sequencing data of the methylation assaying has been deposited into Sequence Read Archive (SRA) under accession number PRJNA752473.

Declarations
Ethics Approval and Consent to Participate Our study was approved by the Medical Ethics Committees of Xi'an Jiaotong University and the Xi'an Mental Health Center, in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki). All study subjects were well informed and provided written informed consent for participating in the study.

Consent for Publication
The authors affirm that all participants provided informed consent for publication.

Competing Interests
The authors declare no competing interests.