FoxO1, A2M, and TGF-β1: three novel genes predicting depression in gene X environment interactions are identified using cross-species and cross-tissues transcriptomic and miRNomic analyses.

To date, gene-environment (GxE) interaction studies in depression have been limited to hypothesis-based candidate genes, since genome-wide (GWAS)-based GxE interaction studies would require enormous datasets with genetics, environmental, and clinical variables. We used a novel, cross-species and cross-tissues "omics" approach to identify genes predicting depression in response to stress in GxE interactions. We integrated the transcriptome and miRNome profiles from the hippocampus of adult rats exposed to prenatal stress (PNS) with transcriptome data obtained from blood mRNA of adult humans exposed to early life trauma, using a stringent statistical analyses pathway. Network analysis of the integrated gene lists identified the Forkhead box protein O1 (FoxO1), Alpha-2-Macroglobulin (A2M), and Transforming Growth Factor Beta 1 (TGF-β1) as candidates to be tested for GxE interactions, in two GWAS samples of adults either with a range of childhood traumatic experiences (Grady Study Project, Atlanta, USA) or with separation from parents in childhood only (Helsinki Birth Cohort Study, Finland). After correction for multiple testing, a meta-analysis across both samples confirmed six FoxO1 SNPs showing significant GxE interactions with early life emotional stress in predicting depressive symptoms. Moreover, in vitro experiments in a human hippocampal progenitor cell line confirmed a functional role of FoxO1 in stress responsivity. In secondary analyses, A2M and TGF-β1 showed significant GxE interactions with emotional, physical, and sexual abuse in the Grady Study. We therefore provide a successful 'hypothesis-free' approach for the identification and prioritization of candidate genes for GxE interaction studies that can be investigated in GWAS datasets.

Introduction

Major Depressive Disorder (MDD) has a well-established genetic contribution, although with a modest (30–40%) heritability [1]. However, we have so far failed to identify vulnerability genes for this disorder. Indeed, while genome wide association studies (GWAS) have identified main genetic effects for schizophrenia [2], autism [3], and bipolar disorder [4], even very large GWAS meta-analyses have failed to identify genome-wide significant associations in depression [5, 6]. Only recently, and using enormous datasets from more than 75 thousand individuals with depression and 200 thousand healthy subjects, 15 genome-wide significant loci for MDD have been identified [7]. One of the potential reasons behind this weak genetic effect is the fact that individual genotypic variations may increase the risk of depression only in the presence of exposure to life stressors and other adverse environmental circumstances, a phenomenon named as ‘gene-environment (GxE) interaction’ [8]. Indeed, environmental factors, and in particular exposure to adversities early in life, have been consistently implicated in the pathophysiology of depression [9]. For example, in a large study conducted by the Center of Disease Control with over 9000 participants, Chapman et al. [10] reported a dose-response relationship between the severity of experienced childhood adversities and the lifetime presence of depressive episodes or chronic depression. However, despite the epidemiological and clinical evidence confirming these effects, there is substantial variability in the outcomes of early life stress, since not all people exposed to early adversities develop depression later in life. GxE interactions may indeed also explain why adverse environmental factors do not increase the risk of depression in each and every individual.

Two of the most established examples of GxE interactions linking stress to depression result from the hypothesis-driven investigation of the serotonin transporter promoter short/long polymorphism (5-HTTLPR) and of a functional single nucleotide polymorphism (SNP) within the FK506-binding protein 51 (FKBP5) gene. Caspi and colleagues were the first to show that individuals with one or two copies of the short allele of the 5-HTTLPR exhibit more depressive symptoms and suicidality in the presence of stressful life events as compared with individuals homozygous for the long allele [11]. Although not all the studies have replicated these findings [12], this GxE interaction has been shown to predict depression in individuals exposed to childhood maltreatment [13, 14]. Similarly, many studies have shown that a functional SNP within FKBP5 (rs1360780) interacts with childhood abuse to predict a multitude of psychiatric phenotypes in adults, including depression and post-traumatic stress disorder [15-19]. For this latter gene, epigenetic mechanisms leading to different cortisol reactivity to stress seem to explain these effects [20, 21]. These two examples of GxE interactions are built on decades of hypothesis-driven research, and have been instrumental in order to identify not only genetic variants that increase the risk of depression, but also the potential biological and molecular mechanisms underlying this increased risk. However, as extensively discussed elsewhere [8], this hypothesis-driven approach can only discover a fraction of such potential existing genetic variants, and it is at odds with the current strategy of using a hypothesis-free approach to identify genes associated with disorders.

Of course, using GWAS data to test GxE interactions would be a theoretically viable option, but it would require enormous datasets with both exposure (such as life events) and outcome (depression) data. To date, only two genome-wide by environment interaction studies (GWEIS) have been performed in this regard. Dunn and colleagues (2016) [22], using data from the SHARe cohort of the Women’s Health Initiative comprising more than 10 thousand African Americans and Hispanics/Latinas, have examined genetic main effects and GxE interactions with stressful life events and social support in the development of depressive symptoms; however, only one interaction signal was genome-wide significant in African Americans (rs4652467 located 14 kb from CEP350), and it was not replicated. The second, subsequent whole genome pilot study, performed in only 320 subjects characterized for recent stressful life events, found no interaction that was genome-wide significant [23].

