Association Between FoxO1, A2M, and TGF-β1, Environmental Factors, and Major Depressive Disorder

Introduction: Investigations of gene-environment (G×E) interactions in major depressive disorder (MDD) have been limited to hypothesis testing of candidate genes while poly-gene-environmental causation has not been adequately address. To this end, the present study analyzed the association between three candidate genes, two environmental factors, and MDD using a hypothesis-free testing approach.

Methods: A logistic regression model was used to analyze interaction effects; a hierarchical regression model was used to evaluate the effects of different genotypes and the dose-response effects of the environment; genetic risk score (GRS) was used to estimate the cumulative contribution of genetic factors to MDD; and protein-protein interaction (PPI) analyses were carried out to evaluate the relationship between candidate genes and top MDD susceptibility genes.

Results: Allelic association analyses revealed significant effects of the interaction between the candidate genes Forkhead box (Fox)O1, α2-macroglobulin (A2M), and transforming growth factor (TGF)-β1 genes and the environment on MDD. Gene-gene (G×G) and gene-gene-environment (G×G×E) interactions in MDD were also included in the model. Hierarchical regression analysis showed that the effect of environmental factors on MDD was greater in homozygous than in heterozygous mutant genotypes of the FoxO1 and TGF-β1 genes; a dose-response effect between environment and MDD on genotypes was also included in this model. Haplotype analyses revealed significant global and individual effects of haplotypes on MDD in the whole sample as well as in subgroups. There was a significant association between GRS and MDD (P = 0.029) and a GRS and environment interaction effect on MDD (P = 0.009). Candidate and top susceptibility genes were connected in PPI networks.

Conclusions: FoxO1, A2M, and TGF-β1 interact with environmental factors and with each other in MDD. Multi-factorial G×E interactions may be responsible for a higher explained variance and may be associated with causal factors and mechanisms that could inform new diagnosis and therapeutic strategies, which can contribute to the personalized medicine of MDD.

Introduction

Major depressive disorder (MDD) is the most common psychiatric disorder and is associated with high morbidity, mortality, costs, and risk of suicide (1–3). The heritability of MDD is approximately 30%–40% (4) and may be greater for recurrent, early-onset disease (5). Establishing the contribution of genetics to MDD can lead to more accurate and timely diagnosis and clinical management, which has significant public health implications.

Studies on main genetic effects in MDD have reported conflicting findings (6). One reason for this is that risk of MDD is polygenic, with many susceptibility genes exerting small effects (7). Additionally, genotypic variations among individuals may increase the risk of MDD only upon exposure to adverse environmental factors, a phenomenon known as gene-environment (G×E) interaction (8). The association between candidate susceptibility genes and disease occurrence has mainly been predicated on the observation that genetic variants not only increase MDD risk but can also explain the underlying biological and molecular mechanisms (9). However, such hypothesis-driven approaches to the study of MDD etiology can only detect a fraction of genetic variants (8). On the other hand, a hypothesis-free approach requires enormous datasets that include both exposure [e.g., stressful life events (10, 11) and childhood adversities (9)] and phenotype (MDD) data. To this end, only three genome-wide environmental interaction studies has been carried out (12–14). In a recent report, a novel omics-based approach was used to identify genetic variants for genome-wide G×E interaction studies based on cross-species and -tissue biological prioritization strategies. However, genomics approaches usually cannot identify causal variants or genes (15). Integrating multiple omics approaches can provide a more comprehensive view of the biology of MDD.

Three candidate MDD susceptibility genes [Forkhead box (Fox)O1, transforming growth factor (TGF)-β1, and α2 macroglobulin (A2M)] that are ideal candidates for G×E interaction analyses were previously identified using an omics-based approach. It has been suggested that environmental and genetic factors contribute equally to the development of mental illness (16). Distal environmental risk factors such as family history (FH) and culture are important because they increase vulnerability to proximal factors; (17) however, the latter are more relevant for G×E investigations since they are more likely to meet the criteria for risk factors and lend themselves to biologically plausible hypotheses regarding their effects on specific neural systems that underlie psychopathological symptoms.

Genome-wide association studies (GWAS) have identified more than 100 genetic variants associated with MDD (18). Each of these has a small effect and the number of associated genetic markers increases with sample size. However, the utility of individual genetic markers for MDD risk prediction is uncertain (19). A multilocus genetic risk score (GRS)-based analysis has been proposed for combining the relatively small effects of single genes to better assess the complex relationship between genetic markers and MDD (20). GRS analyses have shown that including more weakly associated genetic variants improves the prediction of mental illness risk (21).

Proteins interact to mediate many cellular processes; disrupting one subunit of a protein complex can lead to direct and indirect functional consequences in various diseases and physiological conditions (22). For a specific illness, interactions between proteins encoded by susceptibility genes tend to be more frequent than those between random proteins, as revealed by protein-protein interaction (PPI) network analyses (22).

In the present study, we investigated the association between FoxO1, TGF-β1, and A2M genes and MDD in a large population. G×E, G×G, G×G×E, and environment × environment (E×E) interactions were analyzed in the context of MDD. We also assessed haplotypes of these three genes and the relationship between haplotype and environment. We used a GRS approach to calculate the combined effects of these three genes on the prediction of MDD and the GRS–environment interaction (GRS×E) in MDD. Finally, we examined FoxO1, TGF-β1, and A2M PPI networks to clarify their interactions with each other and with other MDD-related proteins.

