Integrative genomic strategies applied to a lymphoblast cell line model reveal specific transcriptomic signatures associated with clozapine response

Clozapine is an important antipsychotic drug. However, its use is often accompanied by metabolic adverse effects and, in rare instances, agranulocytosis. The molecular mechanisms underlying these adverse events are unclear. To gain more insights into the response to clozapine at the molecular level, we exposed lymphoblastoid cell lines (LCLs) to increasing concentrations of clozapine and measured genome-wide gene expression and DNA methylation profiles. We observed robust and significant changes in gene expression levels due to clozapine (n = 463 genes at FDR < 0.05) affecting cholesterol and cell cycle pathways. At the level of DNA methylation, we find significant changes upstream of the LDL receptor, in addition to global enrichments of regulatory, immune and developmental pathways. By integrating these data with human tissue gene expression levels obtained from the Genotype-Tissue Expression project (GTEx), we identified specific tissues, including liver and several tissues involved in immune, endocrine and metabolic functions, that clozapine treatment may disproportionately affect. Notably, differentially expressed genes were not enriched for genome-wide disease risk of schizophrenia or for known psychotropic drug targets. However, we did observe a nominally significant association of genetic signals related to total cholesterol and low-density lipoprotein levels. Together, these results shed light on the biological mechanisms through which clozapine functions. The observed associations with cholesterol pathways, its genetic architecture and specific tissue effects may be indicative of the metabolic adverse effects observed in clozapine users. LCLs may thus serve as a useful tool to study these molecular mechanisms further.


Introduction 54
Antipsychotic drugs (APs) play an important role in the treatment of psychotic disorders such as 55 schizophrenia (SCZ). Clozapine is one of the most effective antipsychotic drugs (AP) (Kane et al., 56 1988; Leucht et al., 2013;Taylor, 2017). However, the decision to prescribe clozapine is 57 complicated by its potential to induce severe adverse effects (Leucht et al., 2013). The most 58 severe adverse effect, with a prevalence of <1%, is clozapine-induced agranulocytosis, a 59 dramatic reduction of white blood cells (Andersohn et al., 2007). More common adverse effects 60 include weight gain, dyslipidemia and type 2 diabetes. These adverse metabolic effects are, in 61 addition to the chance of developing agranulocytosis, the primary reasons for patient 62 noncompliance and discontinuation of treatment (Cohen, 2014;Weiden et al., 2004). Here, we exposed LCLs to increasing doses of clozapine and collected both expression and 97 methylation profiles. Through integrative genomic analyses, we aim to extrapolate in vitro 98 molecular signatures of clozapine towards relevance of in vivo function and study clozapine 99 response without the need to assemble a large cohort of patients. We identified significant 100 changes in transcriptomic signatures associated with cholesterol and cell cycle pathways. By then 101 integrating these molecular profiles with genome-wide association study (GWAS) summary 102 statistics of different traits and diseases, we show that clozapine-associated genes also overlap 103 with genetic signals related to total cholesterol and low-density lipoprotein (LDL) levels but not 104 We processed raw gene expression values using the "Lumi" R-package (Du et al., 2008). We log2 137 transformed and quantile normalized the raw data, keeping only expressed gene transcripts 138 (detection p < 0.01) for further analysis (22,926 probes). We processed DNA methylation values 139 using the "WateRmelon" Bioconductor package (Pidsley et al., 2013), removing probes that were 140 known to cross-hybridize, probes containing SNPs in target CpG regions, probes with detection 141 p-value greater than 0.01 in 5% of samples, and probes with beadcounts > 3 (n=86,068 probes 142 in total) Price et al., 2013). We normalized the data using the dasen function 143 and computed β-values, defined as the ratio of the methylated probe intensity and the overall 144 intensity (sum of methylated and unmethylated probe intensities), to measure methylation levels. 145 To limit the effect of heteroskedasticity, we included only variable probes with β-values between 146 0.2 -0.8 in our analyses (165,014 probes) (Du et al., 2010). 147

