The UKB envirome of depression: from interactions to synergistic effects

Major depressive disorder is a result of the complex interplay between a large number of environmental and genetic factors but the comprehensive analysis of contributing environmental factors is still an open challenge. The primary aim of this work was to create a Bayesian dependency map of environmental factors of depression, including life stress, social and lifestyle factors, using the UK Biobank data to determine direct dependencies and to characterize mediating or interacting effects of other mental health, metabolic or pain conditions. As a complementary approach, we also investigated the non-linear, synergistic multi-factorial risk of the UKB envirome on depression using deep neural network architectures. Our results showed that a surprisingly small number of core factors mediate the effects of the envirome on lifetime depression: neuroticism, current depressive symptoms, parental depression, body fat, while life stress and household income have weak direct effects. Current depressive symptom showed strong or moderate direct relationships with life stress, pain conditions, falls, age, insomnia, weight change, satisfaction, confiding in someone, exercise, sports and Townsend index. In conclusion, the majority of envirome exerts their effects in a dynamic network via transitive, interactive and synergistic relationships explaining why environmental effects may be obscured in studies which consider them individually.


Results
In the following sections, we present the results of multiple analyses performed using specific tools. First, we provide a detailed description of the envirome map based on the variable dependency structure estimated by the BMLA method. Results are shown both as an undirected graph and as edge probabilities (Fig. 1). Second, we introduce the strong relevance measure provided by the BMLA method and compare relevant factors with respect to lifetime depression (reported by the participants and elaborated by trained research nurses) and probable depression diagnosis (derived from the Mental Health Questionnaire with the method described by Smith et al. 26 ). Third, we analyze structural interactions revealed by the BMLA method, and perform parametric analysis using multivariate odds ratios. Fourth, we investigate synergistic effects which differ from structural interactions as all involved variables have an individual effect. Finally, we present results of deep neural networks assessing predictive power. In addition to comparing strongly associated variables in terms of predictive capabilities, we also investigate the predictive power of variable groups. the map of depression envirome. The primary aim of our study was to identify factors that may influence susceptibility to lifetime depression or the severity of its symptoms. This requires multiple analyses that investigate possible relationships qualitatively and quantitatively. We applied the BMLA method for structural analysis, which provides posteriors for arbitrary dependency, relevance and causal patterns based on the multivariate dependency structure learned from data 21,22,24 . Figure 1 presents the network of variable relationships with respect to reported lifetime depression up to its second neighbors (sex and age and their connections are omitted for visibility purposes). An edge between variables X and Y represents that based on Bayesian relevance analysis a dependency relationship exists between X and Y with a probability not less than 0.5. Note that this edge probability takes into account both possible directions of a directed edge assuming an underlying Bayesian network structure representing dependencies.
One of the remarkable features of this graph is that reported lifetime depression has only a few direct relationships, namely current depressive symptoms, neuroticism, parental depression, and body fat percentage. All other factors shown on the graph are in a non-direct relationship with reported lifetime depression, most of them mediated by current depressive symptoms. Table 1 shows the probabilities of various relationship types such as direct and transitive relationships concerning lifetime depression and current depressive symptoms. A direct relationship between two variables X and Y means there are no intermediary variables between them, i.e., they are connected by a directed edge (X → Y or Y → X) in the graph representing dependency relationships. On the other hand, a transitive relationship between X and Y means that there are one or more intermediary variables (Z) in www.nature.com/scientificreports www.nature.com/scientificreports/ between such that there is a path formed between X and Y by directed edges in the dependency graph (X → Z → Y or Y → Z → X).
Life stress and mental health factors. According to our results, the probability that recent negative life events of the past two years (denoted as Life stress) is in a direct relationship with lifetime depression is relatively low (p Dir = 0.2), whereas the probability of a direct relationship with current depressive symptoms is high (p Dir = 1.0). Conversely, the probability of a transitive relationship is relatively higher in the former case (p Trans = 0.4), and low in the latter case (p Trans = 0.0).
The neuroticism personality trait (denoted as Neuroticism) is in a direct relationship with both lifetime depression and current depressive symptoms with high probability (p Dir = 0.999). In addition, it is directly related to parental depression and bipolar disorder. The sleep quality descriptor (denoted as Insomnia) is not related directly to lifetime depression, but it is in a direct relationship with current depressive symptoms (p Dir = 0.999) and neuroticism (p Dir = 0.999) with a high probability.
Social factors. Mental health related social factors such as being able to confide in someone (denoted as Confide), and the sense of satisfaction (Satisfaction) are only transitively connected to lifetime depression (p Trans = 0.2 and 0.6 respectively) as no direct connection was detected. Regarding current depressive symptoms, the probability of a direct relationship with satisfaction is high (p Dir = 1.0), and moderately high in case of Confide (p Dir = 0.6).
Pain descriptors. Variables describing short or long-term presence of pain (i.e. Headache, Back pain, Stomach or abdominal pain, Neck or shoulder pain, Knee pain, and Pain allover) are not in a direct relationship with lifetime depression (p Dir = 0.0). Instead, results indicate transitive relationships with lifetime depression such that most pain related variables are mediated by current depressive symptoms. For example, the probability of a direct relationship between neck or shoulder pain and lifetime depression is negligible (p Dir = 0.0), whereas a transitive relationship is more probable having a moderately high probability (p Trans = 0.629). In contrast to this, the probability of a direct relationship with current depressive symptoms is remarkably high (p Dir = 1.0). In addition, the transitive relationship of headache with lifetime depression (p Trans = 0.8) is not only mediated by current depressive symptoms (p Dir = 1.0), but also mediated by neuroticism with which it has a direct relationship with high probability (p Dir = 0.8).
Dietary change and metabolic descriptors. Variables related to dietary change and metabolism have multiple relations with lifetime depression and current depressive symptoms. Body fat percentage (Body fat), which is highly correlated with the obesity descriptor (Obesity), is in a direct relationship with lifetime depression with Life stress -light brown, (K) Falls -pink, and reported lifetime depression -red. An edge between two nodes represents a direct relationship, and its width is proportional to the Bayesian edge probability which takes into account both possible edge directions assuming an underlying Bayesian network. Edges with a probability lower than 0.5 are omitted.
high probability (p Dir = 0.8). In contrast with body fat percentage, obesity is not in a direct relationship with lifetime depression (p Dir = 0.0), and the probability of direct connection with current depressive symptoms is also low (p Dir = 0.2). On the other hand, a transitive relationship is more probable in both cases (p Trans = 0.6 and 0.8 respectively). Furthermore, weight change is directly connected only to current depressive symptoms (p Dir = 1.0). In addition, the metabolism descriptor (denoted as Metabolic rate) is not in a direct relationship with either lifetime depression or current depressive symptoms (p Dir = 0.0). Similarly to the obesity descriptor, the variable indicating a substantial change in diet (Dietary change) is not in a direct relationship with lifetime depression (p Dir = 0.0), and the probability of a transitive relationship with current depressive symptoms is higher (p Trans = 0.8) than that of a direct (p Dir = 0.2).
Sports and physical activity. Physical activity related variables such as strenuous sports, exercises, walking, vigorous physical activity, moderate physical activity are not in a direct relationship with lifetime depression, instead corresponding transitive relationships are of moderate probability (strenuous sports p Trans = 0.6, exercises p Trans = 0.4, vigorous physical activity p Trans = 0.4, walking p Trans = 0.2) mostly mediated by body fat. Concerning current depressive symptoms, the probability of direct relationships with physical activity descriptors is moderate (strenuous sports p Dir = 0.6, exercises p Dir = 0.6, vigorous physical activity p Dir = 0.4, walking p Dir = 0.2).
Financial background and qualification. Among the investigated socioeconomic status descriptors, the Townsend deprivation index (Townsend) is not connected directly to lifetime depression (p Dir = 0.0), but there is a moderate probability of a transitive relationship (p Trans = 0.6) mediated by current depressive symptoms. The probability of a direct relationship between household income and lifetime depression is relatively low (p Dir = 0.2), although it plays a remarkable role in some of the interactions detailed later. Regarding current depressive symptoms, a transitive relationship with household income is more probable (p Trans = 0.6) than a direct one (p Dir = 0.4). Qualification (that is whether the subject has a college or university degree) is not connected directly to either www.nature.com/scientificreports www.nature.com/scientificreports/ lifetime depression (p Dir = 0.001) or to current depressive symptoms (p Dir = 0.0). However, the probability that there is a transitive relationship between current depressive symptoms and qualification is high (p Dir = 1.0).

