Introduction
Major depression (MD) is a common psychiatric disorder that affects more than 300 million people worldwide (World Health Organization, 2017). The Diagnostic and Statistical Manual of Mental Health Disorders-5 (DSM-5) (American Psychiatric Association, 2013) and the International Classification of Disease 11 (ICD-11) (World Health Organization, 2018) implicitly assume a categorical model for MD, whereby a collection of symptoms reflects a common dysfunction not present in healthy individuals (Helzer, Kraemer, & Krueger, Reference Helzer, Kraemer and Krueger2006; Kraemer, Noda, & O'Hara, Reference Kraemer, Noda and O'Hara2004). An alternative view is that MD is dimensional, existing along a continuum. Under this assumption, diagnostic boundaries using DSM-5 and ICD-11 criteria represent an arbitrary threshold along the distribution to partition ‘affected’ from ‘unaffected’ individuals (Vares, Salum, Spanemberg, Caldieraro, & Fleck, Reference Vares, Salum, Spanemberg, Caldieraro and Fleck2015). A categorical diagnosis is useful clinically as it allows for practical decisions regarding who to treat. In research, however, in case MD was dimensional, dichotomisation would reduce the power of a study to characterise the phenotype. Simulations have showed that categorising a phenotype with a right-skewed distribution – which would be typical of depressive symptoms – according to a clinical cut-off reduced the power to detect a genetic effect explaining 2% of the variance from 90.9% to 9.2% (van der Sluis, Posthuma, Nivard, Verhage, & Dolan, Reference van der Sluis, Posthuma, Nivard, Verhage and Dolan2013).
To explore dimensionality in MD, studies have focused on clinical correlates between sub-threshold and clinical MD (reviewed in Rodríguez, Nuevo, Chatterji, & Ayuso-Mateos, Reference Rodríguez, Nuevo, Chatterji and Ayuso-Mateos2012; Solomon, Haaga, & Arnow, Reference Solomon, Haaga and Arnow2001). These findings support a dimensional classification of MD, showing an association between subthreshold MD and increased risk for disability, impairment, comorbidities and health care use (Cuijpers, de Graaf, & van Dorsselaer, Reference Cuijpers, de Graaf and van Dorsselaer2004; Hybels, Blazer, & Pieper, Reference Hybels, Blazer and Pieper2001; Rucci et al., Reference Rucci, Gherardi, Tansella, Piccinelli, Berardi, Bisoffi and Pini2003) and a linear trend between MD severity and risk of a future episode (Kendler & Gardner, Reference Kendler and Gardner1998; Kessler, Zhao, Blazer, & Swartz, Reference Kessler, Zhao, Blazer and Swartz1997).
There are three key limitations in the literature. First, by grouping individuals into sub-threshold and clinical MD, the full continuum of the phenotype is not explored. Second, grouping individuals in this way fails to account for symptomatic heterogeneity. According to the DSM-5, one of the two ‘core’ symptoms – low mood or anhedonia – in combination with any other combination of four symptoms meets the criteria for an MD diagnosis (American Psychiatric Association, 2013). This gives rise to a possible 227 symptom profiles, with evidence that many of these profiles present clinically (Olbert, Gala, & Tupler, Reference Olbert, Gala and Tupler2014; Zimmerman, Ellison, Young, Chelminski, & Dalrymple, Reference Zimmerman, Ellison, Young, Chelminski and Dalrymple2015). Finally, the symptoms assessed are limited to the nine symptoms specified within the DSM or ICD criteria for MD. While these do reflect the core symptoms of MD, many other symptoms are commonly present. Restricting the analysis to the core symptoms therefore omits potentially relevant information to MD; for example, many patients with MD also show symptoms related to anxiety (Moffitt et al., Reference Moffitt, Harrington, Caspi, Kim-Cohen, Goldberg, Gregory and Poulton2007).
A meta-analysis of twin studies involving 21 000 individuals estimated the heritability of MD to be 37% (95% CI 31–42%) (Sullivan, Neale, & Kendler, Reference Sullivan, Neale and Kendler2000). A recent genome-wide association study (GWAS) using a broad definition of depression has identified 102 independent loci (Howard et al., Reference Howard, Adams, Clarke, Hafferty, Gibson, Shirali and McIntosh2019). These studies confirm a polygenic architecture for MD, where genetic predisposition comprises many common genetic variants of small effect that additively increase risk (Sullivan, Daly, & O'Donovan, Reference Sullivan, Daly and O'Donovan2012). GWAS results can be used to determine a single measure of genetic liability from common variants to MD, the polygenic risk score (PRS), calculated by summing the number of risk alleles carried by an individual, weighted by their effect size (Wray et al., Reference Wray, Lee, Mehta, Vinkhuyzen, Dudbridge and Middeldorp2014). How PRS associates with various characterisations of MD can be used to deepen our nosological understanding of the disorder.
This study aimed to explore the dimensionality of MD using symptom and genetic data. Exploratory and confirmatory factor analyses were used to construct multiple dimensional phenotypes for MD within the UK Biobank, a volunteer-based, national health resource. These phenotypes were tested for association with MD PRS – calculated according to a categorical definition – to understand if dimensional phenotypes contain information that would not be identified from a categorical diagnosis. Individuals were stratified by self-reported MD status to explore how the association differed between MD cases and controls.
Methods
Participants
Participants included in this study were recruited in the UK Biobank, a national health resource of 502 655 individuals aged between 37 and 73 at the time of recruitment (2006–2010) (Bycroft et al., Reference Bycroft, Freeman, Petkova, Band, Elliott, Sharp and Marchini2018). In 2016, participants were offered an online Mental Health Questionnaire (MHQ) which 157 336 participants voluntarily completed (Davis et al., Reference Davis, Coleman, Adams, Allen, Breen, Cullen and Hotopf2020).
Item selection for the dimensional phenotypes
Eighteen items were selected from the MHQ to construct the dimensional phenotypes for MD. These included all items from the Patient Health Questionnaire 9 (PHQ-9) – a measure of depressive symptoms over the last 2 weeks which correspond to the DSM criteria for MD (Kroenke, Spitzer, & Williams, Reference Kroenke, Spitzer and Williams2001); all items from the General Anxiety Disorder 7 (GAD-7) (Spitzer, Kroenke, Williams, & Löwe, Reference Spitzer, Kroenke, Williams and Löwe2006), relating to symptoms of anxiety over the last 2 weeks; and two items from the subjective well-being questionnaire, on general happiness and the meaning of life. For a detailed list of these items, see Table S1.
In total, 148 957 MHQ participants provided a non-missing response to all items above. This sample was randomly and equally split into ‘training’ (N = 74 478) and ‘test’ (N = 74 479) sub-samples. The training sample was used to develop a factor model with good fit using exploratory factor analysis (EFA), and the test sample was used to internally validate this model using confirmatory factor analysis (CFA).
Exploratory factor analysis
Polychoric correlations were computed for the 18 ordinal items in the training sample (Carroll, Reference Carroll1961). So that all correlations were positive, two items from the subjective well-being questionnaire – relating to a participant's general happiness and belief in how meaningful their life is – were reverse coded. Ordinal α (Gadermann, Guhn, & Zumbo, Reference Gadermann, Guhn and Zumbo2012) Keiser–Meyer–Olkin (KMO) (Kaiser, Reference Kaiser1974) and Bartlett's test of sphericity (Bartlett, Reference Bartlett1950) were computed. Factor analysis was considered an acceptable method if ordinal α and KMO were >0.80, and Bartlett's test with p < 0.05 in line with previous guidelines (Table S2) (Beavers et al., Reference Beavers, Lounsbury, Richards, Huck, Skolits and Esquivel2013; Gadermann et al., Reference Gadermann, Guhn and Zumbo2012).
EFA was performed using the weighted least squares method and factors were allowed to correlate using the ‘geominQ’ form of oblique rotation. Parallel analysis (Horn, Reference Horn1965) and Velicer's Minimum Average Partial test (Velicer, Reference Velicer1976) were used to produce an upper and lower bound for the appropriate number of factors, respectively. Based on these bounds, factor models were fit iteratively and compared using the criteria: Tucker–Lewis index (TLI) ⩾0.95, root mean square error approximation (RMSEA) ⩽0.05 and a smaller Bayesian information criteria (BIC) relative to other models. The model with the best fit was retained for further testing.
Items were removed according to a series of post-hoc tests using Thurstone's analytical method for simple structure criteria (Thurstone, Reference Thurstone1947) to attain a good model fit. Where multiple models demonstrated good fit, the model retaining the largest number of items was chosen. All analysis steps for EFA were performed using the psych package in R 3.4.1 (Revelle, Reference Revelle2017).
Confirmatory factor analysis
To internally validate the EFA-derived model, CFA was performed in the ‘test’ sub-sample using the lavaan package in R 3.4.1 (Rosseel, Reference Rosseel2012). Factor loadings >0.3 from the EFA model were used to specify the relationships between latent variables and items within the CFA model. As in EFA, factors were allowed to correlate. In addition to the model fit metrics used in EFA, the comparative fit index (CFI) and standardised root mean square residual (SRMR) were calculated. Models with CFI ⩾0.95 and SRMR ⩽0.05 were considered a good fit.
Deriving the phenotype: factor score calculation
Following confirmation of the proposed model, the selected model was fit to the full sample using CFA, and factor scores were computed for each factor (Fig. 1). Distributions of all factor scores are displayed in Figs S1–S3.
Flow chart displaying the methodology of the study along with the sample sizes for each section.
Polygenic risk score calculation
MD PRS were constructed using PRSice v2 (Choi & O'Reilly, Reference Choi and O'Reilly2019; Euesden, Lewis, & O'Reilly, Reference Euesden, Lewis and O'Reilly2014) in unrelated individuals of European ancestry (N = 119 692) using genotype data and quality control procedures previously described (online Supplementary Methods) (Bycroft et al., Reference Bycroft, Freeman, Petkova, Band, Elliott, Sharp and Marchini2018; Coleman et al., Reference Coleman, Peyrot, Purves, Davis, Rayner, Choi and Breen2020). This sample therefore represents a subset of the full sample for which factor scores had been calculated (N = 148 957). Summary statistics from Wray et al. (Reference Wray, Ripke, Mattheisen, Trzaskowski, Byrne and Abdellaoui2018) with 23andMe and UK Biobank samples removed (N cases = 45 591, N controls = 97 674) were used as the base dataset. To account for linkage disequilibrium, clumping was performed so that single nucleotide polymorphisms (SNPs) had an r 2 < 0.1 and a 500 kb window from other SNPs. MD PRS were then calculated across 11 p value thresholds (p < 5 × 10−8, p < 1 × 10−5, p < 0.001, p < 0.01, p < 0.05, p < 0.1, p < 0.2, p < 0.3, p < 0.4, p < 0.5, p < 1). PRS for height were used as a negative control, and were calculated in the same way as MD PRS, using summary statistics from Wood et al. (Reference Wood, Esko, Yang, Vedantam, Pers, Gustafsson and Frayling2014) (N = 253 228) as the base dataset.
Association testing
Factor scores and MD PRS were standardised, and linear regression was performed using the 119 692 MHQ participants. Each factor score was regressed on MD and height PRS with the first six genetic principal components, genotyping batch and assessment centre fitted as covariates.
Stratification by major depression diagnosis
MD cases and controls were identified using the Composite International Diagnostic Interview – Short Form (CIDI-SF), a structured self-report questionnaire focusing on depressive symptoms during an individual's worst episode of depression (online Supplementary Methods) (Kessler, Andrews, Mroczek, Ustun, & Wittchen, Reference Kessler, Andrews, Mroczek, Ustun and Wittchen1998). Individuals who did not complete the CIDI-SF, or who were removed due to the exclusion criteria (online Supplementary Methods), were assigned an ‘unknown’ status, resulting in 60 760 controls, 27 692 cases and 31 240 unknowns.
Participants were stratified into cases and controls and the same set of linear regressions repeated in each group. To test for group differences, cases and controls were combined (N = 88 452) and an interaction term for MD PRS and case–control status was included in the linear regression model, in addition to interaction terms between each of these variables and all other covariates (Keller, Reference Keller2014). Cases were further stratified into single-episode and recurrent cases with recurrence defined as having self-reported two or more lifetime depressive episodes and the interaction test was re-performed. In all interaction tests, controls were set as the reference for MD status. As a sensitivity analysis, the same tests were performed in a group comprising of both controls and unknowns.
As many of the tests are highly correlated, a method proposed by Nyholt (Reference Nyholt2004) was used to correct for multiple testing. The ‘poolR’ package (Cinar & Viechtbauer, Reference Cinar and Viechtbauer2016) in R 3.6.1 showed the effective number of independent tests for the factor scores and MD PRS was 14, giving a Bonferroni-corrected α < 0.0036 (0.05/14).
Results
Exploratory factor analysis
Initially, a six-factor model produced the best fit for the 18 items from the PHQ-9, GAD-7 and the two subjective well-being questions (Table S3). Model fit statistics did not surpass pre-defined thresholds (TLI = 0.973, RMSEA = 0.054), suggesting either poor discrimination or a lack of shared variance for particular items.
To improve model fit, items were removed according to a series of post-hoc tests based on the loadings from this six-factor model (Table S4) (Thurstone, Reference Thurstone1947). The following criteria produced a model with a good fit to the data that retained the greatest number of items: remove items with no factor loadings above 0.4, remove items with loadings to multiple factors above 0.3, retain MD symptoms from the PHQ-9 that would otherwise qualify for removal except for an item for feelings of inadequacy (test 8; Table S4). Retaining all depressive symptoms previously removed according to the cross-loading and factor loading criteria did not produce a good fit (test 5; Table S4). As such, a leave-one-out scenario was undertaken whereby depressive symptoms that had been previously removed due to the cross-loading and factor loading criteria in test 4 (psychomotor retardation/agitation, impaired concentration, feelings of inadequacy and suicidal thoughts) were removed in ascending order of the items communality – a measure of the shared variance of an item with all other items in the model (tests 6–8; Table S4). The first item removed that resulted in a good model fit was feelings of inadequacy (and so a test removing suicidal thoughts was not performed). These criteria retained 15 of the original 18 items and a five-factor model produced a model with good fit (TLI = 0.981, RMSEA = 0.048). Items relating to trouble relaxing, irritability and feelings of inadequacy were excluded.
Confirmatory factor analysis
In contrast to EFA which allows all items to load on all factors, CFA specifies the relationship between items and factors directly and provides a more stringent test for the proposed model (Thompson, Reference Thompson2004). CFA within the test sub-sample (N = 74 479) confirmed a good fit to the observed covariance matrix (CFI = 0.995, TLI = 0.993, RMSEA = 0.033, SRMR = 0.025). As correlations between the five factors were moderate to high (r = 0.574–0.848), a second-order latent variable was included in the model; the model fit remained within the pre-specified thresholds (CFI = 0.992, TLI = 0.99, RMSEA = 0.038, SRMR = 0.031). See Tables S5a–S5f for detailed results on the model fit, factor loadings and correlations. Fit statistics did not change when fitting the model into the full sample (N = 148 957) using CFA (Table S5g).
In the final model, the five first-order factors reflected feelings of anxiety, psychomotor-cognitive impairment, neurovegetative states, mood and subjective well-being. The second-order factor, representing the correlation between these factors, was termed the ‘internalising’ factor (Fig. 2). Calculated as the average squared factor loading, the variance explained by each first-order factor ranged from 49% to 70%. The internalising factor explained 72% of the variance for all first-order factors. The high degree of variance explained supports the use of these factors as phenotypes for dimensional MD.
Factor model used to derive the dimensional phenotypes. As is customary in structural equation modelling graphs, circles are factors and squares are the self-reported symptoms. Shaded areas relate to either core MD symptoms or factors containing a majority of MD symptoms. Arrows pointing from either one factor to a symptom or a factor to another factor represent the factor loadings. *The items ‘General Happiness’ and ‘Meaningful Life’ have been reverse coded such that they explore ‘general unhappiness’ or ‘lack of meaning in one's life’. Subjective well-being, therefore, also corresponds to a ‘subjective lack of well-being’. Nomenclature has been retained for the brevity of the labelling.
Major depression polygenic risk score associates with an individual's factor score
MD PRS were associated with each of the factor scores (β range: 0.056–0.064; p value range: 2.57 × 10−82–1.89 × 10−107; MD PRS p value threshold: p T < 0.3). While the effect sizes were similar across factors (Table S6; Figs S4–S9), the two factors with the lowest association with MD PRS related to feelings of anxiety and subjective well-being (Fig. 3). In contrast, height PRS was not associated with any factor score (p > 0.05) at any p value threshold (Table S7).
Association of MD PRS on each factor in the full sample and when stratified by MD Case/Control Status. The MD PRS used was calculated at the p value threshold of p T < 0.3. Both MD PRS and Factor scores were standardised to have a mean of 0 and variance of 1 using the full sample.
Stratification by case/control status
MD PRS remained associated with factor scores following stratification into case- and control-only subgroups (controls: β range: 0.022–0.025; p value range: 1.89 × 10−13–1.43 × 10−17; cases: β range: 0.042–0.055; p value range: 7.54 × 10−10–8.95 × 10−17; MD PRS p value threshold: p T < 0.3). Figure 3 shows the β coefficients for MD PRS on each factor for each group (cases, controls and unknowns) (Tables S8–S10; Figs S4–S9). The p value thresholds of MD PRS with the largest effect size differed across factors between cases and controls (case threshold: p T < 0.2; control threshold: p T < 1) suggesting SNPs associated with MD contained greater signal for dimensions within cases whereas dimensions within controls required additional SNPs with weaker evidence for association in the base data.
MD PRS had an attenuated effect in controls relative to cases (Fig. 3). To formally test for a differential genetic effect, factor scores were regressed on MD PRS and MD diagnostic status (case or control; N = 88 452) with an interaction term. An interaction was detected for all factors (β range: 0.020–0.032; p value range: 2.23 × 10−3–3.94 × 10−7; MD PRS p value threshold: p T < 0.3; Table 1) indicating the genetic contribution is potentiated in cases relative to controls. The effect size of MD PRS with both controls and unknowns in a single group was greater than the association of the two individual groups (Fig. S10). This aggregate group showed no interaction effect of MD PRS with MD case status (β range: 0.001–0.01; p value >0.05; MD PRS p value threshold: p T < 0.3).
Main and interaction effects of MD PRS and MD diagnostic status on an individual's standardised factor score
The full sample has been subset to only include individuals with MD diagnostic status (N = 27 692 cases; N = 60 760 controls). Cases have been coded as 1 and controls as 0 (reference). Main effect models for MD PRS and MD diagnostic status include the relevant variable and covariates. The main effects are shown to evidence their statistical significance, often a pre-requisite to performing a statistical interaction.
* p value not possible to determine through lm() function in R as it is below the floor limit of the software. As such it is specified to simply be under a specific value given in the summary results from R (p < 2 × 10−16).
Cases were stratified into single-episode (N = 10 590) and recurrent cases, reporting two or more lifetime depressive episodes (N = 10 726). The interaction test was repeated with controls as the reference category. No interaction effect was detected between controls and single-episode cases (p > 0.05) for any factors. A nominally significant interaction was found between controls and recurrent cases for all factors, except anxiety and subjective well-being. No results survived correction for multiple testing (Tables 2 and 3).
Main effects of MD PRS and MD diagnostic status split into single-episode and recurrent cases
The full sample has been subset to only include individuals with MD diagnostic status who report the number of lifetime episodes of depression (N = 10 590 single-episode cases; N = 10 726 recurrent cases; N = 60 760 controls). Single-episode cases have been coded as 1, recurrent cases as 2 and controls were coded as 0 and set as the reference category. Main effect models for MD PRS and MD diagnostic status include the relevant variable and covariates. The main effects are shown to provide evidence for their statistical significance, often a pre-requisite to performing a statistical interaction.
* p value not possible to determine through lm() function in R as it is below the floor limit of the software. As such it is specified to simply be under a specific value (p < 2 × 10−16) given in the summary results from R.
Interaction effects of MD PRS and MD diagnostic status split into single-episode and recurrent cases
The full sample has been subset to only include individuals with MD diagnostic status who report the number of lifetime episodes of depression (N = 10 590 single-episode cases; N = 10 726 recurrent cases; N = 60 760 controls).
Discussion
The aim of this study was to explore polygenic associations across the continuum of MD and to test if this association differed between cases and controls. Using the UK Biobank, this study shows that the polygenic liability for a categorical MD phenotype (Wray et al., Reference Wray, Ripke, Mattheisen, Trzaskowski, Byrne and Abdellaoui2018) also associates with a dimensional model of depressive symptoms. Moreover, this finding holds in analyses stratified by case–control status. This suggests PRS contains information over and above the risk of becoming a case and may also be used to indicate severity across the continuum.
To account for symptom-level heterogeneity within the dimensional phenotype, this study used factor analysis to derive a five-factor hierarchical structure for MD. Previous studies investigating the latent structure of MD have produced multiple factor solutions. One- and two-factor models have been shown to produce the best fit for the PHQ-9, depending on the sample selected, i.e. case-only or population cohort (Elhai et al., Reference Elhai, Contractor, Tamburrino, Fine, Prescott, Shirley and Calabrese2012; Kocalevent, Hinz, & Brähler, Reference Kocalevent, Hinz and Brähler2013). The five-factor hierarchical model derived in this study has a high level of agreement for depressive symptoms with a model proposed by Kendler, Aggen, and Neale (Reference Kendler, Aggen and Neale2013) in a population-based sample of 7500 twins, which showed three uncorrelated genetic factors best decompose the phenotypic variance in lifetime MD symptoms. Our model differed in two ways, firstly the symptom of feelings of worthlessness or excessive guilt was not included; and suicidal thoughts loaded on the mood factor, whereas it loaded on the psychomotor/cognitive factor in the model proposed by Kendler et al. (Reference Kendler, Aggen and Neale2013). Our study may therefore be considered a quasi-replication that used phenotypic rather than genetic covariance, in a substantially larger dataset (N = 119 692), with current in contrast to lifetime symptoms, to provide support for the multidimensionality of MD. Multidimensionality also indicates a deeper level of complexity in MD that is left unexplored when using either sum-scores of symptoms or a case–control design. A logical extension would be to use this multidimensionality to determine which genetic variants influence each factor; however, the high loadings of the five first-order factors onto the second-order ‘internalising’ factor suggest a significant proportion of this complexity is shared. As a result, it is likely that larger sample sizes will be required to identify genetic variants specific to a given factor.
The separation between symptoms of MD, anxiety and subjective well-being is noteworthy. The GAD-7 has previously been shown to possess a unidimensional factor structure (Löwe et al., Reference Löwe, Decker, Müller, Brähler, Schellberg, Herzog and Herzberg2008) distinct from symptoms of MD (Spitzer et al., Reference Spitzer, Kroenke, Williams and Löwe2006). Compared with the ‘core MD factors’, an interaction between controls and recurrent MD with MD PRS was not detected for the factors relating to anxiety and subjective well-being at the level of nominal significance. This may suggest that although these symptom dimensions contain a highly pleiotropic genetic component (Purves et al., Reference Purves, Coleman, Meier, Rayner, Davis, Cheesman and Eley2019), a diagnosis of MD contains a degree of specificity which reflects the structure suggested by the DSM and ICD. However, as no factors survived correction for multiple testing for the recurrent case interaction, this conclusion warrants further investigation.
When stratified into cases and controls for MD, an attenuated association between MD PRS and the continuous phenotypes was evident for controls relative to cases. This is, perhaps, not surprising as controls were screened for the presence of any psychiatric disorders and high levels of current depressive symptoms (PHQ-9 sum-score <14). As such, compared to cases, controls are expected to form a more homogenous group of healthy individuals, limiting the power for MD PRS to associate with the phenotype. As can be seen in Figs S4–S9, the variance in factor scores in controls is approximately half of that found in cases as scores are condensed heavily towards the left tail. However, the variance in the MD PRS remains approximately equal across all three sub-groups. Sampling from a smaller subset of the total variance within the factor scores will drive the association towards the null. This conclusion is further supported by the removal of any interaction effect of MD diagnostic status when controls and unknowns were combined into a single group (Fig. S10). Nevertheless, a significant association was still evident, suggesting that even in what would typically be considered a healthy, ‘super-control’, group, MD PRS can still differentiate the subtle differences of the continuous phenotype. The larger control sample size limits comparisons of significance between cases and controls; however, the core finding is of significance within groups, not a comparison between groups.
In contrast, the effect size within cases was similar to that of the entire sample, indicating cases contain the majority of the signal for the dimensions. This finding has important implications under the assumption of a purely dimensional phenotype as it suggests ignoring the variation within cases, also ignores a substantial proportion of the association. As highlighted above, the increased association within cases may be primarily due to the increased variability in factor scores within the group. Alternatively, questionnaires may be being interpreted differently between cases and controls for MD, perhaps due to a greater degree of familiarity with the questionnaires by people who have been diagnosed with depression. Familiarity may increase the validity of the responses as individuals are more ‘in tune’ with the symptoms, reducing measurement error and improving power for the study to detect an association.
The association of MD PRS with the factor scores increases as the three subgroups are coalesced, showing the largest association when all participants are included (Fig. 3; Fig. S10). This is consistent with a model in which the mean MD PRS and factor score differ by diagnostic subgroup and the variance of the factor scores within the groups is smaller than that of the full sample. A simulated example to demonstrate the reasoning is shown graphically in Fig. S11. Importantly, the different associations between cases and controls do not provide evidence either for or against a categorical phenotype for MD. Further studies are warranted to explore this within diagnostic sub-group variation for differential genetic risk.
Whilst the presence of an association within cases and controls for MD appears to contradict the purely categorical phenotype of MD – differences within cases and controls are not associated with risk factors – this study cannot exclude the possibility that MD is a categorical phenotype characterised by continuous variation within cases and controls. Two taxometric studies, designed specifically to detect the presence of such groups or ‘taxons’, provide support for this finding (Ruscio, Brown, & Ruscio, Reference Ruscio, Brown and Ruscio2009; Ruscio, Zimmerman, McGlinchey, Chelminski, & Young, Reference Ruscio, Zimmerman, McGlinchey, Chelminski and Young2007). However, this field of research has consistently supported a dimensional classification, and the results appear to depend on the measurement instrument used, whether the symptoms were self-report or clinically ascertained and age of the sample (Hankin, Fraley, Lahey, & Waldman, Reference Hankin, Fraley, Lahey and Waldman2005; Liu, Reference Liu2016).
Limitations
This study has many strengths including large sample size from a volunteer-based, national health resource and its accounting for the heterogeneity inherent to MD; however, important limitations remain. The items used to create the dimensional phenotypes were self-reported, increasing the risk of misclassification and sampling bias. It has been shown that participants who responded to the MHQ have a higher level of education and fewer hospital diagnoses inclusive of mental disorders compared with other UK Biobank participants (Adams et al., Reference Adams, Hill, Howard, Dashti, Davis, Campbell and McIntosh2019). This ‘healthier and wealthier’ bias may hamper our ability to appropriately represent cases at the most severe end of the spectrum. The same study showed that MHQ responders were more likely to have a family history of severe depression relative to non-responders (Adams et al., Reference Adams, Hill, Howard, Dashti, Davis, Campbell and McIntosh2019). Similarly, the influence of personal interest in mental health could limit the generalisability of this sample to the general population.
A combination of theory on the depression phenotype and minimum threshold criteria for fit statistics was imposed to choose the factor model. Thresholds are driven by best practice; ultimately, however, they are arbitrary. With a different set of thresholds, another model may have been chosen during the exploratory phase of the analysis. The authors recognise this but believe this model to make theoretical sense and have a strong fit to the observed covariance matrix as strict thresholds have been chosen.
This study assessed current depressive symptoms at a single time point, when the MHQ was completed. Studies investigating latent class trajectory of depressive symptoms have supported the dynamic nature of MD, with trajectories including persistently low, persistently high, increasing and decreasing symptoms (Byers et al., Reference Byers, Vittinghoff, Lui, Hoang, Blazer, Covinsky and Yaffe2012; Kuchibhatla, Fillenbaum, Hybels, & Blazer, Reference Kuchibhatla, Fillenbaum, Hybels and Blazer2012). Future studies should seek to confirm this model in longitudinal settings to test if it is robust to temporal invariance (Widaman, Ferrer, & Conger, Reference Widaman, Ferrer and Conger2010). Furthermore, the distinction between factors might reflect an artefact of time as the items are from three questionnaires (PHQ-9, GAD-7 and subjective well-being). For example, if the PHQ-9 was taken at the start of the MHQ and the GAD-7 at the end, subtle shifts in mood may occur while completing the MHQ that could artificially reduce the correlations between depressive and anxiety symptoms. This limitation is somewhat mitigated from the fact that the MHQ was completed by most participants in a short period with 82% having completed the questionnaire within 25 min (Davis et al., Reference Davis, Coleman, Adams, Allen, Breen, Cullen and Hotopf2020).
Conclusions
MD PRS supports a multi-dimensional model of MD, indicating that information is contained within cases and controls that would otherwise be omitted using a categorical phenotype. Much of this additional information is held within cases; possibly due to the greater variation in factor scores found within this group. Considering this additional variation in future study design may enhance the power to detect genetic associations, elevate our current aetiological understanding and improve prediction through more accurate PRS.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291720001828.
Acknowledgements
CML is funded by the Medical Research Council (N015746/1). SPH is funded by the Medical Research Council (MR/S0151132). DMH is supported by a Sir Henry Wellcome Postdoctoral Fellowship (Reference 213674/Z/18/Z) and a 2018 NARSAD Young Investigator Grant from the Brain & Behavior Research Foundation (Ref: 27404). This study represents independent research funded by the National Institute for Health Research (NIHR) Biomedical Research Centre at South London and Maudsley NHS Foundation Trust and King's College London. The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care. We thank participants and scientists involved in making the UK Biobank resource available (http://www.ukbiobank.ac.uk/). UK Biobank data used in this study were obtained under approved application 18177.
Conflict of interest
Cathryn M Lewis reports having received fees from Myriad Neuroscience. Bradley S Jermy, Saskia p Hagenaars, Kylie p Glanville, Jonathan RI Coleman, David M Howard, Gerome Breen and Evangelos Vassos reported no biomedical financial interests or potential conflicts of interest.