Gene ontology analysis 167
We performed functional gene ontology analysis using DAVID (Database for Annotation, 168

Weighted gene co-expression network analysis 180
We performed a gene expression network analysis using weighted gene co-expression network 181 analysis (WGCNA) in R. Briefly, WCGNA identifies distinct modules using the shared variation in 182 gene expression based on pairwise correlation. To account for the biases related to differing 183 probe numbers between genes assayed on the array, we provided as input the mean probe 184 expression of genes residing within nominally significant differentially expressed genes (p<0.05), 185 considering 5,708 probes within 4,897 genes (Horvath, 2011;Langfelder and Horvath, 2008;186 Zhang and Horvath, 2005). To assign biological function to each WCGNA module, we performed 187 gene ontology analysis using DAVID we performed. 188 Additional DNA methylation analyses 189 We ran a candidate gene study for CpG sites in close proximity of genes with evidence of 190 clozapine-induced differential gene expression. Methylation probes within the gene body, in the 191 3' and 5' untranslated regions and up to 1,500 nucleotides upstream of the transcription start site 192 of the 463 'top genes' (p<0.05) were selected for a post-hoc analysis (n = 1,004). These results 193 can be found in Table 1  (1 + RPKM value). We then tested 1) if clozapine-associated genes have higher or lower average 200 expression within each tissue and 2) if between tissue distance is different for clozapine-201 associated genes compared to the expected based on chance. To calculate between tissue 202 distance, we used the top half most variable genes across GTEx samples that were also 203 significantly detected in our in vitro assay (n=7,025). We then performed multidimensional scaling

Results 255
Clozapine exposure induces widespread gene expression changes 256 We exposed lymphoblast cell lines to different concentrations of clozapine and tested for 257 subsequent changes in gene expression levels (Supplemental Figure 1A). In total, we tested 258 22,926 gene expression probes, of which 5,708 showed nominal significance (p<0.05) and 518 259 probes exceeded a Bonferroni-corrected p < 2.18 * 10 -6 . These 518 probes consisted of 234 up-260 regulated probes and 284 down-regulated probes, representing a total of 463 unique genes. The 261 top 10 up-regulated and down-regulated probes are shown in Table 1 Gene ontology enrichment analysis of up-regulated transcripts showed significant functional 264 enrichment for cholesterol metabolism (13 genes, p = 4.15 * 10 -15 ) and steroid biosynthesis (11 265 genes, p = 1.01 * 10 -8 ) ( Table 2), while analysis of down-regulated transcripts were enriched for cell 266 division processes and related annotations, such as mitosis (44 genes, p = 1.87 * 10 -39 ), 267 chromosome (49 genes p = 3.03 * 10 -35 ) and nucleosome (16 genes, p = 3.12 * 10 -14 ) and other cell 268 cycle pathways (Table 2) genes), which were nominally significant. The M14 co-expression module was enriched for genes 275 involved in cholesterol metabolism (13 genes, p = 4.8 * 10 -15 ), the M1 co-expression module was 276 enriched for genes involved in cell cycle (133 genes, p = 7.3 * 10 -32 ) and the M9 and M3 co-277 expression module was enriched for mitochondrial genes (15 genes, p = 6.7 * 10 -6 and 46 genes p 278 = 3.8 * 10 -7 respectively). The M10 co-expression module was enriched for genes involved in the 279 nucleosome (9 genes, p = 1.4 * 10 -7 ). 280 Minimal changes in DNA methylation at single CpG sites after clozapine exposure 281 We then performed DNA methylation profiling to assess whether these effects were due to 282 epigenetic changes. For statistical analysis, we applied an analytical approach similar to our gene 283 expression analyses. Targeted analysis on the 1,004 CpG sites near the 463 differentially-284 expressed genes revealed 3 probes exhibiting significant changes in DNA methylation (p < 285 1.08 * 10 -4 ), including a probe upstream of the low-density lipoprotein receptor (LDL-R) gene 286 (cg22971501, p=4.75*10 -5 ) after 24h; one probe upstream of the cyclin F (CCNF) gene showed 287 a significant change in DNA methylation after 96h (Table S3) Table 2. Main functional annotations associated with clozapine exposure based on gene expression. Output of pathway enrichment analyses using DAVID is shown. The enrichment score is used as the main metric of importance. It is defined as the geometric mean of all enrichment pvalues of each annotation term within the group. It is expressed as the minus log of the p-value, an enrichment score of 1.3 is nominally significant. This table shows clusters with an enrichment score > 4.6, corresponding to a p-value < 0.01. Fold enrichment is a measure to express the enrichment of this particular group of genes in comparison with the genes in the human genome. A list of all pathways is available in supplemental information.

