Immunogenetics of lithium response and psychiatric phenotypes in patients with bipolar disorder

The link between bipolar disorder (BP) and immune dysfunction remains controversial. While epidemiological studies have long suggested an association, recent research has found only limited evidence of such a relationship. To clarify this, we investigated the contributions of immune-relevant genetic factors to the response to lithium (Li) treatment and the clinical presentation of BP. First, we assessed the association of a large collection of immune-related genes (4,925) with Li response, defined by the Retrospective Assessment of the Lithium Response Phenotype Scale (Alda scale), and clinical characteristics in patients with BP from the International Consortium on Lithium Genetics (ConLi+Gen, N = 2,374). Second, we calculated here previously published polygenic scores (PGSs) for immune-related traits and evaluated their associations with Li response and clinical features. We found several genes associated with Li response at p < 1×10− 4 values, including HAS3, CNTNAP5 and NFIB. Network and functional enrichment analyses uncovered an overrepresentation of pathways involved in cell adhesion and intercellular communication, which appear to converge on the well-known Li-induced inhibition of GSK-3β. We also found various genes associated with BP’s age-at-onset, number of mood episodes, and presence of psychosis, substance abuse and/or suicidal ideation at the exploratory threshold. These included RTN4, XKR4, NRXN1, NRG1/3 and GRK5. Additionally, PGS analyses suggested serum FAS, ECP, TRANCE and cytokine ligands, amongst others, might represent potential circulating biomarkers of Li response and clinical presentation. Taken together, our results support the notion of a relatively weak association between immunity and clinically relevant features of BP at the genetic level.


INTRODUCTION
Bipolar disorder (BP) has been associated with some degree of immune dysfunction. Epidemiological data has linked immune-related medical comorbidities, including autoimmune and metabolic diseases, and chronic low-grade in ammation with BP. In particular, increases in pro-in ammatory cytokines are observed during affective episodes in patients with BP [1]. In addition, genomic studies have revealed weak yet signi cant genetic correlation between BP and immune-related diseases [2]. Nevertheless, as a number of these observations originated from underpowered studies [3], further investigations are required to elucidate the proposed relationships.
Lithium (Li), mainly used in the treatment of BP, is an effective pharmacological agent in the treatment of an array of psychiatric conditions [4,5]. In addition to its mood-stabilizing effects, Li shows anti-viral and immune cell regulatory properties [6,7]. The immune regulatory activity of Li has been partially attributed to the modulation of pro-in ammatory cytokines and GSK-3β. Therefore, it has been suggested that the mechanism through which Li improves symptom progression may be via anti-in ammatory effects [8,9]. The Retrospective Assessment of the Lithium Response Phenotype Scale (Alda scale) is the most widely used clinical measure of Li response. Most often, it is dichotomized such that individuals with scores ≥ 7 are classi ed as "responders" and those with scores < 7 as "non-responders" [10,11]. Using this metric, previous genetic studies have implicated human leukocyte antigen (HLA) and in ammatory cytokine genes in the response to Li treatment in BP [12,13]. Therefore, we hypothesized that single nucleotide polymorphisms (SNPs) in immune-related genes contribute, to some extent, to Li response and further, may impact speci c clinical features within BP. To test our hypothesis, we performed association studies of a comprehensive collection of immune-related genes in 2,374 patients with BP from the International Consortium on Lithium Genetics (ConLi + Gen) [14]. Additionally, we tested associations with published polygenic scores (PGSs) for immune-relevant traits.

METHODS
Since our study follows a candidate approach to selected genes, pathways and networks, a diagram summarizing the methodology employed can be found in the Supplementary Figures: Figure S1.

