Gene expression insights: Chronic stress and bipolar disorder: A bioinformatics investigation

1.   Background

Bipolar disorder (BD) is a psychiatric disorder that affects around 2% of the world's population ^[1,2]. It is a severe and chronic disorder that includes bipolar I disorder and bipolar II disorder, with a series of clinical manifestations like manic episodes, severe emotional disorders, neuropsychological deficits, physiological and immune changes, dysfunction, etc. Lifestyle often contributes to the development of the disease, including diet, physical activity (PA), sleep, stress management, social relationships (SR) and others ^[3]. As epidemiological investigations show, BD commonly starts in young adults or adolescence, leading to disability and premature mortality. With family genetic characteristics, it brings a heavy burden to families and society ^[4,5]. In addition, BD has been reported to be associated with high a risk of suicide, self-harm and other autotomy behaviors. Research has shown that the rate of suicide among patients with bipolar disorder (BD) is approximately 10–30 times higher than the rate in the general population. More than one fifth of BD patients, especially those untreated, end their lives by suicide. Equally discouraging is that 20–60% of BD patients attempt suicide at least once in their lifetime. In a meta-analysis, suicide attempte rates averaged for teenagers with BD were 7.44% per year ^[6].

Mental disorders have been considered a major social and medical burden. It is hard to prevent and treat because we do not fully understand how it affects the brain. However, recent discoveries in psychiatric genomics have revealed that many common genetic variants can influence the risk of developing these conditions, but each one has a small effect. By connecting these genetic findings to the ways in which they change brain functions and behaviors, we might be able to figure out how these conditions really work. This may help us develop better ways to diagnose and treat them ^[7].

The causes and etiology of BD are unclear. Some studies have suggested that BD may be related to some genetic factors, while others have shown that BD patients have widespread changes in their brain structure, such as reduced gray matter volume, enlarged ventricles and a damaged corpus callosum. These changes make BD difficult to treat clinically ^[8]. Therefore, it is crucial to identify the key genes that influence BD and how they affect the brain and behavior of BD patients. This may help improve the diagnosis and treatment of this disorder, which can have diverse and variable symptoms across different individuals and stages of life ^[9].

Chronic stress (CS) is a typical psychophysiological response that the body generates in response to undesirable, challenging and difficult circumstances or stressors. It is a key cause of affective disorders by influencing the hypothalamic-pituitary-adrenal axis (HPA-Axis) ^[10], Autonomic Nervous System (ANS), immune system, adrenocorticotropic hormone (ACTH), brain-derived neurotropic factor (BDNF), metabolic biomarkers, antioxidants, glucose, hemoglobin, C-reactive protein (CRP), cytokines, pro-inflammatory cytokines, anti-inflammatory cytokines and tumor necrosis factor (TNF) ^[11].

Chronic stress can impair one's coping ability and have harmful effects on their behavior and physiology, leading to a greater amount of allostatic overload. The brain and developing nervous system have a remarkable ability to adapt their structure and function to stressful and other experiences, such as by replacing neurons, reshaping dendrites and changing synapses. Stress can disrupt the balance of neural circuits that support cognition, decision-making, anxiety and mood, leading to altered behaviors and emotional states ^[12]. This disruption affects the whole-body physiology through neuroendocrine, autonomic, immune and metabolic pathways. These changes may be helpful in the short term, but if the stressor is gone and the emotional state remains along with the neural changes, such maladaptation needs to be treated with a combination of drugs and behavioral therapies.

Numerous studies have shown that chronic stress may increase the risk of developing or relapsing from various diseases and impair the human brain. Being exposed to stress can weaken the higher cognitive functions of the prefrontal cortex (PFC) ^[12], hippocampus ^[13], amygdala and other brain regions ^[14]. It can also reduce its top-down regulation of behavior.

It was found that chronic mild stress was related to the pathogenesis of BD. Previous studies have reported that patients with BD experience more negative and stressful life events than the general population. In addition, factors affecting stress activation can be eliminated through appropriate stress management, thereby improving the treatment and condition of BD ^[15]. To study the brain mechanisms of BD in animal models, researchers often use stress exposure to trigger BD-like behavior. There are different types of stress models, such as chronic unpredictable mild stress (CUMS), chronic social defeat stress (CSDS) ^[16,17] and chronic restraint stress (CRS) ^[7]. These models are affected by different kinds of stress, and they can induce different aspects of BD symptoms, such as depression, mania or anxiety. However, the influence and role of CS on BD need to be further clarified.