In the current paper, we propose a different ‘omics-based’ approach to identify candidate genes for GxE interactions studies in GWAS datasets, using cross-species and cross-tissue biological prioritization strategies that limits the number of investigated genes and thus enhances the statistical power to identify significant findings. In particular, we have first analyzed the transcriptome and miRNome data from the hippocampus of adult rats exposed to prenatal stress (PNS), a well-established model of early life stress leading to depressive behavior and hypothalamic-pituitary-adrenal (HPA) axis hyperactivity in adulthood [24], in order to obtain a list of genes that are both modulated by PNS and targeted by the miRNAs that are modulated by PNS. We have then integrated the resulting genes list with transcriptome data obtained from blood mRNA of human adults who had suffered from childhood trauma, and analyzed the overlapping genes list using network analysis, which identifies genes and molecules that interact with each other because of physical interaction or co-expression, or because they are related in common signaling pathways. Finally, we have tested the top network cluster of genes for GxE interactions in order to examine the effects of early life stress on depressive symptoms in adulthood, in two different clinical samples from the Grady Trauma Project in Atlanta (USA) [25] and the Helsinki Birth Cohort Study (Finland) [26, 27]. In vitro experiments in a human hippocampal progenitor cell line were also conducted to confirm the role of the identified genes in stress responsivity.

Material and methods

Animal model and clinical samples

Prenatal stress model (for transcriptomics and miRNomics analyses)

PNS procedure was performed as already published [28-30]; briefly, pregnant dams in the last week of gestation were restrained in a transparent Plexiglas cylinder, under bright light, for 45 min, three times a day for 1 week. Control pregnant females were left undisturbed in their home cages. Male offspring from control and PNS groups were killed at postnatal day (PND) 62 (early adulthood) for whole hippocampal dissection (for further details see Supplementary Material). Rat handling and experimental procedures were performed in accordance with the EC guidelines (EC Council Directive 86/609 1987) and with the Italian legislation on animal experimentation (D.L. 116/92), in accordance with the National Institute of Health Guide for the Care and Use of Laboratory Animals.

Subjects exposed to early life trauma (for transcriptomics analysis)

We recruited volunteer subjects from the local population living in the area served by the South London and Maudsley NHS Trust, in south-east London. Childhood traumatic events were assessed using the Childhood Experience of Care Abuse Questionnaire (CECA-Q; [31]). All the volunteers were screened using the Psychosis Screening Questionnaire (to exclude any psychotic and psychotic-like symptoms). We performed transcriptomic analyses in 20 subjects who reported at least two types of abuse (separation and physical abuse = 40%; separation and sexual abuse = 20%; physical and sexual abuse = 5%; physical and emotional abuse = 5%) or one type of severe abuse (physical abuse only = 20%; emotional abuse only = 5%; sexual abuse only = 5%); and in 20 subjects matched for age, gender and body mass index (BMI) with no history of early life trauma. Mean age ± SD was 27 ± 1.6 and 25 ± 0.9 (df = 15; ℵ2 = 0.222, p-value = 0.3), and percentage of females was 45% (7 F/13 M) and 35% (9 F/11 M) (ℵ2 = 0.417, p-value = 0.5), respectively, in the subjects with and without childhood trauma. We collected information about current medications (one subject was taking oral contraception and one was in treatment with mefenamic acid), general health (asthma, eczema, frequent headaches, recent infections), education, employment, nicotine, and cannabis smoking; the two groups were not different for any of these variables (data not showed).

Blood samples were collected by using PaxGene Blood Tubes. After collection, blood samples were kept at room temperature for 2 h, then at −20 °C for 2 days and then at −80 °C until their processing. The project was approved by the Local Ethical Committee and subjects signed a consent form (for further details see Supplementary Material).

Grady trauma project (for GxE analysis)

In the Grady Trauma Project, recruited in Atlanta (Georgia, USA), detailed trauma interviews were collected in a large sample of 4791 adults (mean age ± SD of 40.1 ± 13.9) with high rates of current and lifetime PTSD. For this study, we used the Childhood Trauma Questionnaire (CTQ), a 28-item, psychometrically validated inventory assessing self-reported levels of sexual, physical and emotional abuse [32]. The presence of moderate-severe sexual, physical and emotional abuse was coded as described before [33]. Emotional abuse was used for the main GxE analyses on early life emotional stress, across both this and the Helsinki Birth Cohort Study. Depressive symptoms were measured using the Beck Depression Inventory (BDI) [34]. The institutional review board at Emory University approved the study procedures. All participants provided written informed consent before participating.

Helsinki birth cohort study (for GxE analysis)

Between 2001 and 2004, 1620 individuals (mean age ± SD of 61.5 ± 2.9) who were born at Helsinki University Central Hospital between 1934 and 1944, and who were alive and living in Finland in 1971, participated into a clinical examination during which blood samples for DNA were used to obtain GWAS data for both directly genotyped and imputed SNPs. Depressive symptoms were measured using the BDI [34]. Early life emotional stress in this samples was defined as separation from biological parents in 1939–1946, during the second world war, when around 80,000 Finnish children were sent to Sweden and Denmark in order to escape the dangers of war; 384 subjects were identified as having been separated in childhood [26, 27, 35]. These children were separated by both parents and their siblings at the median age of 3.8 years (Range 0.6–10.1 years) for a mean duration of 1.4 years (Range 0.2–4.7 years), thus severely disrupting the attachment to significant caregivers. Previous studies in this cohort have clearly shown that these subjects show worse mental health in adulthood. For example, they report almost double the risk of experiencing depressive symptoms over time [36] and have higher risks of mental disorders, including substance abuse and personality disorders [37, 38]. Interestingly for the context of our paper, they also show higher cortisol levels and more pronounced cortisol reactivity [39]. The project was approved by the Local Ethical Committee (for further details see Supplementary Material).

Biological assessments

mRNA and miRNA isolation from hippocampus of animals and from the blood of subjects

Total RNA, including miRNAs, was isolated from rat brain tissues by using RNeasy mini kit (Qiagen, Italy) and from human blood samples using PaxGene miRNA kit (Qiagen, Italy), according to the manufacturer’s instructions. Samples were also treated with DNase, and RNA quantity and quality was assessed by evaluation of the A260/280 and A260/230 ratios using a Nanodrop spectrometer (NanoDrop Technologies).

Whole Genome expression microarray analyses in the hippocampus of animals and in the blood of subjects