Study sample
The ConLi + Gen cohort has been previously described in detail [15]. Brie y, peripheral blood samples from individuals with a diagnosis of a bipolar spectrum disorder (in accordance with the criteria established in the Diagnostic and Statistical Manual of Mental Disorders -DSM-versions III or IV) that had taken Li for a minimum of six months (with no additional mood stabilizers), were collected from 2003 to 2013. The isolated DNA was genotyped in two phases. This resulted in two sample batches originally referred to as "GWAS1" and "GWAS2", comprising 1,162 and 1,401 individuals, respectively. Long-term responses to Li treatment were assessed in both sample batches using the Alda scale. Here, the A subscale rates the degree of response on a 10-point scale, and the B subscale re ects the relationship between improvement and treatment. A total score, ranging from 0-10, is obtained by subtracting the B score from the A score of these subscales. Negative scores are set to 0. Data on age-at-onset (AAO), age (at sample collection and phenotyping), sex and diagnostic subtype were available for both sample batches. Diagnoses included bipolar disorder type I and type II, schizoaffective bipolar disorder and bipolar disorder not otherwise speci ed. Additionally, information on psychiatric features, namely the number of episodes of depression, mania and hypomania, the presence of psychosis, alcohol and substance abuse, and suicidal ideation, were available for patients in the "GWAS1" batch.
The Ethics Committee at the University of Heidelberg provided central approval for the ConLi + Gen Consortium. Written informed consent from all participants was obtained according to the study protocols of each of the participating sites and their institutions. All procedures were performed in accordance with the guidelines of the Declaration of Helsinki.

Immune gene collection
A comprehensive set of immune-related genes was collated from gene lists available in the online databases MSigDB [16] and InnateDB [17]. From MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/), the following gene sets contained in the C2 "curated gene sets" collection were retrieved: M1036: Reactome-innate immune system, M1058: Reactome-adaptive immune system, M39895: WikiPathways (WP)-neuroin ammation, M39711: WP-cytokines and in ammatory response, and M39641: WPin ammatory response pathway. From InnateDB (https://www.innatedb.com/index.jsp), the curated gene lists derived from the Immunology Database and Analysis Portal (ImmPort), the Immunogenetic Related Information Source (IRIS) and the Immunome Database, were downloaded. Chromosomal locations were annotated from Ensembl using the hg19 build. Herein, the combined collection is referred to as the ImmuneSet and contained 4,925 autosomal genes to be included in association analyses.
Genotype data Schubert et al. (2021) [18] have previously described the creation of the genotype dataset used herein.
Brie y, DNA samples were originally genotyped using either Affymetrix or Illumina SNP arrays. These genotype data from multiple cohorts were separately imputed using the 1000 Genomes Project reference panel phase 3 v5. Each imputed dataset underwent a basic quality control (QC) step to keep variants with minor allele frequency (MAF) > 0.01, Hardy-Weinberg equilibrium p-value (HWE) ≥ 1x10 − 6 and imputation quality score (Rsq) ≥ 0.6. Genotype calls were derived from the imputed dosage scores and all datasets were merged by retaining only common sets of SNPs. To update this dataset and obtain a higher number of good quality variants, we re-imputed the genotype data via the Michigan Imputation Server [19] using the Haplotype Reference Consortium (HRC) panel for European ancestry. The re-imputed genotypes underwent a QC step to keep variants with Rsq ≥ 0.8, MAF ≥ 0.01 and HWE ≥ 1x10 − 6 . Additionally, individuals were removed if they failed the heterozygosity test and/or showed relatedness, according to the tests performed using the plinkQC R package [20]. In the latter case, one individual from each pair of related individuals (PI-HAT > 0.25) was removed. For analysis of the ImmuneSet, SNPs within each gene's boundaries (± 0 kb) were retained. The nal ImmuneSet genotype datasets contained 701,031 SNPs from 1,024 and 1,350 individuals in "GWAS1" and "GWAS2", respectively.

Polygenic Scores
A set of 32 published PGSs available at the PGS Catalog [21] were used to approximate markers of in ammation and immune-related phenotypes that were not experimentally measured in the "GWAS1" ConLi + Gen sample. These PGSs, created and evaluated in large samples of predominantly European ancestry, stemmed from three recent publications and corresponded to the following traits: autoimmune disease [22], lymphocyte / monocyte / eosinophil / neutrophil / basophil percentage of white (blood) cells [23], and serum levels of 26 markers of in ammation [24]. After downloading and harmonizing weight les, we performed allelic scoring in ConLi + Gen using the sum method applied in Plink 1.9 [25].

Association analyses
The "GWAS1" (N = 853) and "GWAS2" (N = 1,258) samples were tested separately for associations of 701,031 SNPs in the ImmuneSet with: 1) Li response (responder / non-responder, de ned by Alda scores ≥ 7 or < 7, respectively), 2) total Alda score, 3) Alda subscale A, and 4) Alda subscale B (total). Because the most reliable continuous Li response phenotype has been previously shown to be the Alda A score, when excluding individuals with Alda B scores > 4 [10], we tested this as the primary continuous phenotype in our study. All association tests were performed applying an additive model in Plink 1.9, and all models were adjusted for age at recruitment, age-at-onset (AAO), sex, diagnosis and the rst principal components (PCs) obtained for each ImmuneSet genotypes dataset. PCA plots were explored to determine the optimal number of PCs to be used as covariates for each sample. Therefore, the rst ve PCs were used as covariates for "GWAS1" while the rst six PCs were used for "GWAS2". Population strati cation due to ancestry was successfully corrected by the selected numbers of PCs (Supplementary Figures: Figure S2). Next, the association results for Li response from "GWAS1" and "GWAS2" samples were meta-analyzed using the weighted-z (METAL) method applied in Plink 1.9. These meta-analysis results were QCed to exclude variants with I 2 heterogeneity index (I) > 40 and p-value for Cochran's Q statistic (Q) < 0.1 (highly heterogeneous). We searched rst for associations at the commonly accepted thresholds for GWASs (genome-wide signi cance, p < 5x10 − 8 , and suggestive signi cance, p < 1x10 − 5 ). However, considering this a candidate gene rather than a genome-wide approach, we chose to look further into ndings with p < 1x10 − 4 , a threshold that has been previously used to select association ndings for follow-up in GWASs [26] and, therefore, represents an acceptable exploratory threshold.
In "GWAS1", the ImmuneSet was further tested for associations with other BP clinical phenotypes (i.e. AAO, the number of episodes of depression, mania and hypomania, as well as the presence of psychosis, alcohol and/or substance abuse, and suicidal ideation). These models were adjusted for age at recruitment, AAO (except when AAO was tested as phenotype), sex, diagnosis and the rst ve PCs. Statistical signi cance was considered as above.
Associations between PGSs and the various BP clinical phenotypes were tested using linear or binomial regression models, as appropriate, adjusted for age at recruitment, AAO (except when AAO was tested as phenotype), sex, diagnosis and population using the robustbase R package. Signi cance was set to false discovery rate (FDR) < 0.05. However, we also looked into the nominally signi cant ndings (p < 0.05) for exploratory purposes.

