Drug Response-Related DNA Methylation Changes in Schizophrenia, Bipolar Disorder, and Major Depressive Disorder

Pharmacotherapy is the most common treatment for schizophrenia (SCZ), bipolar disorder (BD), and major depressive disorder (MDD). Pharmacogenetic studies have achieved results with limited clinical utility. DNA methylation (DNAm), an epigenetic modification, has been proposed to be involved in both the pathology and drug treatment of these disorders. Emerging data indicates that DNAm could be used as a predictor of drug response for psychiatric disorders. In this study, we performed a systematic review to evaluate the reproducibility of published changes of drug response-related DNAm in SCZ, BD and MDD. A total of 37 publications were included. Since the studies involved patients of different treatment stages, we partitioned them into three groups based on their primary focuses: (1) medication-induced DNAm changes (n = 8); (2) the relationship between DNAm and clinical improvement (n = 24); and (3) comparison of DNAm status across different medications (n = 14). We found that only BDNF was consistent with the DNAm changes detected in four independent studies for MDD. It was positively correlated with clinical improvement in MDD. To develop better predictive DNAm factors for drug response, we also discussed future research strategies, including experimental, analytical procedures and statistical criteria. Our review shows promising possibilities for using BDNF DNAm as a predictor of antidepressant treatment response for MDD, while more pharmacoepigenetic studies are needed for treatments of various diseases. Future research should take advantage of a system-wide analysis with a strict and standard analytical procedure.

Therapeutics for SCZ, BD, and MDD are generally based on similar classes of molecules, targeting similar pathways, with distinct doses and proper combinations. However, drug selection is clinically subjective and treatment typically requires weeks of symptom evaluation to determine treatment efficacy. A large proportion of patients fail to respond to first-line drug treatment. For example, the Sequenced Treatment Alternatives to Relieve Depression trial (STAR * D) study reported that only 35% of patients remit after their primary antidepressant treatment (Rush et al., 2006;Katz et al., 2012). This illustrates the importance of developing biomarkers that can support decisions regarding optimal drug choice, and identify likely poor responders as quickly as possible. This emphasizes the need for a better biomarker-based stratification of patients that could facilitate treatment planning.
Pharmacogenetic approaches for guiding the treatment of psychiatric disorders has been a rapidly expanding area of research in the last two decades (Nelson et al., 2016). The main hypothesis of these studies was that genetic variants could predict the influences of drug treatment. Numerous pharmacogenetic studies have investigated the genetic contribution to the treatment response of these disorders (Baum et al., 2008;Brandl et al., 2016;Shi et al., 2017;Yu et al., 2018). Clinical practice for pharmacogenetic findings is emerging, as a few of the Food and Drug Administration (FDA) approved drug labels changed for SCZ and MDD. Although psychiatric disorders are thought to be highly heritable, gene-environment interactions are relevant, and results of pharmacogenetic studies have been inconsistent, limiting the clinical utilization of this approach.
Emerging evidence suggests that epigenetic marks could be used as a predictor of drug response for psychiatric disorders (Chan and Baylin, 2012;Heerboth et al., 2014). Dysregulation of epigenetic events can be pathological, and associated with SCZ, BD, and MDD (Mill et al., 2008;Sabunciyan et al., 2012;Guintivano et al., 2013;Chen et al., 2014;Liu et al., 2018). This implies that biological pathways and cellular processes are under the impact of epigenome status. Unlike genetics, epigenetic status is dynamic and could better reflect various environmental events during disease progress and drug treatment. Because of the reversibility of epigenetic events, we postulate that modulation of epigenetic regulators could be valuable for therapeutic potential. Pharmacoepigenetics may be more promising for guiding the treatment of psychiatric disorders than pharmacogenetics.
The recent advent of interest in pharmacoepigenetics of psychiatric disorders is an alternative research avenue. The epigenetic modification known as DNA methylation (DNAm), one of major epigenetic modifications, is the product of the interaction between genetic variants and environmental influence, and for this reason might be a better predictor of treatment outcomes (Kubota et al., 2012;Moore et al., 2013). Measurable in peripheral blood, alterations of DNAm have been reported to be involved in both the pathology and drug treatment of SCZ, BD, and MDD (Schroeder et al., 2012;Vialou et al., 2013;Jaffe et al., 2016). Drugs may exert their effects by reversing these DNAm deregulations. Moreover, some FDA-approved drugs, like clozapine and sulpiride, have demonstrated the ability to activate brain DNA demethylation (Dong et al., 2008). Therefore, DNAm is a promising molecular approach to study mechanisms and prediction of drug response in psychiatric disorders (Goud Alladi et al., 2018;Lisoway et al., 2018). However, the reproducibility of those results, crucial for their clinical application, has not been thoroughly evaluated.
We systematically reviewed the literature on DNAm and drug response in SCZ, MDD, and BD to determine if epigenetic markers have produced any reproducible finding as predictors of drug response. In this review, we focused on studies providing information about the effect of medications on DNAm and the utility of DNAm as a predictor of drug treatment for major psychiatric disorders (Figure 1). First, we evaluated the reproducibility of reported changes of DNAm-related genes associated with drug response in each disorder. We also evaluated whether the results of candidate gene studies could be reproduced in genome-wide studies. Finally, we discuss limitations of existing studies and provide recommendations for future research.