DNA methylation is affected at the pathway level 293
To investigate if our top associated DNAm probes aggregated to changes at the pathway level, 294 we performed pathway enrichment analysis using GREAT, which incorporates functional 295 annotation from various databases to predict cis-regulatory function of genomic regions of 296 interest. When considering probes exhibiting FDR < 10%, we observed enrichment for protein 297 binding and regulation of cellular processes after 24h of clozapine treatment. When considering 298 the top 1000 probes, we also observed enrichment of immune-related functions, such as the 299 Major Histocompatibility Complex (MHC) class II protein complex and Graft-versus-host disease 300 after 24h. We found significant enrichment for estradiol regulation and various embryonic 301 developmental processes (Figure 1 and Table S5)

Clozapine transcriptomic profiles highlight multi-tissue effects in GTEx 330
We then asked whether our in vitro derived transcriptomic signatures could be used to help 331 translate the function of clozapine in humans. For this purpose, we used gene expression data 332 from the GTEx Project, including LCLs (n = 22 GTEx tissues, Figure 3A). First, we overlapped 333 preferentially-expressed genes of each GTEx tissue, as previously reported (Melé et al., 2015), 334 with the clozapine-associated genes detected in our assay. Preferentially expressed genes in 335 GTEx-LCLs exhibited the most overlap with the genes identified by our experiment followed by 336 whole blood (Supplementary Figure 2). We then investigated if the mean tissue expression of 337 clozapine-associated genes was significantly different from the mean expression of all the genes 338 tested in our experiment. We found that genes downregulated by clozapine have lower expression 339 in all tissues except testis and LCLs ( Figure 3B and Table S6). Downregulated genes are enriched 340 for cell cycle processes and their higher expression in the testis and LCLs likely represent the 341 proliferative nature of these tissues. Genes upregulated by clozapine have significantly higher 342 average expression in liver, muscle, lung, and fibroblasts, and lower average expression in testis 343 tissue. We then asked whether clozapine-associated genes have different between-tissue 344 distances as a proxy to investigate possible functional links with other tissues. Distance is 345 calculated using the Euclidean distance measure. A lower value indicates that two tissues are 346 more similar while a higher distance value indicates that these tissues are more dissimilar. Of the 347 406 tissue pairs, we found that genes upregulated by clozapine have significantly deviating 348 between-tissue distances across 31 pairs of tissues ( Figure 3C and Table S7). Spleen tissue 349 stands out with having the most different distance with other tissues (N=20). We also found 350 multiple differences for adipose, lung, and breast tissue pairs. The distance between cervix uteri 351 and ovary tissue, adrenal gland and thyroid tissue, muscle and nerve tissue are significantly more 352 dissimilar as well for upregulated clozapine genes compared to all genes detected in our assay. 353 For downregulated clozapine genes, we detected significantly dissimilar distances for 27 tissue-354 pairs, all involving testis tissue ( Figure 3C).  (n=7,025), that were also significantly detected in our experiment, across a random subset of 2,000 GTEx samples were used as input in the MDS analysis. (B) A forest plot visualizing mean tissue gene expression of clozapine-associated genes across several tissues. Clozapine genes were divided into groups that are up regulated and down regulated after drug exposure. The expected mean gene expression, based on 10,000 weighted samplings, is shown as well (dotted vertical line). Within each tissue, the observed mean tissue gene expression (x-axis) is normalized by subtracting the expected mean gene expression. P-values indicate whether the mean gene expression of clozapine up regulated genes significantly deviate from the expected mean expression within a tissue. P-values are corrected for the number of tissues tested. Downregulated genes were not tested. See Supplementary Tables for results of all tissues (C) Between tissue distance across all GTEx tissue pairs. Each point represents one tissue pair. The y-axis shows the Euclidean distance between tissues computed using only clozapine-associated genes. The x-axis shows the expected mean Euclidean distance across 50,000 weighted samplings. Tissue pairs for which the between tissue distance, based on clozapine genes only, significantly deviates from mean expected distance (Padjusted <0.05), are color-coded. P-values are adjusted for multiple testing by Bonferroni correction (n test = 465 pairs x tests = 930).

