# Focal adhesion is associated with lithium response in bipolar disorder: evidence from a network-based multi-omics analysis

Vipavee Niemsiri, Sara Brin Rosenthal, Caroline M Nievergelt, Adam X

Maihofer, Maria C Marchetto, Renata Santos, Tatyana Shekhtman, Ney

Alliey-Rodriguez, Amit Anand, Yokesh Balaraman, et al.

► To cite this version:

Vipavee Niemsiri, Sara Brin Rosenthal, Caroline M Nievergelt, Adam X Maihofer, Maria C Marchetto, et al.. Focal adhesion is associated with lithium response in bipolar disorder: evidence from a network-based multi-omics analysis. Molecular Psychiatry, 2023, 10.1038/s41380-022-01909-9. inserm-04118355

HAL Id: inserm-04118355

https://inserm.hal.science/inserm-04118355

Submitted on 6 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE **OPEN**

<figure><img src="image_2.png" /></figure>

Check for updates

# Focal adhesion is associated with lithium response in bipolar disorder: evidence from a network-based multi-omics analysis

Vipavee Niemsiri<fnref n="1" />, Sara Brin Rosenthal², Caroline M. Nievergelt ID <fnref n="1" />, Adam X. Maihofer<fnref n="1" />, Maria C. Marchetto³,⁴, Renata Santos ID ⁴,<fnref n="5" />, Tatyana Shekhtman<fnref n="1" />, Ney Alliey-Rodriguez ID <fnref n="6" />,⁷, Amit Anand⁸,⁹, Yokesh Balaraman¹⁰, Wade H. Berrettini¹¹, Holli Bertram¹², Katherine E. Burdick⁹, Joseph R. Calabrese¹³,¹⁴, Cynthia V. Calkin ID ¹⁵, Carla Conroy¹³,¹⁴, William H. Coryell¹⁶, Anna DeModena¹⁷, Lisa T. Eyler<fnref n="1" />, Scott Feeder¹⁸, Carrie Fisher¹⁰, Nicole Frazier¹², Mark A. Frye ID ¹⁸, Keming Gao¹³,¹⁴, Julie Garnham ID ¹⁹, Elliot S. Gershon<fnref n="6" />, Fernando S. Goes ID ²⁰, Toyomi Goto¹³, Gloria J. Harrington¹², Petter Jakobsen²¹, Masoud Kamali ID ⁸,¹², Marisa Kelly ID ¹², Susan G. Leckband¹⁷, Falk W. Lohoff ID ¹¹, Michael J. McCarthy<fnref n="1" />, Melvin G. McInnis ID ¹², David Craig²², Caitlin E. Millett⁹, Francis Mondimore²⁰, Gunnar Morken<fnref n="23" />, John I. Nurnberger¹⁰,<fnref n="24" />, Claire O' Donovan¹⁹, Ketil J. Øedegaard²¹, Kelly Ryan¹², Martha Schinagle¹³, Paul D. Shilling ID <fnref n="1" />, Claire Slaney¹⁹, Emma K. Stapp ID ²⁵, Andrea Stautland ID ²⁶, Bruce Tarwater¹⁶, Peter P. Zandi²⁰, Martin Alda ID ¹⁹,²⁷, Kathleen M. Fisch²,²⁸, Fred H. Gage ID ⁴ and John R. Kelsoe<fnref n="1" />,<fnref n="29" />

© The Author(s) 2023

Lithium (Li) is one of the most effective drugs for treating bipolar disorder (BD), however, there is presently no way to predict response to guide treatment. The aim of this study is to identify functional genes and pathways that distinguish BD Li responders (LR) from BD Li non-responders (NR). An initial Pharmacogenomics of Bipolar Disorder study (PGBD) GWAS of lithium response did not provide any significant results. As a result, we then employed network-based integrative analysis of transcriptomic and genomic data. In transcriptomic study of iPSC-derived neurons, 41 significantly differentially expressed (DE) genes were identified in LR vs NR regardless of lithium exposure. In the PGBD, post-GWAS gene prioritization using the GWA-boosting (GWAB) approach identified 1119 candidate genes. Following DE-derived network propagation, there was a highly significant overlap of genes between the top 500- and top 2000-proximal gene networks and the GWAB gene list ($P_{\text{hypergeometric}} = 1.28E-09$ and 4.10E-18, respectively). Functional enrichment analyses of the top 500 proximal network genes identified focal adhesion and the extracellular matrix (ECM) as the most significant functions. Our findings suggest that the difference between LR and NR was a much greater effect than that of lithium. The direct impact of dysregulation of focal adhesion on axon guidance and neuronal circuits could underpin mechanisms of response to lithium, as well as underlying BD. It also highlights the power of integrative multi-omics analysis of transcriptomic and genomic profiling to gain molecular insights into lithium response in BD.

Molecular Psychiatry; https://doi.org/10.1038/s41380-022-01909-9

## INTRODUCTION

Bipolar disorder (BD) is a major psychiatric disorder characterized by recurrent episodes of mania and depression, and a high risk of suicide. Approximately 50% of BD patients suffer

psychosis, and, if left untreated, up to about 17% will complete suicide [1]. Though effective treatments exist, little is understood regarding etiology to guide clinical drug selection or drug design.

<footnote n="1">¹Department of Psychiatry, University of California, San Diego, La Jolla, CA, USA. ²Center for Computational Biology and Bioinformatics, University of California, San Diego, La Jolla, CA, USA. ³Department of Anthropology, University of California, San Diego, La Jolla, CA, USA. ⁴Laboratory of Genetics, Salk Institute for Biological Studies, La Jolla, CA, USA.</footnote>

<footnote n="5">⁵University of Paris, Institute of Psychiatry and Neuroscience of Paris (IPNP), INSERM U1261266, Laboratory of Dynamics of Neuronal Structure in Health and Disease, Paris, France.</footnote>

<footnote n="6">⁶Department of Psychiatry and Behavioral Neuroscience, University of Chicago, Chicago, IL, USA. ⁷Department of Psychiatry and Behavioral Neuroscience, Northwestern University, Chicago, IL, USA. ⁸Department of Psychiatry, Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA. ⁹Department of Psychiatry, Brigham and Women's Hospital, Harvard Medical School, Boston, MA, USA. ¹⁰Department of Psychiatry, Indiana University School of Medicine, Indianapolis, IN, USA. ¹¹Department of Psychiatry, University of Pennsylvania, Philadelphia, PA, USA. ¹²Department of Psychiatry, University of Michigan, Ann Arbor, MI, USA. ¹³Mood Disorders Program, Case Western Reserve University School of Medicine, Cleveland, OH, USA. ¹⁴Mood Disorders Program, University Hospitals Cleveland Medical Center, Cleveland, OH, USA. ¹⁵Department of Psychiatry and Medical Neuroscience, Dalhousie University, Halifax, NS, Canada. ¹⁶Department of Psychiatry, University of Iowa, Iowa City, IA, USA. ¹⁷Psychiatry Service, VA San Diego Healthcare System, San Diego, CA, USA. ¹⁸Department of Psychiatry, The Mayo Clinic, Rochester, MN, USA. ¹⁹Department of Psychiatry, Dalhousie University, Halifax, NS, Canada. ²⁰Department of Psychiatry and Behavioral Sciences, Johns Hopkins University, Baltimore, MD, USA. ²¹Norment, Division of Psychiatry, Haukeland University Hospital and Department of Clinical medicine, University of Bergen, Bergen, Norway. ²²Department of Translational Genomics, University of Southern California, Los Angeles, CA, USA.</footnote>

<footnote n="23">²³Division of Mental Health Care, St Olavs University Hospital, and Department of Mental Health, Norwegian University of Science and Technology (NTNU), Trondheim, Norway.</footnote>

<footnote n="24">²⁴Medical and Molecular Genetics, Stark Neurosciences Research Institute, Indiana University School of Medicine, Indianapolis, IN, USA. ²⁵Department of Mental Health, Johns Hopkins Bloomberg School of Public Health, Johns Hopkins University, Baltimore, MD, USA. ²⁶Department of Clinical Medicine, University of Bergen, Bergen, Norway. ²⁷National Institute of Mental Health, Klecany, Czech Republic. ²⁸Department of Obstetrics, Gynecology & Reproductive Sciences, University of California, San Diego, La Jolla, CA, USA.</footnote>

<footnote n="29">²⁹Institute for Genomic Medicine, University of California, San Diego, La Jolla, CA, USA. ³⁰email: jkelsoe@health.ucsd.edu</footnote>
Lithium (Li) is the first and remains the best mood-stabilizing medication for BD [2, 3]. The mechanism of action of lithium has been studied for over six decades and multiple effects on cellular signaling processes have been identified such as: regulation of GSK3/Akt, G proteins and PKA signaling, inositol turnover, neuronal excitability (via Na⁺-K⁺-ATPase), or neurotransmitters [4]. Lithium is clinically effective in treating both mania and depression, but primarily used for prophylaxis. Approximately 30% of patients with BD enjoy a very robust response to lithium with almost complete elimination of symptoms [5, 6]. However, after onset, most patients go through multiple medication trials often over several years during which time they suffer and are at risk for suicide [1]. Many who would be excellent lithium responders never receive a trial of lithium. For these reasons, there is a great need for a predictor of lithium response to guide clinicians in prescribing lithium. Genetics may provide such a predictor as lithium responders have been shown to have a stronger family history of both BD and lithium responsiveness [5].

One of the challenges of pharmacogenomics is the labor-
intensive task of phenotyping. The gold standard for assessing
drug response is the prospective clinical trial, but sample sizes for
such studies are orders of magnitude smaller (500 vs 50 000) than
those currently successful for genome-wide association studies
(GWASs). Furthermore, there are few GWASs focusing on lithium
response, most suffering from lack of power or failure in
replication [7]. Though more data is being accumulated [8–14],
GWAS has so far not had the power to consistently detect
reproducible genes.

Human induced pluripotent stem cells (iPSCs) provide an alternative and complementary approach to identifying genes and mechanisms of lithium response. It is a revolutionary set of methods enabling access to living neurons from specific individuals and in part overcoming a major hurdle in neuropsychiatric research, the inability to readily access living brain tissue [15]. iPSC methods are now being developed to derive a variety of specific neurons, which in turn can model diseases and drug response [16, 17]. For instance, we have previously demonstrated a differential response to lithium in vitro between iPSC neurons derived from lithium responders vs non-responders [18, 19]. These data are consistent with the notion that there are two different BD sub-populations with different pathophysiologies defined by lithium response [5].

Network-based analysis is a powerful bioinformatic approach that employs the fundamental connectivity of gene networks and genetic data to derive models of disease or drug response. Such network models may implicate biological functions associated with complex traits and presumably serve specific cellular processes [20]. Network-based methods require comprehensive information but have been shown to produce promising insights in studies of various diseases, including psychiatric disorders [21, 22]. Recently, network construction has advanced by integrating multiple sources of data, (e.g., genomic, transcriptomic, and proteomic) that, in turn, improves performance [23, 24]. Therefore, a multi-omics network-based approach can be employed in psychiatric research for efficiently uncovering the mechanisms of complex traits and the goal of precision psychiatry [25, 26].

