Susanne Edelmann1
Ariane Wiegand1,2,3
Thomas Hentrich4,5
Sarah Pasche1
Julia Maria Schulze-Hentrich4,5
Matthias H. J. Munk1
Andreas J. Fallgatter1
Benjamin Kreifelts1
Vanessa Nieratschker1*- 1Department of Psychiatry and Psychotherapy, University Hospital of Tuebingen, Tuebingen Center for Mental Health (TüCMH), Eberhard Karls University of Tuebingen, Tuebingen, Germany
- 2Max Planck Fellow Group Precision Psychiatry, Max Planck Institute of Psychiatry, Munich, Germany
- 3Department of Psychiatry and Psychotherapy, Ludwig-Maximilians-University Munich, Munich, Germany
- 4Institute for Medical Genetics and Applied Genomics, Eberhard Karls University of Tuebingen, Tuebingen, Germany
- 5Department of Genetics and Epigenetics, Faculty NT, Saarland University, Saarbrücken, Germany
Social anxiety disorder (SAD) is a psychiatric disorder characterized by severe fear in social situations and avoidance of these. Multiple genetic as well as environmental factors contribute to the etiopathology of SAD. One of the main risk factors for SAD is stress, especially during early periods of life (early life adversity; ELA). ELA leads to structural and regulatory alterations contributing to disease vulnerability. This includes the dysregulation of the immune response. However, the molecular link between ELA and the risk for SAD in adulthood remains largely unclear. Evidence is emerging that long-lasting changes of gene expression patterns play an important role in the biological mechanisms linking ELA and SAD. Therefore, we conducted a transcriptome study of SAD and ELA performing RNA sequencing in peripheral blood samples. Analyzing differential gene expression between individuals suffering from SAD with high or low levels of ELA and healthy individuals with high or low levels of ELA, 13 significantly differentially expressed genes (DEGs) were identified with respect to SAD while no significant differences in expression were identified with respect to ELA. The most significantly expressed gene was MAPK3 (p = 0.003) being upregulated in the SAD group compared to control individuals. In contrary, weighted gene co-expression network analysis (WGCNA) identified only modules significantly associated with ELA (p ≤ 0.05), not with SAD. Furthermore, analyzing interaction networks of the genes from the ELA-associated modules and the SAD-related MAPK3 revealed complex interactions of those genes. Gene functional enrichment analyses indicate a role of signal transduction pathways as well as inflammatory responses supporting an involvement of the immune system in the association of ELA and SAD. In conclusion, we did not identify a direct molecular link between ELA and adult SAD by transcriptional changes. However, our data indicate an indirect association of ELA and SAD mediated by the interaction of genes involved in immune-related signal transduction.
1. Introduction
Anxiety disorders are common and highly comorbid with other psychiatric disorders (1). A distinctive form is Social Anxiety Disorder (SAD) with an estimated worldwide lifetime prevalence of 4% (2). SAD is described by severe fear and avoidance behavior in social situations, such as fear of being the center of attention or fear of negative social evaluation, which can have a detrimental impact on daily life (3). The etiology of SAD is influenced by genetic (4, 5) as well as environmental factors. One of the most relevant environmental influences is stress. Stressful experiences in critical periods of life, especially during childhood and adolescence, such as physical or emotional neglect, sexual or physical abuse, or exposure to other types of violence as well as limited family resources (6), can lead to structural and regulatory alterations—such as disturbed programming of the hypothalamic–pituitary–adrenal (HPA) axis—contributing to disease vulnerability (7–9). Furthermore, dysregulation of the (inflammatory) immune response through childhood stress exposure can affect brain development and neurophysiology, cognition, stress reactivity and resilience and, hence, the risk for psychopathology later in life (10–14). Early life adversity (ELA) therefore represents one of the main environmental factors contributing to an increased risk for SAD (15, 16). However, the molecular link between an early stressor, such as adverse events during childhood, and the risk for SAD in adulthood remains unclear.
Changes of gene expression patterns following ELA have been identified in different organisms (17–19). In humans, monocytes of individuals exposed to early childhood maltreatment showed altered HPA axis responses to stress, evidenced by lower blood adreno-corticotropic hormone and cortisol levels. Moreover, the analysis of transcriptome-wide gene expression patterns in the same samples showed that stress-responsive transcripts were enriched for genes involved in cytokine- and inflammation-related pathways (20). In addition, co-expression network analysis identified an association of ELA with inflammation-related pathways (21). Furthermore, RNA sequencing (RNA-seq) in brain tissue revealed enrichment of differentially expressed genes in immune and GTPase function in individuals with a history of ELA as compared to control individuals without the experience of ELA (22).
Aberrant gene expression patterns of various genes, with some of them involved in the immune system, have also been identified in humans across different social environments such as social isolation or low socioeconomic status (23–25). Moreover, the expression of genes involved in immune response as well as transcriptional regulation and cell proliferation has been shown to be sensitive to social regulation [more precisely, the level of loneliness, (26)]. Therefore, not only ELA, but also acute social stress is likely to impact gene expression in humans as it has already been proven in mice, in which the vascular system and inflammatory pathways were mainly affected (27). Furthermore, several studies have indicated an association between expression changes of diverse genes in mouse brain and social fear (28) as well as anxiety (29, 30). In humans, an investigation of the blood transcriptome has suggested altered immune function in generalized anxiety disorder (31). In addition, expression differences of α-endomannosidase (MANEA) are associated with SAD and panic disorder in human blood (32). Moreover, RNA-seq (33) has identified higher ITM2B gene expression levels associated with higher anxiety scores in a cohort of 25 monozygotic (MZ) twins, which has been validated in a second cohort of 22 MZ twins (33).
A molecular link between ELA and the sensitivity to social stress on the transcriptome level has been shown in mouse brain tissue, where distinct transcriptional patterns depending on ELA in socially stressed adult mice have been revealed: Several genes involved in, among others, cell differentiation and nervous system development have been identified as differentially expressed in mice with social stress in adulthood and ELA compared to controls without ELA. Interestingly, their gene expression levels have not been altered when exposed to either ELA or adult social stress alone (19). However, in humans the association between transcriptional changes induced by ELA and adult SAD still remain elusive.
In the current study, we aimed to identify gene expression patterns associated with SAD, ELA, and their interaction on a transcriptome-wide level in order to identify in a hypothesis-free approach molecular pathways induced by ELA which could contribute to the development of SAD later in life.
2. Materials and methods
2.1. Study population
In total, 159 participants of German ancestry (not genetically verified) between 19 and 50 years of age took part in the study. Assuming a medium effect size of 0.25, the sample size is sufficient to detect a significant effect at p < 0.05 with a power > 0.8. Power calculation has been carried out using GPower 3.1.9.2 (34) and applying an ANOVA test (fixed effects, special, main effects and interactions, no. of groups = 4). All participants were assessed using the Structured Clinical Interview for DSM-IV (SCID) and 70 participants were found to be suffering from SAD as a primary diagnosis. The severity of social anxiety was evaluated using the Liebowitz Social Anxiety Scale [LSAS, (35)]. ELA was assessed using the Childhood Trauma Questionnaire (CTQ) that measures five dimensions (further referred to as subscales) of maltreatment: emotional and physical neglect and emotional, physical, and sexual abuse (36, 37). Responses are measured on a 5-point Likert scale (1 = never true, 2 = rarely true, 3 = sometimes true, 4 = often true, 5 = very often true). Each subscale is represented by five questions with a score range from 5 to 25. Participants with at least a moderate score in one of the five categories (sexual abuse: >8; physical abuse: >10; physical neglect: >10; emotional abuse: >13; emotional neglect: >15) (Supplementary Table S1) were classified as participants with high levels of ELA (37, 38). Thus, four groups emerged: (1) control participants without SAD and low levels of ELA (n = 62), (2) control participants without SAD and high levels of ELA (n = 27), (3) participants suffering from SAD with low levels of ELA (n = 43), and (4) participants suffering from SAD with high levels of ELA (n = 27). All participants gave written informed consent to the experimental procedure prior to inclusion in the study. The study was performed in accordance with the Declaration of Helsinki and approved by the University of Tübingen local ethics committee.
2.2. RNA extraction, library preparation, and 3’RNA-sequencing
Total RNA from whole blood stored in PAXgene Blood RNA tubes was extracted using the PaxGene Blood miRNA kit (Qiagen, Hilden, Germany). Quality of RNA was assessed using a Bioanalyzer (Agilent, Santa Clara, USA). Only samples with an RNA integrity number (RIN) of 7 and higher were used for sequencing library preparation. Libraries for 3′ RNA-seq were prepared using the 3′ method by Lexogen (39) as used in the NGS Competence Center Tübingen (NCCT) where both library preparation and sequencing in randomized batches was performed. First strand synthesis of polyA-tailed RNA from total RNA using oligo dT primers was followed by degradation of the RNA template, second strand synthesis with random primers containing 5′ Illumina-compatible linker sequences, and amplification using random primers that add barcodes and cluster generation sequences (39). The libraries were sequenced on the NCCT Nova sequencing platform at a depth of about 10 million reads with 100 bp in length.
2.3. RNA-seq reads preprocessing
Read preprocessing was performed using the Lexogen pipeline (40) implementing the bbduk tool from the BBTools suite1 for quality trimming and the STAR aligner (41) that has been shown to consistently provide a well-balanced ratio between a high accuracy and sensitivity compared to runtime as well as computational resources (42) for mapping to the reference genome (vGRCh38.104). Principal component analysis (PCA) was performed using normalized reads (standard DESeq2 read normalization) to detect sample outliers by DESeq2 (43). The R package OUTRIDER (44) was used to identify gene count outliers that were excluded from further analyses (Supplementary Table S2).
To control for the effect of blood cell type composition variability on gene expression, blood cell type proportions were estimated using the granulator package in R using TPM (transcripts per million) normalized counts. Benchmarking in granulator was performed using reference cell type counts of a subset of the cohort (Supplementary Figure S2). The R package variancePartition (45) was used to calculate the variance explained by differential cell type composition and covariates. The package implements a linear mixed model method to characterize the contribution of selected variables to transcriptional variability. As deconvolution results showed a minor contribution of most cell types to the variance between the samples (Supplementary Figure S3), we used an adjustment approach of the gene counts to all cell type ratios resulting from the deconvolution approach based on a linear model adapted from Jones et al. (46) instead of using the cell type ratios as covariates in the later analyses.
2.4. Data analyses
2.4.1. Statistical analysis
All statistical analyses were performed using the software environment R. Statistical tests, that are available within the R package ggpubr, were used depending on the analysis specified in the following sections.
2.4.2. Demographic and clinical information
Normality of data was tested using Shapiro–Wilk test. The test revealed non-normal distributions for all variables (Supplementary Table S3). Therefore, the comparison of the trait medians between the independent groups was performed using the Wilcoxon Mann–Whitney rank sum test.
2.4.3. Differential gene expression
Differential gene expression (DGE) analysis was performed using the R package DESeq2 (43), which analyzes differences in gene expression based on a negative binomial generalized linear model. Cell type adjusted read counts were normalized by DESeq2’s median of ratios method (47). Genes with low counts were removed and only those with at least 20 counts in all samples were kept, as huge on/off changes were not expected due to the research question. A linear model with the factors of interest SAD, ELA and covariates age and sex was fitted. The Benjamini–Hochberg (FDR) procedure was used to correct for multiple testing. Differentially expressed genes (DEG) were considered statistically significant at FDR < 0.1 (DESeq2 default setting, as we did not expect a large number (≤ 100) of genes to be differentially expressed) with further visual inspection of the count distribution of the DEGs. Results were filtered for DEGs with an absolute log2 fold-change larger than 0.3.
2.4.4. Weighted gene co-expression network analysis
The R package WGCNA uses a validated principle called guilt by association, which relies on the assumption that associated or interacting genes share expression patterns and are likely to function together. Scale-free co-expression networks were constructed using WGCNA that defines modules using a dynamic tree-cutting algorithm based on hierarchical clustering of expression values (minimum module size = 100, cutting height = 0.99). WGCNA was performed using filtered (≥ 20 counts per sample) and variance stabilized count data (generated from the read count matrix using DESeq2’s getVarianceStabilizedData function). The network was constructed at a soft power of 10 at which the scale-free topology fit index reached 0.9. The module eigenvalue was used to perform the correlation analysis with the variables (i.e., questionnaire scores of LSAS and its subcategories as well as CTQ and its subcategories; covariates age and sex) with each whole module. Modules additionally significantly correlating with sex and age were discarded from further analyses.
2.4.5. Gene functional enrichment analysis
Gene list functional enrichment analysis was performed using the R package gProfiler2 (48, 49) by using the Gene Ontology (GO) resource [vOBO 1.4, (50, 51)], the Kyoto Encyclopedia Genes and Genomes (KEGG) pathways database [v103.0, (52)] and Reactome database [v81, (53)]. Terms with FDR-corrected p values of <0.05 were considered significantly enriched within modules.
2.4.6. Network analysis and visualization
MAPK3, the top hit of the DGE analysis, was imported into the online Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database v11.52 (54) for known and predicted protein–protein interactions (ppi). We used the following conditions for network generation: medium confidence (0.4), maximum 50 interactors for the first shell and 10 for the second shell.
The interactome of MAPK3 together with all genes of significant WGCNA modules (in total 1815) was generated using the STRING database (v11.5) starting with a full network (edges indicating both functional and physical protein associations) and then filtering for interaction scores >0.9, thereby increasing confidence. For the final interactome, all direct neighbors of MAPK3 were selected. The interactome was visualized using Cytoscape [v3.9.1., (55)].
3. Results
3.1. Demographic and clinical information
Table 1 shows the sample characteristics with respect to the four groups emerging from the factors SAD and ELA in more detail. While there was neither a significant group difference in age (Wilcoxon test, n = 159, W = 3,370, p = 0.38, with a mean age of 26 ± 7 in the SAD group and 26 ± 5 in the control group) nor sex (Pearson’s Chi-square test, χ2 = 1.62, p = 0.20, with 50 female and 20 male participants in the SAD group and 55 female and 34 male participants in the control group) with respect to SAD, significant differences in age (Wilcoxon test, n = 159, W = 3,398, p = 0.040, with a mean age of 28 ± 8 in the group with high levels of ELA and 25 ± 5 in the group with low levels of ELA), but not sex (Pearson’s Chi-square test, χ2 = 0.05, p = 0.82, 35 female and 19 male participants in the group with high levels of ELA and 70 female and 35 male participants in the group with low levels of ELA) emerged with respect to ELA. Additionally, Levene’s test revealed variance heterogeneity of the age data among the ELA groups (DF = 1, F = 14.814, p < 0.001).
Table 1. Sample characteristics for the four groups emerging from the factors SAD and ELA.
The total score of the CTQ of our cohort was not significantly different with respect to sex (Wilcoxon test, n = 159, W = 2373.5, p = 0.093), but it correlated positively with age (r = 0.117, p < 0.001). Nevertheless, we included age (in addition to sex) as covariate in the differential gene expression (DGE) analysis and tested each candidate gene expression count for correlation with age to exclude age effects on the expression of the respective gene (Supplementary Figure S3). For the total score of the LSAS, there was no significant difference with respect to sex (Wilcoxon test, n = 159, W = 3,337, p = 0.068) and no correlation with age (r = −0.001, p = 0.83). Finally, there was a significant correlation between the total scores of the LSAS and CTQ (r = 0.104, p < 0.001). This correlation was mainly due to the highly significant correlation of the emotional CTQ subscales emotional abuse (r = 0.107, p < 0.001, n = 29) and neglect (r = 0.092, p < 0.001, n = 35) with the LSAS score (Supplementary Figure S4), whereas the other subscales of the CTQ did not or less significantly correlate with the LSAS total score (physical abuse: r = 0.018, p = 0.05; sexual abuse: r = 0.018, p = 0.06; physical neglect: r = 0.004, p = 0.010). Importantly, in our cohort, we have eight cases of physical and five cases of sexual abuse only (Supplementary Table S1), which needs to be kept in mind when interpreting the results.
3.2. Differentially expressed genes were observed in SAD, but not in ELA condition
Investigating gene expression of all participants, visual inspection of the PCA revealed no obvious grouping of samples (Supplementary Figure S5). This is in line with rather subtle gene expression changes that one may expect in blood in the context of mental disorders (56, 57). Analyzing DGE using DESeq2 (43), 13 significantly (FDR-corrected p ≤ 0.1 for multiple correction) differentially expressed genes (DEGs) were identified which had a |l2fc| ≥ 0.3 in SAD condition compared to healthy individuals without SAD (Figure 1A; Supplementary Table S4), while no significant differences in expression were identified in ELA condition compared to individuals without the experience of ELA (Supplementary Table S5).
Figure 1. DGE analysis between SAD and control group using normalized gene counts adjusted to cell type ratio. (A) Volcano plot displaying the significantly (FDR-corrected p ≥ 0.1) differentially (|l2fc| ≥ 0.3) expressed genes with downregulated genes marked in blue and upregulated genes marked in red. (B) The expression patterns of the most significantly differentially expressed gene MAPK3 displayed in gene counts with respect to SAD and ELA show increased gene counts in the SAD groups without an influence of ELA. Wilcoxon rank sum test was applied and p values were adjusted for multiple testing using Benjamini–Hochberg correction. Kruskal-Wallis-Test additionally shows significantly different MAPK3 expression between the SAD/ no SAD (with respect to ELA) groups.
Visualizing the count distribution of all SAD associated DEGs, eight of the candidates exhibited an expression pattern rather caused by extreme values or other effects than being a true DEG. After removing the sample(s) exhibiting gene-specific extreme values (e.g., for MX1, Supplementary Figure S6; Supplementary Table S4), those genes no longer appeared as a DEG. Therefore, we excluded these genes from further analyses (marked in Supplementary Table S4). The remaining DEGs include (in order of significance) MAPK3 (Mitogen-Activated Protein Kinase 3), ANAPC1 (Anaphase Promoting Complex Subunit 1), PFKL (Phosphofructokinase, Liver Type), FGFBP2 (Fibroblast Growth Factor Binding Protein 2), and AC008937.2 (long non-coding (lnc) RNA, Supplementary Figure S7). The most significantly expressed gene MAPK3 (FDR-corrected p = 0.003, l2fc = 0.33) was upregulated in the SAD group compared to control individuals without SAD (Figure 1B). Figure 1B shows the counts of MAPK3 for each experimental group revealing that there is no ELA specific expression pattern within the SAD group. SAD groups with high and low levels of ELA displayed an equally high MAPK3 count (SAD/ELA: adj. Mean gene count = 132.59 ± 44.71, SAD/no ELA: adj. Mean gene count = 135.62 ± 38.59), whereas the groups of individuals without SAD showed a significantly lower MAPK3 mean count, no matter whether ELA levels were high or low (no SAD/ELA: adj. Mean gene count = 101.21 ± 20.60, no SAD/no ELA: adj. Mean gene count = 110.49 ± 38.06).
In addition to the comparison of ELA groups, the CTQ subscales were classified according to Bernstein and Fink (58), with a score moderate and higher indicating the respective trauma (Supplementary Table S1). DGE analysis was carried out for each CTQ subscale (Supplementary Table S5). Furthermore, we analyzed DGE in the SAD group only with respect to ELA and each subscale as well as in the ELA group with respect to SAD to examine the potential transcriptomic association of ELA and SAD (Supplementary Table S5). The subscale DGE analyses resulted in two genes differentially expressed in individuals with or without the experience of childhood sexual abuse in the entire cohort (Supplementary Table S6) and 197 genes with respect to physical abuse within the SAD group only (Supplementary Table S7). However, as there were only five cases of sexual abuse in the entire cohort and four cases of physical abuse in the SAD group only (Supplementary Table S1), the DGE analysis in those subgroups has to be interpreted with caution. In summary, there were no relevant significant DEGs identified for ELA as well as the CTQ subscales neither in the entire cohort nor in the SAD group only.
3.3. Gene co-expression clusters correlated with emotional ELA, but not SAD
WGCNA was performed on gene counts matching the same filters as for the DGE analysis and a soft threshold power of 10 to identify ELA and/or SAD specific gene co-expression (for more details on the analysis, see Supplementary Figure S8). We identified 11 gene co-expression modules correlating with any of the variables available for the cohort (Supplementary Table S8) with sizes ranging from 73 to 1750 genes. 1559 genes were assigned as not correlated (module gray).
Interestingly, whereas the DGE analysis revealed only associations of gene expression and SAD, the WGCNA resulted only in modules significantly correlated with ELA but not SAD or LSAS scores, respectively. In more detail, the modules red (Supplementary Table S9), greenyellow (Supplementary Table S10), and turquoise (Supplementary Table S11) were significantly correlated with the ELA groups (turquoise), the subscales emotional abuse and emotional neglect as well as total CTQ score (red and greenyellow), respectively, but not with any of the non-disease-related variables (age, sex, size and weight, Table 2). None of the ELA co-expression cluster top hub genes (Supplementary Table S12) overlapped with SAD DEGs. Furthermore, MAPK3 was found in the gray module containing the genes not correlated with any variable. The other relevant DEGs were found in the following co-expression modules: the green module that was not significantly correlated with any variable (Supplementary Figure S9) included ANAPC1 and FGFBP2. PFKL was found in the blue module which is associated with, among others, ELA, sex, and size (Supplementary Figure S9). Finally, AC008937.2 was co-expressed with genes in the red module.
Table 2. WGCNA module gene count and correlation coefficients (r) and corresponding p-values (p) for each variable.
3.4. Gene functional enrichment analysis revealed relevance of signal transduction pathways and immune system
A protein–protein-interaction (ppi) network of MAPK3 was generated using the STRING database to find pathways in which co-expressed genes within ELA-specific WGCNA modules and MAPK3, associated with SAD, potentially interact. A functional enrichment analysis was performed to examine the enrichment of annotated terms within the three modules significantly correlating with ELA and/or the respective CTQ (sub-)scale modules (turquoise, red, and greenyellow) and the MAPK3 ppi network. The MAPK3 ppi network was enriched mainly for MAPK signal transduction pathways and NTRK (neurotrophin receptor) signaling. The red and the turquoise WGCNA modules were enriched for cellular structural processes/compartments (Figures 2A,B). The greenyellow co-expression cluster contained genes particularly involved in immune-related pathways (especially interleukin regulation and production) and JAK–STAT signaling (Figure 2C).
Figure 2. Gene functional enrichment of the (A) WGCNA module red, (B) WGCNA module greenyellow, and (C) WGCNA module turquoise. Significance values are color-coded. Abbreviations indicate the database with G: KEGG database, BP: Biological process (GO term), CC: Cellular component (GO term) and MF: Molecular function (GO term).
The enrichment analysis did not reveal any shared or overlapping pathways between the modules. The results therefore do not indicate a direct molecular mediation of ELA on adult SAD by one single process.
3.5. Network analysis identified common genes between SAD-related MAPK3 and ELA associated co-expression modules
The gene lists of each co-expression module correlating with ELA and/or CTQ (sub)scales, i.e., of the turquoise, red, and greenyellow modules, were compared with the genes contained in the MAPK3 ppi network that plays a role in SAD to identify overlapping and thus potentially interacting genes. We identified PTPN7 (Tyrosine-protein phosphatase non-receptor type 7) as a gene present in the MAPK3 ppi network and the red module (Supplementary Table S13). PTPN7 is a member of the phosphatase family and a known negative regulator of MAPK signal cascade activation (59). The turquoise module and the MAPK3 ppi network share a set of 20 genes (Supplementary Table S13). Gene set enrichment analysis identified an involvement of most of those genes especially in signal transduction (MAPK, NTRK, neurotrophin, Supplementary Figure S10). The hub gene of the greenyellow module STAT3 (Signal transducer and activator of transcription 3) as well as RAF1 (RAF Proto-Oncogene Serine/Threonine-Protein Kinase) were also found in the MAPK3 ppi network (Supplementary Table S13). RAF1 activation initiates a mitogen-activated protein kinase cascade and is in part regulated by cytokine signaling (60) and STAT3 mediates cellular responses to interleukins and other growth factors (61–67) as well as inflammatory responses by regulating differentiation of naive CD4+ T-cells into T-helper Th17 or regulatory T-cells (68). Therefore, both genes are involved in the immune response and are linked to the mitogen-activated signaling cascade (60, 69, 70), where MAPK3 plays a central role (60). Therefore, an interaction of STAT3, RAF1, and MAPK3 in immune signaling is likely.
To verify the interaction of MAPK3 and genes from the ELA-correlated modules, we used the STRING database to extract information on interaction scores of the 1815 genes (MAPK3 + genes from the three modules). After removing all genes with an interaction score ≤ 0.9 indicating the highest confidence, 51 genes remained. We generated an interactome of those genes (Figure 3). The interactome highlighted the interaction of MAPK3 and the before-mentioned STAT3, RAF1, and PTPN7. However, the interactome also revealed complex interrelations between the genes with several to many interaction partners of each gene (Figure 3).
Figure 3. Network visualization of the interactome of MAPK3 and the genes from the turquoise, red, and greenyellow module. The nodes were colored by module, node size displays module membership score from WGCNA, the node sort and node transparency were set by STRING degree and the edge transparency was set by STRING score.
The functional enrichment of the gene list was performed by using the Reactome database only. This enabled focusing on interaction of the genes to form a biologically relevant network. The analysis revealed mainly enrichment of the genes in terms related to immune-related signaling (Figure 4, Signaling by Receptor Tyrosine Kinase, Signaling by NTRKs, etc.). However, the most significantly enriched pathway was Signal Transduction (Figure 4).
Figure 4. Gene functional enrichment of the MAPK3 interactome associated with genes of the turquoise, red, and greenyellow module. Significance values are color-coded.
4. Discussion
Transcriptome analyses have become highly relevant over recent years to investigate the molecular basis of psychiatric disorders. In particular, RNA-seq has been widely used to analyze psychiatric disorders and interrelations (71–73). In the study presented here, we focused on social anxiety disorder and the molecular connection with a potential environmental trigger—early life adversity.
DGE analysis revealed genes associated with SAD, with MAPK3 being the most significantly upregulated in individuals with SAD compared to control individuals. No DEGs were identified between individuals with and without a history of ELA. MAPK3 is a serine/threonine kinase which acts as an essential component of the MAP kinase signal transduction pathway. MAPK1/ERK2 and MAPK3/ERK1 play an important role in the MAPK/ERK cascade (extracellular signal-regulated kinase-dependent cascade). Depending on the cellular context, the MAPK/ERK cascade mediates diverse biological functions such as cell growth, adhesion, survival, and differentiation through the regulation of transcription, translation, and cytoskeletal rearrangements. The MAPK/ERK cascade also plays a role in initiation and regulation of meiosis, mitosis, and postmitotic functions in differentiated cells by phosphorylating a number of transcription factors that, for example, promote breast cancer (74). Interestingly, another finding was the differential expression of the lncRNA AC0008937.2 in the context of SAD, which—as an antisense lncRNA to MAP3K1—has a potentially regulative role (75) in the MAPK signaling cascade. Another important gene that can regulate and is regulated by MAPK signaling is the brain-derived neurotrophic factor (BDNF) (76, 77), which is an important player in neurogenesis (78). The neurotrophin BDNF, which is known to be affected by neuroinflammation (79), is essential for dendritic development in peripheral and central nervous system and regulates dendritic growth (80). Moreover, BDNF level changes in serum are associated with anxiety disorders (81, 82). Furthermore, BDNF protein levels were shown to be reduced in brain tissue of rats with a history of ELA which were subsequently exposed to stress (83). However, in our data derived from whole blood, we did not identify a SAD-specific expression pattern of BDNF. MAPK3 may function as a transporter in blood to regulate the expression of BDNF in brain tissue which, in turn, may lead to altered structural brain plasticity playing a role in SAD. Therefore, the analysis of BDNF expression levels would be of interest in different brain areas in the context of SAD. MAPK3 expression might also be altered in the brain of patients as MAPK phosphorylation levels in the amygdala were directly associated with anxiety symptoms in a previous study (84). The authors demonstrate that the rate of extracellular signal-related kinase phosphorylation in the amygdala is negatively and independently associated with anxiety symptoms (84). These findings further support our results of an involvement of MAPK signaling in SAD. Nevertheless, as BDNF is associated with several mental disorders, MAPK3 differential expression may also not be restricted to social anxiety. An upregulation of MAPK-related genes is also found in Major Depressive Disorder [MDD, (85)] especially in cases with a history of ELA. Altered MAPK signaling involved in the development of mental disorders may be a result of stress on the mental or even cellular (e.g., infection) level. Indeed, psychological stressors such as early childhood adversities or other factors such as infections are shown to severely affect the immune system resulting in an inflammatory phenotype and increase the risk for adult psychiatric disorders (10, 86, 87). Therefore, ELA may cause (sterile) inflammation (88) and might lead to molecular alterations mainly involved in inflammatory processes and those changes may induce an aberrant synaptic development in SAD transferred by MAPK signaling (89, 90).
Gene co-expression analysis revealed gene clusters significantly associated with the emotional aspects of ELA (emotional abuse and neglect). Although no direct link on the level of differential gene expression was identified, this is an interesting finding as social anxiety and especially the emotional forms of childhood maltreatment are shown to significantly correlate (91, 92), which is supported by our data as well (Supplementary Figure S3). Therefore, a connection of emotional ELA and adult SAD on the molecular level seems likely. The MAPK3 interaction network shared one gene, PTPN7, with the red co-expression module, that significantly correlated with the scores of the CTQ subscales emotional abuse and emotional neglect as well as the total CTQ score. PTPN7 is a member of the phosphatase family and specifically inactivates MAPKs. Nothing is known about the regulation of PTPN7 in the context of ELA or SAD so far, although Schwieck et al. (2020) did not identify differential expression in MDD cases together with suicide risk and ELA history (85). However, it acts as a regulator of MAPK signaling activity (93, 94). The immune system is a plausible pathway how ELA could be molecularly involved in adult psychiatric disorders as ELA is known to cause inflammatory mimicking effects—or sterile inflammation—that are still measurable in adults, e.g., higher levels of typical markers of inflammation such as white blood cell count, circulating proinflammatory cytokine levels, and the acute phase molecule C-Reactive Protein (CRP) and lower NK cell activity (10, 95). Although we did not measure cytokine levels or NK cell activity, none of the mentioned marker genes associated with emotional neglect and/or abuse were differing in the context of ELA and/or SAD in our sample (data not shown). Nevertheless, PTPN7 may be a promising candidate connecting ELA, the immune response, and MAPK signaling as a potential regulator of anxiety disorders. Furthermore, blood cell type ratios were estimated and gene counts were adjusted to the ratios reflecting the immune activity. However, cell type ratios did not differ with respect to ELA or SAD in our sample.
The greenyellow co-expression module significantly correlated with emotional abuse and neglect in early childhood and was enriched for immune-related terms and JAK–STAT signaling. STAT3 and RAF1 were shared between this module and the MAPK3 ppi network. RAF1 is a known upstream regulator of MAPK signaling [Raf/MEK/ERK cascade, (96)]. RAF1 expression increases upon infection, which is mediated by interleukin 2 [IL-2, (97)], whereas inhibition of RAF1 affects production of IL-6 and IL-8 in cultured human corneal epithelial cells (98). STAT3—a transcription factor—regulates processes involved in inflammation and tumorigenesis by regulating cell proliferation, differentiation and metabolism (99). STAT3 is a member of the JAK–STAT signaling pathway whose canonical mode is based on cytokine release followed by MAPK signaling activation (100). In the non-canonical signaling pathway, unphosphorylated STATs are localized on heterochromatin in the nucleus in association with proteins regulating the maintenance of heterochromatin state (101). Therefore, STAT3—like MAPK3 and RAF1, respectively—is involved in a complex crosstalk of signaling pathways and may be involved in epigenetic regulation of downstream processes. A STAT3 knockout in mice leads to reduced negative behavioral reactivity (102) Additionally, STAT3 is involved in alcohol withdrawal (103) and depressive symptoms in rats (104, 105). Therefore, RAF1 and STAT3 are potential candidates connecting ELA, immune response, and SAD.
Gene set functional enrichment of the genes overlapping in the turquoise module, that significantly correlated with ELA, and the MAPK3 ppi network revealed mainly terms related to signal transduction pathways (i.e., MAPK and NTRK and neurotrophin signaling, Supplementary Figure S10). Among others, NTRK signaling was enriched, pointing toward a role of those genes in BDNF-related processes. Moreover, the gene set enrichment of the MAPK3 interactome substantiates the role of NTRK signaling and therefore in (neuro)inflammation in the association of ELA and SAD.
4.1. Summary and model
As stated above, ELA is known to have an effect on neuronal structures by affecting the immune system (e.g., cytokine levels) via (neuro)inflammation (106), which is linked to structural changes involved in the development of mental disorders. Our findings support this assumption on the molecular level: The gene set enrichment of the co-expression clusters revealed terms especially involved in the cellular structure, signal transduction, and immune response. Genes co-expressed in clusters associated especially with emotional ELA are potential interactors of MAPK3, which is significantly differentially expressed in individuals suffering from SAD compared to controls. Especially STAT3, the hub gene of such a co-expression cluster, may be regulated by the emotional ELA-dependent release of interleukins like IL-6 (107, 108) and thus may be involved in the cell type-specific regulation of more growth factor and cytokine release (109–111) which for their part increase MAPK3 expression (112). MAPK3 may be further involved in the expression of genes shaping synaptic plasticity, e.g., BDNF.
In a recent study of our group, the blood DNA methylome was analyzed in the same cohort presented here. Differentially methylated regions (DMRs) specifically associated with SAD, ELA, or the interaction of SAD and ELA were identified (113). None of these regions were overlapping with the DEGs found in the current gene expression analyses. However, STAT3 is shown to interact with the DNA methyltransferase DNMT1 (114, 115). Therefore, STAT3-directed DNA methylation is a possible step in the signal transduction cascade transferring ELA to SAD. In the study of Camilo et al. (2019), a gene network approach revealed a direct association of MAPK3 methylation and cocaine use disorder (116). Therefore, in a follow-up study, we aim to conduct a multivariate machine learning-based analysis to integrate DNA methylation and gene expression data in the context of SAD following ELA.
The approach presented here has several limitations: First, the expression of genes can vary between blood cell types and therefore, differential cell type composition can affect results. Immune responsive cell types are known to play a role in mental disorders such as MDD or panic disorder (117–119), which make them unsuitable for usage as covariate in statistical tests and DGE analysis. Therefore, we adjusted the gene expression data to the estimated cell type composition with the help of reference cell counts of a subset of our cohort and showed that the number of real counts and estimated ratios correlated for several relevant cell types (Supplementary Figure S1). Furthermore, we have to be aware of the fact that transcriptomic profiles are not only cell type- but also tissue-specific, and that we therefore cannot assume that the differences we observe in blood directly reflect the situation in brain (as mentioned for BDNF earlier). In psychiatric transcriptomics, we are faced with the problem, that the tissue of interest—the brain—is not easily available for transcriptomic analyses in living individuals. However, we can assume that there is some overlap of genes expressed similar in blood and brain, as human whole blood tissue showed a significant similarity in gene expression to multiple brain tissues with a median correlation of 0.5 as revealed by microarray analysis (120). Moreover, in rat brain and blood tissue more than half of the 29,215 genes analyzed by microarray were co-expressed (121). Furthermore, age differed significantly between the groups with high and low levels of ELA. Age effects were not identified for the expression of the candidate genes (Supplementary Figure S3) and co-expression clusters correlating with age were excluded from further analyses. An approach in a larger cohort would be needed to decipher whether age-dependent gene expression has an effect on the results presented here, as well as to identify further environmental factors shared by the analyzed groups’ influencing gene expression. An enlarged sample size as well as less heterogeneity could also increase statistical power of the analyses as we neither identified significant DEGs in an ELA condition nor gene co-expression clusters significantly correlating with SAD. Moreover, the phosphorylation levels of MAPK3 and MAPK signaling (and therefore the activation of the signal transduction cascade) mark an important step in the signal cascade that should be included in future experiments to clarify whether gene product or phosphorylation abundance are the potential drivers behind the molecular development of SAD. Also, genetic risk variants may confer susceptibility to anxiety traits by influencing gene expression. Therefore, in future studies, genetic data should be included in the analyses in order to determine to what extent expression of a certain gene is influenced by genetics or shaped by an environmental stressor such as ELA, especially as gene x environment (G x E) interactions are shown to play a major role in anxiety disorders [(e.g., 122)].
In summary, by investigating gene expression in the context of SAD and its relation to ELA on a transcriptome level, we were able to identify DEGs associated with SAD—with MAPK3 being the most significant DEG—as well as co-expression clusters correlating with ELA and/or its subclasses. Interestingly, functional enrichment of MAPK3 protein–protein interaction network and ELA associated gene co-expression modules pointed toward signal transduction pathways and the immune system. Additionally, shared genes are involved in JAK–STAT and ERK signaling as well as DNA methylation. Although a direct molecular link of ELA leading to adult SAD by gene expression changes was not identified, the data indicate an indirect relation of emotional ELA and SAD mediated by the interaction of genes involved in immune-related signal transduction. Further studies will be needed to replicate our findings in independent, larger cohorts and to investigate the potential effect of the immune responsive gene expression pattern caused by ELA on adult anxiety disorders in more detail.
Data availability statement
The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (123) and are accessible through GEO Series accession number GSE228702 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE228702).
Ethics statement
The studies involving human participants were reviewed and approved by Ethics committee of the University of Tübingen (process number 526/2018BO2). The participants provided their written informed consent to participate in this study.
Author contributions
AF, BK, and VN received funding for the study (DFG NI 1332/6-1, FA 361/27-1, KR 4398/5-1). SE, BK, and VN were involved in designing the experiments. AW and MM recruited and phenotyped the participants. SE and SP ran the experiments. SE and TH analyzed the data. SE and VN wrote the manuscript. AF and JS-H commented on the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was funded by the German Research Foundation (DFG NI 1332/6-1, FA 361/27-1, KR 4398/5-1).
Acknowledgments
We acknowledge the support by the NGS Competence Center Tübingen and the Open Access Publishing Fund of University of Tuebingen. The authors also wish to express their appreciation to all participants.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpsyt.2023.1125553/full#supplementary-material
Footnotes
References
1. Kessler, RC, McGonagle, K, Zhao, S, Nelson, CB, Hughes, M, Eshleman, S, et al. Lifetime and 12-month prevalence of DSM-III-R psychiatric disorders in the United States: results from the National Comorbidity Survey. Arch Gen Psychiatry. (1994) 51:8–19. doi: 10.1001/archpsyc.1994.03950010008002
2. WHO World Mental Health Survey CollaboratorsStein, DJ, Lim, CCW, Roest, AM, de Jonge, P, Aguilar-Gaxiola, S, et al. The cross-national epidemiology of social anxiety disorder: data from the world mental health survey initiative. BMC Med. (2017) 15:143. doi: 10.1186/s12916-017-0889-2
3. American Psychiatric Association. Diagnostic and Statistical Manual of Mental Disorders (DSM-IV). Washington, DC: American Psychiatric Association. (1994).
4. Stein, MB, Chen, CY, Jain, S, Jensen, KP, He, F, Heeringa, SG, et al. Genetic risk variants for social anxiety. Am J Med Genet B Neuropsychiatr Genet. (2017) 174:470–82. doi: 10.1002/ajmg.b.32541
5. Meier, SM, and Deckert, J. Genetics of anxiety disorders. Curr Psychiatry Rep. (2019) 21:16. doi: 10.1007/s11920-019-1002-7
6. Smith, KE, and Pollak, SD. Rethinking concepts and categories for understanding the neurodevelopmental effects of childhood adversity. Perspect Psychol Sci. (2021) 16:67–93. doi: 10.1177/1745691620920725
7. Hettema, JM, Neale, MC, and Kendler, KS. A review and meta-analysis of the genetic epidemiology of anxiety disorders. Am J Psychiatry. (2001) 158:1568–78. doi: 10.1176/appi.ajp.158.10.1568
8. Klauke, B, Deckert, J, Reif, A, Pauli, P, and Domschke, K. Life events in panic disorder-an update on "candidate stressors". Depress Anxiety. (2010) 27:716–30. doi: 10.1002/da.20667
9. McGrath, LM, Weill, S, Robinson, EB, Macrae, R, and Smoller, JW. Bringing a developmental perspective to anxiety genetics. Dev Psychopathol. (2012) 24:1179–93. doi: 10.1017/S0954579412000636
10. Elwenspoek, MMC, Hengesch, X, Leenen, FAD, Schritz, A, Sias, K, Schaan, VK, et al. Proinflammatory T cell status associated with early life adversity. J Immunol. (2017) 199:4046–55. doi: 10.4049/jimmunol.1701082
11. Heim, C, Newport, DJ, Mletzko, T, Miller, AH, and Nemeroff, CB. The link between childhood trauma and depression: insights from HPA axis studies in humans. Psychoneuroendocrinology. (2008) 33:693–710. doi: 10.1016/j.psyneuen.2008.03.008
12. Kuhlman, KR, Chiang, JJ, Horn, S, and Bower, JE. Developmental psychoneuroendocrine and psychoneuroimmune pathways from childhood adversity to disease. Neurosci Biobehav Rev. (2017) 80:166–84. doi: 10.1016/j.neubiorev.2017.05.020
13. Danese, A, and Baldwin, JR. Hidden wounds? Inflammatory links between childhood trauma and psychopathology. Annu Rev Psychol. (2017) 68:517–44. doi: 10.1146/annurev-psych-010416-044208
14. Dantzer, R. Neuroimmune interactions: from the brain to the immune system and vice versa. Physiol Rev. (2018) 98:477–504. doi: 10.1152/physrev.00039.2016
15. Maccari, S, Krugers, HJ, Morley-Fletcher, S, Szyf, M, and Brunton, PJ. The consequences of early-life adversity: neurobiological, behavioural and epigenetic adaptations. J Neuroendocrinol. (2014) 26:707–23. doi: 10.1111/jne.12175
16. Nemeroff, CB. Early-life adversity, CRF dysregulation, and vulnerability to mood and anxiety disorders. Psychopharmacol Bull. (2004) 38:14–20.
17. Cole, SW, Conti, G, Arevalo, JMG, Ruggiero, AM, Heckman, JJ, and Suomi, SJ. Transcriptional modulation of the developing immune system by early life social adversity. Proc Natl Acad Sci U S A. (2012) 109:20578–83. doi: 10.1073/pnas.1218253109
18. Levine, ME, Cole, SW, Weir, DR, and Crimmins, EM. Childhood and later life stressors and increased inflammatory gene expression at older ages. Soc Sci Med. (2015) 130:16–22. doi: 10.1016/j.socscimed.2015.01.030
19. Peña, CJ, Smith, M, Ramakrishnan, A, Cates, HM, Bagot, RC, Kronman, HG, et al. Early life stress alters transcriptomic patterning across reward circuitry in male and female mice. Nat Commun. (2019) 10:5098. doi: 10.1038/s41467-019-13085-6
20. Schwaiger, M, Grinberg, M, Moser, D, Zang, JCS, Heinrichs, M, Hengstler, JG, et al. Altered stress-induced regulation of genes in monocytes in adults with a history of childhood adversity. Neuropsychopharmacology. (2016) 41:2530–40. doi: 10.1038/npp.2016.57
21. Dieckmann, L, Cole, S, and Kumsta, R. Stress genomics revisited: gene co-expression analysis identifies molecular signatures associated with childhood adversity. Transl Psychiatry. (2020) 10:34. doi: 10.1038/s41398-020-0730-0
22. Lutz, P-E, Chay, M-A, Pacis, A, Chen, GG, Aouabed, Z, Maffioletti, E, et al. Non-CG methylation and multiple histone profiles associate child abuse with immune and small GTPase dysregulation. Nat Commun. (2021) 12:1132. doi: 10.1038/s41467-021-21365-3
23. Cole, SW. Social regulation of human gene expression: mechanisms and implications for public health. Am J Public Health. (2013) 103 Suppl 1:S84–92. doi: 10.2105/AJPH.2012.301183
24. Miller, GE, Chen, E, Fok, AK, Walker, H, Lim, A, Nicholls, EF, et al. Low early-life social class leaves a biological residue manifested by decreased glucocorticoid and increased proinflammatory signaling. Proc Natl Acad Sci U S A. (2009) 106:14716–21. doi: 10.1073/pnas.0902971106
25. Chen, E, Miller, GE, Walker, HA, Arevalo, JM, Sung, CY, and Cole, SW. Genome-wide transcriptional profiling linked to social class in asthma. Thorax. (2009) 64:38–43. doi: 10.1136/thx.2007.095091
26. Cole, SW, Hawkley, LC, Arevalo, JM, Sung, CY, Rose, RM, and Cacioppo, JT. Social regulation of gene expression in human leukocytes. Genome Biol. (2007) 8:R189. doi: 10.1186/gb-2007-8-9-r189
27. Stankiewicz, AM, Goscik, J, Majewska, A, Swiergiel, AH, and Juszczak, GR. The effect of acute and chronic social stress on the hippocampal Transcriptome in mice. PLoS One. (2015) 10:e0142195. doi: 10.1371/journal.pone.0142195
28. Royer, M, Pai, B, Menon, R, Bludau, A, Gryksa, K, Perry, RBT, et al. Transcriptome and chromatin alterations in social fear indicate association of MEG3 with successful extinction of fear. Mol Psychiatry. (2022) 27:4064–76. doi: 10.1038/s41380-022-01481-2
29. Malan-Müller, S, Fairbairn, L, Daniels, WMU, Dashti, MJS, Oakeley, EJ, Altorfer, M, et al. Molecular mechanisms of D-cycloserine in facilitating fear extinction: insights from RNAseq. Metab Brain Dis. (2016) 31:135–56. doi: 10.1007/s11011-015-9727-4
30. Czibere, L, Baur, LA, Wittmann, A, Gemmeke, K, Steiner, A, Weber, P, et al. Profiling trait anxiety: transcriptome analysis reveals cathepsin B (Ctsb) as a novel candidate gene for emotionality in mice. PLoS One. (2011) 6:e23604. doi: 10.1371/journal.pone.0023604
31. Wingo, AP, and Gibson, G. Blood gene expression profiles suggest altered immune function associated with symptoms of generalized anxiety disorder. Brain Behav Immun. (2015) 43:184–91. doi: 10.1016/j.bbi.2014.09.016
32. Jensen, KP, Stein, MB, Kranzler, HR, Yang, BZ, Farrer, LA, and Gelernter, J. The alpha-endomannosidase gene (MANEA) is associated with panic disorder and social anxiety disorder. Transl Psychiatry. (2014) 4:e353. doi: 10.1038/tp.2013.122
33. Davies, MN, Verdi, S, Burri, A, Trzaskowski, M, Lee, M, Hettema, JM, et al. Generalised anxiety disorder:a twin study of genetic architecture, genome-wide association and differential Gene expression. PLoS One. (2015) 10:e0134865. doi: 10.1371/journal.pone.0134865
34. Erdfelder, E, Faul, F, and Buchner, A. GPOWER: a general power analysis program. Behav Res Methods. (1996) 28:11.
36. Bernstein, DP, et al. Development and validation of a brief screening version of the childhood trauma questionnaire. Child Abuse Negl. (2003) 27:169–90. doi: 10.1016/s0145-2134(02)00541-0
37. Wingenfeld, K, et al. The German version of the childhood trauma questionnaire (CTQ): Preliminary psychometric properties. Psychother Psychosom Med Psychol. (2010) 60:442–50. doi: 10.1055/s-0030-1247564
38. Bernstein, DP, Stein, JA, Newcomb, MD, Walker, E, Pogge, D, Ahluvalia, T, et al. Development and validation of a brief screening version of the childhood trauma questionnaire. Child Abuse Negl. (2003) 27:169–90. doi: 10.1016/S0145-2134(02)00541-0
39. Ma, F, Fuqua, BK, Hasin, Y, Yukhtman, C, Vulpe, CD, Lusis, AJ, et al. A comparison between whole transcript and 3' RNA sequencing methods using Kapa and Lexogen library preparation methods. BMC Genomics. (2019) 20:9. doi: 10.1186/s12864-018-5393-3
40. Lexogen. QuantSeq 3‘ mRNA-Seq Integrated Data Analysis Pipelines on BlueBee® Genomics Platform. Lexogen. (2021).
41. Dobin, A, Davis, CA, Schlesinger, F, Drenkow, J, Zaleski, C, Jha, S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635
42. Musich, R, Cadle-Davidson, L, and Osier, MV. Comparison of short-read sequence aligners indicates strengths and weaknesses for biologists to consider. Front Plant Sci. (2021) 12:657240. doi: 10.3389/fpls.2021.657240
43. Love, MI, Huber, W, and Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
44. Brechtmann, F, Mertes, C, Matusevičiūtė, A, Yépez, VA, Avsec, Ž, Herzog, M, et al. OUTRIDER: a statistical method for detecting aberrantly expressed genes in RNA sequencing data. Am J Hum Genet. (2018) 103:907–17. doi: 10.1016/j.ajhg.2018.10.025
45. Hoffman, GE, and Schadt, EE. variancePartition: interpreting drivers of variation in complex gene expression studies. BMC Bioinformatics. (2016) 17:483. doi: 10.1186/s12859-016-1323-z
46. Jones, MJ, Islam, SA, Edgar, RD, and Kobor, MS. Adjusting for cell type composition in DNA methylation data using a regression-based approach. Methods Mol Biol. (2017) 1589:99–106. doi: 10.1007/7651_2015_262
47. Anders, S, and Huber, W. Differential expression analysis for sequence count data. Genome Biol. (2010) 11:R106. doi: 10.1186/gb-2010-11-10-r106
48. Kolberg, L, Raudvere, U, Kuzmin, I, Vilo, J, and Peterson, H. gprofiler2 -- an R package for gene list functional enrichment analysis and namespace conversion toolset g:profiler. F1000Res. (2020):9. doi: 10.12688/f1000research.24956.2
49. Raudvere, U, Kolberg, L, Kuzmin, I, Arak, T, Adler, P, Peterson, H, et al. G:profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. (2019) 47:W191–8. doi: 10.1093/nar/gkz369
50. Ashburner, M, Ball, CA, Blake, JA, Botstein, D, Butler, H, Cherry, JM, et al. Gene ontology: tool for the unification of biology. Gene Ontol Consort Nat Genet. (2000) 25:25–9. doi: 10.1038/75556
51. Gene Ontology, C. The Gene Ontology resource: enriching a GOld mine. Nucleic Acids Res. (2021) 49:D325–34. doi: 10.1093/nar/gkaa1113
52. Kanehisa, M, and Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27
53. Gillespie, M, Jassal, B, Stephan, R, Milacic, M, Rothfels, K, Senff-Ribeiro, A, et al. The reactome pathway knowledgebase 2022. Nucleic Acids Res. (2022) 50:D687–92. doi: 10.1093/nar/gkab1028
54. Szklarczyk, D, Gable, AL, Nastou, KC, Lyon, D, Kirsch, R, Pyysalo, S, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. (2021) 49:D605–12. doi: 10.1093/nar/gkaa1074
55. Shannon, P, Markiel, A, Ozier, O, Baliga, NS, Wang, JT, Ramage, D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
56. Cole, JJ, McColl, A, Shaw, R, Lynall, ME, Cowen, PJ, de Boer, P, et al. No evidence for differential gene expression in major depressive disorder PBMCs, but robust evidence of elevated biological ageing. Transl Psychiatry. (2021) 11:404. doi: 10.1038/s41398-021-01506-4
57. Nøhr, AK, Lindow, M, Forsingdal, A, Demharter, S, Nielsen, T, Buller, R, et al. A large-scale genome-wide gene expression analysis in peripheral blood identifies very few differentially expressed genes related to antidepressant treatment and response in patients with major depressive disorder. Neuropsychopharmacology. (2021) 46:1324–32. doi: 10.1038/s41386-021-01002-9
58. Bernstein, DP, and Fink, L. Childhood Trauma Questionnaire: A Retrospective Self-Report. San Antonio, TX: The Psychological Corporation (1998).
59. Inamdar, VV, Reddy, H, Dangelmaier, C, Kostyak, JC, and Kunapuli, SP. The protein tyrosine phosphatase PTPN7 is a negative regulator of ERK activation and thromboxane generation in platelets. J Biol Chem. (2019) 294:12547–54. doi: 10.1074/jbc.RA119.007735
60. Santarpia, L, Lippman, SM, and El-Naggar, AK. Targeting the MAPK-RAS-RAF signaling pathway in cancer therapy. Expert Opin Ther Targets. (2012) 16:103–19. doi: 10.1517/14728222.2011.645805
61. Shen, S, Niso-Santano, M, Adjemian, S, Takehara, T, Malik, SA, Minoux, H, et al. Cytoplasmic STAT3 represses autophagy by inhibiting PKR activity. Mol Cell. (2012) 48:667–80. doi: 10.1016/j.molcel.2012.09.013
62. Jiang, H, Yu, J, Guo, H, Song, H, and Chen, S. Upregulation of survivin by leptin/STAT3 signaling in MCF-7 cells. Biochem Biophys Res Commun. (2008) 368:1–5. doi: 10.1016/j.bbrc.2007.04.004
63. Saxena, NK, Vertino, PM, Anania, FA, and Sharma, D. Leptin-induced growth stimulation of breast cancer cells involves recruitment of histone acetyltransferases and mediator complex to CYCLIN D1 promoter via activation of Stat3. J Biol Chem. (2007) 282:13316–25. doi: 10.1074/jbc.M609798200
64. Dreuw, A, Radtke, S, Pflanz, S, Lippok, BE, Heinrich, PC, and Hermanns, HM. Characterization of the signaling capacities of the novel gp130-like cytokine receptor. J Biol Chem. (2004) 279:36112–20. doi: 10.1074/jbc.M401122200
65. Shao, H, Cheng, HY, Cook, RG, and Tweardy, DJ. Identification and characterization of signal transducer and activator of transcription 3 recruitment sites within the epidermal growth factor receptor. Cancer Res. (2003) 63:3923–30.
66. Yamamoto, T, Sekine, Y, Kashima, K, Kubota, A, Sato, N, Aoki, N, et al. The nuclear isoform of protein-tyrosine phosphatase TC-PTP regulates interleukin-6-mediated signaling pathway through STAT3 dephosphorylation. Biochem Biophys Res Commun. (2002) 297:811–7. doi: 10.1016/S0006-291X(02)02291-X
67. Tsai, Y, Su, YH, Fang, SS, Huang, TN, Qiu, Y, Jou, YS, et al. Etk, a Btk family tyrosine kinase, mediates cellular transformation by linking Src to STAT3 activation. Mol Cell Biol. (2000) 20:2043–54. doi: 10.1128/MCB.20.6.2043-2054.2000
68. Ma, L, Huang, C, Wang, XJ, Xin, DE, Wang, LS, Zou, QC, et al. Lysyl oxidase 3 is a dual-specificity enzyme involved in STAT3 Deacetylation and Deacetylimination modulation. Mol Cell. (2017) 65:296–309. doi: 10.1016/j.molcel.2016.12.002
69. Jain, N, et al. Repression of Stat3 activity by activation of mitogen-activated protein kinase (MAPK). Oncogene. (1998) 17:10. doi: 10.1038/sj.onc.1202238
70. Yang, L, Lin, S, Xu, L, Lin, J, Zhao, C, and Huang, X. Novel activators and small-molecule inhibitors of STAT3 in cancer. Cytokine Growth Factor Rev. (2019) 49:10–22. doi: 10.1016/j.cytogfr.2019.10.005
71. Ramaker, RC, Bowling, KM, Lasseigne, BN, Hagenauer, MH, Hardigan, AA, Davis, NS, et al. Post-mortem molecular profiling of three psychiatric disorders. Genome Med. (2017) 9:72. doi: 10.1186/s13073-017-0458-5
72. Li, X, and Teng, S. RNA sequencing in schizophrenia. Bioinform Biol Insights. (2015) 9:53–60. doi: 10.4137/BBI.S28992
73. Hemmings, SMJ, Swart, P, Womersely, JS, Ovenden, ES, van den Heuvel, LL, McGregor, NW, et al. RNA-seq analysis of gene expression profiles in posttraumatic stress disorder, Parkinson’s disease and schizophrenia identifies roles for common and distinct biological pathways. Discover Mental Health. (2022) 2:9. doi: 10.1007/s44192-022-00009-y
74. Deng, R, Zhang, HL, Huang, JH, Cai, RZ, Wang, Y, Chen, YH, et al. MAPK1/3 kinase-dependent ULK1 degradation attenuates mitophagy and promotes breast cancer bone metastasis. Autophagy. (2021) 17:3011–29. doi: 10.1080/15548627.2020.1850609
75. Liu, B, Xiang, W, Liu, J, Tang, J, Wang, J, Liu, B, et al. The regulatory role of antisense lncRNAs in cancer. Cancer Cell Int. (2021) 21:459. doi: 10.1186/s12935-021-02168-4
76. Medina, JH, and Viola, H. ERK1/2: a key cellular component for the formation, retrieval, reconsolidation and persistence of memory. Front Mol Neurosci. (2018) 11:361. doi: 10.3389/fnmol.2018.00361
77. Genheden, M, Kenney, JW, Johnston, HE, Manousopoulou, A, Garbis, SD, and Proud, CG. BDNF stimulation of protein synthesis in cortical neurons requires the MAP kinase-interacting kinase MNK1. J Neurosci. (2015) 35:972–84. doi: 10.1523/JNEUROSCI.2641-14.2015
78. Li, Y, Li, F, Qin, D, Chen, H, Wang, J, Wang, J, et al. The role of brain derived neurotrophic factor in central nervous system. Front Aging Neurosci. (2022) 14:986443. doi: 10.3389/fnagi.2022.986443
79. Lima Giacobbo, B, Doorduin, J, Klein, HC, Dierckx, RAJO, Bromberg, E, and de Vries, EFJ. Brain-derived Neurotrophic factor in brain disorders: focus on Neuroinflammation. Mol Neurobiol. (2019) 56:3295–312. doi: 10.1007/s12035-018-1283-6
80. Nourbakhsh, K, and Yadav, S. Kinase signaling in dendritic development and disease. Front Cell Neurosci. (2021) 15:624648. doi: 10.3389/fncel.2021.624648
81. Suliman, S, Hemmings, SM, and Seedat, S. Brain-derived Neurotrophic factor (BDNF) protein levels in anxiety disorders: systematic review and meta-regression analysis. Front Integr Neurosci. (2013) 7:55. doi: 10.3389/fnint.2013.00055
82. Farias, T, Cerqueira, RA, Sousa, DF, Freire, JVC, Lopes, ACT, and Lima De Moura Cal, SF. BDNF Protein and Anxiety Disorders In: K Shad and K Dogan, editors. Neurological and Mental Disorders. London: intechOpen (2020)
83. Zhang, X, Li, H, Sun, H, Jiang, Y, Wang, A, Kong, Y, et al. Effects of BDNF signaling on anxiety-related behavior and spatial memory of adolescent rats in different length of maternal separation. Front Psych. (2020) 11:709. doi: 10.3389/fpsyt.2020.00709
84. de Carvalho, CR, Lopes, MW, Constantino, LC, Hoeller, AA, de Melo, HM, Guarnieri, R, et al. The ERK phosphorylation levels in the amygdala predict anxiety symptoms in humans and MEK/ERK inhibition dissociates innate and learned defensive behaviors in rats. Mol Psychiatry. (2021) 26:7257–69. doi: 10.1038/s41380-021-01203-0
85. Schiweck, C, Claes, S, van Oudenhove, L, Lafit, G, Vaessen, T, de Beeck, GO, et al. Childhood trauma, suicide risk and inflammatory phenotypes of depression: insights from monocyte gene expression. Transl Psychiatry. (2020) 10:296. doi: 10.1038/s41398-020-00979-z
86. Elwenspoek, MMC, Kuehn, A, Muller, CP, and Turner, JD. The effects of early life adversity on the immune system. Psychoneuroendocrinology. (2017) 82:140–54. doi: 10.1016/j.psyneuen.2017.05.012
87. Costello, H, Gould, RL, Abrol, E, and Howard, R. Systematic review and meta-analysis of the association between peripheral inflammatory cytokines and generalised anxiety disorder. BMJ Open. (2019) 9:e027925. doi: 10.1136/bmjopen-2018-027925
88. Kuhlman, KR, Horn, SR, Chiang, JJ, and Bower, JE. Early life adversity exposure and circulating markers of inflammation in children and adolescents: a systematic review and meta-analysis. Brain Behav Immun. (2020) 86:30–42. doi: 10.1016/j.bbi.2019.04.028
89. Jadhav, SP, Kamath, SP, Choolani, M, Lu, J, and Dheen, ST. microRNA-200b modulates microglia-mediated neuroinflammation via the cJun/MAPK pathway. J Neurochem. (2014) 130:388–401. doi: 10.1111/jnc.12731
90. Kaminska, B, Gozdz, A, Zawadzka, M, Ellert-Miklaszewska, A, and Lipko, M. MAPK signal transduction underlying brain inflammation and gliosis as therapeutic target. Anat Rec (Hoboken). (2009) 292:1902–13. doi: 10.1002/ar.21047
91. Kuo, JR, Goldin, PR, Werner, K, Heimberg, RG, and Gross, JJ. Childhood trauma and current psychological functioning in adults with social anxiety disorder. J Anxiety Disord. (2011) 25:467–73. doi: 10.1016/j.janxdis.2010.11.011
92. Lochner, C, Seedat, S, Allgulander, C, Kidd, M, Stein, D, and Gerdner, A. Childhood trauma in adults with social anxiety disorder and panic disorder: a cross-national study. Afr J Psychiatry (Johannesbg). (2010) 13:376–81. doi: 10.4314/ajpsy.v13i5.63103
93. Seo, H, Lee, IS, Park, JE, Park, SG, Lee, DH, Park, BC, et al. Role of protein tyrosine phosphatase non-receptor type 7 in the regulation of TNF-alpha production in RAW 264.7 macrophages. PLoS One. (2013) 8:e78776. doi: 10.1371/journal.pone.0078776
94. Sabio, G, and Davis, RJ. TNF and MAP kinase signalling pathways. Semin Immunol. (2014) 26:237–45. doi: 10.1016/j.smim.2014.02.009
95. Witek Janusek, L, Tell, D, Albuquerque, K, and Mathews, HL. Childhood adversity increases vulnerability for behavioral symptoms and immune dysregulation in women with breast cancer. Brain Behav Immun. (2013) 30:S149–62. doi: 10.1016/j.bbi.2012.05.014
96. Krzyzowska, M, Swiatek, W, Fijalkowska, B, Niemialtowski, M, and Schollenberger, A. The role of map kinases in immune response. Adv Cell Biol. (2010) 2:125–38. doi: 10.2478/v10052-010-0007-5
97. Turner, B, Tonks, NK, Rapp, UR, and Reed, JC. Interleukin 2 regulates Raf-1 kinase activity through a tyrosine phosphorylation-dependent mechanism in a T-cell line. Proc Natl Acad Sci U S A. (1993) 90:4.
98. Zhao, GQ, Lin, J, Hu, LT, Yin, XN, Wang, Q, Xu, Q, et al. The role of Dectin-1/Raf-1 signal cascade in innate immune of human corneal epithelial cells against Aspergillus fumigatus infection. Int J Ophthalmol. (2016) 9:1371–5. doi: 10.18240/ijo.2016.10.01
99. Qi, QR, and Yang, ZM. Regulation and function of signal transducer and activator of transcription 3. World J Biol Chem. (2014) 5:231–9. doi: 10.4331/wjbc.v5.i2.231
100. Li, WX. Canonical and non-canonical JAK-STAT signaling. Trends Cell Biol. (2008) 18:545–51. doi: 10.1016/j.tcb.2008.08.008
101. Seif, F, Khoshmirsafa, M, Aazami, H, Mohsenzadegan, M, Sedighi, G, and Bahar, M. The role of JAK-STAT signaling pathway and its regulators in the fate of T helper cells. Cell Commun Signal. (2017) 15:23. doi: 10.1186/s12964-017-0177-y
102. Reisinger, SN, Sideromenos, S, Horvath, O, Derdak, S, Cicvaric, A, Monje, FJ, et al. STAT3 in the dorsal raphe gates behavioural reactivity and regulates gene networks associated with psychopathology. Mol Psychiatry. (2021) 26:2886–99. doi: 10.1038/s41380-020-00904-2
103. Chen, WY, Chen, H, Hamada, K, Gatta, E, Chen, Y, Zhang, H, et al. Transcriptomics identifies STAT3 as a key regulator of hippocampal gene expression and anhedonia during withdrawal from chronic alcohol exposure. Transl Psychiatry. (2021) 11:298. doi: 10.1038/s41398-021-01421-8
104. Kong, E, Sucic, S, Monje, FJ, Reisinger, SN, Savalli, G, Diao, W, et al. STAT3 controls IL6-dependent regulation of serotonin transporter function and depression-like behavior. Sci Rep. (2015) 5:9009. doi: 10.1038/srep09009
105. Kwon, SH, Han, JK, Choi, M, Kwon, YJ, Kim, SJ, Yi, EH, et al. Dysfunction of microglial STAT3 alleviates depressive behavior via neuron-microglia interactions. Neuropsychopharmacology. (2017) 42:2072–86. doi: 10.1038/npp.2017.93
106. Andersen, SL. Neuroinflammation, early-life adversity, and brain development. Harv Rev Psychiatry. (2022) 30:24–39. doi: 10.1097/HRP.0000000000000325
107. Carpenter, LL, Gawuga, CE, Tyrka, AR, Lee, JK, Anderson, GM, and Price, LH. Association between plasma IL-6 response to acute stress and early-life adversity in healthy adults. Neuropsychopharmacology. (2010) 35:2617–23. doi: 10.1038/npp.2010.159
108. Garbers, C, Aparicio-Siegmund, S, and Rose-John, S. The IL-6/gp130/STAT3 signaling axis: recent advances towards specific inhibition. Curr Opin Immunol. (2015) 34:75–82. doi: 10.1016/j.coi.2015.02.008
109. Eddahri, F, Denanglaire, S, Bureau, F, Spolski, R, Leonard, WJ, Leo, O, et al. Interleukin-6/STAT3 signaling regulates the ability of naive T cells to acquire B-cell help capacities. Blood. (2009) 113:2426–33. doi: 10.1182/blood-2008-04-154682
110. Huang, W, Yeh, HH, Lin, CC, Lai, W-W, Chang, J-Y, Chang, W-T, et al. Signal transducer and activator of transcription 3 activation up-regulates interleukin-6 autocrine production: a biochemical and genetic study of established cancer cell lines and clinical isolated human cancer cells. Mol Cancer. (2010) 9:16. doi: 10.1186/1476-4598-9-309
111. Agerholm, R, Rizk, J, Viñals, MT, and Bekiaris, V. STAT3 but not STAT4 is critical for gammadeltaT17 cell responses and skin inflammation. EMBO Rep. (2019) 20:e48647. doi: 10.15252/embr.201948647
112. Zhang, W, and Liu, H. MAPK signal pathways in the regulation of cell proliferation in mammalian cells. Cell Res. (2002) 12:9. doi: 10.1038/sj.cr.7290105
113. Wiegand, A, Kreifelts, B, Munk, MHJ, Geiselhart, N, Ramadori, KE, MacIsaac, JL, et al. DNA methylation differences associated with social anxiety disorder and early life adversity. Transl Psychiatry. (2021) 11:104. doi: 10.1038/s41398-021-01225-w
115. Zhang, Q, et al. MAPK signal pathways in the regulation of cell proliferation in mammalian cells. Proc Natl Acad Sci USA. (2005) 102:5.
116. Camilo, C, Maschietto, M, Vieira, HC, Tahira, AC, Gouveia, GR, Feio dos Santos, AC, et al. Genome-wide DNA methylation profile in the peripheral blood of cocaine and crack dependents. Braz J Psychiatry. (2019) 41:485–93. doi: 10.1590/1516-4446-2018-0092
117. Baek, JH, Kim, HJ, Fava, M, Mischoulon, D, Papakostas, GI, Nierenberg, A, et al. Reduced venous blood basophil count and anxious depression in patients with major depressive disorder. Psychiatry Investig. (2016) 13:321–6. doi: 10.4306/pi.2016.13.3.321
118. Schleifer, SJ, Keller, SE, and Bartlett, JA. Panic disorder and immunity: few effects on circulating lymphocytes, mitogen response, and NK cell activity. Brain Behav Immun. (2002) 16:698–705. doi: 10.1016/S0889-1591(02)00022-3
119. Medina-Rodriguez, EM, Lowell, JA, Worthen, RJ, Syed, SA, and Beurel, E. Involvement of innate and adaptive immune systems alterations in the pathophysiology and treatment of depression. Front Neurosci. (2018) 12:547. doi: 10.3389/fnins.2018.00547
120. Sullivan, PF, Fan, C, and Perou, CM. Evaluating the comparability of gene expression in blood and brain. Am J Med Genet B Neuropsychiatr Genet. (2006) 141B:261–8. doi: 10.1002/ajmg.b.30272
121. Witt, SH, Sommer, WH, Hansson, AC, Sticht, C, Rietschel, M, and Witt, CC. Comparison of gene expression profiles in the blood, hippocampus and prefrontal cortex of rats. Silico Pharmacol. (2013) 1:15. doi: 10.1186/2193-9616-1-15
122. Su, X, Li, W, Lv, L, Li, X, Yang, J, Luo, XJ, et al. Transcriptome-wide association study provides insights into the genetic component of Gene expression in anxiety. Front Genet. (2021) 12:740134. doi: 10.3389/fgene.2021.740134
Keywords: RNA sequencing, gene expression, social anxiety disorder, early life adversity, transcriptome, immune system, signal transduction
Citation: Edelmann S, Wiegand A, Hentrich T, Pasche S, Schulze-Hentrich JM, Munk MHJ, Fallgatter AJ, Kreifelts B and Nieratschker V (2023) Blood transcriptome analysis suggests an indirect molecular association of early life adversities and adult social anxiety disorder by immune-related signal transduction. Front. Psychiatry. 14:1125553. doi: 10.3389/fpsyt.2023.1125553
Edited by:
Gabriela Xavier, Harvard Medical School, United StatesReviewed by:
Vanessa Kiyomi Ota, Federal University of São Paulo, BrazilGustavo Satoru Kajitani, Federal University of São Paulo, Brazil
Carolina Muniz Carvalho, Universidade Federal de São Paulo, Brazil
Copyright © 2023 Edelmann, Wiegand, Hentrich, Pasche, Schulze-Hentrich, Munk, Fallgatter, Kreifelts and Nieratschker. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Vanessa Nieratschker, vanessa.nieratschker@med.uni-tuebingen.de