METHODS
We systematically reviewed published studies that examined DNAm in relation to drug response in SCZ, BD, or MDD, using the PubMed, Web of Science, and MEDLINE with the following search builder: (((((((pharmacogenetic) OR pharmacogenomic) OR pharmacoepigenetic) OR drug response) OR treatment) OR drug efficiency) OR #DRUG#) AND ((((DNA methylation) OR methylation) OR epigenetic) OR EWAS) AND (#DISORDER#). #DRUG# is the major drug for each disorder: antipsychotic, mood stabilizer, lithium, antidepressant, selective serotonin reuptake inhibitor (SSRI) and selective norepinephrine reuptake inhibitor (SNRI). #DISORDER# is one of the major psychiatric disorders: SCZ, BD and MDD. The inclusion criteria were the following: (1) published in English before April 2020; (2) assessed drug response on DNAm levels; (3) sample size larger than 10. Studies on animal and cell culture models were excluded. We discarded one study testing for drug-related co-methylation modules for BD, since it did not directly compare the efficacy of the individual drugs and was incompatible with the other studies. The literature selection followed the standard by Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) (Liberati et al., 2009).
We included 37 publications about drug response-related DNAm in SCZ, BD, and MDD (Figure 2). Studies included SCZ (number of studies, n = 7), BD (n = 6), MDD (n = 18), and cross-disorder studies (BD and MDD, n = 2; SCZ and BD, n = 4). Since the studies involved patients of different treatment stages, we partitioned them into three groups based on their primary focuses (Figure 3): (1) medication-induced DNAm changes; (2) the relationship between DNAm and clinical improvement; and (3) comparison of DNAm status across different medications. To evaluate the reproducibility of published changes of DNAmrelated genes associated with treatment response in psychiatric disorders, we examined the stability of the reported changes in each disorder. Here, we defined stability as the presence of significant results (p < 0.05) that were consistent in more than one independent study. We also evaluated whether the FIGURE 1 | Schematic representation of traditional medicines vs. personalized medicine using pharmacoepigenetic approaches.
results of candidate gene studies could be reproduced in genomewide studies.

KEY CONCEPTS, DEFINITIONS, AND PRINCIPLES IN PHARMACOEPIGENETICS
Pharmacoepigenetics research in psychiatry focuses on the effects of treatments on DNAm and their potential influences on treatment response. The main aim of pharmacoepigenetics is the identification of DNAm of specific genes that are associated with treatment outcomes of psychiatric disorders, with the ultimate goal of translating this information into strategies that can support personalized medicine.

DNA Methylation
DNA methylation is the epigenetic modification by which methyl groups are added to DNA nucleotides, primarily cytosine, and adenine. DNAm is thought to have a prominent influence on the structure and functions of DNA, particularly in the regulation of gene expression. However, DNAm regulates gene expression differently according to its genomic context. Previous studies revealed that DNAm levels of CpG sites near the transcription start sites could repress gene expression levels, while DNAm levels of gene body CpGs could activate the expression of their target genes (Holliday and Pugh, 1975;Hellman and Chess, 2007;Jones, 2012;Moore et al., 2013). Conflicting results were reported regarding intergenic CpGs: some claimed that they might regulate gene expression less frequently (Chen et al., 2014;Wagner et al., 2014;Liu et al., 2018), while others suggested that intergenic CpGs could repress the gene expression (Gaudet et al., 2004;Hutnick et al., 2010;Moore et al., 2013). Though CpGs are the primary sites for DNAm, it has also been observed at non-CpG sites (CpA, CpT, and CpC) in mammals (Pinney, 2014). Non-CpG methylation sites are common in several tissue types, including skeletal muscle (Yan et al., 2011) and brain (Lister et al., 2013;Guo et al., 2014). However, the functions of non-CpG methylation are still largely unknown and information on their function is just beginning to emerge. Our review includes only studies on the DNAm at CpG sites.
Tissues from postmortem samples, preserved clinical specimens from humans, can be used to study DNAm. For study pharmacoepigenetics for psychiatric disorders, peripheral tissues, such as blood, saliva, and buccal cells, are useful for developing clinically informative biomarkers. Psychiatric disorders are considered primarily a brain dysregulation. Therefore, we prefer to investigate genomic or epigenomic features in peripheral tissues that are similar to brain. Given that DNAm signals can be tissue-specific, it is important to establish the pros and cons of different types of tissue for investigation (Mill and Petronis, 2007).
Several studies investigated correlations of DNAm in blood and multiple regions of the brain, with inconsistent results (Davies et al., 2012;Horvath et al., 2012;Kaminsky et al., 2012;Masliah et al., 2013;Walton et al., 2016). Studies also indicated that DNAm in saliva and buccal cells were more informative surrogate tissues than blood (Lowe et al., 2013;Smith et al., 2015). To address this issue, Braun et al. (2019) performed genome-wide DNAm comparison across the human brain, blood, saliva and buccal cells. They proposed that to determine the optimal surrogate tissue for representing brain FIGURE 2 | Flowchart of data selection using the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA).
DNAm, the DNAm patterns specific to the genomic region of interest between these two tissues must be considered (Braun et al., 2019). More importantly, the proxy tissue should have similar responsive changes to environmental influences such as drugs. Building a reliable peripheral-brain relationship will be critical in proving the biological relevance and illustrating the underlying biological mechanisms of available in vivo human tissue for clinical application.
Laboratory techniques have been developed for measuring DNAm not only at global methylation levels but also individual nucleotide levels. In this review, the research strategies used to assess drug response-related DNAm included target analyses for candidate genes, the discovery of genome-wide DNAm with microarrays, and assessment for global methylation. Microarray-based technologies such as Illumina Infinium HumanMethylation27 BeadChip R Array (HM27K) (Bibikova et al., 2009), Infinium R HumanMethylation450 BeadChips (HM450K) (Dedeurwaerder et al., 2011(Dedeurwaerder et al., , 2014 and Infinium R MethylationEPIC BeadChip (EPIC) (Moran et al., 2016), have been widely used for methylation profiling. Methylated DNA immunoprecipitation sequencing (MeDIP), Reduced Representation Bisulfite Sequencing (RRBS) and others are also used for DNAm measurement (Song et al., 2005;Couldrey and Cave, 2014). Whole-genome bisulfite sequencing is the gold standard for measuring CpG and non-CpG methylation, but due to its cost, few published pharmacoepigenetic studies have used it.

Research Strategies
The most widely used approach for studying drug responserelated DNAm is still a "candidate gene" method. We found 31 studies of 26 candidate genes completed in the three disorders of interest to us, and findings from 19 of them (73%) reached nominal significance (p < 0.05). These candidate gene studies only focus on a limited number of genes, based on their established association with a disorders' pathogenesis or hypothesized relation to response to drug treatment. Although FIGURE 3 | Current research types of pharmacoepigenetics in psychiatric disorders. The red number represents the number of "candidate gene studies," the green color means the number of "genome-wide studies," and the purple one means the number of "global methylation studies." previous candidate gene studies have identified target genes that may predict or reflect the clinical effect of drug response, most studies had limited statistical power, a concern especially for dense array data. Because of this limitation, the candidate gene method can provide little comprehensive understanding regarding the mechanism of action of drugs, but can only confirm previous knowledge or test relevance of a particular gene methylation (Harrison, 2015).
Along with next-generation sequencing, array-based analyses are commonly used in DNAm measurement to assess sitespecific methylation in the genome. These genome-wide analyses offer a non-biased experimental approach to identify novel candidates (Kurdyukov and Bullock, 2016), but we found only 3 studies that that implemented genome-wide analysis, using microarrays to evaluate the drug response. These studies identified several loci associated with clinical improvement, almost all of which were not previously predicted to be relevant. Unfortunately, genome-wide DNAm measurement faces two challenges, first its cost (added to costs of a clinical trial) and second, the challenge of having large samples in a study of patients receiving controlled treatment regimens. These problems limit the size of data sets, and therefore statistical power, for array analysis of drugrelated treatment effects. However, genome-wide approaches using a microarray do provide more comprehensive data and could identify novel and unexpected DNAm effects (Barros-Silva et al., 2018).
Additionally, three studies evaluated global DNAm in a specific tissue for a cross-medications comparison. Global DNAm reflects the DNAm status of total genomic content within a sample. Several methods exist to assess global DNAm status, including enzyme-linked immunosorbent assays (ELISA), the use of a previously validated protocol of Linear Interspersed Nuclear Element 1 (LINE-1) and the use of restriction enzymes (Burghardt et al., 2020). However, this approach only assesses DNAm changes at the genome level, and thus does not further provide data on region-specific DNAm changes to characterize critical regulatory regions (e.g., CpG islands and promoter regions).

PHARMACOEPIGENETIC FINDINGS IN SCZ, BD, AND MDD
The included pharmacoepigenetic studies were divided into three types: (1) medication-induced DNAm changes; (2) relationships between DNAm and clinical improvement, and (3) comparison of DNAm status across different medications (Figure 3). In this section, we summarize recent findings from pharmacoepigenetic studies in SCZ, BD and MDD to evaluate the reproducibility of drug response-related DNAm changes in these disorders.

Medication-Induced DNAm Changes
The changes of DNAm after a period of drug treatment were examined in studies comparing DNAm changes between preand post-treatment, which could reflect the effects of drug treatment directly. There were eight studies ( Table 1, SCZ = 2, BD = 1, and MDD = 5) that examined medication-induced DNAm changes, including seven that were candidate gene studies and one genome-wide study.
For candidate gene studies, there were 9 candidate genes included and only 5 reached nominal significance (p < 0.05). Of those five, only brain-derived neurotrophic factor (BDNF) was analyzed for MDD patients in 2 independent studies (Tadić et al., 2014;Wang et al., 2018b). A significant increase of DNAm in BDNF after 8 weeks of escitalopram treatment for 44 MDD patients based on blood samples was reported by Wang and collaborators (Wang et al., 2018b). The other study did not find significant DNAm changes in leukocytes after 6 weeks of monoaminergic drug treatment for 39 MDD patients (Tadić et al., 2014).
One explanation for this discrepancy is that escitalopram belongs to the SSRI class which has distinctly different mechanisms of action from monoaminergic drugs. The treatment duration and tissue types were also different between these two studies.
The genome-wide study identified 29,134 sites (p < 0.05) that were significantly differentially methylated after 1 year of treatment with clozapine for 21 SCZ patients based on blood samples, using HM450K assays . This relatively small study has not been replicated to date. We further evaluated whether the results of candidate gene studies for antipsychotics treatment could be reproduced in genome-wide studies. We found none of those significant candidate genes tested reached a nominal significance in genome-wide analysis.

The Relationship Between DNAm and Clinical Improvement
More studies explored the relationship between DNAm levels at baseline and clinical improvement after a period of drug treatment. The most frequently used method to assess the clinical improvement is still the scale score (e.g., PANSS, HAM-D). Patients are separated into remission, non-remission, response or non-response based on changes of scale scores. Remission, the symptom of virtual absence of disease, is the goal of treatment (Paykel, 1998). Response is defined as a clinically meaningful reduction in symptoms (e.g., a reduction of at least 50% in baseline symptom levels) (Aaronson et al., 2017). In this section, there were 24 studies included ( Table 2). Examining the correlation between baseline DNAm and treatment outcomes for patients, these studies address the question of whether pre-treatment DNAm features could predict clinical improvement.
Three studies used a genome-wide approach to evaluate the correlation between baseline DNAm and clinical improvement. Takeuchi et al. (2017) analyzed the changes of baseline DNAm between responders and non-responders of 20 paroxetine-treated MDD patients based on blood samples, using HM450K assays. Responders and non-responders had two CpG sites on genes PPFIA4 and HS3ST1 that were differently methylated (q < 0.05). Another study compared the DNAm changes in blood between responders and non-responders with 8 weeks of escitalopram treatment for 177 MDD patients, using the EPIC assays (Ju et al., 2019). They identified 303 (p < 0.05) sites with nominally significant differences between responders and non-responders, but none were significantly different after correction for multiple comparisons. Since these two drugs, paroxetine and escitalopram, are both SSRIs, we compared the reproducibility of the nominally significant results between these two studies. We found that none of the significant loci found to be associated with paroxetine could be replicated in the escitalopram treatment analysis. This failure at replication has to be considered in the context of differences in the platforms used, limitations in statistical power, and the 303 escitalopram-related sites detected by EPIC analysis were not included in HM450K assays. In evaluating whether results of candidate gene studies for MDD could be reproduced in genome-wide studies, we did not detect any overlap.
We did find one study involving schizophrenia, which evaluated the correlation between DNAm changes in blood after 1-year of clozapine treatment and changes of PANSS scores for 21 SCZ patients . They found that DNAm in a site located at the CREBBP gene was negatively correlated with clinical improvement after multiple testing correction. However, this finding remains to be replicated. Further analysis showed that none of these significant candidate genes tested for antipsychotics (e.g., HTR1A, CYP3A4, COMT), especially clozapine, reached a nominal significance in genome-wide analysis.

Comparison of DNAm Status Across Different Medications
Major medications for SCZ, BD, and MDD included antipsychotics, mood stabilizers and antidepressants. Antipsychotic medications are one of the primary treatments for the acute phases and the long-term prevention of recurrence for SCZ patients (Leucht et al., 2009;Hasan et al., 2017). Mood stabilizers, especially lithium, are a cornerstone in the long-term treatment of BD (Rybakowski, 2011;Geddes and Miklowitz, 2013). MDD patients are often treated with antidepressants (Drozda et al., 2014). However, the drug medications for each disorder are not mutually exclusive, as antipsychotic drugs may also be used in BD and MDD patients (Cruz et al., 2010;Patkar and Pae, 2013;Poo and Agius, 2014). The therapeutic efficacy of individual drugs is quite variable, and treatment trials to establish efficacy are laborious, so establishing biomarkers to improve guidance for treatment planning of individual patients is an important aim of psychiatric research. Predicting response to specific therapies could promote personalized medicine.
In this review, we evaluated 14 studies that compared DNAm after treatment with different medications, including 11 studies that used a candidate gene method and 3 that analyzed global methylation ( Table 3). We found that antidepressant medication increased the DNAm status of promoter regions of the BDNF gene in both BD and MDD patients beyond the level seen in antidepressant-free patients. With advantages and disadvantages, it is worth noting that drug selection and treatment duration varied across studies.

CHALLENGES AND FUTURE DIRECTIONS Underpowered Studies -Large Cohort Collaborations Are Needed
Insufficient statistical power is a major challenge for pharmacoepigenetic studies in psychiatric disorders. The cost of using methylation panels, the high cost of clinical trials where drug therapy is controlled along with diverse drugs used across studies, and the large number of sites for methylation investigation, together conspire to limit sample sizes and statistical power to detect effects of interest. These factors limit our ability to evaluate reproducibility of current pharmacoepigenetic findings. Although several interesting associations between DNAm and drug treatment outcomes occurred in relatively small sample sizes, findings from these studies should be considered preliminary until replication in a larger sample size is pursued. So, given that collecting clinical data with a large-sample is challenging and costly, small-sample studies remain the focus of most investigators in this area via hypothesis-generating studies.
To fill this gap, multi-site collaborations with ethnically diverse groups are needed in studies of the pharmacoepigenetics of psychiatric disorders. Prominent consortia have been established which looked at pharmacogenomics in psychiatric disorders, like the STAR * D project (Rush et al., 2004), the Genome-based Therapeutic Drugs for Depression (GENDEP) project (Uher et al., 2010), the Chinese Antipsychotics Pharmacogenomics Consortium (CAPC) (Yu et al., 2018), and others. These consortia benefit ongoing and future pharmacoepigenetic studies by providing preliminary findings that can be examined in new samples. For example, Powell et al. (2013) found that DNAm in IL11 could predict the clinical response to antidepressants using samples from the GENDEP project. Large international consortia, as have been developed for genetic association studies, are needed to study adequate numbers of patients with samples collected using a uniform procedure and with standardized drug treatment.

Recommendation for Clinical Design: Drug Selection, Treatment Duration, Evaluation, and Tissue Selection for Pharmacoepigenetic Studies
Neuropsychiatric drugs are grouped into various classes with slightly different mechanisms of action. It is clear that from an efficacy perspective, first line treatments for all these disorders leaves considerable room for improvement, potentially in new drug development and in strategies for personalized medicine (Leucht et al., 2012;Cipriani et al., 2018;Huhn et al., 2019;Pillinger et al., 2020). Identification of new biomarkers for these disorders is difficult, primarily because of the lack of knowledge about disease pathophysiology and mechanisms of drug action.
A major aim of psychopharmacology research is to develop valid strategies for selecting the right drug with the right dose for the "right" patients to optimize outcomes with low risks of side effects and relapse. Prior work is limited in supporting progress in this area. Across the published studies in this area, the types of medications were not uniform for each disorder, making comparison of studies difficult. Only 8% of studies used monotherapy for patients. While 27% of studies claimed to use monotherapy for patients, but we noticed that besides the same  primary medication, they temporarily permitted secondary medications (e.g., benzodiazepines, flunitrazepam) for insomnia or some side effects, which themselves can also exert epigenetic effects. Choice of duration for drug treatment trials is also important, yet variable across prior work. Treatment-resistance, nonadherence, late-emerging side effects and relapse are frequent events in clinical practice, some of which are seen only in longer term protocol participation. Aside from these factors, DNAm status is dynamic and can be affected by various environmental events during treatment that are unrelated to drug effects (Bjornsson et al., 2008).
Definition of treatment outcomes for psychiatric disorders is also a challenge, as different approaches may be used across studies. A recent study suggested that integrating genomics and phenotypic measures data could increase the accuracy of prediction of drug response (Kauppi et al., 2018). Thus, future research should consider an integrated examination of genetic, DNAm and dense phenotype assessment of drug effects not only with behavioral ratings but with direct in vivo assessment of brain physiology and treatment . Notably, all of these prediction models need to be covered across the relevant age span, ethnic groups, and sexes.
Predictors of drug response for psychiatric disorders ideally need to have biological relevance based on a solid mechanistic understanding of the pathophysiology of brain dysfunction. Most of the pharmacoepigenetic studies utilized clinically accessible tissues such as blood, saliva, and others. However, whether biomarkers identified in these tissues could reflect similar DNAm changes in the brain needs further study. Several brain banks with post-mortem tissue from patients that had neuropsychiatric disorders have been established, which generated big data sets at multiple regulatory levels (e.g., DNAm, gene expression, and other omics) . These multidimensional data sets could facilitate the exploration of the peripheral-brain relationship. Since epigenetic modifications are known to fluctuate over the lifespan, resolving the causal relationships in the epigenetic study of drug response using postmortem brains will be challenging. Cellular models will likely fill the gaps where both peripheral tissue and postmortem brain will not be able to deliver the answers.

System-Wide Analysis Could Benefit Pharmacoepigenetic Studies
Candidate gene studies contributed to answering specific questions in biology, but they do not provide system-wide information. In this review, we noticed that the assessed candidate genes were generally not the top signals in genomewide studies, and could not even reach a genome-wide significance. The genes that change the most between pre-and post-treatment, or are most strongly associated with clinical improvement, remain to be discovered. Also, pathway-based research rather than gene-based approaches may be more efficient. Different genes that converge on the same pathway may each contribute modestly, but together may robustly co-influence the functional abnormality. Therefore, systemwide research could provide better biomarkers than individual candidate genes.
Regulation of drug responses occurs at various levels including genetics, epigenetics, transcriptional, and protein modification, and also involves many functional pathways (Amare et al., 2017). To date, the underlying system of drug response-related DNAm remains poorly understood. Some have proposed that one or more components involved in a pathway, even if not risk genes, may be a better drug target than risk genes themselves (Wang et al., 2010;Harrison, 2015). One study assessed the interaction between genetic variates with DNAm for the IL-11 gene, which identified a possible regulatory relationship that could be used to predict antidepressant response for MDD patients (Powell et al., 2013). However, studies to date only scratched the surface of the complex regulatory system of drug response. Systematic integration of genomic data and genomewide DNAm data will be more effective and unbiased for learning about these regulatory processes, hopefully leading to the discovery of underlying drug response networks.

Data Analysis for Pharmacoepigenetic Studies
Strict statistical criteria and sufficient attention to covariances are critical for reducing false-positive rates. DNAm status is not a constant like DNA sequence are therefore is subject to influences from more confounders. Further, approaches for DNAm assessment can be easily confounded by technical artifacts. DNAm levels are also affected by factors including age (Horvath, 2013;Jaffe et al., 2016;Merid et al., 2020), sex (Maschietto et al., 2017;Xia et al., 2019), circadian rhythms (Lim et al., 2014), smoking (Shenker et al., 2013), drinking (Philibert et al., 2012), and others. All of these considerations for study design and data analysis need to meet standard requirements for drug response-related DNAm studies as in most genomic studies.
Methylation profiling could be easily confounded by batch effects. Systematic error can be introduced when samples are processed in multiple batches (e.g., the same sample measured at different times), which cannot be eliminated unless all samples are run in a single batch (Chen et al., 2011). It is especially important that cases and controls not be put in separate batches. Except for batch effects, positional effects also exist in the microarray and bias analysis . Positional effects are emerging when the same sample in different physical positions on the array and could bias methylation levels and lead to false findings. From our review, we observed that only one study did the batch effects correction (Ju et al., 2019), and none of these studies corrected for positional effects according to the method description in those papers. We cannot rule out the possibility that the data were properly processed but failed to be reported in the papers, but not reporting such details at least indicated the lack of attention to the serious issues.
Considering the influence of demographic information (e.g., age, sex, BMI) is also important. In this review, nearly 49% of studies did not control for covariates in their analyses. When there was a correction for covariates, typically only sex and age was controlled. Since DNAm is highly cell-typespecific, cellular heterogeneity of blood may skew DNAm patterns, influencing findings in drug response for psychiatric disorders (You et al., 2020). All studies did not control cell-type compositions. Collecting data about potentially useful covariates, and controlling for them in analyses, should be standard practice in future research for drug response studies. However, we note that large samples will be needed to develop appropriate modeling for covariate effects, as a few outliers or complex covariate interactions can exert effects that are challenging to deal with in smaller sample studies.
Additionally, strict and consistent statistical criteria for defining significance need to be applied that correct for multiple comparisons. Most published studies considered a nominal p < 0.05 to be the threshold for statistical significance, even when multiple genes or CpG sites were tested. We noticed that only 59% of the previous studies did the correction for multiple testing. False-positive results likely occurred in published work (Joober et al., 2012).

Summary
In this review, we found that only DNAm of BDNF was consistently related to clinical outcomes, an effect observed in four independent studies of MDD patients. Although pharmacoepigenetic studies to date are still limited in number and sample size, preliminary studies to date that evaluated changes in drug response-related DNAm across major psychiatric disorders provide evidence for an intriguing possibility of predictive DNAm biomarkers for drug treatment. A system-wide study using strict analytical and statistical procedures in large multisite studies may best advance identification of useful DNAm predictors for drug response in psychiatric disorders.

AUTHOR CONTRIBUTIONS
JZ, ML, XW, and YH did the manuscript collection. JZ did the manuscript writing. JZ and ML designed and generated the manuscript figures and tables. YX, JS, and RK helped with manuscript writing. CC and CL helped with study design and manuscript writing. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank Zhi Xu for her manuscript collection contribution, which greatly improved the review.