The cause and etiology of the disorder are unknown and likely involve multiple factors. Therefore, studying the mechanisms and genetic influences of bipolar disorder is important. First, it can help us identify biomarkers and subtypes of the disorder, which can improve diagnosis and prognosis, as well as helping us understand the genetic overlap and differences between bipolar disorder and other psychiatric disorders, such as schizophrenia and major depression. Furthermore, it can even address the health disparities and diversity issues that affect people with bipolar disorder from different ancestral backgrounds. Moreover, it is an important way to develop novel treatments that target the underlying biological pathways and causal mechanisms of the disorder.

Recently, machine learning has transformed the field of biomedicine ^[18,19,20] by providing predictive models for the diagnosis ^[21,22,23], prognosis ^[24] and prediction of therapy response in diseases ^[25]. In drug response prediction, it enables the assessment of individual patient characteristics to forecast how they will respond to specific medications, advancing precision medicine ^[26,27]. It also has a crucial function in bioinformatics by examining extensive biological datasets ^[28,29,30], assisting in gene expression research and unraveling complex patterns within genomic data ^[31]. Furthermore, machine learning algorithms have accelerated the process of discovering biomarkers ^[32,33] by effectively analyzing multi-omics data, detecting patterns that are suggestive of diseases ^[34] and facilitating the identification of possible biomarkers for early detection and accurate disease monitoring ^[35,36].

In this paper, we employed machine learning and bioinformatics to examine how chronic stress affects BD in several ways. The gene expression data from the postmortem brains of BD patients and healthy controls in two datasets revealed 11 chronic stress-related genes that were differently elevated in BD. We then analyzed five biomarkers using machine learning. Their expression and diagnostic value were verified in other datasets and used for functional enrichment, regulatory network analysis and medication prediction. In our work, chronic stress affected monoamine oxidase synthesis and breakdown, neuroinflammation, ion permeability, pain perception and other features of BD. This work underlined the importance of exploring the genetic influences of chronic stress on BD and other mental disorders and suggested that chronic stress biomarkers may be diagnostic and therapeutic targets for BD.

2.   Materials and methods

Figure 1 presents the flow of our study. Details are described in the following subsections.

Figure 1.  Flow diagram of the study.

DownLoad: Full-Size Img PowerPoint

2.1.   Database of human genes