Gene expression microarray assays were performed using Rat Gene 2.1st Array Strips (which covers 27,147 coding transcripts) or Human Gene 2.1st Array Strips (which covers 31,650 coding transcripts) on GeneAtlas platform (Affymetrix), following the WT Expression Kit protocol described in the Affymetrix GeneChip Expression Analysis Technical Manual, and as we have done before [28, 40] (for further details see Supplementary Material). Of note, this is the first report of the rat hippocampal transcriptome following PNS using the novel Rat Gene 2.1st Array Strips, hence we have reported here all the relevant data and pathways analyses.

Real Time validation analyses of transcriptomic findings in the hippocampus of animals and in the blood of subjects

All the genes that were significantly modulated by stress both in the hippocampus of animals and in the blood of subjects (n = 22 genes, see Results Section) were validated using Real-Time PCR by using Biorad qPCR Mix (BioRad, Milan) and Taqman Assays on a 96 wells Real Time PCR System (One Step TaqMan Real Time PCR). Each sample was assayed in triplicate and each target gene was normalized to β-actin and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) in rodents, and to β−2-microglobulin (B2M) and GAPDH in humans. The Pfaffl Method was used to determine relative target gene expression.

MiRNome expression analysis in the hippocampus of animals

Five hundred nanogram of total RNA (including miRNAs) were processed with the FlashTag Biotin HSR RNA Labeling kit (Affymetrix, Santa Clara, CA, USA) and subsequently hybridized onto the GeneChip miRNA 4.1 Array Strip on a GeneAtlas platform (Affymetrix, Santa Clara, CA, USA). The GeneChip miRNA 4.1 Array Strip shows a comprehensive coverage as they are designed to interrogate all mature miRNA sequences in miRBase Release 20. Washing/staining and scanning procedures were conducted on the Fluidics station following the manufacturer’s instructions.

Genotyping of polymorphisms within the Grady Trauma Study

DNA was extracted from saliva in Oragene collection vials (DNA Genotek, Ottawa, ON, Canada) using the DNAdvance kit (Beckman Coulter Genomics, Danvers, MA, USA), while DNA from blood was extracted using either the E.Z.N.A. Mag-Bind Blood DNA Kit (Omega Bio-Tek, Norcross, GA, USA) or ArchivePure DNA Blood Kit (5 Prime, Gaithersburg, MD, USA). N = 4,791 subjects were genotyped on the HumanOmniExpress (6%) and the Omni1-Quad BeadChip (94%) (Illumina Inc), and genotypes were called with Illumina’s Genome Studio. The HumanOmniExpress interrogates 730,525 individual SNPs per sample, whereas the Omni1- Quad BeadChip interrogates 1,011,219 individual SNPs (for further details see Supplementary Material).

Genotyping of polymorphisms in the helsinki birth cohort study

DNA was extracted from blood samples using Gentra Kit. Participants were genotyped with the modified Illumina 610 k array, at the Wellcome Trust Sanger Institute, Cambridge, UK according to standard protocols [38] (for further details see Supplementary Material).

Statistical and bioinformatics analysis

Transcriptomic and miRNome statistical analysis

Raw data were imported and analyzed with the software Partek Genomic Suite 6.6 (Partek, St. Louis, MO, USA). We checked for quality control and batch effects using the Principal Component Analyses (PCA) in the Partek Genomic Suite Software, and we did not observe any outlier or any batch effect, thus we did not correct for this. After quality control of the data, Analysis of variance (ANOVA) test was performed to assess the effects of PNS on genes expression or miRNAs expression in rats, and of history of early life adversities on gene expression in humans. A maximum filter of p < 0.05 (FDR corrected) and a minimum absolute fold change cut-off of 1.4 in animals and of 1.2 in humans was applied to select the lists of significant genes. The miRNAs-mRNAs combining analyses was performed using a specific sub-feature in Partek Genomic Suite 6.6. For blood transcriptome, possible differences in cell type between the two groups were evaluated by using CellMix package builds on R/BioConductor and available from http://web.cbio.uct.ac.za/~renaud/CRAN/web/CellMix, and as we have not detected differences in childhood trauma exposed subjects as compared to non-exposed individuals we did not use it as covariate.

Pathway analyses

List of significant genes were analyzed for pathway analyses by using Ingenuity Pathway Analyses Software (IPA) where, as a background, we used gene lists that we obtained applying the minimum absolute fold change cut-off of 1.4 in animals and of 1.2 in humans and q-value <0.05.

Animal and human integration data

The genes lists deriving from rats (following miRNAs-mRNAs combining analyses) and humans (mRNAs only) were cross-integrated, and the list of common genes was tested for random probability by Hyper-Geometric Distribution test [41] in R.

Network analyses

The list of common genes from animal and human data were used for network analyses using specific tools in IPA.

GXE interactions

Grady trauma project

We transformed imputed genotypes to best guessed genotypes using a probability threshold of 90%. To assess GxE interaction, we individually tested each SNP within the three genes of interest resulting from the final network analysis and also available in the Helsinki Birth cohort study, for interactions between emotional stress (emotional abuse) and BDI scores in adulthood, as a similar phenotype was also present in the Helsinki Birth Cohort Study. Age, gender, and the first two principal components of the identity-by-state (IBS) matrix were used as covariates in the linear regression analyses in R. We made further adjustments for SNP X covariate and child abuse X covariate interactions, as suggested by Keller [42]. In addition, secondary analyses were conducted testing GxE interaction between each type of childhood trauma (emotional, sexual and physical abuse) and BDI scores in adulthood, in this sample only.

Helsinki birth cohort study

We tested those SNPs within the three genes of interest resulting from the final network analysis that were also available in the Grady Trauma Project, for interactions between emotional stress (whether the individual was separated from parents during WWII, or not) and BDI scores in adulthood, using linear regression analysis. Each SNP was tested in a separate model; age at testing, gender, and the first three principal components of the IBS matrix, main effects of the gene and the environment were used as covariates. We made further adjustments for SNP × covariate and separation status × covariate interactions, as suggested by Keller [42].

