From 2021 to 2023, a total of 14 Caucasian patients (seven male and seven female) age 18 or older with acute suicidal ideation were recruited from the Clinic for Psychiatry and Psychotherapy II at Ulm University, located at the district hospital in Günzburg, Bavaria, Germany. The study was conducted in accordance with the Declaration of Helsinki and received approval from the ethical review board of the Bavarian State Medical Association (Nr. 21007). Patient recruitment and blood collection were performed within the first 24 h of hospital admission by two independent doctoral students. For the purposes of this study, “acute suicidal ideation” was defined as the presence of active suicidal thoughts accompanied by an acute intent to commit suicide, with symptom severity comparable to ICD-10 diagnosis F32.2 and requiring inpatient treatment. Patients presenting with such an acute symptom severity were identified and referred to the doctoral students by on-duty clinical staff.

Inclusion criteria comprised ICD-10 diagnoses F32, F33, F34, F38, or F39, as well as sufficient proficiency in German.

Exclusion criteria included any diagnosis within the schizophrenia spectrum, intellectual disability, bipolar disorder, eating disorders, substance use disorders (SUDs), as well as cancer, chronic inflammatory diseases, epilepsy, diabetes, coronary heart disease, obesity, and pregnancy. Although posttraumatic stress disorder (PTSD) was not explicitly listed as an exclusion criterion, none of the participants had a documented PTSD diagnosis in their medical records. With regard to personality disorders, only one male and one female participant had a diagnosis of borderline personality disorder. Given the very low numbers and the equal distribution across sexes, it is unlikely that borderline personality disorder would account for the observed differences, and it was therefore not included as a covariate in the analysis.

Additionally, 14 sex- and ethnicity-matched healthy controls were recruited via public notices, social media advertisements, and personal contacts, following the same exclusion criteria. For specific sociodemographic details, see Table 1. In Table 1, the term “smoker” refers to participants who answered “yes” when asked whether they consume cigarettes. The exact number of cigarettes smoked per day was not assessed.

Blood samples were collected in PAXgene^® Blood RNA tubes and stored according to manufacturer’s protocol at −80 °C. RNA was then extracted using the MagNA Pure 96 Cellular RNA Large Volume Kit following the manufacturer’s instructions. RNA concentration and quality were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States), and only samples with an A260/A280 ratio greater than 1.8 were included in the experiment. Complementary DNA (cDNA) was synthesized using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, United States). Quantitative PCR was performed using the TaqMan Gene Expression Master Mix (Applied Biosystems) with the following TaqMan gene expression assays: Interleukin 1 beta (Il1b; Hs01555410_m1), Interferon gamma receptor 2 (Ifgnr2; Hs00194264_m1), and TNF alpha-induced protein 6 (Tnfai6; Hs00200180_m1). Glyceraldehyde 3-phosphate dehydrogenase (Gapdh; Hs02786624_g1) was used as an endogenous control. All three genes were selected based on their ranking among the most differentially expressed genes. Reactions were conducted using the Real-Time 7900 Fast system (Applied Biosystems). Relative quantification was performed using the 2^-ΔΔCT method, and gene expression changes were expressed as fold change relative to the corresponding control group (male or female healthy volunteers).

