NR3C1 and NR3C2 Genes Increase the Risk of Suicide Attempt in Psychiatric Disorder Patients with History of Childhood Trauma

Introduction

Suicidal behavior is a major public health issue. An estimated 800,000 deaths annually are attributed to suicide; these show a high prevalence of psychiatric disorders, supported by psychological autopsy studies.^1–3 Family, twin, and adoption studies show genetic and environmental components to suicidal behavior, estimating a heritability of 45%, which suggests that genetic factors could be interacting with environmental factors in the risk factors for suicide attempt (SA).⁴ In particular, childhood trauma (CT) has been found to affect the stress response in individuals with genetic vulnerability to suicidal behavior.⁵

The hypothalamic-pituitary-adrenal (HPA) axis is considered the main regulator of the stress system stimulating forward and feedback inhibition loops involving the brain, pituitary, and adrenal glands, which regulates cortisol production.⁶ During acute stress, cortisol levels increase and pulsatility remains intact, which is beneficial for promoting the survival of the fittest as part of the fight-or-flight response.⁷ However, chronic exposure to stress can result in the reversal of beneficial effects, generating long-term high cortisol levels, which can become a maladaptive stress response, which has been associated with major depressive disorder, bipolar disorder patients with a late age-of-onset, and drug-naïve first-episode patients with psychosis.^7–9 Cortisol is released from the adrenal glands and exerts its effects in the hippocampus by binding with the glucocorticoid receptor (GR) and the mineralocorticoid receptor (MR). The GR has a relatively low affinity to cortisol but is gradually occupied by higher levels of the steroid hormone under acute stress conditions. Intracellular MR, however, has ten times its affinity to cortisol, leading to a higher occupancy, even under baseline conditions. These two hippocampal steroid receptors exert opposite actions: the MR mediates the tonic influences of basal cortisol levels, and the GR facilitates negative feedback actions in response to higher cortisol levels to maintain steroid homeostasis.¹⁰ The two genes are thus essential for activating and modulating the HPA axis.

Glucocorticoid (NR3C1) and mineralocorticoid (NR3C2) receptor genes codify to proteins containing a variable amino-terminal transactivation domain, a small well-conserved DNA-binding domain that comprises two repeats of a zinc finger protein motif, and a relatively conserved carboxy-terminal domain responsible for hormone binding.¹¹ NR3C1 is located at 5q31-32 and the NR3C2 gene at 4q31.1–31.2. NR3C1 is expressed in two main isoforms due to alternative splicing of exon 9, generating the GRα and GRβ proteins. GRα diverges from GRβ only at the carboxyl-terminal end, with 15 nonhomologous amino acids. This difference impairs GRβ from binding to steroids such as cortisol. The GRα isoform is thus considered a biologically active receptor isoform.¹²

The NR3C1 gene contains a functional polymorphism, rs6198 (A → G), which alters the rates of GRα/GRβ expression, where the G allele favors the expression of the GRβ isoform. G allele carriers have also shown significantly higher cortisol responses following the psychological stressor and dexamethasone suppression test.¹¹ Another genetic variant is the rs6191, located in exon 9β in the 3’UTR region, in linkage disequilibrium with rs6198, a block related to mRNA stability and increased protein expression of the receptor.¹³ Moreover, rs33388 is in the region of intron 2, which has been suggested as a target sequence involved in regulating the NR3C1 gene or as a binding site for specific miRNAs regulating the network of genes involved in the regulation of the HPA axis.¹³

The expression and activity of the MR have also been associated with functional genetic variants of NR3C2. One of them is rs5522, a G → A substitution in exon 2, resulting in an amino acid change (Ile180Val), where a carrier of the Val allele was associated with a loss of cortisol function when used as a ligand, and with increased saliva and plasmatic cortisol responses to psychosocial stressors and the dexamethasone suppression test.^11,14 The single-nucleotide polymorphism (SNP) rs2070951 (C/G) of the NR3C1 gene has also been associated with modulation of MR activity through the altering of expression, transactivation, and activity. In vitro studies have demonstrated that the C allele is associated with increased expression of MR.¹⁰

