Educational attainment and mental health conditions: a within-sibship Mendelian randomization study

Importance: Observational studies have demonstrated consistent protective effects of higher educational attainment (EA) on the risk of suffering mental health conditions (MHC). Determining whether these beneficial effects are causal is challenging given the potential role of dynastic effects and demographic factors (assortative mating and population structure) in this association. Objective: To evaluate to what extent the relationship between EA and various MHC is independent from dynastic effects and demographic factors. Design: Within-sibship Mendelian randomization (MR) study. Setting: One-sample MR analyses included participants’ data from the Trøndelag Health Study (HUNT, Norway) and UK Biobank (United Kingdom). For two-sample MR analyses we used summary statistics from publicly available genome-wide-association-studies. Participants: 61 880 siblings (27 507 sibships). Exposure: Years of education. Main outcomes: Scores for symptoms of anxiety, depression and neuroticism using the Hospital Anxiety Depression Scale (HADS), the 7-item Generalized Anxiety Disorder Scale (GAD-7), the 9-item Patient Health Questionnaire (PHQ-9), and the Eysenck Personality Questionnaire, as well as self-reported consumption of psychotropic medication. Results: One standard deviation (SD) unit increase in years of education was associated with a lower symptom score of anxiety (−0.20 SD [95%CI: −0.26, −0.14]), depression (−0.11 SD [−0.43, 0.22]), neuroticism (−0.30 SD [−0.53, −0.06]), and lower odds of psychotropic medication consumption (OR: 0.60 [0.52, 0.69]). Estimates from the within-sibship MR analyses showed some attenuation, which however were suggestive of a causal association (anxiety: −0.17 SD [−0.33, −0.00]; depression: −0.18 SD [−1.26, 0.89]; neuroticism: −0.29 SD [−0.43, −0.15]); psychotropic medication consumption: OR, 0.52 [0.34, 0.82]). Conclusions and Relevance: Associations between EA and MHC in adulthood, although to some extend explained by dynastic effects and demographic factors, overall remain robust, indicative of a causal effect. However, larger studies are warranted to improve statistical power and further validate our conclusions.


Introduction
Mental health conditions (MHC) are mental disorders, psycho-social disabilities as well as mental states associated with significant distress, impaired functioning, or risk of self-harm 1 .Common MHC including anxiety and mood disorders (e.g., depression) are among the top causes of years lived with disability across all age groups 2 .Moreover, people with MHC have a mortality rate that is at least twice higher than that of the general population or people without such conditions 3 .
Socioeconomic circumstances including those related to level of educational attainment (EA) have been related to the risk of suffering MHC 4 .
Higher EA is associated with a lower risk of suffering MHC including, anxiety, depression, personality disorders, and substance abuse [5][6][7] .This inverse relationship might be explained by the wider availability of mental/physical resources in adulthood to cope with daily hassles among educated individuals, which could assist preventing MHC 8 .Nevertheless, it has been shown that family-related factors influence both educational level and mental health status, which suggests that alternative causal pathways in adulthood may have a role [9][10][11][12] .
Mendelian randomization (MR) is a method based on instrumental variable analysis used to investigate relationships between exposures and outcomes by using genetic variants or polygenic scores (PGS) as instruments for the modifiable factor 13 .MR relies on the premise that there is a causal pathway from an individual's genotype to the individual's phenotype (e.g., a genotypephenotype association).If its three core assumptions hold, the approach can overcome both exposure-outcome confounding and reverse causation.
To provide evidence of causality i) the genetic instrument must be associated to the exposure of interest and this relationship is required to be reasonably strong (relevance); ii) there mustn't be common factors associated to the genetic instrument or the outcome of interest (independence); iii) the genetic instrument must influence the outcome only through the exposure of interest, which 6 .implies that the effect of the genetic instrument on an outcome is fully mediated by the exposure (exclusion) 13 .Two recent population-based MR studies have provided evidence that supports the role of EA in the aetiology of MHC like anxiety and depression in adulthood 14,15 .However, since the phenotypic variation of complex traits cannot be fully attributed to direct genetic effects, such MR estimates might be biased 16,17 .In fact, violations of the independence assumption resulting from indirect genetic effects (i.e., dynastic effects) and confounding due demographic factors (i.e., assortative mating and population stratification) are now well documented for various traits including education [16][17][18] .
Dynastic effects are the result of genetic nurturing, namely the effect of parents' or other relatives' genotypes on an individual's phenotype through their contributions to the environment 17,18 .Assortative mating is when partners are not selected at random but based on particular characteristics 19 .Population structure is the presence of subpopulation differences in allele frequencies that correlate with both the phenotype and the genotype 16 .In the context of MR all three mechanisms have the potential of amplifying the effect of a given genotype-phenotype association 17,20 .
The objective of the present study is to evaluate if the association between EA and symptoms of anxiety, depression, and neuroticism as well as the consumption of psychotropic medication are explained by indirect genetic effects and/or are confounded by demographic factors.
We applied a within-sibship design coupled with genomic data to perform MR to account for dynastic effects, assortative mating, and population structure.

