Introduction
Bipolar disorder (BD) is a chronic and recurrent psychiatric disorder affecting ~1% of the population worldwide and presenting a major public health burden (Weissman et al., Reference Weissman, Bland, Canino, Faravelli, Greenwald, Hwu, Joyce, Karam, Lee, Lellouch, Lepine, Newman, Rubio-Stipec, Wells, Wickramaratne, Wittchen and Yeh1996; Eaton et al., Reference Eaton, Martins, Nestadt, Bienvenu, Clarke and Alexandre2008). It is characterized clinically by instability in mood resulting in manic and depressive episodes interspersed between neutral and euthymic states (Eaton et al., Reference Eaton, Martins, Nestadt, Bienvenu, Clarke and Alexandre2008). Risk for BD is highly genetic, with heritability estimates as high as 85% (McGuffin et al., Reference Mcguffin, Rijsdijk, Andrew, Sham, Katz and Cardno2003) and common genetic variation explaining up to a third of that (Cross-Disorder Group of the Psychiatric Genomics et al., Reference Lee, Ripke, Neale, Faraone, Purcell, Perlis, Mowry, Thapar, Goddard, Witte, Absher, Agartz, Akil, Amin, Andreassen, Anjorin, Anney, Anttila, Arking, Asherson, Azevedo, Backlund, Badner, Bailey, Banaschewski, Barchas, Barnes, Barrett, Bass, Battaglia, Bauer, Bayes, Bellivier, Bergen, Berrettini, Betancur, Bettecken, Biederman, Binder, Black, Blackwood, Bloss, Boehnke, Boomsma, Breen, Breuer, Bruggeman, Cormican, Buccola, Buitelaar, Bunney, Buxbaum, Byerley, Byrne, Caesar, Cahn, Cantor, Casas, Chakravarti, Chambert, Choudhury, Cichon, Cloninger, Collier, Cook, Coon, Cormand, Corvin, Coryell, Craig, Craig, Crosbie, Cuccaro, Curtis, Czamara, Datta, Dawson, Day, De Geus, Degenhardt, Djurovic, Donohoe, Doyle, Duan, Dudbridge, Duketis, Ebstein, Edenberg, Elia, Ennis, Etain, Fanous, Farmer, Ferrier, Flickinger, Fombonne, Foroud, Frank, Franke, Fraser, Freedman, Freimer, Freitag, Friedl, Frisen, Gallagher, Gejman, Georgieva, Gershon, Geschwind, Giegling, Gill, Gordon, Gordon-Smith, Green, Greenwood, Grice, Gross, Grozeva, Guan, Gurling, De Haan, Haines, Hakonarson, Hallmayer, Hamilton, Hamshere, Hansen, Hartmann, Hautzinger, Heath, Henders, Herms, Hickie, Hipolito, Hoefels, Holmans, Holsboer, Hoogendijk, Hottenga, Hultman, Hus, Ingason, Ising, Jamain, Jones, Jones, Jones, Tzeng, Kahler, Kahn, Kandaswamy, Keller, Kennedy, Kenny, Kent, Kim, Kirov, Klauck, Klei, Knowles, Kohli, Koller, Konte, Korszun, Krabbendam, Krasucki, Kuntsi, Kwan, Landen, Langstrom, Lathrop, Lawrence, Lawson, Leboyer, Ledbetter, Lee, Lencz, Lesch, Levinson, Lewis, Li, Lichtenstein, Lieberman, Lin, Linszen, Liu, Lohoff, Loo, Lord, Lowe, Lucae, MacIntyre, Madden, Maestrini, Magnusson, Mahon, Maier, Malhotra, Mane, Martin, Martin, Mattheisen, Matthews, Mattingsdal, McCarroll, McGhee, McGough, McGrath, McGuffin, McInnis, McIntosh, McKinney, McLean, McMahon, McMahon, McQuillin, Medeiros, Medland, Meier, Melle, Meng, Meyer, Middeldorp, Middleton, Milanova, Miranda, Monaco, Montgomery, Moran, Moreno-De-Luca, Morken, Morris, Morrow, Moskvina, Muglia, Muhleisen, Muir, Muller-Myhsok, Murtha, Myers, Myin-Germeys, Neale, Nelson, Nievergelt, Nikolov, Nimgaonkar, Nolen, Nothen, Nurnberger, Nwulia, Nyholt, O'Dushlaine, Oades, Olincy, Oliveira, Olsen, Ophoff, Osby, Owen, Palotie, Parr, Paterson, Pato, Pato, Penninx, Pergadia, Pericak-Vance, Pickard, Pimm, Piven, Posthuma, Potash, Poustka, Propping, Puri, Quested, Quinn, Ramos-Quiroga, Rasmussen, Raychaudhuri, Rehnstrom, Reif, Ribases, Rice, Rietschel, Roeder, Roeyers, Rossin, Rothenberger, Rouleau, Ruderfer, Rujescu, Sanders, Sanders, Santangelo, Sergeant, Schachar, Schalling, Schatzberg, Scheftner, Schellenberg, Scherer, Schork, Schulze, Schumacher, Schwarz, Scolnick, Scott, Shi, Shilling, Shyn, Silverman, Slager, Smalley, Smit, Smith, Sonuga-Barke, St Clair, State, Steffens, Steinhausen, Strauss, Strohmaier, Stroup, Sutcliffe, Szatmari, Szelinger, Thirumalai, Thompson, Todorov, Tozzi, Treutlein, Uhr, van den Oord, Van Grootheest, Van Os, Vicente, Vieland, Vincent, Visscher, Walsh, Wassink, Watson, Weissman, Werge, Wienker, Wijsman, Willemsen, Williams, Willsey, Witt, Xu, Young, Yu, Zammit, Zandi, Zhang, Zitman, Zollner, Devlin, Kelsoe, Sklar, Daly, O'Donovan, Craddock, Sullivan, Smoller, Kendler and Wray2013). Still, however, the pathophysiological characteristics of BD are not well understood. Investigating molecular phenotypes such as gene expression as intermediate measures between genetic variation and clinical variation is a viable strategy for uncovering disease mechanisms. Many such studies have been carried out for BD, and we present a summary that reveals a lack of consistency between findings likely owing to clinical heterogeneity, differing study designs, and the low numbers of samples investigated (N 62 BD subjects; online Supplementary Table S1) (Elashoff et al., Reference Elashoff, Higgs, Yolken, Knable, Weis, Webster, Barci and Torrey2007; Matigian et al., Reference Matigian, Windus, Smith, Filippich, Pantelis, Mcgrath, Mowry and Hayward2007; Choi et al., Reference Choi, Higgs, Wendland, Song, Mcmahon and Webster2011; Akula et al., Reference Akula, Barb, Jiang, Wendland, Choi, Sen, Hou, Chen, Laje, Johnson, Lipska, Kleinman, Corrada-Bravo, Detera-Wadleigh, Munson and Mcmahon2014; Beech et al., Reference Beech, Leffert, Lin, Sylvia, Umlauf, Mane, Zhao, Bowden, Calabrese, Friedman, Ketter, Iosifescu, Reilly-Harrington, Ostacher, Thase and Nierenberg2014; Mostafavi et al., Reference Mostafavi, Battle, Zhu, Potash, Weissman, Shi, Beckman, Haudenschild, Mccormick, Mei, Gameroff, Gindes, Adams, Goes, Mondimore, Mackinnon, Notes, Schweizer, Furman, Montgomery, Urban, Koller and Levinson2014; Witt et al., Reference Witt, Juraeva, Sticht, Strohmaier, Meier, Treutlein, Dukal, Frank, Lang, Deuschle, Schulze, Degenhardt, Mattheisen, Brors, Cichon, Nothen, Witt and Rietschel2014; Xiao et al., Reference Xiao, Camarillo, Ping, Arana, Zhao, Thompson, Xu, Su, Fan, Ordonez, Wang, Mao, Zhang, Cruz, Escamilla, Li and Xu2014; Cruceanu et al., Reference Cruceanu, Tan, Rogic, Lopez, Torres-Platas, Gigek, Alda, Rouleau, Pavlidis and Turecki2015; Madison et al., Reference Madison, Zhou, Nigam, Hussain, Barker, Nehme, Van Der Ven, Hsu, Wolf, Fleishman, O'dushlaine, Rose, Chambert, Lau, Ahfeldt, Rueckert, Sheridan, Fass, Nemesh, Mullen, Daheron, Mccarroll, Sklar, Perlis and Haggarty2015; Mertens et al., Reference Mertens, Wang, Kim, Yu, Pham, Yang, Zheng, Diffenderfer, Zhang, Soltani, Eames, Schafer, Boyer, Marchetto, Nurnberger, Calabrese, Odegaard, Mccarthy, Zandi, Alda, Nievergelt, Mi, Brennand, Kelsoe, Gage and Yao2015; van Eijk et al., Reference Van Eijk, De Jong, Strengman, Buizer-Voskamp, Kahn, Boks, Horvath and Ophoff2015; Zhao et al., Reference Zhao, Xu, Chen, Kim, Reimers, Bacanu, Yu, Liu, Sun, Wang, Jia, Xu, Zhang, Kendler, Peng and Chen2015; Anand et al., Reference Anand, Mcclintick, Murrell, Karne, Nurnberger and Edenberg2016; Breen et al., Reference Breen, White, Shekhtman, Lin, Looney, Woelk and Kelsoe2016; Fromer et al., Reference Fromer, Roussos, Sieberts, Johnson, Kavanagh, Perumal, Ruderfer, Oh, Topol, Shah, Klei, Kramer, Pinto, Gumus, Cicek, Dang, Browne, Lu, Xie, Readhead, Stahl, Xiao, Parvizi, Hamamsy, Fullard, Wang, Mahajan, Derry, Dudley, Hemby, Logsdon, Talbot, Raj, Bennett, De Jager, Zhu, Zhang, Sullivan, Chess, Purcell, Shinobu, Mangravite, Toyoshiba, Gur, Hahn, Lewis, Haroutunian, Peters, Lipska, Buxbaum, Schadt, Hirai, Roeder, Brennand, Katsanis, Domenici, Devlin and Sklar2016; Hess et al., Reference Hess, Tylee, Barve, De Jong, Ophoff, Kumarasinghe, Tooney, Schall, Gardiner, Beveridge, Scott, Yasawardene, Perera, Mendis, Carr, Kelly, Cairns, Tsuang and Glatt2016; Jansen et al., Reference Jansen, Penninx, Madar, Xia, Milaneschi, Hottenga, Hammerschlag, Beekman, Van Der Wee, Smit, Brooks, Tischfield, Posthuma, Schoevers, Van Grootheest, Willemsen, De Geus, Boomsma, Wright, Zou, Sun and Sullivan2016; Peterson et al., Reference Peterson, Service, Jasinska, Gao, Zelaya, Teshiba, Bearden, Cantor, Reus, Macaya, Lopez-Jaramillo, Bogomolov, Benjamini, Eskin, Coppola, Freimer and Sabatti2016; Fries et al., Reference Fries, Colpo, Monroy-Jaramillo, Zhao, Zhao, Arnold, Bowden and Walss-Bass2017; Kittel-Schneider et al., Reference Kittel-Schneider, Hilscher, Scholz, Weber, Grunewald, Schwarz, Chiocchetti and Reif2017; Pacifico and Davis, Reference Pacifico and Davis2017; Vizlin-Hodzic et al., Reference Vizlin-Hodzic, Zhai, Illes, Sodersten, Truve, Parris, Sobhan, Salmela, Kosalai, Kanduri, Strandberg, Seth, Bontell, Hanse, Agren and Funa2017). Moreover, there are many potential confounds that impact gene expression, including medication.
Therefore, to explore gene expression changes associated with BD, we generated RNA sequencing data from peripheral whole blood collected in a large case-control cohort from The Netherlands. We examined gene expression differences between groups both at the individual gene level and at the level of gene co-expression to shed light on disease-relevant molecular profiles. We also investigated the effects of lithium use, the most widely used prescription drug in our cohort, and other variables on gene expression reflect the complex, polygenic nature of bipolar disorder as measured by genetic risk identified through genome-wide association studies (GWASs) of psychiatric disorders. The main findings suggest that there are nominal BD-related gene expression effects in blood but numerous effects related to lithium treatment, and that these changes are independent of genetic risk of psychiatric disorders. This work highlights the importance of incorporating medication use in psychiatric transcriptomic studies and provides insight into blood-based gene expression effects related to BD and the effects of lithium.
Methods
Participant recruitment
Data were generated according to protocols approved by the respective local ethics committees: the Medical Ethical Review Board at University Medical Center Utrecht and the Institutional Review Board at University of California Los Angeles. Informed consent was obtained from all subjects. Participants were included upon the criteria of having at least three Dutch grandparents and being older than 18 years of age. Patients were recruited via clinicians, the Dutch patients' association, pharmacies, and advertisements as previously described (Abramovic et al., Reference Abramovic, Boks, Vreeker, Bouter, Kruiper, Verkooijen, Van Bergen, Ophoff, Kahn and Van Haren2016; Abramovic et al., Reference Abramovic, Boks, Vreeker, Verkooijen, Van Bergen, Ophoff, Kahn and Van Haren2018). Case inclusion criteria included having a diagnosis of BD-I or BD-II and a current euthymic state. Diagnosis was confirmed via assessment as previously described (Abramovic et al., Reference Abramovic, Boks, Vreeker, Bouter, Kruiper, Verkooijen, Van Bergen, Ophoff, Kahn and Van Haren2016) with the Structured Clinical Interview for DSM-IV (http://www.scid4.org). A portion of the controls (ascertainment group A) was recruited via advertisements and involvement in previously studies after having agreed to be re-contacted for new research as previously described (Abramovic et al., Reference Abramovic, Boks, Vreeker, Bouter, Kruiper, Verkooijen, Van Bergen, Ophoff, Kahn and Van Haren2016; Abramovic et al., Reference Abramovic, Boks, Vreeker, Verkooijen, Van Bergen, Ophoff, Kahn and Van Haren2018). Another portion of the controls (ascertainment group B) was recruited at outpatient preoperative screening services in four hospitals in and around Utrecht, Netherlands as previously described (Luykx et al., Reference Luykx, Bakker, Lentjes, Neeleman, Strengman, Mentink, Deyoung, De Jong, Sul, Eskin, Van Eijk, Van Setten, Buizer-Voskamp, Cantor, Lu, Van Amerongen, Van Dongen, Keijzers, Kappen, Borgdorff, Bruins, Derks, Kahn and Ophoff2014). Control subjects did not have a diagnosis of BD or any psychotic or neurological disorder and had no first-degree relative with a diagnosis of BD or any psychotic disorder. In addition to diagnostic assessments, subjects were assessed for medication and tobacco use.
Lithium use assessment
Information about patients' lithium use was gathered as previously described (Abramovic et al., Reference Abramovic, Boks, Vreeker, Bouter, Kruiper, Verkooijen, Van Bergen, Ophoff, Kahn and Van Haren2016). In brief, patients lithium treatment prescribed and monitored by their own physician was self-reported on in three ways: (1) in an online medical questionnaire that inquired about medication use, (2) during the on-site assessment where a list of current and lifetime medication use was discussed, and (3) in an assessment of a lithium satisfaction questionnaire. The data of these three measures was combined to accurately determine the current use of lithium in subjects with BD, which constituted the lithium use phenotype in subsequent analyses. Although data were missing for seven individuals, information regarding response to lithium showed that a majority of subjects being treated with lithium had experienced less frequent (N = 104, 71.7%) and less severe (N = 113, 77.9%) mood episodes, and were satisfied or very satisfied with the use of lithium (N = 113, 77.9%) since starting the medication. Additionally, although data were missing for 16 individuals, information regarding past lithium use indicated that only ten subjects had never used lithium (4.5%). Information regarding subjects' blood lithium levels was available for 82 of the 152 lithium users via self-report based on their most recent medication check with their doctor (not at the time of blood draw for the current study). Levels ranged between 0.34 and 1.2 mmol/l with a mean of 0.78 mmol/l, and 73 subjects had blood levels within the therapeutic range (0.6 to 1.2 mmol/l). Of the 152 subjects using lithium, 59 were using lithium with another mood stabilizer and 93 were using lithium as their only mood stabilizer, however, information regarding which other mood stabilizer used was unavailable. Data about the use of antipsychotics was too sparse to be used in subsequent analyses.
Gene expression quantification and differential expression analysis
Peripheral whole blood was drawn and processed for genotyping and RNA sequencing from 240 controls and 240 cases yielding an average of 24.9 million paired-end reads per sample, which were then mapped to human reference genome hg19 using TopHat2 (Kim et al., Reference Kim, Pertea, Trapnell, Pimentel, Kelley and Salzberg2013). Additional details about sample preparation, RNA sequencing, and read alignment can be found in the online Supplementary Methods. Known Ensembl gene levels were quantified using HTSeq in the union mode to obtain integral counts of reads that intersect the union of all transcripts of genes. Principal component analysis of gene expression quantification was used for data visualization and additional QC, after which four samples were removed for apparent mix-up (online Supplementary Methods). Thirty-two additional samples were excluded due to incomplete demographic information. Differential expression and co-expression analyses were therefore limited to a set of 444 subjects (240 cases and 204 controls).
Gene expression counts from HTSeq were filtered for genes having >10 counts in 90% of samples, yielding 12 344 genes for subsequent analyses, of which 1796 were non-coding. A standard differential expression analysis for gene expression counts generated from RNA sequencing was performed using limma voom (Law et al., Reference Law, Chen, Shi and Smyth2014). Details of this analysis can be found in the online Supplementary Methods. p Values were corrected for multiple testing using the Benjamini–Hochberg false discovery rate (FDR) estimation, and a gene was considered to be differentially expressed if it had an FDR-corrected p value < 0.05. An overview of covariates can be found in Table 1 and a description of the final covariate models can be found in the online Supplementary Methods. Differentially expressed genes (DEGs) were checked for overlap and concordance with other datasets (online Supplementary Methods). Fold changes (FCs) reported are in log 2 FC units. Functional annotation of DEGs was performed as described in the online Supplementary Methods.
Demographic and technical variables
PC, principal component; s.d., standard deviation.
p Values computed by Fisher's exact test (categorical variable) or t test (continuous variable).
Co-expression network analysis
To prepare the data for co-expression network analysis, the 12 344 filtered and normalized genes were residualized adjusting for covariates. To determine networks of genes with correlated expression, weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, Reference Langfelder and Horvath2008) was performed using the WGCNA package in R. WGCNA defines a network of genes as nodes with edges between genes based on pairwise correlations between genes, and separates the network into modules of gene clusters with highly coordinated expression. The gene expression profiles of each module were summarized by calculating the module eigengene, which is defined as the first principal component of the expression matrix of that module, and each gene was then assigned a measure of module membership for each module. To determine biologically significant modules, gene significance measures were assigned to each gene for each of our traits of interest, including BD diagnosis and lithium use, by calculating the absolute correlation between the trait and the expression profiles. Then a measure of module-trait significance was calculated by correlating module membership values with gene significance values. An association was considered significant if its p value surpassed Bonferroni correction for testing multiple modules (p < = 0.05/N modules). Additional details of this analysis and functional annotation of modules are described in the online Supplementary Methods.
Cell-type proportion estimation and enrichment
To estimate cell-type composition in our sample we employed the CIBERSORT online software (cibersort.stanford.edu) (Newman et al., Reference Newman, Liu, Green, Gentles, Feng, Xu, Hoang, Diehn and Alizadeh2015). Details of this analysis can be found in the online Supplementary Methods. The resulting estimated cell-type proportions were regressed on covariates, and the residuals were used to predict lithium use in a stepwise linear regression using the stepAIC function in the MASS package in R. The estimated cell-type proportions were also appended to the table of technical and biological covariates and then used to re-run the differential expression analysis while accounting for cell-type heterogeneity in the sample. The enrichment of LM22 cell types in gene co-expression modules determined from WGCNA was calculated in two ways: first via hypergeometric overlap between modules and cell type signature genes, and second by using binary cell type signatures to predict module membership values in a linear model. Details of these analyses are described in online Supplementary Methods.
Integration of GWAS data with transcriptomic signatures
Prior to gene-set analyses, heritability and genetic correlation of traits of interest were estimated to confirm significant non-zero single-nucleotide polymorphism (SNP)-based heritability (online Supplementary Methods). Analyses were performed across three psychiatric GWAS traits from publicly available datasets (bipolar disorder, schizophrenia (SCZ), and self-reported depression) and two sets of DEGs (BD at FDR < 0.2 and lithium-use at FDR < 0.05). Differential expression log 2 FCs and FDR-corrected p values for each of the 12 344 genes expressed at >10 counts in 90% of samples were obtained from limma to integrate whole-blood gene expression signatures with GWAS data using Multi-marker Analysis of GenoMic Annotation (MAGMA v1.06) (de Leeuw et al., Reference De Leeuw, Mooij, Heskes and Posthuma2015).
GWAS summary statistics were obtained for the following three GWAS traits:
(1) SCZ (Schizophrenia Working Group of the Psychiatric Genomics, 2014): 36 989 cases and 113 075 controls;
(3) Depression (Hyde et al., Reference Hyde, Nagle, Tian, Chen, Paciga, Wendland, Tung, Hinds, Perlis and Winslow2016): 75 607 cases and 231 747 controls;
The 1000 Genomes Project Phase 3 release European reference panel (N = 503) was used to model linkage disequilibrium in all analyses (Genomes Project et al., Reference Auton, Brooks, Durbin, Garrison, Kang, Korbel, Marchini, Mccarthy, Mcvean and Abecasis2015). Eight gene lists were used from two different DEG models along with a positive and negative control:
(1) Lithium-use DEGs at FDR < 0.05: N = 897 genes;
(2) Up-regulated lithium-use DEGs at FDR < 0.05: N = 680 genes;
(3) Down-regulated lithium-use DEGs at FDR < 0.05: N = 217 genes;
(4) BD DEGs at FDR < 0.2: N = 630 genes;
(5) Up-regulated BD DEGs at FDR < 0.2: N = 389 genes;
(6) Down-regulated BD DEGs at FDR < 0.2: N = 241 genes;
(7) Positive control gene-set: the top 100 most significant genes from a random draw of N = 1000 using the BD GWAS gene-level test statistics;
(8) Negative control gene-set: a random draw of N = 1000 genes using the BD GWAS gene-level test-statistics.
The MAGMA software settings used for this analysis can be found in online Supplementary Methods. Secondary gene-set analyses were run on a limited number of DEG gene sets and additional, sleep-related GWAS traits (online Supplementary Methods).
Results
Minimal changes in bipolar disorder gene expression
To explore the transcriptomic signatures of BD, we first evaluated whether subjects with BD harbored transcriptional differences on a per gene level compared with controls. Of the 12 344 genes tested, only six were differentially expressed in BD after correcting for multiple testing (FDR < 0.05; Fig. 1a). The differences in expression were very small, with absolute FC ranging from 0.12 to 0.44. While the number of identified DEGs was too small to perform functional enrichment analysis, we did find that three of the six genes (COG4, DOCK3, and BBS9) were expressed in GTEx frontal cortex tissue (median TPM > 1) and show relatively stable expression across brain cell types except for DOCK3, which is enriched in neurons (FC relative to other cell types = 6.82; online Supplementary Table S4). All six genes were present in at least one of the three BD cortical gene expression datasets examined, and one of the genes was significantly up-regulated in BD cases v. controls, PVT1 (Gandal et al., Reference Gandal, Haney, Parikshak, Leppa, Ramaswami, Hartl, Schork, Appadurai, Buil, Werge, Liu, White, Horvath and Geschwind2018a, Reference Gandal, Zhang, Hadjimichael, Walker, Chen, Liu, Won, Van Bakel, Varghese, Wang, Shieh, Haney, Parhami, Belmont, Kim, Moran Losada, Khan, Mleczko, Xia, Dai, Wang, Yang, Xu, Fish, Hof, Warrell, Fitzgerald, White, Jaffe, Psych, Peters, Gerstein, Liu, Iakoucheva, Pinto and Geschwind2018b). Another gene, COG4, was differentially expressed in the same direction as the current study in BD individuals in the Stanley Genomics brain gene expression database, and was reported as differentially expressed in a schizophrenia mega-analysis of nine whole blood microarray datasets (Hess et al., Reference Hess, Tylee, Barve, De Jong, Ophoff, Kumarasinghe, Tooney, Schall, Gardiner, Beveridge, Scott, Yasawardene, Perera, Mendis, Carr, Kelly, Cairns, Tsuang and Glatt2016). Using polygenic risk scores (PRS) for BD as the differential expression trait of interest rather than the dichotomous case-control phenotype did not yield any significant genes, even though PRS did significantly differ between BD cases and controls (t = −3.42, p = 6.88 × 10−4; online Supplementary Fig. S1 and Supplementary Methods).
Differentially expressed genes. (a) Six BD DEGs. FC, log 2 FC; q, FDR-adjusted p < 0.05. (b) 976 genes differentially expressed between lithium users and non-lithium users (shown as blue triangles, FDR-adjusted p < 0.05; all other genes tested shown as light gray circles). (c) DAVID (Da Huang et al., Reference Da Huang, Sherman and Lempicki2009) functional annotation cluster enrichment of all 976 DEGs (upper) and 754 up-regulated DEGs (lower). Enrichment scores increase when the gene list is limited to up-regulated genes only. Clusters were considered significant if the enrichment score >1 and at least one term in the cluster survived Bonferroni correction for multiple testing.
Widespread subtle gene expression changes in lithium users
Following the same differential expression pipeline as above, we found 976 genes with small differences in gene expression between lithium users and non-lithium users (|FC| mean = 0.20, max = 0.82, s.d. = 0.10; Fig. 1b, online Supplementary File S1). These genes were enriched for biological terms related to calcium signaling and other signaling pathways, and immunity (Fig. 1c). To distinguish between up- and down-regulated gene pathways, we stratified genes by their direction of change in expression. The 754 up-regulated genes were annotated for many of the same terms as the full set but with greater enrichment scores, indicating that the up-regulated genes are driving the enrichment scores in the full set (Fig. 1c). Of the 976 lithium-use DEGs, 804 were expressed in GTEx frontal cortex samples (TPM > 1), and 488, 553, 503, 478, 512, and 403 were expressed in neurons, fetal astrocytes, mature astrocytes, oligodendrocytes, microglia/macrophages, and endothelia, respectively (FPKM > 1). However, none of these gene sets were significantly enriched (hypergeometric p > 0.05).
To validate our results, the 976 lithium-use DEGs were tested for overlap with lists of DEGs from similar studies found in the literature (online Supplementary Table S5). Although none of these studies has the same design as ours, we did find a significant overlap between our 976 lithium-use DEGs and the lists from two studies. In the first study (Anand et al., Reference Anand, Mcclintick, Murrell, Karne, Nurnberger and Edenberg2016), 35 DEGs were detected by comparing peripheral monocyte gene expression in subjects before and after lithium monotherapy. Of these 35 DEGs, 18 were shared with the current study [hypergeometric odd ratio (OR) 13.57, p = 4.66 × 10−12], and all 18 were concordant in direction (online Supplementary Fig. S4A). In the second study (Breen et al., Reference Breen, White, Shekhtman, Lin, Looney, Woelk and Kelsoe2016), DEGs were detected by comparing LCL gene expression before and after lithium treatment in vitro. Of the 1504 DEGs discovered, 134 were shared with our study (hypergeometric OR 1.27, p = 9.23 × 10−3), and 84.6% of these were concordant in direction (online Supplementary Fig. S4B). There were two genes shared between all three lists, RFX2 and SLC29A1. We report genes in these overlapping lists as high-confidence lithium-associated genes (online Supplementary File S1).
Next, in search of genes with differential co-expression, we constructed a gene expression network in the entire sample using WGCNA and assessed the detected modules for association with traits of interest. This network consisted of 27 modules ranging in size from 48 to 2760 genes (mean N genes = 441, online Supplementary File S2), five of which were significantly associated with lithium-use. No modules were associated with BD or any other clinical or technical variable (online Supplementary Table S6). Of the five modules associated with lithium use, three shared significant overlap with lithium-use DEGs (Table 2). M26 was most significantly associated with lithium (p = 2.00 × 10−4; online Supplementary Fig. S6A) but was not significantly enriched for lithium DEGs. M1 was also associated with lithium (p = 9.04 × 10−4; online Supplementary Fig. S6B) and had the most significant enrichment of DEGs (431 of 2092 genes in the module were DEGs; hypergeometric OR 4.62, p = 2.03 × 10−97). Functional annotation clustering of the genes in M1 showed an enrichment of terms related to cell signaling, immunity, and glycophosphatidylinositol anchor. Module preservation analysis was also performed to assess differences in network density and connectivity between groups, but showed full preservation indicating that networks constructed in separate groups maintain their underlying structure (online Supplementary Methods and Fig. S7).
Co-expression module association with lithium use
DEGs, differentially expressed genes; GPI, glycophosphatidylinositol.
Functional annotation cluster enrichment determined using DAVID (Da Huang et al., Reference Da Huang, Sherman and Lempicki2009). Correlation with lithium use calculated by correlating gene module membership values with gene significance values for lithium use. Overlap was calculated by testing for hypergeometric overlap between the list of lithium-use DEGs and the list of genes within each module.
Estimated neutrophil proportions are increased in lithium users
We then sought to determine if variation in our sample could be explained by differences in blood cell-type composition, which might represent a biologically meaningful effect of lithium. Because white blood cell counts were unavailable, to deconvolve cellular heterogeneity, we utilized CIBERSORT, a method that has been shown to accurately characterize blood cell composition based on gene expression profiles (Newman et al., Reference Newman, Liu, Green, Gentles, Feng, Xu, Hoang, Diehn and Alizadeh2015). Using our gene expression quantifications and a reference panel of 22 blood cell-type signatures, we estimated cell-type proportions (Fig. 2a) and examined their relationship with lithium use in BD cases only. Each cell type was regressed on demographic and technical variables then the residuals were used to predict lithium use in a stepwise linear model. Neutrophils are the one cell type that significantly predicted lithium use within the BD cases (β = 0.63, p = 0.024), with elevated proportions in individuals being treated with lithium (Fig. 2b). Indeed, 16 of 60 signature neutrophil genes were also lithium-use DEGs (hypergeometric OR 4.64, p = 4.45 × 10−6).
Estimated neutrophil composition association with lithium use. (a) Leukocyte cell-type proportions per sample as estimated from gene expression, sorted by neutrophil proportions. Mean proportion across samples shown in parentheses. Lithium users, shown in the bar on the bottom, cluster on the right where neutrophil proportions are higher. (b) Lithium users have higher estimated neutrophil proportions (β = 0.63, p = 0.024).
The number of genes showing differential expression in subjects undergoing lithium treatment decreased from 976 in the model without cell-type estimates to 233 in the model with cell-type estimates (FDR < 0.05; online Supplementary Fig. S8A and File S1), of which 194 (83.2%) were significant in the original model and concordant in the direction of effect (online Supplementary Fig. S8B). No functional annotation cluster terms remained significant after correcting for multiple testing. The number of genes differentially expressed between BD cases and controls decreased to zero after accounting for estimated cell-type proportions.
We then sought to determine if the various lithium-associated modules of co-expressed genes reflected biologic signatures of distinct populations of blood cell types. We did this in two ways. First, a hypergeometric overlap between lithium-associated module gene lists and cell-type signature gene lists revealed a significant overlap between module M1 with monocyte and neutrophil signature genes and M9 with eosinophil and activated mast cell signature genes (Fig. 3a, left). Second, the expression of cell-type signature genes was used to predict module membership values in a linear model for each of the five lithium-associated modules. Neutrophils, monocytes, and eosinophils were again implicated (Fig. 3a, right). In both of these analyses, the most significant cell type-module relationship was M1 with neutrophil estimates (hypergeometric p = 5.68 × 10−21, linear model p < 2.20 × 10−16). Indeed, neutrophil signature genes had higher M1 membership values (Fig. 3b).
Lithium-associated co-expression module M1 enrichment for neutrophil gene expression signatures. (a) Lithium-associated module enrichment for leukocyte cell types. Left, Hypergeometric overlap between leukocyte cell type signature genes and genes in each module. Right, Linear regression of leukocyte cell type signature genes to predict module membership values. (b) Neutrophil signature genes have higher module membership values for M1 than other leukocyte signature genes (β = 0.60, p < 2.20 × 10−16).
Genes with altered expression are not enriched for genes with common psychiatric risk alleles
To evaluate if BD and lithium-use DEG sets were associated with a higher burden of psychiatric risk alleles, we performed gene-set analyses using MAGMA (de Leeuw et al., Reference De Leeuw, Mooij, Heskes and Posthuma2015). Analyses were performed across three psychiatric GWAS traits: BD (Stahl et al., Reference Stahl, Breen, Forstner, Mcquillin, Ripke, Trubetskoy, Mattheisen, Wang, Coleman, Gaspar, De Leeuw, Steinberg, Pavlides, Trzaskowski, Byrne, Pers, Holmans, Richards, Abbott, Agerbo, Akil, Albani, Alliey-Rodriguez, Als, Anjorin, Antilla, Awasthi, Badner, Baekvad-Hansen, Barchas, Bass, Bauer, Belliveau, Bergen, Pedersen, Boen, Boks, Boocock, Budde, Bunney, Burmeister, Bybjerg-Grauholm, Byerley, Casas, Cerrato, Cervantes, Chambert, Charney, Chen, Churchhouse, Clarke, Coryell, Craig, Cruceanu, Curtis, Czerski, Dale, De Jong, Degenhardt, Del-Favero, Depaulo, Djurovic, Dobbyn, Dumont, Elvsashagen, Escott-Price, Fan, Fischer, Flickinger, Foroud, Forty, Frank, Fraser, Freimer, Frisen, Gade, Gage, Garnham, Giambartolomei, Pedersen, Goldstein, Gordon, Gordon-Smith, Green, Green, Greenwood, Grove, Guan, Guzman-Parra, Hamshere, Hautzinger, Heilbronner, Herms, Hipolito, Hoffmann, Holland, Huckins, Jamain, Johnson, Jureus, Kandaswamy, Karlsson, Kennedy, Kittel-Schneider, Knowles, Kogevinas, Koller, Kupka, Lavebratt, Lawrence, Lawson, Leber, Lee, Levy, Li, Liu, Lucae, Maaser, MacIntyre, Mahon, Maier, Martinsson, McCarroll, McGuffin, McInnis, McKay, Medeiros, Medland, Meng, Milani, Montgomery, Morris, Muhleisen, Mullins, Nguyen, Nievergelt, Adolfsson, Nwulia, O'Donovan, Loohuis, Ori, Oruc, Osby, Perlis, Perry, Pfennig, Potash, Purcell, Regeer, Reif, Reinbold, Rice, Rivas, Rivera, Roussos, Ruderfer, Ryu, Sanchez-Mora, Schatzberg, Scheftner, Schork, Shannon Weickert, Shehktman, Shilling, Sigurdsson, Slaney, Smeland, Sobell, Soholm Hansen, Spijker, St Clair, Steffens, Strauss, Streit, Strohmaier, Szelinger, Thompson, Thorgeirsson, Treutlein, Vedder, Wang, Watson, Weickert, Witt, Xi, Xu, Young, Zandi, Zhang, Zollner, e, Consortium, Adolfsson, Agartz, Alda, Backlund, Baune, Bellivier, Berrettini, Biernacka, Blackwood, Boehnke, Borglum, Corvin, Craddock, Daly, Dannlowski, Esko, Etain, Frye, Fullerton, Gershon, Gill, Goes, Grigoroiu-Serbanescu, Hauser, Hougaard, Hultman, Jones, Jones, Kahn, Kirov, Landen, Leboyer, Lewis, Li, Lissowska, Martin, Mayoral, McElroy, McIntosh, McMahon, Melle, Metspalu, Mitchell, Morken, Mors, Mortensen, Muller-Myhsok, Myers, Neale, Nimgaonkar, Nordentoft, Nothen, O'Donovan, Oedegaard, Owen, Paciga, Pato, Pato, Posthuma, Ramos-Quiroga, Ribases, Rietschel, Rouleau, Schalling, Schofield, Schulze, Serretti, Smoller, Stefansson, Stefansson, Stordal, Sullivan, Turecki, Vaaler, Vieta, Vincent, Werge, Nurnberger, Wray, Di Florio, Edenberg, Cichon, Ophoff, Scott, Andreassen, Kelsoe and Sklar2019), SCZ (Schizophrenia Working Group of the Psychiatric Genomics, 2014), and self-reported depression (Hyde et al., Reference Hyde, Nagle, Tian, Chen, Paciga, Wendland, Tung, Hinds, Perlis and Winslow2016). SCZ and depression were used because of their high degree of overlap in SNP-based heritability with BD (Cross-Disorder Group of the Psychiatric Genomics et al., Reference Lee, Ripke, Neale, Faraone, Purcell, Perlis, Mowry, Thapar, Goddard, Witte, Absher, Agartz, Akil, Amin, Andreassen, Anjorin, Anney, Anttila, Arking, Asherson, Azevedo, Backlund, Badner, Bailey, Banaschewski, Barchas, Barnes, Barrett, Bass, Battaglia, Bauer, Bayes, Bellivier, Bergen, Berrettini, Betancur, Bettecken, Biederman, Binder, Black, Blackwood, Bloss, Boehnke, Boomsma, Breen, Breuer, Bruggeman, Cormican, Buccola, Buitelaar, Bunney, Buxbaum, Byerley, Byrne, Caesar, Cahn, Cantor, Casas, Chakravarti, Chambert, Choudhury, Cichon, Cloninger, Collier, Cook, Coon, Cormand, Corvin, Coryell, Craig, Craig, Crosbie, Cuccaro, Curtis, Czamara, Datta, Dawson, Day, De Geus, Degenhardt, Djurovic, Donohoe, Doyle, Duan, Dudbridge, Duketis, Ebstein, Edenberg, Elia, Ennis, Etain, Fanous, Farmer, Ferrier, Flickinger, Fombonne, Foroud, Frank, Franke, Fraser, Freedman, Freimer, Freitag, Friedl, Frisen, Gallagher, Gejman, Georgieva, Gershon, Geschwind, Giegling, Gill, Gordon, Gordon-Smith, Green, Greenwood, Grice, Gross, Grozeva, Guan, Gurling, De Haan, Haines, Hakonarson, Hallmayer, Hamilton, Hamshere, Hansen, Hartmann, Hautzinger, Heath, Henders, Herms, Hickie, Hipolito, Hoefels, Holmans, Holsboer, Hoogendijk, Hottenga, Hultman, Hus, Ingason, Ising, Jamain, Jones, Jones, Jones, Tzeng, Kahler, Kahn, Kandaswamy, Keller, Kennedy, Kenny, Kent, Kim, Kirov, Klauck, Klei, Knowles, Kohli, Koller, Konte, Korszun, Krabbendam, Krasucki, Kuntsi, Kwan, Landen, Langstrom, Lathrop, Lawrence, Lawson, Leboyer, Ledbetter, Lee, Lencz, Lesch, Levinson, Lewis, Li, Lichtenstein, Lieberman, Lin, Linszen, Liu, Lohoff, Loo, Lord, Lowe, Lucae, MacIntyre, Madden, Maestrini, Magnusson, Mahon, Maier, Malhotra, Mane, Martin, Martin, Mattheisen, Matthews, Mattingsdal, McCarroll, McGhee, McGough, McGrath, McGuffin, McInnis, McIntosh, McKinney, McLean, McMahon, McMahon, McQuillin, Medeiros, Medland, Meier, Melle, Meng, Meyer, Middeldorp, Middleton, Milanova, Miranda, Monaco, Montgomery, Moran, Moreno-De-Luca, Morken, Morris, Morrow, Moskvina, Muglia, Muhleisen, Muir, Muller-Myhsok, Murtha, Myers, Myin-Germeys, Neale, Nelson, Nievergelt, Nikolov, Nimgaonkar, Nolen, Nothen, Nurnberger, Nwulia, Nyholt, O'Dushlaine, Oades, Olincy, Oliveira, Olsen, Ophoff, Osby, Owen, Palotie, Parr, Paterson, Pato, Pato, Penninx, Pergadia, Pericak-Vance, Pickard, Pimm, Piven, Posthuma, Potash, Poustka, Propping, Puri, Quested, Quinn, Ramos-Quiroga, Rasmussen, Raychaudhuri, Rehnstrom, Reif, Ribases, Rice, Rietschel, Roeder, Roeyers, Rossin, Rothenberger, Rouleau, Ruderfer, Rujescu, Sanders, Sanders, Santangelo, Sergeant, Schachar, Schalling, Schatzberg, Scheftner, Schellenberg, Scherer, Schork, Schulze, Schumacher, Schwarz, Scolnick, Scott, Shi, Shilling, Shyn, Silverman, Slager, Smalley, Smit, Smith, Sonuga-Barke, St Clair, State, Steffens, Steinhausen, Strauss, Strohmaier, Stroup, Sutcliffe, Szatmari, Szelinger, Thirumalai, Thompson, Todorov, Tozzi, Treutlein, Uhr, van den Oord, Van Grootheest, Van Os, Vicente, Vieland, Vincent, Visscher, Walsh, Wassink, Watson, Weissman, Werge, Wienker, Wijsman, Willemsen, Williams, Willsey, Witt, Xu, Young, Yu, Zammit, Zandi, Zhang, Zitman, Zollner, Devlin, Kelsoe, Sklar, Daly, O'Donovan, Craddock, Sullivan, Smoller, Kendler and Wray2013) (online Supplementary Table S7). The 23andMe self-reported depression GWAS was used instead of MDD GWAS because of the large sample size and successful findings of this study. A lithium-response GWAS was not used because the SNP-based heritability estimate for this trait is not different from zero (personal communication with Drs. Thomas G. Schulze and Francis McMahon). Because the set of BD DEGs at FDR < 0.05 was too small to test, we used a more lenient significance threshold of FDR < 0.2 for this analysis instead. None of the comparisons demonstrated an association with genetic risk across the genes identified in the current study (except for the positive control gene set), even after stratifying by up- and down-regulated genes (online Supplementary Fig. S9 and Table S8). Because sleep disturbances are a hallmark of BD (Harvey, Reference Harvey2008), and due to the genetic correlation of sleep-related phenotypes with BD (Jones et al., Reference Jones, Tyrrell, Wood, Beaumont, Ruth, Tuke, Yaghootkar, Hu, Teder-Laving, Hayward, Roenneberg, Wilson, Del Greco, Hicks, Shin, Yun, Lee, Metspalu, Byrne, Gehrman, Tiemeier, Allebrandt, Freathy, Murray, Hinds, Frayling and Weedon2016), we performed a secondary gene-set analysis with genes implicated from chronotype, sleep duration, oversleeping, and undersleeping GWAS, which failed to demonstrate association with genes identified in the current study (online Supplementary Table S9).
Discussion
In our whole blood BD case-control gene expression study we observed widespread subtle changes in gene expression in subjects undergoing lithium treatment but few transcriptomic differences in euthymic BD cases compared to controls. These effects were partially driven by variation in leukocyte cell type composition, and we find no evidence for a link with genetic risk for BD. Upon validation of our findings with previous in vivo and in vitro lithium treatment gene expression studies, we present a high-confidence list of genes that display altered expression associated with lithium treatment.
Lithium is the first-line treatment for BD, not only for the treatment of acute episodes but also for maintenance and suicide prevention (Cipriani et al., Reference Cipriani, Hawton, Stockton and Geddes2013; Malhi et al., Reference Malhi, Gessler and Outhred2017). However, only about 30% of BD patients fully respond to lithium, it has several adverse side effects, and its mechanisms of action are not well understood (Alda, Reference Alda2015; Gitlin, Reference Gitlin2016; Pickard, Reference Pickard2017). One probable reason for this lack of understanding is the magnitude of lithium's physiological interactions (Roux and Dosseto, Reference Roux and Dosseto2017). In pharmacological terms, lithium is a small molecule (the third smallest element in fact) without a defined target (Pickard, Reference Pickard2017). This lack of specificity makes it difficult to discern therapeutic mechanisms from off-target effects, which likely lead to many of lithium's undesirable side effects and even its toxicity in serum concentrations that are above therapeutic levels. Lithium ions (Li+) have a single positive charge and are hypothesized to mimic and disrupt the actions and targets of more ubiquitous metal ions such as magnesium (Mg2+) (Pickard, Reference Pickard2017). Theorized therapeutic mechanisms of lithium include its inhibition of the protein GS3Kβ, and its effect on intracellular signaling cascades such as those involving protein kinases and phosphatidylinositol (Luykx et al., Reference Luykx, Boks, Terwindt, Bakker, Kahn and Ophoff2010; Brown and Tracy, Reference Brown and Tracy2013). It is not clear how these mechanisms relate to higher order properties thought to be involved in BD etiology like neuronal function, chronobiology, and brain structure. Examining lithium mechanisms at high-biological resolution is therefore not only crucial for understanding the high rates of non-response and non-adherence to prophylactic lithium treatment in BD patients but also for understanding BD etiology itself.
The widespread but subtle gene expression changes observed in lithium users are in line with lithium's broad scope of physiological effects (Roux and Dosseto, Reference Roux and Dosseto2017) and with the complex genetic architecture of BD (Fromer et al., Reference Fromer, Roussos, Sieberts, Johnson, Kavanagh, Perumal, Ruderfer, Oh, Topol, Shah, Klei, Kramer, Pinto, Gumus, Cicek, Dang, Browne, Lu, Xie, Readhead, Stahl, Xiao, Parvizi, Hamamsy, Fullard, Wang, Mahajan, Derry, Dudley, Hemby, Logsdon, Talbot, Raj, Bennett, De Jager, Zhu, Zhang, Sullivan, Chess, Purcell, Shinobu, Mangravite, Toyoshiba, Gur, Hahn, Lewis, Haroutunian, Peters, Lipska, Buxbaum, Schadt, Hirai, Roeder, Brennand, Katsanis, Domenici, Devlin and Sklar2016). These genes were enriched for functional annotations related to transmembrane, cell signaling, protein kinase, and immunity. These pathways have been implicated in previous BD transcriptome studies (Xiao et al., Reference Xiao, Camarillo, Ping, Arana, Zhao, Thompson, Xu, Su, Fan, Ordonez, Wang, Mao, Zhang, Cruz, Escamilla, Li and Xu2014; Cruceanu et al., Reference Cruceanu, Tan, Rogic, Lopez, Torres-Platas, Gigek, Alda, Rouleau, Pavlidis and Turecki2015; Mertens et al., Reference Mertens, Wang, Kim, Yu, Pham, Yang, Zheng, Diffenderfer, Zhang, Soltani, Eames, Schafer, Boyer, Marchetto, Nurnberger, Calabrese, Odegaard, Mccarthy, Zandi, Alda, Nievergelt, Mi, Brennand, Kelsoe, Gage and Yao2015; Kittel-Schneider et al., Reference Kittel-Schneider, Hilscher, Scholz, Weber, Grunewald, Schwarz, Chiocchetti and Reif2017; Pacifico and Davis, Reference Pacifico and Davis2017) and are known targets of lithium (Alda, Reference Alda2015; Maddu and Raghavendra, Reference Maddu and Raghavendra2015). The elevated levels of neutrophil proportions we observed are in line with lithium-induced neutrophilia, which has been described since the medication's early use in psychiatry (Maddu and Raghavendra, Reference Maddu and Raghavendra2015). Lithium is thought to induce neutrophilia through a complex pathway involving GSK3 and immune-related transcription factors and genes (Kast, Reference Kast2008). Increased levels of neutrophils are typically associated with anti-inflammatory or infection-fighting immune responses (Rosales et al., Reference Rosales, Demaurex, Lowell and Uribe-Querol2016). Whether these immunity-related mechanisms play a role in the mood stabilizing effects of lithium remains to be determined. Immune components of psychiatric illness including BD (Rosenblat and McIntyre, Reference Rosenblat and Mcintyre2017) have long been recognized, but it remains unclear if they represent a causal pathway, a property of the disease state, or a consequence of environmental factors like body mass index or smoking. These results contribute to the understanding of the genomics of lithium action, which may be essential for the future of personalized psychiatric medicine for patients with BD. Future studies with larger sample sizes and independent replication datasets will be needed to confirm our findings, and whether these genes and pathways play a role in the mood-stabilizing mechanisms of lithium remains to be determined.
The lack of enrichment of the genetic signal from common alleles associated with BD, schizophrenia, or self-reported depression suggests that genes transcriptionally associated with lithium treatment in peripheral blood most likely represent secondary effects of treatment that are independent from disease susceptibility. The lack of genetic enrichment could also indicate that our gene expression study is underpowered for this purpose, or that the transcriptomic mechanisms of genetic risk for BD are not present in whole blood. In addition, the currently available GWAS may still be underpowered thereby impacting our ability to detect a significant enrichment. With the expected rapidly increasing sample sizes of these GWAS studies we will be able to test this hypothesis more fully in the near future. We did explore the opportunity to examine enrichment of genetic susceptibility of lithium response, but because this phenotype has a SNP-based heritability not different from zero, this specific analysis is not meaningful. In this regard, it is important to distinguish between lithium use, the phenotype we used in our study, and lithium response. Self-reported answers to a lithium questionnaire by participants in our study show that the majority of subjects being treated with lithium had a positive response to the treatment and the majority of non-users have been treated with lithium in the past (online Supplementary Methods). We therefore consider that the lithium use phenotype partially captures lithium response, but disentangling the complex interplay between these phenotypes is an avenue for further exploration.
Because lithium use is a trait only present in BD subjects and therefore confounded with BD diagnosis, it is a confounder by indication and likely eliminated most of the observable BD effects. Our results highlight the importance of correcting for cell type composition as well as medication use in BD transcriptome studies. A lithium-naive study design is warranted to optimize BD transcriptomic signal that is independent of lithium use. Nevertheless, investigating the BD transcriptome in whole blood remains valuable for the following reasons. It is an accessible tissue, it has the potential for biomarker discovery, and it can be used in longitudinal study designs, which are appealing due to the episodic nature of BD. It may also be a choice tissue to observe the suggested immune component of BD etiology. In addition, peripheral tissues such as blood partially recapitulate gene expression signatures of the brain (Cai et al., Reference Cai, Langfelder, Fuller, Oldham, Luo, Van Den Berg, Ophoff and Horvath2010), and compared to post-mortem tissues are less subject to poor quality due to rapid degradation upon death (Popova et al., Reference Popova, Mennerich, Weith and Quast2008). However, studies involving post-mortem tissue or in vitro neuronal cells will still be needed to determine the therapeutic effect of lithium on BD-associated brain-related function.
We recognize some limitations of this study. For example, designed to investigate the gene expression signatures of bipolar disorder, the setup was not intended to measure the effects of lithium on gene expression. Lithium use was overlapped partially with BD diagnosis, which we attempted to disentangle, and likely with the use of other medications and comorbid disorders. We did not consider comorbid disorders in addition to BD in the current analysis, although it is possible that other psychiatric and somatic comorbidities might influence gene expression in this sample. Additionally, the lithium use phenotype used was self-reported and may be biased as such, and did not necessarily reflect a therapeutic response to the medication. A longitudinal design in a more homogenous, well-characterized cohort of lithium-naïve cases may be better suited to capture the gene expression signatures of lithium.
In summary, our findings suggest that there are minimal bipolar disorder-associated gene expression changes in whole blood independent of medication use and underline the importance of accounting for such confounders in psychiatric genomic studies. While limited in their ability to uncover mechanisms associated with genetic risk, blood-based transcriptome analyses of BD may still be informative with larger sample sizes and careful designs. Lastly, our findings provide molecular insights into the potential therapeutic actions of lithium, including cell signaling and immunity-related functions.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291719002745.
Data
Gene expression data is available under the GEO accession number GSE124326.
Acknowledgements
We thank the study subjects for their willingness to provide specimens and clinical data. We thank GWAS research participants and researchers involved in making summary statistics available, including 23andMe research participants and employees, for making this work possible. We thank Drs. Francis McMahon and Thomas Schulze for correspondence regarding lithium-response GWAS.
Financial support
This work was supported by National Institutes of Health grants R01MH090553 (RAO), R01MH114152 (RAO), 25602-NARSAD Distinguished Investigator Award (RAO), K99MH116115 (LMOL), and T32HG002536 (CEK).
Conflict of interest
None.
Ethical standards
The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees of human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.