Clozapine transcriptome signatures are not enriched for schizophrenia disease risk 356
Previous studies have found associations between schizophrenia genetic susceptibility and 357 antipsychotic drug targets (Gaspar and Breen, 2017;Skene et al., 2018). To investigate whether 358 clozapine-induced in vitro gene expression signatures are also associated with genetic 359 susceptibility of schizophrenia, we conducted gene-set enrichment analysis using all differentially 360 expressed genes (n=463). We did not observe a significant enrichment (p = 0.91). We then 361 considered upregulated genes, and did not observe significant enrichment (p = 0.74), nor for 362 downregulated genes (p = 0.64). Gene sets defined according to <1% and <5% FDR did not 363 change these findings (Table S8). 364

365
To further understand these findings, we investigated whether SCZ genetic susceptibility 366 aggregates to antipsychotic drug target genes with and without clozapine targets. As SCZ genetic 367 risk has been associated with antipsychotic drug target genes, it could be that this signal is 368 primarily driven by non-clozapine genes. To examine this, we used drug target gene lists from 369 drug-gene interaction databases as previously reported (Roth et al., 2000;Wagner et al., 2016). 370 We were able to extract 53 drug target gene-sets belonging to the N05A class of drugs with 371 antipsychotic actions, including clozapine. Across drug target gene-sets, we mapped all targets 372 to 104 unique genes of which 41 were detectable in our gene expression data but none 373 overlapped with differentially expressed genes identified. We found no evidence for SCZ risk to 374 be enriched in antipsychotic drug target genes overall (n = 96 genes, p = 0.96) nor with clozapine 375 targets excluded (n = 52, genes, p = 0.60). We did observe a strong concordance between the p-376 values of individual drug gene-sets tested in our analysis and the p-values reported by the 377 previous study (Gaspar and Breen, 2017) (n = 50 drug gene sets, rho = 0.85, p = 1.49 * 10 -15 ), 378 indicating our analysis framework was able to reproduce previous findings at the level of individual 379 drug target gene-sets. Our analysis however does not observe an association between 380 schizophrenia genetic risk and clozapine-associated genes nor antipsychotic drug target genes. 381