Despite to the limited literature, there are some studies about NR3C1, NR3C2, CT, and SA. It was observed an interaction between NR3C1 and physical neglect increasing the risk of dependence of crack/cocaine.¹⁵ Another study shows that NR3C1 and emotional trauma (abuse or neglect) in childhood interact inducing reductions in left hippocampal volumes suggesting vulnerability to mood disorders in adolescent girls.¹⁶ It is important to mention that substance use disorder and mood disorders are highly implicated with SA.¹⁷ For other hand, it was observed significant association between NR3C1 and SA in patients with schizophrenia.¹⁸ Two studies failed to find an association between suicide behavior and the NR3C2 gene.^19,20 To date, few studies analyzing NR3C1, NR3C2, and CT have been conducted on patients with SA. A gene x environment (GxE) study found no interaction between 19 gene variants involved in the HPA axis, including NR3C1, and early childhood abuse in patients with SA and bipolar disorder.²¹ A previous study failed to identify interaction of NR3C1, NR3C2, and CT in the development of SA in bipolar disorder patients.²⁰ The only gene expression study has been of NR3C1, showing an under-expression (UE) of GRα in patients with suicidal behavior and childhood abuse.²² To date, the expression of NR3C2 has not been evaluated in SA. Findings concerning the interaction of NR3C1 and NR3C2 with CT in the development of SA have not been conclusive.

Aim

The aim of this study was to analyze GxE interaction of NR3C1, NR3C2, and CT, and the gene expression of NR3C1 and NR3C2 in the development of SA in patients with CT.

Hypothesis

Genetic variants of the NR3C1 and NR3C2 genes are associated with differential expression levels, which interact with a history of childhood trauma, contributing to an increase of at least one SA.

Materials and Methods

Study Design

The study design was comparative, observational, and transversal.

Sample

We included a sample of 516 Mexican patients (404 female and 112 male), who met DSM-5 criteria for bipolar disorder (n = 231), major depressive disorder (n = 221), borderline personality disorder (n = 34), and other psychiatric disorders (n = 30). They were recruited from the affective disorder’s clinic, the outpatient service, the emergency department, and the inpatient population of the Instituto Nacional de Psiquiatría Ramón de la Fuente Muñiz from October 2016 to November 2021. The sample included 274 patients (53%) who had made at least one suicide attempt and 242 (47%) who had not attempted suicide. SA was defined as self-destructive behavior with the intent to end one’s life, regardless of the resulting damage. The information about CT was systematically obtained by psychiatrists in the clinical interview and registered in the patients’ medical records. The CT was evaluated by a psychiatrist who reviewed the clinical information and classified the different subtypes of CT according to the criteria of Bernstein and Fink (1998) for neglect (emotional and physical) and abuse (emotional, physical, and sexual).²³

In addition, we analyzed the gene expression in a subsample of 96 patients (86 female and 10 male) that included 63 patients with major depressive disorder (65.6%), 24 with borderline personality disorder (25%), and nine with other psychiatric disorders (9.4%). The assessment of CT for these patients was assessed with the Childhood Trauma Questionnaire (CTQ) short form, which evaluates five types of traumas: sexual, physical, and emotional abuse, and physical and emotional neglect.²³

Inclusion criteria were psychiatric patients diagnosed with DSM-5 with and without SA; female and male patients over 18 years of age; Mexican Mestizos with a family background of three generations born in México; and capacity to provide written informed consent. Additionally, the patients implicated in the gene expression analyzes did not presented medical conditions that generate an imbalance in the HPA axis. Whereas exclusion criteria were psychiatric patients not competent to consent to the study, and insufficient sample to genetic analysis.

Ethical Approval

This study was conducted in agreement with the Declaration of Helsinki standards. The study protocol was approved by the Ethics Committee of the Instituto Nacional de Psiquiatría Ramón de la Fuente Muñiz (Approval No. CEI/C/015/2019). All participants in the study provided written informed consent.

Genotyping

