Combined study of genetic and epigenetic biomarker risperidone treatment efficacy in Chinese Han schizophrenia patients.

Nowadays, risperidone is an atypical antipsychotic drug that has been increasingly used for treatment and maintenance therapy in schizophrenia. However, partially affected by genetic or environmental factors, there is significant difference in treatment outcomes among patients. In this study, we aimed to interpret the difference between good and poor responders treated with risperidone in both genetic and epigenetic levels in 288 mainland Chinese patients. We recruited a Henan cohort including 98 patients as initial discovery group and then confirmed our results in Shanghai cohort. In genetic studies, we found 10 candidate single-nucleotide polymorphisms (SNPs) and 2 rare variants in Henan cohort by next-generation sequencing of 100 risperidone-response-related genes. After replication in Shanghai cohort by massarray platform, ultimately, rs6706232 and rs4818 were significantly associated with risperidone response in the two cohort meta-analysis (P=0.024 and 0.04, respectively). Besides, we also selected another reported 17 candidate SNPs associated with risperidone drug response to replicate in our mainland Chinese samples, while, we found no significant SNPs after Bonferroni correction. In epigenetic studies, we investigated the methylation status in promoters or gene-coding region of risperidone drug response-related genes including CYP3A4, CYP2D6, ABCB1, HTR2A, DRD2. Totally we found seven significant CpG sites in the meta-analysis with Bonferroni-corrected PCYP3A4_CpG_-36=0.0014, PCYP3A4_CpG_-258=0.0013, PCYP3A4_CpG_-296=0.0014, PCYP3A4_CpG_-367:-372:-374=0.028, PCYP2D6_CpG_193=0.012, PCYP2D6_CpG_242:244:250=0.00076 and PCYP2D6_CpG_284=0.034, respectively. As genetic and epigenetic factors may interactively affect drug response, we finally carried out a multivariant interaction analysis with multifactor dimensionality reduction and discovered a significant four-locus model (CYP3A4_CpG_-82:-86 +rs6280+rs1800497+rs6265, P=0.038) affecting drug response. These findings could partially explain different risperidone response outcome in Chinese population in a systematic level.

Introduction

Risperidone is an antipsychotic drug that has been increasingly used for treatment and maintenance therapy in schizophrenia and related psychotic disorders.[1] Risperidone is metabolized to its active metabolite, 9-hydroxyrisperidone (9-OH-risperidone with potent antagonistic properties for the dopamine D2 and serotonin 5-hydroxytryptamine-2 (5-HT2) receptors). However, there are significant interindividual differences in clinical response and side effects, meanwhile, optimizing drug treatment for patient is often by trial and error which costs a lot of time and money, so it is crucial to identify more novel drug-response-related markers to predict drug response.

The variability in the risperidone response can be caused by genetic, epigenetic, physiologic and environmental factors. Genetic factors are mostly assumed to have a close relationship with drug treatment response,[2] and on this basis, a number of pharmacokinetic studies have been performed. To date, most studies have typically focused on candidate genes, mainly selected drug metabolizing enzyme genes, transport genes, neurotransmitter receptors genes, such as dopamine or serotonin receptors. Several studies have shown positive associations between genetic variation and risperidone response, for example, CYP2D6, ABCB1, dopamine receptors and serotonin receptors have ever been reported to be significant associated with risperidone’s efficiency and risperidone-induced adverse effects.[3, 4, 5] However, these studies mainly focused on a few genes through candidate gene association study method and most were performed in small sample size without independent replication, because of which these results could not been confirmed by different groups and used in clinical practices. In recent years, 14 candidate genes have been identified in relation to risperidone treatment response in a genome-wide association study.[6] Although this study was also performed in relatively small sample size, it has the low power to uncover the typical genome-wide association study significant variants and sometimes the genome-wide association study results are not replicated across other studies or populations.

In addition, epigenetic factors can also affect drug treatment by modulating the expression of key genes involved in the metabolism and distribution of drugs as well as drug targets, thereby contributing to interindividual variation in drug response,[7] which mainly derived from DNA methylation modification changes. There have been some reports that DNA methylation status may serve as a pharmacogenomics biomarker. More recently, DNA methylation in drug absorption, distribution, metabolism, excretion (ADME) genes such as GSTP1, GPx3, ABCB1, ABCG2 and the nuclear receptor vitamin D receptor have been reported to be associated with drug response,[8, 9, 10, 11] which reflects a strong potential of epigenetic marks to serve as predictors of antipsychotic drug response. However, epigenetic studies of risperidone response have rarely been reported. Furthermore, the drug response of risperidone involves complex drug ADME-related molecular networks and pathways; as a result, the common variants associated with risperidone may have smaller effect sizes, and they may predict a response only when combined variants in genes in a known molecular pathway test whether the pathway is associated with the phenotype. This pathway-based association approach provides a more powerful strategy for pharmacogenomics study.

In this study, to comprehensively discover the predictor of risperidone response, we conducted the pharmacogenomics study using target sequencing technology and epigenetic study using massARRAY technology. Also, we performed combined analysis of different markers and replication studies in independent subjects, which laid the foundation for comprehensively discovering the predictor of risperidone response in Chinese population.

Materials and methods

Subjects

Two cohorts of Chinese Han in-patients with schizophrenia were enrolled in the present study. Figure 1 is research workflow. One hundred and ninety patients were recruited from Shanghai Mental Health Center, 98 patients were recruited from Henan Provincial Mental Health Center. All patients recruited met the following criteria: (1) they satisfied the Diagnostic and Statistical Manual of Mental Disorders-IV (DSM-IV) criteria for schizophrenia; (2) they had no physical complication or other substance abuse; (3) they had no history suggesting resistance to antipsychotic drug treatment; (4) they had not received any medication for 4 weeks; and (5) they had not previously received second-generation antipsychotics. The study protocol was drawn up according to the principles of the Helsinki Accord and was reviewed and approved by the Shanghai Ethical Committee of Human Genetic Resources. The statement of informed consent was obtained from all the subjects after full explanation of the procedure. Genomic DNA was extracted by QIAamp DNA Blood Midi Kit (Qiagen, Hilden, Germany) and quantified on three platforms including NanoDrop 2000 (Thermo, Wilmington, DE, USA), Qubit Fluorimeter (Invitrogen, Eugene, OR, USA), 2100 Bioanalyzer (Agilent, Waldbronn, Germany) according to manufacturer’s protocol.

Figure 1

Research workflow of this study. ADME, absorption, distribution, metabolism, excretion; MDR, multifactor dimensionality reduction; NGS, next-generation sequencing; SNP, single-nucleotide polymorphism.

Clinical assessment

For the recruited subjects, the dosage of risperidone was 2 mg per day initially and then gradually increased to 4 mg per day within the first week, which was maintained until the end of week 2. After that, the dosage was adjusted according to individual tolerance. For all the participants, medication compliance was closely monitored and confirmed by nursing staff, and no other medication was given except bedridden for extrapyramidal side effects, flunitrazepam for insomnia and sennoside for constipation during the study period.

