The Potential Role of Regulatory Genes (DNMT3A, HDAC5, and HDAC9) in Antipsychotic Treatment Response in South African Schizophrenia Patients

Despite advances in pharmacogenetics, the majority of heritability for treatment response cannot be explained by common variation, suggesting that factors such as epigenetics may play a key role. Regulatory genes, such as those involved in DNA methylation and transcriptional repression, are therefore excellent candidates for investigating antipsychotic treatment response. This study explored the differential expression of regulatory genes between patients with schizophrenia (chronic and antipsychotic-naïve first-episode patients) and healthy controls in order to identify candidate genes for association with antipsychotic treatment response. Seven candidate differentially expressed genes were identified, and four variants within these genes were found to be significantly associated with treatment response (DNMT3A rs2304429, HDAC5 rs11079983, and HDAC9 rs1178119 and rs11764843). Further analyses revealed that two of these variants (rs2304429 and rs11079983) are predicted to alter the expression of specific genes (DNMT3A, ASB16, and ASB16-AS1) in brain regions previously implicated in schizophrenia and treatment response. These results may aid in the development of biomarkers for antipsychotic treatment response, as well as novel drug targets.


INTRODUCTION
The onset of schizophrenia is marked by a first psychotic episode, typically followed by subsequent relapse episodes, separated by intervals of remission (Lieberman et al., 2001). Diagnosis remains difficult due to a heterogeneity of symptoms as well as symptom overlap with other disorders (Tandon et al., 2009). Furthermore, treatment strategies are not optimal (Brandl et al., 2014), and it is estimated that approximately half of all patients with schizophrenia will not respond satisfactorily to antipsychotics (Lohoff and Ferraro, 2010), which are the mainstay of treatment and, as such, widely used. They are effective for positive symptoms (such as delusions and hallucinations); however, their efficacy for negative symptoms (such as apathy, anhedonia, and social withdrawal) is limited (Leucht and Davis, 2017). Moreover, antipsychotics are known to result in a number of adverse drug reactions (ADRs), including motor abnormalities (Chowdhury et al., 2011) and metabolic deficits (Tandon et al., 2010). These ADRs are often severe and long lasting, resulting in reduced compliance and diminished positive outcomes (Brandl et al., 2014). Considering the high rate of nonresponders to treatment and the potential severe side effects of treatment, there is a clear need to improve our understanding of antipsychotic treatment response.
Pharmacogenetics, the study of the effects of genetic variation on treatment outcomes, has been moderately successful in explaining variability in inter-individual treatment response. Variation within the dopaminergic pathway has been extensively investigated, and several variants within dopamine receptor genes are associated with treatment response (Lencz et al., 2010) and ADRs (Bakker et al., 2008). In addition, variation within genes encoding drug-metabolizing enzymes has yielded similar findings of association (Lohoff and Ferraro, 2010;Brandl et al., 2014). Larger hypothesis-free-driven genome-wide association studies (GWAS) have also associated common variation with antipsychotic treatment response (Liou et al., 2012;Zhang and Malhotra, 2013); however, there has been little validation or replication of these associations, and their biological relevance remains to be determined (Liou et al., 2012;Zhang and Malhotra, 2013). In addition to these challenges, the majority of treatment response heritability is not explained by common variation, suggesting that other factors must also play a role (Manolio et al., 2009). This underscores the complexity and multi-factorial nature of treatment response since common and rare genetic factors, environmental factors, and gene-environment interactions need to be considered (Manolio et al., 2009;Majchrzak-Celińska and Baer-Dubowska, 2017).
Epigenetics refers to molecular mechanisms that determine inherited cellular phenotypes without alteration of the genotype (Wu et al., 2012). These mechanisms include various molecular processes, such as histone modification, nucleosome remodeling, non-coding RNAs, and DNA methylation (Wu et al., 2012). Both unique and overlapping altered epigenetic modifications are associated with schizophrenia etiology and pathogenesis as well as antipsychotic treatment (Swathy and Banerjee, 2017). As such, an extremely complex, multi-directional relationship needs to be considered in studying the role of epigenetics in treatment response (pharmacoepigenetics), since genes that are implicated may be regulated by epigenetic modification independent of disease etiology, treatment outcomes, and/or ADRs (Kurita et al., 2012;Melas et al., 2012;Tang et al., 2014;Majchrzak-Celińska and Baer-Dubowska, 2017). For example, alterations of DNA methylation profiles (Melas et al., 2012;Tang et al., 2014) and altered chromatin structure (Kurita et al., 2012) may influence treatment response.
Although specific gene-environment interactions are required for epigenetic modifications to occur, it is important to note that specific genes directly or indirectly produce proteins and other products necessary for these modifications. For example, the DNMT1 gene encodes for Dnmt1, which is responsible for the maintenance of existing methylation patterns during cell division (Bostick et al., 2007), while DNMT3A and DNMT3B encode for the enzymes responsible for establishing de novo methylation patterns (Okano et al., 1999). Regulatory genes, such as those involved in DNA methylation and transcriptional repression, are therefore excellent candidates for investigation with regards to antipsychotic treatment response.
The aims of this study were, therefore, to identify candidate regulatory genes that may be involved in antipsychotic treatment response, and then to determine if variation within these genes is associated with treatment outcome. In order to address these aims, we first sought to identify regulatory genes that were differentially expressed between patients with schizophrenia (including chronically medicated patients and drug-naïve first-episode patients) and healthy controls. By identifying genes differentially expressed between the chronically medicated patients and both the drug-naïve first-episode patients and healthy controls, we identified a subset of candidate genes that may be implicated in treatment response independent of schizophrenia pathogenesis. Variation within these candidate genes was then investigated in an independent cohort for association with treatment response over time.