Parental illnesses.
Parental depression is the only parental illness descriptor that is directly related to lifetime depression with high probability (p Dir = 0.999), all other descriptors are only transitively related with various degrees of probability. Results also indicate a transitive relationship with current depressive symptoms with high probability for all such descriptors (p Trans > 0.999).
Other factors. According to our results, alcohol intake and tobacco consumption are not in a direct relationship with lifetime depression (p Dir = 0.0), the probability of a transitive relationship is also low in both cases (p Trans = 0.4 and 0.2 respectively). However, both variables are in a transitive relationship with current depressive symptoms with high probability (p Trans = 0.8).
Age and sex are both not directly connected to lifetime depression, although there is a high probability of a transitive relationship. Regarding current depressive symptoms, the probability of a direct relationship with age is high (p Dir = 0.999), while sex is transitively related with high probability (p Trans = 1). Among the investigated childhood descriptors none of the variables are in a direct relationship with lifetime depression, however it should be noted that childhood trauma items were not available for this analysis.

Relevance of environmental factors.
Identifying direct relationships is a major step towards discovering relevant factors, however relevance can be interpreted in multiple ways. Here we utilize the strong relevance concept according to which strongly relevant variables of a selected target variable consist of (1) direct relationships and (2) interaction terms that have a joint effect on the target involving another variable. This requires the analysis of relevant sets of variables with respect to the target variable. In a Bayesian structural approach, strong relevance (or relevance for short) of a variable is quantified by the posterior probability of the occurrence of the variable in possible models as a direct relationship or as an interaction term with respect to the target (see Strong relevance section of methods for details). Table 2 shows posterior probabilities of strong relevance (p Rel ) for relevant variables with respect to lifetime depression using a cutoff value of 0.2 including both direct and interaction type relations. These results indicate that in addition to the previously investigated direct relationships, i.e. current depressive symptoms (p Rel = 1.0), neuroticism trait (p Rel = 0.999), parental depression (p Rel = 1.0), and body fat percentage (p Rel = 0.8), there are several other variables that are relevant with respect to lifetime depression to some extent due to multivariate interactions. In other words, there are interaction terms forming multivariate interaction patterns involving lifetime depression. For example, sex, risk taking, parental Alzheimer's disease and parental bronchitis are such factors that have a moderate probability of being interaction terms, and thus they can be considered as strongly relevant variables to a certain degree (p Rel > 0.3). On the other hand, household income and life stress are in a direct relationship with lifetime depression with a low but non-negligible probability (p Rel > 0.2) and consequently can be considered worthy of further investigation.  Table 2. Probability of strong relevance and dependency types with respect to lifetime depression. Displayed relationship types include direct relations and interaction terms. Associated posterior probabilities reflect the likeliness that a variable is in a given type of relationship with lifetime depression. The probability of strong relevance is the sum of these probabilities.
www.nature.com/scientificreports www.nature.com/scientificreports/ In addition, we investigated relevant relationships with respect to probable depression diagnosis variables (single depressive episode, moderate depression, severe depression) created by Smith et al. 26 and compared it to lifetime depression as a validation (see the Validation section of methods for details). Results indicate that the envirome map of lifetime depression is similar to that of probable depression diagnosis, such that it consists of similar patterns regarding current depressive symptoms, neuroticism, parental depression, and several moderately relevant factors. A notable difference between the two relationship maps is that while lifetime depression is directly connected only to body fat percentage among diet and metabolism related variables, in case of probable depression diagnosis this relationship is partially replaced by connections with obesity and dietary change.