Clinical effect was assessed on the Positive and Negative Syndrome Scale (PANSS), including the positive, negative and general psychopathology subscales. For the recruited patients, clinical assessments were conducted on the day of admission, as well as at the end of week 4. In each cohort, all PANSS ratings were conducted independently by two qualified psychiatrists, who were blind to the genotype of patients. And the inter-rater reliability between the two psychiatrists is good. Risperidone treatment efficacy was measured in terms of the reduction in PANSS scores. For the independent samples, risperidone response was classed into four groups, cured (⩾75%), significant progress (50–74%), progress (25~49%) and ineffective (<25%) based on PANSS scores.

Methylation analysis

We screened CpG-rich spot in the upstream of promoter region or within candidate genes including CYP3A4, CYP2D6, HTR2A, ABCB1 and DRD2 by UCSC Human Genome Browser Gateway (http://genome.ucsc.edu/cgi-bin/hgGateway), and designed specific PCR primers for bisulfate treatment amplification by EpiDesigner software (www.epidesigner.com). We detected CpG methylation status on MassARRAY Analyzer 4 platform and analyzed using software Epityper1.0.5 (Sequenom, San Diego, CA, USA).

Next-generation sequencing

The systematic association study between exon polymorphisms of candidate genes and response to risperidone in Chinese Han schizophrenia patients was carried out with Miseq pair-end sequencing technology after HaloPlex Target Enrichment system (HaloPlex Custom Kits, 1–500 kb, ILMFST, 96, G9901B). The primers were designed on SureDesign software (http://earray.chem.agilent.com/suredesign/, Agilent) and data analysis was performed on SureCall (www.genomics.agilent.com) software. The result of next-generation sequencing (NGS) was validated by iPLEX Gold SNP method on MassARRAY Analyzer 4 platform.

Statistical analysis

Statistical analysis was carried out as described previously[12] with minor modifications. Cohort characteristics including gender, age and weight differences of demographic and clinical variables were examined first to affirm the homogeneity of the good responders and poor responders in our samples with Student’s t-test. The reduction of the total and subscale scores of PANSS was used as a measure of clinical improvement of risperidone treatment. SPSS for Windows, version 11.0 (SPSS, Chicago, IL, USA) was used for statistical analysis.

To substantiate the results, allele and genotype frequencies of each polymorphism were compared between good-responders and poor-responders groups using the χ2 test on the online software SHEsis (http://analysis.bio-x.cn).[13] And clinical good responders were defined as patients with 50% or even higher reduction in PANSS scores than the average level of all subjects; correspondingly, poor responsers were defined as patients with lower reduction than 50% in PANSS scores. All tests were two-tailed and statistical significance was assumed at P<0.05.

The mutual interactions between methylation sites and single-nucleotide polymorphisms (SNPs) were analyzed on multifactor dimensionality reduction (MDR) software as previously described.[14] In the configuration file, 10-fold cross-validation was defined and the threshold ratio was set at 1.0. We ran the analysis 10 times using constant random number seeds and the results were averaged to avoid spurious outcomes due to chance divisions of the data. This MDR procedure can be carried out for each possible model size, if computationally feasible. Due to computation restrictions, we considered two-locus interactions through four-locus interactions. We determined statistical significance by comparing the average prediction error from the observed data with the distribution of average prediction errors under the null hypothesis of no association derived empirically from 1000 permutations. The null hypothesis was rejected when the upper-tail Monte Carlo P-value derived from the permutation test was <0.05.

Results

Patients characteristics

Detailed information about patients and research workflow could be found in Table 1 and Figure 1, respectively. There is no gender, age and weight difference between risperidone good-response group and poor-response group either in the Henan cohort or in the Shanghai cohort or in the whole sample. A power calculation indicated that we had the power of 80.9% to detect many of the SNPs with effect size =0.167 and d.f. =1 in the combined cohort of Henan and Shanghai.

Table 1

Descriptive statistics for patient-related variables with regard to good and poor responders

Cohorts and variables	Henan			Shanghai			Meta-analysis
	Good responders	Poor responders	P	Good responders	Poor responders	P	Good responders	Poor responders	P
Male	14	22	0.860	74	32	0.643	88	54	0.176
Female	23	39		56	28		79	67
Age	33.57±10.34	29.67±9.04	0.053	37.30±15.26	41.12±19.35	0.182	36.46±14.36	35.35±16.06	0.538
Weight (kg)	58.41±10.55	58.54±10.36	0.95	61.24±13.18	60.57±15.22	0.773	60.58±12.64	59.51±12.90	0.507

Age and weight difference between good and poor responders were tested by Student's t-test and gender difference was based on χ2 test. Value: mean±s.d.

The discovery and validation of genetic variant biomarkers associated with risperidone response

A hundred candidate genes that related with risperidone response were introduced in our study (Supplementary Table 1). In target sequencing, average read depth was 63.7 ×, and 1 ×, 8 ×, 20 × coverage with Q30 were 94.35%, 87.70%, 74.34%, respectively. A total of 330 SNPs were found by plink software with filtration conditions such as MAF (minor allele frequency) >0.01, Hardy–Weinberg equilibrium P>0.001, call rate >95% and then 10 SNPs were significantly associated with risperidone treatment response by target sequencing after data quality control and tested with generalized linear regression model with non-genetic confounding factors (age, onset age, weight) as covariants (Table 2). However, after P-value correction, all the 10 SNPs were no longer significant statistically (P>0.05). Second, to validate the technical accuracy of Miseq sequencing platform, we performed MassARRAY analysis in 80 Henan samples with specific primers (Supplementary Table 2) and the average accordance of the two platforms reached 94.88%, which guaranteed a high credibility of our data (data not shown). Besides, when we neglected the non-genetic confounding factors (age, onset age, weight), we found two mutations with low frequency (1%< MAF <5%), which were significantly different in risperidone treatment response (Supplementary Table 3). Third, in order to validate whether the 10 candidate SNPs detected by NGS remain significant, based on MassARRAY platform, we first tested the results in a smaller Henan cohort and then we found that five SNPs and four genotypes were significantly associated with risperidone treatment response without P-value correction; then we used a larger Shanghai cohort for further validation. A meta-analysis was carried out in the whole cohort, two SNPs and one genotype showed statistical significance before P-value correction ultimately (Table 3). For comparison, we included another 17 reported risperidone-response-related SNPs for validation in our cohort, while we only validated two significant SNPs in Henan cohort and two genotypes in meta-analysis without P-value correction (Table 4).

Table 2

10 SNPs were significantly related with risperidone treatment response by NGS

Gene	Chr	SNP	Pos/Mb	MA	PR_freq	R_freq	OR (95% CI)	P
UTG1A3	2	rs6706232	234.6	A	0.25	0.50	0.41 (0.20–0.85)	0.017
CYP1B1	2	rs1056827	38.3	A	0.21	0.1	4.37 (1.20–15.91)	0.026
CYP1B1	2	rs10012	38.3	C	0.18	0.1	3.16 (1.00–9.93)	0.049
DRD3	3	rs6280	113.9	C	0.42	0.24	2.60 (1.09–6.20)	0.031
HTT	4	rs362267	3.2	T	0.44	0.24	2.50 (1.10–5.70)	0.029
HTT	4	rs362306	3.2	A	0.49	0.31	2.22 (1.01–4.90)	0.048
CYP2E1	10	rs2515641	135.4	T	0.07	0.24	0.22 (0.07–0.66)	0.007
COMT	22	rs4633	20.0	T	0.28	0.13	3.56 (1.24–10.21)	0.018
COMT	22	rs4680	20.0	A	0.28	0.15	3.06 (1.13–8.28)	0.028
COMT	22	rs4818	20.0	G	0.31	0.48	0.42 (0.18–0.96)	0.039

Abbreviations: Chr, chromosome; CI, confidence interval; MA, minor allele; NGS, next-generation sequencing; OR, odds ratio; PR, poor responder; R, responder; SNP, single-nucleotide polymorphism. P and OR were tested using logistic regression model with non-genetic confounding factors as covariants.

Table 3

Meta-analysis of candidate SNPs associated with risperidone treatment response in two cohorts with two factors adjustment (age and gender)

Cohort	Gene	SNPs	Allele		OR (95% CI)	Chi-square	P1	P1a	Genotype			P2	P2a
Henan	UGT1A3	rs6706232	G	A					GG	AA	GA
		Good responders	28 (0.666)	14 (0.333)	2.681 (1.055–6.810)	3.578	0.038	0.260	9 (0.428)	2 (0.095)	10 (0.476)	0.119	0.610
		Poor responders	75 (0.815)	17 (0.184)					30 (0.652)	1 (0.021)	15 (0.326)
	HTT	rs362267	G	A					GG	AA	GA
		Good responders	46 (0.741)	16 (0.258)	0.469 (0.229–0.960)	5.584	0.018	0.264	18 (0.58)	3 (0.096)	10 (0.322)	0.035	0.293
		Poor responders	64 (0.561)	50 (0.438)					17 (0.298)	10 (0.175)	30 (0.526)
	CYP2E1	rs2515641	C	T					CT	CC
		Good responders	48 (0.799)	12 (0.199)	2.933 (1.027–8.379)	2.967	0.044	0.260	18 (0.6)	12 (0.399)		0.06	0.293
		Poor responders	102 (0.894)	12 (0.105)					45 (0.789)	12 (0.21)
	COMT	rs4818	G	C					GG	CG	CC
		Good responders	26 (0.565)	20 (0.434)	2.753 (1.133–6.692)	6.772	0.025	0.260	8 (0.347)	10 (0.434)	5 (0.217)	0.010	0.293
		Poor responders	24 (0.324)	50 (0.675)					2 (0.054)	20 (0.540)	15 (0.405)
	DRD3	rs6280	C	T					CT	TT	CC
		Good responders	11 (0.323)	23 (0.676)	0.121 (0.022–0.669)	3.145	0.015	0.260	11 (0.647)	6 (0.352)	0 (0)	0.046	0.293
		Poor responders	25 (0.520)	23 (0.479)					19 (0.791)	2 (0.083)	3 (0.125)
Shanghai	UGT1A3	rs6706232	G	A					GG	AA	GA
		Good responders	446 (0.753)	146 (0.246)	1.395 (0.921–2.115)	2.763	0.116	1	154 (0.52)	4 (0.013)	138 (0.466)	0.114	1
		Poor responders	211 (0.805)	51 (0.194)					80 (0.61)	0 (0)	51 (0.389)
	HTT	rs362267	G	A					GG	AA	GA
		Good responders	420 (0.648)	228 (0.351)	0.956 (0.711–1.287)	0.12	0.767	1	136 (0.419)	40 (0.123)	148 (0.456)	0.91	1
		Poor responders	182 (0.636)	104 (0.363)					57 (0.398)	18 (0.125)	68 (0.475)
	CYP2E1	rs2515641	C	T					CT	CC	TT
		Good responders	528 (0.832)	106 (0.167)	0.860 (0.598–1.236)	0.641	0.416	1	88 (0.277)	220 (0.694)	9 (0.028)	0.348	1
		Poor responders	232 (0.811)	54 (0.188)					38 (0.265)	97 (0.678)	8 (0.055)
	COMT	rs4818	G	C					GG	CG	CC
		Good responders	237 (0.375)	395 (0.624)	1.260 (0.932–1.703)	1.747	0.132	1	44 (0.139)	149 (0.471)	123 (0.389)	0.371	1
		Poor responders	95 (0.329)	193 (0.67)					17 (0.118)	61 (0.423)	66 (0.458)
	DRD3	rs6280	C	T					CT	TT	CC
		Good responders	176 (0.286)	438 (0.713)	1.067 (0.768–1.483)	3.28 × 10−4	0.699	1	132 (0.429)	152 (0.500)	22 (0.071)	0.956	1
		Poor responders	81 (0.287)	201 (0.712)					59 (0.418)	70 (0.500)	11 (0.078)
Meta-analysis	UGT1A3	rs6706232	G	A					GG	AA	GA
		Good responders	474 (0.747)	160 (0.252)	1.511 (1.054–2.166)	4.649	0.024	0.591	163 (0.514)	6 (0.018)	148 (0.466)	0.047	1
		Poor responders	286 (0.807)	68 (0.192)					110 (0.621)	1 (0.005)	665 (0.372)
	HTT	rs362267	G	A					GA	GG	AA
		Good responders	466 (0.656)	244 (0.343)	0.829 (0.641–1.072)	1.9	0.152	0.895	158 (0.445)	154 (0.433)	43 (0.121)	0.335	1
		Poor responders	246 (0.615)	154 (0.385)					98 (0.49)	74 (0.37)	28 (0.139)
	CYP2E1	rs2515641	C	T					CC	CT	TT
		Good responders	576 (0.829)	118 (0.17)	1.016 (0.730–1.415)	0.045	0.924	1	238 (0.685)	100 (0.288)	9 (0.025)	0.451	1
		Poor responders	334 (0.835)	66 (0.165)					142 (0.71)	50 (0.25)	8 (0.039)
	COMT	rs4818	G	C					GG	CG	CC
		Good responders	263 (0.387)	415 (0.612)	1.324 (1.012–1.732)	3.556	0.04	0.591	52 (0.153)	159 (0.469)	128 (0.377)	0.166	1
		Poor responders	119 (0.328)	243 (0.671)					19 (0.104)	81 (0.447)	81 (0.447)
	DRD3	rs6280	C	T					CT	TT	CC
		Good responders	187 (0.288)	461 (0.711)	0.858 (0.637–1.156)	1.109	0.314	1	143 (0.441)	159 (0.49)	22 (0.067)	0.548	1
		Poor responders	106 (0.321)	224 (0.678)					78 (0.472)	73 (0.442)	14 (0.084)

Abbreviations: CI, confidence interval; OR, odds ratio; SNP, single-nucleotide polymorphism.

P: Bonferroni-corrected P-value. Bold values signify P-value<0.05.

Table 4

Validation of 17 SNPs associated with risperidone treatment response found in other literatures (no non-genetic factor adjustment)

Cohort	Gene	SNPs	Allele		OR (95% CI)	Chi-square	P1	P1a	Genotype			P2	P2a
Henan	ABCB1	rs1128503	T	C					TT	TC	CC
		Good responders	44 (0.733)	16 (0.266)	0.557 (0.281–1.105)	2.831	0.092	1	16 (0.533)	12 (0.399)	2 (0.066)	0.21	1
		Poor responders	69 (0.605)	45 (0.394)					24 (0.421)	21 (0.368)	12 (0.21)
	HTR3A	rs1176713	C	T					CC	TC	TT
		Good responders	19 (0.741)	41 (0.258)	0.819 (0.419–1.603)	0.337	0.561	1	2 (0.066)	15 (0.5)	13 (0.433)	0.518	1
		Poor responders	39 (0.561)	69 (0.438)					8 (0.148)	23 (0.425)	23 (0.425)
	DRD2	rs1799978	G	A					GA	AA	GG
		Good responders	8 (0.133)	52 (0.866)	0.547 (0.23–1.302)	1.89	0.169	1	8 (0.266)	22 (0.733)	0 (0)	0.327	1
		Poor responders	25 (0.219)	89 (0.78)					21 (0.368)	34 (0.596)	2 (0.035)
		rs1800497	T	C					TT	CC	TC
		Good responders	19 (0.475)	21 (0.525)	0.845 (0.378–1.887)	0.168	0.681	1	9 (0.45)	10 (0.5)	1 (0.05)	0.926	1
		Poor responders	26 (0.433)	34 (0.566)					12 (0.399)	16 (0.533)	2 (0.066)
	MTHFR	rs1801133	T	C					TC	CC	TT
		Good responders	34 (0.566)	26 (0.433)	1.13 (0.6–2.128)	0.144	0.704	1	20 (0.666)	3 (0.099)	7 (0.233)	0.631	1
		Poor responders	68 (0.596)	46 (0.403)					32 (0.561)	7 (0.122)	18 (0.315)
	5HTR6	rs1805054	C	T					CC	TC	TT
		Good responders	40 (0.666)	20 (0.333)	0.973 (0.501–1.891)	0.006	0.937	1	11 (0.366)	18 (0.6)	1 (0.033)	0.995	1
		Poor responders	74 (0.66)	38 (0.339)					20 (0.357)	34 (0.607)	2 (0.035)
	AKT1	rs2494732	C	T					CC	TC	TT
		Good responders	45 (0.725)	17 (0.274)	0.707 (0.358–1.395)	1.001	0.316	1	19 (0.612)	7 (0.225)	5 (0.161)	0.063	0.929
		Poor responders	73 (0.651)	39 (0.384)					23 (0.41)	27 (0.482)	6 (0.107)
		rs3803300	A	G					AA	GG	GA
		Good responders	43 (0.716)	17 (0.283)	0.628 (0.319–1.237)	1.818	0.177	1	15 (0.5)	2 (0.066)	13 (0.433)	0.388	1
		Poor responders	70 (0.614)	44 (0.385)					21 (0.368)	8 (0.14)	28 (0.491)
	HRH3	rs3787429	C	T					CC	TC	TT
		Good responders	42 (0.7)	18 (0.3)	0.756 (0.377–1.514)	0.623	0.429	1	15 (0.5)	12 (0.399)	3 (0.099)	0.707	1
		Poor responders	60 (0.638)	34 (0.361)					19 (0.404)	22 (0.468)	6 (0.127)
	5HTR2C	rs3813928	G	A					GG	AA	GA
		Good responders	51 (0.85)	9 (0.15)	1.006 (0.419–2.418)	2.38 × 10−4	0.987	1	23 (0.766)	2 (0.066)	5 (0.166)	1	1
		Poor responders	97 (0.85)	17 (0.149)					44 (0.771)	4 (0.07)	9 (0.157)
		rs3818929	C	T					CC	TC	TT
		Good responders	49 (0.816)	11 (0.183)	0.642 (0.294–1.399)	1.25	0.263	1	20 (0.666)	9 (0.3)	1 (0.033)	0.542	1
		Poor responders	83 (0.741)	29 (0.258)					31 (0.553)	21 (0.375)	4 (0.071)
	HRH4	rs4483927	T	G					TT	GG	TG
		Good responders	50 (0.806)	12 (0.193)	0.357 (0.171–0.744)	7.856	0.005	0.274	19 (0.612)	0 (0)	12 (0.387)	0.015	0.603
		Poor responders	67 (0.598)	45 (0.401)					20 (0.357)	9 (0.16)	27 (0.482)
	BDNF	rs6265	G	A					GG	AA	GA
		Good responders	28 (0.466)	32 (0.533)	1.183 (0.632–2.213)	0.278	0.597	1	4 (0.133)	6 (0.199)	20 (0.666)	0.079	0.929
		Poor responders	58 (0.508)	56 (0.491)					17 (0.298)	16 (0.28)	24 (0.421)
	5HTR1A	rs6295	C	G					CG	GG	CC
		Good responders	35 (0.583)	25 (0.416)	2.417 (1.231–4.745)	6.748	0.009	0.274	17 (0.566)	4 (0.133)	9 (0.3)	0.02	0.603
		Poor responders	88 (0.771)	26 (0.228)					18 (0.315)	4 (0.07)	35 (0.614)
	5HTR2A	rs6311	T	C					TT	CC	TC
		Good responders	29 (0.467)	33 (0.532)	1.411 (0.757–2.629)	1.178	0.277	1	5 (0.161)	7 (0.225)	19 (0.612)	0.423	1
		Poor responders	62 (0.553)	50 (0.446)					16 (0.285)	10 (0.178)	30 (0.535)
		rs6313	T	C					TT	CC	TC
		Good responders	29 (0.483)	31 (0.516)	1.274 (0.681–2.384)	0.577	0.447	1	5 (0.166)	6 (0.199)	19 (0.633)	0.481	1
		Poor responders	62 (0.543)	52 (0.456)					16 (0.28)	11 (0.192)	30 (0.526)
	COMT	rs9606186	G	C					GG	CG	CC
		Good responders	50 (0.862)	8 (0.137)	0.387 (0.164–0.91)	4.938	0.026	0.511	22 (0.758)	6 (0.206)	1 (0.034)	0.041	0.811
		Poor responders	75 (0.707)	31 (0.292)					25 (0.471)	25 (0.471)	3 (0.056)
Shanghai	ABCB1	rs1128503	T	C					TC	CC	TT
		Good responders	392 (0.626)	234 (0.373)	1.025 (0.766–1.373)	0.029	0.864	1	158 (0.504)	38 (0.121)	117 (0.373)	0.134	1
		Poor responders	177 (0.632)	103 (0.367)					57 (0.407)	23 (0.164)	60 (0.428)
	HTR3A	rs1176713	T	C					TC	CC	TT
		Good responders	479 (0.746)	163 (0.253)	1.165 (0.84–1.615)	0.842	0.358	1	123 (0.383)	20 (0.062)	178 (0.554)	0.51	1
		Poor responders	226 (0.773)	66 (0.226)					48 (0.328)	9 (0.061)	89 (0.609)
	DRD2	rs1799978	A	G					AA	GG	AG
		Good responders	505 (0.786)	137 (0.213)	1.182 (0.83–1.683)	0.865	0.352	1	199 (0.619)	15 (0.046)	107 (0.333)	0.577	1
		Poor responders	231 (0.813)	53 (0.186)					93 (0.654)	4 (0.028)	45 (0.316)
		rs1800497	T	C					TT	CC	TC
		Good responders	216 (0.395)	330 (0.604)	0.974 (0.709–1.338)	0.025	0.872	1	52 (0.19)	109 (0.399)	112 (0.41)	0.804	1
		Poor responders	88 (0.389)	138 (0.61)					19 (0.168)	44 (0.389)	50 (0.442)
	MTHFR	rs1801133	T	C					TC	CC	TT
		Good responders	256 (0.408)	370 (0.591)	1.184 (0.886–1.582)	1.309	0.252	1	158 (0.504)	106 (0.338)	49 (0.156)	0.491	1
		Poor responders	104 (0.368)	178 (0.631)					70 (0.496)	54 (0.382)	17 (0.12)
	5HTR6	rs1805054	C	T					CC	TC	TT
		Good responders	454 (0.727)	170 (0.272)	1.05 (0.761–1.449)	0.09	0.763	1	159 (0.509)	136 (0.435)	17 (0.054)	0.892	1
		Poor responders	202 (0.737)	72 (0.262)					71 (0.518)	60 (0.437)	6 (0.043)
	AKT1	rs2494732	C	T					CC	TC	TT
		Good responders	475 (0.73)	175 (0.269)	1.038 (0.759–1.419)	0.055	0.813	1	174 (0.535)	127 (0.39)	24 (0.073)	0.485	1
		Poor responders	217 (0.738)	77 (0.261)					77 (0.523)	63 (0.428)	7 (0.047)
		rs3803300	A	G					AG	AA	GG
		Good responders	398 (0.633)	230 (0.366)	1.207 (0.895–1.628)	1.523	0.217	1	152 (0.484)	123 (0.391)	39 (0.124)	0.39	1
		Poor responders	188 (0.676)	90 (0.323)					60 (0.431)	64 (0.46)	15 (0.107)
	HRH3	rs3787429	C	T					CC	TC	TT
		Good responders	355 (0.618)	219 (0.381)	1.015 (0.748–1.376)	0.009	0.922	1	113 (0.393)	129 (0.449)	45 (0.156)	0.533	1
		Poor responders	158 (0.622)	96 (0.377)					47 (0.37)	64 (0.503)	16 (0.125)
	5HTR2C	rs3813928	G	A					GG	AA	AG
		Good responders	544 (0.855)	92 (0.144)	1.296 (0.847–1.982)	1.442	0.229	1	254 (0.798)	28 (0.088)	36 (0.113)	0.059	1
		Poor responders	253 (0.884)	33 (0.115)					115 (0.804)	5 (0.034)	23 (0.16)
		rs3818929	C	T					CC	TC	TT
		Good responders	555 (0.851)	97 (0.148)	1.326 (0.873–2.013)	1.764	0.184	1	255 (0.782)	45 (0.138)	26 (0.079)	0.065	1
		Poor responders	258 (0.883)	34 (0.116)					116 (0.794)	26 (0.178)	4 (0.027)
	HRH4	rs4483927	T	G					TT	TG	GG
		Good responders	404 (0.627)	240 (0.372)	0.968 (0.728–1.288)	0.047	0.826	1	131 (0.406)	142 (0.44)	49 (0.152)	0.977	1
		Poor responders	181 (0.619)	111 (0.38)					58 (0.397)	65 (0.445)	23 (0.157)
	BDNF	rs6265	A	G					AG	AA	GG
		Good responders	320 (0.503)	316 (0.496)	0.959 (0.725–1.27)	0.081	0.774	1	156 (0.49)	82 (0.257)	80 (0.251)	0.66	1
		Poor responders	139 (0.492)	143 (0.507)					63 (0.446)	38 (0.269)	40 (0.283)
	5HTR1A	rs6295	C	G					CG	CC	GG
		Good responders	482 (0.75)	160 (0.249)	1.005 (0.727~1.39)	0.001	0.974	1	118 (0.367)	182 (0.566)	21 (0.065)	0.91	1
		Poor responders	212 (0.751)	70 (0.248)					54 (0.382)	79 (0.56)	8 (0.056)
	5HTR2A	rs6311	C	T					CC	TC	TT
		Good responders	279 (0.431)	367 (0.568)	1.064 (0.803–1.41)	0.188	0.664	1	65 (0.201)	149 (0.461)	109 (0.337)	0.872	1
		Poor responders	120 (0.416)	168 (0.583)					26 (0.18)	68 (0.472)	50 (0.347)
		rs6313	C	T					CC	TC	TT
		Good responders	277 (0.436)	357 (0.563)	1.123 (0.845–1.492)	0.648	0.42	1	65 (0.205)	147 (0.463)	105 (0.331)	0.659	1
		Poor responders	116 (0.408)	168 (0.591)					24 (0.169)	68 (0.478)	50 (0.352)
	COMT	rs9606186	C	G					CG	GG	CC
		Good responders	159 (0.258)	457 (0.741)	0.869 (0.634–1.192)	0.749	0.386	1	111 (0.36)	173 (0.561)	24 (0.077)	0.069	1
		Poor responders	80 (0.285)	200 (0.714)					66 (0.471)	67 (0.478)	7 (0.05)
Meta-analysis	ABCB1	rs1128503	T	C					TT	TC	CC
		Good responders	436 (0.635)	250 (0.364)	0.953 (0.737–1.231)	0.134	0.713	1	133 (0.387)	170 (0.495)	40 (0.116)	0.038	0.744
		Poor responders	246 (0.624)	148 (0.375)					84 (0.426)	78 (0.395)	35 (0.177)
	HTR3A	rs1176713	C	T					CC	TC	TT
		Good responders	182 (0.259)	520 (0.74)	0.983 (0.743–1.3)	0.013	0.906	1	22 (0.062)	138 (0.393)	191 (0.544)	0.485	1
		Poor responders	105 (0.262)	295 (0.737)					17 (0.085)	71 (0.355)	112 (0.559)
	DRD2	rs1799978	G	A					GA	AA	GG
		Good responders	145 (0.206)	557 (0.793)	1.067 (0.785–1.452)	0.175	0.675	1	115 (0.327)	221 (0.629)	15 (0.042)	0.76	1
		Poor responders	78 (0.195)	320 (0.804)					66 (0.331)	127 (0.638)	6 (0.03)
		rs1800497	T	C					TT	CC	TC
		Good responders	235 (0.401)	351 (0.598)	0.989 (0.741–1.321)	0.004	0.945	1	61 (0.208)	119 (0.406)	113 (0.385)	0.905	1
		Poor responders	114 (0.398)	172 (0.601)					31 (0.216)	60 (0.419)	52 (0.363)
	MTHFR	rs1801133	T	C					TC	CC	TT
		Good responders	290 (0.422)	396 (0.577)	0.953 (0.742–1.224)	0.138	0.71	1	178 (0.518)	109 (0.317)	56 (0.163)	0.915	1
		Poor responders	172 (0.434)	224 (0.565)					102 (0.515)	61 (0.308)	35 (0.176)
	5HTR6	rs1805054	C	T					CC	TC	TT
		Good responders	494 (0.722)	190 (0.277)	0.965 (0.731–1.273)	0.063	0.801	1	170 (0.497)	154 (0.45)	18 (0.052)	0.656	1
		Poor responders	276 (0.715)	110 (0.284)					91 (0.471)	94 (0.487)	8 (0.041)
	AKT1	rs2494732	C	T					CC	TC	TT
		Good responders	520 (0.73)	192 (0.269)	0.923 (0.703–1.211)	0.333	0.563	1	193 (0.542)	134 (0.376)	29 (0.081)	0.276	1
		Poor responders	290 (0.714)	116 (0.285)					100 (0.492)	90 (0.443)	13 (0.064)
		rs3803300	A	G					AA	GG	GA
		Good responders	441 (0.64)	247 (0.359)	1.078 (0.831–1.399)	0.322	0.57	1	138 (0.401)	41 (0.119)	165 (0.479)	0.75	1
		Poor responders	258 (0.658)	134 (0.341)					85 (0.433)	23 (0.117)	88 (0.448)
	HRH3	rs3787429	C	T					CC	TC	TT
		Good responders	397 (0.626)	237 (0.373)	1.001 (0.763–1.311)	6.18 × 10−5	0.993	1	128 (0.403)	141 (0.444)	48 (0.151)	0.533	1
		Poor responders	218 (0.626)	130 (0.373)					66 (0.379)	86 (0.494)	22 (0.126)
	5HTR2C	rs3813928	G	A					GG	AA	GA
		Good responders	595 (0.854)	101 (0.145)	1.188 (0.825–1.709)	0.865	0.352	1	277 (0.795)	30 (0.086)	41 (0.117)	0.094	1
		Poor responders	350 (0.874)	50 (0.125)					159 (0.794)	9 (0.045)	32 (0.159)
		rs3818929	C	T					CC	TC	TT
		Good responders	604 (0.848)	108 (0.151)	0.967 (0.69–1.356)	0.035	0.849	1	275 (0.772)	54 (0.151)	27 (0.075)	0.021	0.619
		Poor responders	341 (0.844)	63 (0.155)					147 (0.727)	47 (0.232)	8 (0.039)
	HRH4	rs4483927	T	G					TT	GG	TG
		Good responders	454 (0.643)	252 (0.356)	0.882 (0.685–1.136)	0.942	0.331	0.274	150 (0.424)	49 (0.138)	154 (0.436)	0.632	0.603
		Poor responders	248 (0.613)	156 (0.386)					78 (0.386)	32 (0.158)	92 (0.455)
	BDNF	rs6265	G	A					GG	AA	GA
		Good responders	344 (0.494)	352 (0.505)	1.054 (0.824–1.349)	0.179	0.672	1	84 (0.241)	88 (0.252)	176 (0.505)	0.3	0.929
		Poor responders	201 (0.507)	195 (0.492)					57 (0.287)	54 (0.272)	87 (0.439)
	5HTR1A	rs6295	C	G					CG	GG	CC
		Good responders	517 (0.736)	185 (0.263)	1.118 (0.841–1.486)	0.592	0.441	0.274	135 (0.384)	25 (0.071)	191 (0.544)	0.745	0.603
		Poor responders	300 (0.757)	96 (0.242)					72 (0.363)	12 (0.06)	114 (0.575)
	5HTR2A	rs6311	T	C					TT	CC	TC
		Good responders	396 (0.559)	312 (0.44)	1.065 (0.832–1.365)	0.255	0.613	1	114 (0.322)	72 (0.203)	168 (0.474)	0.799	1
		Poor responders	230 (0.575)	170 (0.425)					66 (0.33)	36 (0.18)	98 (0.49)
		rs6313	T	C					TT	CC	TC
		Good responders	386 (0.556)	308 (0.443)	1.092 (0.851–1.401)	0.484	0.486	1	110 (0.317)	71 (0.204)	166 (0.478)	0.713	1
		Poor responders	230 (0.577)	168 (0.422)					66 (0.331)	35 (0.175)	98 (0.492)
	COMT	rs9606186	G	C					GG	CG	CC
		Good responders	507 (0.752)	167 (0.247)	0.816 (0.615–1.081)	2.008	0.156	0.511	195 (0.578)	117 (0.347)	25 (0.074)	0.017	0.619
		Poor responders	275 (0.712)	111 (0.287)					92 (0.476)	91 (0.471)	10 (0.051)

Abbreviations: CI, confidence interval; OR, odds ratio; SNP, single-nucleotide polymorphism.

P: Bonferroni-corrected P-value. Bold values signify P-value<0.05.

The discovery and validation of DNA methylation biomarkers associated with risperidone response

The PCR primers for bisulfate-treated amplification were listed (Supplementary Table 4). Significant CpG sites within CYP3A4 gene promoter region or CYP2D6 gene body were found associated with risperidone treatment (Figures 2, 3, 4), whereas no significant CpG sites in HTR2A, ABCB1, DRD2 gene promoters, respectively (data not shown).

Figure 2

Methylation rate of CYP3A4 and CYP2D6 gene promoter and gene-coding region in Henan Cohort. P*: P-value after Bonferroni correction; CYP3A4_CpG_-36: 36 bases ahead of CYP3A4 transcriptional start site; CYP2D6_CpG_36: 36 bases after CYP2D6 transcriptional start site.

Figure 3

Methylation rate of CYP3A4 and CYP2D6 gene promoter and gene-coding region in Shanghai Cohort as a replication group for validation of the results found in Henan cohort.

Figure 4

Meta-analysis about the methylation rate of CYP3A4 and CYP2D6 gene promoter and gene-coding region in Henan and Shanghai Cohort. Ten CpG sites were significantly related with risperidone treatment response; whereas after Bonferroni correction, there remained seven significant CpG sites between good responders and poor responders.

Multivariate interactions analysis of genetic variant and DNA methylation biomarkers associated with risperidone response

MDR software was used for detecting factor–factor interactions (including all SNPs, CpG rates, sex, age, weight and so on) in genetic case–control studies as having great advantages vs the conventional statistical approaches. We combined all the SNPs and CpG sites into MDR and built two-, three- and four-level interaction models with MDR software, then we found a best four-locus model (CYP3A4_CpG_-82:-86+rs6280+rs1800497+rs6265)with a testing-balanced accuracy of 58.94% and a cross-validation consistency of 7/10. The permutation testing showed that the four-locus model was significant with P=0.044 (Table 5).

Table 5

Multivariate interactions analysis results

Antipsychotics	Model	Training balance accuracy	Testing balance accuracy	P-valuea	Cross-validation consistency
Risperidone	rs6280+rs3787429	0.6733	0.5283	0.32	5/10
	CYP2D6_CpG_242:244:250+CYP3A4_CpG_-82:-86 +rs1800497	0.7583	0.5489	0.16	3/10
	CYP3A4_CpG_-82:-86+rs6280+rs1800497+rs6265	0.8685	0.5894	0.038	7/10

P-value based on 1000 permutations. Bold values signify P-value<0.05.

Discussion

To provide the most effective and safe treatment for each schizophrenic patient, understanding the internal and external factors affecting response to antipsychotic drugs is important. Most studies put emphasis on either genetic or environmental aspect alone instead of combining both pharmacogenetics and pharmacoepigenetics. This prompted us to perform a more systematic study using two methods: next-generation sequence of 100 risperidone treatment response-related genes and epigenetic research five representative of these genes, including drug metabolic enzyme genes CYP450s (phase I), Uridine 5′-diphosphate-glucuronosyltransferases (UGTs; phase II), mostly studied receptor genes, other potentially related genes and five risperidone tightly associated ADME gene promoters or coding regions, which might optimize clinical treatment in antipsychotic drugs treatment. More importantly, most of the reported findings remain inconclusive and require either clinical examination or a larger clinical sample to validate, so we selected 17 other SNPs from literatures for validation and comparison of our NGS results based on Chinese Han population background.

NGS study

In NGS study, we found 10 candidate SNPs significantly related with risperidone treatment response in six genes (Table 2). We also confirmed NGS results by MassARRAY genotyping platform with an accuracy of 94.88% accordance in average (data not shown), whereas the significance of 10 candidate SNPs discovered by NGS platform no longer demonstrated significance after Bonferroni correction. As Henan cohort only contained less than 100 cases, we recruited another Shanghai cohort for the validation of NGS results and implemented a meta-analysis by combining the two cohorts (Table 3). Ultimately, we discovered two rare variant SNPs with 1%< MAF <5%, though both of them had limited impact in downstream protein functions predicted by protein online databases. We could not ignore the high effect score predicted by PolyPhen2; in addition, we found two SNPs, rs6706232 in UGT1A3 and rs4818 in catechol-O-methyltransferase (COMT) were significantly associated with risperidone treatment in meta-analysis (P=0.024 and P=0.04, respectively).

UGT1A3

UGTs are phase II drug metabolism enzymes in human tissues, containing nine functional isoforms including (UGT1A1, UGT1A3–UGT1A10) and four pseudogenes (UGT1A2, UGT1A11–UGT1A13). Among these enzymes, UGT1A3 has a critical role in endobiotic and xenobiotic compounds metabolism by catalyzing the glucuronidation of endogenous compounds such as bilirubin, bile acids, thyroid hormone, steroid hormones and substantial exogenous substrates such as many therapeutic drugs, heterocyclic and polycyclic hydrocarbons, and heterocyclic amines,[15, 16, 17] Polymorphisms in UGT1A3 genes (http://www.pharmacogenomics.pha.ulaval.ca/cms/ugt_alleles/) could have significant influence on metabolism of endogenous compounds like bilirubin or variability in response to irinotecan among other drugs.[18, 19] For example, effect of rs6706232 on OTS167 (a novel synthetic anticancer agent molecule undergoing clinical development) glucuronidation formation rates was merely modest suggesting this SNP may not significantly contribute to OTS167 clearance.[20] In our study, rs6706232 was first found related with risperidone treatment in mainland Chinese Han population, allele A accounted for a higher proportion in good responders compared with poor responders at 33.3%, 24.6%, 25.2% in the Henan, Shanghai and total cohort, respectively, which means the G>A mutation might have certain impact on risperidone therapy. Recent research found that rs6706232 was significantly associated with UGT1A3 gene transcription and activity.[21] Whether enzyme UGT1A3 directly, or its derivatives, affects risperidone response remains elusive although aromatic heteropolycyclic risperidone was cataloged a potential target catalyzed by UGT1A3 enzyme.

COMT polymorphisms

COMT, first discovered in 1957,[22] has always been regarded as susceptibility gene in schizophrenia for its function of degrading catecholamine such as dopamine, epinephrine and norepinephrine.[23, 24] Polymorphisms in COMT affected enzyme activity, for example, well-studied variant rs4680 G>A (substitution of valine to methionine at codon 158) resulted in decreased enzyme activity, which contributed to inefficient catalyzation and accumulation of dopamine, thus inducing positive symptoms in schizophrenia or even drug treatment failure.[25, 26] Gupta et al.[27] found that haplotype rs4818 and rs4680 in exon 4 was related with risperidone treatment response (P=0.028) in 398 schizophrenia patients and 241 healthy individuals from a homogeneous South Indian population. Recently, our group found that rs4818 was significantly associated with quetiapine drug response (P=0.00081) in the Shanghai cohort, but there was no correlation with risperidone drug response (P=0.515),[28] whereas in this study, we confirmed the significance of rs4818 in 288 risperidone-treated schizophrenia patients (P=0.04) with similar G/C allele frequency in good responders or poor responders. The difference might partially result from different sample size or drug metabolism pathways and indications. Risperidone has a higher pituitary D2 receptor occupancy, which is suitable for positive symptoms than quetiapine, whereas quetiapine has various affinities for cerebral serotonergic(5-HT), D2 dopaminergic, histaminergic(H1) and adrenergic receptors, including much wilder indications than risperidone that involved more genetic pathways.[29, 30] As a synonymous variant, rs4818 may have a milder influence on COMT enzyme function compared with rs4860. Another interesting result is that most significant polymorphisms were simultaneously found in COMT including rs4633, rs4680 and rs4818 (based on NGS), rs9606186 (based on MassARRAY) with P=0.018, 0.028, 0.039, 0.041 in the Henan cohort. The discovery of rs9606186, rs4680 and rs4818 also confirmed previous findings in Chinese patients carried out by our group.[31] Taken together, these results suggest that COMT polymorphisms contribute greatly to risperidone treatment.

Replication of 17 SNPs

Replication of other reported 17 SNPs was carried out preliminarily in the Henan cohort and we found three candidate SNPs including rs4483927 in HRH4, rs6295 in 5HTR1A, rs9606186 in COMT with P=0.005, 0.009, 0.026, respectively, while rejected by multiple testing corrections. However, after validation in the Shanghai cohort for replication and meta-analysis, none of the three SNPs remained significant. This might partially attribute to different population background such as geography, nutrition[32, 33] and so on, and the fact of this study using a small sample size also accounts for the inconsistent results in our study and others. This promotes us to further expand our sample size to nearly 600 qualified subjects for discovering the effect markers.

Comparison of previous results by our group

We compared this study with other pharmacogenomics results found by our group in the last decade to investigate the accordance of different studies, which might make complement or new clues to pharmacogenomics research. Previously, our group found no significant SNPs in CYP2D6, CYP3A4, CYP2E1 associated with risperidone treatment response,[34, 35, 36] however, SNPs in CYP2D6, CYP3A4 and CYP2E1 in this study were not overlapped and rs2515641 in CYP2E1 was found significantly related with risperidone treatment response in the Henan cohort using the NGS method (P=0.007). Rs2515641, located in the eighth exon of CYP2E1 as a synonymous mutation (Phe421Phe), was first uncovered in the scanning of CYP2E1 in Chinese mainland Han population in 2010 (ref. 37; MAF=15.1%) and might relate to acute rejection in kidney transplantation recipients.[38] However, the significance was not observed in the Henan and Shanghai cohort or meta-analysis when conducted in MassARRAY platform (P>0.05, data not shown). Besides, Xing, et al.[39] found that the wild-type TT of rs1128503 in ABCB1 (synonymous mutation) carriers had better risperidone treatment response than other genotype carriers in 130 Chinese schizophrenia patients (F=3.976, P=0.021). But we failed to validate this SNP in meta-analysis, which might reflect the limited role in risperidone treatment response though the brain entry of risperidone and 9-OH-risperidone is greatly limited by ABCB1 product: p-glycoprotein. Meanwhile, loci rs6280 (Ser9Gly) of DRD3 gene has shown significant association with risperidone treatment response in the Henan cohort by the NGS method (adjusted P=0.031) or MassARRAY platform (adjusted P=0.015), and more allele C carriers in poor responders group than good responders with 41.8% vs 24.2% (Henan cohort, NGS platform), 52% vs 32% (Henan cohort, MassARRAY), which might be related with therapy response, but this result was inconsistent with Xuan et al.[40]

Epigenetic study

The evidence that epigenetic events can have an important role in regulating the expression of drug ADME genes, drug transporters, nuclear receptors and drug targets strongly implied that interindividual differences in their epigenetic status can contribute to the clinically observed variability in drug response,[41] which could not be explained by genetic polymorphisms alone. We chose gene promoter or gene body region[42] of five genes (CYP3A4, CYP2D6, HTR2A, DRD2, ABCB1) to investigate the correlation between their methylation status and risperidone treatment response in the Henan and Shanghai cohorts and found seven significant CpG sites related with drug treatment after Bonferroni correction (Figure 4). The CpG sites in CYP3A4 promoter had a methylation rate ranging from 20% to approximately 100%. We also observed that every single CpG site has slight difference between the Henan and Shanghai cohort excepting CYP3A4_CpG_-258 suggesting that the geography factor might have little effect in methylation status during risperidone treatment therapy period in both good and poor responders. The initial two positive CpG sites (CYP3A4_CpG_-36, CYP3A4_CpG_-296, before Bonferroni correction) found in the Henan cohort failed to be validated in Shanghai replication group while another novel positive CpG site (CYP3A4_CpG_-390, P=0.048, P=0.78, before and after Bonferroni correction, respectively) was found. However, in the meta-analysis, CYP3A4_CpG_-36 and CYP3A4_CpG_-296 were validated ultimately and another two novel CpG sites were found (P=0.0014, P=0.0013, P=0.0014, P=0.028, after Bonferroni correction). Kacevska et al.[43] investigated the methylation status in ~12 kb CYP3A4 regulation region including the proximal promoter, XREM and CLEM4 and in separate C/EBP and HNF4α-binding regions in only 79 subjects and found CYP3A4_CpG_-383 showing significant Spearman’s rank coefficient between adjacent CpG sites. Their research covered a wider region of CYP3A4 than ours, but we excavated much more thoroughly in the CYP3A4 promoter region, which contains 16 sites vs 2 sites for analysis. Besides, the four candidate CpG sites were all located before CYP3A4_CpG_-383 and the methylation rates of CYP3A4_CpG_-36, CYP3A4_CpG_-258, CYP3A4_CpG_-296 were higher in poor responders than good responders in both the Henan and Shanghai cohort, possibly suggesting that the inhibition of CYP3A4 protein expression by high methylation rate might result in low efficiency in metabolizing risperidone. Though CYP2D6 enzyme accounts for only 1.3%–4.3% of all hepatic CYPs but metabolites ~20% medications in human liver,[44] little knowledge is accessible about methylation regulation in CYP2D6 promoter and gene body regions. Park and colleagues found the methylation frequency in the gene body region of CYP2D6 containing 32 CpG islands, the methylation frequency was 45.5% and 90.3% in human embryonic stem cell-derived hepatocytes and primary hepatocytes, respectively, which exhibits a dynamic methylation pattern change. Our specific CpG sites were not overlapped with that found by Park and colleagues, and we found three positive CYP2D6_CpG sites (P =0.012, P =0.00076, P =0.034) in meta-analysis in 15 CpG sites. The methylation status of CYP2D6 gene body was highly varied and irregular across CpG site. We also observed that the methylation rate in poor responders was relatively higher than good responders in both the Henan and Shanghai cohort for all the three significant CpG sites, methylation rate of CYP2D6_242:244:250 in poor responders reached 90% which means CYP2D6 enzyme might be repressed and induced undesired risperidone treatment response. Although these findings may still be preliminary and still require either clinical examination or a larger clinical sample group, such translatable pharmacological effects demonstrate the potential of epigenetic phenomenon in explaining the interindividual differences in drug treatment outcome. Characterizing such epigenetic marks and developing noninvasive approaches to examine them holds a promising tool in the effective treatment of schizophrenia.

MDR analysis

Many factors affect risperidone treatment response. Genetic polymorphism or epigenetic regulation alone may be incapable of explaining drug response. A previous study observed this combination effect in drug response.[45] In this study, we conducted a combined analysis including all the 27 SNPs and the methylation status of genes that encoded drug metabolism enzyme, drug transporters and neurotransmitter receptors and so on using MDR software. The four-locus model (CYP3A4_CpG_-82:-86+rs6280 (CYP1B1) + rs1800497 (DRD2) + rs6265 (BDNF)) was regarded as the best significant model (P=0.038), which might provide new clues in predicting drug treatment response. Description of the best four-locus model could be found in the Supplementary Materials. This model revealed that mutual interactions of genetic and epigenetic factors might be related with risperidone treatment efficacy in Chinese Han schizophrenia patients.

Conclusion

In conclusion, pharmacogenetics studies of antipsychotic drugs are promising despite many challenges. Our results may push the field closer to routine clinical utilization of pharmacogenetics testing to maximize therapeutic effects and minimize adverse effects. We found genetic and epigenetic biomarkers in risperidone treatment efficacy due to our systematical study design, but there are still some shortcomings in our research and more samples are suggested to be recruited to strengthen our results in future researchs. Progress in genomic technology and bioinformatics, larger sample sizes, better phenotype characterization and precise study design will help to promote antipsychotic pharmacogenetics to its next level.