Genetic Markers for Bipolar Disorder in the Taiwanese Popu-lation: Insights from a Genetic Association Study

1. Introduction

Bipolar disorder (BD) is a clinically severe mood disorder with a lifetime prevalence of 4%[1]. From 1990 to 2017, the incidence of BD increased by 47.74% and the disability-adjusted life years (DALYs) associated with the condition increased by 54.4% during this period [2], taking away years of healthy functioning from individuals with the illness.

Currently, BDs are primarily diagnosed by careful assessment of behavior combined with subjective reports of abnormal experiences to group patients into disease categories standardized in The Diagnostic and Statistical Manual of Mental Disorders (DSM) and The International Statistical Classification of Diseases and Related Health Problems (ICD). However, BD types I and II are difficult to diagnose accurately in clinical practice, particularly in early stages. Besides, nearly a quarter of adults (22.5%) and adolescents with major depressive disorder (MDD), followed up for a mean length of 12-18 years developed BD [3]. Underestimated prevalence and difficulties in diagnosis have resulted in an overarching need for more precise and early diagnostic methods.

Functional studies of candidate genes and their locations within genomic regions linked to schizophrenia (SCZ) or bipolar disorder (BD) through linkage studies were prevalent until around 2006, when genome-disease association studies began to explore potential risk gene variants. This is particularly important in psychiatry, given the unconvincing and inconsistent evidence from candidate gene studies and the genetic architecture for most diseases that seem to be polygenic [4]. Since the first genome-disease association study of BD in 2007, a handful of risk loci have been identified, notably ANK3 [5], NCAN [6], CACNA1C [7] and ODZ4 [8].In recent years, an increasing number of studies have reported hundreds of genes and proteins related to BD, many of which have been suggested as potential biomarkers for this disease, including BDNF [9], RELN [10] and ANK [5]. More importantly, many genes have been reported to influence the pathogenesis of BD through multiple possible mechanisms [4].

In this study, we aimed to demonstrate the capabilities of genome-disease association analysis using the TPM array to detect genetic variations of bipolar disorder (BD) in the Taiwanese population. This analysis aids in predicting specific Taiwanese risk-relevant single nucleotide polymorphisms (SNPs), which could be applied in the diagnosis of BD, by investigating individuals who participated in the Tri-Service General Hospital (TSGH) genetic study project. Additionally, we focused on identifying the potential mechanisms of the most influential characterized variants, as well as related variants in ABCD1 and MTHFSD, which may further inform treatment strategies for BD.

2. Materials and Methods

2.1. Study Participants and Ethical Approval

The data was collected from the TPM array, which aims to create a database containing detailed clinical information and genetic profiles of one million individuals from the Taiwanese Han population. Participants in this study were recruited from medical centers and genotyped by Academia Sinica. DNA extraction and SNPs identification were performed as follows: First, each participant’s genomic DNA was extracted and purified from 3 mL of peripheral blood collected in EDTA vacutainers using the QIAsymphony SP system (QIAGEN, Hilden, Germany). Second, the purified DNA from each participant was loaded onto the TPM array chip, and genome-type signals were detected using the Axiom GeneTitan system (Thermo Fisher Scientific, Sunnyvale, CA, USA). Conversion and quality control of SNP calling and sample annotation was performed using the Axiom Analysis Suite (Thermo Fisher Scientific, Sunnyvale, CA, USA).

2.2. Disease Association Analysis

Figure 1 illustrates the steps involved in the disease association analysis. Initially, 128 individuals diagnosed with BD (ICD-10 code F31), determined through diagnostic interviews conducted by clinical psychiatrists according to the DSM, Fifth Edition (DSM-5) criteria, were recruited from the TSGH. The control group comprised 26,122 individuals without BD. Comprehensive participant information is presented in Table 1.

The case group (128 bipolar disorder patients) and the control group (26,122 participants). Genotyping was performed using the Taiwan Precision Medicine (TPM) Array Chip. Data from 493,852 SNPs were filtered, and 280,177 SNPs were passed through the chi-squared test to detect risk factors.