Materials and Methods

Study Population

From November 2014 to December 2018, 800 MDD patients (564 women and 236 men) were recruited for the study (mean age: 45.64 ± 14.10 years) along with 800 age- and sex-matched control subjects from the same geographic area in Northern China. All subjects were of Chinese Han ethnicity and provided written, informed consent before participating in the study. The study protocol was approved by the ethics committee of Harbin Medical University.

Independent Measures

Participants completed three questionnaires: a socio-demographic questionnaire, the Chinese version of the 24-item Hamilton Rating Scale for Depression (HRSD-24), and the Life Events Scale (LES). The self-rating socio-demographic questionnaire was adapted from the version developed by the Epidemiology Department of Harbin Medical University and was used to collect detailed background information on family psychiatric history, childhood trauma history, and socioeconomic background. The HRSD-24 (23) is widely used to measure depression symptoms; (24–26) a number of patients above the threshold (21 points) were included in the study. The LES questionnaire for measuring negative life events contained 48 items in three dimensions: family life (28 items), work-related problems (13 items), and social and other aspects (seven items). The LES was scored based on the occurrence/absence (1 and 1, respectively) and frequency (0 or more) of SLEs (26, 27).

Genotyping

Single nucleotide polymorphisms (SNPs) in the FoxO1 (rs2297626, rs7319021, rs28553411, rs17592468, and rs17592371), A2M (rs10492115, rs226415, rs10842849, rs11048839, rs10842847, and rs669), and TGF-β1 (rs2317130, rs1800469, rs12983775, rs12462166, and rs2241715) genes were selected for genotyping, the candidate SNPs comprised most of the allelic variants with r² > 0.8 in the Asian population. Genomic DNA was extracted from venous blood samples using the AxyPrep Blood Genomic DNA Miniprep kit (Axygen, Union City, CA, USA). Polymerase chain reaction (PCR) was performed in a reaction volume of 2 μl DNA, 7.5 μl 2× PCR mix, 2 μl primer mix, 0.2 μl ExoI enzyme (Fermentas, Burlington, ON, Canada), 0.7 μl ExoI buffer (Fermentas), and 0.8 μl FastAP enzyme (Fermentas). Amplification conditions were as follows: at 95°C for 3 min; 35 cycles of 94°C for 15 s, 55°C for 15 s, and 72°C for 30 s; and 72°C for 3 min. PCR products were purified by incubation with ExoI and FastAP at 37°C for 15 min and at 80°C for 15 min. Extension conditions were as follows: 96°C for 1 min, and 30 cycles of 96°C for 10 s, 52°C for 5 s, and 60°C for 30 s. Extension products after denaturation at 95°C for 3 min were analyzed by DNA sequencing (Applied Biosystems, Foster City, CA, USA).

Statistical Analysis

