Novel ABCD1 and MTHFSD Variants in Taiwanese Bipolar Disorder: A Genetic Association Study

TL;DR

This study identifies new genetic variants linked to bipolar disorder in the Taiwanese population, offering insights into potential diagnostic and treatment approaches.

Contribution

The study discovers novel ABCD1 and MTHFSD gene variants associated with bipolar disorder in the Taiwanese Han population.

Findings

The variant rs11156606 in the ABCD1 gene is strongly associated with bipolar disorder and is linked to regulatory gene expression.

The MTHFSD gene variant rs3829533 is in strong linkage with missense variants, suggesting possible protein function alterations.

Abstract

Background and Objectives: In recent years, bipolar disorder (BD), a multifaceted mood disorder marked by severe episodic mood fluctuations, has been shown to have an impact on disability-adjusted life years (DALYs). The increasing prevalence of BD highlights the need for better diagnostic tools, particularly those involving genetic insights. Genetic association studies can play a crucial role in identifying variations linked to BD, shedding light on its genetic underpinnings and potential therapeutic targets. This study aimed to identify novel genetic variants associated with BD in the Taiwanese Han population and to elucidate their potential roles in disease pathogenesis. Materials and Methods: Genotyping was conducted using the Taiwan Precision Medicine Array (TPM Array) on 128 BD patients and 26,122 control subjects. Following quality control, 280,177 single nucleotide polymorphisms (SNPs) were analyzed via chi-square tests, and linkage disequilibrium (LD) analyses were employed to examine the associations among key SNPs. Results: Eleven SNPs reached significance (p < 10−5), with the variant rs11156606 in the ABCD1 gene—implicated in fatty acid metabolism—emerging as a prominent finding. LD analysis revealed that rs11156606 is strongly linked with rs73640819, located in the 3′ untranslated region, suggesting a regulatory role in gene expression. Additionally, rs3829533 in the MTHFSD gene was found to be in strong LD with the missense variants rs3751800 and rs3751801, indicating potential alterations in protein function. Conclusion: These findings enhance the genetic understanding of BD within a Taiwanese cohort by identifying novel risk-associated variants and support the potential for using these markers in early diagnosis and targeted therapeutic strategies.

Figures

Full text

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, in which phenomenology can be extremely similar. When a true bipolar patient is assessed in a depressive phase, there is a relative risk of more than 40% of being mistakenly diagnosed as MDD [3,4]. Moreover, 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 [5]. Underestimated prevalence and difficulties in diagnosis have resulted in an overarching need for more precise and early diagnostic methods.

The first GWAS was published in 2005. Since then, researchers have become optimistic about the prospect of the genome–disease association approach [6]. 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 [7]. Since the first genome–disease association study of BD in 2007, a handful of risk loci have been identified, notably ANK3 [8], NCAN [9], CACNA1C [10], and ODZ4 [11]. 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 [12], RELN [13], and ANK [8]. More importantly, BD exhibits a polygenic architecture, in which numerous genes can collectively influence its pathogenesis through diverse mechanisms—such as gene regulation, epigenetic modifications, and gene–environment interactions—where the cumulative burden of many small-effect variants contributes to individual risk [7].

Findings from genome-wide association studies (GWAS) in Western populations often have limited transferability to Taiwanese and other non-Western groups due to differences in genetic architecture, allelic frequency distributions, linkage disequilibrium patterns, and distinct environmental and lifestyle factors [14]. 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

Our study is part of the TSGH clinical genetics project under the Taiwan Precision Medicine Initiative (TPMI), spearheaded by Academia Sinica in collaboration with 16 medical centers nationwide. The TPMI has established a comprehensive biobank and research infrastructure, enabling the collection of DNA samples for genetic analyses and facilitating both retrospective and prospective access to participants’ disease phenotypes for longitudinal monitoring. A distinctive feature of this program is its reciprocal approach, providing participants with their genetic risk profiles and opportunities to join validation studies aimed at refining disease risk-prediction models and implementing risk-stratified healthcare protocols.

We recruited patients diagnosed with BD (ICD-10 code F31) based on diagnostic interviews conducted by clinical psychiatrists according to the DSM, Fifth Edition (DSM-5) criteria. Genotyping data were collected from the TPM array, the third-generation Taiwan Biobank SNP array designed by Academia Sinica and Thermo Fisher Scientific specifically for the Taiwan Precision Medicine Initiative (TPMI) project. This array facilitates large-scale genetic studies within the Taiwanese population by capturing genetic variation unique to the Taiwan Han Chinese. It comprises genome-wide SNPs, markers for complex diseases, pharmacogenomic indicators, and variants exclusive to Taiwanese genetics. This emphasis on Taiwanese diversity ensures relevance to East Asian populations while highlighting rare variations not commonly observed elsewhere [15].

Participants in this study were recruited from medical centers and genotyped by Academia Sinica. DNA extraction and SNP 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 were 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.