In order to determine the specific genes associated with chronic stress in individuals with BD, we initially identified the genes that are expressed differently in individuals with BD compared to healthy individuals using publicly available data. We then compared these genes with the genes known to be related to chronic stress. The data utilized for this analysis was obtained from datasets in the GEO database (https://www.ncbi.nlm.nih.gov/geo), specifically GSE12649 ^[37], GSE53987 ^[38], GSE5388 ^[39] and GSE78936 ^[40]. The initial two datasets used in this study were the training sets. GSE12649 consisted of 67 samples, with a breakdown of 33 samples in the BD group and 34 samples in the control group. GSE53987 consisted of 36 samples, with 17 samples in the BD group and 19 samples in the control group. We utilized GSE5388, which consisted of 61 samples (BD: Control = 30: 31), and GSE78936, which consisted of 54 samples (BD: Control = 30: 24), as validation sets to confirm the expression levels of biomarkers.

2.2.   Database of chronic stress-related genes

We obtained 895 genes associated with chronic stress (referred to as chronic stress-related genes or CSRGs) from the Genecards database (https://www.genecards.org). This was achieved by conducting a search using the phrase "chronic stress-related genes" and setting the relevance score threshold to be greater than 5.

2.3.   Analysis of differentially expressed genes

Differential expression analysis was performed using the $\texttt{limma}$ package in the R 4.2.0 software. Volcano/difference ranking plots and heatmap plots were generated using $\texttt{ggplot2}$ and $\texttt{pHeatmap}$, respectively. To visualize the distribution of DEGs on chromosomes, we used the R packages $\texttt{RCircos}$ and $\texttt{RIdeogram}$. The DEGs were screened using the criteria of $|\log_2(FC)| > 1$ and $p$-value $<$ 0.05.

The KEGG pathway and GO enrichment of DEGs were analyzed using the R packages $\texttt{clusterProfiler}$ and $\texttt{GOplot}$.

To obtain the 11 CSRGs and then the top 5 biomarkers, we used the Cytohubba plugin software in Cytoscape (Version 3.9.0) ^[41]. The enriched pathways of GSEA were screened using the criteria of FDR $<$ 0.25 and $p$-value $<$ 0.05. To explore whether there is interaction between these 11 genes and their functions, we conducted correlation analysis and plotted them using $\texttt{corrplot}$. The GSEA software (Version 4.1.0) was used to clarify the potential mechanism in BD.

2.4.   Machine learning models

In order to identify the most informative CSRGs as biomarkers for BD diagnosis, we initially employed the "Least Absolute Shrinkage and Selection Operator Regression" technique using the $\texttt{glmnet}$ package in R for LASSO analysis. Subsequently, we employed machine learning algorithms including Gradient Boosting Machine (GBM), Random Forest (RF), Support Vector Machine (SVM) and Logistic Regression (LR) using the $\texttt{randomForest}$ and $\texttt{library(e1071)}$ libraries. A $p$-value less than 0.05 was deemed to be statistically significant.

2.5.   Identification of diagnostic genes

To screen the diagnostic genes, the expression levels of hub genes between BD patients and healthy controls were visually displayed in the form of scatter plots and boxplots. ROC curve analysis was performed, and the AUCs were calculated using the $\texttt{pROC}$ package in R to determine the predicted values of the hub genes. Diagnostic genes were selected from the training set and validation set using the criterion of AUC $>$ 0.700 ^[42].

2.6.   Construction of miRNA-gene regulatory network

The miRNet database is a web-based platform that provides a comprehensive analysis of miRNA-target interactions and functional annotations. We used the miRNet database to predict the interaction between diagnostic genes. To screen the diagnostic genes, the expression levels of hub genes between BD patients and healthy controls were visually displayed in the form of scatter plots and boxplots. ROC curve analysis was performed, and the AUCs were calculated using the $\texttt{pROC}$ package in R to determine the predicted values of the hub genes. Diagnostic genes were selected from the training set and validation set using the criterion of AUC $>$ 0.500.

2.7.   Potential therapeutic drug prediction

The gene-drug interaction network is a graphical tool that helps predict potential new targets for drugs. In the study, the Network Analyst (https://www.networkanalyst.ca) and DGIdb databases (https://dgidb.genome.wustl.edu) were used to analyze and identify biomarkers–drug interactions. The results of these two databases were used to predict drugs that may interact with biomarkers.

2.8.   Statistical analysis

GraphPad Prism Version 8.2.0 for Windows (https://www.graphpad.com) was used to perform statistical tests. The number of independent experiments is indicated in the plots. Bar graphs represent the mean ± SD, unless stated otherwise. Boxplots represent the median (box: Two quantiles around the median; whiskers: Minimum and maximum value; points superimposed on the graph: Individual values). The Shapiro-Wilk test was used to test normality. One or two sample student's $t$-tests (always two-sided) or ANOVA followed by a post hoc Tukey test were used for normally distributed data, and a nonparametric Mann-Whitney U-test was used for non-normal data. Correlations were calculated using the Spearman correlation coefficient with a two-tailed analysis. Pairwise comparison using a nonparametric Wilcoxon rank sum test was used to determine significant changes in BD patients. A $p$-value $<$ 0.05 was considered statistically significant.

3.   Results

3.1.   Identification and analysis of differentially expressed CSRGs in BD

3.1.1.   Identification of differentially expressed CSRGs

Using the R package $\texttt{limma}$, we performed differential analysis between the BD group and the Control group in the two datasets, GSE12649 and GSE53987, respectively. The results showed that in the GSE12649 dataset, there were 1311 differentially expressed genes in the BD-vs-Control comparison group, of which 477 were upregulated and 834 were downregulated (Figure 2A). In the GSE53987 dataset, there were 2597 differentially expressed genes in the BD-vs-Control comparison group, of which 1431 were upregulated and 1166 were downregulated (Figure 2B). The criteria for selecting differentially expressed genes were $|\log_2($fold change$)| > 0$ and $p$-value $<$ 0.05. In order to screen out differentially expressed genes with consistent expression trends from the two datasets GSE12649 and GSE53987, we performed Venn analysis on the upregulated and downregulated genes of these two datasets. The results showed that there were 61 common upregulated genes and 39 common downregulated genes (Figure 2C, D). Then we enriched the GO analysis (Figure 2E) and KEGG pathways (Figure 2F) of 100 differentially expressed genes.

Figure 2.  (A) The Volcano plot shows the 1311 differentially expressed genes in the GSE13649 dataset, including 477 upregulated genes and 834 downregulated genes; (B) The Volcano plot shows the 2597 differentially expressed genes in the GSE53987 dataset, with 1431 upregulated genes and 1611 downregulated genes; (C) Venn analysis of the upregulated genes in the two datasets; (D) Venn analysis of the downregulated genes in the two datasets; (E) GO enrichment analysis of 100 differentially expressed genes in the datasets; and (F) KEGG pathway analysis of 100 differentially expressed genes in the datasets.

DownLoad: Full-Size Img PowerPoint

In order to screen out differentially expressed genes related to CSRGs from the above common differentially expressed genes, we performed Venn analysis on the three groups and found that there were 11 CSRGs (Figure 3A) that were differentially expressed between the BD and Control samples, which are MAOA, TCF7L2, SOX9, HIF1A, NTRK2, LGALS3, TRPM3, IGFBP6, AIF1, GPX4 and ALOX5AP.

Figure 3.  (A) Venn analysis of the upregulated and downregulated genes in the two datasets; (B) Heatmap of the 11 differentially expressed CSRGs in the GSE12649 dataset; and (C) Heatmap of the 11 differentially expressed CSRGs in the GSE53987 dataset. In (B) and (C), each small square in the heatmaps represents a differentially expressed CSRG, and its color indicates the gene expression level. The larger the expression level, the darker the color (red for high expression, blue for low expression).

DownLoad: Full-Size Img PowerPoint

The heatmaps of the 11 differentially expressed CSRGs in the two datasets GSE12649 and GSE53987 were displayed, as shown in Figure 3B, C. The first row indicates the sample grouping: green for Control samples and orange for BD samples. Each row represents the expression level of each gene in different samples, and each column represents the expression level of all differentially expressed CSRGs in each sample. The tree analysis on the left shows the results of cluster analyses of different genes from different samples.

3.1.2.   Differential distribution and correlation analysis of CSRGs

To further identify the differential distribution of 11 differential chronic stress-related genes in BD, we used the R packages $\texttt{RCircos}$ and $\texttt{RIdeogram}$ to visualize the distribution on their chromosomes (Figure 4A, B). In order to explore whether there is interaction between these 11 genes and their functions, we used the $\texttt{corrplot}$ package to conduct correlation analysis and plot them (Figure 4C). We can see that there are significant correlations between most of the genes.

Figure 4.  (A) Differential distribution of 11 CSRGs on chromosomes; (B) Distribution of 11 different CSRGs on chromosomes. Note: Blue represents a lower density of genes on the chromosome, and red represents a higher density of genes on the chromosome; (C) Correlation analysis results. Note: The numbers in the box represent the correlation coefficient; $\times$ means that there is no significant correlation between the two genes on the horizontal and vertical axes; and (D) Functional enrichment analysis of differential CSRGs.

DownLoad: Full-Size Img PowerPoint

3.1.3.   Functional enrichment analysis of differential CSRGs

In addition, we utilized the metascape database (https://metascape.org/gp/index.html#/main/step1) to condense the functional analysis of 11 CSRGs. The enrichment analysis revealed a significant association between these genes and specific factors, including the inhibition of the canonical Wnt signaling pathway and the response to oxidative stress (Figure 4D).

3.2.   Enrichment and evaluation of biomarkers

3.2.1.   Screening biomarkers by machine learning

Biomarkers, also referred to as biological markers, are quantifiable indicators of the biological state and condition of diseases. Additionally, they have the capability to assess and analyze biological processes, pathogenic processes or pharmacologic responses to a therapeutic intervention. From this perspective, we believe that biomarkers can facilitate a deeper investigation into the underlying mechanism of chronic stress in individuals with BD. Biomarkers are quantifiable indicators of a biological state or condition that can be utilized to investigate regular biological processes, disease-causing processes or the effects of a therapeutic intervention. From this perspective, we employed biomarkers to assess the levels of gene expression across different datasets.

The LASSO (least absolute shrinkage and selection operator) method is a compression estimation technique that reduces the variable set (order reduction). It compresses the coefficient of the variable by constructing a penalty function, which makes some regression coefficients become 0, thereby achieving the purpose of variable selection.

We used the $\texttt{glmnet}$ package to conduct LASSO regression analysis on 11 differential CSRGs in the GSE12649 dataset and obtained two common graphs in LASSO regression: One is the graph of gene coefficient (Figure 5A), and the other is the error graph of cross-validation (Figure 5B). After that, we screened out five characteristic genes, including IGFBP6, ALOX5AP, MAOA, AIF1 and TRPM3.

Figure 5.  (A) Graph of gene coefficient in LASSO regression analysis on 11 differential CSRGs, conducted by using the $\texttt{glmnet}$ package in R; (B) Graph of error graph of cross-validation in LASSO regression analysis on 11 differential CSRGs, conducted using the $\texttt{glmnet}$ package; (C) The ROC curve of machine learning models determines the value of these signature genes in diagnosing BD. The LR model had the best AUC value of 0.81; (D) The distribution difference of risk scores in the two samples of the BD group and the Control group in GSE12649. ($N_{BD} = 33$, $N_{Con} = 34$; two-tailed student's $t-$test, ^****: $p < 0.0001$); and (F) The distribution difference of risk scores in the two samples of the BD group and the Control group in GSE53987. ($N_{BD} = 33$, $N_{Con} = 34$; two-tailed student's $t-$test, ^****: $p < 0.0001$).

DownLoad: Full-Size Img PowerPoint

We employed machine learning algorithms, specifically GBM, RF, SVM and LR models, to assess the diagnostic significance of these signature genes in BD (Figure 5C). The findings indicate that all the area under the Receiver Operating Characteristic (ROC) curve values in the models are greater than 0.5, demonstrating the reliability of our biomarkers in diagnosing BD with a favorable diagnostic value. The LR model outperformed other models with an AUC value of 0.81. Thus, we utilized the coefficient of characteristic genes in the logistic regression model to compute the risk score. The risk score is calculated using the following formula: (-3.8673 $\times$ IGFBP6) + (-1.7361 $\times$ ALOX5AP) + (0.3371 $\times$ MAOA) + (-1.225 $\times$ AIF1) + (1.7401 $\times$ TRPM3). Subsequently, we utilized the risk scores to compute the disparity in distribution between the BD and Control groups in GSE 12649. The disparity between the groups is notably evident (Figure 5D, E), demonstrating that the amalgamation of our biomarkers could serve as an effective means of assessing and diagnosing BD.

3.2.2.   Construction and evaluation of diagnostic nomogram in BD

The diagnostic nomogram calculates the value levels of influence factors according to the contribution degree of outcome (the size or size of the regression coefficient) and then adds all scores to get the total score. It is a graphical tool that helps practitioners calculate the probability of having a particular condition based on the results of diagnostic tests. The nomogram calculates the value levels of influence factors according to the contribution degree of outcome (the size or size of the regression coefficient) and then adds all scores to get the total score. Therefore, the predicted value of the outcome event for the individual is calculated through the function conversion relationship between the total score and the probability of the outcome event.

The regression coefficient is a measure of how much a dependent variable changes when an independent variable changes. In linear regression, coefficients are the values that multiply the predictor values. The aim of linear regression is to find the regression coefficients that produce the best-fitted line.

Based on the samples from the GSE12649 dataset, we use the R language package $\texttt{rms}$ to construct a diagnostic nomogram of the biomarkers. We mapped all the biomarkers, and each factor corresponds to a score. We added up the scores of each factor to the total score and then predicted the survival rate of patients according to it. The higher the score, the lower the survival rate (Figure 6A). The calibration curve is drawn based on the above prediction model to verify the diagnostic nomogram (Figure 6B). The slope of the calibration curve distribution is used because the closer it is to 1 on both sides of the diagonal, the more accurate the prediction will be. According to the analysis, the Area Under the Curve (AUC) is 0.565 and the $c-$index/$c-$statistic of the ROC is 0.783 (Figure 6C), proving that the predictive power of this model is sufficient.

Figure 6.  (A) Diagnostic nomogram for the survival rate of patients; (B) The Calibration Curve of the diagnostic nomogram; (C) ROC Curve of biomarkers, Nomogram AUC is 0.783; (D) Standardized Net Curve Analysis: the established model and the nomogram we performed before. The curve is plotted with the horizontal black line indicating that all samples were negative (Pi $<$ Pt), and the benefit of none being treated was 0. The oblique gray curve means that all samples are positive and all people are treated; and (E) Impact Curve of prediction. Under High Risk Threshold values, the Number high risk curve represents the Number of people classified by the model as positive (high risk), and the Number high risk with event curve represents the number of people who are truly positive.

DownLoad: Full-Size Img PowerPoint

We compared the prediction efficiency of the pure model and the model with biomarkers by performing decision curve analysis on the established model and the nomogram.

The decision curve analysis showed that the benefit rate of the nomogram model with biomarkers included was higher than that of the single biomarker feature (Figure 6D). The decision curve and influence curve of biomarkers, respectively, indicate whether the predicted results of the model (red curve) and the real situation (blue curve) are consistent as the probability threshold increases. We further plotted the biomarker Impact Curve based on the biomarker decision curve (Figure 6E), which is a graphical representation of the nomogram model's predicted risk stratification by assuming 1000 people. Under High Risk Threshold values, the Number high risk curve represents the Number of people classified by the model as positive (high risk), and the Number high risk with event curve represents the number of people who are truly positive. From 0 to 1, the curve of "Number high risk" and the curve of "Number high risk with event" under the high risk threshold mostly overlap, indicating that the nomogram model has more accurate forecasting ability.

3.3.   Verification of biomarkers in BD

We used the Wilcox test assay to analyze the expression levels of biomarkers in GSE12649 and GSE53987 (Figure 7A, B), showing all biomarkers were differentially expressed and MAOA had a significant difference in both sets. A Wilcox test is a statistical test used to compare two groups of data that are not normally distributed. It is a non-parametric test that does not assume any particular distribution for the data. The test is used to determine whether two groups of data have different medians. As the result showed, in GSE12649 and GSE53987, all biomarkers were found to be differentially expressed, and MAOA showed significant differences in both sets.

Figure 7.  (A) The expression levels of biomarkers in GSE12649. All five biomarkers were differentially expressed; AIF1 and MAOA have significant differences; (B) The expression levels of biomarkers in GSE53987. All five biomarkers were differentially expressed; IGFBP6 and MAOA have significant differences; (C) The expression levels of biomarkers in GSE5388. The trend of difference in all five biomarkers was the same as the result in GSE12649 and GSE53987; among them, ALOX5AP and MAOA have significant differences; and (D) The expression levels of biomarkers in GSE78936. The trend of difference in all five biomarkers was the same as the result in GSE12649 and GSE53987, but only IGFBP6 had a significant difference.

DownLoad: Full-Size Img PowerPoint

Then, we verified the expression of biomarkers in the validation datasets GSE5388 and GSE78936 (Figure 7C, D). The expression trend of all biomarkers in the GSE5388 dataset was completely consistent with the training set, and the expression of ALOX5AP and MAOA was significantly different in the BD group and the Control group. The expression trend of some biomarkers in the GSE5389 dataset was consistent with that in the training set. The expression trend of all biomarkers in the GSE78936 dataset was completely consistent with that in the training set, and IGFBP6 expression was significantly different between the BD group and the Control group.

3.4.   Functions analysis of biomarkers

3.4.1.   Biological functions of biomarkers

We used Uniprot (https://www.uniprot.org) and Genescard (https://www.genecards.org) to search for the biological functions of these five biomarkers (IGFBP6, ALOX5AP, MAOA, AIF1 and TRPM3). The full name and partial function are shown in Table 1.

Table 1.  Name, partial function and related articles of the biomarkers.

Gene	Full name	Function	PubMed
IGFBP6	Insulin-like growth factor-binding protein 6	IGF-binding proteins prolong the half-life of the IGFs and have been shown to either inhibit or stimulate the growth promoting effects of the IGFs on cell culture	24003225
ALOX5AP	Arachidonate 5-lipoxygenase-activating protein	Required for leukotriene biosynthesis by ALOX5 (5-lipoxygenase)	2300173; 8440384
MAOA	Amine oxidase [flavin-containing] A	Encoding mitochondrial enzymes which catalyze the oxidative deamination of amines, such as dopamine, norepinephrine, serotonin and some neurotransmitters. It plays important roles in the metabolism of neuroactive and vasoactive amines in both the central nervous system and peripheral tissues	20493079; 8316221; 18391214; 24169519
AIF1	Allograft inflammatory factor 1	Actin-binding protein that enhances membrane ruffling and RAC activation	15117732; 16049345; 18699778
TRPM3	Transient receptor potential cation channel subfamily M member 3	Calcium channel mediating constitutive calcium ion entry	21278253

---

 | Show Table

DownLoad: CSV

3.4.2.   GSEA functional enrichment analysis of biomarkers

We conducted Gene Set Enrichment Analysis (GSEA) on the five biomarkers (IGFBP6, ALOX5AP, MAOA, AIF1 and TRPM3) to study the main GO (Gene Ontology) enrichment and pathway analysis for possible biological processes and related molecular functions. First, we sequenced biomarkers according to the correlation with all gene expression; next, we set SIZE $>$ 20 and NOM.p.val $<$ 0.05 as the significantly enriched pathways; then we drew the top 5 pathways. The results showed that IGFBP6 was significantly enriched in mitochondrion organization, mitochondrial inner membrane, and mitochondrial protein-containing complexes (Figure 8A). ALOX5AP was significantly enriched in cell killing, leukocyte-mediated cytotoxicity, and regulation of leukocyte-mediated cytotoxicity (Figure 8B). MAOA was significantly enriched in DNA-binding transcription activator activity, proteasome complexes, RNA splicing and articles such as transesterification reactions (Figure 8C). AIF1S was significantly enriched in cell killing, leukocyte-mediated cytotoxicity, microglial cell activation, positive regulation of cytokine production and regulation of leukocyte-mediated cytotoxicity (Figure 8D). TRPM3 was significantly enriched for mRNA splicing via spliceosome, proteasome complex and protein Polybiquitination (Figure 8E).

Figure 8.  Gene set enrichment analysis of (A) IGFBP6, (B) ALOX5AP, (C) MAOA, (D) AIF1 and (E) TRPM3.

DownLoad: Full-Size Img PowerPoint

3.4.3.   Transcriptional and post-transcriptional level regulation of biomarkers

Using the hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget#!), we searched five biomarkers in turn, set Tissue as "Brain" and obtained the transcription factor (TF) of each of them (Figure 9A). Next, we used the miRDB database (https://mirdb.org) to predict the miRNAs that are likely to be regulated at the post-transcriptional level of the biomarkers and visualize the mRNA-miRNA network in Cytoscape (Figure 9B). Transcriptional and post-transcriptional level regulation of biomarkers can affect the expression of genes and their corresponding proteins. Transcriptional regulation refers to the process by which gene expression is controlled at the level of transcription, which involves the binding of transcription factors to specific DNA sequences in the promoter region of a gene. Post-transcriptional regulation, on the other hand, refers to the regulation of gene expression after transcription has occurred, which includes processes such as RNA splicing, RNA editing and RNA stability. These processes can affect the amount and type of protein that is produced from a given gene, which can have downstream effects on cellular processes and disease states. For example, we found that transcription factors like REST, also known as Repressor element-I silencing transcription factor or neuron-restrictive silencer factor, may work a lot in BD through the biomarkers MAOA, TRPM3 and IGFBP6.

Figure 9.  (A) mRNA-TF network of biomarkers (34 nodes, 60 edges). Note: The blue rectangle represents the biomarker, and the pink diamond represents the TF; (B) mRNA-miRNA network of biomarkers (240 nodes, 246 edges). Note: The blue rectangle represents the biomarker, and the orange hexagon represents the miRNA. The blue rectangle represents the biomarker, and the pink diamond represents the TF; (C) Analysis and identification of biomarkers---drugs in Network Analyst (25 nodes, 24 edges); and (D) Analysis and identification of biomarkers---drugs in DGIdb (30 nodes, 28 edges).

DownLoad: Full-Size Img PowerPoint

3.4.4.   Analysis and identification of biomarkers–drugs

We used the Network Analyst (www.networkanalyst.ca) and DGIdb (https://dgidb.genome.wustl.edu) databases to analyze and identify biomarkers and drug interactions, as biomarkers may interact with drugs. The results of these two databases for predicting drugs (Figure 9C, D). The analysis and identification of biomarkers-drugs can help identify potential drug targets for a particular disease. The information we concluded can be used to develop new drugs or repurpose existing ones for the treatment of BD.

4.   Discussion and conclusions

In this study, we aimed to analyze the relationship between chronic stress and bipolar disorder from a genetic perspective, especially the expression of related genes in the prefrontal cortex.

Therefore, we explored the possible mechanisms of chronic stress affecting bipolar disorder (BD) through bioinformatics analysis of gene expression data from the postmortem brains of BD patients and healthy controls. We chose datasets GSE12649 and GSE53987 for enrichment analysis of the bioinformatics of BD, which contain gene expression after large-scale DNA microarray analysis in the prefrontal cortex of BD patients' postmortem brains. The differences between the BD group and the Control group in the GSE12649 and GSE53987 datasets, respectively, The results showed that a total of 1311 differentially expressed genes were obtained in the BD-VS-control comparison group in the GSE12649 dataset. There were 477 up-regulated genes and 834 down-regulated genes. A total of 2597 differentially expressed genes were obtained in the BD-vs-Control comparison group in the GSE53987 dataset, including 1431 up-regulated genes and 1166 down-regulated genes. Venn analysis was performed on up-down-regulated genes in GSE12649 and GSE53987. The results showed that 61 up-regulated genes and 39 down-regulated genes were obtained, respectively. In order to screen out CSRGS-related differential genes from the above shared differential genes, analysis was performed on the three genes, and the results showed that 11 CSRGs were differentially expressed in BD and Control samples. We then visualized the distribution of 11 differential CSRGs on chromosomes and performed PPI analysis and functional enrichment analysis of these genes. Subsequently, we performed LASSO regression analysis on 11 differential CSRGs in the GSE12649 dataset and screened out 5 characteristic genes, namely IGFBP6, ALOX5AP, MAOA, AIF1 and TRPM3. Subsequently, GBM, RF, SVM and LR machine learning models were used to determine the value of these feature genes in the diagnosis of BD. The results showed that the LR model had the best AUC value, and then the coefficient of feature genes in the LR model was used to calculate the risk score, and the distribution difference of risk scores between the BD group and the Control group was shown. Moving forward, there is potential for further exploration of alternative machine learning methods ^[43,44], alongside the utilization of data-centric approaches ^[45].

In order to enhance credibility, relevant biological markers were verified in the dataset containing human gene expression, GSE5388 and GSE78936. Subsequently, to further analyze the effects of chronic stress on BD through diversity gene expression, we enriched and analyzed their functions. We queried the biological functions of these five biomarkers, and based on the biomarkers, we constructed diagnostic nomographs, decision curves, GSEA functional enrichment analysis, regulatory networks, drug prediction and other analyses. Finally, we showed the expression levels of biomarkers in the datasets: GSE12649, GSE53987, GSE5388 and GSE78936.

We screened a total of five biomarkers: IGFBP6, ALOX5AP, MAOA, AIF1 and TRPM3. We predicted the therapeutic drugs of the biomarkers and constructed the regulatory network based on the biomarkers. In subsequent studies, the specific molecular mechanisms of these five biomarkers in BD can be further explored.

We can further discuss the relationship between chronic stress and BD by analyzing biomarkers' related functions. MAOA is one of the genes most commonly associated with a variety of psychiatric disorders, including major depressive disorder, bipolar disorder and antisocial behavior ^[46,47]. Monoamine oxidase A (MAOA) plays an important role in balancing the secretion of serotonin in the brain, which affects areas of the brain that are used to control anger. Therefore, the level of MAOA gene expression directly affects the emotional level of individuals and has an impact on their behavior. The MAOA gene regulates the activity of an enzyme in the brain that breaks down dopamine and other neurotransmitters, similar to the "feel-good" chemicals found in antidepressants. Low expression of the MAOA gene increases levels of monoamines, which boost mood by increasing the amount of these neurotransmitters in the brain. Gene Ontology (GO) annotations also relate it to oxidoreductase activity and primary amine oxidase activity. ALOX5A is a protein that localizes to the plasma membrane. Inhibitors of its function impede the translocation of 5-lipoxygenase from the cytoplasm to the cell membrane and inhibit 5-lipoxygenase activation. It is related to pathways such as neurotransmitter clearance and oxidation by cytochrome ^[48]. IGFBP6 is one of the Insulin-like Growth Factor family ^[5], which is a group of growth-promoting peptides. These peptides are secreted by cells that are widely distributed in the human liver, kidney, lung, heart, brain and intestine. IGFBP6 has a strong growth-promoting effect and is often considered an important autocrine and paracrine signaling molecule during the development of the central nervous system. Aif1, also known as Iba1, is a calcium-binding protein that is widely used as a marker of microglia. It is produced by activated monocytes and microglia and selectively expressed in microglia and macrophages ^[49]. TRPM3 belongs to the TRP family and is a pain-sensing pathway ^[4]. TRPM3 is activated by thermal and chemical ligands, such as the neurosteroid pregnenolone sulfate (PregS) and the synthetic ligand CIM0216. Upon activation, it is highly permeable to calcium ions. Some studies have shown that TRPM3 is strongly associated with inflammatory and neuropathic pain and plays a key role in nociceptive temperature perception, particularly mechanical and thermal pain.

The study found that chronic stress has a substantial effect on the frequency and advancement of BD in various aspects, including the production and breakdown of monoamine oxidase, neuroinflammation, ion permeability and pain perception. Additionally, it identified five potential biomarkers that could function as both diagnostic indicators and targets for the treatment of BD. In this study, we emphasize the crucial necessity of examining the impact of genetic factors in chronic stress on both BD and other mental health disorders. From this standpoint, biomarkers linked to persistent stress show potential as valuable instruments for diagnosing BD and influencing therapeutic strategies. To gain a deeper understanding of the connection between chronic stress and BD, it would be beneficial to verify the accuracy of these biomarkers in both animal models and clinical trials. In addition, it would be advantageous to develop novel treatments that focus on these biomarkers.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was funded by the General Program of National Natural Science Foundation of China (82174278).

Conflict of interest

The author has no conflict of interest to declare. The funder did not participate in the study's design, data collection, analysis, interpretation, manuscript writing or the decision to publish the results.