Downstream analyses
All variants under the p < 1x10 − 4 threshold were annotated for known regulatory effects on gene expression (i.e. expression quantitative trait loci, eQTLs) in all human brain, blood, spleen and thyroid tissues, as well as in immune cells (e.g. monocytes and macrophages) using Qtlizer [27].
A protein-protein interaction (PPI) network to explore the functional relevance of the identi ed genes associated with Li response was created using the ReactomeFIViz app [28] for Cytoscape 3.7 [29]. This analysis used as input a list composed of the ImmuneSet genes showing associations at the p < 1x10 − 4 threshold with the dichotomous and continuous Li response phenotypes. The network also incorporated "linker" genes (i.e. genes not in the input gene list that create indirect connections between input genes) to increase biological interpretability. Moreover, pathway overrepresentation analysis was performed on the PPI network (including linker genes) using the pathway enrichment network function of the app. Because the linker genes were not drawn from the ImmuneSet collection, we used the standard background genes of the ReactomeFiViz app for this analysis. The resulting overrepresented pathways were ltered to exclude terms that: 1) had FDR > 0.05, 2) corresponded to a speci c disease (e.g. bladder cancer, herpes virus infection), 3) had less than two genes overlapping between the pathway set and the network set, and/or 4) the overlap with the pathway set represented less than 3% of genes in the set. Additionally, we repeated the pathway overrepresentation analysis including not only the variant mapped genes, but also the annotated eQTL genes.
For associations with clinical phenotypes in the "GWAS1" sample, functional analyses were performed using the GENE2FUNC tool of the Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA-GWAS) platform [30]. The input gene lists included mapped and eQTL genes annotated for variants below the p < 1x10 − 4 threshold for each studied phenotype. Because eQTL genes were not drawn from our ImmuneSet collection, we used all protein-coding genes as background for these analyses. Overrepresented gene sets were those that showed FDR < 0.05, following a hypergeometric test, and a minimum of two overlapping genes. Curated gene sets from pathway databases in the "canonical pathways" category were preferred when available. Otherwise, Gene Ontology biological processes (GO_BPs) or any other available category (including GWAS Catalog trait associations) were taken. For GTEx-based enriched tissues of expression, as our focus is on immune-brain relationships, we kept only those enrichments corresponding to brain expression, as these are the most relevant tissues for the analysis of Li response and clinical features of BP.
In addition, the relative importance for (dichotomous) Li response of the calculated PGSs in "GWAS1" was assessed through a machine learning (ML) screening approach using the Auto Model extension of RapidMiner Studio. This applied various classi cation algorithms to the raw PGS data. Auto Model provides the following models: Naïve Bayes, Generalized Linear Model, Logistic Regression, Fast Large Margin, Deep Learning, Decision Tree, Random Forest, Gradient Boosted Trees and Support Vector Machines. Because ML algorithms are sensitive to class imbalance, an equal number of responder and non-responder individuals were randomly selected for the analysis (N = 657) using the sample method of the Python's Pandas library. The resulting le with balanced classes was used as input for the ML screening in Auto Model, where the sample was randomly divided into training (60%) and test (40%) sets. Default parameters for all algorithms were applied. Given that different types of ML algorithms can differ in their feature selection procedure due to their inherent characteristics, here, features were considered important for Li response, with either a positive (i.e. favoring response) or a negative (i.e. favoring nonresponse) effect, when at least two algorithms selected the same feature with the same effect direction as important for the classi cation task.