DNA was extracted from peripherical leukocytes using the Flexigene DNA kit (Qiagen; Minneapolis, MN). NR3C1 (rs6198, rs6191, and rs33388) and NR3C2 (rs5522 and rs2070951) polymorphisms were selected, given that they had been previously studied in people with SA and CT.^20,21 Genotyping was analyzed using TaqMan assays for rs6198 (C_8951023_10), rs6191 (C_3234245_20), rs33388 (C_1046426_10), rs5522 (C_12007869_20), and rs2070951 (C_1594392_10). Genotypes were identified with the TaqMan allele-specific assay method using the ABI Prism^®7500 sequence detection system according to the manufacturer´s protocols (Applied Biosystems Inc.; Foster City, CA). The final volume of the reaction was 7μL, consisting of 100ng of genomic DNA, 1X TaqMan Universal Master Mix, and 0.71X SNP Genotyping Assay Mix (Applied Biosystems Inc.). After denaturing at 95°C for 10 min, 40 cycles of PCR were performed under the following conditions: denaturing at 95°C for 15s and annealing at 60°C for 1 min. The genotyping was performed without knowing the clinical status of the subjects.

Gene Expression

Total RNA was isolated from the blood samples collected within the first 24 hours after hospitalization (7:00–9:00 am), using the reagent TRIzol (Roche; Germany). The first-strand cDNA was synthesized with SuperScript™ IV VILO™ Master Mix (Thermofisher Scientific) using 500 ng of total RNA. Gene expression was analyzed with custom TaqMan assays for NR3C1 (Hs00353740_m1) and NR3C2 (Hs01031809_m1), with Glyceraldehyde-3-phosphate dehydrogenase (GAPDH; Hs02786624_g1) as a reference gene transcript, using the ABI Prism^®7500 Sequence Detection System according to the manufacturer´s protocols (Applied Biosystems Inc.; Foster City, CA). Real-time PCR was performed using a final reaction volume of 10μL, consisting of 62.5 ng of resulting cDNA, 1X TaqMan Universal Master Mix, and 1X Gene Expression Assay Mix (Applied Biosystems Inc.). cDNA was amplified for 40 cycles at 95°C for 15s and 60°C for 1 min. All samples were performed in duplicate and information about the study sample blinded.

Statistical Analyses

Demographic and clinical characteristics were analyzed using chi-square and Student’s t-tests. The Shapiro–Wilk test was used to analyze the normality of the total and subscale scores of the CTQ short form. Total CTQ, emotional neglect, and emotional abuse scores were divided into low and high trauma groups defined by the upper and lower 27% of the total scores, due to their normal distribution, based on previous grouping methods.^24,25 Physical neglect, physical abuse, and sexual abuse were not normally distributed; a median-split approach was therefore used to dichotomize these scores between high and low trauma severity.²⁶

The statistical analyses were performed using the program RStudio version 4.2.0.²⁷ The power analysis was performed with the R package “gap” version 1.2.3–6,²⁸ showing a sample power of 0.99, assuming an additive genetic model, a risk allele frequency of 0.09, a population prevalence of SA of 2.7%, an α level of 0.05, and a portion of cases of 0.53 in a sample of 516 patients. Analysis of genotype and allele distributions in patients with and without SA were performed using chi-square tests with Epidat version 3.1.²⁹ Bonferroni correction for multiple testing was applied (for five polymorphisms corrected at p < 0.01). A linkage disequilibrium (LD) matrix based on the patients’ D’ parameter was estimated using Haploview version 4.2.³⁰ THESIAS software was used to analyze haplotype effects. The results were expressed as a haplotypic odds ratio in reference to the most frequent haplotype.³¹

We also analyzed the GxE interaction using Multifactor Dimensionality Reduction (MDR) software version 3.0.2 and MDR permutation testing software version 1.0 beta.^32,33 The MDR method is used for detecting high-order nonlinear or non-additive interactions in case control studies with small sample size to improve the effect of multiple genetic loci and environment in the development of a disease.^33,34 MDR collapses multi-locus and environment data into a single-dimensional variable with two levels: high and low risk.³² The interaction models are evaluated using testing balanced accuracy (TBA), a measure of how often individuals are correctly classified with respect to their case/control status, and cross-validation consistency (CVC), which evaluates the consistency with which individuals are classified. The best final model is defined as one with values of TBA between 0.55 and 0.69 and a CVC of 10, which is considered strong evidence of a multifactor interaction, and a significant p-value derived from 1000 permutations.^34,35