Meta-analyses

We ran a fixed-effects meta-analysis combining the interaction results testing emotional stress in both cohorts, using the R-package Rmeta. Meta-analysis p-values were corrected for multiple testing using Bonferroni-correction over all tested SNPs.

Cell culture, gene expression and neurogenesis assay

The immortalized, multipotent human fetal hippocampal progenitor cell line, HPC0A07/03C (propriety of ReNeuron), was used for all cellular and molecular analyses (mRNA levels, neuronal proliferation). Details of this cell line, the proliferation protocol, gene expression and immunocytochemistry can be found in Supplementary Information Materials and Methods.

Results

Transcriptomic changes in the hippocampus of adult male rats following PNS

Our first aim was to identify gene expression changes in the hippocampus of male adult rats (at PND 62) that had been exposed to PNS. We analyzed the entire transcriptome and we found a significant modulation of 916 genes in the comparison with control animals (using 1.4 < FC < −1.4, q-value <0.05). We have visualized the most significant genes in a hierarchical clustering in Fig. 1, and listed the 916 significant genes in Supplementary Table 1.

Fig. 1

Hierarchical Clustering of gene expression changes by prenatal stress (PNS) in the hippocampus of rats, in comparison with control animals. The red bar on the horizontal axis indicates control animals and the orange bar indicates PNS-exposed animals. The blue and the red squares in each group (CTRL or PNS) indicate the modulation of the gene expression, with red squares indicating genes that are up-regulated and blue squares gene that are down regulated

We then performed a pathway analysis on the significantly modulated genes, in order to identify the main biological processes involved. We identified 49 biological processes that were altered by the PNS exposure, including some involved in neurodevelopment (Axonal Guidance Signaling, Protein Kinase A Signaling, Glutamate Receptor Signaling, Wnt/β-catenin Signaling, CREB Signaling in Neurons, Synaptic Long Term Depression, Synaptic Long Term Potentiation, Glucocorticoid Receptor Signaling) and inflammation (IL-8 Signaling, STAT3 Pathway, CDK5 Signaling, IL-1 Signaling, IL-6 Signaling, BMP signaling pathway and TGF-β Signaling). All the significant pathways are summarized by the pie-chart detailing the relevant functional areas (Fig. 2a), while the entire list of significant pathways is presented in Supplementary Table 2. We conducted the subsequent target prioritization steps using the full list of 916 differentially regulated genes.

Fig. 2

a Pathways pie chart in prenatally stressed animals. The pie chart represents the functional relevance of the significant pathways found modulated in the hippocampus of PNS-exposed adult rats. b Pathways pie chart in adults with a history of childhood trauma. The pie chart represents the functional relevance of the significant pathways found modulated in the blood of adult subjects exposed to childhood trauma

MiRNome changes in the hippocampus of PNS rats, and combined analyses with transcriptomic changes

We subsequently investigated whether exposure to PNS causes changes in miRNAs levels, and whether the identified miRNAs could target transcripts identified in the transcriptomics analysis. After the selection of the species (Rattus Norvegicus), we identified a total number of 1218 miRNAs (including pre-miRNAs and mature miRNAs); out of these, 68 miRNAs (47 of them were mature miRNAs) were significantly modulated (all with q-value <0.05) in animals exposed to prenatal stress, with sixty-five downregulated and only 3 up-regulated (see Table 1).

Table 1

List of miRNAs (both mature and non-mature miRNAs) significantly modulated (all with q-value <0.05) in animals exposed to prenatal stress, compared with control animals

	miRNA name	Fold-change
1	rno-miR-3473	–1,7
2	rno-mir-3593	–1,2
3	rno-miR-339–5p	–1,4
4	rno-mir-1839	–1,2
5	rno-mir-3556b	–1,2
6	rno-miR-133b-5p	–1,3
7	rno-miR-1839–5p	–1,2
8	rno-mir-322	–1,3
9	rno-mir-425	–1,2
10	rno-mir-18a	–1,3
11	rno-mir-3085	–1,3
12	rno-mir-329	–1,2
13	rno-miR-370–3p	–1,4
14	rno-miR-181c-5p	–1,3
15	rno-miR-872–3p	–1,5
16	rno-miR-423–3p	–1,6
17	rno-miR-93–3p	–1,6
18	rno-mir-872	–1,3
19	rno-miR-152–3p	–1,4
20	rno-mir-448	–1,3
21	rno-miR-3068–3p	–1,4
22	rno-mir-532	–1,3
23	rno-mir-290	–1,2
24	rno-miR-106b-3p	–1,4
25	rno-mir-27b	–1,2
26	rno-mir-203a	1,3
27	rno-miR-362–5p	–1,5
28	rno-miR-877	–1,4
29	rno-miR-7a-1–3p	–1,5
30	rno-miR-135b-3p	–1,5
31	rno-miR-6328	–1,2
32	rno-miR-540–3p	–1,4
33	rno-mir-376a	–1,2
34	rno-miR-10b-5p	–1,3
35	rno-miR-425–3p	–1,5
36	rno-miR-330–5p	–1,4
37	rno-miR-6326	–1,3
38	rno-mir-764	–1,3
39	rno-miR-200a-5p	–1,2
40	rno-mir-494	–1,2
41	rno-mir-379	–1,3
42	rno-miR-151–3p	–1,2
43	rno-miR-28–3p	–1,5
44	rno-miR-15b-5p	–1,4
45	rno-miR-324–3p	–1,3
46	rno-miR-871–3p	1,2
47	rno-miR-214–3p	–1,4
48	rno-mir-101b	–1,3
49	rno-miR-325–5p	–1,3
50	rno-miR-339–3p	–1,4
51	rno-mir-196c	–1,4
52	rno-miR-6324	–1,6
53	rno-miR-19b-3p	–1,4
54	rno-miR-20b-5p	–1,6
55	rno-miR-666–3p	–1,5
56	rno-miR-351–5p	–1,2
57	rno-miR-99b-3p	–1,3
58	rno-mir-135a	–1,3
59	rno-miR-139–3p	–1,4
60	rno-miR-342–5p	–1,3
61	rno-miR-376b-3p	–1,3
62	rno-miR-450a-5p	1,3
63	rno-miR-140–5p	–1,6
64	rno-miR-493–3p	–1,5
65	rno-miR-532–5p	–1,3
66	rno-miR-124–5p	–1,3
67	rno-miR-1843–5p	–1,3
68	rno-miR-6215	–1,3