MATERIALS AND METHODS
An outline of the approach used to perform this study is provided in Figure 1.

Participants for Gene Expression Analyses
The cohort consisted of 20 unrelated, age-matched, and male South African participants of South African colored (SAC) descent [10 patients with schizophrenia (SCZ), 26.4 ± 7.9 years old, and 10 healthy controls (CON), 26.5 ± 7.6 years old]. Patients with SCZ were further divided into two equal groups consisting of five first-episode patients (FES) (24.8 ± 10.9 years old) and five chronic patients (CHR) (28.0 ± 4.1 years old). The FES patients were recruited and sampled within 1 week of their first episode of psychosis and were subsequently administered flupenthixol decanoate (Fluanxol, Lundbeck, Copenhagen, Denmark). The CHR patients were recruited 6.2 ± 0.4 years after their first episode, and all were treated with flupenthixol decanoate (Fluanxol, Lundbeck, Copenhagen, Denmark)

Total RNA Isolation and cDNA Synthesis
Whole blood was collected from all participants by venipuncture of a forearm vein into PAXgene Blood RNA tubes (Qiagen, California, USA), which were stored at −20°C until processed.
Total RNA was extracted using the PAXgene Blood RNA Kit IVD according to the manufacturer's instructions (Qiagen, California, USA). All samples were eluted in 80 µl of elution buffer and stored at −80°C until further analysis. RNA yield and quality were assessed using an Agilent Model 2100 Bioanalyzer (Agilent Technologies, California, USA) and a DropSense 16 spectrophotometer (TRINEAN, Belgium). All samples had 260/280 > 2.0 and RIN > 7.0.
Reverse transcription was performed using the High-Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (Applied Biosystems, California, USA), according to the manufacturer's specifications. Briefly, for each sample, 100 ng of RNA was added to 2 μl of random and oligo (dT) primers in a final reaction of 20 μl. These reaction tubes were then placed in the GeneAmp ® PCR Systems 2700 (Applied Biosystems, California, USA) thermocycler at 25°C for 10 min, 37°C for 120 min, and 85°C for 5 min to inactivate the reverse transcriptase enzyme. The cDNA samples were then stored at −20°C until analyzed.