RESULTS
Immune-related genes associated with response to Li treatment in BP After excluding individuals with missing phenotypic data (age and/or AAO), the effective sample sizes for the association analyses in ConLi + Gen were 853 and 1,258 in "GWAS1" and "GWAS2", respectively. A basic description of both samples is shown in Table 1. In general, there were more female than male patients in both samples and there were minimal differences in the mean ages at recruitment and disease onset between "GWAS1" and "GWAS2". Therefore, the total sample size of our meta-analyses of Li response was 2,111, including 606 (28.7%) responders and 1,505 non-responders for the dichotomized variable, which included 1,224 (58%) females and 887 males, with mean age 47 (± 14) years and mean AAO 25 (± 11) years. For the continuous Li response phenotype (i.e. Alda A score, excluding individuals with Alda B score > 4), the effective sample sizes were 828 for "GWAS1" and 1,044 for "GWAS2". After the post-meta-analysis exclusion of variants with heterogeneous effects between both ConLi + Gen samples, a mean of 556,196 SNPs remained in each set of summary statistics. This was higher for the continuous phenotype, in which 625,818 SNPs remained. Age-at-onset (mean ± SD) 25 ± 10 25 ± 11 25 ± 11 We found no associations with Li response at the genome-wide GWAS threshold (p < 5x10 − 8 ). At the suggestive threshold for GWAS (p < 1x10 − 5 ), the dichotomous Li response phenotype and Alda B (total) showed associations with FAT3 (best SNP: rs4313539, p = 2.1x10 − 6 , z = 4.744), and with ADAMTS5 (best SNP: rs162501, p = 9.18x10 − 7 , z = 4.909) and GRID2 (best SNP: rs62312225, p = 2.71x10 − 6 , z = 4.692), respectively (Supplementary File 1: Table 1). At the exploratory threshold (p < 1x10 − 4 ), when considering linkage disequilibrium (LD), we identi ed between 9 and 12 genomic loci in relation with different aspects of Li response (Supplementary File 1: Table 1). The most signi cant SNPs from the analyses of the dichotomous and continuous phenotypes mapped to FAT3 and BMPR1A, respectively (Table 2). In total, 42 genes were implicated in the response to Li in patients with BP from our exploratory analyses in the ConLi + Gen cohort (Supplementary File 2: Table 1). As expected, there was a number of gene-based overlaps between different aspects of Li response, particularly between the dichotomous variable and total Alda score, and between the continuous variable and the other Alda variables. Twenty-four ImmuneSet genes were implicated in the primary Li response phenotypes (i.e. dichotomous and continuous) in our exploratory analysis (Fig. 1A). These were used as input for a network analysis to facilitate biological interpretation of the ndings. This network analysis provided known and predicted functional interactions between a subset of 21 input genes from our association results and 16 linkers drawn from the total of protein-coding genes in the background reference of the ReactomeFiViz app ( Fig. 1B). Functional analysis of the network uncovered an overrepresentation of crucial developmental pathways and regulatory networks, suggesting the involvement of processes such as assembly and stability of the cell-cell signaling machinery (e.g. adherens junction, E-cadherin signaling, focal adhesion, integrin signaling, L1 cell adhesion molecule signaling), neuronal development and function (e.g. neurotrophic signaling, lysophosphatidic acid receptor mediated events, regulation of pluripotency, Wnt signaling), as well as activation of in ammatory (e.g. sphingolipid signaling, S1P pathways, toll-like receptor signaling) and adaptive immune pathways (e.g. T and B cell receptor signaling) in the biological response to Li. Interestingly, processes such as angiogenesis, long-term potentiation, thyroid hormone signaling pathways, melanogenesis and sensory processing were also overrepresented in our network analysis (Supplementary File 2: Table 2). These observations suggest that variation in immune-related genes and its effects that expand beyond in ammatory and immune responses from early developmental stages contribute to determine the extent of the organism's response to Li treatment in patients with BP. Moreover, when eQTL genes were incorporated into the analysis, there was a marked overrepresentation of in ammatory and autoimmune disease pathways (e.g. asthma, type 1 diabetes mellitus, autoimmune thyroid disease, in ammatory bowel disease) and of vitamin D metabolism. The Wnt signaling, cell adhesion and adaptive immune pathways remained signi cantly overrepresented (data not shown). All protein-coding eQTL genes annotated for Li response phenotypes are shown in Fig. 2A.
Immune-related genes associated with clinical phenotypes in BP The association analyses of the ImmuneSet with speci c clinical features that were available for the "GWAS1" sample showed no associations at the genome-wide GWAS threshold. However, at the suggestive GWAS threshold there were, collectively, 100 associations between the ImmuneSet and AAO (3), number of depressive (54) and manic (14) episodes, and the presence of psychosis (1), substance use disorder (15) and/or suicidal ideation (13). These implicated 17 genes that associated to speci c clinical features (i.e. no overlaps were observed at this threshold Table 3). When we moved forward to the exploratory analysis, we identi ed 786 SNP-phenotype associations in total (Supplementary File 1: Tables 2-9). These implicated 166 immune-related genes ( Fig. 1A; Supplementary File 2: Table 1) mostly involved in adaptive immunity and in ammation. Beyond their immune functions, however, these genes play important roles in the development of the nervous system, signal transduction, synaptic processes and cell adhesion (Supplementary File 2: Tables 3-9). In particular, large numbers of associations were found for AAO and mood episodes (Table 3). In addition, 156 eQTL genes were collectively annotated for these phenotypes. Figure 2A shows the protein-coding eQTL genes annotated for each clinical feature.
Even when these showed no overlaps among the clinical phenotypes, there were some overlaps with the ImmuneSet genes meaning that, in some instances, the eQTL gene corresponded to the gene mapped to the variant while, in others, the eQTL gene was different from the mapped gene. A summary of exploratory ndings for each clinical phenotype studied is shown below. Age-at-onset. Fifty-four associations in 21 ImmuneSet genes were found for AAO in our study (Table 3, Supplementary File 1: Table 2). These genes were enriched for negative regulation of cell death and synaptic transmission, as well as expression in the cerebellum (Supplementary File 2: Table 3). The top variant, rs1248079 (p = 3.9x10 − 6 , beta = 2.75), mapped to an intronic region in GRK5 (G Protein-Coupled Receptor Kinase 5). Other important genes included PLD3, AKT2 and IL1B. The addition of eQTL genes to the functional enrichment analysis resulted in an additional overrepresentation of processes related to cellular stress responses.
Depression. With an effective sample of 692 individuals, 107 associations in 31 ImmuneSet genes for the number of depressive episodes were found (Table 3, Supplementary File 1: Table 3). These genes were enriched for synaptic processes, as well as expression in frontal and anterior cingulate cortices (Supplementary File 2: Table 4). While the top variant, rs55975329 (p = 1.7x10 − 7 , beta = 4.54), localized to an intron in BLNK (B cell linker), other important genes included PHLPP1, ZCCHC11 (TUT4) and CPPED1.
The addition of eQTL genes to the functional enrichment analysis resulted in an additional overrepresentation of axonal and synaptic components.  Table 5).
Mania. With an effective sample of 665 individuals, 116 associations in 32 ImmuneSet genes for the number of manic episodes were found (Table 3, Supplementary File 1: Table 5). These genes were enriched for neuronal development and differentiation, as well as expression in spinal cord and frontal cortex (Supplementary File 2: Table 6). The top variant, rs59134172 (p = 2.4x10 − 6 , beta = 1.62), localized to an intron in CNTN6 (Contactin 6). Other important genes included KALRN, PCDH9, PTK2B and CNTNAP2. Interestingly, we also found enrichment for response to serotonin re-uptake inhibitors in major depressive disorder (FDR = 0.042) and serum thyroid-stimulating hormone levels (FDR = 0.0028), from the GWAS Catalog trait associations, in mania-associated ImmuneSet genes. The addition of eQTL genes to the functional enrichment analysis had no impact on the overrepresented gene set classes.
Psychosis. The effective sample for our analysis of the presence of psychosis in BP was 692 individuals.
The addition of eQTL genes to the functional enrichment analysis had no impact on the overrepresented processes. However, this resulted in a considerable increase in overepresented brain tissues of expression, including the hippocampus, amygdala, hipothalamus, anterior cingulate cortex, putamen and substantia nigra.
Alcohol and substance abuse. The effective sample sizes for alcohol and substance use disorders in ConLi + Gen were 835 and 832, respectively. Twenty-nine SNPs in nine genes were associated with alcohol use, with the rs7698751 SNP in SCD5 (Stearoyl-CoA Desaturase 5) being the top association (p = 1.8x10 − 5 , beta = 2.61). Although no gene set enrichments were found for associations with alcohol abuse, other implicated genes included the BP-associated DGKH, NRG3 and RIMS1 (Supplementary File 1: Table 7). Interestingly, when incorporating the eQTLs genes into this functional analysis, overrepresentation of genes associated with mood swings, loneliness and anxious behaviors in the GWAS Catalog was observed (Supplementary File 2: Table 8). For substance use disorder, 78 associations implicating 17 genes were found (Supplementary File 1: Table 8). Genes were overrepresented in neurogenesis-related processes, with expression in cerebellum as well as frontal and anterior cingulate cortices (Supplementary File 2: Table 8). The top variant, rs7814474 (p = 4.5x10 − 7 , beta = 6.7), was mapped to an intron in TPD52 (Tumor Protein D52). Other genes included the BP-associated NRG1, the schizophrenia-associated PTPRM, as well as NOD1 and PRKCQ. In addition, the GWAS Catalog traits chronotype (FDR = 0.011) and serum thyroid-stimulating hormone levels (FDR = 0.047) were also overrepresented in substance abuse-associated ImmuneSet genes. The addition of eQTL genes to the functional enrichment analysis resulted in an additional overrepresentation of the phosphatidylinositol signaling system and expression in the amygdala, hippocampus and basal ganglia.
Suicidal ideation. Information on the presence of suicidal thoughts was available for 660 ConLi + Gen individuals. Based on these, 30 variants in seven genes were associated with suicidal ideation in the "GWAS1" sample (Table 3, Supplementary File 1: Table 9). The top SNP was rs2327882 (p = 3.9x10 − 6 , beta = 0.55), located in an intron of the JARID2 (Jumonji and AT-Rich Interaction Domain Containing 2) gene. Gene set enrichment analysis found overrepresentation of functions of nuclear receptors (Supplementary File 2: Table 9). Indeed, these terms related to RARB and THRB. The addition of eQTL genes to the functional enrichment analysis had no impact on the overrepresented gene set classes.