The gene expression analyses were performed using GraphPad Prism software version 5.0.0.³⁶ All of the results were normalized to GAPDH. Cycle threshold values were averaged and relative NR3C1 and NR3C2 expression levels were determined using the comparative cycle threshold method.³⁷ Differences between NR3C1 and NR3C2 expression of patients with and without SA were calculated using a t-test. Changes in gene expression in relationship to low or high total CT score and the subscale scores were analyzed using ANOVA and Tukey’s multiple comparison post-hoc test or the Kruskal–Wallis test and Dunn’s multiple comparison post-hoc test.

Results

Analysis of Demographic and Clinical Characteristics

Analysis of demographic characteristics revealed significant age and sex differences between patients with and without SA (Table 1). There were significant differences in primary diagnoses, with a high level of major depressive disorder and borderline personality disorder in patients with SA (Table 1). The SA patients showed greater comorbidity with borderline personality disorder, eating disorder, and post-traumatic stress disorder than those without SA (Table 1). There were differences in all CT subtypes, except for physical neglect, between patients with and without SA.

Table 1 Demographic and Clinical Characteristics of the Sample

Association Analysis

Genotype and allele distributions of NR3C1 (rs6198, rs6191, and rs33388) and NR3C2 (rs5522 and rs2070951) gene polymorphisms are shown in Table 2. There were no differences in NR3C1 or NR3C2 gene polymorphisms between patients with and without SA (Table 2).

Table 2 Genotype and Allele Frequencies of NR3C1 and NR3C2 SNPs in Patients with and without Suicide Attempt

Haplotype Analysis

The NR3C1 and NR3C2 haplotype structures are shown in Figure 1. NR3C1 showed a block composed of rs6198, rs6191, and rs33388 (D’ = 0.83, r² = 0.14), and NR3C2 showed a block composed of rs5522 and rs2070951 (D’ = 0.91, r² = 0.35). The haplotype analysis did not show any statistical differences between subjects with and without SA (data not shown).

Figure 1 NR3C1 and NR3C2 linkage disequilibrium structure. The numbers in the squares refer to pairwise D´. Haploview analysis showed a block in NR3C1 and a block in NR3C2 using a setting of average pairwise D´ within-block of ≤ 0.80, while the color intensity correlates to the pairwise r2 values. The direction of gene transcription is from left to right.

Gene Expression

There were no significant differences in NR3C1 and NR3C2 gene expression levels between patients with and without SA (t = 0.7, p = 0.48; t = 1.4, p = 0.15, respectively). Analysis of the gene expression in relationship to CT and its subtypes showed overexpression (OE) of NR3C1 associated with high scores for physical abuse (p < 0.0001), sexual abuse (p = 0.0006), and emotional neglect (p = 0.014) in patients with SA. It also detected NR3C1 OE in subjects with high scores for total CTQ (p = 0.0012) and physical neglect (p < 0.0001) compared with subjects with low scores. No gene expression changes were detected between the groups after the post-hoc test (p > 0.05) for emotional abuse. The analysis of NR3C2 gene expression and CT showed UE in psychiatric patients with high scores for total CTQ (p < 0.0001), emotional abuse (p < 0.0001), physical abuse (p < 0.0001), sexual abuse (p < 0.0001), and emotional neglect (p < 0.0001) compared with patients that presented low scores.

Analysis of Association Between SNPs and Gene Expression

Patients were subgrouped into carriers and non-carriers of alleles related to previous reports of altered function or differential levels of gene expression of NR3C1 and NR3C2;^10,11,13 the NR3C1/rs6198 were subgrouped into AA and AG+GG carriers, the rs6191 into CC and CA+AA carriers, and rs33388 into AA y AT/TT, the NR3C2/rs5522 into AA and AG+GG carriers, and the rs2070951 into CC and GG+CG carriers. An OE of NR3C1 was found in patients with SA who were A allele carriers of rs6191 (p = 0.0015); also, an OE of NR3C1 in carriers of the G allele of rs6198 (p < 0.0001) and an UE of NR3C2 in carriers of the G allele of rs5522 (p < 0.0001) in the sample as a whole. No differences in gene expression were found between those with and without SA in NR3C1/rs33388 (p = 0.43) or NR3C2/rs2070951 (p = 0.2).