As both transcriptomic and miRNome analyses were conducted in the hippocampus of the same animals, we then performed a mRNA-miRNAs combining analysis that allows the identification of a panel of top-hit genes that were both modulated by PNS exposure and targeted by the miRNAs that were modulated by PNS (see Methods). These analyses identified 528 significant genes, presented in Supplementary Table 3. The specific pathways analysis shows that these were genes involved in neurodevelopment (Axonal Guidance, Protein Kinase-A Signaling, Glucocorticoid Receptor Signaling, TGF-β Signaling) and inflammation (STAT3 Pathway, PTEN Signaling, ILK Signaling, IL-8 signaling); the full list of the 42 pathways is presented in Supplementary Table 4.

Blood mRNA transcriptomics of early life trauma in humans

Comparing subjects with and without exposure to early life trauma (see Methods for sociodemographic and clinical information), we identified 250 genes that were differentially modulated (FC > 1.2, FDR q-value <0.05) (see Supplementary Table 5), involved in the modulation of 41 significant pathways. Of note, these were, again, pathways involved in neurodevelopment (Wnt/Ca + pathway, cAMP signaling, CREB signaling) and inflammation (eNOS signaling, chemokine signaling, B cell activation), similar to what we found in the PNS model. All the significant pathways are summarized by the pie-chart detailing the relevant functional areas in Fig. 2b, and listed in Supplementary Table 6.

Integrating data from PNS in animals and from early life trauma in humans

In the next step, we integrated the 528 genes obtained by combined mRNA/miRNA analyses in the hippocampus of rats exposed to PNS with the 250 genes significantly modulated in the blood of adults exposed to early life trauma. We identified 22 common genes that were present in both lists; according to the hypergeometric distribution test, the probability that these genes were not overlapping due to random probability was p = 1.8 × 10–8.

Finally, we validated all of the 22 genes by Real–Time PCR, both in the hippocampus of animals and in the blood of adults exposed to childhood trauma. In line with the microarray results, 16 genes were modulated in the same direction both in the rat hippocampus and in the human blood, and 6 were modulated in the opposite direction. Specifically, 15 were up-regulated in both: Alpha-2-Macroglobulin (A2M), AT-Rich Interaction Domain 5B (ARID5B), Arrestin Domain Containing 4 (ARRDC4), EPH Receptor A4 (EPHA4), F-Box Protein 32 (FBXO32), Forkhead box protein O1 (FoxO1), Heat Shock Transcription Factor 2 (HSF2), Isochorismatase Domain Containing 1 (ISOC1), Low Density Lipoprotein Receptor Adapter Protein 1 (LDLRAP1), Leucine-rich repeat neuronal protein 3 (LRRN3), Myosin ID (MYO1D), Phosphatidylinositol-4-Phosphate 3-Kinase, Catalytic Subunit Type 2 Beta (PIK3C2B), Phosphatidic Acid Phosphatase Type 2A (PPAP2A), Sterile Alpha Motif Domain Containing 12 (SAMD12), and Serine Incorporator 5 (SERINC5); one was down-regulated in both: Transforming Growth Factor, Beta 1 (TGF-β1); one was up-regulated in rats and down-regulated in humans: Solute Carrier Family 24 (Sodium/Potassium/Calcium Exchanger) Member 3 (SLC24A3); and 5 were down-regulated in rats and up-regulated in humans: B-Cell CLL/Lymphoma 2 (BCL2), B-Cell CLL/Lymphoma 9 (BCL9), Lymphoid Enhancer-Binding Factor 1 (LEF1), Lin-54 DREAM MuvB Core Complex Component (LIN54), and Post-GPI Attachment to Proteins 1 (PGAP1). These genes are presented in Supplementary Table 7, where we also present the FCs obtained from the microarray studies and the FCs from Real Time PCR analyses. For all the above-mentioned molecular and biochemical analyses, all data met the assumptions of normal distribution and equality of variance.

Gene network analysis and selection of candidates for GxE interaction

We focused on the 16 genes (A2M, ARID5B, ARRDC4, EPHA4, FBXO32, FoxO1, HSF2, ISOC1, LDLRAP1, LRRN3, MYO1D, PIK3C2B, PPAP2A, SAMD12, SERINC5, and TGF-β1) that were modulated in the same direction both in animals and in humans, and we applied a gene network analysis to identify possible interactions among these genes. Indeed, our main theoretical framework was based on the notion that depression is a complex disorder, where multiple genes interact with each other as belonging to common or overlapping signaling systems. Therefore, also to minimize any chance findings, we decided a priori to select, for the subsequent GxE interaction analyses, only those genes that interacted with each other through physical interaction, co-expression or involvement in common pathways.

Using IPA Software, we observed only one cluster of interacting genes, represented by A2M, FoxO1 and TGF-β1 (Fig. 3). Using a gene enrichment/pathway analysis, we then confirmed that this cluster is involved in cytokines signaling, TGF- β1 signaling and glucocorticoid receptor signaling. As we found that A2M, FoxO1, and TGF-β1 form a single individual cluster of genes interacting with each other, we focused on these three genes for the subsequent GxE studies.

Fig. 3

Network Analyses on the 16 common genes between rats and humans. Network Analyses of the 16 common genes modulated in the same direction both in the hippocampus of rats exposed to prenatal stress and in the blood of adults exposed to childhood trauma. The blue lines indicate direct interactions between molecules (A2M, FoxO1, and TGF-β1). The dot lines indicate the involvement of molecules within pathways