Study design and data sources
The present is the report of a within-sibship MR study.One-sample and two-sample MR methods were applied using individual-level data, and summary statistics from genome wide association studies (GWAS).Individual-level data came from the Trøndelag Health Study (HUNT) and UK Biobank (UKB) 21,22 .A brief description of contributing GWAS can be found in eTable 1.

Setting and participants
HUNT is a population-based cohort study that is held in the Trøndelag County in Norway and started in 1984 21 .We used data from all participants of the second (HUNT2) and third HUNT wave (HUNT3) who had been genotyped.From those, we selected all individuals who were > 30 years of age when they participated in the survey, and had at least one sibling.The final sample from HUNT2 and HUNT3 included 26 770 (10 428 sibships) and 16 718 siblings (7010 sibships), respectively (eTable 2).
The UKB is a prospective cohort study that began in 2006.UKB is following nearly 500 000, 40 to 69 years old participants from across the UK, who volunteered to be part of the study and provided consent for follow-up through linkage to their health records 22 .We included all participants who had been genotyped.After restricting the sample to sibships with two or more individuals, our analysis sample included 35 118 participants from 17 079 sibships (eTable 2).

Genetic variants
For HUNT participants we used a weighted PGS as an instrumental variable for EA (PGS-edu).The calculation of the PGS-edu was based on the genetic variants reported as significantly associated to years of education at the genome-wide level (p < 5 x 10 -8 ) in a recent GWAS 23 .From the 3952 independent genetic variants reported by Okbay et al., we included those that were well imputed in the target population (eTable 3) 24 .For further information see eMethods.

Exposure and outcomes
EA was the exposure of interest and was assessed trough the question: What is your highest level of education?for HUNT participants, or Which of the following qualification do you have?for UKB participants.A number of years of attained education was assigned for each of the answers based on The International Standard Classification of Education (ISCED) mapping 1997 (eTable 4).
Symptoms of anxiety, depression and neuroticism as well as consumption of psychotropic medication were the outcomes of interest.In HUNT symptoms of anxiety and depression were assessed by the Hospital Anxiety and Depression Scale (HADS), and in UKB by the 7-item Generalized Anxiety Disorder Scale (GAD 7) and 9-item Patient Health Questionnaire (PHQ-9) -25- 27 .In both HUNT and UKB, neuroticism was measured using the Eysenck Personality Questionnaire and consumption of psychotropic medication was self-reported 28 .A detailed description on how the exposure and outcomes were processed can be found in eMethods.

Ethics approval and informed consent
See eMethods.

