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.

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.

Read preprocessing was performed using the Lexogen pipeline (40) implementing the bbduk tool from the BBTools suite¹ 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.

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, χ² = 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, χ² = 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).

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.

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).

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.

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.

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).

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.

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).

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).

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.