Gene X environment interaction studies

For our main analyses, we focused on all the SNPs located in these three genes and available in both the Grady Trauma Project and Helsinki Birth cohort (FoxO1, n = 132 SNPs; A2M, n = 91 SNPs; and TGF-β1, n = 14 SNPs), and we tested their interaction with emotional stress in childhood in predicting depressive symptoms in adulthood (BDI scores) for both cohorts (Table 2). Then, we assessed the combined effect over both cohorts and performed a meta-analysis (See Supplementary Table 8). The meta-analyses indicated that 6 SNPs, all located within the FoxO1 gene, were significantly associated with emotional stress in predicting depressive symptoms, and also survived multiple testing correction for GxE interactions, with the SNP rs17592371 showing the strongest significance (p = 0.003, Bonferroni correction). For all the 6 SNPs, p-values for heterogeneity of effects were not significant, indicating that fixed-effects meta-analysis was justified. The GxE interaction effect for the SNP rs17592371 in both cohorts is presented in Fig. 4, where we can see that individuals with “at risk” genotypes (CC) developed more depressive symptoms in the presence of early separation as compared to individuals with the “low risk” genotypes (CT and TT).

Table 2

Gene X Environment interactions in the two clinical samples for the six genes identified to be associated with emotional stress in the Grady Trauma Project and in the Helsinki Birth Study Cohort

GENE	SNPs FOXO1	Grady Trauma Project		Helsinki Birth Study Cohort
		β score	p-value	β score	p-value
FoxO1	rs17592371	–3,035	<0,0001	–1,780	0,033
FoxO1	rs2297626	–2,601	0,0002	–1,808	0,037
FoxO1	rs17592468	–2,549	0,0003	–1,763	0,034
FoxO1	rs28553411	–2,549	0,0003	–1,763	0,034
FoxO1	rs7319021	–2,549	0,0003	–1,763	0,034
FoxO1	rs12585452	–2,416	0,0007	–1,796	0,032

Fig. 4

Interaction effect between rs17592371 and stress on BDI symptoms. Estimated marginal means of BDI scores for different genotypes of the FoxO1 single nucleotide polymorphism rs17592371 in subjects who were or were not exposed to emotional abuse in the Grady Trauma Project (on the left) and in the Helsinki Birth Study Cohort (on the right). Values are adjusted for gender and age

We also conducted additional, secondary analyses on these three genes, in the Grady Trauma Project sample, because it has a wider range of abuse phenotypes. Specifically, we tested the interaction between childhood sexual, physical or emotional abuse, and depressive symptoms in adulthood (BDI scores). Within FoxO1, 6 SNPs showed significant GxE interactions with sexual abuse, 40 SNPs with physical abuse and 40 SNPS with emotional abuse (although with a large overlap between SNPs); within A2M, 7 SNPs showed significant GxE interactions with emotional abuse; and within TGF-β1, 4 SNPs showed significant GxE interactions with sexual abuse and 1 SNP with emotional abuse (See Supplementary Table 9).

FoxO1 is activated in hippocampal progenitor stem cells following cortisol exposure, and it mediates the negative effect of stress on neurogenesis

As FoxO1 was up-regulated both in the hippocampus of animals exposed to PNS and in the blood of subjects exposed to childhood trauma, and it was also the most significant gene from our GxE analyses across the two clinical cohorts, we decided to test the functional relevance of this gene. Specifically, we wanted to test: 1) whether the stress-induced FoxO1 upregulation could be replicated in vitro upon treatment with the human stress hormone, cortisol, using our established model of ‘depression in a dish’, the immortalized, multipotent human fetal hippocampal progenitor cell line (HPC0A07/03C); and 2) whether FoxO1 is involved in the cortisol-induced reduction of neuronal proliferation that we had previously described in this model [28].

Consistent with our hypothesis, FoxO1 mRNA levels were significantly increased following 24 h of treatment with 100 μM cortisol (FC = + 3.9, p < 0.0001, Fig. 5a). Moreover, the cell-permeable inhibitor of FoxO1, AS1842856 (100 nM), was able to completely inhibit the cortisol-induced upregulation of FoxO1 (Fig. 5a) and the cortisol-induced reduction in neuronal proliferation (BrdU + cells; one-way ANOVA, p = 0.01 versus cortisol, Fig. 5b, c). AS1842856 alone did not exert any effects on proliferation (p = 0.9), indicating that FoxO1 mediates the effect of stress, but it is not active under baseline conditions.

Fig. 5

Cortisol treatment (100 μM) increases FoxO1 and SGK1 mRNA levels and decreases neurogenesis via a FOXO1-dependent effect. a Gene expression levels of FoxO1 and of SGK1. b BrdU immunocytochemistry was used to assess proliferation. c The FOXO1 inhibitor, AS1842856 (100 nM) counteracted the CORT-induced reduction in proliferation (n = 5). Data are mean ± SEM **p < 0.01; ***p < 0.001

Finally, in order to investigate the molecular pathways by which FoxO1 signaling regulates the effects of stress on neuronal proliferation, we measured also the mRNA levels of a stress- and cortisol-induced gene, the Serum Glucocorticoid Kinase 1 (SGK-1), again following cortisol with or without AS1842856. SGK1 is a stress responsive gene that mediates some of the glucocorticoid effects on brain function [43], including the negative effects of cortisol on neurogenesis [28], and we found that SGK1 is elevated in both the hippocampus of stressed rats and in the blood of depressed patients [44]. We found that 24 h of cortisol (100 uM) treatment up-regulates SGK1 mRNA levels (FC = + 3.7, p < 0.0001, Fig. 5a), replicating again our previous findings [44], and that AS1842856 was able to counteract such changes (Fig. 5a).

Discussion