The mechanisms underlying differential response to lithium in BD remain elusive. In this study, we aim to identify genes associated with lithium response in BD by combining genetic data from a GWAS of lithium response and data from a transcriptome study of iPSC-derived neurons challenged in vitro with lithium. The idea behind this integrative analysis is to improve our power to detect genes for lithium response by combining two different independent sources of data and examining the overlap in derived networks.

To our knowledge, this is the first integrative analysis of multi-
omics data for lithium response. Here, we describe the results of
combining data between a GWAS for lithium response and 41
differentially expressed genes that were identified in iPSC-derived

BD neurons from lithium responders (LR) and lithium non-
responders (NR). GWAS genes showed a highly significant overlap
with the expression-derived network. The functional enrichment
analyses identified focal adhesion and the extracellular matrix
(ECM) as the most significant biological functions.

METHODS

The methods are summarized here and detailed in Supplementary Methods. The overall study design is illustrated in Fig. 1a.

All subjects provided written informed consent according to their institution's approved procedures. Subjects for the GWAS were recruited as part of two studies of lithium response and BD: the multi-site Pharmacogenomics of Bipolar Disorder (PGBD) study and an identical study of veterans recruited from Veterans Affairs San Diego Healthcare System (VA). The PGBD and veteran studies had a relapse prevention design where subjects were followed prospectively for up to 2.5 years [12]. All subjects had diagnoses confirmed using the Diagnostic Interview for Genetic Studies (DIGS) [27] and all were of European American (EA) ancestry.

Subjects for the iPSC studies were selected from the PGBD/VA and Halifax samples. The Halifax sample from Dalhousie University was assessed retrospectively using the Alda scale [28]. Both lithium responders and non-responders were selected from the ends of the distribution of response, either time in study (PGBD/VA sample) or Alda score (Halifax sample). Control (CT) subjects were recruited by advertising and screened for psychiatric diagnoses using the DIGS.

GWAS was conducted using the Illumina Human Psychchip on the Illumina Infinium platform (Illumina, San Diego, CA). Genotypes were called using Genome Studio (Illumina). The GWAS analysis employed the "entered_maintenance" phenotype, representing stabilization on lithium monotherapy after 4 months. Analysis began with quality control (QC), followed by alignment and imputation. Association testing used logistic regression in PLINK [29] with age, sex, and three population principal component covariates.

Following the initial analysis of all single nucleotide polymorphisms (SNPs), GWAS results were analyzed using Versatile Gene-Based Association Study (VEGAS) test [30], which obtained a single empirical *P*-value for each gene. The Genome-Wide Association Boosting (GWAB) algorithm [31] was employed in order to use network information to rank order genes.

Skin (PGBD/VA) or blood (Halifax) samples were obtained, and either fibroblasts or lymphoblasts, respectively were reprogrammed to iPSCs (two clones per subject) and differentiated to *prox1+* hippocampal dentate gyrus glutamatergic granule cells (DG) as described previously [18] and in Supplementary Methods. The five PGBD/VA bipolar and two control cell lines used in this study were identical to the cells reprogrammed as reported previously [18]. The six bipolar and four control subjects in the Halifax sample were also identical to those previously reported [19]. Evidence of pluripotency, normal karyotypes and neural induction had also been previously reported for these lines.

iPSC-derived DG-like neurons from both clinically validated lithium responders and non-responders were treated in culture both with and without lithium. Under the treatment condition, cells were treated for one week with 1 mM lithium chloride (LiCl) [18], a clinically effective blood concentration. RNA-sequencing (RNA-seq) was performed on all samples. Ribo-depleted libraries were constructed and cDNA was sequenced as paired ends on an Illumina HiSeq 2500. QC and RNA-seq analysis are detailed in Supplementary Methods. Validation of selected genes was performed using reverse transcription quantitative real-time PCR (RT-qPCR).

Initial functional analysis of gene expression was conducted using WebGestalt [32] and g:Pro lifer [33]. Network propagation was used to identify the network regions proximal to the RNA-seq differentially expressed genes [20]. A hypergeometric test was used to test for over-representation of the GWAB genes in the RNA-seq network. Clusters were derived using the Louvain graph-based clustering algorithm [34] and functional analysis of the derived clusters was performed using ToppGene [35] and g:Profiler [33].

RESULTS

Subjects

For the GWAS, out of total 728 enrolled subjects, 256 were
selected based on completeness of data, Hardy-Weinberg
a

I. Patient-Derived iPSC Neurons

<figure><img src="image_3.png" /></figure>

IV. Network Propagation Analysis

<figure><img src="image_4.png" /></figure>

V. Functional Enrichment Analysis

* Using the Louvain modularity maximization algorithm for cluster analysis
* Using the ToppGene analysis suite and the g-Profiler tool for functional analysis

The top 500-proximal network genes

<figure><img src="image_5.png" /></figure>

b

<figure><img src="image_6.png" /></figure>

c

45 Significance DE Genes in Across-Group Analysis

<figure><img src="image_7.png" /></figure>

d

Untreated (Li-) Condition

<figure><img src="image_8.png" /></figure>

Li-Treated (Li+) Condition

e

<figure><img src="image_9.png" /></figure>

equilibrium, and EA ancestry. Characteristics of the selected GWAS subjects are summarized in Supplementary Table 1. NR were significantly more likely to have rapid cycling illness. For the RNA-seq study, overall, there were no significant demographic or clinical differences between LR, NR, and CT groups selected for generation of iPSC-derived neurons (Supplementary Table 2).

GWAS analysis yielded no significant results

As shown in Supplementary Fig. 1, no SNP was genome-wide significant. The quantile-quantile plot is consistent with inadequate power. A gene-based VEGAS analysis showed similarly negative results (Supplementary Table 3).

RNA-seq analysis

Overall, a total of 13 691 genes were expressed and included in downstream analyses. Filtering out low expression transcripts and using a trimmed mean of *M*-values transformation successfully normalized the expression levels (Supplementary Fig. 2). We applied two analysis strategies: 'within-group' and 'across-group' in the RNA-seq analysis, described in detail in Supplementary Methods.

The largest difference in gene expression was between LR and NR without lithium

For the within-group analysis, gene expression in neurons treated with lithium (Li+) was compared to those without lithium (Li-). Employing a significance threshold of *P*-value < 0.05 and log₂fold-change ≥ |1|, 14 genes showed nominal significance, but none were significant (Benjamini and Hochberg (B-H) *q* ≤ 0.20) after multiple testing correction (Fig. 1b; Supplementary Fig. 3a; Supplementary Table 4). Therefore, lithium had a limited effect on gene expression.