Quantitative Real-Time PCR
Quantitative real-time PCR (qRT-PCR), to determine the relative mRNA abundance, was performed using the StepOnePlus Real-Time PCR System and SDS Software version 2.3 (Applied Biosystems, California, USA). The commercially available fluorescence-based TaqMan ® Human DNA Methylation and Transcriptional Repression microarray plates (Applied Biosystems, California, USA) were used to assess the relative mRNA content of 27 regulatory genes ( Table 1). Each sample was analyzed in two independent experiments to control for experimental bias using the following PCR conditions: a 10-min heat activation step (95°C) followed by 40 cycles of 15 s at 94°C and 1 min at 63°C. Fluorescence data, indicative of the amount of PCR product, was captured at each cycle. The relative mRNA concentrations were then calculated based on the cycle number that the threshold quantity of PCR product is obtained (Ct). The RefFinder online tool was used to assess the stability of potential housekeeping genes (https://www.heartcure.com. au/reffinder) (Xie et al., 2012), and GAPDH was identified as the most stable. Ct values were therefore normalized to values of GAPDH and expressed relative to this control (Livak and Schmittgen, 2001).

Gene Expression Analyses
Differential gene expression between the CON and SCZ groups was determined using unpaired t-tests or Mann-Whitney U tests where appropriate. Additionally, one-way analysis of variance (ANOVA) and Tukey's post hoc tests were used to determine differences between the CON, FES, and CHR groups. The false discovery rate (FDR) correction (Benjamini and Hochberg, 1995) was used to correct for multiple testing in the CON vs. SCZ and CON vs. FES vs. CHR analyses (27 genes; FDR < 0.01). A significance threshold of p < 0.05 was used for the Tukey's post hoc tests since these were only performed in the case of a significant (FDR < 0.01) ANOVA result. Candidate genes for genetic association analyses with treatment response were selected if they met the following three criteria: i) significant FIGURE 1 | A flow diagram outlining the approach used in this study. (A) The first arm of the study involves identifying differentially expressed genes between healthy controls and schizophrenia patients, further divided into drug-naïve first-episode patients and chronic patients. (B) Differentially expressed genes were then used as candidates for genetic association with antipsychotic treatment response in a drug-naïve FES cohort. FES, first episode schizophrenia; CHR, chronic schizophrenia; SCZ, schizophrenia; CON, healthy control; eQTLs, expression quantitative trait loci. differences in gene expression between the CHR and CON groups, ii) no significant differences in gene expression between the FES and CON groups, and iii) significant differences in gene expression between the CHR and FES groups.