To analyze the genotyping data from the TPM array, we first excluded SNPs with a typing call rate below 80%. Next, we removed variants with a minor allele frequency under 0.05 or a Hardy–Weinberg equilibrium p-value below 1 × 10^−6^. We then performed a chi-squared test for association, setting the significance level to p < 0.05, using PLINK 1.9 software (https://zzz.bwh.harvard.edu/plink/ (v1.90, accessed on 7 October 2024)) [16]. To address potential inflation of test statistics, we implemented Genomic Control (GC), developed by Devlin and Roeder [17], as our primary p-value adjustment method.

2.3. Variant Annotations and Functional Analysis

Genes were identified for variant annotations utilizing the RefSeq Database (https://www.ncbi.nlm.nih.gov/refseq/ (accessed on 5 December 2024)), as described in wANNOVAR (https://wannovar.wglab.org/ (accessed on 5 December 2024)) [18]. To compare the allele frequency in different racial populations, the public domain databases 1000 Genomes [19], Genome Aggregation Database (gnomAD) [20], and Taiwan BioBank (https://taiwanview.twbiobank.org.tw/index (accessed on 20 January 2025)) [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 [21], were elucidated using Enrichr (https://maayanlab.cloud/Enrichr/ (accessed on 20 January 2025)) [22].

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 [23]. In our study, we uploaded the variants rs11156606 and rs3829533 obtained from the association study using the LDproxy Tool, which is part of LDlink (https://ldlink.nih.gov/?tab=home (accessed on 14 December 2024)) [24,25]. 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 R^2^ (R-squared) were calculated based on allele frequencies. The measured values of D′ and R^2^ 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 [26].

3. Results

3.1. Demographics of the Selected Patients

Table 1 lists the demographic composition of the participants enrolled in this 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 gene–disease association study pipeline implemented in our investigation to identify genetic susceptibilities associated with BD. This study differentiated between a cohort case 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 a robust association with BD.

We conducted a search of the GWAS Catalog (https://www.ebi.ac.uk/gwas/ (accessed on 20 January 2025)) 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 Chromosome X, was recorded in the GWAS Catalog (accession ID GCST000821) as being associated with both BD and SCZ [27]. 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. 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.

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, R^2^ = 0.987) and rs3751801 (D′ = 1, R^2^ = 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 S1 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.

4. Discussion

This study represents an advancement in our understanding of the genetic underpinnings of BD, particularly in the Taiwanese population. The identified 11 novel genetic variants associated with BD could contribute to earlier and more accurate diagnosis by helping clinicians identify individuals with a heightened genetic risk. If these variants are confirmed through larger replication studies and integrated into a polygenic risk score (PRS), clinicians might better distinguish BD from unipolar depression, especially during ambiguous early phases. This stratification could lead to more targeted monitoring of at-risk individuals, potentially reducing the risk of misdiagnosis and facilitating timely intervention.

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.

For the analysis, sometimes a Bonferroni-corrected threshold is applied to account for the large number of tests (e.g., p < ~1.8 × 10^−7^ for 280,177 SNPs in our study). However, due to our modest sample size and the exploratory nature of this work, we adopted p < 1 × 10^−5^ as a suggestive significance cutoff. Conventional methods like Bonferroni and FDR, though effective at controlling for multiple comparisons, do not address GWAS-specific confounders such as linkage disequilibrium (LD) and population stratification. In particular, the Bonferroni method assumes complete independence among tests, which is rarely the case in GWAS because of extensive LD, leading to overly conservative thresholds that can mask true associations. To mitigate these issues, we implemented Genomic Control (GC) as our primary p-value adjustment method [17]. GC directly accounts for inflation in test statistics due to population structure by calculating and applying a genomic inflation factor (λ). This approach corrects for ancestry-related allele frequency differences, a critical concern in GWAS of genetically diverse cohorts. A comparative analysis of Bonferroni, FDR, and GC is provided in Revised Supplementary Table S2, demonstrating that GC maintains robust control of type I error while better accommodating population stratification [28,29].

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 [30]. The variant rs11156606, located intronically in the 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, and 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 [31]. It also plays a crucial role in cellular processes such as membrane transport [32,33,34], circulatory system stability [35], lipid transport, and homeostasis [35]. Dysregulation of ABC transporters has been implicated in various diseases, including Mendelian disorders [32], cancer, and drug resistance [34]. 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 [36]. 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 [37]. 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 [38]. 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 [27]. 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 [39]. 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 [40]. 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 [41]. 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 [42]. 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 [43], echoing Ozdogan et al.’s findings from measuring homocysteine levels in patients [44]. Animal models such as those of Akahoshi et al. suggest that a high-methionine diet could induce BD-like behaviors [45], 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 [46]. A 2016 study on methamphetamine-associated psychosis (MAP) identified MTHFSD as a candidate biomarker for RNA degradation [47], 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 [48], 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 in essential biological processes and their impact on human health, warranting further exploration of its mechanisms and therapeutic potential.

A limitation of our study is the small sample size, which restricts statistical power and the robustness of our findings, particularly affecting the analysis of high-OR SNPs. Future studies with larger cohorts are needed to better validate our results and enhance the predictive accuracy of the genetic markers. The small sample size also affects the statistical power to thoroughly validate the rs11156606 variant and its potential X-linked inheritance pattern, where recruitment of additional patients and their family members is required to further investigate rs11156606 and elucidate its functional implications in disease pathogenesis. In a meta-analysis of individuals of East Asian (EAS) ancestry, two BD-associated loci—rs117130410 (p = 3.68 × 10^−8^, OR = 1.31) and rs174576 (p = 7.78 × 10^−9^, OR = 0.86)—were identified [49]. These findings differ from our results, which may be attributed to differences in demographics, analytical methods, or threshold criteria. Future studies could aim to validate the variants identified in both studies. Another limitation of our study is the 80% call-rate threshold, which helps maintain reliable detection but excludes lower call-rate SNPs that might otherwise increase genetic diversity. We plan to enroll additional participants in a new clinical trial to increase sample size and use genetic imputation strategies, thereby enhancing both the robustness of our findings and the overall SNP diversity.

5. Conclusions

This study successfully identified novel genetic variants associated with BD in the 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, potentially paving the way for improved diagnostic and targeted therapeutic approaches.