For the across-group analysis, we examined the differential expression of genes between groups of neurons exposed to the same treatment condition. Six comparisons were made: LR vs NR, LR vs CT, and NR vs CT each treated with (Li+) and without lithium (Li-) (Supplementary Fig. 3b; Supplementary Tables 5, 6). Using the same significance criteria, a total of 45 genes (43 protein-coding) from all six comparisons were significantly differentially expressed (*P* < 0.05, log₂fold-change ≥ |1|, and B-H *q* ≤ 0.20; 'DE' genes; Fig. 1b). About half of the DE genes (25 of 45) were expressed in more than one comparison, suggesting genetic heterogeneity with small genetic effects across the comparisons. The majority of DE genes (43 of 45, 41 protein-coding) were found
**Fig. 1** Overview of the study design and summary results of RNA-sequencing analysis from iPSC-derived neurons. **a** Study design and analysis workflow. Figure 1a shows a summary of our study workflow, including major findings. The study materials and methods are detailed in Supplementary Methods. Our comprehensive findings are described in the 'Results' section. **Step I:** iPSC-derived neurons of three sample groups (6 controls [CT], 6 BD Li responders [BD.LR], and 5 BD Li non-responders [BD.NR]; all samples were obtained from male subjects) were tested with in vitro lithium. Part of step I is modified from Welham, et al. (2015) [81]. **Step II:** RNA-seq pipeline was used for analyzing transcriptomic profiles of all three sample groups under Li-treated (Li+) and untreated (Li-) conditions. The RNA-seq analyses were classified into 'within-group' and 'across-group'. Forty-one significantly DE genes were identified in LR vs NR regardless treatment conditions. **Step III:** The PGBD [12]/NA GWAS for lithium response in BD (n = 256) was used in a GWA-boosting (GWAB) analysis [31], which identified the top 5% of GWAB-prioritized genes (1119 genes). **Step IV:** Network propagation analysis [20] of the top 500-proximal gene network derived from the 41 DE genes revealed a highly significant overlap with the 1119 GWAB-prioritized genes, containing 103 genes (73 GWAB and 30 DE/seed) genes that were significantly proximal to the 41 seed genes with a hypergeometric *P* of 1.28E-09. **Step V:** Functional enrichment analysis of gene clusters and KEGG pathways in the top 500-proximal gene network identified 189 terms significantly enriched in three clusters and 37 significantly enriched KEGG pathways for lithium response in BD. **b** Distribution of RNA-seq genes identified in iPSC-derived neurons. Bar plots of serial filtering RNA-seq genes for 'within-group' (**left panel**) and 'across-group' (**right panel**) analyses. Each analysis was categorized into subgroup comparisons (total *n* = 9; see comparison details in Supplementary Methods). A total of 13 691 RNA-seq genes for each subgroup comparison were removed by serial filters. The RNA-seq genes that passed the serial filtering threshold of CPM > 1, *P*-value < 0.05, log₂FC ≥ 1], and B-H *q*-value ≤ 0.20 were considered as 'differentially expressed' (DE) genes. Each subplot shows a plot of the serial filters (y-axis) against the number of genes (x-axis). RNA-seq genes identified in subgroup comparisons with detailed expression data are presented in Supplementary Tables 4–7. **c, d** Distribution and expression profiles of 45 DE genes identified in across-group analysis of iPSC-derived neurons. A total of 45 DE genes (*P* < 0.05, log₂FC of ≥ 1], and B-H *q*-value ≤ 0.20) were identified in 'across-group' analysis: 39 genes in untreated (Li-) condition and 31 genes in Li-treated (Li+) condition. Note that out of 45 DE genes, almost all (n = 43) were found in LR vs NR comparisons: 37 genes in untreated condition and 28 genes in Li-treated condition, including 22 genes that overlapped between them. Lists of DE genes with detailed expression data under two treatment conditions are presented in Supplementary Table 8. **c** Venn diagrams of 45 DE genes comparisons among CT, BD LR, and BD NR under two treatment conditions. **d** Expression of 45 DE genes classified into subgroup comparisons under untreated (Li-) and Li-treated (Li+) conditions. Each subgroup comparison displays each graph comprising one heatmap and one bar plot. Heatmaps (**left panel**) display hierarchical clustering of gene expression levels for sets of DE genes. Gene symbols (rows) are listed and indicate direction of regulation (down-regulated, light red; up-regulated, light blue). The color scale (**top left**) represents the degree of differential expression (low, blue; high, red). The color boxes above the heatmaps represent sample groups (CT, green; BD.LR, blue; BD.NR, red). Texts (columns) below the heatmap represent samples and are colored by treatment conditions (Li-, dark blue; Li+, orange). Bar plots (**right panel**) present FC expression values (log₂ transformed) and significance of gene expression (*P* < 0.05). The color scale (**right**) represents the degree of significance in expression for each gene (low, dark purple; high, yellow), displayed as nominal *P*-values (-log transformed). **e** RNA-seq validation of selected DE genes using RT-qPCR. Bar plot and scatter plot show a FC expression comparison of RNA-seq and RT-qPCR results for the four selected DE genes (out of total 37; Fig. 1d; Supplementary Table 8a) in Li.-LR vs Li.-NR. Bar plot (**left panel**) represents relative FC expression values (y-axis) of the selected four DE genes (x-axis), measured by RNA-seq (grey bars) and RT-qPCR (yellow bars; Supplementary Fig. 4): *HEY1*, *KLF10*, *PTP4A3* were down-regulated; and *POU3F1* was up-regulated. FC values were presented as in log₂ unit. RT-qPCR were calculated using the 2⁻⁰⁻²ΔΔCt method. Error bars represent the mean ± SEM of triplicate RT-qPCR data. Scatter plot (**right panel**) shows a high correlation (Pearson *R*² = 0.971, *P* = 0.015) of log₂FC expression between RNA-seq (x-axis) and RT-qPCR (y-axis) methods. B-H Benjamini and Hochberg, CPM counts-per-million, CT controls, DE differentially expressed, FC fold-change, GWAS genome-wide association study, iPSC induced pluripotent stem cell, KEGG Kyoto Encyclopedia of Genes and Genomes pathways (https://www.genome.jp/kegg/pathway.html), LR BD Li responders, NR BD Li non-responders, RNA-seq RNA-sequencing, RT-qPCR reverse transcription quantitative real-time PCR, SEM the standard error of the mean.

in LR vs NR comparisons (Fig. 1c, d). Specifically, out of the 43 genes, the most was 37 DE genes in Li.-LR vs Li.-NR; next was 28 DE genes in Li+.LR vs Li+.NR. Moreover, 22 out of 43 genes (51.16%) were common to LR vs NR between the two treatments. Note that the direction of changes in expression of all 45 DE genes were similar among each comparison regardless of the treatment conditions. No interaction effect was significant (Supplementary Fig. 3c; Supplementary Table 7). The entire list and detailed differential expression of 45 DE genes are shown in Supplementary Table 8.

### RT-qPCR validation of selected differential expressed genes

RT-qPCR was performed for technical validation of four selected genes (*HEY1*, *KLF10*, *POU3F1*, and *PTP4A3*) from the 41 protein-coding DE genes in LR vs NR comparisons. These four genes were selected because they were identified in our previous gene expression study of lithium response [36]. RT-qPCR quantification correlated well with RNA-seq for the four genes examined (Fig. 1e; Supplementary Fig. 4).

### There was a highly significant overlap of genes between GWAB and DE-derived networks

Using standard GWAS analytic methods, we failed to identify genes significantly associated with lithium response in BD not only in SNP-based, but also in gene-based analysis (Supplementary Fig. 1; Supplementary Table 3).

In order to prioritize potential candidate genes in the GWAS dataset, we first used a gene-based VEGAS to conduct analysis of

the GWAS data. A total of 1180 genes were identified in the top 5% of the VEGAS-prioritized genes, including three genes (*APCDD1*, *DSP*, and *PTP4A3*) shared with the 41 protein-coding DE gene list. We then employed GWAB to boost the GWAS results and obtain prioritized genes utilizing network information [31]. GWAB reprioritizes GWAS genes by boosting "not quite significant" genes that are near other more significant genes in the network. This method reprioritizes genes but does not provide an association statistic for each gene, nor does it provide weights for edges indicating the functional interaction of genes. We identified a total of 1119 genes in the top 5% of the GWAB-prioritized gene list, including four genes (*FLNC*, *LEF1*, *MBP*, and *PRKD1*) that were shared with the DE gene list (n = 41). The top 5% prioritized genes obtained by VEGAS and GWAB are listed in Supplementary Table 9.

We performed network propagation of 41 DE genes using the GIANT brain interactome database [37] to construct a 500-proximal gene network with 25 020 edges (Fig. 2a-c). We further boosted the network to a 2000-proximal gene network with 157 688 edges (Fig. 2d). Out of 41 DE genes, 34 were present in both networks. For the GWAB-prioritized genes, 73 and 241 were detected in the top 500- and top 2000-gene networks, respectively. For the VEGAS-prioritized genes, 36 and 118 were detected in the top 500- and top 2000-gene networks, respectively. The details and genes included in the two networks are presented in Supplementary Table 10. Next, we examined the significance of the gene overlap between the DE-derived top 500- or top 2000-proximal gene networks and the 1119 GWAB-prioritized genes. This overlap was striking and highly significant
a

<figure><img src="image_10.png" /></figure>

<table><thead><tr><th rowspan="3">Proximal Network Genes</th><th colspan="4">Hypergeometric P</th></tr><tr><th colspan="2">GWAB</th><th colspan="2">VEGAS</th></tr><tr><th>P</th><th>-log P</th><th>P</th><th>-log P</th></tr></thead><tbody><tr><td>Top 500</td><td>1.28E-09</td><td>8.8928</td><td>0.006</td><td>2.2218</td></tr><tr><td>Top 2000</td><td>4.10E-18</td><td>17.3872</td><td>0.007</td><td>2.1549</td></tr></tbody></table>

b

Top 500 Proximal Network Genes
N = 500 genes

Top 2000 Proximal Network Genes
N = 2000 genes

<figure><img src="image_11.png" /></figure>

c Network Propagation: Top 500 Proximal Network Genes

<figure><img src="image_12.png" /></figure>

d

Network Propagation: Top 2000 Proximal Network Genes

<figure><img src="image_13.png" /></figure>

e Cluster Analysis: Top 500 Proximal Network Genes

<figure><img src="image_14.png" /></figure>

($P_{\text{hypergeometric}}$ = 1.28E-09 for the top 500 network genes; 4.10E-18 for the top 2000 network genes). In contrast, the overlap between the DE-derived gene networks and 1180 VEGAS-prioritized genes showed nominal significance ($P_{\text{hypergeometric}}$ = 0.006 for the top

500 network genes; 0.007 for the top 2000 network genes) (Fig. 2a; Supplementary Table 11). The highly significant overlap from independent data sources suggests convergence on valid common biological functions.
**Fig. 2** Summary results of the top 500- and top 2000-proximal gene networks. **a** Results of hypergeometric test between the top 500 or top 2000 proximal network genes and the top 5% prioritized genes. A line graph (**left** panel) plots the -log (hypergeometric *P*-value; y-axis) against the number of proximal network genes (x-axis). Each line corresponds to the top 5% prioritized gene sets for GWAB [31] (orange) and VEGAS [30] (green). Vertical dashed grey lines at 500 and 2000 (on the x-axis) represent the top 500 and top 2000 proximal network genes, respectively. A table (**right** panel) shows summary of hypergeometric test results. A hypergeometric *P*-value of < 0.05 shows the significance of overlap in genes between the top 500 or top 2000 proximal network genes and the top 5% prioritized genes obtained by either GWAB (1119 genes) or VEGAS (1180 genes). See the lists of the top 5% prioritized genes and the top 500 and top 2000 proximal network genes in Supplementary Tables 9 and 10, respectively. See the detailed results of the hypergeometric test in Supplementary Table 11. **b** Distribution of genes in the top 500- and top 2000-proximal gene networks derived from the 41 protein-coding DE genes. Pie charts show the number and types of genes in the proximal gene networks: the top 500-proximal gene network (**left**) and the top 2000-proximal gene network (**right**). Each pie chart represents the proportion of gene types as both numbers and percentages. Colors represent gene types (DE, orange; DE/GWAB, yellow; GWAB, green; network, purple). The list of 41 protein-coding DE genes is presented in Fig. 1c, d, and Supplementary Table 8. **c, d** Network propagation of the top 500 and top 2000 proximal network genes. The network used 41 protein-coding DE genes as the seeds and GIANT interactome [37] (Tissue-specific gene networks from HumanBase; https://hb flatironinstitute.org/about) as the background. Color scale on genes (**nodes**) represents the degree of proximity from the seeds (low/further, yellow; high/nearer, red). Shapes represent gene types. **c** Network propagation of the top 500 proximal network genes. The network resulted in 500 genes (**nodes**) connected with 25 025 edges, and contained 30 DE, 4 DE/GWAB, 69 GWAB, and 397 network genes. A subset of 103 (69 GWAB, 4 GWAB/seed, and 30 seed) genes showed statistical significance of the overlap ($P_{\text{hypergeometric}} = 1.28E-09$) in genes between the top 500-proximal gene network and top 5% GWAB-prioritized data, two independent sources of data. The interactive graph of the top 500-proximal gene network can be accessed at https://www.ndexbio.org/viewer/networks/bc0c278d-aea2-11eb-9e72-0ac135e8bacf. **d** Network propagation of the top 2000 proximal network genes. The network resulted in 2000 genes (**nodes**) connected with 157 688 edges, and contained 30 DE, 4 DE/GWAB, 237 GWAB, and 1729 network genes. Note that a subset of 271 (237 GWAB, 4 GWAB/seed, and 30 seed) genes with 3135 edges are shown here. These relevant 271 genes showed statistically significant overlap ($P_{\text{hypergeometric}} = 4.10E-18$) between the top 2000-proximal gene network and top 5% GWAB-prioritized data. The interactive graph of the top 2000-proximal gene network can be accessed at https://ndexbio.org/viewer/networks/8ddb8cc6-aea3-11eb-9e72-0ac135e8bacf. **e** Cluster enrichment analysis of the top 500-proximal gene network derived from the 41 protein-coding DE genes. Three clusters were identified based on biological functions. Colors on genes (**nodes**) indicate for each cluster (0, blue; 1, orange; 2, green). Shapes represent gene types. The basic network prior cluster enrichment analysis is displayed in Fig. 2c. See the list and distribution of genes for each cluster in Supplementary Table 10 and Fig. 3, respectively. The interactive graph of cluster enrichment for the top 500-proximal gene network can be publicly accessed at https://ndexbio.org/viewer/networks/43073550-aea3-11eb-9e72-0ac135e8bacf. For Fig. 2c-e, gene types are indicated by shapes: DE, filled triangle; DE/GWAB, filled outlined triangle; GWAB, filled outlined square; network, filled circle. DE differentially expressed, GWAB genome-wide association boosting, VEGAS versatile gene-based association study.

## Focal adhesion and the extracellular matrix were the most enriched biological functions

To gain the biological insights into the 41 DE genes, we initially explored the functional enrichment of the 41 gene set using the WebGestalt [32] and g:Profiler [33] analysis tools. The results in Supplementary Fig. 5 show that the highest ranked functions were related to 'focal adhesion' (KEGG:hsa04510; WebGestalt) and the 'extracellular matrix' (REAC:R-HSA-1474244; g:Profiler). However, no significant enrichment ($P < 0.05$ with B-H *FDR* < 1.0E-05) was observed in either test.

We further evaluated the functional enrichment of the genes in the top 500-gene network proximal to the 41 DE gene seed set. Cluster enrichment analysis of the top 500 proximal network genes identified three clusters (Fig. 2e) with 189 enriched terms (B-H $q \le 0.05$, range = 0.0498-1.01E-08), comprising clusters 0, 1, and 2, with 82, 67, and 40 functional terms, respectively. The enrichment for each of the three clusters, including the top 10 terms and gene contribution are summarized in Fig. 3. The details of each functional term for each cluster are described in Supplementary Fig. 6 and Supplementary Tables 12–14.

Overall, cluster 0 had the most enriched terms. Strikingly, almost all of the top 10 enriched terms in cluster 0 involved in the ECM, including the term with the greatest enrichment—'extracellular matrix' (GO:0031012; B-H $q = 1.01E-08$), suggesting that cluster 0 mainly represented ECM-related sub-networks (Fig. 3c). Whereas, the top ranked terms in cluster 1, such as G-protein-coupled receptors (GPCR) signaling and neurotransmission, yielded sub-networks relevant to neuronal systems. The terms with the greatest enrichment in cluster 1 (B-H $q = 1.57E-04$) were 'GPCR downstream signaling' (REAC:R-HSA-388396), 'neuronal system' (REAC:R-HSA-112316), and 'transmission across chemical synapses' (REAC:R-HSA-112315). In contrast, most of the enriched terms in cluster 2 were related to tissue-specific functions, but not strongly significant. The greatest enrichment in cluster 2 (B-H $q = 0.00541$) included 'soft tissue 1; fibroblasts' (HPA:047030_22) and 'breast; myoepithelial cells' (HPA:004030_22). Only one term (HPA:007010_22, cerebellum; Purkinje cells; B-H $q = 0.04$) was

related to neurons. Thus, the major findings of cluster analysis implicate the ECM in lithium response.

The KEGG pathway enrichment analysis of the top 500 proximal network genes revealed a total of 37 KEGG pathways that were significantly enriched (B-H $q \le 0.05$, range = 0.02409-1.05E-21) among 196 genes (Fig. 4a, b; Supplementary Fig. 7; Supplementary Table 15). The 37 significant KEGG pathways involved four main KEGG categories: cellular processes, environmental information processing, human diseases, and organismal systems, of which the one with the greatest enrichment was 'pathways in cancers' (hsa05200; B-H $q = 1.05E-21$), followed by 'focal adhesion' (hsa04510; B-H $q = 8.04E-20$) as the second highest enrichment.

Of these 37 KEGG pathways, 17 were selected a priori based on their involvement in brain function and relevant involvement in pathways modulated by lithium, such as PI3K-Akt signaling, Ras signaling, MAPK signaling, and Wnt signaling (Fig. 4a, b). These 17 relevant pathways with 130 network genes represented 26.0% of the top 500 network genes, including 7 DE, 3 DE/GWAB, and 24 GWAB genes. Of the 17 relevant KEGG pathways, the top 3 were 'focal adhesion' (hsa04510; B-H $q = 8.04E-20$), 'ECM-receptor interaction' (hsa04512; B-H $q = 1.57E-13$), and 'PI3K-Akt signaling pathway' (hsa04151; B-H $q = 9.62E-13$), respectively (Fig. 4b). These top 3 pathways were significantly overrepresented in a subset of 54 network genes (10.8% of the top 500 network genes). Note that 20 out of 54 genes overlapped among the top 3 relevant KEGG pathways (Fig. 4c, d). The top relevant KEGG pathway—'focal adhesion' with 42 network genes, is illustrated in Fig. 4e, which also shows the connection between focal adhesion and other enriched pathways.

Notably, the enrichment of the top 500-proximal gene network was consistent with the preliminary enrichment of the 41 DE gene set obtained via WebGestalt [32] and g:Profiler [33] (Supplementary Fig. 5). Likewise, Fig. 5a summarizes the UniProtKB functions [38] of the 41 DE genes, which were associated with cell/focal adhesion, the ECM, neurogenesis, including development of axon and synapse. Similarly, the UniProtKB function [38] of 22 genes
a

Enriched Terms and Gene Types in Clusters 0, 1, & 2

<figure><img src="image_15.png" /></figure>

b

Term Categories in Clusters 0, 1, & 2
Total 189 Enriched Terms

<figure><img src="image_16.png" /></figure>

c

Top 10 Enriched Terms of Clusters 0, 1, & 2
Enriched Terms and Gene Types

<figure><img src="image_17.png" /></figure>

(Fig. 5b) shared by gene sets of the top 3 terms of cluster 0 (91 genes) and the top 3 relevant KEGG pathways (54 genes) were involved in the same functions.

Hub-like genes are defined as genes that have a high degree of connectivity among pathways. Among a total of 90 genes from the top 3 of 37 KEGG pathways, seven genes (EGFR, TGFB2, CAMK2B, FGFR1, TIMP3, LAMA4, and COL1A1) showed the highest

degree of connectivity (degree > 200, range = 202–258; Supplementary Table 16), which were considered ‘hub-like’ genes. Out of these seven, four genes were also found in a 54-gene subset of the top 3 of 17 KEGG pathways (Fig. 4c, d). The top scoring gene (EGFR) was involved in focal adhesion and PI3K-Akt signaling pathway. The remaining three genes were involved in one (FGFR1, PI3K-Akt signaling) or all three pathways (LAMA4 and COL1A1).
**Fig. 3** Cluster enrichment analysis result of the top 500-proximal gene network. **a–c** Significantly enriched terms (B-H *q*-value ≤ 0.05) identified in cluster enrichment analysis of the top 500-proximal gene network derived from the 41 protein-coding DE genes. The cluster analysis identified three clusters: 0 (blue), 1 (orange), and 2 (green), containing 239, 168, and 93 genes, respectively. A total 189 functional terms were significantly enriched (B-H *q*-value ≤ 0.05) among the three clusters. Overall, cluster 0 had the most enriched terms, of these, the term with the greatest enrichment was 'extracellular matrix' (GO:0031012; B-H *q* = 1.01E-08). Details of corresponding terms for clusters 0, 1, and 2 are listed in Supplementary Tables 12, 13, and 14, respectively. The distribution of genes for each term is summarized in Supplementary Fig. 6. **a** Distribution of significantly enriched terms and gene sets in three clusters of the top 500-proximal gene network. The sunburst diagram (*central*) shows the proportion of genes involved in each of the three clusters. Each subordinate pie chart (*marginal*) shows the proportion of genes involved in significantly enriched terms for each cluster. Colors on the subordinate pie charts represent gene types. Solid and stripe patterns indicate functional (participating) and non-functional (non-participating) of the genes in each cluster, respectively. The size of each pie slice corresponds to the number of genes. The proportion is presented as both numbers and percentages. **b** Categories of significantly enriched terms in three clusters of the top 500-proximal gene network. A sunburst diagram corresponds to a total of 189 significantly enriched terms, comprising of 82, 67, and 40 terms enriched in clusters 0, 1, and 2, respectively. Among six various categories, the Gene Ontology (GO) and Human Phenotype Ontology (HP) were the majority of enrichment in cluster 0; while GO and Human Proteome Atlas (HPA) were the most terms enriched in clusters 1 and 2, respectively. The size of each slice of the pie chart corresponds to the number of enriched terms. Each proportion is displayed as both numbers and percentages. The *inner ring* shows the proportion of enriched terms for each of the three clusters, colored by clusters. The *outer ring* shows the proportion of biological term categories for each cluster, colored by categories of enriched terms. **c** Identification and summarizing features of the top 10 significantly enriched terms of clusters 0, 1, and 2. In brief, cluster 0 had the terms with the greatest enrichment, and almost all the top 10 terms mainly involved in the extracellular matrix (ECM). The top 10 significantly enriched terms in cluster 1 tended to involve in neuronal systems and synapses. Whereas, those in cluster 2 showed various tissue-specific association other than neuronal tissues. Bar plot (**left panel**) shows the number and types of genes in each of the top 10 enriched terms for each cluster. The *x-axis* represents the number of genes (gene counts) and gene types distributing to each corresponding term. The *y-axis* corresponds to significantly enriched terms grouped by clusters, shown in text colors. Colors on bars indicate gene types. Bubble plot (**right panel**) shows the statistics of enrichment for each of the enriched terms. The size of bubbles represents the number of genes (gene counts) in each corresponding term. The bubble coordinate on the *x-axis* represents the degree of enrichment for DE genes known as 'DE gene rich factor'. The DE gene rich factor is the ratio of DE genes in each term to total genes in each term. The larger rich factor represents the greater enrichment. The color scale indicates the degree of significance (B-H *q*-value ≤ 0.05) in enrichment for each corresponding term (low, dark purple; high, yellow). The significance of enrichment is presented as the -log transformed B-H *q*-value. For Fig. 3a, *c*, gene types are indicated by colors: DE, orange; DE/GWAB, light green (**a**) or yellow (**c**); GWAB, green; network, purple. B-H Benjamini and Hochberg, DE differentially expressed, GWAB genome-wide association boosting.

Noticeably, all seven hub-like genes were not DE; but rather GWAB or network genes, since only 34 DE genes (6.80%) were part of the top 500-proximal gene network.

## DISCUSSION

In this study, we have combined GWAS and transcriptomic data to improve our overall power to detect genes and biological functions associated with lithium response. Analysis of the RNA-seq data indicated that the largest overall difference was between BD LR and BD NR, particularly in the absence of lithium, implying that inter-individual differences are a larger effect than the effect of lithium. In addition, the effect of lithium was similar between LR and NR subjects. We further demonstrated a highly significant overlap of the independent gene sets—the DE gene-propagated network and GWAB-prioritized genes ($P_{\text{hypergeometric}}$ = 1.28E-09 for the top 500-proximal gene network; 4.10E-18 for the top 2000-proximal gene network), indicating that GWAB genes are in the same gene neighborhoods as the 41 DE genes. Functional enrichment analyses of the 500-gene network identified more than 180 functions and revealed 'focal adhesion' (KEGG), 'ECM-related functions' (cluster 0 and KEGG), and 'PI3K-Akt signaling' (KEGG) as the functions influencing lithium response. PI3K-Akt signaling has long been implicated in lithium's action, though the role of focal adhesion and the ECM in lithium response in BD are relatively novel results.

Several limitations temper the interpretation of our results. Most prominent is the small sample size of each iPSC-derived neuron group (6 LR; 5 NR; 6 controls) and the GWAS (n = 256). Not surprisingly, no significant SNPs in the GWAS were detected due to low power. Only EA subjects were included, limiting generalization to other racial/ethnic groups. Our iPSC model using very young neurons may not reflect the behavior of mature neurons interacting with other surrounding cells in vivo from BD patients. Furthermore, our iPSC model focuses on only one neuronal type, the *prox1+* DG hippocampal neuron. This type of neurons was selected as it previously demonstrated lithium rescue of hyperexcitability specifically in BD LR [18, 19]. For iPSC generation, we used different tissues of origin, which

were included as a covariate in our RNA-seq analysis. In addition, lithium has been shown to facilitate reprogramming [39], possibly leading to confounds. Another concern is about "epigenetic memory", in which epigenetic marks and aspects of cell identity may be preserved through reprogramming. This issue is difficult to address, as it is not well understood. Nonetheless, it has been shown that the epigenetic effect on iPSCs diminishes during differentiation and stabilizes over time with increasing passage number [40, 41]. Also, genetic variation in iPSCs seems to be mainly driven by genetic background in individual donors rather by the induced reprogramming process [42–44]. Moreover, though incorporating network information substantially increased the power, this must be interpreted in the light of the limitations of the datasets which may not be complete. Extended multi-omics data would provide further support. Lastly, our results are all correlative and the biology not rigorously tested.

There are limited reports with which to compare our results. Recently, studies of lithium response in BD have examined gene expression but utilized different tissues or different comparisons (BD vs controls; BD LR vs BD NR; Li-exposure vs non-exposure). The majority of studies utilized lymphocytes [36, 45–47], whose expression pattern is of limited relevance to brain; while few studies utilized either iPSC-derived neurons [18, 48] or post-mortem brains [49]. In the current study, six of our 41 DE genes from LR vs NR comparisons were also reported in previous transcriptomic studies of lithium response in BD. For five genes (*HEY1*, *KLF10*, *PDGFA*, *POU3F1*, and *PTP4A3*), expression in lymphoblasts was shown to be modulated by lithium [36]. Another gene (*LEF1*) was reported to be responsible for resistance to lithium in BD NR [48]. Surprisingly, only *PDGFA* appeared to be involved in both focal adhesion and the ECM (Fig. 5a). Our study failed to detect genes previously reported for lithium response (i.e., *GADL1* [10], *GRIA2* [8], *SESTD1* [13], IncRNAs [11], and HLA antigen genes [11]). Nonetheless, we performed a *post-hoc* comparison of our 37 DE genes (Li-LR vs Li-NR) with genes identified in our PGBD/VA and two other GWA studies: the ConLiGen consortium GWAS of lithium response [9] (n = 2563)
a

Top 500 Proximal Network Genes
500 genes

<figure><img src="image_18.png" /></figure>

37 Enriched KEGG Pathways
196 genes

<figure><img src="image_19.png" /></figure>

17 Enriched KEGG Pathways
Relevant to BD/Neuronal System
130 genes

<figure><img src="image_20.png" /></figure>

b

KEGG Classification in KEGG Pathway Enrichment
37 Enriched KEGG Pathways

<figure><img src="image_21.png" /></figure>

Enriched KEGG Pathways

<figure><img src="image_22.png" /></figure>

Enriched KEGG Pathways and Gene Types

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

Statistics of Enrichment

<figure><img src="image_25.png" /></figure>

c

Network Propagation: Top 3 Relevant Enriched KEGG Pathways

<figure><img src="image_26.png" /></figure>

d

Top 3 Relevant Enriched KEGG Pathways

<figure><img src="image_27.png" /></figure>

e

<figure><img src="image_28.png" /></figure>

* DE genes
* DE/GWAB genes
* GWAB genes
* Network genes
* Hub-like genes

and the Psychiatric Genomics Consortium Bipolar Disorder Working Group (PGC-BD) GWAS of BD [50] (n = 51 710). Of the 35 DE genes for which GWAS results were available, 14 genes showed significant SNP association after Bonferroni correction (12 in our PGBD/VA study; DSP and *LMX1B* in the ConLiGen study; *ADAMTS14*

and *FOXO6* in the PGC-BD study; Supplementary Table 17). *LMX1B* and *ADAMTS14* were significant in two studies. Of these 14 significantly associated genes, 12 were part of the top 500-proximal gene network and achieved statistically significant functional enrichment. Taken together, our *post-hoc* findings
**Fig. 4** KEGG pathway enrichment analysis result of the top 500-proximal gene network. **a**, **b** Significantly enriched KEGG pathways (B-H *q*-value ≤ 0.05) identified in KEGG pathway enrichment analysis of the top 500-proximal gene network derived from the 41 protein-coding DE genes. The KEGG pathway analysis (containing at least one seed gene and one GWAB gene) identified a total of 37 KEGG pathways, including a subset of 17 KEGG pathways relevant to BD/neuronal system, which were significantly enriched (B-H *q*-value ≤ 0.05) in the top 500-proximal gene network. Details for each pathway are listed in Supplementary Table 15. The distribution of genes for each pathway is summarized in Supplementary Fig. 7. **a** Summary of significantly enriched KEGG pathways and gene sets in the top 500-proximal gene network. Pie charts display the number of genes from each source that comprise the top 500-proximal gene network (**left**), a set of 196 genes of the 37 enriched KEGG pathways (**middle**), and a set of 130 genes of the 17 relevant enriched KEGG pathways (**right**). Each pie represents the proportion of gene types as both numbers and percentages. Colors represent gene types. Solid and stripe patterns indicate function (participating) and non-function (non-participating) of genes, respectively, in each gene set of either 37 or 17 enriched KEGG pathways. **b** KEGG pathway classification and enrichment statistics of significantly enriched KEGG pathways for the top 500-proximal gene network. A total of 37 enriched KEGG pathways is on the **left** panel, and a subset of 17 relevant enriched KEGG pathways is on the **right** panel. These pathways were classified into four classes: cellular processes (yellow), environmental information processing (blue), human diseases (grey), and organismal systems (green). Of the 37 enriched KEGG pathways, the most significant enrichment was 'pathways in cancer' (hsa05200; B-H *q* = 1.05E-21; see red box). The second most significant enrichment was 'focal adhesion' (hsa04510; B-H *q* = 8.04E-20; see red box), which was also the most significant one among the 17 relevant enriched KEGG pathways. Of the 17 relevant enriched KEGG pathways, the top 3 significant enrichment were 'focal adhesion' (as mentioned), 'ECM-receptor interaction' (hsa04512; B-H *q* = 1.57E-13), and 'PI3K-Akt signaling pathway' (hsa04151; B-H *q* = 9.62E-13), respectively. Pie charts (**top panel**) show the distribution of KEGG pathway classification. Each pie displays the proportion of KEGG classes as both numbers and percentages, colored by KEGG classes. Bar and bubble plots (**bottom panel**) show the distribution of genes and/or the statistical enrichment for each of the enriched KEGG pathways, respectively. Bar plots show the distribution of genes for each KEGG pathway. The x-axis represents the number of genes contributing to each corresponding pathway. The y-axis corresponds to significantly enriched KEGG pathways categorized into KEGG classes and shown in text colors. For the 37 enriched KEGG pathways (**bottom left panel**), the color scale on bars indicates the degree of significance (B-H *q*-value ≤ 0.05) of enrichment (low, blue; high, red). An asterisk (*) on bars specifies the KEGG pathways that are relevant to BD/neuronal system. Whereas, for a subset of the 17 relevant enriched KEGG pathways ranked by the degree of significance in enrichment (**bottom middle panel**), colors on bars indicate gene types. Bubble plot (**bottom right panel**) shows the statistical enrichment of a subset 17 relevant enriched KEGG pathways. The size of bubbles represents the number of genes (gene counts) in each corresponding KEGG pathway. The bubble coordinate on the x-axis represents the degree of enrichment for DE genes known as 'DE gene rich factor'. The DE gene rich factor is the ratio of DE genes in each pathway to total genes in each pathway. The larger rich factor represents the greater enrichment. The color scale indicates the degree of significance (B-H *q*-value ≤ 0.05) in enrichment for each corresponding pathway (low, dark purple; high, yellow). The significance of enrichment is presented as the -log transformed B-H *q*-value. **c-e** Identification and summarizing features of the top 3 significantly enriched KEGG pathways relevant to BD/neuronal system (B-H *q*-value ≤ 0.05). The top 3 of the 17 relevant significantly enriched KEGG pathways are 'top 1' - focal adhesion (hsa04510; light blue), 'top 2' - ECM-receptor interaction (hsa04512; light yellow), and 'top 3' - PI3K-Akt signaling pathway (hsa04151; light pink), respectively. A subset of 54 genes (out of the top 500 proximal network genes) are present in the top 3 relevant enriched KEGG pathways. **c** Sub-networks among the 54 genes in the top 3 significantly enriched KEGG pathways relevant to BD/neuronal system. A propagation-based network illustrates sub-networks among a set of 54 genes (out of the top 500 proximal network genes) present in the top 3 relevant enriched KEGG pathways. This network included 54 genes (**nodes**) connected with 506 edges, containing 3 DE, 1 DE/GWAB, 12 GWAB, and 38 network genes. Shapes represent gene types (DE, filled triangle; DE/GWAB, filled outlined triangle; GWAB, filled outlined square; network, filled circle). Color scale indicates the degree of connectivity calculated from the top 500 proximal network genes (low degree/poorly connected, yellow; high degree/ highly connected, dark purple). Background colors are specified for the top 3 relevant enriched KEGG pathways. The connectivity of the 54 proximal network genes of the top 3 relevant enriched KEGG pathways is presented in Supplementary Table 16. The interactive graph exclusively displayed the connectivity degree of the top 500 proximal network genes can be accessed at https://www.ndxbio.org/viewer/ networks/16809bed-82ed-11ec-b3be-0ac135e8bacf. **d** Venn diagram of the 54 genes in the top 3 significantly enriched KEGG pathways relevant to BD/neuronal system. A total of 54 gene names are listed and categorized into gene types. There are 42, 23, and 44 genes involved in the top 1, top 2, and top 3 relevant enriched KEGG pathways, respectively, including 20 overlapping genes and four 'hub-like' genes. Venn diagram colors indicate the top 3 relevant enriched KEGG pathways. **e** Pathview of focal adhesion (hsa04510), the top significantly enriched KEGG pathway relevant to BD/neuronal system. The top relevant significantly enriched KEGG pathway, 'focal adhesion' (B-H *q*-value = 8.04E-20), contained a total of 42 genes (including three 'hub-like' genes), which are listed and categorized into gene types. This figure also illustrates several KEGG pathways connected with focal adhesion. Those relevant significantly enriched pathways (B-H *q*-value ≤ 0.05) identified in our study are highlighted in blue. The color scale indicates the degree of network proximity to seed genes (more proximal/farther, yellow; less proximal/nearer, red). Pathview map created by the R-based Pathview [82] software from Bioconductor (https:// www.bioconductor.org/). For Fig. 4a, b, d, and e, gene types are indicated by colors: DE, orange; DE/GWAB, yellow; GWAB, green; network, purple. An asterisk (*) in dark red represents the 'hub-like' genes. B-H Benjamini and Hochberg, DE differentially expressed, ECM the extracellular matrix, GWAB genome-wide association boosting, KEGG Kyoto Encyclopedia of Genes and Genomes pathways (https:// www.genome.jp/kegg/pathway.html).

indicate a striking level of correspondence between the RNA-seq and three independent GWAS results.

Focal adhesion and the ECM were novel and unexpected findings. In neurons, focal adhesions are a complex of proteins that bind multiple ECM proteins [51–53]. Integrins on the cell surface are activated by mechanical or chemical signals from the ECM, and in turn form the focal adhesion complex, which promotes actin-microtubule polymerization [52, 54, 55]. Axon guidance occurs as filopodia and lamellipodia of the growth cone detect and respond to axon guidance signals from the ECM proteins (known as guidance cues), resulting in growth-cone motility and turning [56–60]. The ECM, comprising various proteins, e.g., collagens and non-collagenous glycoproteins [61], has been shown to be a dynamic structure that provides not just structural support for neurons and glia cells but has an important

role in axon guidance and regulation of axonal growth [62, 63]. Thus, in neurons, focal adhesions and the ECM together form a "motor" that propels growth-cone movement and steering via downstream regulation of actin cytoskeleton organization [53, 55, 58–60, 64, 65].

Figure 4e is the KEGG pathway diagram for focal adhesion (hsa04510), the most enriched among the 17 pathways relevant to BD/neuronal system, annotated with gene involvement from our analyses. Among these multiple genes (n = 42), of particular interest is *EGFR*, our top hub-like gene, that appears to be essential for neuronal development, including neurite outgrowth and axonal regeneration [66]. A recent study also demonstrates that *EGFR* can modulate integrin tension and focal adhesion formation [67]. Altogether, it supports the role of *EGFR* in the focal adhesion process and the mechanism of lithium response in BD. It also
a

<figure><img src="image_29.png" /></figure>

b

<figure><img src="image_30.png" /></figure>

c

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

shows an example of biological interconnections contributing to lithium response in BD, which were successfully identified by integrative multi-omics approaches.

A variety of studies indicate a strong effect of lithium on neurons in culture and in animal models [68, 69]. Lithium has been shown at therapeutic concentrations to regulate axon morphology, i.e., promote axon growth, enlarge the growth cone, and increase neurite branching [70, 71]. Inhibition of GSK3β results in a similar phenotype of elongated axons and increased branching consistent with lithium's action being mediated by its inhibition of GSK3β [72, 73]. Lithium and two other mood stabilizing drugs, carbamazepine and valproate, all prevent growth-cone collapse, increase growth-cone area and axonal branching [74].

BD neurons are shown to have morphopathological changes in a variety of measures, e.g., reductions in number, size, density, and/or dendrite lengths [75-77], which suggests that BD neurons tend to be smaller with short dendrites and axons, as compared with normal/healthy neurons. Together with our findings, we hypothesize that BD LR inherit genetic defects in 'focal adhesion' and the ECM including the integrin-ECM interactions that cause disruption of focal adhesion function. This results in poorly branched and shortened axons with malformed growth cones that convey susceptibility to BD

(Fig. 5c). Lithium, through its actions (likely upon inhibition of GSK3β) of facilitating branching, re-arborization, and supporting the growth cone, is effective for BD LR by correcting these morphological and functional defects (Fig. 5c"). BD NR, on the other hand, have BD for reasons other than dysfunctional focal adhesion, and therefore, lithium does not rescue the relevant mechanisms and they fail to respond (Fig. 5c"). In this model, BD LR and BD NR result from different disease mechanisms as has been previously proposed [5, 18, 19, 36]. Consistent with this, another line of evidence shows a significant difference in cellular adhesion between BD LR and BD NR in induced neuron-like cells [78]. Moreover, emerging proteomic evidence reveals a fundamental role of CRMP2, a cytoskeleton modulator, in lithium response in BD neurons, via cytoskeletal dynamics, mainly at dendritic spines [79, 80]. Taken together, these thus highlight the role of the biological cell adhesion-ECM process in the underlying mechanism of lithium response in BD.

In conclusion, to our knowledge, our study is the first to report the significant role of 'focal adhesion' and the 'ECM' influencing lithium response in BD. Our results, along with the known effect of lithium on axonal growth and extension, suggest that both genetic and functional studies of lithium response and/or BD should focus efforts on the pathways of focal adhesion and the
**Fig. 5** Functional annotation of 41 DE genes from RNA-sequencing analysis, summary of gene sets—from the top 3 cluster 0 enriched terms and the top 3 17 relevant enriched KEGG pathways, and schematic illustration of a proposed model of focal adhesion in lithium response in BD. **a** Summary of functional annotation for 41 protein-coding DE genes: individual- and multi-gene functions. The bar plots show the functions of each protein-coding DE gene (n = 41; Supplementary Table 8) based on the functional enrichment analysis of the 500-proximal gene network from this study. The 41 DE genes were annotated in three clusters (0, blue; 1, orange; 2, green) (Supplementary Tables 12, 13, and 14, respectively), shown on the left bar plot, and annotated in the 37 significantly enriched KEGG pathways (17 relevant enriched KEGG, light blue; the remaining 20 KEGG, dark blue) (Supplementary Table 15), shown on the right bar plot. Note that, out of 41, seven genes (ADAMTS14, ADAMTS4, APCDD1, BTBD17, CPEN5, FOXO6, and IAH1; texts in grey) were excluded from the top 500-proximal gene network prior the functional enrichment analysis. Other six (CHST6, MBP, POU3F1, PTP4A3, RWDD2B, and TPD52) out of 34 remaining genes were not in either clusters or KEGG functions. The x-axis indicates the number of participating terms in clusters or participating pathways in KEGG functions. The y-axis indicates a list of 41 protein-coding DE gene names in alphabetical order. The tubular texts (left box) show the single and co-occurrence gene functional modules based on the UniProt Knowledgebase evidence [38] (UniProtKB; https://www.uniprot.org/). The functions for each DE gene and co-occurrence genes are shown in right and left tableau, respectively. Among the single-gene functions, the majority were 'cell/focal adhesion', 'cell migration', and 'extracellular matrix' (ECM). While the co-occurrence functions in multi-genes also tended to cluster for 'focal adhesion', 'nervous system development', including 'neurogenesis' and 'axon development'. **b** Distribution and summary of gene sets in the top 3 cluster 0 significantly enriched terms (91 genes) and the top 3 17 relevant significantly enriched KEGG pathways (54 genes) (B-H *q*-value ≤ 0.05). Venn diagrams (top panel) show the numbers of genes participating in the top 3 cluster 0 significantly enriched terms (91 genes, tan; Supplementary Fig. 6d; Supplementary Table 12) and the top 3 17 relevant significantly enriched KEGG pathways (54 genes, grey; Fig. 4c, d; Supplementary Fig. 7; Supplementary Table 15), including 22 genes overlapping (dark red bold) between them. The size of the Venn diagrams is shown as proportion to a size of 500 genes. Doughnut charts (bottom panel) exhibit the distribution of gene types and functions for the top 3 cluster 0 significantly enriched terms (left) and the top 3 17 relevant significantly enriched KEGG pathways (right). (i) **Inner ring** represents gene types: DE (orange), DE/GWAB (yellow), GWAB (green), and network (purple) genes. (ii) **Middle ring** represents numbers of enriched terms/KEGG pathways in which each gene was participating, shown as bar scale (range from 1 to 3). (iii) **Outer ring** represents gene names that are listed in alphabetical order and colored by gene types. Overlapping genes are highlighted in dark red bold. A total of 22 overlapping genes including their single-gene functions based on evidence in the UniProt Knowledgebase [38] (UniProtKB; https://www.uniprot.org/) are listed (top right box). The functions among 22 genes were mainly involved in 'cell/focal adhesion', 'cell migration', 'ECM', 'neurogenesis', including 'axon guidance/extension', which are similar to the functions of 41 protein-coding DE genes (Fig. 5a, left box). For Fig. 5a, b, the UniProtKB was accessed on 31 May 2021. **c** Schematic of a proposed model of focal adhesion and its role in lithium response in BD. Boxes are magnified for physiological cytoskeleton dynamics (5c') and a proposed model (5c''–5c''' at the neuronal growth cones. As shown in Fig. 5c, based on our KEGG pathway enrichment result, we hypothesize that a defect in focal adhesion and the ECM, including the integrin-ECM interactions, accounts for an underlying mechanism of lithium response in BD. However, to date, the assessment of focal adhesion function in BD, including the cytoskeleton dynamics and morphology of growth cones in BD neurons, has not yet been well studied. In addition, lithium has been reported to have positive effects on axonal and growth cone morphology of normal neurons [70–74] (in grey) (top right panel). Moreover, without drug treatments, collective evidence favors the pathological morphology [75–77] in BD neurons (in purple) (bottom left panel) compared with healthy neurons (in grey) (top left panel). Altogether, we propose that, lithium treatment (bottom right panel) results in substantial improvement of axon guidance and synaptic connectivity in BD responder neurons (with blue axon), thereby correcting an underlying causative defect and resulting in successful treatment in BD responders. Whereas, BD non-responder neurons (with red axon) may have BD not only due to focal adhesion but also due to other different molecular etiologies or functional pathways not involving focal adhesion, and therefore not responsive to lithium. The proposed model at growth cones is illustrated in Fig. 5c'' (BD responder neurons) and c''' (BD non-responder neurons). 5c' Physiological cytoskeleton dynamics. The neuronal growth cone is a motile structure at the peripheral tip of the axon, which is enriched in two cytoskeletal filaments—filamentous actin (F-actin) and microtubules (MTs). The growth cone is essential in axon outgrowth, guidance and pathfinding in order to form proper synaptic connections. The interaction between actin filaments and MTs is a dynamic process resulting in protrusion or retraction along the edge of growth cones [56–60]. In brief, the protrusion of filopodia (F-actin parallel bundles) and lamellipodia (F-actin meshwork) are responsive to extrinsic attractive cues in the ECM. This induces actin polymerization and subsequent coupling of F-actin and MTs for MT polymerization and stabilization, resulting in growth cone turning and axon outgrowth. Whereas, retraction occurs in response to repulsive guidance cues in the ECM by F-actin severing, disassembly of F-actin and MTs, and MT depolymerization/destabilization. To develop proper axon outgrowth, focal adhesion is required as an adhesive linkage between the cytoskeleton and the ECM that regulates actin cytoskeletal organization (see 5c''–5c'''), thereby enabling growth cone behavior and axon guidance [64, 65]. This figure is adapted from Vitriol and Zheng (2012) [59] and Cammarata, Bearce, and Lowery (2016) [60] with permission. 5c''–5c''' Focal adhesion at growth cones of BD neurons in response to lithium in the proposed model. Focal adhesion plays a pivotal role in linking the cytoskeleton to ECM via an interaction between transmembrane integrin receptors and ECM proteins (i.e., laminins, collagens, and fibronectin) [52–55]. Focal adhesions and integrin-ECM interactions occur in neuronal cells, mainly at the tip of growth cones, named point contact (PC). At PC, several signaling and adhesion molecules such as FAK-Src, paxillin, vinculin, talin, and Cas-csk are recruited and form the focal adhesion complex [52–55], which regulates actin cytoskeleton organization downstream of the focal adhesion complex [56–60]. Thus, focal adhesion has complex, dynamic functions and is essential for regulating axon outgrowth and guidance in response to axon guidance cues (see 5c'), which leads to growth cone turning and synapse formation in neuronal cells [64, 65]. Our hypothesized model is described in detail here. In BD responder neurons (5c'') treated with lithium, focal adhesion dysfunction is corrected, resulting in strong and well-organized actin filaments and MTs. Vigorous filopodia and lamellipodia are restored at the growth cone, and then axon elongation and growth cone dynamics (protrusion and retraction, see 5c') are normalized. This suggests that focal adhesions can be rescued and preserved by lithium effects in BD responder neurons, leading to response to lithium treatment in BD responders. In contrast, for BD non-responder neurons (5c'''), when treated with lithium, focal adhesion function remains defective perhaps due to inefficient effects of lithium on those neurons, resulting in persistently weak and poorly-organized actin and microtubule structures. Unwell filopodia and lamellipodia at the growth cone cause poor neuronal axon guidance and disrupted synaptic formation. Therefore, lithium has no major effects on axonal outgrowth and growth cone morphology in BD non-responder neurons, leading to no response to lithium in BD non-responders. More details and references are provided in the 'Discussion' section. B-H Benjamini and Hochberg, Cas Crk-associated substrate, csk c-terminal Src kinase, DE differentially expressed, ECM the extracellular matrix, FAK focal adhesion kinase, GWAB genome-wide association boosting, KEGG Kyoto Encyclopedia of Genes and Genomes pathways (https://www.genome.jp/kegg/pathway.html), Src steroid receptor coactivator.

ECM as well as regulation of axonal growth/extension and synaptic connectivity. Distinguishing two distinct forms of BD would advance our understanding of disease mechanism and facilitate the development of novel therapeutics or a clinical test

for lithium response. In addition, our study demonstrates the power of applying network methods to multi-omics data, whereas advanced rigorous validation awaits additional multi-omics data and future examination not only the biology of axon guidance,

<figure><img src="image_33.png" /></figure>
axon extension and branching, but also the function of pathway-
related genes in BD neurons.

DATA AVAILABILITY

Our study was registered on the ClinicalTrials.gov with identifiers NCT01272531 (the PGBD study) and NCT00252577 (the VA study). The raw RNA-sequencing data and metadata are archived in a public repository: the NCBI's *Gene Expression Omnibus (GEO) repository* with accession number GSE205422.

CODE AVAILABILITY

Analytical source codes created for this study are openly accessed via a **GitHub repository: "Bipolar Disorder and Lithium Response: Pharmacogenomic Study"** at
https://github.com/vniems/BD-Lithium.

REFERENCES

1. Miller JN, Black DW. Bipolar disorder and suicide: a review. Curr Psychiatry Rep. 2020;22:6.

2. Balance investigators collaborators, Geddes JR, Goodwin GM, Rendell J, Azorin JM, Cipriani A, et al. Lithium plus valproate combination therapy versus mono-therapy for relapse prevention in bipolar I disorder (BALANCE): a randomised open-label trial. Lancet. 2010;375:385-95.

3. Kessing LV. Lithium as the drug of choice for maintenance treatment in bipolar disorder. Acta Psychiatr Scand. 2019;140:91-93.

4. Alda M. Lithium in the treatment of bipolar disorder: pharmacology and pharmacogenetics. Mol Psychiatry. 2015;20:661-70.

5. Grof P, Duffy A, Cavazzoni P, Grof E, Garnham J, MacDougall M, et al. Is response to prophylactic lithium a familial trait? J Clin Psychiatry. 2002;63:942-7.

6. Rybakowski JK. Response to lithium in bipolar disorder: clinical and genetic findings. ACS Chem Neurosci. 2014;5:413-21.

7. Vecera CM, Fries GR, Shahani LR, Soares JC, Machado-Vieira R. Pharmacogenomics of lithium response in bipolar disorder. Pharmaceuticals 2021;14:287.

8. Perlis RH, Smoller JW, Ferreira MA, McQuillin A, Bass N, Lawrence J, et al. A genomewide association study of response to lithium for prevention of recurrence in bipolar disorder. Am J Psychiatry. 2009;166:718-25.

9. Schulze TG, Alda M, Adli M, Akula N, Ardau R, Bui ET, et al. The International Consortium on Lithium Genetics (ConLiGen): an initiative by the NIMH and IGSLI to study the genetic basis of response to lithium treatment. Neuropsychobiology. 2010;62:72-8.

10. Chen CH, Lee CS, Lee MT, Ouyang WC, Chen CC, Chong MY, et al. Variant GADL1 and response to lithium therapy in bipolar I disorder. N. Engl J Med. 2014;370:119-28.

11. Hou L, Heilbronner U, Degenhardt F, Adli M, Akiyama K, Akula N, et al. Genetic variants associated with response to lithium treatment in bipolar disorder: a genome-wide association study. Lancet. 2016;387:1085-93.

12. Oedegaard KJ, Alda M, Anand A, Andreassen OA, Balaraman Y, Berrettini WH, et al. The Pharmacogenomics of Bipolar Disorder study (PGBD): identification of genes for lithium response in a prospective sample. BMC Psychiatry. 2016;16:129.

13. Song J, Bergen SE, Di Florio A, Karlsson R, Charney A, Ruderfer DM, et al. Genome-wide association study identifies SESTD1 as a novel risk gene for lithium-responsive bipolar disorder. Mol Psychiatry. 2016;21:1290-7.

14. Amare AT, Schubert KO, Hou L, Clark SR, Papiol S, Cearns M, et al. Association of polygenic score for major depression with response to lithium in patients with bipolar disorder. Mol Psychiatry. 2020. https://doi.org/10.1038/s41380-020-0689-5.

15. Soliman MA, Aboharb F, Zeltner N, Studer L. Pluripotent stem cells in neuropsychiatric disorders. Mol Psychiatry. 2017;22:1241-9.

16. Silva MC, Haggarty SJ. Human pluripotent stem cell-derived models and drug screening in CNS precision medicine. Ann N. Y Acad Sci. 2020;1471:18-56.

17. Hoffmann A, Sportelli V, Ziller M, Spengler D. From the Psychiatrist's couch to induced pluripotent stem cells: bipolar disease in a dish. Int J Mol Sci. 2018;19:770.

18. Mertens J, Wang QW, Kim Y, Yu DX, Pham S, Yang B, et al. Differential responses to lithium in hyperexcitable neurons from patients with bipolar disorder. Nature. 2015;527:95-9.

19. Stern S, Santos R, Marchetto MC, Mendes APD, Rouleau GA, Biesmans S, et al. Neurons derived from patients with bipolar disorder divide into intrinsically different sub-populations of neurons, predicting the patients' responsiveness to lithium. Mol Psychiatry. 2018;23:1453-65.

20. Cowen L, Ideker T, Raphael BJ, Sharan R. Network propagation: a universal amplifier of genetic associations. Nat Rev Genet. 2017;18:551-62.

21. Li J, Shi M, Ma Z, Zhao S, Euskirchen G, Ziskin J, et al. Integrated systems analysis reveals a molecular network underlying autism spectrum disorders. Mol Syst Biol. 2014;10:774.

22. Li J, Ma Z, Shi M, Malty RH, Aoki H, Minic Z, et al. Identification of Human Neuronal Protein Complexes Reveals Biochemical Activities and Convergent Mechanisms of Action in Autism Spectrum Disorders. Cell Syst. 2015;1:361-74.

23. Karczewski KJ, Snyder MP. Integrative omics for health and disease. Nat Rev Genet. 2018;19:299-310.

24. Di Nanni N, Bersanelli M, Milanesi L, Mosca E. Network diffusion promotes the integrative analysis of multiple omics. Front Genet. 2020;11:106.

25. Sullivan PF, Geschwind DH. Defining the genetic, genomic, cellular, and diagnostic architectures of psychiatric disorders. Cell. 2019;177:162-83.

26. Zuo Y, Wei D, Zhu C, Naveed O, Hong W, Yang X. Unveiling the pathogenesis of psychiatric disorders using network models. Genes 2021;12:1101.

27. Nurnberger JI Jr., Blehar MC, Kaufmann CA, York-Cooler C, Simpson SG, Harkavy-Friedman J, et al. Diagnostic interview for genetic studies. Rationale, unique features, and training. NIMH Genetics Initiative. Arch Gen Psychiatry. 1994;51:849-59. discussion 863-844.

28. Manchia M, Adli M, Akula N, Ardau R, Aubry JM, Backlund L, et al. Assessment of response to lithium maintenance treatment in bipolar disorder: a Consortium on Lithium Genetics (ConLiGen) Report. PLoS One. 2013;8:e65636.

29. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559-75.

30. Mishra A, Macgregor S. VEGAS2: Software for more flexible gene-based testing. Twin Res Hum Genet. 2015;18:86-91.

31. Shim JE, Bang C, Yang S, Lee T, Hwang S, Kim CY, et al. GWAB: a web server for the network-based boosting of human genome-wide association data. Nucleic Acids Res. 2017;45:W154-61.

32. Wang J, Vasaikar S, Shi Z, Greer M, Zhang B. WebGestalt 2017: a more comprehensive, powerful, flexible and interactive gene set enrichment analysis toolkit. Nucleic Acids Res. 2017;45:W130-7.

33. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, et al. g:Profiler—a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res. 2016;44:W83-9.

34. Blondel VD, Guillaume JL, Hendrickx JM, de Kerchove C, Lambiotte R. Local leaders in random networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2008;77: 036114.

35. Chen J, Bardes EE, Aronow BJ, Jegga AG. ToppGene Suite for gene list enrichment analysis and candidate gene prioritization. Nucleic Acids Res. 2009;37:W305-11. (Web Server issue).

36. Breen MS, White CH, Shekhtman T, Lin K, Looney D, Woelk CH, et al. Lithium-responsive genes and gene networks in bipolar disorder patient-derived lymphoblastoid cell lines. Pharmacogenomics J. 2016;16:446-53.

37. Greene CS, Krishnan A, Wong AK, Ricciotti E, Zelaya RA, Himmelstein DS, et al. Understanding multicellular function and disease with human tissue-specific networks. Nat Genet. 2015;47:569-76.

38. UniProt Consortium. UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Res. 2021;49:D480-9.

39. Wang Q, Xu X, Li J, Liu J, Gu H, Zhang R, et al. Lithium, an anti-psychotic drug, greatly enhances the generation of induced pluripotent stem cells. Cell Res. 2011;21:1424-35.

40. Ankam S, Rovini A, Baheti S, Hrstka R, Wu Y, Schmidt K, et al. DNA methylation patterns in human iPSC-derived sensory neuronal differentiation. Epigenetics. 2019;14:927-37.

41. Bar S, Benvenisty N. Epigenetic aberrations in human pluripotent stem cells. EMBO J. 2019;38:e101033.

42. Kilpinen H, Goncalves A, Leha A, Afzal V, Alasoo K, Ashford S, et al. Common genetic variation drives molecular heterogeneity in human iPSCs. Nature. 2017; 546:370-5.

43. de Boni L, Gasparoni G, Haubenreich C, Tierling S, Schmitt I, Peitz M, et al. DNA methylation alterations in iPSC- and hESC-derived neurons: potential implications for neurological disease modeling. Clin Epigenetics. 2018;10:13.

44. Volpato V, Webber C. Addressing variability in iPSC-derived models of human disease: guidelines to promote reproducibility. Dis Model Mech. 2020;13.

45. Anand A, McClintick JN, Murrell J, Karne H, Nurnberger JI, Edenberg HJ. Effects of lithium monotherapy for bipolar disorder on gene expression in peripheral lymphocytes. Mol Neuropsychiatry. 2016;2:115-23.

46. Fries GR, Colpo GD, Monroy-Jaramillo N, Zhao J, Zhao Z, Arnold JG, et al. Distinct lithium-induced gene expression effects in lymphoblastoid cell lines from patients with bipolar disorder. Eur Neuropsychopharmacol. 2017;27:1110-9.

47. Paul P, Iyer S, Nadella RK, Nayak R, Chellappa AS, Ambardar S, et al. Lithium response in bipolar disorder correlates with improved cell viability of patient-derived cell lines. Sci Rep. 2020;10:7428.

48. Santos R, Linker SB, Stern S, Mendes APD, Shokhirev MN, Erikson G, et al. Deficient LEF1 expression is associated with lithium resistance and hyperexcitability in neurons derived from bipolar disorder patients. Mol Psychiatry. 2021. https://doi.org/10.1038/s41380-020-00981-3.
49. Akkouh IA, Skrede S, Holmgren A, Ersland KM, Hansson L, Bahrami S, et al. Exploring lithium's transcriptional mechanisms of action in bipolar disorder: a multi-step study. Neuropsychopharmacology. 2020;45:947-55.

50. Stahl EA, Breen G, Forstner AJ, McQuillin A, Ripke S, Trubetskoy V, et al. Genome-wide association study identifies 30 loci associated with bipolar disorder. Nat Genet. 2019;51:793-803.

51. Robles E, Gomez TM. Focal adhesion kinase signaling at sites of integrin-mediated adhesion controls axon pathfinding. Nat Neurosci. 2006;9:1274-83.

52. Myers JP, Santiago-Medina M, Gomez TM. Regulation of axonal outgrowth and pathfinding by integrin-ECM interactions. Dev Neurobiol. 2011;71:901-23.

53. Short CA, Suarez-Zayas EA, Gomez TM. Cell adhesion and invasion mechanisms that guide developing axons. Curr Opin Neurobiol. 2016;39:77-85.

54. Chacon MR, Navarro AI, Cuesto G, del Pino I, Scott R, Morales M, et al. Focal adhesion kinase regulates actin nucleation and neuronal filopodia formation during axonal growth. Development. 2012;139:3200-10.

55. Lilja J, Ivaska J. Integrin activity in neuronal connectivity. J Cell Sci. 2018;131:jcs212803.

56. Geraldo S, Gordon-Weeks PR. Cytoskeletal dynamics in growth-cone steering. J Cell Sci. 2009;122:3595-604.

57. Lowery LA, Van, Vactor D. The trip of the tip: understanding the growth cone machinery. Nat Rev Mol Cell Biol. 2009;10:332-43.

58. Dent EW, Gupton SL, Gertler FB. The growth cone cytoskeleton in axon outgrowth and guidance. Cold Spring Harb Perspect Biol. 2011;3:a001800.

59. Vitriol EA, Zheng JQ. Growth cone travel in space and time: the cellular ensemble of cytoskeleton, adhesion, and membrane. Neuron. 2012;73:1068-81.

60. Cammarata GM, Bearce EA, Lowery LA. Cytoskeletal social networking in the growth cone: How +TIPs mediate microtubule-actin cross-linking to drive axon outgrowth and guidance. Cytoskeleton (Hoboken). 2016;73:461-76.

61. Roumazeilles L, Dokalis N, Kaulich E, Lelievre V. It is all about the support—The role of the extracellular matrix in regenerating axon guidance. Cell Adh Migr. 2018;12:87-92.

62. Barros CS, Franco SJ, Muller U. Extracellular matrix: functions in the nervous system. Cold Spring Harb Perspect Biol. 2011;3:a005108.

63. Benarroch EE. Extracellular matrix in the CNS: Dynamic structure and clinical correlations. Neurology. 2015;85:1417-27.

64. Kalil K, Dent EW. Branch management: mechanisms of axon branching in the developing vertebrate CNS. Nat Rev Neurosci. 2014;15:7-18.

65. Pacheco A, Gallo G. Actin filament-microtubule interactions in axon initiation and branching. Brain Res Bull. 2016;126:300-10.

66. Romano R, Bucci C. Role of EGFR in the Nervous System. Cells 2020;9:1887.

67. Rao TC, Ma VP, Blanchard A, Urner TM, Grandhi S, Salaita K, et al. EGFR activation attenuates the mechanical threshold for integrin tension and focal adhesion formation. J Cell Sci. 2020;133:jcs238840.

68. Salcedo-Tello P, Ortiz-Matamoros A, Arias C. GSK3 function in the brain during development, neuronal plasticity, and neurodegeneration. Int J Alzheimers Dis. 2011;2011:189728.

69. Can A, Schulze TG, Gould TD. Molecular actions and clinical pharmacogenetics of lithium therapy. Pharm Biochem Behav. 2014;123:3-16.

70. Kim WY, Zhou FQ, Zhou J, Yokota Y, Wang YM, Yoshimura T, et al. Essential roles for GSK-3s and GSK-3-primed substrates in neurotrophin-induced and hippocampal axon growth. Neuron. 2006;52:981-96.

71. Shah SM, Patel CH, Feng AS, Kollmar R. Lithium alters the morphology of neurites regenerating from cultured adult spiral ganglion neurons. Hear Res. 2013;304:137-44.

72. Owen R, Gordon-Weeks PR. Inhibition of glycogen synthase kinase 3beta in sensory neurons in culture alters filopodia dynamics and microtubule distribution in growth cones. Mol Cell Neurosci. 2003;23:626-37.

73. Kim YT, Hur EM, Snider WD, Zhou FQ. Role of GSK3 signaling in neuronal morphogenesis. Front Mol Neurosci. 2011;4:48.

74. Williams RS, Cheng L, Mudge AW, Harwood AJ. A common mechanism of action for three mood-stabilizing drugs. Nature. 2002;417:292-5.

75. Rajkowska G. Cell pathology in bipolar disorder. Bipolar Disord. 2002;4:105-16.

76. Gigante AD, Young LT, Yatham LN, Andreazza AC, Nery FG, Grinberg LT, et al. Morphometric post-mortem studies in bipolar disorder: possible association with oxidative stress and apoptosis. Int J Neuropsychopharmacol. 2011;14:1075-89.

77. Konradi C, Zimmerman EI, Yang CK, Lohmann KM, Gresch P, Pantazopoulos H, et al. Hippocampal interneurons in bipolar disorder. Arch Gen Psychiatry. 2011;68:340-50.

78. Wang JL, Shamah SM, Sun AX, Waldman ID, Haggarty SJ, Perlis RH. Label-free, live optical imaging of reprogrammed bipolar disorder patient-derived cells reveals a functional correlate of lithium responsiveness. Transl Psychiatry. 2014;4:e428.

79. Tobe BTD, Crain AM, Winquist AM, Calabrese B, Makihara H, Zhao WN, et al. Probing the lithium-response pathway in hiPSCs implicates the phosphoregulatory set-point for a cytoskeletal modulator in bipolar pathogenesis. Proc Natl Acad Sci USA. 2017;114:E4462-71.

80. Zhao WN, Tobe BTD, Udeshi ND, Xuan LL, Pernia CD, Zolg DP, et al. Discovery of suppressors of CRMP2 phosphorylation reveals compounds that mimic the

behavioral effects of lithium on amphetamine-induced hyperlocomotion. Transl Psychiatry. 2020;10:76.

81. Welham NV, Ling C, Dawson JA, Kendziorski C, Thibeault SL, Yamashita M. Microarray-based characterization of differential gene expression during vocal fold wound healing in rats. Dis Model Mech. 2015;8:311-21.

82. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics. 2013;29:1830-1.

## ACKNOWLEDGEMENTS

We thank the patients who participated. The study was primarily supported by a grant to JRK from the NIMH (U01 MH92758) as part of the Pharmacogenomics Research Network (PGRN) and a grant from the Department of Veterans Affairs. VN was supported by the UCSD Fellowship in Biological Psychiatry and Neuroscience (T32MH018399). KMF is supported by a grant from the National Center for Advancing Translational Sciences, part of the NIH, as part of support for the UCSD CTSA (UL1TR001442). FHG and MCM are supported by U19MH106434, part of the National Cooperative Reprogrammed Cell Research Groups (NCRCRG) to Study Mental Illness. FHG is also supported by grants from AHA-Allen Initiative in Brain Health and Cognitive Impairment Award made jointly through the American Heart Association and The Paul G. Allen Frontiers Group (19PABH134610000), The JPB Foundation, Bob and Mary Jane Engman, Annette C. Merle-Smith, R01 MH095741, and Lynn and Edward Streim. The Halifax group (MA, CVC, JG, CO, and CS) is supported by grants from Canadian Institutes of Health Research (#166098), ERA PerMed project PLOT-BD, Research Nova Scotia, Genome Atlantic, Nova Scotia Health Authority and Dalhousie Medical Research Foundation (Lindsay Family Fund).

## AUTHOR CONTRIBUTIONS

Designed the study: KEB, JRC, WHC, ESG, SGL, MJM, MGM, DC, JIN, KJØ, PPZ, MA, and JRK. Collected the data: VN, MCM, RS, TS, NA, AA, YB, WHB, HB, KEB, JRC, CVC, CC, WHC, AD, LTE, SF, CF, NF, MAF, KG, JG, ESG, FSG, TG, GJH, PJ, M Kamali, M Kelly, SGL, FWL, MJM, MGM, DC, CEM, FM, GM, JIN, CO, KJØ, KR, MS, PDS, CS, EKS, AS, BT, PPZ, MA, FHG, and JRK. Analysis and interpretation of the data: VN, SBR, CMN, AXM, DC, PDS, KMF, and JRK. Wrote the manuscript: VN, SBR, CMN, JRC, WHC, ESG, MJM, MGM, DC, JIN, KJØ, PDS, PPZ, MA, KMF, FHG, and JRK. Securing funding: KEB, JRC, WHC, ESG, SGL, MJM, MGM, DC, JIN, KJØ, PPZ, MA, and JRK. Overall principal investigator: JRK. All authors have read and approved the final version of the manuscript.

## COMPETING INTERESTS

JIN is an investigator for Janssen Pharmaceuticals, Inc. (US Headquarters, Raritan, NJ). No other author declares a conflict of interest.

## ADDITIONAL INFORMATION

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41380-022-01909-9.

**Correspondence** and requests for materials should be addressed to John R. Kelseo.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_34.png" /></figure>

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.