Statistical analysis
Before running the analyses using individual-level data we standardized all numerical variables so that they had a mean of 0 and standard deviation (SD) of 1. EA, as well as symptoms of anxiety, depression and neuroticism were analysed as continuos, while consumption of psychotropic medication as categorical (yes/no).The association between EA and the outcomes was analysed by applying one-sample MR (two-stage least squares regression), and ordinary least squares (OLS) or logistic (LOG) regression for comparison.
All models were adjusted by sex and age.However, when the PGS-edu was included as predictor, the model was also adjusted by the first 10 principal components of ancestry (PCA) and genotyping batch to account for population structure.In all cases we assumed that standard errors were correlated within sibships and therefore clustered standard errors were computed.OLS, LOG and one-sample MR analyses were performed using the "feols" (for continuos) and "felgm" (for categorical) functions of the "fixest" package in R 29,30 .
Any difference between families due to indirect genetic effects (i.e.dynastic effects) or bias due to assortative mating was accounted for by using family fixed effects (anxiety, depression and neuroticism) or the sibling difference method (consumption of psychotropic medication) 16 .All analyses were performed using the "fixest" package in R (see eMethods for a detailed description) 30 .
Individual-level data for each of the two cohorts were analysed separately using the same model specifications and R packages.Then, results were meta-analysed using the "rma" function from the "metafor" R package 31 .We applied fixed effect models, except when heterogeneity between HUNT and UK Biobank estimates was detected, that is Cochran's Q Chi 2 < 0.05 and I 2 > 50% (eTable 5).
Two-sample MR analyses were performed using summary statistics from GWAS described in eTable 1 and the R package "TwoSampleMR" 32 .The inverse-variance estimator weighted (IVW) estimator and its 95% confidence intervals (CI) are reported in the present work.Robust estimators including weighted median, weighted mode, and MR-Egger coefficients were calculated to investigate pleiotropy and reported in supplementary tables.The directionality of the effect was evaluated using Steiger test of directionality.MR-Egger intercept test was performed to assess pleiotropy.Within-sibship two-sample MR was performed for depressive symptoms and neuroticism, as only for them we had access to summary statistics from within-sibship GWAS (eTable 1) 33 .

Handling of missing information
We imputed missing data for various questions of the HADS score as well as education and consumption of psychotropic medication in HUNT (see eTable 6 for details on missing information).We applied multivariate imputation to each HUNT survey, using fully conditional specification implemented by the MICE algorithm 34 .For each imputed data set, we then calculated the corresponding score.For further information see eMethods.

Descriptive statistics
In both cohorts participation rate among females was higher than that among males.Participants from HUNT were younger, but had similar level of education than those from the UKB (Table 1).
Symptoms of anxiety were more often reported by HUNT than from UKB participants (8% vs. 3%).
About 5% and 9% of HUNT and UKB participants reported symptoms of depression and consumption of psychotropic medication, respectively (Table 1).However, the consumption of psychotropic medication was slightly higher among HUNT than UKB participants.Differences between females and males can be seen in Table 1.For each of the cohorts, the total number of siblings varied depending on the outcome (eTable 2).

Main analyses
The direction of the regression and MR estimates was consistent among both cohorts and across all analyses.There were however some differences in the strength of the associations between HUNT and UKB (Table 2).These differences were more evident for the depression and neuroticism MR estimates.
The results of the MR analyses suggest that EA has a beneficial effect on anxiety, depression and neuroticism symptomatology (Figure 1, 3 and 4).A higher EA was also associated with lower consumption of psychotropic medication (Table 2).A SD unit increase in years of education was associated with lower scores of anxiety (-0.20 SD [-0.26, -0.14]), and neuroticism (-0.30SD [-0.53, -0.06]) as well as lower odds of psychotropic medication consumption (OR: 0.60 [0.52, 0.69]).
There was no association between EA and depression (-0.11SD [-0.43, 0.22]) (Figure 2).Estimates from one-sample MR were larger than those of the regression analyses, except from depression, but were consistent in direction between each other (Table 2).
Most of the estimates of the within-sibling MR analyses showed some attenuation in comparison to the MR estimates, the exception was depression (Table 2 and Figure 1 to 4).The effect of EA attenuated from -0.20 to -0.17 [-0.33, -0.00] for anxiety, and from -0.30 to -0.29 [-0.43, -0.15] for neuroticism.The odds ration of psychotropic medication decreased from -0.60 to -0.52 [0.33, 0.82].Some attenuation was also observed among the estimates of the within-sibling regression analyses, except from consumption of psychotropic medication (Table 2).
Two-sample MR analyses were performed for depression and neuroticism and in both cases, results were consistent with those from one-sample MR and likewise suggest a protective role of education (eTable9).Each SD increase in EA caused a -0.22 [-0.25, -0.18] and -0.21 [-0.25, -0.17] SD decrease in the scores of depression and neuroticism.Weighted median, weighted mode and MR-Egger estimators were consistent with IVW.In both cases, the directionality test indicated that the causal direction is likely to be correct, and that horizontal pleiotropy is unlikely (eTable 9).
Within-sibship two-sample MR estimates suggested an attenuation of both associations, which were however less precise, and overlapped with those of the population-based MR and the null.