Immune-related genes showed pleiotropy for BP phenotypes
Not surprisingly, a number of genes showed shared associations with the different phenotypes included in our exploratory ImmuneSet analyses in ConLi + Gen (Fig. 1A). Therefore, these genes can be prioritized for follow-up studies due to their pleiotropic effects in Li response, clinical features, or both. In this way, we prioritized 21 genes that showed associations with more than one BP phenotype (Table 4). Here, we excluded genes associated with Alda A, total Alda B and/or Total Alda when there was no overlap with the dichotomous and/or continuous Li response phenotypes. However, a complete list of corresponding gene-phenotype associations (197 in total) is shown in Supplementary File 2: Table 1. Five of the prioritized genes (CNTNAP5, DSP, NFIB, BMPR1A, and HAS3) were associated with multiple Li response phenotypes, 10 of them (XKR4, NRXN1, RTN4, NRG1/3, ALK, GRK5, LRP1B, NPSR1-AS1 and CPPED1) were associated with multiple clinical features, and another six genes (BANK1, ROBO2, CNTNAP2, PCDH9, CDH12 and FAT3) were associated with Li response as well as with clinical features.
Polygenic scores for immune-related traits associated with BP phenotypes In addition to testing associations of the ImmuneSet with BP phenotypes, we calculated a set of 32 previously published immune-related PGSs, namely for: 1) (general) autoimmune disease, 2) the proportions of white blood cell populations and 3) in ammatory marker levels in serum (Supplementary File 2: Table 10). The overlap of variants between the PGS weight les obtained from PGS Catalog and the SNPs available in our ConLi + Gen "GWAS1" sample was, in general, better for the serum levels of in ammatory markers (80.6% in average) than for the other PGSs. The lowest valid SNP overlap was found for the PGS of general autoimmune disease (38.8%), suggesting that this PGS may relatively poorly index autoimmune disease signatures. For the proportion of white blood cell populations, the valid SNP overlap was also not fully satisfactory (49.4% in average). Nevertheless, we were able to obtain some valuable insights by using these scores (Fig. 2B). We found associations (FDR < 0.05) of the PGSs for CXCL1 (C-X-C Motif Chemokine Ligand 1; z = 14.9, p < 2x10 − 6 ), ECP (Eosinophilic Cationic Protein; z =

DISCUSSION
There is apparent mounting evidence of immune dysregulation in BP and other major psychiatric diseases. Nevertheless, some observations have originated from underpowered studies, resulting in a lack of reproducibility [3]. Therefore, it becomes crucial to gain a better understanding of the relationships between the immune and central nervous systems, and to discern between causes and consequences of disease. In very general terms, it can be assumed that the existence of genetic associations with a given phenotype indicates causal contributions of the associated loci to the studied phenotype. With this in mind, we sought to investigate how genetic factors relevant to immune activity relate to disease phenotypes, such as response to Li treatment, AAO and psychiatric symptoms, in patients with BP. Using an exploratory and extensive candidate gene approach, our study prioritized various genes and in ammatory markers that appeared to represent pleiotropic factors contributing to multiple phenotypes. However, we should note that, here, we refer to pleiotropy as the (suggested) association with multiple traits in the ConLi + Gen cohort and, by no means, have we implied that these features are independent from each other. In fact, it should be expected that, given the important correlation between psychiatric disorders, the features that we have studied in ConLi + Gen are, to some extent, also correlated with each other. Also, because genes can play different roles in different tissues and cell types, we observed widespread enrichments of biological pathways participating in the development and function of the brain. This suggested that the genes implicated in BP phenotypes might affect in parallel both the immune and nervous systems. Nevertheless, because we found no associations at the genome-wide GWAS threshold of signi cance, and those observed at the suggestive GWAS threshold were limited, our ndings are also consistent with a relatively weak effect of immune genetic factors over BP phenotypes.
The results of our exploratory assessment of associations of genetic polymorphisms in immune-related genes with Li response in the ConLi + Gen cohort suggest that variations in in ammatory and adaptive immune processes might contribute to the e ciency of the response to Li treatment in patients with BP.
Importantly, our network and gene set enrichment analysis uncovered an involvement of numerous biological pathways that participate in cell adhesion, migration and intercellular communication, helping in the development and maintenance of the central and peripheral nervous systems, as well as of the immune and vascular systems. Interestingly, many of these appear to converge in the participation of GSK-3β (glycogen synthase kinase-3 beta), as assessed through comparative overlap analysis of the enriched KEGG and Reactome gene sets. GSK-3β is involved in multiple major developmental pathways, such as the Wnt, Notch and Hedgehog signaling pathways. Genetic manipulation in mouse models has shown an antidepressant-like behavior upon GSK-3β knockdown in hippocampus, as well as cognitive, behavioral and biochemical changes associated with psychiatric disorders, including Alzheimer's disease, BP and schizophrenia, upon GSK-3β overexpression [31]. Li possesses a well-known inhibitory effect over GSK-3β [32]. Therefore, our analysis suggests that GSK-3β might be highly relevant for the response to Li treatment in BP. Our ndings are also in agreement with other epidemiological and molecular investigations of Li effects. For example, we repeatedly observed overrepresentation of gene sets related to thyroid function, such as thyroid-stimulating hormone signaling and autoimmune thyroid disease. This is in line with the reports of a reversible association of Li treatment with hypothyroidism, particularly in women [33,34].
Taken together, the encouraging literature supporting our ndings sparked our interest in exploring how a genetic measure (PGS) for in ammatory marker levels in serum might associate with Li response in ConLi + Gen. Our results supported a relationship between Li treatment response and the genetic in uence on the levels of various in ammatory markers circulating in serum. These included CXCL1, ECP, FAS, TRANCE and CCL4, molecules that regulate the activation and recruitment of immune cells (T and B lymphocytes, monocytes, macrophages, neutrophils and eosinophils).
The results of our exploratory assessment of associations of genetic polymorphisms in the ImmuneSet genes with BP's AAO, numbers of mood episodes and psychiatric comorbidities in ConLi + Gen identi ed various candidate genes particularly contributing to mood episodes and substance abuse, including XKR4, NRXN1, GRK5 and NRG1/3. These genes, besides their immune-related functionalities, seem to play important roles in neuronal development and function, according to our gene set enrichment analyses. Indeed, this could be corroborated by the literature in many instances. For example, NRXN1, a cell surface protein involved in cell-cell interactions, exocytosis of secretory granules and regulation of signal transmission, has been associated with autism, schizophrenia and nicotine dependence [35]. GRK5 has a role in the regulation of motility in polymorphonuclear leukocytes and in ammation [36,37]. It also regulates the activity of various G-protein coupled receptors, including neurotransmitter receptors [38]. XKR4, a phospholipid scramblase strongly expressed in brain tissue and activated by caspases, has been suggested to participate in the remodeling of neural networks by triggering microglial responses to the exposure of phosphatidylserine on axons, dendrites and synapses [39].
It is worth noting that three of our candidate genes for depressive episodes have previously shown associations with either major depression at the gene-level (DCC) [40] or mapped to schizophrenia loci (CR1L and DGKI) [41] in large GWASs. Moreover, one of our candidate genes for manic episodes, ESR2, and one for substance use disorder, BTN3A2, were previously associated at the gene-level with major depression in the large GWAS, while one of our candidates for hypomanic episodes and alcohol dependence, the BP-associated RIMS1 [26], mapped previously to schizophrenia loci. This, together with the overlaps with reported BP-associated genes that we observed, particularly concerning known eQTL genes, provides another level of support to the validity of our ndings.
Additionally, the results of our exploratory assessment of PGS associations led to interesting observations. For example, that activation of macrophages and neutrophils, re ected by the associations with the PGSs for various cytokines and chemokines produced by or targeting these cells (e.g. interleukin-8, CX3CL1, CXCL6/16), might contribute to disease AAO in ConLi + Gen. Indeed, these observations are supported by studies that have found increases in neutrophil counts in psychiatric disorders, including BP [42,43], as well as association of genetic polymorphisms in interleukin-1β, a pro-in ammatory cytokine produced by activated immune cells, including neutrophils and macrophages, with age of onset of depression in geriatric patients [44]. In addition, if we consider that macrophage/neutrophil activation was also a suggested mechanism of Li response in our study, it would be easy to speculate that activation of these cells might be linked with some aspect of the disease onset.
In conclusion, we performed an exploratory study that indicates a relationship between immunity and clinically relevant BP phenotypes at the genetic level, and pinpoints various interesting candidates for follow-up studies. We acknowledge that our study was limited by a relatively small sample size, particularly for the episodes of hypomania, and by incomplete overlap between the variants in the PGSs and our ConLi + Gen dataset. The latter, which likely resulted from a limited overlap among the different SNP arrays initially used to genotype samples in different collection centers, caused an incomplete indexing of the immune phenotypes of interest, and a relatively low rate of survival of correction for multiple comparisons in the PGS-BP phenotype relationships. Nevertheless, despite inherent limitations, we believe that our study provides valuable insight and furthers the understanding of the immune implications in BP. These results complement the evidence coming from epidemiological data and previous ndings in ConLi + Gen, and support the hypothesis that, to some extent, immune regulation might represent a feasible strategy to improve the symptomatology and treatment response in patients with BP.  Exploratory ndings for the immune-related polygenic scores calculated in ConLi + Gen. Signi cance of