Genes upregulated after clozapine exposure show an association with GWAS signal of total 382
cholesterol and low-density lipoprotein. 383 Next, we explored whether differentially expressed genes were enriched for genetic signal of 384 cholesterol-and cardiovascular disease-related traits. Using summary statistics of large GWASs 385 for each trait, we observe significant SNP-h 2 and a rich correlation structure between these traits 386 ( Figure 4A and 4B). We then integrated the observed SNP-h2 with our detected clozapine 387 transcriptomic signatures. While no association remained significant after correction of multiple 388 testing (72 tests, p < 6.9 * 10 -4 ), we did observe an increasing association between total cholesterol 389 and LDL heritability and clozapine gene-sets across more stringent thresholds of differentially 390 expressed genes (i.e. <FDR 5%, <FDR 1%, and Bonferroni correction ( Figure 4C and Table S9). 391 Such trends were not observed for HDL or any of the other traits tested, including SCZ. To further explore our gene expression findings, we investigated their association with genetic 470 susceptibility of schizophrenia. We did not observe significant enrichment of schizophrenia 471 heritability across clozapine-associated genes. Two previous studies did report an association 472 between schizophrenia heritability and antipsychotic target genes (Gaspar and Breen, 2017;473 Skene et al., 2018). Their approach, however, differed from ours. While they used lists of 474 antipsychotic target genes originating from pharmacological databases, we used a list of 475 experimentally derived gene sets after in vitro exposure to clozapine in LCLs. The genes we found 476 to be differentially expressed did not overlap with antipsychotic drug target genes used in these 477 two previous studies, which may explain the discrepancy in findings. If we assume that the 478 clozapine-induced gene expression differences are a proxy for adverse effects, the lack of 479 evidence of the enrichment analysis suggests that the genetic architecture of disease 480 susceptibility is likely independent from susceptibility to adverse effects. We explored this further 481 by examining the enrichment of genetic signal linked to cholesterol and cardiovascular disease-482 related traits such as body mass index, coronary artery disease, type 2 diabetes, and triglyceride 483 levels. We observed a nominally significant increased enrichment for total cholesterol and LDL 484 genetic signal in upregulated genes in response to clozapine exposure. While these findings did 485 not survive multiple testing correction and should be interpreted with caution, the observed 486 cholesterol and LDL heritability enrichment increased as we narrowed down on the most strongly 487 associated clozapine genes suggesting that this association warrants further investigation. in one tissue may thus be informative for other tissues. We observed that preferentially expressed 506 genes in GTEx-LCL tissue overlap the most with clozapine genes identified in LCLs in vitro, with 507 whole blood ranked as second highest tissue. Preferentially expressed genes however, represent 508 only a subset of the clozapine genes. Using all associated genes, we demonstrated that 509 upregulated clozapine genes have significantly different average expression in liver, muscle, lung, 510 and testis tissue. This suggests multi-tissue downstream effects through clozapine treatment, 511 which fits its clinical presentation of a diverse set of adverse effects that have been reported (Iqbal 512 et al., 2003). Our finding of higher expression in the liver is in line with the observed cholesterol 513 gene ontology signature of upregulated clozapine genes in our assay. The liver is a central player 514 in the function of cholesterol in the body and clozapine-related hepatotoxicity has furthermore 515 been reported in cases of treatment with the drug (Keane et al., 2009;Kellner et al., 1993). We 516 observe no differences for any of the GTEx brain tissues. 517 518 On top of within tissue effects, we also performed a more explorative analysis to examine 519 clozapine gene expression-based similarity between tissues and found specific tissue pairs that 520 were significantly different in tissue similarities. Here, we highlight some interesting observations. 521 We observe a decrease in similarity for tissues involved in immune (spleen), endocrine (testis, 522 ovary, adrenal and thyroid gland) and metabolic (adipose) functioning. The spleen and testis in 523 particular stand out with a significant change in dissimilarity with many tissues. The spleen is the 524 largest secondary lymphoid organ in the body and hosts a wide range of immunological functions, 525 including the storage of leukocytes (Lewis et al., 2019). While it remains speculative if splenic 526 dysfunction could for example lead to agranulocytosis, our findings do solicit for more research 527 on the mechanism between clozapine and splenic function. An intriguing, and perhaps also 528 surprising, finding is the implication of the testis among analysis of genes downregulated by 529 clozapine. The testis has two primary functions; to produce sperm and to produce hormones, in 530 particular testosterone, which is a sex steroid synthesized from cholesterol (Eacker et al., 2008). 531 Genes downregulated by clozapine are enriched for cell division processes and cell proliferation 532 is an important function of testicular cells (Sohni et al., 2019). While little is currently known about 533 how clozapine or other antipsychotics may affect testicular function in adult humans, our finding 534 does again point to cholesterol-related biology. In addition to the testis, we also observe deviation 535 in between-tissue similarity for several other endocrine tissues. As endocrine abnormalities are 536 known causes of human obesity, the identified transcriptomic profile of clozapine response 537 together with downstream tissue effects may point to new avenues to study clozapine-induced 538 weight gain (Baptista, 1999). Together, these findings suggest that clozapine disproportionately 539 affects the function of specific tissues. Our results furthermore demonstrate how experimental 540 molecular signatures can be integrated with external genomic datasets, such as the GTEx project, 541 to help translate in vitro findings to human biology.