Microarray analyses were performed using 200 ng total RNA as starting material and 5.5 µg ssDNA per hybridization (GeneChip Fluidics Station 450; Affymetrix, Santa Clara, CA). The total RNAs were amplified and labeled following the Whole Transcript (WT) Sense Target Labeling Assay (http://www.affymetrix.com). Labeled ssDNA was hybridized to Human Clariom S Affymetrix GeneChip arrays (Affymetrix, Santa Clara, CA). The chips were scanned with an Affymetrix GeneChip Scanner 3000 and subsequent images analyzed using Affymetrix^® Expression Console™ Software (Affymetrix). A transcriptome analyses was performed using BRB-ArrayTools developed by Dr. Richard Simon and BRB-ArrayTools Development Team (http://linus.nci.nih.gov/BRB-ArrayTools.html), as published previously (Zoller et al., 2017). Raw feature data were normalized and log₂ intensity expression summary values for each probe set were calculated using robust multiarray average (Irizarry et al., 2003).

Genes showing minimal variation across the set of arrays were excluded from the analysis. Genes whose expression differed by at least 1.5-fold from the median in at least 20% of the arrays were retained.

We identified genes that were differentially expressed among the two classes using a two-sample t-test. Genes were considered statistically significant if their p value was less than 0.05 and displayed a fold change between the two groups of at least 1.5-fold. We used the Benjamini and Hochberg correction to provide 90% confidence that the false discovery rate was less than 10% (Benjamini and Hochberg, 1995). The false discovery rate is the proportion of the list of genes claimed to be differentially expressed that are false positives.

We developed models for utilizing gene expression profile to predict the class of future samples. We developed models based on the Compound Covariate Predictor (Radmacher et al., 2002), Diagonal Linear Discriminant Analysis (Dudoit et al., 2002), Nearest Neighbor Classification (Dudoit et al., 2002), and Support Vector Machines with linear kernel (Ramaswamy et al., 2001). The models incorporated genes that were differentially expressed at the 0.001 significance level as assessed by the random variance t-test (Wright and Simon, 2003). We estimated the prediction error of each model using leave-one-out cross-validation (LOOCV) as described by (Simon et al., 2003). For each LOOCV training set, the entire model building process was repeated, including the gene selection process. We also evaluated whether the cross-validated error rate estimate for a model was significantly less than one would expect from random prediction. The class labels were randomly permuted and the entire LOOCV process was repeated. The significance level is the proportion of the random permutations that gave a cross-validated error rate no greater than the cross-validated error rate obtained with the real data. 1000 random permutations were used.

To identify associated biological processes, as defined by Gene Ontology annotation, we used the GoMINER analysis tool (Zeeberg et al., 2003). This package allows the automatic analysis of multiple microarrays and then integrates the results, across all of them, to find the GO categories that were significantly over- or under-represented.

Differential gene expression data were generated by downloading count data from the dataset GSE247998 (Sun et al., 2024), downloaded from Gene Expression Omnibus, and analyzed using the R-package Desq2 (Love et al., 2014) version 1.44 under R version 4.4. Genes with read counts below 10 were excluded.

To investigate the biological pathways and functional gene sets enriched in the blood transcriptomic profiles of female patients with MDD and suicidal ideation, we performed Gene Set Enrichment Analysis (GSEA) using GSEA Version 4.3.3 (https://www.gsea-msigdb.org/gsea/index.jsp) and genesets Hallmark (H), curated genesets (C2) and immunologic genesets (C7) (Subramanian et al., 2005; Mootha et al., 2003) on the preprocessed gene expression data.

GSEA was run using the following parameters: gene set size filter of 15–500 genes, 1,000 permutations, and the weighted enrichment statistic. Statistical significance was determined based on a false discovery rate (FDR) q-value <0.25, as recommend in the documentation (https://docs.gsea-msigdb.org/#GSEA/GSEA_User_Guide/#interpreting-gsea-results).

In recent years, a powerful method for systematic analysis called weighted gene co-expression network analysis (WGCNA) has been widely applied in bioinformatic analysis of various diseases (Zhang and Horvath, 2005). The robust co-expression network is able to cluster genes with similar expression patterns into modules to identify the underlying biological functions. To identify gene co-expression modules associated with suicidal ideation, we used this method on the transcriptomic data using the WGCNA R-package (R version 4.3) developed by P. Langfelder, and S. Horvath (Langfelder and Horvath, 2008).

STRING database 12.0 (Szklarczyk et al., 2023), https://string-db.org/) was used to identify possible protein-protein interactions between differentially regulated genes as well as between WGCNA hub genes. Default settings (Full STRING network; network edges: evidence; all active interaction sources and medium confidence) were used for both analyses.

We first analyzed male and female acutely suicidal patients with MDD together and identified 87 differentially expressed genes (91 probesets) (p-value <0.05, FDR < 0.1, |FC| > 1.5x), a threshold applied across all subsequent comparisons. The complete gene list is provided in Supplementary Table S1.

To determine whether transcriptional profiles could effectively distinguish suicidal patients with MDD from healthy controls, we applied several class prediction methods and observed a striking difference between male and female patients. In the male group, only two of seven algorithms (Compound Covariate Predictor and Diagonal Linear Discriminant Analysis) correctly predicted a maximum of 9 out of 14 samples. In contrast, three algorithms in the female group correctly predicted all samples, while the remaining four accurately classified 13 out of 14 samples. This result was further supported by a PCA analysis of the gene expression data.

As shown in Figure 1, there is only a clear separation between the female suicidal ideation group and their corresponding healthy control group, but not between the corresponding male subgroups. Based on these findings, we proceeded with sex-specific differential gene expression analyses.

This approach revealed no differentially expressed genes in suicidal men with MDD compared to healthy male controls. However, in women, we identified 671 differentially expressed probesets corresponding to 665 genes, with 304 downregulated and 367 upregulated. A list of the up- and downregulated genes is available in Supplementary Table S2.

To ensure these findings were not driven by differences between the male and female healthy control groups, we compared their gene expression profiles. This analysis yielded only 23 differentially expressed genes, all located on the sex chromosomes. Additionally, when comparing suicidal male and female MDD patients, we identified only 40 differentially expressed genes, 26 of which were XY-chromosomal (provided as Supplementary Table S3). These results strongly indicate that the observed gene expression differences are specifically associated with MDD-driven suicidality and sex-dependent to their respective healthy controls.

Gene Ontology enrichment analysis of the differentially expressed genes in suicidal women with MDD revealed significant enrichment in immune system and inflammatory pathways. The top 20 GO terms are presented in Supplementary Table S4.

To identify and visualize potential protein-protein interactions within the set of differentially regulated genes in the female MDD patient group with acute suicidal ideations, we used the STRING database. The network as identified by STRING [Figure 2; Supplementary Figure S1 (high resolution)] revealed a complex network of interactions mainly centered around two clusters of genes consisting of immunological genes such as SYK, FCGR1A or LCK or genes involved in ribosomal biosynthetic processes such as RPS24, RPL27a, NSA2 or BMS1.

Gene Set Enrichment Analysis (GSEA) of the female gene expression data identified 12 significantly enriched gene sets (FDR < 0.25) when compared to the Hallmark gene sets (Table 2). The top three enriched pathways were TNFA signaling via NFKB, IL6-JAK-STAT3 signaling, and inflammatory response. Among the identified pathways, TNF-alpha Induced Protein 6 (TNFAIP6) emerged as the top-ranked gene in both TNFA signaling via NFKB and inflammatory response pathways. Other notably differentially expressed genes were Bestrophin 1 (BEST1), Protein Tyrosine Phosphatase Receptor Type r (PTPRE), Plasminogen Activator/Urokinase Receptor (PLAUR), Aquaporin 9 (AQP9), the inward-rectifying potassium channel KCJN2 and alpha-2-macroglobulin (A2M).

Since TNFAIP6 was found to be significantly overexpressed in the suicidal group, it was together with Interleukin 1-beta (IL1B) and Interferon Gamma Receptor 2 (IFNgR2) (Fritz et al., 2018) selected for validation using quantitative PCR. Consistent with the initial array findings, all three were not upregulated in the male population. However, IL1B, IFNgR2, and TNFAIP6 exhibited significant differential expression between suicidal women with MDD and their respective healthy controls (see Figure 3).

Given the overlap between genes in the top five GSEA-enriched pathways, we constructed a Venn diagram to illustrate shared gene expression patterns (Figure 4A; Table 3). Interleukin-6 (IL6) emerged as the top shared gene, present in all five pathways.

To further validate our sex-specific findings, we repeated the same analysis with the male dataset, especially since GSEA analysis is not limited to a set of differentially expressed genes that could be biased by filtering conditions on fold-change and p-value. In support of our sex-specific findings, not a single pathway was identified (FDR < of 0.25) in both the hallmark and curated dataset collection (C2) in contrast to 12 and 61 genesets in the female dataset.

To further explore the regulatory networks within our dataset, we performed weighted gene co-expression network analysis (WGCNA). Similar to GSEA, WGCNA is not restricted to a predefined set of differentially expressed genes but instead utilizes the entire gene expression dataset. After testing multiple candidate powers, we selected a soft threshold power of 12 based on the approximate scale-free topology criterion. This resulted in the identification of 13 distinct WGCNA modules (Figure 4B), with module sizes ranging from 52 to 2,076 genes (Figure 4C).

To identify modules significantly associated with suicidal ideations (correlation >0.5, p < 0.001), we correlated eigengenes with both sex and suicide attempt status. The resulting correlation plot (Figure 4D) again revealed a striking difference between the male and female MDD-patients with suicidal ideation. In contrast to the male group, where no highly significant correlations were observed, the female group exhibited a highly significant module (turquoise) with a strong correlation (r = 0.67, p < 0.001). Interestingly, when male and female data were combined in the analysis, the correlation strength for this module decreased (from r = 0.67 to r = 0.51) alongside its statistical significance (from p < 0.001 to p < 0.05). This further supports the necessity of analyzing male and female data separately in studies of gene expression related to suicidality.

Next, we identified key regulatory genes, or “hub genes,” within the turquoise module using gene significance (GS > 0.3) and module membership (MM > 0.8) criteria (Liang et al., 2020). A total of 363 hub genes (Supplementary Table S8) were identified, which were then subjected to gene set enrichment analysis using the Hallmark collection. The results closely mirrored those from the GSEA of differentially expressed genes, with significant enrichment in TNF alpha signaling via NFKB, inflammatory response, and complement pathways (Table 4). To further investigate interactions between these hub genes, we conducted a STRING analysis. The resulting gene interaction network (Figure 5; Supplementary Figure S5) revealed two major central nodes: IL1B (inflammatory response) and GRB2 (RAS-ERK signaling pathway).

To test whether we could reproduce our results on a previously published dataset, we used data from Sun et al. (2024), who compared brain and blood transcriptome profiles with respect to common genetic pathways in suicidal ideation and suicide. Since the authors of this study did not publish differential expression data on their samples, we analyzed the count data as provided in the corresponding GEO dataset using the R-package DeSeq 2. Surprisingly, we see only a small set of six genes (one pseudogene, one lncRNA and four genes related to hemoglobin biosynthesis) when comparing males with suicide attempt (SA, eight samples) to male healthy controls (11 samples) and not a single gene (|FC| > 1.5x, FDR < 0.05) when comparing females with SA (15 samples) to female healthy controls (15 samples). On the other hand, we found 439 differentially expressed genes when comparing the female suicide attempter group with the male suicide attempter group. This striking difference is also visible in a PCA plot of the samples (Figure 6), which clearly separates these two groups but not the sex-specific samples with or without SA. Since this study by Sun et al. was focused on a gene network approach we extracted the genes of the 18 significantly enriched modules from their Ingenuity Pathway analysis (IPA, data available in the supplement to this study) and used this set of 450 genes for a GSEA to identify possible overlaps with our data from the analysis of the DE genes within the HALLMARK collection. As it turned out, the overlap consisted of three genesets, namely, INFLAMMATORY_RESPONSE, COMPLEMENT and APOPTOSIS.

When we searched for corresponding GO groups within the IPA data using DAVID (Huang et al., 2009; Sherman et al., 2022), the results were similar to our own data, with the highest scoring GO groups containing “positive regulation of NF-kappaB transcription factor activity”, “positive regulation of inflammatory response”, or “positive regulation of interleukin 1 beta production”. Thus, despite the fact that we could not reproduce the sex-specific differences in gene expression between suicidal patients and healthy individuals, the differential gene expression between male and female subjects with suicide attempts leads to the same biological pathways.

It is not possible to rule out the possibility that our findings are influenced by the relatively small sample size in this study. We anticipate that differentially expressed genes will also be identified in male participants; however, these differences may be less pronounced than those observed in females. Therefore, we strongly encourage researchers conducting similar studies to assess the reproducibility of these findings by incorporating sex-specific analyses within their datasets.

Alternatively, the observed discrepancy between male and female patients may be attributable to the composition of our sample group and could stem from genetic variability among individuals, which may obscure the detection of statistically significant expression changes. This increased variability likely has a comparatively smaller impact on gene set enrichment analysis, potentially explaining the consistent findings in pathway-level analyses. Nevertheless, if this is the case, our results underscore not only the importance of considering sex-specific differences in psychiatric research but also highlight the need for a more personalized approach to patient care, moving away from the current one-size-fits-all model.

Additionally, because adverse childhood experiences (ACEs) are associated with depression and have been shown to influence gene expression (Parel and Peña, 2022), we cannot rule out their potential impact on our findings, as ACEs were not assessed in this study.

Finally, this study refrains from reporting specific biomarkers for suicidality, as it does not include a comparison group of MDD patients who are not acutely suicidal, which represents a notable limitation. Nevertheless, its primary and well-founded contribution lies in highlighting the striking sex differences observed, thereby encouraging future research to confirm whether some of the identified genes may serve as genuine sex-specific biomarkers of suicidality.

Supplementary Tables S1–S3 present lists of differentially expressed genes (DEGs) identified through transcriptomic analysis. Table 1 includes DEGs (p < 0.05, FDR < 0.1, |FC| > 1.5x) in the combined male and female cohort compared to healthy controls, while Table 2 focuses specifically on DEGs within the female subset and Table 3 lists the DEGs between the male and female samples with suicide attempts. Supplementary Table S4 lists Gene Ontology (GO) categories enriched among DEGs in the female subgroup, as identified via GoMiner analysis. Supplementary Tables S5–S7 detail the leading-edge genes from gene set enrichment analysis (GSEA) for three hallmark pathways: TNFA_SIGNALING_VIA_NFKB, IL6_JAK_STAT3_SIGNALING, and INFLAMMATORY_RESPONSE. Supplementary Table S8 lists hub genes (Gene Significance >0.6, Module Membership >0.8) derived from Weighted Gene Co-expression Network Analysis (WGCNA) of the female sample group. Supplementary Figures S1 and S5 display high resolution STRING network analyses of DEGs in female suicide attempters (vs controls) and hub genes within the highly significant WGCNA “turquoise” module, respectively. Supplementary Figures S2–S4 provide enrichment plots for the hallmark pathways examined in GSEA: TNFA_SIGNALING_VIA_NFKB, IL6_JAK_STAT3_SIGNALING, and INFLAMMATORY_RESPONSE.