Participants for Genetic Association Analyses
The patient cohort included 103 unrelated South African FES patients meeting DSM-IV TR (American Psychiatric Association, 2000) diagnostic criteria for SCZ, schizophreniform disorder, or schizo-affective disorder (80% SAC, 12% Xhosa, and 8% European descent) and have been described previously (Drogemöller et al., 2014;Chiliza et al., 2015;Ovenden et al., 2017;O'Connell et al., 2018). All patients received treatment with flupenthixol decanoate (Fluanxol, Lundbeck, Copenhagen, Denmark), a long-acting injectable antipsychotic, according to a fixed protocol. Treatment response was assessed using the Positive and Negative Syndrome Scale (PANSS) (Kay et al., 1987) over a period of 12 months, with measurements taken biweekly for the first 6 weeks, and every 3 months thereafter. Written informed consent was obtained from all patients, or their caregivers, prior to the study, and ethical approval was granted by HREC, Faculty of Medicine and Health Sciences, Stellenbosch University (N06/08/148).

Genetic Association Analyses
Variants within these candidate genes, obtained from each National Center for Biotechnology Information gene page (https://www.ncbi.nlm.nih.gov/gene/), were mined from available genome-wide genotype data. All 103 FES patients described above were previously genotyped with the Infinium OmniExpressExome-8 Kit (Illumina, California, USA) in accordance with the standard Illumina protocol. Variants were excluded from downstream analysis if they exhibited high rates of missing genotype data (> 5%), if their minor allele frequency (MAF) was < 1% or if they showed departure from the Hardy Weinberg equilibrium (p < 1×10 -4 ). For this study, only variants with a minor allele frequency greater than 5% and a call rate greater than 98% were included for the association analyses. Furthermore, only one representative variant was included when two or more variants were shown to be in linkage disequilibrium (LD, r 2 > 0.6) with one another. In total, 274 variants were included for analyses. All processing of genetic data was performed using Plink v1.9 (Purcell et al., 2007;Purcell and Chang).
Association analyses were conducted in R (R Core Team, 2017) using R packages lme4 (Bates et al., 2014) and lmerTest (Kuznetsova et al., 2016). Linear mixed-effects models were used to investigate the effect of genetic variants on change in PANSS scores for each subscale (positive, negative, and general) and total over the 12-month period, adjusting for age, gender, proportion ancestry, and baseline PANSS scores. Multiple modes of inheritance were investigated and the Bonferroni correction method was used to correct for multiple testing (274 variants, four modes of inheritance, four PANSS domains; threshold p = 1.14 × 10 -5 ).

Relative Gene Expression Results
After the initial analysis, comparing gene expression between the CON and SCZ groups, 14 genes were shown to have significantly different expression between these groups. Specifically, the CH4, DNMT1, DNMT3A, DNMT3B, HDAC3, HDAC5, HDAC6, HDAC7, HDAC9, HDAC11, MBD3, RBBP7, SAP30, and SIN3A genes showed increased expression in the SCZ group when compared to the CON group (Table 1). When comparing the CON group to the FES and CHR groups, 19 genes were shown to be differentially expressed ( Table 1). Post hoc analyses identified seven genes that were significantly over-expressed in the CHR group when compared to both the CON and FES groups.
For follow-up investigation as candidates for association with antipsychotic treatment response, seven genes were selected. Specifically, the B2M, DNMT3A, HDAC5, HDAC9, MBD2, MBD3, and RPLP0 genes were selected since the expression levels of these genes were significantly different between the CHR and CON groups and CHR and FES groups, respectively (Table 1, Figure 2). Furthermore, no significant differences in gene expression were identified for these genes when comparing the FES and CON groups (Table 1, Figure 2).

Genetic Association Results
Of the 274 variants investigated within the seven candidate genes, four were significantly associated with treatment response, as described by a change in PANSS scores over time, when considering correction for multiple testing (p < 1 × 10 -5 ) ( Table 2). All significant associations were identified with the PANSS-negative domain. Specifically, the DNMT3A rs2304429 CC genotype was significantly associated with an improved treatment response (greater reduction in PANSS-negative scores per month) when compared to the TC genotype. The rate of PANSS-N score reduction was significantly faster over time (an additional 1.92% per month) in patients with the DNMT3A rs2304429 CC genotype when compared to patients with the TC genotype. The HDAC5 rs11079983 TT genotype was significantly associated with a poorer treatment response (less reduction in PANSS-negative scores per month) when compared to the CC genotype. Patients with the HDAC5 rs11079983 TT genotype had significantly slower rate of PANSS-N score reduction over time (less by 2.36% per month) when compared to patients with the CC genotype. Two variants within HDAC9 (rs1178119 and rs11764843) were also significantly associated with poorer treatment response (less improved PANSS-negative treatment trajectory scores) as shown in Table 2. For HDAC9 rs1178119, patients with the GA genotype had a significantly slower rate of reduction in PANSS-N scores over time (less by 2.04% per month) than patients with the AA genotype. Presence of the HDAC9 rs1178119 G allele was also significantly associated with a slower rate of reduction in PANSS-N scores over time (less by 1.52% per month per G allele). Similarly, a significantly slower rate of reduction in PANSS-N scores over time (less by 1.92% per month) was also identified for patients with the HDAC9 rs11764843 CA genotype when compared to those with the AA genotype.

Bioinformatics
The variants significantly associated with change in PANSS scores were assessed using the BRAINEAC (Ramasamy et al., 2014) and GTex (GTEx Consortium, 2013, GTEx Consortium, 2015 databases. Two variants (rs2304429 and rs11079983) were identified as potential brain-specific eQTLs. The rs2304429 variant was suggested by BRAINEAC to alter the expression of DNMT3A in particular brain regions-namely, the putamen (PUTM), the cerebellar cortex (CRBL), the temporal cortex (TCTX), and the medulla (MEDU). Specifically, the TT genotype was associated with reduced expression of DNMT3A in these brain regions. In addition, the HDAC5 rs11079983 variant was shown to alter the expression of ASB16 and ASB16-AS1 in the cerebellum when considering the GTex database. Specifically, the TT genotype is associated with reduced expression of ASB16 and increased expression of ASB16-AS1 when compared to the CC genotype, respectively. None of the other variants were identified as brain-specific eQTLs in the BRAINEAC or GTex databases.

DISCUSSION
We identified a number of DNA methylation and transcriptional repression genes to be significantly over-expressed in patients with SCZ, in first-episode and CHR, when compared to healthy controls (Table 1). Specifically, seven candidate regulatory genes were identified, including B2M, DNMT3A, HDAC5, HDAC9, MBD2, MBD3, and RPLP0, and variation within these genes was assessed for association with antipsychotic treatment response. Four variants were found to be significantly associated with poorer treatment trajectory in the PANSS-negative domain.
In this study, significant increases in the expression of DNMT1 and DNMT3A were identified between all patients with SCZ and controls; however, further analyses revealed that this increase was only present when considering the CHR patients and not in the FES patients. Increased DNMT1 and DNMT3A gene expression in the GABAergic neurons of SCZ patients has been previously identified (Zhubi et al., 2009), while the increased expression of DNMT1 was also established in the peripheral blood lymphocytes of patients with SCZ (Auta et al., 2013). When considering these previous studies, it is interesting to note that the mean ages of their study cohorts were 57 ± 11 and 56 ± 18 (Zhubi et al., 2009) and43.6 ± 10.3 (Auta et al., 2013), respectively. Given the young age of onset of SCZ, it is likely that the patients in these cohorts are more indicative of chronic SCZ and that the results of this study therefore replicate these previous findings.
The expression of a number of HDAC genes (1-4, 6, and 9) was previously investigated in the prefrontal cortex of patients FIGURE 2 | The relative expression levels of B2M, DNMT3A,HDAC5,HDAC9,MBD2,MBD3,and RPLP0 in the SCZ, FES, and CHR groups when compared to the CON group. ***Significant after correction, versus CON group. # Significant after correction, versus FES group. Data is presented as mean ± standard error. Additional details available in Table 1.
with SCZ (Sharma et al., 2008). Only HDAC1 was shown to have increased expression in patients with SCZ when compared to controls, while no significant differences in expression were identified for the other HDAC genes (Sharma et al., 2008). One possible explanation for the differences between these results and those of our study is that these patients were subject to a range of different medications, including typical and atypical antipsychotics, mood stabilizers (including valproic acid), antidepressants, stimulants, and sedatives which all may have an effect on the gene expression observed (Majchrzak-Celińska and Baer-Dubowska, 2017). Of the other differentially expressed genes in this study, SAP30 expression was previously investigated and no significant changes were identified (Vawter et al., 2006). These results highlight the need for well-defined and deep phenotyping when investigating the molecular etiologies of neuropsychiatric disorders since their molecular architecture is malleable to disorder progression as well as treatment and other environmental factors (Gurwitz and Pirmohamed, 2010). Due to the nature of the genes investigated in this study, changes in gene expression that were identified are indicative of altered regulatory mechanisms in chronically medicated patients with SCZ ( Table 1). The alteration of these regulatory mechanisms is likely the result of a combination of disease progression, antipsychotic medication (Swathy and Banerjee, 2017), and the influence of other environmental factors through the process of epigenetics (Feil and Fraga, 2012). To further elucidate these complex interactions, gene expression differences between first-episode and CHR, as well as healthy controls, were assessed. Significantly increased expression of seven genes (B2M, DNMT3A, HDAC5, HDAC9, MBD2, MBD3, and RPLP0) was observed in patients with chronic SCZ, when compared to FES and healthy controls (Figure 2). These seven genes were selected as candidates for association with antipsychotic treatment response.
Novel associations were identified between variants within DNMT3A (rs2304429), HDAC5 (rs11079983), and HDAC9 (rs1178119, rs11764843) and antipsychotic treatment response as defined by a change in PANSS scores over time (Table 2). Specifically, all four of these variants were associated with a significantly worse treatment trajectory in the negative PANSS symptom domain. Bioinformatics analyses of these variants revealed that two variants are predicted to exert functional changes as eQTLs. The DNMT3A rs2304429 and HDAC5 rs11079983 variants were predicted to alter expression of specific genes in particular brain regions. Specifically, individuals with the rs2304429 CC genotype have increased DNMT3A gene expression in the PUTM, CRBL, TCTX, and MEDU when compared to individuals with the TC or TT genotypes. In addition, individuals with the HDAC5 rs11079983 TT genotype have reduced expression of ASB16 and increased expression of ASB16-AS1 in the cerebellum when compared to the CC genotype. In this study, the rs2304429 CC and rs11079983 CC genotypes infer improved treatment response indicating that the brain region-specific gene expression changes associated with these variants may play a role. From these results, it may be hypothesized that increased expression of DNMT3A, as result of the rs2304429 CC genotype, in the abovementioned brain regions may result in de novo methylation patterns (Okano et al., 1999) that result in increased efficacy of antipsychotic medication. Similarly, reduced expression of ASB16 and ASB16-AS1 in the cerebellum, in the presence of the rs11079983 CC genotype, may result in altered ubiquitin-mediated pathways (Kohroki et al., 2005) and cytokine signaling (Babon et al., 2009) with beneficial effects when considering antipsychotic response. These brain regions have previously been implicated in SCZ (Hokama et al., 1995;Turetsky et al., 1995;Harrison, 2004;Yeganeh-Doost et al., 2011;Williams et al., 2014;Tohid et al., 2015) and suggested as targets for treatment (Buchsbaum et al., 2003;Hugdahl et al., 2009;Mitelman et al., 2009;Parker et al., 2014). These brain region-specific eQTLs should therefore be investigated as biomarkers and potential targets for antipsychotic treatment response.
The potential underlying mechanisms of action associated with the remaining two significant findings are not clear. Further studies investigating the exact role that these variants may have in antipsychotic treatment response is warranted. Furthermore, the novel associations identified in this study should be replicated and validated in additional independent cohorts. Moreover, the eQTL data presented were not generated from Sub-Saharan African-ancestry-related individuals and further studies are required to confirm whether the variants presented here have similar eQTL effects in populations of African-ancestry. Due to the small sample sizes used for the gene expression analyses in this study, these results should be interpreted with caution and require replication. Functional studies incorporating in situ and in vivo assays should also be considered for validation of these results. Moreover, the gene expression results presented in this study should be replicated in female patients with SCZ and controls.
In conclusion, this study identified significant differential expression of DNA methylation and transcriptional repression genes between FES with SCZ, chronically medicated patients with SCZ and healthy controls. Variants within specific differentially expressed genes were significantly associated with antipsychotic treatment response, and highlighted particular brain regions in which altered expression of specific genes may play a role in treatment outcome. These results may aid in the development of biomarkers for antipsychotic treatment response, as well as, novel drug targets, and treatment strategies.

ETHICS STATEMENT
Written and informed consent was obtained from all patients, or their caregivers, prior to the study and ethical approval was granted by the Human Research and Ethics Committee (HREC), Faculty of Medicine and Health Sciences, Stellenbosch University (N13/08/115 and N06/08/148).

AUTHOR CONTRIBUTIONS
KO'C NM, and LW conceived the study. KO'C performed the laboratory work and statistical analyses. RE and SS were involved