Gene x Environment Interaction

Finally, we analyzed the GxE interaction for the risk of developing a SA using the different subgroups of carriers or non-carriers of alleles related to differential levels of function. In the analysis of NR3C1, NR3C2, and CT, we found a prediction of SA risk (TBA = 0.6027) with an OR of 2.76 (1.93–3.94). In addition, the analysis of CT subtypes showed three best interaction models in relation to neglect (TBA = 0.5621; OR 2.12; 95% CI [1.45–3.11]), emotional abuse (TBA = 0.5616; OR 2.13; 95% CI [1.48–3.05]), and sexual abuse (TBA = 0.5673; OR 2.38; 95% CI [1.56–3.64]) for the risk of developing SA. All GxE interaction models showed a CVC of 10 out of 10 and p < 0.0001 (Table 3).

Table 3 Best Models for the Risk of Suicide Attempt Between NR3C1, NR3C2 and Childhood Trauma

Discussion

This study investigated the effect of two candidate genes involved in the activation and modulation of the HPA axis and CT in the development of SA in a sample of Mexican patients with psychiatric disorders. Consistent with previous studies, we observed a high proportion of women, primary diagnosis of major depressive disorder, and comorbidity with borderline personality disorder.^38–40 Our CT assessment was carried out by psychiatrists using clinical interviews and medical records with the operative definition of neglect and abuse proposed by Bernstein and Fink.²³ CT frequencies were found to be similar to those reported in previous studies using the CTQ in SA patients.^41–43 This study did not used a healthy control group and to our knowledge, there are not Mexican studies that analyzed the CT in a healthy group. A study in healthy Colombian population obtained the information of the subtypes of CT, showing a similar frequency of sexual abuse in our patients without SA (21.74–19%); also, observed similar frequencies of emotional (36.84–42.3%) and physical abuse (29.63–29.9%) and emotional neglect (28.57–31.7%) with our patients with SA, and higher frequency of physical neglect (30.56%) compared with the groups with (14.2%) and without SA (10.7%).⁴⁴ Could be interesting to analyze the frequency of the CT in a healthy Mexican population to know more directly their impact in the development of SA in psychiatric patients.

We did not find an association between NR3C1 and NR3C2 gene polymorphisms and SA. Yin et al⁴⁵ likewise found no association of an NR3C1 SNP in patients with SA.

SA is a multifactorial behavior caused by the interaction of various genes and environmental factors. Our study showed that the patients with SA with a history of emotional neglect or physical or sexual abuse presented OE of the NR3C1 gene. A previous study observed an UE of NR3C1 in a sample of patients with suicidal behavior and a background of childhood abuse, as compared with a control group.²² The difference in the results may be explained by the inclusion in their sample of patients with SA and suicidal ideation, and the measurement of gene expression measure of 9β isoform of NR3C1. In contrast, our study included only patients with SA and measured the total gene expression of NR3C1. We also analyzed the gene expression of NR3C2, observing that both genes presented changes in the level of gene expression related to high scores for neglect and abuse. Our findings suggest that a history of CT in psychiatric patients affects the gene expression of both steroid hormone receptors. Several studies in psychiatric disorder patients with CT reported changes in the methylation patterns altering the expression levels of NR3C1.⁴⁶ Regarding NR3C2, a study showed an under-expression of NR3C2 in hippocampus of suicide victims compared with non- suicide subjects.⁴⁷ To our knowledge, the effects of CT in the NR3C2 expression levels have not been investigated in SA. Therefore, it is important to replicate our findings in an independent sample.

We also observed GxE interaction of NR3C1, NR3C2, and CT in psychiatric patients with SA. In contrast, two previous studies did not show interaction. The first study showed no interaction of NR3C1, NR3C2, and CT in Caucasian bipolar disorder patients with SA.²⁰ The second analyzed 235 SNPs of genes involved in the HPA axis and reported no interaction between NR3C1, childhood physical or sexual abuse, and SA in patients of European-American ancestry with bipolar disorder.²¹