Discussion
The objective of the present study was to assess to what extent the relationships between EA and anxiety, depression, neuroticism and consumption of psychotropic medication are explained by family and/or demographic factors in adulthood.Regression and population-based MR analyses indicate that higher EA is associated with lower levels of anxiety, neuroticism, and consumption of psychotropic medication, but not depression.The protective effect of EA attenuated after accounting for family and/or demographic factors.The observed attenuation was however modest, which is indicative of a causal association between EA and the studied mental health outcomes.
Various potential mechanisms that could explain a causal relationship between EA and the propensity to develop mental health symptoms have been described, and include availability of resources, cognition, knowledge and social integration 35 .Higher educated people are more likely to have higher income than less educated people, and this difference in earning can affect health by providing easier access to health care services, better nutrition and opportunities for being physical active.Similarly, EA can influence cognition, which can affect health outcomes by leading people to adopt healthier lifestyles and avoiding risk factors associated with mental health conditions such as stress.
However, as family is a driving force of adult health, the existence of other pathways than those in adulthood should be considered [9][10][11][12]36 . Suc pathways might involve environment-mediated genetic factors during early life (specially pertained to parents' characteristics such as geneticnurturing), similar but reasonably less strong genetic nurturing effect might be expected from the interaction of the individual and other relatives (dynastic-effects) 17 .
EA might also associate with mental health outcomes because of parental cross-trait assortative-mating on the exposure and outcome or on variables genetically correlated with them (e.g., highly educated individuals are more likely to select partners with lower neuroticism symptomatology) 20,37 .In some situations parental single-trait assortative-mating on a phenotype genetically correlated with both educational level and the mental health outcome might also explain their association (e.g., partner selection based on cognitive ability that is correlated with EA and the mental health outcome) 20,38 .Nonetheless, with the present design we are unable to estimate the relative influence of dynastic effects isolated from that of assortative mating or that of population structure.
Since we observed that the level of change was dissimilar among the examined MHC, it is possible that the contribution of indirect genetic effects may vary depending on the outcome of interest and its clinical severity.Accordingly, MHC whose aetiology cannot be fully attributed to direct genetic effects may be influenced by genetic nurturing to a higher extent in comparison to others with a stronger genetic basis.For example, a less pronounced change may be expected for the association between EA and attention-deficit/hyperactivity disorder or schizophrenia, since both seem to rely strongly on direct genetic effects 39,40 .

Strengths and limitations
A major strength of our work is that we applied a robust method to evaluate the impact dynastic effects, assortative mating, and population structure may have on the association between EA and the risk of suffering a MHC 16,33 .As genetic variants are randomly assigned at conception the design is generally likely to avoid bias due to confounding or reverse causality.Moreover, the application of a within-family design has allowed us to overcome the limitations of population-based MR such as non-compliance with the independence criterion 16 .Importantly, our results were replicated in two different cohorts and across different analyses: Results of observational and MR models were consistent in both cohorts.Nevertheless, our approach has limitations that must be considered when interpreting our findings.
For highly polygenic phenotypes such as EA, limitations arise from the uncertainty regarding the compliance of core assumptions including relevance and exclusion criteria 36 .Even though first-stage F-statistic indicated our analyses did not suffer from weak instrument bias, our results might still be underpowered to detect an effect of interest, as the PGS-edu only explained nearly 2% of the variation on EA.Similarly, we found little evidence of horizontal pleiotropy.
Importantly, the phenotypical expression of the genetic architecture of EA is contingent on social, historical, and cultural contexts, which might limit the informativeness of the genetic instruments in real scenarios 36 .EA in the context of MR can be thought of as propensity to achieve a given level of education, which is environmental-dependent, and hence the impact of EA on MHC may vary depending on the environmental context 36 .Moreover, MR may not be well suited to inform either novel or actionable policy recommendations, as it does not provide much insight on how or why EA is associated to any of the specific MHC we analysed 36 .
Our within-family analyses had limited power, and could not rule out a causal effect EA on none of the analysed MHC in adulthood.Additionally, as MR relies on available GWAS data (e.g., to calculate PGS-edu) our study suffers from lack of representativeness of populations from other geographical contexts and ancestry which might limit the generalizability of our findings.Our results may be affected by selection bias, as participants of large studies such as UK Biobank were non-randomly sampled 41 .However, simulations have shown that in the context of MR the impact of selection bias is likely to be less than that of other biases such as pleiotropy or population stratification 42 .

Conclusions
Our findings suggest that previously observed associations between EA and MHC in adulthood may to some extent originate and be explained by family-related factors (dynastic effects) and/or assortative mating, and population structure.However, overall all estimates remain robust, indicative of a causal protective effect of EA on the risk of suffering MHC.Future larger family based studies are needed to provide more precise evidence about these effects.
The funders had no role in the design, analysis, interpretation of the data, drafting of the manuscript, review, or final approval of the version to be published, nor in the decision to submit the manuscript for publication.

Figure 1 .
Figure 1.Educational attainment and symptoms of anxiety.Standard deviation (SD) changes in the anxiety score and its 95% confidence interval per SD increase in years of education are shown.Estimated associations are displayed for ordinary least square regression (OLS) and Mendelian randomization models.Abbreviations: SD, standard deviation unit; OLS EA, ordinary least square regression model with educational attainment as exposure; OLS PGS-edu, ordinary least square regression model with the educational attainment polygenic score as exposure; 1SMR, one-sample Mendelian randomization.

Figure 2 .
Figure 2. Educational attainment and symptoms of depression.Standard deviation (SD) changes in the depression score and its 95% confidence interval per SD increase in years of education are shown.Estimated associations are displayed for ordinary least square regression (OLS) and Mendelian randomization models.Abbreviations: SD, standard deviation unit; OLS EA, ordinary least square regression model with educational attainment as exposure; OLS PGS-edu, ordinary least square regression model with the educational attainment polygenic score as exposure; 1SMR, one-sample Mendelian randomization; 2SMR, two-sample Mendelian randomization.

Figure 3 .
Figure 3. Educational attainment and neuroticism.Standard deviation (SD) changes in the neuroticism score and its 95% confidence interval per SD increase in years of education are shown.Estimated associations are displayed for ordinary least square regression (OLS) and Mendelian randomization models.Abbreviations: SD, standard deviation unit; OLS EA, ordinary least square regression model with educational attainment as exposure; OLS PGS-edu, ordinary least square regression model with the educational attainment polygenic score as exposure; 1SMR, one-sample Mendelian randomization; 2SMR, two-sample Mendelian randomization.

Figure 4 .
Figure 4. Educational attainment and consumption of psychotropic medication.Log odds changes in consumption of psychotropic medication and its 95% confidence interval per SD increase in years of education are shown.Estimated associations are displayed for logistic regression (LOG) and Mendelian randomization models.Abbreviations: SD, standard deviation unit; LOG EA, logistic regression model with educational attainment as exposure; LOG PGS-edu, logistic regression model with the educational attainment polygenic score as exposure; 1SMR, one-sample Mendelian randomization.

Tables able 1. General characteristics of HUNT and UKB samples
. In HUNT symptoms of anxiety and depression were assessed by the Hospital Anxiety and Depression Scale (HADS), and in UK Biobank by the 7-item Generalized Anxiety Disorder Scale (GAD-7) and 9-item Patient Health Questionnaire (PHQ-9), respectively.Neuroticism was assessed using a six-item and a 12-item Eysenck Personality Questionnaire in HUNT and UK Biobank, respectively (see eMethods).