We provide a novel ‘hypothesis-free’ approach for the identification and prioritization of candidate genes that can be investigated in GWAS datasets to test GxE interaction studies in depression. In particular, by employing “omics” approaches in animals (mRNAs and miRNAs) and in humans (mRNAs), and cross-species validation, we have identified one cluster of genes, comprising FoxO1, A2M, and TGF-β1, that show significant GxE interactions predicting adult depression in the context of early life trauma in our study, and have also been previously involved in the regulation of mechanisms relevant to stress and depression [45-49]. Our main findings, confirmed by a meta-analysis run on both cohorts, indicate that several SNPs within FoxO1 interact with emotional stress during childhood in predicting adult depression, both in the Grady Trauma Project, as well as in the Helsinki Birth Cohort Study. Importantly, we also describe a functional role of FoxO1 in stress responsivity in our in vitro model of human hippocampal progenitor cells. In addition, secondary analyses show that several SNPs within FoxO1, A2M, and TGF-β1 interact with emotional, physical and sexual abuse in predicting adult depression in the Grady Trauma Project.

A large body of evidence has rapidly accumulated over the past years describing interactive effects of genetic factors and early life stressors in determining the risk of depression [8]. These ‘hypothesis-driven’ studies have mostly focused on genetic variations in candidate genes acting in neurobiological systems that have been previously implicated in the pathophysiology of major depression. As mentioned in the Introduction, the most consistent examples of this approach are studies examining the 5-HTTLPR SNP in the serotonin transporter [13, 50–53] and several functional SNPs within FKBP5 [17, 18]. However, numerous studies, and some meta-analyses, have produced discrepant results [50, 54–56]. New gene-environment interactions continue to be reported in psychiatry, including in genes as biologically varied as BDNF [57, 58], CD38 [59], ADCYAP1R1 [60], DRD2 [61] and GRIN2B [62], but always from hypothesis-driven research, which has a high risk of non-replication. Indeed, considering the current emphasis on omics-based approach in genetic-research, and especially in the context of a paucity of significant genetic effects for depression in GWAS studies, other statistical and biological prioritization strategies need to be developed to facilitate a systematic discovery of novel gene-environment interactions in genome-wide analyses. To our knowledge, the present study is the first report of a hypothesis-free, omics-, cross-tissues and cross-species approach to identify candidate genes for GxE interaction analyses. Interestingly, Mamdani and colleagues (2015) [63] recently used a similar approach to detect genes of interest in the context of alcohol dependence, combining data from mRNA and miRNA expression patterns in the brain of 18 patients with alcohol dependence and 18 matched controls; although they did not test for GxE interactions, they identified a network of hub genes that was enriched for genes identified by the Alcohol Dependence Genome Wide Association Study.

Our strongest findings, replicated in both clinical cohorts, pertain to FoxO1. FoxO1 belongs to the forkhead family of transcription factors, which are characterized by a distinct forkhead domain; it is the main target of insulin signaling, it regulates metabolic homeostasis in response to oxidative stress, it is modulated by glucocorticoids and it enhances inflammation [64]. It is also an important regulator of cell death acting downstream of Cyclin-dependent kinase 1 (CDK1), AKT Serine/Threonine Kinase 1 (AKT1) and macrophage stimulating 1 (MST1), and indeed several studies have suggested that activation of FoxO1 contributes to glucocorticoid-induced cell death [65]. Importantly, FoxO1 is also involved in the maintenance of human embryonic stem cells pluripotency, a function that is mediated through direct control by FoxO1 of OCT4 and SOX2 gene expression, through occupation and activation of their respective promoters.

The six FoxO1 SNPs that were found associated with emotional stress in both cohorts are intronic SNPs, and therefore the mechanisms of action remain unclear, although likely affecting gene expression and not the function of the gene product. It is also possible that these SNPs can influence the protein conformation or the accessibility to phosphorylation sites. Indeed, a number of kinases can phosphorylate and regulate FoxO1 proteins either positively or negatively [66-68], including AKT serine/threonine protein kinase downstream of the PI(3)K (phosphatidylinositol-3-OH kinase) signaling pathway [69, 70], as well as SGK1, a kinase that is specifically activated by stress, glucocorticoids and depression. In particular, SGK1 phosphorylates FoxO1 at different sites, leading to subcellular redistribution of FoxO1 from the nucleus to the cytosol [71, 72]. Of note, we have previously described that SGK1 mediates the inhibitory effect of glucocorticoids on neurogenesis, and that its blood mRNA expression levels are elevated in patients with major depression [28]. Moreover, in the present paper we show that the glucocorticoid, cortisol, induces the mRNA levels of both FoxO1 and SGK1, and that an inhibition of FoxO1 prevents the SGK1 upregulation, indicating that this could be dependent on FoxO1, and thus that SGK1 may represent a novel target gene for FoxO1. Interestingly, Kaiser et al., [65] had shown similar findings in mouse pancreatic islet cells, although their paper only identified FoxO1 as a target of SGK1, and not vice versa. Taking into account also this last study, it is possible to speculate that FoxO1 activation during stress stimulates SGK1 activation, which could then further enhance FoxO1 phosphorylation.

There is some limited evidence already for a role of FoxOs in psychiatric disorders. Indeed FoxOs, especially FoxO6, are mainly expressed in hippocampus, amygdala and nucleus accumbens, which are important brain areas involved in aversive and rewarding responses to emotional stimuli [73, 74], suggesting that FoxOs may be involved in the regulation of mood and emotion. Moreover, FoxOs are regulated by antidepressants, as well as serotonin and norepinephrine receptor signaling [75]. Indeed, FoxO3a has emerged as a candidate gene involved in bipolar and unipolar disorders [76]. Moreover, in a learned helplessness animal model, inescapable shocks significantly reduce the phosphorylation of FoxO3a and induces the nuclear location of FoxO3a in the cerebral cortex; moreover, learned helplessness behavior is markedly diminished in FoxO3a-deficient mice [77]. These results are consistent with the report that FoxO3a influences behavioral processes linked to anxiety and depression, and that FoxO3a- deficient mice show reduced depressive behavior [78] while FoxO1 KO mice displayed reduced anxiety [78] behavior. This is in line with our finding in the present paper, as we describe higher levels of FoxO1 in the brain of rats exposed to prenatal stress, in the blood of subjects exposed to childhood trauma, and in cells treated with the stress hormone, cortisol; moreover, in vitro FoxO1 mediates the cortisol-induced reduction of neuronal proliferation. Together, these data indicate that enhanced levels of FoxO1 due to stress exposure are associated with enhanced vulnerability to develop stress related disorders.