environment-environment interactions. Structural interactions.
In order to analyze interactions first we investigated strongly relevant sets (with respect to lifetime depression) provided by the applied BMLA method (see the methods section for details). These sets can also be called as structurally relevant sets of variables as they are based on the dependency structure of variables. Table 3 presents the top 4 most probable structurally relevant variable sets detailing the components of possible structural interaction patterns.
Generally, structural interaction patterns have at least one term that has an individual main effect (direct relationship) with respect to the target variable, while the other terms typically have minor or negligible effects individually. The key feature of interactions is the multivariate context, in which a particular set of variables have a considerable effect on the target variable. This context is provided either by the variable with the individual main effect or by additional variables. Figure 2 shows possible structural interactions among elements of each relevant variable set using differently colored markers. Variables representing the neuroticism trait, current depressive symptoms, and parental depression are present in all sets as direct relationships with individual main effects. In addition, parental depression plays a central role in several structural interactions by providing context for interaction terms (this assessment requires the analysis of possible dependency structures of variables not discussed here in details). Similarly, body fat percentage (p Dir = 0.800), life stress (p Dir = 0.200) and household income (p Dir = 0.200) are in a direct relationship with lifetime depression and also play roles in interaction patterns. Bipolar disorder on the other hand is present in almost all relevant sets, but only as an interaction term (p Int = 0.999).
According to Table 3 the most probable relevant set consists of risk taking, sex, parental depression, body fat percentage, bipolar disorder, current depressive symptoms and neuroticism. Based on Table 2, variables such as risk taking (p Int = 0.328), sex (p Int = 0.395) and bipolar disorder (p Int = 0.999) are potential interaction terms, whereas parental depression, body fat percentage, neuroticism and current depressive symptoms are potential main effects of these structural interactions. The second sets includes qualification (p Int = 0.200), parental bronchitis (p Int = 0.477), and bipolar disorder as potential interaction terms. The third set contains several interaction terms out of which exercises (p Int = 0. These results (and structural interactions in general) indicate that based on the dependency structure these variables have a multivariate effect on the target. In most cases the individual relevance of these variables is  www.nature.com/scientificreports www.nature.com/scientificreports/ moderate or low with respect to the target, but they have a higher relevance as a pattern. Although structural interactions do not provide information on the parametric nature of these effects, they can be utilized to direct effect size analysis efforts.
Parametric interactions. In most cases when interactions are considered, the parametric aspect of relevance is investigated by applying various effect size measures such as the odds ratio for a binary target. The challenge is that the individual effect of interaction terms tends to be moderate or small, whereas their joint effect is considerably larger. The latter requires a multivariate effect size measure that is capable of computing an odds ratio for value configurations of multiple variables. Since selecting the base configuration of values -against which all other configurations are compared -is non-trivial, we utilized a value configuration relative odds ratio, i.e. a given value configuration is compared against all other possible configurations (see the Multivariate effect size measure section of methods for details).
Structural interaction results indicate higher-order interactions among members of relevant variable sets, however subsets of these variables can also be of interest. Furthermore, note that these interactions are interpreted on the level of variables and provide no further insight neither on the value level nor on the parametric level, i.e. variable value configurations. In order to investigate interactions on a parametric level additional analysis is required which involves the arbitrary selection of a subset of variables from a relevant set. For example, we can investigate parametric aspects of the interaction between body fat and sporting activities with respect to lifetime depression based on a relevant set of variables. Note that this selection is arbitrary as any subset of this relevant set could be technically investigated. Table 4 shows the parametric aspect of this interaction involving strenuous sports and exercises. In general, both physical activity types provide a protective effect with respect to lifetime depression, whereas higher body fat percentage presents a risk (OR(high versus normal) = 1.56). Regarding the joint effects of body fat and physical activity descriptors, the protective effect of doing sports or exercises is larger in case of subjects with high body fat percentage (e.g. CR-OR(Sports: No, Body fat: High) = 1.62, CR-OR(Sports:   Table 4. Parametric interactions of Body fat, Exercises and Sports with respect to lifetime depression. CR-OR and CI 95% denotes the configuration relative odds ratio and its 95% confidence interval respectively. www.nature.com/scientificreports www.nature.com/scientificreports/  Fig. 3).
Additionally, a detailed analysis of several interaction patterns is provided as additional information.

Synergistic effects.
Interactions and synergistic effects can be distinguished based on the constraints they impose on the dependency structure. Whereas interactions can be related to specific dependency structures, synergistic effects are more general in the sense that there is no such hierarchy implied concerning the dependency structure of variables. Although it is reasonable to assume that variables with distinct individual main effects are among the first to be investigated for additional synergistic effects. In our case Neuroticism, Parental depression, Current depressive symptoms, and Body fat variables can be considered as a set of relevant variables with synergistic effects. The BMLA method identified each of these variables as strongly relevant and directly connected to lifetime depression. As a set of relevant variables, they are part of the majority of possible strongly relevant variable sets. To quantify synergistic effects from a parametric aspect, a multivariate effect size analysis can be performed similarly to that which was applied in case of parametric interactions. Table 5 shows configuration relative odds ratios for various Neuroticism -Parental depression -Body fat -Current depressive symptoms subgroups. Both the neuroticism score and current depressive symptoms affect In general, the presence of parental depression entails higher risk for lifetime depression, as does high body fat percentage. These effects are more pronounced in subgroups with high neuroticism scores. In addition, the current depressive symptoms variable influences the effect size of a particular configuration to the largest extent. The before-mentioned effects of risk factors are observable especially in case of high current depressive symptom scores. The two extreme points of variable configurations in terms of effect size are: (1) subjects with high neuroticism score, high current depressive symptoms score, high body fat percentage, and with parental depression (CR-OR: 10.11), and (2) subjects with low neuroticism score, low current depressive symptoms score, normal body fat percentage, and no parental depression (CR-OR: 0.16). In terms of ordinary odds ratio, this means that subjects with the former traits are 52.61 times more likely to suffer from depression than subjects with the latter traits.     Table 6 displays joint effect sizes regarding body fat percentage and weight change. Results indicate that subjects with high body fat percentage have higher risk of lifetime depression compared to similar configurations with normal body fat percentage (see Fig. 5). Concerning weight change, gaining weight conveys a larger risk with respect to lifetime depression predictive power. In addition to relevance, the predictive power of variables can also be of interest when building a predictive model. In most cases variables that are directly relevant to a target are among those that have the highest predictive power with respect to the target. However, it should be also taken into account that (1) some strongly relevant variables are highly predictive only in a multivariate context, i.e. as part of a set of predictor variables, and (2) there are variables that are only transitively relevant but are highly predictive.
In order to investigate the predictive nature of specific variable sets and also individual variables we utilized a deep neural network based classifier having lifetime depression as the target (class) variable. For the purposes of evaluation we used an information theoretic measure (cross-entropy), residual variance, and a predictive performance measure. The reduction of residual variance due to a feature was also computed by comparing the measured residual variance to that of a random classifier. Furthermore, the predictive performance of a feature was also compared to the performance of a saturated model containing all variables (for further details see the Deep neural network based modeling section in methods). Table 7 displays measures related to some of the most significant variables in terms of association with lifetime depression. Results indicate that the neuroticism trait and current depressive symptoms are the most relevant variables as they have the best scores both in terms of reduction of residual variance (19.02% and 16.67% respectively) and of predictive performance (97.17% and 88.05% relative to the saturated model). Other variables achieve remarkably lower reduction in residual variance. Parental depression is also among the highly predictive variables (predictive rank: 3). Interestingly, the variable describing satisfaction (association rank: 3) has one of the lowest scores among these variables in terms of residual relevance reduction (1.21%, rank: 10) and also in terms of predictive performance (75.4%, rank: 11). Life stress and household income have smaller predictive power (79.41%, rank: 4 and 78.02%, rank: 7 respectively) compared to parental depression. Similarly, the BMLA method detected both variables having a direct relationship with lifetime depression but with a considerably lower probability than parental depression. To the contrary, Headache was not detected as a relevant variable even though it achieved better results in residual variance reduction (5.46%, rank: 3) than life stress (2.93%, rank: 7) and slightly better predictive power (78.53%, predictive rank: 6) than household income (78.02%, predictive rank:  Table 6. Synergistic effects of Body fat and Weight change with respect to lifetime depression. CR-OR and CI 95% denotes the configuration relative odds ratio and its 95% confidence interval respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ 7). According to the envirome map, this lack of a direct effect of headache on lifetime depression is probably due to the mediatory role of neuroticism and current depressive symptoms. In other words, despite the fact that the presence of headache has a highly significant association with lifetime depression (association rank: 5) and has considerable predictive power, these properties do not entail strong relevance.
Furthermore, the body fat descriptor is among the less predictive variables (individually) on this list (predictive rank: 9). In a multivariate context however, body fat along with neuroticism, current depressive symptoms and parental depression forms one of the most highly predictive variable sets with respect to lifetime depression. Supplementary Table S19 shows the predictive performance of this set and that of the most relevant sets of variables identified by the BMLA.
In addition, the predictive power of variable groups (e.g. social factors, financial background, etc.) was also investigated using a deep neural network classifier. In these cases only the variables related to a selected group formed the input layer of the network. Table 8 displays cross-entropy, residual variance and predictive performance measures and corresponding ranks for each variable group. Cross-entropy measured the remaining uncertainty between predicted and actual outputs of the classifier. According to results, mental health descriptors have the highest predictive power with respect to lifetime depression, followed by pain descriptors, diet and metabolism factors, life stress and parental illnesses. Note that childhood descriptors only consisted of general factors such as body size, height, and maternal smoking. Childhood trauma and maltreatment indicators (with the exception of first sexual intercourse) were not included in this analysis.   Table 8. The predictive power of variable groups with respect to lifetime depression. Reduction denotes the residual variance reduction compared to the random classifier. Ratio stands for the ratio of the predictive performance score and to that of a saturated model.
www.nature.com/scientificreports www.nature.com/scientificreports/ The superior predictive performance of the mental health group (0.713) was expected as it contained two of the most relevant variables (neuroticism and current depressive symptoms) which also had the highest predictive power. In addition, the sleep quality descriptor (Insomnia) also contributed to the predictive power.
Pain descriptors appear to be the second most predictive group (0.594) containing all pain related variables such as headache. This result indicates that various types of pain descriptors can be correlated with lifetime depression, however in the multivariate context of the environmental factors, these relationships are mediated by other factors.
The third most predictive group is diet and metabolism factors (0.582) which consists of body fat which was found directly relevant with respect to lifetime depression, and several transitively relevant factors such as weight change and the obesity descriptor. In the context of the envirome, this variable group is relevant as some of its effects directly influence lifetime depression.
Following life stress, the parental illnesses variable group is the fifth in the predictive ranking (0.570). Previous analyses revealed that parental depression is the only directly relevant variable within this group, and parental Alzheimer's disease and bronchitis may play roles in interactions of moderate relevance. As this group represents the hereditary aspect of lifetime depression, these results (i.e. several groups of environmental factors are better predictors) also confirm that investigating environmental factors is essential in predicting depression.

Discussion
Our study investigating the effects of environmental, social, lifestyle, metabolic and mental health factors on lifetime depression have shown that a surprisingly small number of core factors mediates the effects of the envirome. That is the majority of the envirome variables do not have an independent direct and relevant effect and they are only indirectly related to lifetime depression, exerting their effects in a dynamic network via transitive, interactive and synergistic relationships. This may also explain why existing environmental effects may be obscured in most studies which consider them individually in an isolated way. This result also implies that a narrowly focused set of factors can shield the effects of the whole envirome, suggesting therapeutic, clinical and pharmacological consequences. This drastic reduction of the set of relevant variables related to MDD by filtering the mediated, non-direct factors was also observed in an earlier study, which investigated the ratio of direct multi-morbidities of MDD among all the statistically associated co-morbidities 25 .
In our present study the only directly related and also highly relevant factors with respect to lifetime depression were neuroticism, current depressive symptoms, parental depression, and body fat percentage, which factors also mediated the majority of other effects in various ways. First, current depressive symptoms play a central role mediating the effects of a wide range of variables on lifetime depression including financial factors, sporting activity (partially), social factors (partially), insomnia (partially) and most pain descriptors. Second, neuroticism mediates the effects of social factors (partially), insomnia (partially), and some of the pain descriptors (headache). Third, body fat mediates the effects of most sporting activity and metabolism-related variables. In addition, body fat and parental depression play significant roles in structural interactions involving variables such as life stress, risk taking, maternal smoking, qualification, sport, alcohol consumption and household income. Furthermore, neural network-based analysis of predictive power also indicated that neuroticism, current depressive symptoms and parental depression have the highest predictive capabilities individually, while body fat had a lower predictive power. However, considered jointly as a set these four variables have the highest predictive power and form a core component of the majority of possible strongly relevant sets of variables.
Parametric analysis, besides confirming the relevance of these four variables also revealed synergistic effects between them. While high neuroticism and high current depressive symptoms had the largest individual quantitative effect on lifetime depression (OR: 9.26 and OR: 10.66 respectively), their synergistic effect coupled with presence of parental depression and high body fat had the largest multivariate odds ratio (CR-OR: 10.11) which indicates a 52.61-fold increased risk for lifetime depression compared to those with low neuroticism, low current depressive symptoms, low body fat and no parental depression (CR-OR: 0.16).
The complex pattern of relationships between the identified core factors and their role in communicating the relationship of several other factors on lifetime depression is novel compared to previous studies 14,15 . These core variables, however, have been previously implicated in association with lifetime depression in individual studies. Nonetheless, a closer look may also reveal their role in influencing the effect of other environmental variables as well. The strong relevance of parental depression with respect to lifetime depression as well as its high predictive capacity found in our study may in part reflect the significant heritability and familial aggregation of major depression 27 . However, while heritability of major depression is estimated from approximately 37% in general population samples 27 up to 75% in severely depressed recurrent depressive samples 28 , the effects of parental depression on mood disorders in off-springs go beyond genetic transmission 29 . This may determine early environmental influences including rearing and financial conditions, as well as possible early neglect and abuse on the one hand, and transmission of coping strategies and shaping character traits and behaviors by model learning on the other [30][31][32] . These factors may impact eliciting and responding to depression-relevant environmental events and stressors lasting into adulthood 33 , and may also influence future social, lifestyle and metabolic status 34,35 , although it must be noted that we had not included information on the timing of occurrence of parental depression within the lifetime of our subjects and whether it was paternal or maternal, and our data does not permit drawing conclusions concerning the direction of effects.
Similarly, neuroticism, besides its impact as a vulnerability factor for risk of depression [36][37][38][39][40][41][42][43] , with possible overlapping genetic susceptibility 4,44-46 is one of the fundamental traits of personality associated with emotional instability, negativity, increased vigilance and reactibility for negative environmental cues and a tendency for maladaptive reactiveness upon stressors 36 . Our present results confirm not only the fundamental role of Neuroticism in lifetime depression but also its central role in mediating effects of other relevant factors within the envirome.
www.nature.com/scientificreports www.nature.com/scientificreports/ Current depressive symptoms similarly emerged as showing a direct relationship with high relevance and high predictive capacity with respect to lifetime depression, an illness with a strong tendency to manifest in recurring episodes 47 . Our results indicated that the majority of effects and factors considered in our study were mediated via current depressive symptoms which can, on the one hand, be provoked and precipitated by environmental and lifestyle factors and on the other hand they are likely to determine the perception of social, lifestyle and metabolic factors, activity and functioning as well as perception of and reaction to environmental events and stressors [48][49][50][51][52] . Given the recurrent nature of major depression it is crucial to understand the influence of current depressive status including symptom profile and severity on factors determining lifetime depression as well as its position within the envirome.
Body fat percentage showed high relevance with respect to lifetime depression and mediated the effects of several related variables including stress, alcohol and tobacco consumption, physical activity and income. The association and complex relationship between obesity and depression is well-known 25,53 . While depression and consequential changes in appetite and decreased activity and motivation, as well as lower quality of life and side effects of pharmacotherapy may all be involved in weight changes in both directions but mainly weight gain in depressed patients, obesity may also be involved in the increased risk 54,55 and development of depression highlighting their reciprocal relationship 56 . Our results showing the core relevance of body fat percentage with respect to lifetime depression may also offer a possible model for understanding the association of depression with several obesity and metabolic-linked disorders such as cardiovascular disorders or type II diabetes which present serious public health challenges. Identifying the relationship between obesity and depression is among the key interests of several research studies 57 . Such analyses are typically based on the body mass index (BMI), and in some cases similar relationships concerning body fat are also explored 58,59 . Results indicated that body fat central distribution measured by waist/hip circumference ratio (WHR) had more significant associations with respect to depression symptoms than BMI. Even though BMI and WHR are closely related, these earlier results suggested that different mechanisms were responsible for the associations with depression symptoms. Our results reflect these observations since body fat was found directly relevant to lifetime depression and interacting with several other factors, whereas obesity was detected as transitively relevant, i.e. its effects on depression are mediated by other variables.
In addition to the above core factors, life stress, which has been explored in relation to depression from several aspects in the past decades 14,15,60 particularly in the context of gene-environment interactions 8,61,62 and household income reflecting financial hardship, strongly associated with depression in previous studies 63-68 , were also in a direct relationship with lifetime depression with a low but non-negligible probability. However, analysis of interactions indicated that life stress and household income had considerable effects on lifetime depression in conjunction with other factors with negligible individual effects. Life stress formed an interaction pattern with a considerable effect with lesser investigated factors including maternal smoking and parental Alzheimer's disease, while household income showed a considerable joint effect in interaction with variables reflecting physical activity and alcohol intake. While considerable research focused on the association with physical activity and alcohol intake in depression [69][70][71] trying to decipher their role as causative or consequential factors and their risk and protective effects, our study provides information on their position in the context of an envirome network of several interconnected variables which, on the longer term, could also inform research and clinical work on how they can be exploited in decreasing risk or managing symptoms of depression. All other environmental, social, lifestyle and metabolic factors investigated in our study appeared to have effects only in conjunction with our core factors.
Our results thus clearly indicate that environmental, lifestyle, social and other health-related factors in the envirome of depression act mostly in an indirect way on lifetime depression through a small number of core factors with direct effects and in transitive, interactive and synergistic relationships significantly influencing the effects of other members of the envirome network. Therefore, to understand their impact on depression, or the impact of depression on them, and exploit them for understanding, preventing, screening or treating depression as well as finding their role in restoring well-being, functionality and quality of life in depressed patients they have to be conceived together with their position and complex dynamic relationships in the envirome of depression.
The main challenge of this study is the presence of complex dependency patterns among environmental factors and large effect sizes with respect to lifetime depression. These effects result in strong associations between environmental factors and also between these factors and lifetime depression. However, not all significant associations between environmental factors and depression mean direct influence on depression. The BMLA method 8,72 is applied in this study to provide the necessary tools to distinguish between such dependency relationship types. In addition, a neural network based classifier was constructed and utilized to enable the assessment of predictive power for selected groups of factors. Our results demonstrated that environment-environment interactions need to be considered when investigating multi-factorial disorders, and multivariate methods are required to explore the relationships involving environmental factors. The presented envirome map of depression is the first comprehensive network of environmental factors that was learned from such a large cohort data. Although it must be noted that the number of included childhood descriptors is limited as childhood trauma related variables were available for only a subset of the investigated cohort and thus were not considered in the analysis.

Methods
UK Biobank data set. We analyzed a subset of the UK Biobank data set with a sample size of 110,599 involving subjects that completed the mental health questionnaire (UK Biobank Resource under Application Number 1602). The analyses were focused on lifetime depression (binary, reported by the participants and elaborated by trained research nurses) as the main phenotype descriptor (i.e. target variable). In addition, probable depression diagnosis based on experience of depressed mood and help-seeking for mental health (by Smith et al. 26 ) was utilized as a secondary target for a comparative analysis. Investigated variables (57)  Since the applied BMLA method only allows categorical variables to be analyzed, categorization was performed when needed on continuous variables. Supplementary Tables S1-S4 display the distribution of possible values for each categorical variable included in the analysis. In addition, Supplementary Tables S5-S16, describe the derivation protocol for each variable, and Supplementary Tables S17 and S18 show their individual effect size with respect to lifetime depression.
Categorization levels were based on standard thresholds when they were available, e.g. in case of systolic and diastolic blood pressure. In case of complex variables which aggregate several items, such as current depressive symptoms and neuroticism, weighted summary scores were computed and then categorized based on experts' advice. In cases where no guidelines were available, the variable levels were selected based on the distribution of values. Note that given the reasonably large sample size available in this study, we aimed to have a detailed categorization of variables where applicable, i.e. we defined multiple categories instead of having dichotomous variables. Since categorization levels were selected to correspond to practically available and semantically appropriate categories, this approach ensures that the number of levels does not contradict some essential property of the variable, and provide robust results.
The need for categorizing variables can be considered as a limitation of this study and its results. However, the main conclusions of the study should not change even if the number of levels were changed within a reasonable interval. In addition, post-hoc effect size analyses have shown that the categorization level of variables can be considered appropriate because the direction of observed effects (based on these categories) reflect previously published results.
Validation of the investigated depression phenotype. In addition to lifetime depression, we investigated relevant relationships with respect to probable depression diagnosis variables: single depressive episode, moderate depression, and severe depression created by Smith et al. 26 . We performed a comparative analysis to identify differences among relevant factors related to various phenotypes. Table 9 displays relevant factors with respect to lifetime depression and probable depression diagnosis variables treated as a group, i.e a single posterior probability of relevance is calculated for each variable with respect to the targets.
Results indicate that current depressive symptoms, neuroticism, and parental depression are highly relevant in both cases. Furthermore, several moderately relevant factors such as risk taking, life stress, age, household income, exercises, and sports have similar probabilities of relevance. A notable difference between the relevant factors of lifetime depression and probable depression diagnosis is that metabolism related variables act differently. In case of lifetime depression only body fat percentage is relevant with a high posterior probability (p Rel = 0.8), whereas in case of probable depression diagnosis that relevance is 'distributed' between body fat percentage (p Rel = 0. 25 Table 9. Relevance of variables with respect to reported lifetime depression and to probable depression diagnosis. www.nature.com/scientificreports www.nature.com/scientificreports/ the low but non-negligible relevance with respect to lifetime depression (p Rel = 0.2). On the other hand, variables such as sex, the level of satisfaction and alcohol intake are more relevant with respect to probable depression diagnosis.
All in all, lifetime depression and probable depression diagnosis appear to be adequately similar in terms of relevant factors considering that the frequency of depression cases is 5.2% for lifetime depression and 18.5% in case of probable depression diagnosis.
A Bayesian framework for identifying dependency relationships. We applied Bayesian networks in the Bayesian statistical framework to explore dependency, causal and relevance relationships between variables based on data. A Bayesian network is a probabilistic graphical model BN(G,θ), in which a directed acyclic graph (DAG) G represents the multivariate dependency relationships of variables, and parametrization θ quantitatively defines the dependency relationships by conditional probability distributions. We applied a DAG-based Markov Chain Monte Carlo (MCMC) method referred to as the Bayesian multilevel analysis of relevance (BMLA) to infer posteriors of complex structural patterns, which represent systematic, scalable repertoires of dependency, causal and relevance relationships 24,73 . In the following sections we introduce some of the corresponding fundamental concepts.
Strong relevance. Structural properties of Bayesian networks express various forms of structural relevance. An essential relevance measure called strong relevance is related to the concept of minimal Markov blanket (a.k.a. Markov Boundary Set) 74 . Assuming that a given graph G properly represents all dependency relationships defined by the data, the Markov blanket set is represented by the neighborhood of a node (variable) Y which consists of those variables X i that are in either direct structural relationship with Y or that are interaction terms. In other words, the Markov blanket set of Y 'isolates' Y from the effects of the rest of the network, such that if the values of the variables within the Markov blanket are known then no further information is required to infer the value of Y. Throughout the paper we refer to the Markov blanket set as strongly relevant set of variables which are represented by structural properties (i.e. nodes) of a Bayesian network. In the Bayesian framework a posterior probability can be induced for several types of structural properties. The posterior of strong relevance of a variable X i can be estimated by model averaging, that is assessing the probability mass of those Markov blanket sets of which X i is a member. The higher is the sum of posterior probabilities related to Markov blanket sets of which X i is part of, the more relevant X i is considered.
Technically, for each variable X i an indicator function called Markov blanket membership can be defined as , which takes the value of 1 if X i is a member of the Markov blanket of Y in DAG structure G. Assuming that the data D admits multiple possible dependency structures (i.e. models) G, and that the indicator function can be evaluated for each such G, then the posterior probability of strong relevance P(MBM(X i , Y, G|D)) can be defined by model averaging 75 based on the posteriors of possible structures 21,76 . Therefore, strong relevance can be considered as an aggregate of possible multivariate dependency models in the form of a pairwise relevance measure between X i and Y. The challenge of this approach is to properly assess the posterior probability of structural features such as Markov blanket sets.
Bayesian statistical framework. Relying on model averaging, the general task would be to compute the posterior of selected structural properties over the space of possible DAG structures G. However, due to the high cardinality of possible DAG structures, the exact computation of such posteriors is generally infeasible. Instead, various approximation methods using a Markov Chain Monte Carlo (MCMC) scheme over the space of DAGs were developed 22,23 and improved 77,78 .
In order to efficiently estimate the posteriors of complex structural patterns representing potential dependency, causal and relevance relations, especially the complex patterns of higher-order interactions, we introduced subtypes of relevances to bridge the gap between the pairwise Markov Blanket Membership and multivariate Markov Blanket Set 24 . The methodology and corresponding toolset is referred to as the Bayesian multilevel analysis of relevance (BMLA) 24 . The BMLA method can perform MCMC over the space of DAGs, PDAGs and orderings. In the current paper we applied the method in the space of DAGs using the unnormalized posterior to aid the random walk. The alteration between various DAG structures (i.e. transition between states) is facilitated by special operators deleting, inverting or inserting edges 77 .
In an optimal case, the MCMC process reaches a stationary state after a number of steps. The collected samples are utilized for the computation of posteriors only after the burn-in period, which is (approximately) the interval preceding convergence. The result of this process is a set of DAG structures which can be used to evaluate various structural properties such as direct edges or more complex properties such as Markov blanket sets. In order to compute strong relevance with respect to a target Y, the Markov blanket set of Y is identified in each sampled G structure and subsequently Markov blanket membership functions are evaluated (I MBM X Y G ( , , ) i ). By applying model averaging, the latter can be used to estimate the posterior probability of strong relevance for each X i with respect to a selected target Y. This means that each instantiation of the investigated structural property is registered and corresponding frequency information is collected, and as a last step a normalization is performed to ensure that the results can be interpreted as probabilities [0,1].
Note that BMLA can handle the redundancy of variables automatically. Redundancy from the perspective of Markov blanket sets means that e.g. variable X k is not present in a Markov blanket set of target Y if variable X i is already a member of that set. If however X i is not present in a Markov blanket set of Y then X k can be a member (assuming that it is relevant to some extent with respect to Y). If variable X k is redundant (with respect to Y) then that is reflected in a relatively lower relevance score since it is present in a smaller number of possible models.
www.nature.com/scientificreports www.nature.com/scientificreports/ Deep neural network based classifier. In order to measure the predictive power of environmental factors, and also relevant sets of environmental factors, we utilized a deep neural network 79,80 based classifier using the Tensorflow framework 81 . The network consisted of three fully connected layers using ReLU (rectified linear unit) activation functions, containing neurons matching the number of variables in the study, and an output layer. Weights and biases were initialized randomly according to the uniform distribution. Binary-cross entropy was utilized as loss function, the ADAM method was used as an optimizer 82 .
The evaluation was based on a weighted accuracy measure taking into account the unbalanced case-control ratio (i.e. subjects with and without depression). In addition, residual variance and cross-entropy was computed between the actual and predicted values. Residual variance values were compared to the residual variance of a random classifier thereby measuring the reduction in its value. Weighted accuracy measures were compared to the performance of a saturated model including all variables. All measures were computed using a k = 10 fold cross-validation, i.e. the data was split into k = 10 partitions and for each phase 9 partitions were used as training data and 1 partition as testing data. This process was repeated n = 10 times.
Multivariate effect size measure. Effect size measures are quantitative and are typically utilized in a pairwise form describing the effect of one discrete variable on another. For example an odds ratio of a variable X with respect to a target trait Y defines the extent that a change in the value of X causes in the conditional distribution of Y. By convention if the target is binary then odds are computed such that the numerator is related to the cases, and the denominator is related to controls, otherwise a base value is selected. Straightforwardly, odds ratios are computed using odds corresponding to various values of X such that one of the values is selected to be the base. If there are multiple values of X then multiple odds ratios can be defined according to the selected base. The challenge arises when the selection of a base value is non-trivial, because for example X has six possible values (x 1 , x 2 , …, x 6 ) and neither the smallest nor the largest value is a good choice for a base value. Then a possible solution is to apply value relative odds ratio which means that each odds related to a certain value of X (e.g. X = x 1 ) is compared against the odds related to all other values (e.g. X = x 2 , x 3 , …, x 6 ): Odds(X = x 1 )/Odds(X ≠ x 1 ). In this case no base value selection is required, and the effect size is more robust as more samples are utilized.
In a multivariate case, the aim is to measure the joint effect of two or more variables (e.g. X and Z) on a single target Y. Considering odds ratios with a binary target, this means that for every configuration of X and Z an odds has to be computed. However, depending on the number and cardinality of involved variables, the number of configurations can be considerably large. Therefore, defining an appropriate base value for odds ratios can be even more challenging than in the pairwise case. Instead, utilizing the previously presented idea, a configuration relative odds ratio (CR-OR) can be defined such that an odds for a given configuration of values (e.g. Odds(X = x 1 , Z = z 1 )) is compared to an odds corresponding to all other possible configurations Odds(X ≠ x 1 , Z ≠ z 1 ), that is CR-OR(X = x 1 , Z = z 1 ) = Odds(X = x 1 , Z = z 1 )/Odds(X ≠ x 1 , Z ≠ z 1 ). This multivariate effect size measure provides a robust view on the joint effect of multiple examined variables on a single target.

Data Availability
This study only utilized the UK Biobank dataset (application number 1602) under standard terms and conditions.