Haploview v.4.0 software was used to assess Hardy-Weinberg Equilibrium (HWE) and pairwise linkage disequilibrium (LD) and calculate minimal allele frequency (MAF) for genotyped polymorphisms (28). The χ² test/Fisher’s exact test was used to analyze differences in the distributions of independent variables. The Bonferroni method was adopted for multiple-testing correction. Associations between phenotype and independent variables were analyzed by multivariable logistic regression, with polymorphisms scored as 0, 1, or 2 depending on the carrier status of the minor allele, sex, family history (FH), occurrence/absence or number of SLEs and CAs, and interaction effects (E×E/G×E/G×G/G×G×E/GRS×E). Allele-carrier status and dose-response effect of environmental factors were also determined. Hierarchical regression analysis was performed using RStudio v.1.1.423 software. Alpha levels were corrected with the number of polymorphisms of candidate genes (five polymorphisms for FoxO1 and TGF-β1; six polymorphisms for A2M), so that the p values correction was 0.05/5 = 0.01 for FoxO1 and TGF-β1, 0.05/6 = 0.0083 for A2M. Study power was calculated with QUANTO 1.2.4 (http://hydra.usc.edu/gxe/).

Haplotype analysis was performed using UNPHASED v.3.0.11 software (29). Maximum likelihood haplotype frequencies in the study population were assessed with the expectation maximization algorithm. Rare haplotypes with a frequency < 1% were excluded from analyses. The global association and effect of individual haplotypes on phenotype were analyzed in the total sample and in four subgroups classified according to positive/negative environmental factors (SLEs/CAs). The significance level of the global association was determined with the likelihood ratio test. Individual effects of each haplotype, i.e., the difference in effect between a haplotype and all others pooled together, were computed with the score test. A permutation analysis was conducted to assess the reliability of the results; 1,000 random permutations were set to generate empirical P values. Permuted P < 0.05 was considered significant in the haplotype analyses.

GRS Analysis

To estimate the cumulative contribution of genetic factors to MDD in an individual by taking into account the reported risk alleles, we used the predictABEL package in RStudio software to compute weighted GRS with the following formula:

$$\text{GRS}={\Sigma }_{i=1}^k{\beta }_i{N}_i$$

where GRS is the sum of the effect estimates, k is the number of independent genetic variants with strong association as risk predictors, β_i is the weighted coefficient from logistic regression analysis, and N_i is the number of risk alleles for each locus. The association between GRS and MDD was evaluated by logistic regression analysis of FH, sex, GRS, SLEs/CAs, and the interaction between GRS and SLEs/CAs.

Population Stratification Analysis

Population stratification analysis was performed to eliminate the possibility of false-positive associations using STRUCTURE v.2.3.4 software (http://pritch.bsd.uchicago.edu/structure.html). FoxO1, A2M, and TGF-β1 genotype data were used for population stratification analysis of the third-stage sample set. Sixteen SNPs in these genotype datasets were obtained from 311 samples derived from the 1000 Genomes Project, including 99 samples from the Yoruba population of Ibadan, Nigeria (YRI), 109 samples from Chinese Han in Beijing (CHB), and 103 samples from northern and western Europe and the United States (CEU) (1000 Genomes Project Phase 3 at www.internationalgenome.org/data-portal/sample). STRUCTURE assumes that there are K populations in the dataset. The admixture and correlated frequency models had a burn-in length of 10,000 and 10,000 Markov chain Monte Carlo repeats and took into consideration immigration and geography-based genetic isolation. The program was run several times at each K value from 2 to 6 to obtain consistent results.

PPI Analysis

We investigated whether FoxO1, A2M, and TGF-β1 genes interact with each other and are involved in the PPI network containing protein products of the top MDD susceptibility genes identified by a recent GWAS (30–32). The PPI network comprises nodes and edges representing protein and physical interactions, respectively. Proteins encoded by MDD susceptibility genes were used as seed proteins. STRING v.11.0 software (https://string-db.org/cgi/input.pl) (33) was used to reconstruct the PPI network.

Results

Descriptive Statistics

Demographic information on the study population is shown in Table 1. Genetic markers were successfully genotyped at a rate > 95%. No significant deviation from HWE was observed, and MAF was > 5% for each polymorphism (Supplemental Table S1). Pairwise LD D’ values are shown in Figure 1.

TABLE 1

Table 1 Summarized frequencies of demographic and independent variables.

FIGURE 1

Figure 1 Position and linkage disequilibrium (LD) map of genotyped polymorphisms in FoxO1, A2M, and TGF-β1gene. Pairwise LD statistics in examined genes were calculated with Haploview. Squares are bright red if the |D’| value is high, that is, LD is strong.

Association Between Independent Variables and MDD

Two polymorphisms of FoxO1 (rs17592371 and rs2297626), two of A2M (rs669 and rs226415), and one of TGF-β1 (rs12462166) were significantly associated with MDD (Table 2). Family history of MDD (χ² = 54.823, P < 0.0001), occurrence (χ² = 406.993, P <0.0001), and number (χ^{2 =} 781.522, P < 0.0001) of SLEs were related to MDD. The results remained robust after Bonferroni correction (Table 2). The power for association study of FoxO1 (rs17592371 and rs2297626) was 87.43% and 90.45% and the power of TGF-β1 (rs12462166) was 86.55%.

TABLE 2

Table 2 Association of Study Variables with MDD.

Interaction Between SLEs and Candidate Gene Genotypes in MDD

The two FoxO1 (rs17592371 and rs28553411), three A2M (rs10842847, rs10842849, and rs226415), and two TGF-β1 (rs12462166 and rs12983775) gene polymorphisms showed significant interactions with environmental factors (Table 3). In the hierarchical regression analysis, the effect of SLEs on MDD was stronger for the TT genotype (β = 0.612, SE = 0.25, z = 2.450, P = 0.014) than for the TC genotype (β = 0.508, SE = 0.176, z = 2.878, P = 0.003) of FoxO1 rs17592371; and for the GG genotype (β = 0.692, SE = 0.291, z = 2.379, P = 0.017) as compared to the GT genotype (β = 0.542, SE = 0.166, z = 3.258, P = 0.001) of TGF-β1 rs12462166 (Figure 2A). A dose-response effect of the interaction between number of SLEs and genotypes in MDD was only observed for A2M rs10842848, which did not interact with an SLE of 0 (β = 0.099, SE = 0.170, z = 0.586, P = 0.558) but showed interactions with an SLE of 1 (β = 0.420, SE = 0.179, z = 2.339, P = 0.019), 2 (β = 0.647, SE = 0.278, z = 2.330, P = 0.019), and 3 or more (β = 0.749, SE = 0.295, z = 2.535, P = 0.011). That is, high exposure to SLEs was related to an increase in MDD risk (Figure 2B).

TABLE 3

Table 3 Effect of Candidate Genes, Environment, and Their Interaction on MDD.

FIGURE 2

Figure 2 (A) Effect of SLEs on MDD in different genotypes of rs17592371 in FoxO1 and rs12462166 in TGF-β1. (B) Effect of rs10842848 in A2M on MDD in different number of SLEs.

We evaluated the G×G interaction between FoxO1 rs17592371 and rs28553411, between the TGF-β1 rs12462166 and rs12983775, and among A2M rs10842847, rs10842849, and rs226415 as well as the three-way interaction between the two genetic markers that were significant in the G×G interaction analysis and LES (Table 3). The interaction between the two polymorphisms of the FoxO1 gene was significant (β = 0.275, SE = 0.108, z = 2.547, P = 0.010), as was the interaction between those of the A2M gene (β = −0.676, SE = −0.289, z = −2.341, P = 0.019). The effects of three-way interaction (G×G×E) on MDD were also significant for FoxO1 (G×G×SLE: β = 0.211, SE = 0.087, z = 2.409, P = 0.016; G×G×CA: β = 0.321, SE = 0.109, z = 2.922, P = 0.003) and A2M (G×G×CA: β = 1.132, SE = 0.378, z = 2.989, P = 0.002) (Table 3). There results remained significant after controlling for family history and sex in the regression model. The effect of E×E interaction on MDD was significant (β = 0.746, SE = 0.312, z = 2.392, P = 0.016).

Haplotype Analysis

Haplotype analysis was performed for the total sample and subgroups. In the former, 12 haplotypes of five FoxO1 gene polymorphisms, 15 haplotypes of five TGF-β1 gene polymorphisms, and 14 haplotypes of A2M gene polymorphisms had a frequency > 1%. The haplotypes of the FoxO1 gene showed a significant global association (χ² = 19.732, df = 13, P_global < 0.001) and the less common haplotype (G-A-G-C-T-A) of the A2M gene was significantly associated with MDD (χ² = 4.931, P_effect = 0.026). In subgroup samples, haplotypes of the A2M gene within the SLE occurrence group (χ² = 21.267, df = 11, P_global = 0.030) and haplotypes of FoxO1 (χ² = 41.792, df = 11, P_global < 0.001) and TGF-β1 (χ² = 34.531, df = 14, P_global = 0.001) genes within the CA occurrence group were significant in the global association. Two less common haplotypes (G-A-A-C-G-G and G-A-G-C-T-A) of the A2M gene within the SLE occurrence group (χ² = 5.815, P_effect = 0.015; χ² = 5.814, P_effect = 0.015) and a less common A2M haplotype (G-A-G-C-T-A) in the CA occurrence group (χ² = 3.947, P_effect = 0.046) were associated with MDD.

GRS Analysis

After controlling for FH and sex, logistic regression revealed that the 16 SNP GRSs comprising risk variants were associated with MDD (β = 0.022, SE = 0.010, z = 2.180, P = 0.029). Interaction analyses indicated that GRSs interacted with SLEs in MDD (β = 0.054, SE = 0.021, z = 2.577, P = 0.009).

Population Stratification Analysis

The combined population of YRI, CHB, and CEU exhibited obvious stratification (Supplemental Figure S1A). In the triangle chart with K = 3, each angle represented a potentially independent ancestry and the different colored dots represented individuals in assumed population components that did not cluster with the same color in the triangle, indicating that there was no significant stratification in our sample (Supplemental Figure S1B). Consistent results were obtained with K values ranging from two to six. Thus, population stratification was unlikely to be a confounding factor in this study.

PPI Analysis

We detected interactions among FoxO1, A2M, and TGF-β1 proteins and the protein products of multiple MDD risk genes. The top susceptibility genes except for the calcium voltage-gated channel subunit α1 C gene, identified in GWASs of MDD (30–32) formed a densely interconnected PPI network (Figure 3) that included FoxO1, A2M, and TGF-β1. Interestingly, A2M was found to directly interact with FoxO1 and TGF-β1. TGF-β1 directly interacted with FOS, JUN, cyclic AMP response element-binding protein (CREB) 1, brain-derived neurotrophic factor, AKT1, and Mothers against decapentaplegic homolog (SMAD) 2; and FoxO1 directly interacted with FOS, JUN, CREB binding protein, SMAD2, AKT1, AKT2, and sirtuin 1. These factors have been implicated in MDD genetic risk studies (34–38).

FIGURE 3

Figure 3 Proteins encoded by FoxO1, A2M, and TGF-β1 in a densely interconnected protein-protein interaction (PPI) network formed by susceptibility genes of MDD. (33) The amaranthine lines indicate known PPIs that have been experimentally determined, and the light blue lines indicate known PPIs from curated databases. The green lines indicate predicted interactions from gene neighborhood, the dark blue lines indicate predicted interactions through gene co-occurrence and the black lines indicate known PPIs from co-expression. The light green lines indicate predicted PPIs through text mining.

Discussion

The interaction between genetic markers and environment determines the risk of MDD (8). However, most hypothesis-driven studies of MDD to date have not produced reproducible findings. In this study, we investigated candidate MDD risk genes using a hypothesis-free, omics-based, cross-tissue and -species approach. We found that FoxO1 rs17592371 and rs2297626, A2M rs669 and rs226415, and TGF-β1 rs12462166 alleles were associated with MDD. We also found that FoxO1 rs17592371 and rs28553411, A2M rs10842847, rs10842849, and rs226415, and TGF-β1 rs12462166 and rs12983775 interacted with the environment in MDD, suggesting a potential relationship between these genes and the etiology of MDD, which could be a bona fide association or the result of linkage. A known allelic association between a gene and a disorder may reflect a G×E interaction, but the absence of such an association does not disqualify a gene as a candidate for disease risk. Some evidences showed that FoxO1 could express in hippocampus and corpus striatum, FoxO1 may contribute to the pathological process of MDD or other psychiatric disorders (18, 39). It was suggested that immune system associated with MDD and A2M played an important role in the immune system, which indicated that A2M may associate with MDD. Some studies found that the expression of A2M was higher in MDD patients comparing healthy people (21, 40). TGF-β1 can restrain autoimmune response and maintain immune system stability, some studies found that TGF-β1 could moderate the imbalance of proinflammatory cytokines and anti-inflammatory cytokines in MDD patients (15, 16). All of these evidences above were consistent with our findings.

The carrier status of the minor allele and number of SLEs were separately used to test different genotypic effects on MDD and the dose-response relationship between SLEs and MDD. The effect of SLEs on MDD was more highly significant for the TT than for the TC genotype of FoxO1 rs17592371 and for the GG as compared to the GT genotype of TGF-β1 rs12462166. These results indicated that mutant homozygous genotypes may be more significant than the heterozygous genotype for the effect of SLEs on MDD, as previously suggested (41, 42). We also found that exposure to SLEs was related to an increased risk for MDD, which is consistent with previous reports (42, 43), indicating that although the effect of a single environmental factor may be quite small, the cumulative effect of multiple factors may be large and that strong effects typically result from a chain of related events rather than a single factor.

Because gene variants can interact, we examined the SNPs that showed significant G×E interactions to investigate G×G interactions. We found that two polymorphisms in the FoxO1 gene and the A2M gene showed significant interactions and that the G×G×E interaction was also significant. This was similar to a previous finding that serotonin transporter-linked polymorphic region (5-HTTLPR) gene and SNP rs140700 both interacted with the environment and with each other in MDD, and that the 5-HTTLPR×rs140700×E interaction was significant (41). Much of the genetic variation associated with a complex trait can be explained by the joint contribution of multiple genetic markers as well as environmental influence (44). G×G and G×G×E interaction analyses are considered as biologically relevant and explain the role of heritability (45–48). We also found that the E×E (CA × SLE) interaction in MDD was significant. CA increases the likelihood of SLE occurrence (49) that is, early negative experiences increase vulnerability to such events later in life, resulting in a sequential E×E interaction. For instance, mistreatment in childhood caused sensitization to the effects of specific types of SLE in adulthood (49). The present study is the first report of an E×E interaction in MDD.

Although direct associations between haplotype and psychiatric disorders have been reported (50, 51), there have been no studies demonstrating an interaction between haplotype and environmental factors. In our study, the global effect of haplotype on MDD was significant for the FoxO1 and TGF-β1 genes in the total sample and in subjects that had experienced CA. A significant individual effect and global effect of haplotype on MDD was found for the A2M gene in both the total and subgroup samples. The global effect of haplotype on MDD was always significant in combination with environmental factors. Interactions existed not only between the environment and genes, but also between the environment and haplotype in MDD.

GRS analyses have shown that the prediction of mental illness is improved by including more weakly associated genetic variants, suggesting that these influence the risk of mental disorder (21, 40). Our results showed that the association between 16 SNP GRS and MDD was significant, with a significant interaction between GRS and SLEs. Some of the 16 SNPs did not show a significant conditional or interaction effect. GRSs may detect the effects of weaker associated genetic variants (52). The results of the PPI analysis showed that FoxO1, A2M, and TGF-β1 were connected, suggesting that a wide variety of cellular processes are involved in MDD. This could explain the reason for the additive effects of FoxO1, A2M, and TGF-β1 gene variants indicated by the GRSs. In addition, several top MDD susceptibility genes were included in this PPI network, although the clinical significance of these associations requires further investigation. The disturbances of certain cellular processes or pathways contributing to the risk of MDD has been emerging and gained evidences supporting (22, 53, 54). Constructing highly interconnected PPI networks between FoxO1, A2M, and TGF-β1 protein and proteins of multiple defined risk genes for MDD may reveal the underlying biological mechanisms. FoxO1, A2M, and TGF-β1 protein could participate in this PPI network, indicating their potential involvement in the common molecular network modulating the pathogenesis of MDD. It is noteworthy that A2M protein directly interacted with TGF-β1 protein, indicating their connection in the biological process relevant to MDD. A recent study found that A2M as an inflammatory fluid proteinase scavenger could bind to a plethora of cytokines, including TGF-β1. (55)

We investigated the potential influence of population stratification in our samples using STRUCTURE software v.2.3.4 and the third-stage sample set, but did not find any evidence of stratification, which confirmed that the results were not affected by this confounding factor.

Our study had three major limitations. First, we did not use standard questionnaires that have confirmed reliability and validity to measure CA. Second, we did not analyze additional SNPs in our samples and the third-stage sample set, which can improve the reliability of population stratification analysis. Lastly, including a gene expression microarray analysis in our study would have allowed us to compare the transcriptomes of subjects exposed to SLEs/CAs and non-exposed individuals, which could provide more information for predicting MDD risk.

In conclusion, our data indicated that the hypothesis-free approach is useful for identifying novel genes contributing to the G×E interaction in MDD. FoxO1, A2M, and TGF-β1 genes not only interact with environmental factors but also associate with multiple other genes in the MDD G×E interaction. FoxO1, A2M, and TGF-β1 genes may serve the clinical diagnose and treatment of MDD by providing biomarkers, which can contribute to personalized medicine. Future studies will need to focus on the complexity of poly-gene–environmental causation to obtain more detailed insight into the etiology of MDD.

Data Availability Statement

All datasets for this study are included in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by ethics committee of Harbin Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

MZ and LC conducted the statistical analyses and wrote the first draft of the manuscript. ZQ and JZ provided expertise in MDD search. TZ, WZ, and SK did the experiment. XZ, EZ, and XS collected the data. XQ and HP revised the manuscript. YY and XY designed this study and provided expertise. All authors were involved in modifying the secondary-analysis design and editing the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China (81773536) to YY, the Fundamental Research Funds for State Universities of Heilongjiang Province (2017JCZX23) to EZ, and China Postdoctoral Science Foundation Grant (2019M651311) to JZ.

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.

The reviewer JM declared a past co-authorship with several of the authors YY, ZQ, EZ, XS, XY to the handling Editor.

Acknowledgments

We want to express our gratitude to all patients and healthy controls parting in the study, as well as to the psychiatrists for their help in the recruitment and identification of patients with major depressive disorder.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpsyt.2020.00675/full#supplementary-material

References

1. Wittchen HU, Jacobi F, Rehm J, Gustavsson A, Svensson M, Jonsson B, et al. The size and burden of mental disorders and other disorders of the brain in Europe 2010. Eur Neuropsychopharmacol (2011) 21(9):655–79. doi: 10.1016/j.euroneuro.2011.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ferrari AJ, Charlson FJ, Norman RE, Patten SB, Freedman G, Murray CJ, et al. Burden of depressive disorders by country, sex, age, and year: findings from the global burden of disease study 2010. PloS Med (2013) 10(11):e1001547. doi: 10.1371/journal.pmed.1001547

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Murray CJ, Vos T, Lozano R, Naghavi M, Flaxman AD, Michaud C, et al. Disability-adjusted life years (DALYs) for 291 diseases and injuries in 21 regions, 1990-2010: a systematic analysis for the Global Burden of Disease Study 2010. Lancet (2012) 380(9859):2197–223. doi: 10.1016/S0140-6736(12)61690-0

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Sullivan PF, Neale MC, Kendler KS. Genetic epidemiology of major depression: review and meta-analysis. Am J Psychiatry (2000) 157(10):1552–62. doi: 10.1176/appi.ajp.157.10.1552

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Viktorin A, Meltzer-Brody S, Kuja-Halkola R, Sullivan PF, Landen M, Lichtenstein P, et al. Heritability of Perinatal Depression and Genetic Overlap With Nonperinatal Depression. Am J Psychiatry (2016) 173(2):158–65. doi: 10.1176/appi.ajp.2015.15010085

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bosker FJ, Hartman CA, Nolte IM, Prins BP, Terpstra P, Posthuma D, et al. Poor replication of candidate genes for major depressive disorder using genome-wide association data. Mol Psychiatry (2011) 16(5):516–32. doi: 10.1038/mp.2010.38

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hyman S. Mental health: depression needs large human-genetics studies. Nature (2014) 515(7526):189–91. doi: 10.1038/515189a

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Heim C, Binder EB. Current research trends in early life stress and depression: review of human studies on sensitive periods, gene-environment interactions, and epigenetics. Exp Neurol (2012) 233(1):102–11. doi: 10.1016/j.expneurol.2011.10.032

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Caspi A, Sugden K, Moffitt TE, Taylor A, Craig IW, Harrington H, et al. Influence of life stress on depression: moderation by a polymorphism in the 5-HTT gene. Science (2003) 301(5631):386–9. doi: 10.1126/science.1083968

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kendler KS, Kuhn JW, Vittum J, Prescott CA, Riley B. The interaction of stressful life events and a serotonin transporter polymorphism in the prediction of episodes of major depression: a replication. Arch Gen Psychiatry (2005) 62(5):529–35. doi: 10.1001/archpsyc.62.5.529

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Stephan R, Benjamin MN, Aiden C, James TRW, Kai-How F, Peter AH, et al. Biological insights from 108 schizophrenia-associated genetic loci. Nature (2014) 511(7510):421–7. doi: 10.1038/nature13595

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Otowa T, Kawamura Y, Tsutsumi A, Kawakami N, Kan C, Shimada T, et al. The First Pilot Genome-Wide Gene-Environment Study of Depression in the Japanese Population. PloS One (2016) 11(8):e0160823. doi: 10.1371/journal.pone.0160823

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Dunn EC, Wiste A, Radmanesh F, Almli LM, Gogarten SM, Sofer T, et al. Genome-Wide Association Study (GWAS) And Genome-Wide By Environment Interaction Study (GWEIS) Of Depressive Symptoms In African American And Hispanic/Latina Women. Depress Anxiety (2016) 33(4):265–80. doi: 10.1002/da.22484

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Cattaneo A, Cattane N, Malpighi C, Czamara D, Suarez A. FoxO1, A2M, and TGF-beta1: three novel genes predicting depression in gene X environment interactions are identified using cross-species and cross-tissues transcriptomic and miRNomic analyses Mol Psychiatry (2018) 23(11):2192–208. doi: 10.1038/s41380-017-0002-4

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Visscher PM, Brown MA, McCarthy MI, Yang J. Five years of GWAS discovery. Am J Hum Genet (2012) 90(1):7–24. doi: 10.1016/j.ajhg.2011.11.029

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Uher R, Zwicker A. Etiology in psychiatry: embracing the reality of poly-gene-environmental causation of mental illness. World Psychiatry (2017) 16(2):121–9. doi: 10.1002/wps.20436

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Costello EJ, Compton SN, Keeler G, Angold A. Relationships between poverty and psychopathology: a natural experiment. JAMA (2003) 290(15):2023–9. doi: 10.1001/jama.290.15.2023

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hyde CL, Nagle MW, Tian C, Chen X, Paciga SA, Wendland JR. Identification of 15 genetic loci associated with risk of major depression in individuals of European descent. Nat Genet (2016) 48(9):1031–6. doi: 10.1038/ng.3623

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Karg K, Burmeister M, Shedden K, Sen S. The serotonin transporter promoter variant (5-HTTLPR), stress, and depression meta-analysis revisited: evidence of genetic moderation. Arch Gen Psychiatry (2011) 68(5):444–54. doi: 10.1001/archgenpsychiatry.2010.189

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Dudbridge F. Polygenic Epidemiology. Genet Epidemiol (2016) 40(4):268–72. doi: 10.1002/gepi.21966

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wray NR, Lee SH, Mehta D, Vinkhuyzen AA, Dudbridge F, Middeldorp CM. Research review: Polygenic methods and their application to psychiatric traits. J Child Psychol Psychiatry (2014) 55(10):1068–87. doi: 10.1111/jcpp.12295

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Luo X, Huang L, Jia P, Li M, Su B, Zhao Z, et al. Protein-protein interaction and pathway analyses of top schizophrenia genes reveal schizophrenia susceptibility genes converge on common molecular networks and enrichment of nucleosome (chromatin) assembly genes in schizophrenia susceptibility loci. Schizophr Bull (2014) 40(1):39–49. doi: 10.1093/schbul/sbt066

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Hamilton M. A rating scale for depression. J Neurol Neurosurg Psychiatry (1960) 23:56–62. doi: 10.1136/jnnp.23.1.56

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Assmann N, Schramm E, Kriston L. Moderating effect of comorbid anxiety disorders on treatment outcome in a randomized controlled psychotherapy trial in early-onset persistently depressed outpatients. Depres Anxiety (2018) 35(10):1001–8. doi: 10.1002/da.22839

CrossRef Full Text | Google Scholar

25. Erkens N, Schramm E, Kriston L, Hautzinger M, Harter M, Schweiger U, et al. Association of comorbid personality disorders with clinical characteristics and outcome in a randomized controlled trial comparing two psychotherapies for early-onset persistent depressive disorder. J Affect Disord (2018) 229:262–8. doi: 10.1016/j.jad.2017.12.091

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Ma J, Wang L, Yang Y, Qiao Z, Fang D, Qiu X, et al. GNB3 and CREB1 gene polymorphisms combined with negative life events increase susceptibility to major depression in a Chinese Han population. PloS One (2017) 12(2):e0170994. doi: 10.1371/journal.pone.0170994

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Han D, Qiao Z, Chen L, Qiu X, Fang D, Yang X, et al. Interactions between the vascular endothelial growth factor gene polymorphism and life events in susceptibility to major depressive disorder in a Chinese population. J Affect Disord (2017) 217:295–8. doi: 10.1016/j.jad.2017.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics (2005) 21(2):263–5. doi: 10.1093/bioinformatics/bth457

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Dudbridge F. Pedigree disequilibrium tests for multilocus haplotypes. Genet Epidemiol (2003) 25(2):115–21. doi: 10.1002/gepi.10252

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Chang LC, Jamain S, Lin CW, Rujescu D, Tseng GC, Sibille E. A conserved BDNF, glutamate- and GABA-enriched gene module related to human depression identified by coexpression meta-analysis and DNA variant genome-wide association studies. PloS One (2014) 9(3):e90980. doi: 10.1371/journal.pone.0090980

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Xiao X, Zhang C, Grigoroiu-Serbanescu M, Wang L, Li L, Zhou D, et al. The cAMP responsive element-binding (CREB)-1 gene increases risk of major psychiatric disorders. Mol Psychiatry (2018) 23(9):1–11. doi: 10.1038/mp.2017.243

CrossRef Full Text | Google Scholar

32. Haenisch B, Herms S, Mattheisen M, Steffens M, Breuer R, Strohmaier J, et al. Genome-wide association data provide further support for an association between 5-HTTLPR and major depressive disorder. J Affect Disord (2013) 146(3):438–40. doi: 10.1016/j.jad.2012.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics (2003) 164(4):1567–87.

PubMed Abstract | Google Scholar

34. Burokas A, Arboleya S, Moloney RD, Peterson VL, Murphy K, Clarke G, et al. Targeting the Microbiota-Gut-Brain Axis: Prebiotics Have Anxiolytic and Antidepressant-like Effects and Reverse the Impact of Chronic Stress in Mice. Biol Psychiatry (2017) 82(7):472–87. doi: 10.1016/j.biopsych.2016.12.031

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Xu F, Yang J, Chen J, Wu Q, Gong W, Zhang J, et al. Differential co-expression and regulation analyses reveal different mechanisms underlying major depressive disorder and subsyndromal symptomatic depression. BMC Bioinf (2015) 16:112. doi: 10.1186/s12859-015-0543-y

CrossRef Full Text | Google Scholar

36. Noto C, Ota VK, Santoro ML, Gouvea ES, Silva PN, Spindola LM, et al. Depression, Cytokine, and Cytokine by Treatment Interactions Modulate Gene Expression in Antipsychotic Naive First Episode Psychosis. Mol Neurobiol (2016) 53(8):5701–9. doi: 10.1007/s12035-015-9489-3

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Dow AL, Russell DS, Duman RS. Regulation of activin mRNA and Smad2 phosphorylation by antidepressant treatment in the rat brain: effects in behavioral models. J Neurosci (2005) 25(20):4908–16. doi: 10.1523/JNEUROSCI.5155-04.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Kim HD, Hesterman J. SIRT1 Mediates Depression-Like Behaviors in the Nucleus Accumbens. J Neuronsci (2016) 36(32):8441–52. doi: 10.1523/JNEUROSCI.0212-16.2016

CrossRef Full Text | Google Scholar

39. Howard DA-O, Adams MA-O, Shirali MA-O, Clarke TK, Marioni RE, Davies G, et al. Genome-wide association study of depression phenotypes in UK Biobank identifies variants in excitatory synaptic pathways. Nat Commun 9(1):1470. doi: 10.1038/s41467-018-03819-3

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Dudbridge F. Power and predictive accuracy of polygenic risk scores. PloS Genet (2013) 9(3):e1003348. doi: 10.1371/journal.pgen.1003348

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lazary J, Lazary A, Gonda X, Benko A, Molnar E, Juhasz G, et al. New evidence for the association of the serotonin transporter gene (SLC6A4) haplotypes, threatening life events, and depressive phenotype. Biol Psychiatry (2008) 64(6):498–504. doi: 10.1016/j.biopsych.2008.03.030

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Evans GW. The environment of childhood poverty. Am Psychol (2004) 59(2):77–92. doi: 10.1037/0003-066X.59.2.77

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Sameroff AJ, Seifer R, Bartko WT. Environmental perspectives on adaptation during childhood and adolescence. In Luthar S. S., Burack J. A., Cicchetti D., Weisz J. R. (Eds.), Developmental Psychopathology: Perspectives on Adjustment, Risk, and Disorder. (Cambridge England; New York, NY, USA: Cambridge University Press) (1997). pp. 507–526.

Google Scholar

44. Li S, Cui Y. Gene-centric gene-gene interaction: A model-based kernel machine method. Ann Appl Stat (2012) 6(2012):1134–61. doi: 10.1214/12-AOAS545

CrossRef Full Text | Google Scholar

45. Lunzer M, Golding GB, Dean AM. Pervasive cryptic epistasis in molecular evolution. PloS Genet (2010) 6(10):e1001162. doi: 10.1371/journal.pgen.1001162

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Huang W, Richards S, Carbone MA, Zhu D, Anholt RR, Ayroles JF, et al. Epistasis dominates the genetic architecture of Drosophila quantitative traits. Proc Natl Acad Sci U.S.A. (2012) 109(39):15553–9. doi: 10.1073/pnas.1213423109

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zuk O, Hechter E, Sunyaev SR, Lander ES. The mystery of missing heritability: Genetic interactions create phantom heritability. Proc Natl Acad Sci U S A. (2012) 109(4):1193–8. doi: 10.1073/pnas.1119675109

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Hemani G, Shakhbazov K, Westra HJ, Esko T, Henders AK, McRae AF, et al. Detection and replication of epistasis influencing transcription in humans. Nature (2014) 508(7495):249–53. doi: 10.1038/nature13005

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Starr LR, Hammen C, Conway CC, Raposa E, Brennan PA. Sensitizing effect of early adversity on depressive reactions to later proximal stress: Moderation by polymorphisms in serotonin transporter and corticotropin releasing hormone receptor genes in a 20-year longitudinal study. Dev Psychopathol (2014) 26(4 Pt 2):1241–54. doi: 10.1017/S0954579414000996

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Brookes KJ, Mill J, Guindalini C, Curran S, Xu X, Knight J, et al. A common haplotype of the dopamine transporter gene associated with attention-deficit/hyperactivity disorder and interacting with maternal use of alcohol during pregnancy. Arch Gen Psychiatry (2006) 63(1):74–81. doi: 10.1001/archpsyc.63.1.74

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Todd RD, Neuman RJ. Gene-environment interactions in the development of combined type ADHD: evidence for a synapse-based model. Am J Med Genet B Neuropsychiatr Genet (2007) 144b(8):971–5. doi: 10.1002/ajmg.b.30640

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Demirkan A, Penninx BW, Hek K, Wray NR, Amin N, Aulchenko YS, et al. Genetic risk profiles for depression and anxiety in adult and elderly cohorts. Mol Psychiatry (2011) 16(7):773–83. doi: 10.1038/mp.2010.65

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Gilman SR, Chang J, Xu B, Bawa TS, Gogos JA, Karayiorgou M, et al. Diverse types of genetic variation converge on functional gene networks involved in schizophrenia. Nat Neurosci (2012) 15(12):1723–8. doi: 10.1038/nn.3261

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Network, Pathway Analysis Subgroup of Psychiatric Genomics C. Psychiatric genome-wide association study analyses implicate neuronal, immune and histone pathways. Nat Neurosci (2015) 18(2):199–209. doi: 10.1038/nn.3922

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Hope C, Mettenburg J, Gonias SL, DeKosky ST, Kamboh MI, Chu CT. Functional analysis of plasma alpha(2)-macroglobulin from Alzheimer’s disease patients with the A2M intronic deletion. Neurobiol Dis (2003) 14(3):504–12. doi: 10.1016/j.nbd.2003.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar