Transcriptome Profiling of the Dorsomedial Prefrontal Cortex in Suicide Victims

The default mode network (DMN) plays an outstanding role in psychiatric disorders. Still, gene expressional changes in its major component, the dorsomedial prefrontal cortex (DMPFC), have not been characterized. We used RNA sequencing in postmortem DMPFC samples to investigate suicide victims compared to control subjects. 1400 genes differed using log2FC > ±1 and adjusted p-value < 0.05 criteria between groups. Genes associated with depressive disorder, schizophrenia and impaired cognition were strongly overexpressed in top differentially expressed genes. Protein–protein interaction and co-expressional networks coupled with gene set enrichment analysis revealed that pathways related to cytokine receptor signaling were enriched in downregulated, while glutamatergic synaptic signaling upregulated genes in suicidal individuals. A validated differentially expressed gene, which is known to be associated with mGluR5, was the N-terminal EF-hand calcium-binding protein 2 (NECAB2). In situ hybridization histochemistry and immunohistochemistry proved that NECAB2 is expressed in two different types of inhibitory neurons located in layers II-IV and VI, respectively. Our results imply extensive gene expressional alterations in the DMPFC related to suicidal behavior. Some of these genes may contribute to the altered mental state and behavior of suicide victims.


Introduction
The dorsomedial prefrontal cortex (DMPFC) represents one of the nodes in the defaultmode network (DMN). It includes the superior frontal and the paracingulate gyrus (also called the dorsal anterior cingulate cortex) on the medial surface of the frontal cortex (see Methods). The DMPFC is involved in a variety of functions, combining sensory signals from the outside world transferred by the lateral prefrontal (cognitive) and the orbitofrontal (motivational) cortices, as well as viscerosensory signals through the insula and the anterior cingulate cortex [1,2]. The DMPFC prepares executive motor programs for the supplementary and presupplementary motor cortices [3][4][5]. The functions of the DMPFC include activation of working memories, decision making and planning [6][7][8][9][10]. The DMPFC is linked to social processes [11], emotional processing [12], novelty processing, dynamic contextual updating [13] and conscious threat appraisal [14,15].
The other node in the DMN is present in the medial posterior parietal cortex, including the precuneus and the posterior cingulate cortex [16][17][18]. The precuneus plays a significant After counting the reads aligning to each gene, we used DESeq2 to analyze gene expression in suicide victims compared to matched controls. We detected expression of 19,692 protein-coding genes in the DMPFC, out of which 1400 were differentially expressed between the two groups, 1262 genes were downregulated and 138 genes were upregulated at the adjusted p-value < 0.05 and log2 FC > ±1 (Supplementary File S1). All samples were compared to one another using Euclidean distance. The correlation matrix based on all gene expression data indicated that samples mostly cluster by individual and diagnostic group with three control samples showing some similarity to the suicide samples ( Figure 1A). We compared the age, sex and PMI between the three controls clustered with the suicide victims and those who were not clustered using a t-test. We found no significant differences between the two groups of control individuals for any of these covariates (age: p = 0.51; PMI: p = 0.39; sex ratio: p = 0.67). Therefore, we do not assume any systematic confounding variables that can distort the results. Hierarchical clustering of suicide and control samples using 1400 differentially expressed genes revealed that these genes also distinguished the two groups ( Figure 1B). Among the DEGs, 1262 genes were downregulated and 138 were upregulated in suicide victims, as shown in the volcano plot ( Figure 2).

Functional Annotation and Classification of the DEGs
To understand how the DEGs from the DMPFC relate to suicidal behavior, Gene Ontology (GO) classification was performed. Figure 3A shows the top 10 down-and upregulated DEGs. The top three downregulated genes were CSF3, IL1R2 and SOCS3. They all play a critical role in inflammation. The top three upregulated genes were P2RX2, ATP4A and BX276092.9. P2RX2 and ATP4A are found in the plasma membrane and both of them have a role in ATP binding and ion transport. BX276092.9 is a barely characterized protein. Additional GO functional enrichment analysis identified significant differences in more than 500 functional pathways of downregulated genes and 37 functional pathways of upregulated genes at adjusted p-value < 0.05 (Supplementary File S2). A comparison with the Allen Human Brain Atlas (AHBA) expression database showed that some of the altered   16 RNA-seq samples based on the Euclidean distance. A heatmap of the distance matrix shows the similarities and dissimilarities between the samples as calculated from the rlog transformation of the count data. (B) Hierarchical clustering heatmap of Pearson correlation coefficients between suicide and control individuals. The heatmap represents the top 500 protein-coding gene expression for the different groups in columns, with a dendrogram presented at the top of the heatmap. Proteins were selected based on their variability between all samples. Counts were log-transformed, normalized, and used for clustering based on similarity in expression patterns. The x-axis represents the sample. The y-axis represents the top 500 differentially expressed genes. The color represents the log10 transformed gene expression level. The scale shows the level of expression: the yellow color means high, the blue color means low, while the black color represents medium expression level. shows the similarities and dissimilarities between the samples as calculated from the rlog transformation of the count data. (B) Hierarchical clustering heatmap of Pearson correlation coefficients between suicide and control individuals. The heatmap represents the top 500 protein-coding gene expression for the different groups in columns, with a dendrogram presented at the top of the heatmap. Proteins were selected based on their variability between all samples. Counts were logtransformed, normalized, and used for clustering based on similarity in expression patterns. The xaxis represents the sample. The y-axis represents the top 500 differentially expressed genes. The color represents the log10 transformed gene expression level. The scale shows the level of expression: the yellow color means high, the blue color means low, while the black color represents medium expression level.
Among the DEGs, 1262 genes were downregulated and 138 were upregulated in suicide victims, as shown in the volcano plot ( Figure 2).

Figure 2.
Visualization of RNA-seq results. Volcano plot showing the log2 fold change (log2FC) of gene expression and the statistical significance of the differential expression (DE) analysis performed between suicide and control individuals. The x-axis represents the log2 fold change of genes, while the y-axis represents the −log10 of the corrected p-values (padj) for the different pairs of conditions. Each dot represents a gene and the colored areas represent the DEGs that met the following selection criteria: log2FC of at least ±1 (log2FC ≥ 1 or ≤−1) and Benjamini and Hochberg-adjusted pvalue < 0.05. Upregulated genes are shown in red, while the downregulated ones are blue. The top 40 significant DEGs are labeled.

Functional Annotation and Classification of the DEGs
To understand how the DEGs from the DMPFC relate to suicidal behavior, Gene Ontology (GO) classification was performed. Figure 3A shows the top 10 down-and upregulated DEGs. The top three downregulated genes were CSF3, IL1R2 and SOCS3. They all play a critical role in inflammation. The top three upregulated genes were P2RX2, ATP4A and BX276092.9. P2RX2 and ATP4A are found in the plasma membrane and both of them have a role in ATP binding and ion transport. BX276092.9 is a barely characterized protein.
Additional GO functional enrichment analysis identified significant differences in more than 500 functional pathways of downregulated genes and 37 functional pathways of upregulated genes at adjusted p-value < 0.05 (Supplementary File S2). A comparison with the Allen Human Brain Atlas (AHBA) expression database showed that some of the altered genes are expressed in specific cell types. For example, one of the top downregulated genes, the MT-ND5, is produced in FEZF2 and RORB expressing excitatory and SST+ Figure 2. Visualization of RNA-seq results. Volcano plot showing the log2 fold change (log2FC) of gene expression and the statistical significance of the differential expression (DE) analysis performed between suicide and control individuals. The x-axis represents the log2 fold change of genes, while the y-axis represents the −log10 of the corrected p-values (padj) for the different pairs of conditions. Each dot represents a gene and the colored areas represent the DEGs that met the following selection criteria: log2FC of at least ±1 (log2FC ≥ 1 or ≤−1) and Benjamini and Hochberg-adjusted p-value < 0.05. Upregulated genes are shown in red, while the downregulated ones are blue. The top 40 significant DEGs are labeled.
The top three enriched GO terms within the biological process, molecular function, and cellular component categories were visualized (Figure 3B,C). A number of pathways related to cell surface receptor signaling (GO:0007166) and growth factor binding (GO:0019838) were enriched in downregulated genes in suicidal individuals. In contrast, pathways pertaining to synaptic signaling (GO:0099536) and sodium channel activity (GO:0005272) were enriched in upregulated genes of suicidal individuals ( Figure 3B,C). To further identify the pathways involved in suicidal behavior, pathway analysis of DEGs was performed using pathway classifications from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome resources (Statistics of Pathway Enrichment). The top three enriched KEGG and Reactome pathways are shown in Figure 3D. A total of 23 significantly enriched downregulated KEGG pathways and 25 significantly enriched downregulated Reactome pathways were found in downregulated genes, while there was one significantly enriched upregulated KEGG pathway and two upregulated Reactome pathways at adjusted p-value < 0.05. The complete list of significantly enriched pathways can be found in Supplementary File S2. shows the terms and the y-axis shows the enrichment p-values on the −log10 scale. Each circle on the plot corresponds to a single term. Circles are colored according to the origin of annotation and size scaled according to the total number of genes annotated to the corresponding term. The locations on the x-axis are fixed. Terms from the same GO subtree are located closer to each other on the x-axis, which helps to highlight different enriched GO sub-branches making plots from different queries comparable. The top 3 terms from each GO category are indicated with numbers on the plot. The corresponding statistics of annotation are listed in Supplementary File S2. (C) GO classification and scatter plot-enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs. Bar plots reporting the top 9 significantly enriched GO terms related to down-and upregulated genes with the adjusted p-value < 0.05, respectively. The table describes, for each GO term, the number of mapped annotated genes in the reference data set and its −log10 p-value. The color code is proportional to the three ontologies. (D) Dot plots of top 3 pathway annotations from KEGG and Reactome pathways illustrate the distributions of gene sets among up-and downregulated differentially expressed genes. The x-axis represents the gene ratio. The y-axis represents the pathway term. Gene ratio refers to the value of enrichment, which is the ratio of DEGs annotated in the pathway to the total gene amount annotated in the pathway. The larger the value, the more significant the enrichment. The circle size indicates the DEG number associated with each significant pathway. The color indicates the adjusted p-value, the lower p-value indicates the more significant enrichment.

Protein-Protein Interaction Analysis of DEGs and Identification of Key Genes
Computational methods analyzing protein protein interaction (PPI) networks are useful for understanding the biological meaning of gene expressional alterations. Therefore, we used the STRING online database (http://string-db.org, (accessed on 11 November 2021)) to construct the PPI networks of down-and upregulated DEGs, which were then transferred into the Cytoscape software. Analysis of top DEGs was achieved using the degree method in cytoHubba where the top 10 down-and upregulated DEGs were identified ( Figure 4A,C). Enrichment analysis of the top 10 downregulated DEGs s ( Figure 4A) through stringApp revealed that these genes are mainly associated with astrocyte differentiation (GO:0048708), regulation of MAPK cascade (GO:0043408) and cell surface receptor signaling pathway (GO:0007166) ( Figure 4B); meanwhile, the top 10 upregulated DEGs ( Figure 4C) were associated with regulation of membrane potential (GO:0042391), nervous system process (GO:0050877) and cell-cell signaling (GO:0007267) ( Figure 4D). The complete list of genes and their associated pathways can be found in Supplementary File S3. Using the AHBA database for comparison, we found that most of the top downregulated genes are abundantly expressed in astrocytes. In turn, the top 10 upregulated genes are expressed in excitatory and inhibitory neurons and can be also found in oligodendrocytes (Supplementary File S6).

Validation of RNA-Seq Data
To verify changes in gene expression associated with suicidal behavior, we performed qRT-PCR. We have chosen fourteen functionally relevant genes (six up-and eight downregulated) from the DEGs, which have been potentially implicated in depressive behavior based on the literature: GRIK1, GRIK2, GRM2, NRGN, SYT5 and NECAB2 as upregulated genes and AQP1, ITPKB, ITGB4, SLCO2B1, GJA1, PRKCH, GLUL and S100B as downregulated genes. As shown in Supplementary Figure S2, results of quantitative PCR analysis of the selected genes confirmed the RNA-seq results. Thus, we verified the increased expression of GRIK1, GRIK2 and GRM2 involved in glutamate signaling, as well as NECAB2, while the reduced expression of astrocyte-related genes, such as S100B, AQP1, GJA1, GLUL, ITPKB and ITGB4 in suicide victims (  On the graph, y-axis represents the significant enrichment terms. Each bar describes the number of mapped annotated genes in the reference data set while the x-axis and color gradient indicate the significance (−log10 p-value) in each category.

Validation of RNA-seq Data
To verify changes in gene expression associated with suicidal behavior, we performed qRT-PCR. We have chosen fourteen functionally relevant genes (six up-and eight downregulated) from the DEGs, which have been potentially implicated in depressive behavior based on the literature: GRIK1, GRIK2, GRM2, NRGN, SYT5 and NECAB2 as upregulated genes and AQP1, ITPKB, ITGB4, SLCO2B1, GJA1, PRKCH, GLUL and S100B as downregulated genes. As shown in Supplementary Figure S2, results of quantitative PCR analysis of the selected genes confirmed the RNA-seq results. Thus, we verified the increased expression of GRIK1, GRIK2 and GRM2 involved in glutamate signaling, as well as NECAB2, while the reduced expression of astrocyte-related genes, such as S100B, AQP1, GJA1, GLUL, ITPKB and ITGB4 in suicide victims (  On the graph, y-axis represents the significant enrichment terms. Each bar describes the number of mapped annotated genes in the reference data set while the x-axis and color gradient indicate the significance (−log10 p-value) in each category.

Depression-Focused Gene Set Enrichment
The top 10 down-and upregulated DEGs based on fold change as well as the PPI and co-expression network analyses were investigated in the DisGeNET database to search for significant enrichments associated with depression and other mental disorders using the disease classification terms "Mental Disorders" and "Behavior and Behavior Mechanisms". We found 95 enriched categories associated with more than one gene. We found 24 genes that have been previously associated with depression, including CSF3, IL1R2, IL6, SERPINA3, MT-ND5, MT-ND6, P2RX2, ATP4, NOTCH1, EGFR, TLR4, STAT3, ERBB2, PROC, CARTPT, GRIK2, CACNA1G, GABRD, BCL2, CHDH, NRGN, SERPINF1, CCK and CRHR1. Genes associated with depressive disorders, schizophrenia and impaired cognition were strongly enriched in the DMPFC. Furthermore, we found that the upregulated genes CCK, GRIK2 and CRHR1 are affected in suicide (Supplementary File S4). Although the database does not contain information specific to DMPFC, our data suggest that these genes may be associated with suicide in this brain region, too.

Co-Expression Network Analysis and Hub Gene Screening in the DMPFC of Suicidal Individuals
We built two independent gene co-expression networks for down-and upregulated genes to identify gene co-expression communities (Supplementary File S5). Analysis was conducted using highly correlated gene pairs with Pearson correlation > 0.9 and Padj < 0.01.
In the networks, the DEGs were represented by nodes and pairwise co-expression relationships between genes by edges. The co-expression network of downregulated DEGs revealed the top 10 highest ranked hub genes as SORBS3, RHOC, S100A16, SZRD1, TRIP6, AHNAK, GRIN2C, CHDH, MAPKAPK2 and FTL with higher degree and betweenness centrality than others, indicating a more critical role played by them in the network. The network was further subdivided into 21 functional clusters composed of a total of 710 downregulated DEGs (nodes) (Supplementary File S5). There was a significant difference in the biological processes as different modules were enriched ( Figure 5A). For example, the biggest module Cluster 1 (1) was significantly enriched using biological terms of the GO classification, such as "gliogenesis" (GO:0042063) and "response to cholesterol" (GO:0070723). The second biggest module, Cluster 2 (2) was significantly enriched using biological terms, such as "regulation of epithelial cell differentiation" (GO:0030856) and "type I interferon signaling pathway" (GO:0060337) ( Figure 5A). The co-expression network of upregulated DEGs revealed the top 10 most ranked hub genes as CALY, SOHLH1, CPNE9, NRGN, MAST1, GRIK1, SYT5, TRPM2, BX276092.9 and MYO15A. The network was further subdivided into six functional clusters of a total of 97 upregulated DEGs (Supplementary File S5), where Cluster 1 (1) was significantly enriched using biological terms "neurotransmitter receptor internalization" (GO:0099590) and "biotin metabolic process" (GO:0006768). Cluster 2 (2) was significantly enriched using biological terms, such as "sodium channel activity" (GO:0005272) and "neuropeptide hormone activity" (GO:0005184) ( Figure 5B). The hub proteins of clusters were specifically identified as most ranked, based on their degree and betweenness centrality (Supplementary File S5). For instance, the hub protein of downregulated co-expression network in Cluster 1 is the SORBS3 encoding vinexin protein, which is implicated in promoting the upregulation of actin stress fiber formation. An additional hub protein of the downregulated network in Cluster 2 is RHOC encoding the Rho-related GTP-binding protein, which is known to regulate the signal transduction pathway linking plasma membrane receptors to actin stress fibers. The hub protein of upregulated co-expression network in Cluster 1 is CALY encoding the neuron-specific vesicular protein calcyon, which interacts with the clathrin light chain A and stimulates clathrin-mediated endocytosis. Another hub protein from the upregulated network in Cluster 2 is SOHLH1 encoding the spermatogenesis-and oogenesis-specific basic helix-loop-helix-containing protein 1, which was implicated in transcriptional regulation of both male and female germline differentiation. The entire coexpression networks and the annotation information of the top three clusters are available in Supplementary File S5.
Using the CytoHubba plugin in Cytoscape, we determined the top 10 hub genes of down-and upregulated co-expression networks according to their highest Maximal Clique Centrality (MCC) scores ( Figure 6A,B). We constructed a protein-protein interaction network (PPI) for all the top 10 down-and upregulated hub genes using the STRING database, out of which two of the downregulated and four of the upregulated genes were connected to each other ( Figure 6C). In GO and Reactome pathway enrichment analyses of hub genes, we identified enriched gene sets related to developmental processes, such as actin filament organization, neuron maturation or anterograde axonal transport ( Figure 6D). We found similar results when analyzing cell type expressing features of the top hub genes by using the AHBA database, where most of the downregulated hub genes are expressed in astrocytes, while the upregulated ones can be found in excitatory and inhibitory neurons and oligodendrocytes (Supplementary File S6). The top 10 hub genes from the co-expression networks and their functional annotation are listed in Supplementary File S5.

Distribution of NECAB2 in the DMPFC
As a functionally intriguing validated DEG, NECAB2 was selected for analysis of its distribution within the DMPFC. In situ hybridization histochemistry was applied in a control and a suicidal individual to localize the expression of NECAB2 mRNA. The hybridization signal was abundant in layers II-VI ( Figure 7A,B). NECAB2 immunolabeled cells showed the same distribution ( Figure 7C). Immunolabeling also revealed that NECAB2 is present in two different interneuron subtypes. Most of the NECAB2-positive interneurons had small size cell bodies, while some of them demonstrated large soma (Figure 7(C3)). The larger interneurons were located in deep layers; meanwhile, the smaller ones were distributed in the upper layers, predominantly in layer II (Figure 7(C1,C3)). Immunolabeling also revealed that NECAB2 is present not only in the somato-dendritic compartments, but also in axon terminals in neurons (Figure 7(C2,C3)). To determine whether NECAB2 protein is also present in glial cells, we performed combined immunolabeling for Iba1 (a microglia marker) and S100B (an astrocyte marker) with NECAB2. The lack of colocalization indicates that NECAB2 is not expressed in glial cells ( Figure 7D-F Each node refers to a specific Gene Ontology (GO) Biological process term and is grouped based on the similarity of their associated genes. The size of the nodes is equivalent to the statistical power of significance of each term: the larger the node size is, the smaller the p-value is. Node colors represent different functional groups. The emphasized term of each functional group is given by the most significant term of the group.
down-and upregulated co-expression networks according to their highest Maximal Clique Centrality (MCC) scores ( Figure 6A,B). We constructed a protein-protein interaction network (PPI) for all the top 10 down-and upregulated hub genes using the STRING database, out of which two of the downregulated and four of the upregulated genes were connected to each other ( Figure 6C). In GO and Reactome pathway enrichment analyses of hub genes, we identified enriched gene sets related to developmental processes, such as actin filament organization, neuron maturation or anterograde axonal transport ( Figure  6D). We found similar results when analyzing cell type expressing features of the top hub genes by using the AHBA database, where most of the downregulated hub genes are expressed in astrocytes, while the upregulated ones can be found in excitatory and inhibitory neurons and oligodendrocytes (Supplementary File S6). The top 10 hub genes from the co-expression networks and their functional annotation are listed in Supplementary File S5. The larger interneurons were located in deep layers; meanwhile, the smaller ones were distributed in the upper layers, predominantly in layer II ( Figure 7C(C1,C3)). Immunolabeling also revealed that NECAB2 is present not only in the somato-dendritic compartments, but also in axon terminals in neurons (Figure 7(C2,C3)). To determine whether NECAB2 protein is also present in glial cells, we performed combined immunolabeling for Iba1 (a microglia marker) and S100B (an astrocyte marker) with NECAB2. The lack of colocalization indicates that NECAB2 is not expressed in glial cells ( Figure 7D-F).

Discussion
The underlying pathophysiological basis of suicidal behavior remains basically obscure. Neuroimaging studies revealed that committing suicide, similar to MDD, may be a neural network-level disturbance [58,59]. Neuroimaging studies revealed abnormal resting-state functional connectivity (RSFC) between distributed brain areas in MDD. The connection between the medial prefrontal cortex and the medial posterior parietal cortex has been named the resting-state network (RSN) since they have temporal correlation in the absence of a variety of tasks in resting condition [60][61][62][63]. The RSN includes distinct networks, such as the default mode, salience and attention networks [3,9,16,[64][65][66][67]. The default mode network (DMN) has been highlighted in neuroimaging studies, since its discovery has attracted increased interest due to its implications for MDD and suicide-related behav-ior. The DMN is involved through convergent findings of increased resting-state functional connectivity between two core DMN regions, namely the DMPFC and precuneus [68].
Even though the DMPFC is a major component of the DMN, no attention has been assigned to examine transcriptome changes in this region in suicide victims or depressed patients [51,53,57]. We present here the results of RNA-seq analyses in the DMPFC from sixteen subjects, eight with fatal suicide action and eight controls. The number of samples in both groups are sufficiently high to conclude significant results and draw proper conclusions [69]. We identified more than 1000 DEGs. The number of DEGs is in line with previous studies investigating suicide victims [53] and schizophrenic individuals [70]. In our study, a high number of genes decreased their expression level. Since our RNA-seq data fit the DESeq2 model, technical issues are not likely. Therefore, we conclude that suicidal behavior is accompanied by a generally reduced gene expression in the DMPFC.

Conclusions Based on Individual DEGs
Among the DEGs the top 10 downregulated genes were involved in immune response, neurodegeneration, mitochondrial electron transport and cell adhesion. The role of immune dysfunction in suicide victims has been reported in the prefrontal cortex [51] and the insula [71]. It has also been proposed that neurodegeneration and impaired structural neuroplasticity [71,72], as well as mitochondrial dysfunction [53,73], are associated with suicide completion. The top 10 upregulated genes were involved in ATP signaling, protein and vesicle-mediated transport, and were components of membrane and cytoskeleton and regulating protein secretion. Changes in protein and vesicle-mediated transport in suicide have also been observed previously [51,71,74].
In addition to the DEGs with the highest fold changes, we also focused on the individually validated DEGs. For qRT-PCR validation, we selected 14 genes whose alterations were reported in human depression and suicide or rodent depression-like behavior [52,54,[75][76][77][78][79][80]. For ITPKB, GJA1 and GLUL genes, the previously reported alteration in the DLPFC [54] was demonstrated to also occur in the DMPFC. For AQP1 and ITGB4, studied only in rodent models of depression [75,78], the present study was the first to implicate them in human suicidal behavior.
Several genes connected to astrocyte function, such as ATP1A2, ALDH1L1, GFAP, S100B, GJA1 and AQP1, showed decreased expression levels in suicide victims. These observations suggest that genes promoting astrocyte functions may be involved in the pathophysiology of depression and suicidal behavior, as shown previously for depression [81][82][83].

Conclusion Based on the Distribution of NECAB2
NECAB2 is one of the validated genes whose level increased in suicide victims. It was shown to bind specifically to the type 5 metabotropic glutamate receptor (mGluR5) to modulate its function [84]. Our findings, consistent with a previous study [85], provided evidence that the glutamatergic signaling system is involved in the pathogenesis of suicidal behavior. Therefore, it is likely that NECAB2 exerts its role on suicidal behavior by increasing the activity of mGluR5. The AHBA suggests that NECAB2 is expressed in interneurons located in layers II-VI in the human cerebral cortex. The majority of NECAB2 expressing cells are GABAergic interneurons belonging to the vasoactive intestinal peptide (VIP), lysosome-associated membrane glycoprotein 5 (LAMP5), paired box 6 (PAX6) and somatostatin (SST)-expressing interneurons (Supplementary File S6). Our in situ hybridization and immunohistochemistry study confirmed that NECAB2 is not expressed in glial cells. We found that NECAB2 was presented in the cell body of two morphologically different interneuron subtypes in the DMPFC. The interneurons with larger bodies were located in deep layers, while the smaller, more abundant cell types were distributed in upper layers, mainly in layer II. Based on the distributional and morphological data, the larger cell type may correspond to SST+ interneurons while the smaller one corresponds to VIP+ interneurons [86,87]. The presence of NECAB2 in somatodendritic compartments, as well as in axon terminals, in DMPFC interneurons is consistent with previous data on the hippocampus [88]. Based on these data, it is likely that at least one of the cell types expressing NECAB2 is implicated in suicidal behavior, possibly due to the regulation of mGluR5.

Functional Implications Based on Pathway Analyses
Our integrated analyses using GO functional annotation, as well as KEGG and Reactome pathways revealed a number of pathways altered in suicide victims. The downregulated DEGs were significantly enriched in the cell surface receptor signaling pathway of biological process and growth factor binding of molecular function, in PI3K-Akt signaling pathway, in TNF signaling pathway and in cytokine-cytokine receptor interaction. Alterations in inflammatory responses have been related to mental illnesses [89][90][91]. However, our study is the first to suggest that suicidal behavior is linked to reduced inflammatory ability in the DMPFC, suggesting a reduced microglial function in suicidal behavior.
Upregulated DEGs were enriched in axons and synapses as cell components, in synaptic function as biological process and in sodium channel activity as molecular function. Pathways found among upregulated DEGs were neuroactive ligand-receptor interactions and activation of Na-permeable kainate receptor pathways, which all suggest increased glutamatergic, specifically kainate, signaling. This change has been implicated in mood disorders [92,93], addictive disorders [94] and in alcohol dependence [95]. Further evidence from postmortem studies revealed higher expression levels of the majority of glutamatergic genes in the DLPFC associated with MDD and suicide [74,96], and now we provide evidence that it also holds for the DMPFC. Target genes of high interest included GRIK2, which likely plays a role in emerging suicidal thoughts after antidepressant treatment [97]. In a subsequent study, polymorphisms in GRIK2 gene were further associated with suicidal ideation in MDD patients following antidepressant treatment [98]. A rodent study showed that GRIK2-deficient mice were more impulsive and aggressive, suggesting that it has a unique role in controlling the behavioral symptoms of mania [99]. In particular, the outcome of completed suicide has been associated with increased expression of GRIK2 [73]. The present data further suggest that elevated kainate receptor signaling in the DMPFC can predispose the development of committing suicide, and hence can be a potential biomarker.
Functional annotation of down-and upregulated DEGs revealed that there are gene sets overlapping common pathways associated with depression. Searching for genes, which share common functions with depression-related pathways, we found that both down-and upregulated genes were present, which participate in neurotransmission and cell metabolism processes taking place in the DMPFC (Table 2).

Functional Cluster Analysis of Gene Expressions in the DMPFC
Co-expression analysis on the DEGs was performed to identify the key modules of highly co-expressed genes. The largest cluster in a downregulated gene network was associated with gliogenesis, while another cluster was associated with type-I interferon signaling pathway, suggesting reduced glial function and the role of interferon gamma in reduced inflammatory response ability. In turn, modules in upregulated gene networks were associated with neurotransmitter receptor internalization, sodium channel activity and peptide neuromodulator function. Alteration in the neurotransmitter receptor internalization pathway is consistent with the assumption that dysfunction of neurotransmitter receptors may contribute to the pathophysiology of MDD [93,100,101].
Intramodular hubs have an important role in the enriched pathways of the modules in co-expression networks. The major hub gene of Cluster 1 in the downregulated gene network was SORBS3 (vinexin). This protein modulates the actin cytoskeleton and negatively regulates autophagy. Its expression increases with age in mouse and human brain tissue, contributing to autophagic decline in mammalian brain aging [102,103]. Moreover, a meta-analysis provided evidence that SORB3 has a strong relationship with MDD [104]. The hub gene from upregulated gene network, calcyon-neuron-specific vesicular protein (CALY) regulates dopamine-related signaling [105] and dopamine D1 receptor internalization of pyramidal cells in the prefrontal cortex and dorsal striatum [106]. Other studies reported that CALY has trafficking functions primarily involved with neural development and synaptic plasticity [107,108]. These results suggest that intramodular hub genes could function as potential biomarkers for future therapeutic interventions.
The top 10 suicide-related down-and upregulated hub genes were screened from the co-expression networks according to the MCC scores, including the downregulated genes MSI2, SNTA1, S100A16, RHOC and SORBS3, and the upregulated genes ACTL6B, CCK, CALY, CPNE9 and SYT5. One of the top downregulated hub genes, the Musashi-2 (MSI2) RNA-binding protein, was found to be expressed in neural progenitor cells [109] and involved in the regulation of cell cycle and development [110]. Another downregulated hub gene, the alpha-1-syntrophin (SNTA1), plays an important role in the organization of the localization of a variety of membrane proteins and establishes a link between the extracellular matrix and receptors [111,112]. Between the upregulated hub genes, the actinlike protein 6B (ACTL6B) belongs to the neuron-specific chromatin remodeling complex, and thus plays a crucial role in neural development and dendritic outgrowth [113]. A defect in neuronal dendritic branching was observed in ACTL6B knockout iPSC-derived neuronal cells [114]. Another upregulated hub gene is the cholecystokinin (CCK) produced in the gastrointestinal tract. It acts as a pancreatic enzyme to regulate digestion; however, it has been also identified as a neuropeptide abundantly expressed in vertebrate brains [115]. The functional and pathway enrichment analysis found integral biological processes enriched in terms related to cell cycle regulation, cell motility and signal transduction.

Functions Supported by Known PPIs of DEGs
Analysis of the top 10 most ranked downregulated genes (PECAM1, ERBB2, ITGB1, EGFR, STAT3, NOTCH1, CD44, IL6, TLR4 and FN1) from the PPI network revealed that these genes were strongly enriched in cell-surface receptor signaling pathway, astrocyte differentiation and regulation of MAPK cascade, suggesting that these pathways might be implicated in suicide behavior. As previously defined, two of the top most ranked genes, the adhesive stress-response protein PECAM1 and the proto-oncogene ERBB2 were reported as decreasing in the DLPFC in suicide victims [54]. Analysis of top 10 most ranked upregulated genes (CCK, GABRD, CRHR1, CALY, COL11A2, COL24A1, SCN3B, GRIK1, GRIK2, CACNA1G) revealed that these genes were associated with regulation of membrane potential, nervous system processes and cell-cell signaling. A pilot study examining the GABAergic system in suicide victims found GABRD variants in the prefrontal cortex of patients associated with depressive disorder [116]. Recently, a strong association between CRHR1 polymorphism and suicide attempts was shown [117]. Löfberg et al. [118] found that MDD patients associated with the number of previous suicide attempts had higher CCK levels. A more recent study showed that CCK levels were higher among first suicide attempters after the first 12 h following the attempt [119]. Likewise, a study investigated postmortem brain tissues of completed suicides observed higher number of CCK receptors in the frontal cortex compared to healthy controls [120].
A multi-locus GWAS study demonstrated that several allelic variations in glutamatergic synaptic signaling, such as GRIK1, GRIK2, GRIK3, GRIN1, GRIN2A, GRIN2C or GRM7, were associated with MDD [77]. A subsequent systematic review revealed that several crucial candidate genes for MDD are involved in glutamate neurotransmission, [45]. Moreover, several studies found decreases in mRNA expression of NMDA and AMPA receptor subunits associated with schizophrenia in the frontal cortex [121][122][123]. The involvement of glutamate in psychiatric and medical conditions has been intensively examined; however, earlier studies mostly focused on the biology and pathophysiology of ionotropic glutamate receptors. In turn, metabotropic receptors (mGluR) can also modify neuronal activity through G-protein coupled signaling. Indeed, strong interactions have been reported between mGluR5 and NMDA receptors, suggesting that mGluR5 might be implicated in mediating neural plasticity as well as learning and memory processes [124]. Two proteins from upregulated DEGs, the CACNG8 and CALY, are considered neurotransmitter receptor regulatory proteins. CALY regulates dopamine D1 receptor internalization in the prefrontal cortex and dorsal striatum [106], and CACNG8 (also known as TARP-γ8) negatively modulates AMPA receptor signaling [125].

Genes Associated with Depression and Comorbidities
Among the top DEGs, 10 genes have been previously associated with depression, including CSF3, IL1R2, IL6, SERPINA3, MT-ND5, MT-ND6, P2RX2, ATP4, PROC and CARTPT, based on the DisGeNET database [126]. Studies of genetic polymorphisms in psychiatric disorders have shown that CARTPT gene mutation exhibits increased anxiety and depression [127]. Furthermore, the peptide encoded by CARTPT is a candidate biomarker for MDD because of its effects on mood regulation [128,129]. Additionally, we identified 12 genes from the most ranked down-and upregulated DEGs associated with depression, such as GABRD, IL6, NOTCH1, EGFR, ERBB2, TLR4, STAT3, GRIK1, GRIK2, CACNA1G, CCK and CRHR1 and six genes from the hub gene list of the co-expression analysis, including BCL2, CHDH, GABRD, NRGN, SERPINF1 and CCK. The medical history of the subjects was obtained from clinical records, interviews with family members and relatives, as well as pathological and neuropathological reports. All personal data are stored in strict ethical control, and samples were coded before the analyses of tissue.

Tissue Preparation
Postmortem human brain tissue samples from the dorsomedial prefrontal cortex (DMPFC) (Brodmann area 9) were acquired from the Human Brain Tissue Bank (Semmelweis University, Budapest, Hungary). There is considerable variability in how the topographical extension of the DMPFC is defined. The medial prefrontal cortex (MPFC), as its name indicates, occupies the medial surface of the frontal lobe over the anterior cingulate gyrus. The MPFC is divided into dorsal and ventral prefrontal cortex with a theoretical horizontal line through the most rostral (genual) point of the corpus callosum. (These two parts are frequently called anterior and posterior prefrontal cortex.) The posterior border of the DMPFC is defined by the precentral sulcus, which separates it from the premotor cortex. The border of the dorsomedial and dorsolateral prefrontal cortex in coronal sections can be defined by the deep superior frontal sulcus (Figure 8).
The white matter also serves as an excellent landmark: the superior longitudinal cortical pathway, as a well visible white matter bundles, are divided into a medial and a lateral portion in the dorsal prefrontal cortex. The cortical samples were collected from 16 subjects: 8 control subjects with no history of psychiatric or neurological diseases (2 females and 6 males, mean age of 65.4 ± 5.6) and 8 suicide victims (3 females and 5 males, mean age of 53.6 ± 4.8). The age of the subjects ranged from 31 to 89 years ( Table 3). The selected control subjects were not diagnosed with any psychiatric disorder, while the suicide victims were without evidence of acute or chronic depression. Brains were removed from the skull with a postmortem delay of 1 to 10 h, frozen rapidly on dry ice and stored at −80 • C until microdissection. Serial coronal sections were cut from the frontal lobe. Sections between 50.0 and 30.0 mm from the origin of AC-PC coordinates (see [130]) (practically at the first appearance of the genu corporis callosi or the first appearance of the anterior pole of the lateral ventricle) were punched out. Special microdissection needles with 8 and 15 mm inside diameters were used [131,132]. Tissue pellets (1-3 per brain) include samples from the superior frontal gyrus (Brodmann area 9) (both gray and white portions within the gyrus) and the paracingulate cortex (also called the dorsal part of the pre-and supragenual portions of the anterior cingulate cortex, Brodmann area 32), the upper bank and the deep portion of the gyrus (Figure 8). In samples, the superior frontal gyrus/paracingulate gyrus ratio was about 80-85%, respectively. The microdissected tissue pellets were collected in 1.5 mL Eppendorf tubes and kept until use at −80 • C. Throughout the microdissection procedure (slicing, micropunch, storage) the tissue samples were kept frozen.

RNA Sequencing Analysis and Data
We performed RNA sequencing on cortical samples to determine RNA expression changes in the dorsomedial prefrontal cortex related to suicidal behavior. Sample preparation, sequencing library construction and RNA sequencing were conducted by the Beijing Genomics Institute's (BGI, Hongkong, Shenzhen, China) standard procedure. To extract total RNA, 20 mg of postmortem brain tissue was subjected to RNA sequencing on the BGISEQ-500 platform (BGI) using the 100-bp paired-end sequencing strategy. RNA size, concentration and integrity were verified using Agilent 2100 Bioanalyzer (Agilent Technologies). All the generated raw sequencing reads were filtered using SOAPnuke (1.5.2), a filter software developed by BGI company, to remove reads containing adapters, reads in which unknown bases are more than 5%, and low-quality reads (>20% of the bases with a quality score lower than 15). After filtering, the remaining clean reads were obtained and stored in FASTQ format [133]. HISAT2 (2.0.4) software in combination with the Bowtie index (2.2.5) [134,135] was used to map clean reads to the human reference genome (hg19), respectively. The average mapping ratio with genes was 76.12%, and the average mapping ratio with reference genome was 91.28% (Supplementary Table S1). The output of HISAT2 was imported into the RStudio environment (R version 4.0.4; RStudio version 1.4.1106) using the Rsubread package (2.4.3), which converted data from the transcript level to the gene level, then aligned to the human genome (GRCh38) and counted into genes focusing on all annotated protein-coding genes. The Rsubread facilitates the RNA-seq read data analyses, producing the quality assessment of sequence reads, read alignment and read summarization, among others [136]. To identify differentially expressed genes (DEGs), the DESeq2 package (1.34.0) as the standard workflow was used for the detection of DEGs [137]. Before using DESeq2, we performed minimal prefiltering to keep the rows that have at least one read total. To increase power, stricter filtering is automatically applied via independent filtering on the mean of normalized counts within the results function. The complete dataset of DEGs is listed in Supplementary File S1. The DEGs were identified by Benjamini and Hochberg-adjusted p-value (padj < 0.05) [138] and log2 FC ratio ≥ ±1 was defined as the thresholds to discriminate significant DEGs. DEGs were visualized in a volcano plot created with the R package ggplot2 (3.3.5). Hierarchical clustering of DEGs was analyzed using the heatmap.2 function in the gplots (3.1.1) package for R. Heatmap created using two color arrays. Pathway-based and co-expression analysis helps to further understand genes' biological functions. The white matter also serves as an excellent landmark: the superior longitudinal cortical pathway, as a well visible white matter bundles, are divided into a medial and a lateral portion in the dorsal prefrontal cortex. The cortical samples were collected from 16 subjects: 8 control subjects with no history of psychiatric or neurological diseases (2 females and 6 males, mean age of 65.4 ± 5.6) and 8 suicide victims (3 females and 5 males, mean age of 53.6 ± 4.8). The age of the subjects ranged from 31 to 89 years (Table 3). The selected control subjects were not diagnosed with any psychiatric disorder, while the suicide victims were without evidence of acute or chronic depression. Brains were removed from the skull with a postmortem delay of 1 to 10 h, frozen rapidly on dry ice and stored at −80 °C until microdissection. Serial coronal sections were cut from the frontal lobe. Sections between 50.0 and 30.0 mm from the origin of AC-PC coordinates (see [130]) (practically at the first appearance of the genu corporis callosi or the first appearance of the an-

Gene Ontology and Pathway Enrichment Analysis of DEGs
To assess the function of the DEGs, functional enrichment analysis was performed using the gprofiler2 R package (version 0.2.0). The Gene Ontology (GO) analysis was queried using the gost function in "ordered" mode [139]. The differentially expressed genes were ranked by the adjusted p-value significance (p.adjust < 0.05). Functional enrichment was performed for upregulated and downregulated DEGs, separately. The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome resources were used for functional annotation and pathway analysis [140,141]. The results of the functional enrichment were visualized in Manhattan plots created with the gprofiler2. The significant GO terms were listed in Supplementary File S2. The functional enrichment analyses of the top 10 down-and upregulated DEGs, and also the top 10 down-and upregulated hub genes identified by gene co-expression networks, were carried out using Database for Annotation, Visualization and Integration Discovery (DAVID) v6.8 online server to annotate gene-related biological mechanisms using standardized gene terminology. Cell types in the brain are diverse with broadly different expressional profile, morphology and function. Thus, there is a need to define the functional assessment of DEGs to localize their expression in molecularly defined subtypes of the DMPFC. To address this need, we used the Allen Human Brain Atlas (AHBA) transcriptomic data derived from adult human brain. We used the "Multiple Cortical Area-SMART-seq (2019) dataset for cell-type specific genes which was downloaded from the Allen Institute for Brain Science website (www.brain-map.org, (accessed on 13 December 2021)).

Protein-Protein Interaction Network Construction
As an alternative approach, STRING version 11.5 online tool [142] was used to construct a protein-protein interaction (PPI) network of selected genes. Using the STRING database, genes with a minimum required interaction score of 0.4 were chosen to build a network model visualized by Cytoscape software version 3.8.2 [143]. CytoHubba [144], a Cytoscape Plugin, was used to evaluate the most ranked genes of PPI network. The top 10 most ranked genes were screened out by node degree from both down-and upregulated gene networks. Red colors denoted the higher degree of a gene. Gene list enrichment of the top 10 most ranked DEGs was identified through the stringApp. Functional annotation terms were considered significantly enriched with an FDR-corrected p-value < 0.05.

Disease-Associated Gene Sets
Human disease-associated gene sets were obtained from the DisGeNET database, which is known to integrate disease-gene links from several sources [126]. The list of top 10 down-and upregulated DEGs and hub genes were imported into the DisGeNET database browser (http://www.disgenet.org/, (accessed on 15 December 2021)), then the curated gene-disease association file was downloaded. Keywords with special emphasis on mood disorders and associated comorbidities were used for filtering categories to gather those genes involved in the pathophysiology of depression. The result was imported and listed in Supplementary File S4.

Co-Expression Network Construction and Functional Annotation
Coexpression correlation was computed using the Pearson correlation coefficient, as described previously [145]. Read counts of each gene and each sample were normalized via median normalization using the EBSeq R package [146]. The correlation (Pearson's) and correlation significance of every pair DEG (for downregulated genes adjusted p-value < 0.01, for upregulated genes adjusted p-value < 0.05) was calculated using a logarithmic matrix of read counts with the psych R package [147]. A network table containing the statistically significant correlations across the whole data set for every pair of DEGs was generated, and to calculate network statistics, igraph R package [148] was used. The network statistics (degree and betweenness centrality for each node) were uploaded to Cytoscape. The gene symbols were designated as the identifiers of nodes. The correlation, degree and betweenness were mapped to the edge color, edge width and node size, respectively. Regarding the topological analysis, GLay community clustering evaluation [149] was performed in order to identify densely connected nodes. The functional annotation analysis of the network was performed using the ClueGO application [150]. Functional annotation terms were considered significantly enriched with a Bonferroni-corrected p-value < 0.05. Further details on the interplay of the reconstructed networks were examined using the stringApp. For screening hub genes from co-expression networks, a Maximal Clique Centrality (MCC) algorithm was used to select the hub genes [151]. The network datasets and functional annotation results of the networks were shown in Supplementary File S5.

Validation of Expression Changes by qRT-PCR
To confirm the RNA-seq data, a subset of differentially expressed genes whose selection was based on the known or potentially probable involvement in depression and neuron-specific functions were validated using quantitative real-time PCR (qRT-PCR). The procedure was carried out as described previously [152]. Briefly, total RNA was isolated from approximately 20 mg of frozen postmortem brain tissue-from the same samples used in RNA-seq-using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) as lysis buffer combined with RNeasy Mini kit (Qiagen, Germany) following the manufacturer's instructions. The quality and quantity of extracted RNA were determined using NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and only those with A260/A280 ratio between 1.8 and 2.1 were used in subsequent experiments. The concentration of RNA was adjusted to 500 ng/µL, and it was treated with Amplification Grade DNase I (Invitrogen, Carlsbad, CA, USA). The isolated RNA concentration was calculated and normalized with RNase-free water and reverse transcribed into cDNA using SuperScript II Reverse Transcriptase Kit (Invitrogen, Carlsbad, CA, USA). After 10-fold dilution, 2.5 µL of the resulting cDNA was used as a template in PCR performed in duplicate using SYBR Green dye (Sigma, St Louis, MO, USA). The PCR reactions were performed on CFX-96 C1000 Touch Real-Time System (Bio-Rad Laboratories, Hercules, CA, USA) with iTaq DNA polymerase (Bio-Rad Laboratories, Hercules, CA, USA) in total volumes of 12.5 µL under the following conditions: 95 • C for 3 min, followed by 35 cycles of 95 • C for 0.5 min, 60 • C for 0.5 min and 72 • C for 1 min. A melting curve was performed at the end of amplification cycles to verify the specificity of the PCR products. The primers used for qRT-PCR were synthesized by Integrated DNA Technologies, Inc., (IDT, Coralville, IA, USA) and used at 300 nM final concentration. Sequences of primers are listed in Supplementary Table S2. Housekeeping genes ACTB, GAPDH and LDHA were used as internal controls and the relative gene expression values were calculated from their averages using the 2 − Ct method.

Preparation of In Situ Hybridization Probes
A mixture of cDNAs prepared with RT-PCR from the human DMPFC was used as a template in PCR reactions with primers for NECAB2 (CAGGATCTTGGTGCCAGCT and TGTGGTCAGTGTGGGTCATG) yielding a probe that corresponds to the GenBank accession number NM_019065.2. The purified PCR products were applied as templates in a PCR reaction with the primer pairs specific for the probe and also contained the T7 RNA polymerase recognition site added to the reverse primers. Finally, the identities of the cDNA probes were verified by sequencing them with T7 primers.

In Situ Hybridization Histochemistry
Two freshly frozen DMPFC brain blocks of subjects were used: one from a 75-year-old female control subject with negative clinical reports for any major diseases and one from a 72-year-old female suicidal individual. Using a cryostat, serial coronal sections (12 µm) were cut and mounted on positively charged slides (SuperfrostPlus ® , Fisher Scientific), dried and stored at −80 • C until use. Further steps were performed according to the procedure described previously by Dobolyi et al. [153]. Briefly, antisense [35S]UTP-labeled riboprobes were generated from the above-described DNA probes using T7 RNA polymerase of the MAXIscript Transcription Kit (Ambion, Austin, TX, USA) and used for hybridization at 1 million DPM (discharges per minute) activity per slide. Washing procedures included a 30 min incubation in RNase A, followed by decreasing concentrations of sodium-citrate buffer (pH 7.4) at room temperature and subsequently at 65 • C. Following hybridization and washes, slides were dipped in NTB nuclear track emulsion (Eastman Kodak) and stored at 4 • C for 3 weeks for autoradiography. Then, the slides were developed and fixed with Kodak Dektol developer and Kodak fixer, respectively, counterstained with Giemsa and coverslipped with Cytoseal 60 (Stephens Scientific, Kalamazoo, MI, USA).

Tissue Collection for Immunolabeling
Immunohistochemistry was used to assess the distribution of NECAB2 protein in the DMPFC. For immunolabeling, one DMPFC brain block from a 62-year-old female control individual was cut into 10 mm thick coronal slice and immersion fixed in 4% paraformaldehyde in 0.1 M phosphate-buffered saline (PBS) for 5 days. Subsequently, the block was transferred to PBS containing 0.1% sodium azide for 2 days to remove excess paraformaldehyde. Then, the block was placed in PBS containing 20% sucrose for 2 days of cryoprotection. The block was frozen and cut into 60 µm thick serial coronal sections on a sliding microtome. Sections were collected in PBS containing 0.1% sodium azide and stored at 4 • C until further processing.

Microscopy and Photography
Sections were examined using an Olympus BX60 light microscope also equipped with fluorescent epi-illumination and a dark-field condenser. Images were captured at 2048 × 2048 pixel resolution with a SPOT Xplorer digital CCD camera (Diagnostic Instruments, Sterling Heights, MI, USA) using a 4× objective for dark-field images, and 4-40X objectives for bright-field and fluorescent images. Confocal images were acquired with a Zeiss LSM 70 Confocal Microscope System using a 40-63X objectives at an optical thickness of 1 µm to count varicosities and 3 µm to count labeled cell bodies. Contrast and sharpness of the images were adjusted using the 'levels' and 'sharpness' commands in Adobe Photoshop CS 8.0. Full resolution was maintained until the photomicrographs were cropped, at which point the images were adjusted to a resolution of at least 300 dpi.

Statistical Analysis
Demographics were contrasted using the Chi-Square test or Welch's unequal variances t-test (p < 0.05). The plot and heatmap were generated using R version 4.0.4 (https://www.r-project.org/, (accessed on 10 January 2022)). Quantitative statistics of qPCR results were performed with the GraphPad Prism version 8.0.1 (GraphPad Software, San Diego, CA, USA). Normality distribution was assessed via Shapiro-Wilk test. For comparisons between values of the two groups, Welch's unequal variances t-test was applied. The nonparametric Mann-Whitney U test was used for data that were not normally distributed. Differences were considered statistically significant when p < 0.05 and plotted as mean ± S.E.M (standard error of the mean).

Conclusions
We performed a comprehensive transcriptome study of DMPFC in suicide victims, using RNA-seq and observed significant associations between specific transcriptome changes in DMPFC and suicidal behavior. They include (but are not limited to) the suppression of genes that regulate the inactivation of immune responses and glial cell differentiation, and the upregulation of genes involved in glutamatergic synaptic transmission. The identified gene expression alterations may reveal dysregulated pathways and functional cascades involved in the pathophysiology of suicidal behavior. Therefore, we conclude that the DEGs and pathways identified in this study could be linked to depression, albeit a common secondary pathology arising from dysregulated network activity is also possible. These results are supported by the large overlap of DEGs between previously reported genes altered in depressive disorders. Furthermore, due to the alteration in pathways relating to gliogenesis, suicidal behavior can affect glial cell numbers and cellular morphology in the medial part of the prefrontal cortex. Additionally, alterations in glutamate neurotransmission including both ionotropic and metabotropic receptors, as well as neuroinflammatory processes potentially contributing to neuroplastic changes, can contribute to the dysfunction of the DMPFC. The DMPFC is one of the nodes in connections with different cortical networks receiving direct viscero-and somatosensory inputs, motivational/emotional inputs from the orbitofrontal cortex and cognitive inputs from the dorsolateral and ventral prefrontal cortex. By processing all this information together using activated working memories, the DMPFC may work out motor programs and projects (through the default mode network) to the precuneus and directly to the premotor and supplementary motor cortical areas. The precuneus is able to interfere or even block different actions of the motor cortex, and through self-reflection may determine suicidal behavior.
Taken together, our findings may contribute to the pathophysiology of depressive or suicidal behavior. The high percentage of altered genes in suicide victims involved in depression is in line with the fact that depression is the leading cause of suicide. We conclude that our study provides new insights into the genomic underpinnings of major depressive disorder and may allow for novel precision therapeutic development strategies targeting depression to be implemented in future studies.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available in the supplementary material of this article. The raw sequencing data files are available in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) under the Bioproject ID: PRJNA828151, with the following sample accession numbers: SAMN27632065, SAMN27632062, SAMN27632064, SAMN27632063, SAMN27632060, SAMN27632061, SAMN27632059, SAMN27632058, SAMN27632071, SAMN27632070, SAMN27632068, SAMN27632069, SAMN27632066, SAMN27632057, SAMN27632056 and SAMN27632067.