These SNPs may alter the binding of transcription factors, creating a new protein-binding motif or causing alterations in the binding affinity of transcription factors. Alteration of the nucleotide sequence in a protein-binding motif may alter the binding of DNA transcription factors, affecting their regulatory function.⁴⁸ To find out how rs6198 and rs6191 of NR3C1 and rs5522 of NR3C2 could affect the regulation of gene expression by altering the dynamics of transcription factors, we searched in JASPAR⁴⁹ for DNA-binding motifs. We found that rs6198/NR3C1 is located near predictive cis-regulatory element (CRE) regions, which could be altering the expression of the neighboring genes. CREs are DNase hypersensitive sites that are commonly shared with histone modifications such as H4K4me3 and H3K27ac, or CTCF sites, where their function is to regulate the transcription of target genes. It is important to note that this SNP is in a region of DNA binding motifs, such as FOXD3, HOXB, and SOX2, that play an important role in the development of the nervous system.⁵⁰ In addition, we observed that rs6198 correlates with the DNA binding motif of the factor associated with RNApolII ZNF354A, with a possible role in the regulation of the transcription of its target genes.⁵¹ ZNF354A belongs to the family of zinc finger proteins. The binding motifs must remain intact for correct binding to DNA, and deregulation in the anchoring of these proteins could impact the transcription of target genes, particularly ZNF354A, because of its association with RNAPolII.⁵¹ The rs6191 SNP of NR3C1 is in the consensus sequence of binding of the KLF9 transcription factor regulating the transcription of genes; we hypothesize that alteration of a base pair can affect the binding of the transcription factor KLF9 to DNA, increasing transcription levels.⁵² Finally, rs5522 of NR3C2 belongs to the consensus sequence described for the binding of the transcription factor HOXD11, which plays a key role in cell differentiation, the maintenance of cell identity, and the DNA-binding motif of the TFAP4 transcription factor, whose function is the positive regulation of its target genes.^53,54 We suggest that a change in its binding motif can be associated with inhibition in the transcription of NR3C2 through modifications in the binding strength of the protein to DNA.

Based on all of these observations, our hypothesis is that the genetic variants in NR3C1 and NR3C2 are implicated in changes in gene expression due to anomalies in the binding of transcription factors to DNA and in the structure of the binding motif. These generate changes in the density of GR and MR in the hippocampus, which in turn produce deficient activation and modulation of the HPA axis in response to stress, increasing the susceptibility to SA in psychiatric patients with a history of CT.

We acknowledge some strengths of this study reside in that the operational definition of CT was defining by consensus of the group of psychiatrists based on the definition of Bernstein and Fink.²³ In addition, the use of MDR program for the GxE interaction minimized the false positive results due to multiple testing, avoid problems developed with the use of parametric statistics in the analysis of high-order interactions, and it is not assumed a genetic model.³⁴

There are several limitations to this study. First, we did not include the CTQ to assess CT in the entire sample, so we cannot perform the GxE interaction dichotomizing between high and low trauma. However, previous studies consider it valid to obtain the information from retrospective reports.^55,56 Second, we analyzed the total expression of the NR3C1 and NR3C2 genes. However, there are many isoforms for GR and MR, and it would have been useful to analyze the gene expression in these different isoforms. Third, we did not include the assessment of personality traits involved in the development of SA, such as impulsivity, aggressiveness, and anxiety.

Conclusion

In summary, our study supports the evidence of interaction between NR3C1, NR3C2, and a history of CT in the development of a SA. Our findings suggest that genetic variants of these two steroid hormone receptors are related to differential expression levels, increasing the risk of a SA in psychiatric patients with a history of neglect, emotional abuse, and sexual abuse. However, these findings should be further analyzed in a larger sample of psychiatric patients with SA, taking into consideration clinical evaluations related to personality traits and the expression profiles of genes involved in the stress regulator system to elucidate the role of the HPA axis and CT in the etiology of SA.