To analyze the genotyping data from the TPM array, SNPs with a low typing call rate (< 80%) were filtered out, and the remaining data was subjected to chi-squared test association analysis using PLINK 1.9 software (https://zzz.bwh.harvard.edu/plink/) [11]. Variants with low quality (minor allele frequency less than 0.05 and Hardy-Weinberg equilibrium less than 1 × 10⁻⁶) were removed, and highly significant p-values (less than 1 × l0-5) were selected for further analysis.

2.3. Variant Annotations and Functional Analysis

Genes were identified for variant annotations utilizing the RefSeq Database (https://www.ncbi.nlm.nih.gov/refseq/), as described in wANNOVAR (https://wannovar.wglab.org/) [12]. To compare the allele frequency in different racial populations, the public domain databases 1000 Genomes[13], Genome Aggregation Database (gnomAD) [14], and Taiwan BioBank (https://taiwanview.twbiobank.org.tw/index) [15] were used. Further characterization of the genes associated with the identified variants involved examining their biological functions and molecular pathways. The biological processes of the genes, as defined by Gene Ontology [16], were elucidated using Enrichr (https://maayanlab.cloud/Enrichr/) [17].

2.4. Linkage Disequilibrium (LD) Analysis

LD analysis helps researchers understand which genes are likely to be functionally relevant to a trait or disease, as well as the number and location of the contributing genes[18]. In our study, we uploaded the variants rs11156606 and rs3829533 obtained from association study using the LDproxy Tool, which is part of LDlink (https://ldlink.nih.gov/?tab=home) [19,20] and analysis parameters were as follows: 1. Genome Build Version: GRCH37 ( rs11156606) and GRCH38 ( rs3829533). 2. Population: CDX (Chinese Dai in Xishuangbanna, China), CHB (Han Chinese in Beijing, China), and CHS (Southern Han Chinese). 3. LD measurements: D’ (D prime) and R2 (R-squared) were calculated based on allele frequencies. The measured values of D' and R² range from 0 to 1, where a value of 1 indicates complete disequilibrium, and a value of 0 indicates complete equilibrium. The base-pair window is 500,000 bp. Regulatory potential prediction was based on FORGEdb. [21]

3. Results

3.1. Demographics of the Selected Patients

Table 1 lists the demographic composition of the participants enrolled in the study stratified by sex and age group. Within the BD cohort, female participants slightly outnumbered male participants, with 69 females and 59 males. The control group was more extensive, comprising 11,793 males and 14,329 females, indicating a demographic balance that mirrors the expected population distributions. The age distribution within the BD cohort revealed the highest prevalence in the 50-59 age demographic, with a progressive decrease observed in both the younger and older cohorts. By contrast, the control group demonstrated substantial inclusivity across all adult age ranges, ensuring adequate age-matched controls for robust comparative analyses.

3.2. Study Workflow

Figure 1 shows the genetic-disease association study pipeline implemented in our investigation to identify genetic susceptibilities associated with BD. This study differentiated between a case cohort of 128 individuals diagnosed with BD and a sizable control cohort consisting of 26,122 participants. Genotyping was performed using the Taiwan Precision Medicine (TPM) array chip, designed to accurately detect a comprehensive spectrum of SNPs. The initial genotypic yield from this high-throughput platform was an extensive array of 493,852 SNPs.

The subsequent stage involved a filtering process to ensure data integrity and relevance, in which a substantial number of SNPs were excluded from the initial pool owing to factors such as low minor allele frequency, deviations from Hardy-Weinberg equilibrium, or missing data. After filtering out the data that could potentially confound the analysis, a curated set of 280,177 high-quality SNPs was obtained.

The selected SNPs were analyzed using the chi-square test to detect significant associations between SNP frequencies and BD incidence. Variants exhibiting a p-value of less than 10-5 were deemed significant and flagged for in-depth investigation to ascertain their potential as risk factor indicators.

3.3. Genetic Variants Identified, Statistical Analysis and Significance

According to our association study results, there were 11 genetic variants that demonstrated significant (p<10-5) associations with the condition (refer to the data points above the red line in the Manhattan plot in Figure 2). The Q-Q plot (Figure 2) shows significant deviations above the diagonal line, suggesting a potential link with BD. The genomic inflation factor, lambda, derived from the median of the observed and expected chi-square values was close to 1. This indicates minimal inflation and suggests that our test statistics are robust. Therefore, we were able to confidently interpret the results, affirming the minimal bias in the findings of our study.

The variants identified with high significance (specifically rs6427761, rs6658422, rs7530898, rs74419649, rs79221549, rs6496485, rs3829533, rs755981, rs6081464, rs7289613, and rs11156606) were distributed across various chromosomes. These SNPs were associated with several coding and non-coding genes, such as ABCD1, C20orf78, DPP10, LOC284930, LPP, MTHFSD, NR5A2, and non-coding RNA genes, including LINC01221 and NTRK3-AS1 (Table 2). Each variant demonstrated a high odds ratio (>1.5), indicating robust association with BD.

We conducted a search of the GWAS Catalog (https://www.ebi.ac.uk/gwas/) to further understand the novelty and relevance of these findings. Our search revealed that all variants except rs11156606 have not been previously reported to be associated with BD or any physical disease. Notably, rs11156606, located at the 153,741,041 bp position on the X chromosome, was recorded in the GWAS Catalog (accession ID GCST000821) as being associated with both BD and SCZ [22]. In our study, this variant showed a particularly high odds ratio of 2.36, suggesting a more than two-fold increase in the risk of BD in carriers of the allele.

The allele frequencies of the identified variants in our study population varied significantly, ranging from 15.7% to 33.2% in cases and 7.3% to 21.6% in controls. These variations highlight the distinct prevalence of genetic markers. Moreover, the allele frequencies in our study were closely aligned with those recorded in the Taiwan Biobank, which profiles the genetics of healthy Taiwanese individuals. This similarity suggests that the genetic makeup in our control group was representative of the general healthy population in Taiwan (Table 3).

Furthermore, a comparison with other global populations revealed that the allele frequencies in our control group were similar to those found in Asian populations. Notably, the frequencies of rs6496485, rs3829533, rs755987, and rs6081464 were higher in our study than in African, American, and European populations (Table 3). This suggests that these four variants may serve as specific biomarkers for BD in the Taiwanese population, offering the potential for more targeted diagnostic approaches.

3.4. Functional Annotations of Risk Genes

As shown in Table 4, we identified 64 significant Gene Ontology (GO) biological process terms with p-values less than 0.05. These terms are predominantly related to biological functions such as fatty acid metabolism, protein localization, potassium ion transport, vascular transport, and viral replication, which may be linked to BD. Notably, fatty acid metabolism emerged as a major significant GO term, particularly associated with the variant rs11156606 located in the ABCD1 gene.

Analysis of similar GO terms revealed three key cluster networks: (1) Fatty Acid Metabolism/Oxidative Stress Regulation, (2) Lipid Homeostasis and Membrane Potential Regulation, and (3) Localization of Lipid Transport and Membrane Proteins, as shown in Figure 3. Additionally, four ABCD1-related GO terms were identified with direct relevance to neuronal functions: GO:0106027 (neuron projection organization), GO:1990535 (neuron projection maintenance), GO:0043217 (myelin maintenance), and GO:0042552 (myelination). These findings suggest that mutations in ABCD1 play a crucial role in the functioning of the nervous system, potentially contributing to the pathology of BD and other neurological conditions.

The BP functions of the genes were identified using the GO database, with terms selected based on p-values < 0.05. The GO terms are listed in Table 4. Similarity networks were constructed using Navier–GO.

3.5. LD Analysis

The significant variant rs11156606, identified in the association study and located in intron 7 of ABCD1, was linked through LDproxy analysis (Figure 4)to a rare missense mutation, rs782720024 (c.436T>A, p.Phe146Ile), in exon 1. This mutation exhibits low allele frequencies, approximately 0.002, specifically within the EAS and Taiwanese populations, indicating its rarity in the Chinese population. This suggests that individuals in populations with this variant may have a higher risk of BD. Additionally, another SNP rs73640819 (c.131G>A), found in the 3' UTR with high LD values, plays a crucial role in gene expression regulation, including transcription, translation, RNA stability, and localization. Its allele frequency was notably the lowest in the Taiwanese population at 0.07 (Figure S1), highlighting population-specific genetic differences and their potential impact on disease susceptibility and gene regulation mechanisms.

Furthermore, LDproxy analysis revealed that the variant rs3829533, located in exon 6 of MTHFSD, is closely associated with two other SNPs, rs3751800 (D’=1, R2=0.987) and rs3751801 (D’=1, R2=1), based on data from the Chinese population. Notably, the associated SNPs located in MTHFSD exon 8 resulted in missense coding changes (rs3751800: c.727C>T/p. Arg243Cys and rs3751801: c.730G>C/p.Ala244Pro) (Supplemental Table S2 and Figure 5). (Figure 5) This finding is significant because it suggests that while the Taiwanese population may possess the synonymous variant rs3829533 at a relatively high frequency (0.238 in cases and 0.138 in controls, as shown in Table 2), which does not alter coding, the associated missense variants rs3751800 and rs3751801 could potentially modify the molecular function of MTHFSD owing to changes in amino acids, thereby increasing the risk of BD. Figure S2 illustrates that the allele frequencies of the two variants in the East Asian and Taiwanese populations (0.16 and 0.15, respectively) were higher than those in the African, American, and European populations.

Additionally, FORGEdb scores, as depicted in Figure 4 and Figure 5, predicted the likelihood of genetic variants functioning as regulatory elements. These scores range from 0 to 10, and are based on various regulatory DNA datasets, including transcription factor (TF) binding and chromatin accessibility. The FORGEdb scores for rs11156606, rs73640819, rs3829533, rs3751800, and rs3751801 are 7, 10, 7, 8, and 8, respectively. These high scores indicated a significant regulatory effect on ABCD1 and MTHFSD, further suggesting that these variants could be specific risk factors for BD in the Taiwanese population.

We investigated the nucleotide polymorphisms in the three variants within the MTHFSD gene (Figure 6). The variants rs3829533 and rs3751801 were synonymous; the former retained the amino acid threonine at position 175 of the MTHFSD protein translated from exon 6. Conversely, rs3751800 represented a missense variant, whereas its LD-associated variant rs3751801 induced amino acid substitutions at positions 243 and 244 of the MTHFSD protein (arginine to cysteine and alanine to proline, respectively). These changes occurred within exon 8 of MTHFSD.

The synonymous variant rs3829533 preserves the amino acid threonine at position 175 of the MTHFSD protein (translated from exon 6). In contrast, the LD-associated variants rs3751800 and rs3751801 result in amino acid substitutions at positions 243 and 244, respectively: arginine to cysteine and alanine to proline. These alterations occur within exon 8 of MTHFSD.

4. Discussion

This study represents an advancement in our understanding of the genetic underpinnings of BD, particularly in the Taiwanese population. By identifying 11 novel genetic variants associated with BD, our research not only contributes to the growing body of genomic knowledge, but also highlights the importance of conducting population-specific genetic studies. The discovery of these variants, most of which have not been previously reported in relation to BD, opens new avenues for understanding the complex genetic architecture of this condition. In particular, the high odds ratios observed for these variants underscore their potential role as robust biomarkers, offering insights for the development of more precise diagnostic tools and personalized treatment strategies. Furthermore, the alignment of allele frequencies with those in the Taiwan Biobank reinforces the relevance of our findings to Taiwanese demographics, providing a foundation for future studies to explore genetic risk factors and their implications in a targeted manner. BD is a complex disease with unknown causes, involving various factors such as demographics, genetics, and environment, some of which have strong evidence supporting their link to BD. Therefore, this study not only enriches our genetic understanding of BD but also underscores the critical need for tailored genetic research in diverse populations to uncover the nuanced nature of complex psychiatric disorders.

There is growing attention on the interaction between specific genes and the environment as well as the suggested involvement of fatty acid metabolism and neuronal function in the disorder [23]. The variant rs11156606, located intronically in ABCD1 gene, presents intriguing possibilities regarding its mechanistic role in BD. The protein encoded by this gene is a member of the ATP-binding cassette (ABC) transporter superfamily. ABC genes are divided into seven distinct subfamilies (ABC1, MDR/TAP, MRP, ALD, OABP, GCN20, White). This protein is a member of the ALD subfamily, which is involved in the peroxisomal import of fatty acids and/or fatty acyl-CoAs in organelles [24]. It also plays a crucial role in cellular processes such as membrane transport [25,26,27], circulatory system stability [28], lipid transport and homeostasis [28]. Dysregulation of ABC transporters has been implicated in various diseases, including Mendelian disorders [25], cancer, and drug resistance [27]. Several studies have indicated that fatty acid metabolism plays a significant role in BD development and manifestation. In this study, Leclercq S. et al. tested PUFA as potent inducers of the ABCD genes. The expression levels of ABCD2 and ABCD3 were significantly higher in n-3-deficient rats than in rats fed ALA- or DHA-supplemented diets, indicating sensitivity towards dietary PUFA [29]. In a study comparing patients with BD and healthy controls, subjects with BD had distinctly lower levels of omega-3 eicosapentaenoic acid (EPA) and higher omega-6 arachidonic acid levels, coupled with increased plasma IL-6 and TNF-α levels [30]. We established a connection between the genotype of the fatty acid desaturase gene cluster and altered polyunsaturated fatty acid levels in BD. In one study, fatty acids in the erythrocyte membranes of patients with bipolar manic disorder and healthy controls were analyzed using thin-layer chromatography and gas chromatography. The results showed lower levels of arachidonic acid (20:4n-6) and docosahexaenoic acid (22:6n-3) in BD patients with BD compared to normal controls [31]. The association of the ABCD1 variant rs11156606 with altered fatty acid metabolism in BD underscores its potential role in the pathophysiology of the condition. This link, highlighted by the distinct PUFA profiles in patients with BD, suggests that dysregulation of ABC transporters, which are key in lipid transport and homeostasis, could be a critical factor in BD development.

The discovery of rs11156606 within intron region 7 of ABCD1 via disease-gene association interpretation exemplifies the effectiveness of such research in identifying genetic variants that are potentially crucial for disease predisposition. This study combined data from people with SCZ and BD to find shared genetic factors, and revealed that SNP rs11156606 is associated with BD and SCZ [22] Notably, the unique presence of the rare allele rs782720024 in East Asian and Taiwanese populations, along with the discovery of rs73640819, which shows high LD and is situated at the 3’ UTR of ABCD1, underlines the significance of population-specific genetic studies. These findings reveal the intricate nature of gene regulation and its influence on disease, highlighting how allele frequency variations across populations, especially their scarcity in Taiwanese individuals, may shed light on the unique regulatory mechanisms that affect gene expression and disease outcomes.

The discovery of the high LD-associated variant rs73640819 in the 3' untranslated region (3’ UTR) of ABCD1 highlights its importance beyond a mere sequence tailing the coding region. The 3’-UTR plays a pivotal role in the regulation of gene expression by influencing mRNA stability, translation efficiency, and subcellular localization, thereby fine-tuning protein synthesis. This region is integral to the complex post-transcriptional gene regulation machinery and acts as a critical regulator of the pathway from DNA transcription to mRNA translation into proteins [32]. Such insights open avenues for further exploration of how variations in the 3’-UTR affect gene function and disease mechanisms, indicating a rich field for ongoing research.

The variant with c.131G>A changes introduces polyadenylation in the 3 ‘-UTR, highlighting the nuanced role of 3’-UTRs in genetic regulation. Unlike coding regions, 3’-UTRs are marked by AT-rich sequences that play a critical role in directing polyadenylation to these regions, thus influencing mRNA stability and cellular localization. This highlights the importance of 3’-UTRs in controlling gene expression through management of the mRNA lifecycle. Additionally, the diversity of mRNA isoforms attributed to alternative cleavage and polyadenylation (APA) allows for the production of multiple 3’ UTR variants from a single gene, each with unique regulatory capabilities [33]. This complexity underscores the significance of APA in the broader context of post-transcriptional regulation, opening paths for further investigation of its role in gene expression diversity and implications for cellular function and disease.

The Methenyltetrahydrofolate Synthetase Domain Containing (MTHFSD) gene plays a pivotal role in diverse physiological processes, including vascular health and cancer susceptibility, through its cytoplasmic activity in RNA binding, which facilitates crucial cellular functions [34]. It encodes a protein that is essential for the conversion of 5,10-methylenetetrahydrofolate to 5-methyltetrahydrofolate, a key step in the re-methylation of homocysteine to methionine, which is vital for cell health. Notably, the potential of S-adenosyl methionine (SAM) to alter lipid raft arrangements and influence receptor and transporter functions has been linked to therapeutic prospects in depression, BD, and SCZ [35]. Recent studies, including those by Lyu, highlight nutrient deficiencies, particularly in the conversion process from homocysteine to methionine, as is common in BD and SCZ [36], echoing Ozdogan et al.'s findings from measuring homocysteine levels in patients [37]. Animal models such as those of Akahoshi et al. suggest that a high-methionine diet could induce BD-like behaviors [38], pointing to a potential connection between high methionine intake and emotional states. Further research on MTHFSD expression identified it as a top differentially expressed gene in SCZ at various illness stages [39]. A 2016 study on methamphetamine-associated psychosis (MAP) identified MTHFSD as a candidate biomarker for RNA degradation [40], highlighting its role in folate catabolism and methionine metabolism. This connection is further supported by a meta-analysis by Hsieh et al., which found that individuals with BD generally have lower serum folate levels than healthy controls [41], underscoring the enzyme's role in maintaining active folate levels for vital biosynthetic pathways.

LD analysis revealed two correlated variants (rs3751800: c.727C>T/p. Arg243Cys and rs3751801: c.730G>C/p.Ala244Pro) within the MTHFSD gene, indicating its potential as a disease marker for BD in the Taiwanese population owing to its effect on MTHFSD function. These findings emphasize MTHFSD's critical involvement of MTHFSD in essential biological processes and their impact on human health, warranting further exploration of its mechanisms and therapeutic potential.

5. Conclusions

This study successfully identified novel genetic variants associated with BD in a Taiwanese Han population. These findings illuminate the intricate genetic basis of BD and highlight the significance of research on diverse populations. The discovered variants implicated pathways related to fatty acid metabolism, lipid homeostasis, and neuronal function, suggesting potential targets for future research. This study advances our understanding of BD's complex genetic architecture of BD, potentially paving the way for improved diagnostic and targeted therapeutic approaches.