First, we studied gene- and pathway-level associations with lifetime depression of 2,376 genes annotated to 20 curated immune pathways (Suppl.Table.1) in 13,309 individuals from GS. This sample contained 1,867 depression cases with a mean age of 47.5 (± 12.8) years and high proportion of females (70.8%), and a control group of 11,442 individuals of similar age (Table 1). Using the multi-model analysis in MAGMA, we identified one immune pathway gene, the growth hormone receptor (GHR; FDR = 0.03, Z = 4.2), and one immune pathway, the Reactome “RUNX1 regulates transcription of genes involved in differentiation of myeloid cells” (FDR = 0.016, beta = 1.2), associated with depression. Three (CSF2, RUNX1 and PRKCB) of the seven genes annotated to this pathway were tagged as potentially responsible (p < 0.05) for the observed association. The SNP-level association results for GHR are shown in Fig. 1A. We also show in Table 2 the top results (FDR < 0.25) from the gene association analysis that were not considered statistically significant. These included 10 genes annotated to the cytokine signaling, cellular responses to stimuli, interferon signaling and adaptive immune system pathways, among others. Results from the pathway-level association analysis with depression are shown in Fig. 1B.
Second, we investigated the association of gene-CM interactions with depression using the CTQ in a subsample of 749 GS participants (mean age: 52.3 ± 10.2 years, 58.3% females), including 99 depression cases, and an independent sample of 509 individuals from the population-based cohort of BD (mean age: 54.1 ± 8 years, 49.7% females), including 96 depression cases (Table 1). The analyses using the different CTQ susbcales in GS identified 122 immune pathway genes that interacted with the total EA scores in association with depression (Suppl.Table.2), 127 that interacted with total PA (Suppl.Table.3), 416 with total SA (Suppl.Table.4), 10 with total EN (Suppl.Table.5) and 42 with total PN (Suppl.Table.6). However, these analyses in BD resulted in almost only nominal associations that did not survive correction for multiple comparisons. Therefore, we decided to look-up the significant findings from GS in BD using a threshold of p < 0.05. With this, we found 18 overlaps for EA, 14 for PA, 34 for SA, 2 for EN and 7 for PN between both samples (Table 3, Fig. 2, Suppl.Tables.2–6). The most significant GS genes interacting with each CTQ subscale that showed nominal significance in BD were: ASB15 in EA (FDRGS=0.0026, pBD=0.0131), TNIK in PA (FDRGS=7.7x10− 4, pBD=0.002), SERPINH1 in SA (FDRGS=2.1x10− 4, pBD=0.03), FYN in EN (FDRGS=2.3x10− 4, pBD=0.0423) and IKBKE in PN (FDRGS=0.0033, pBD=0.0041).
After, we combined the p-values obtained for the CTQ subscales into one p-value that should reflect a general association of gene-CM interactions with depression, regardless of the type of trauma. We found 653 significant (FDR < 0.05) genes in GS (Suppl.Table.7). Because we had adopted the look-up approach using a p < 0.05 criterium for BD before, we applied the same strategy here and found 167 overlaps with GS (Suppl.Table.7). Importantly, we found that 56 from these genes were statistically significant (FDR < 0.05) also in BD (Tables 3–4, Fig. 2) and, therefore, these gene-CM interactions in association with depression were considered replicated in our study.
In addition, we used PPI network analysis to better inform the biological context of the 56 replicated gene-CM interactions favoring depression (i.e. query proteins). The resulting network (Fig. 1C) had a PPI enrichment p = 2.18x10− 7, 138 PPI pairs, and 37 connected proteins from a total of 66 nodes. The hub nodes in the network were all non-query proteins: IFNG and STAT1 (with 20 interactions each), followed by RBX1 (18 interactions) and MAPK1 (16 interactions). Enrichment analysis identified 78 GO_BPs and 31 tissues of expression that were overrepresented (FDR < 0.05, strength > 1) in our network (Suppl.Table.8). Here, we were interested in observing enrichments of biological processes that were different from those pathways initially selected for our study, as well as of specific tissues and cell types expressing the identified genes. As expected, most of the GO_BP terms overrepresented in our network analysis were related to the original Reactome pathways included in the study, such as cellular responses to biotic stimulus, response to interferons, regulation of innate immune response and chemokine production, etc. However, we found some GO_BPs that provided more defined insights and some hints of effects beyond peripheral immunity including, for example, Astrocyte differentiation (FDR = 6.11x10− 5), Cellular response to amyloid-beta (FDR = 0.0026), Regulation of long-term synaptic potentiation (FDR = 0.0035), Response to L-glutamine (FDR = 0.014), and various terms related to oxidative stress (e.g. Response to reactive oxygen species, Positive regulation of nitric-oxide synthase biosynthetic process, Positive regulation of oxidative stress-induced cell death, Positive regulation of reactive oxygen species metabolic process, etc.), the production of particular interleukins (IL-6, 1β, 17, -12 and − 23), and stem cell differentiation, including Regulation of hematopoietic progenitor cell differentiation (FDR = 0.028) and Regulation of pro-B cell differentiation (FDR = 0.0305). In general, the network proteins showed enriched expression in peripheral macrophages and microglial cells (Suppl.Table.8). Nevertheless, we also found overrepresentation of other types of white blood cells, including B-lymphocytes, natural killer cells, helper and memory T-lymphocytes, monocytes and neutrophils. Interestingly, there was also enrichment for expression in the blood-brain barrier (BBB; FDR = 0.0124).
Finally, we performed pathway-level association analyses using the CTQ interaction gene results. However, we found only nominal associations in GS that did not survive correction for multiple comparisons (Fig. 1B). These were with the “Inflammasomes” (p = 0.0044, beta = 0.6) and “Cellular responses to stimuli” (p = 0.04, beta = 0.084) pathways in the SA analysis, the “Class I MHC-mediated antigen processing & presentation” pathway in the EN (p = 0.019, beta = 0.13) and EA (p = 0.024, beta = 0.12) analyses, the “Adaptive immune system” (p = 0.038, beta = 0.084) and Signaling by BCR (p = 0.042, beta = 0.18) pathways in the PN analysis, and the “Cellular responses to stimuli” (p = 0.048, beta = 0.078) pathway in the PA analysis. After meta-analysis of CTQ subscales, the top pathway for this CM analysis in GS was the “Class I MHC-mediated antigen processing & presentation” (p = 5.8x10− 4, FDR = 0.058). This was followed by the “Adaptive immune system” (p = 0.0168), Signaling by BCR (p = 0.0215) and “Cellular responses to stimuli” (p = 0.0476) pathways. In BD, there were no significant findings from the pathway-level analyses.