The findings that both A2M and TGF-β1 have SNPs with significant GxE interactions in the Grady Trauma Project, albeit limited by the secondary nature of the analyses and the lack of a replication sample, still are of some interest, as both genes have been implicated in biological processes relevant to depression. A2M has traditionally been viewed as an inflammatory fluid proteinase scavenger, and recent studies have demonstrated the ability of A2M to bind to a plethora of cytokines, including TGF-β1 [79]. TGF-β1, in turn, is a member of TGF-β superfamily, which regulates neuronal survival, neurogenesis, synaptogenesis and gliogenesis [80-83]. Similar to the present study, where we found increased A2M mRNA levels following stress both in the rodent hippocampus and in the human blood, A2M is elevated in clinical samples and in experimental models relevant to depression; for example, higher levels of blood A2M are present in patients with total gastrectomy that develop depression, indicating that A2M elevation may be implicated in depression pathogenesis in the context of a pro-inflammatory status [84]; and fishes subjected to the stress of repeated emersions show higher blood levels of A2M [85]. High levels of A2M have indeed been associated with systemic inflammation [86], and therefore it is possible to speculate that the high levels of A2M described in this and other studies may translate into an increased risk of depression by activating the inflammatory system [40, 87, 88]. TGF-β1, in contrast, is reduced in conditions related to stress and depression, in the present study, as well as in previous studies. For example, TGF-β1 levels are reduced in the blood of depressed patients [89-91], and we have previously shown that the TGF-β1-SMAD signaling is one of the pathways that is most significantly down-regulated in neurons by in vitro exposures to high concentrations of cortisol, which mimics depression in vitro and negatively affects neurogenesis [28]. Finally, is of note that several SNPs within A2M and TGF-β1, including those that we have identified through our GxE interaction analyses, are associated with the development of several Complex Diseases and Disorders (see genetic association database from complex diseases and disorders https://geneticassociationdb.nih.gov). Moreover, genetic variants in TGF-β1 have been associated with other brain disorders, such as multiple sclerosis [92] and Alzheimer disease [93, 94].

It is important to emphasize that GxE interactions, while based on the effects of specific changes in DNA sequence, are likely to involve epigenetic mechanisms [95-97]. For example, recent studies have shown that the above-mentioned functional polymorphism in the FKBP5 gene increases the risk of developing depression and PTSD in adulthood by allele-specific, childhood trauma–dependent DNA demethylation in functional glucocorticoid response elements of FKBP5, followed by a long-term dysregulation of the stress hormone system and a global effect on the function of immune cells and of brain areas associated with stress regulation [98]. Besides DNA methylation, miRNAs have also recently emerged as important in the long-term regulation of gene expression associated with stress early in life [99-103]. The ability of miRNAs to selectively and reversibly silence the mRNAs of their target genes [104, 105], together with their involvement in biological processes modulated by stressful life events [106], make miRNAs well-suited to serve as fine regulators of the complex and extensive molecular network involved in stress response; hence, our decision in this paper of merging transcriptomics and miRNomics data in the selection of top-hit genes.

There are two limitations in our study that should be mentioned. First, only interactions with emotional stress were tested in both samples, and six FoxO1 SNPs survived multiple-testing correction and were confirmed by meta-analysis; the other SNPs in FoxO1, A2M and TGF-β1 that were associated with multiple types of abuse were only tested in the Grady Trauma sample, in secondary analyses. We trust that the use of a different clinical population (the South-East London sample) to generate transcriptomics, together with the cross-species validation, support the relevance of all of our findings, including the evidence for A2M and TGF-β1, but clearly these secondary analyses are only exploratory at this stage. Second, the samples size of the cohort of adult subjects with a history of childhood trauma used for transcriptomic analyses is relatively small, with n = 20 subjects per group. However, to the best of our knowledge, this is the first gene expression microarray study comparing the transcriptome in subjects exposed to childhood trauma versus non-exposed individuals, and thus provides informative findings. Moreover, the validity of our transcriptomics results is strengthened by the fact that we compare subjects not exposed to early life trauma with subjects exposed to the most severe types of trauma or to more than one traumatic event, thus maximizing the biological differences between the two groups. Finally, only the genes regulated also in the animal brains were selected for the further genetic association analyses, again strengthening the validity of our findings. All together, these lines of evidence defend the appropriateness of our sample size for the transcriptomics analyses.

In conclusion, our data provide a novel approach for the selection of novel susceptibility genes for GxE interaction in depression, alternative to ‘hypothesis-driven research’ and resulting in the prioritization of candidate genes starting from a ‘hypothesis-free’ approach. The evidence that the three identified genes all interact with environmental stress in inducing depressive symptoms, and that the statistically strongest gene, FoxO1, has functional relevance in stress-induced reduction of neurogenesis, support our theoretical framework. We propose that our prioritization strategy will limit the number of investigated genes in GWAS-based GxE interaction studies and thus will enhance the statistical power to identify significant findings. Moreover, we propose that FoxO1 may represent a target for the development of novel pharmacological therapies to prevent or treat stress-induced mental disorders.

Supplementary Table 1

Supplementary Table 2

Supplementary Table 3

Supplementary Table 4

Supplementary Table 5

Supplementary Table 6

Supplementary Table 7

Supplementary Table 8

Supplementary Table 9

Supplementary Material