{'fraction_non_alphanumeric': 0.04133046756632959, 'fraction_numerical': 0.024699622987209752, 'mean_word_length': 5.058057630736393, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Depressive symptom trajectories and polygenic risk scores in individuals with an immune-mediated inflammatory disease

Objective: To develop group-based trajectories of depressive symptoms in immune-mediated inflammatory disease (IMID) to understand their evolution and identify any associated factors, with the overall goal of identifying those at highest risk of higher depressive symptom burden. Method: 922 participants had an IMID or anxiety/depression. The PHQ-9 was administered at four visits, and polygenic risk scores (PRS) for major depressive disorder, depressive symptoms, and body mass index (BMI) were generated. Group-based trajectory modelling of PHQ-9 scores estimated distinct trajectories. Regression tested whether specific factors were associated with the trajectories. Mediation analyses assessed whether IMID mediated the association between BMI PRS and trajectories. Results: Three trajectories were identified. Regression demonstrated those in Group 3 (‘high symptoms’) had significantly higher PRS for the three traits, compared to Group 1 (‘minimal symptoms’) (OR: 1.34–1.66, P < 0.01). Stratified analyses in the IMID subgroup revealed an increased effect for BMI PRS in Group 3 (OR: 2.31, P < 0.001), in contrast, BMI PRS was no longer associated in the non-IMID sample. No significant indirect effect of BMI PRS on depressive symptoms trajectories was identified via IMID. Conclusions: A significant association between polygenicity and PHQ-9 trajectories supports a role for genetic inheritance in the variability in depressive symptoms in IMID.


Introduction
Major depressive disorder (MDD) affects approximately 5.0% of adults globally [1]. A two-to three-fold increase in the prevalence of MDD is observed in individuals with an immune-mediated inflammatory disease (IMID), such as multiple sclerosis, inflammatory bowel disease or rheumatoid arthritis [2]. Individuals with comorbid IMID and MDD have an increased hazard for all-cause mortality, compared to those with an IMID and no MDD [3]. In a cohort study of 2312 people with multiple sclerosis, significantly greater annual disability progression was found in those with a mood or anxiety disorder (beta = 0.28, P = 0.0002) [4]. In those with inflammatory bowel disease, the co-occurrence of MDD is associated with increased disease relapses [5], treatment failure [2,6] and reduced quality of life [7]. Comorbid MDD is also associated with lower cognitive function, as evidenced by slower processing speed, and lower verbal learning and memory in individuals with IMID [8].
Several studies have also suggested that anxiety and depressive symptoms are more likely to be elevated during periods of increased IMID disease activity and improve when disease activity is reduced [9]. Elevated depressive symptoms, even in the absence of meeting formal MDD diagnostic criteria, can also adversely affect outcomes in the general, non-IMID population [10]. The Patient Health Questionnaire-9 (PHQ-9) comprises nine questions based on the Diagnostic and Statistical Manual of Mental Disorders (DSM-IV) diagnostic criteria, and is a validated tool for assessing depressive symptoms severity both in the general population [11] and in persons with IMID [12][13][14]. A qualitative study showed that PHQ-9 patterns reflected patients' experience of depression severity over a 6-month period [15]. Changes in PHQ-9 scores are consistent with changes in the clinical diagnosis of depression over time [16]. The factors that influence depressive symptoms in IMID are poorly understood, but may include genetic factors as is the case in the general population [17]. Polygenicity refers to a trait in which a large number of genetic variants are involved in explaining its heritability [18]. In many common diseases, the effect of polygenic inheritance is comparable to that of rare monogenic mutations [19]. MDD is no exception, with an estimated heritability from twin studies of 37% [18] and of single nucleotide polymorphism (SNP)-based heritability of 8.7% [20]. The polygenic inheritance can be quantified using a polygenic risk score (PRS), representing the weighted sum of the number of risk alleles that an individual possesses for a trait [21]. Such a score can be utilized to predict the disease status of an individual.
The relationship between MDD and the immune system is complex [22] and is further complicated by a risk factor common to both MDD [23] and IMID [24][25][26]: body mass index. Studies employing Mendelian randomization, a method used to examine causal effects of a modifiable exposure on an outcome, have shown BMI, but not MDD, play a causal role in determining the likelihood of developing two IMID: multiple sclerosis [24] and rheumatoid arthritis [27]. Other Mendelian randomization studies have shown both BMI [26] and MDD as risk factors for inflammatory bowel disease, but the investigation into MDD did not account for BMI [28]. There is also some distinction between depressive symptoms as measured by the individual PHQ-9 questions and MDD genetic architecture as noted by a genetic correlation of <1 [29].
In the general population, high depressive symptoms trajectories are associated with poor social and psychiatric outcomes [30]. Therefore, identifying those in a higher depressive symptom trajectory is important for targeting individuals who may benefit from more intensive interventions to mitigate poor outcomes. In the non-IMID population, we know depressive symptom trajectories can be heterogeneous [30][31][32][33], but to the best of our knowledge, there are no previous longitudinal studies investigating depressive symptoms in IMID. Factors influencing depressive symptom trajectories in IMID have also not been investigated previously, but in the general population, gender, income/education and genetic factors have been associated with greater depressive symptoms [17,30]. Therefore, our objective was to develop trajectory models of depressive symptoms in individuals with an IMID, or depression or anxiety history and no IMID to determine the pattern of trajectories. We then aimed to investigate whether non-genetic and genetic (i.e. PRS for MDD, depressive symptoms and BMI) factors were associated with depressive symptom trajectories, with the overall goal of identifying factors that could be utilized to identify those at highest risk of higher depressive symptom burden.

Study population
As previously reported, participants residing in Manitoba, Canada [2] were recruited between November 2014 and July 2016 for a large cohort study of psychiatric comorbidity in IMID. Individuals were grouped as: IMID (one of multiple sclerosis; rheumatoid arthritis; or inflammatory bowel disease, including Crohn's disease and ulcerative colitis), or had a lifetime history of anxiety or depression disorders but no IMID. We used multiple recruitment methods such as poster placement in hospitals, private medical clinics, and educational institutions. Other forms of recruitment included in person or mail recruitment of patients from community-based and tertiary care clinics. Participants were ≥ 18 years of age, provide informed consent, and were sufficiently proficient in English to complete questionnaires. The University of Manitoba Health Research Ethics Board (HS17706/ H2014:201) and Winnipeg Regional Health Authority (2015-051) approved the study.

Measures
2.2.1. Sociodemographic and clinical characteristics-Self-report questionnaires captured the following characteristics at the baseline visit: sex, age, BMI, annual household income, highest level of education attained, race, and smoking history. To ensure reasonable cell sizes, the following were categorized: annual household income in Canadian dollars (<$50,000, ≥$50,000, or "Decline to answer"), educational attainment (high school or below, above high school), and self-reported race (white, non-white). Non-white selfreported race was further classified into the following categories if N > 5: Asian/Filipino/ Latin American/African/Canadian Indigenous. Participants who reported having smoked ≥100 cigarettes in their lifetime were categorized as ever smokers [2]. Participants also reported their age at symptom onset from which we calculated disease duration for either the IMID or anxiety or depression.

Psychiatric morbidity-
The PHQ-9 is a validated tool of 9 self-report items used to assess depressive symptom severity over the last two weeks. Each of the nine items has response options from 0 (not at all) to 3 (nearly every day). The total score ranges from 0 to 27, with a total score ≥ 10 indicating the presence of depression [34]. The PHQ-9 tool was administered at the baseline study visit and then annually for 3 years, for a total of 4 visits. Two PHQ-9 cut-offs that have been validated for the presence of depression in IMIDs were also utilized: PHQ-9 ≥ 11 and 12 [12][13][14]. Lifetime MDD status was determined using the Structured Clinical Interview for DSM-IV (SCID), as the DSM-IV was the prevailing diagnostic criteria at the time the study was designed [2].

Genotyping and polygenic risk scores-Peripheral blood was collected for
DNA extraction at one of the study visits. Samples were genotyped on the Illumina GSA-MD SNP array (GSAMD-24v2) at the Genome Québec Innovation Centre (McGill University, Montreal, Quebec). Quality control and imputation details can be found in the Supplement.
Polygenic risk scores (PRS) were calculated using summary statistics from recent genome wide association study for each phenotype: MDD (N = 170,756 MDD, 329,443 controls) [40], depressive symptoms (N = 180,866) [17], and BMI (N = ~681,275) [41]. These traits were selected for PRS generation given the trait under study (depressive symptoms) and because BMI is known to be phenotypically and causally associated with depression [23,42,43]. PRS were generated for the samples as the sum of the risk allele scores, weighted by their effect size in the discovery GWAS samples [40]. We performed linkage disequilibrium clumping (r 2 < 0.1 in 1-Mb window) on any overlapping SNPs with the 1000 Genomes Project European samples for the reference. PRS were calculated using PLINK (v1.9) for eight different scores based on different p-value thresholds (p ≤ 5 × 10 −8 , p ≤ 1 × 10 −5 , p ≤ 1 × 10 −3 , p ≤ 0.01, p ≤ 0.05, p ≤ 0.1, p ≤ 0.5, and p ≤ 1). PRS were standardized as mean of 0 and standard deviation of 1 for interpretability, with higher values indicating increased genetic burden of the respective trait.

Statistical analyses
The population was described regarding their characteristics using either the median (interquartile range, IQR), mean (standard deviation, SD) or frequency (%).

Depressive symptoms trajectory analysis-Trajectories
were estimated for each patient by group-based modelling analysis implemented in SAS 9.4 (PROC TRAJ, SAS Institute Inc.) [44]. This semiparametric technique was used to identify subgroups of individuals with distinct longitudinal trajectories for PHQ-9 scores. This method has been used to model longitudinal depressive symptoms in primary care and in the older population [31,45,46]. Individuals were included if they had ≥2 PHQ-9 scores to capture the change in depressive symptoms over time. There were four study visits and for a study visit to be included, we required ≥50% of the cohort to have completed the PHQ-9 at that visit. A censored normal (CNORM) model was used for PHQ-9 to estimate the depressive symptoms trajectories. Model selection was based on the differences of Bayesian Information Criterion (BIC) between the complex and null models. As suggested by Nagin [47], model fitting began with one group with a second-order (quadratic) pattern and then the number of groups and polynomial orders were added until the BIC was minimized. The polynomial orders reflect the shape of the pattern of change for each group over time (0 = zero-order; 1 = linear; 2 = quadratic, 3 = cubic). The polynomial orders were considered based on whether the highest polynomial's coefficient for each trajectory group was significantly different from 0. The number of groups was limited to a minimum class membership of 10% of the total population. Following modelling, the posterior probability threshold was set to 0.7 to determine a correct group membership assignment [47], with any participants with <0.7 excluded from further analyses. We compared the characteristics between those participants who were included in the analyses (i.e. those with a posterior probability ≥0.7) to those excluded (i.e. those with posterior probability <0.7) using the appropriate statistical tests (categorical: chi-square test, continuous: Student's t-test or Wilcoxon rank sum test).

Multinomial regression modelling-
To test whether various genetic and nongenetic factors were associated with differing depressive symptom group-based trajectories, we used univariate and multivariable multinomial regression. Although MDD and depressive symptoms PRS capture similar dimensions of depression, they were both included in the modelling because of low multicollinearity (variance inflation factor ~ 1.2 for both PRS).
To determine which p-value threshold of the PRS to include in the regression modelling, we computed the variance in group-based trajectories [Nagelkerke's pseudo-R 2 [48]] explained by PRS as the difference in R 2 from the full model including the PRS and a baseline model with ancestry principal components only (Table S1). For the multivariable model, covariates were tested and added into the model when significant either by univariate analyses (P ≤ 0.1) or based on the literature. We also investigated the association between PRS quartiles and group-based trajectories to compare those at the extremes of the PRS distribution. Quartile modelling included the same covariates as the multivariable model and designating the first quartile as the reference for each PRS. To determine the effect of disease status (IMID vs. non-IMID) on the regression, we performed a complementary regression stratified by disease status (IMID vs. non-IMID).

Causal mediation analysis-We performed complementary mediation analysis
to test whether IMID mediates the association between BMI PRS and depressive symptoms trajectories. We expected that we would not observe mediation based on prior Mendelian randomization studies suggesting BMI is causally related to IMID and to MDD, but IMID is not causally associated to MDD [24][25][26][27][28] (Fig. S2). However, as MDD and depressive symptoms can capture various dimensions of depression, the previous Mendelian randomization studies in MDD may differ for that of depressive symptoms. We set the statistical significance level at P ≤ 0.05. Missing data were not imputed. All analyses were performed using SAS (v.9.4) or R for Statistical Computing (v.4.1.2) [49] with the following R packages used: tidyverse [50], data.table [51], nnet [52], DescTools [53], mediation [54], and cowplot [55].

Results
Of the 922 study participants, we included 767 European genetic ancestry study participants with available genetic data and ≥ 2 PHQ-9 scores (Fig. 1). At each time point, patient completion of the PHQ-9 scores were as follows: baseline (97.8%, n = 17 missing), year 1 (98.8%, n = 9 missing), year 2 (94.1%, n = 45 missing), and year 3 (90.0%, n = 77 missing); thus all 4 visits fulfilled the 50% completeness requirement. After fitting various group-based trajectory models to minimize the BIC, a 3-group model was selected (Table S2). From the 3-group model, we then excluded 70 individuals because of a low posterior probability (<0.7) of being assigned to one of the three groups. These 70 excluded participants did not differ largely from the 697 participants included in the subsequent analyses (Table S3). The remaining 697 individuals were assigned to one of three PHQ-9 trajectory groups as follows: Group 1 represents 52.6% of participants, Group 2 with 35.3% and Group 3 with 12.1% (Fig. 2A). The mean posterior probabilities for each group were: 0.997 (Group 1), 0.976 (Group 2), 0.998 (Group 3), with a Wald test demonstrating that the three groups were significantly different (X 2 = 111.56, P = 1.64 × 10 −24 ). The three groups were characterized by the level of depressive symptoms, which included minimal (Group 1), moderate (Group 2), and high (Group 3). The three groups were also characterized by stability over time. From here on, we will refer to these three trajectories by the Group number (1, 2 or 3).
Of the 697 individuals included in the further modelling, there was a preponderance of women (76.3%), self-reported white race (96.4%), and higher education (69.3%), with a median baseline age of 51.4 years, and range of 18.3-84.0 years (Table 1). Within each trajectory group, there were similar patterns of these socio-demographic factors. The mean PHQ-9 scores of the three groups were similar by the presence of an IMID or non-IMID (Fig. 2B), with increasing rates of lifetime major depression across the three groups (Table  1).
We used multinomial regression modelling to test whether the non-genetic and genetic factors were associated with depressive symptom trajectories (Table S4). In the univariate analyses, younger age, lower income, lower education, and absence of an IMID were significantly associated with membership in either groups 2 or 3 compared with group 1 (Table S4). Significant findings associated with only group 2 membership included women and ever-smoking.
In the multivariable analyses, the following were significantly associated with higher odds of membership in groups 2 or 3, compared with group 1: younger age, ever-smoking, lower income and non-IMID disease status (Table S4). Those in group 3 had significantly higher odds of reporting lower education, compared with group 1. The only PRS that was significantly associated with a depressive symptom trajectory was that of BMI (Table S4). We then examined the effects of individual factors on the association between BMI and MDD PRS and the depressive symptom trajectory groups, by separately adjusting for the various genetic and non-genetic factors (Table S5A). For BMI PRS, adjusting for other genetic and non-genetic factors made little difference to its effect. However, adjusting the MDD PRS for non-genetic factors, rendered its association with depressive symptom trajectories non-significant.
We also examined the association between quartiles of the three PRS and PHQ-9 trajectories (Fig. S1). Those with the top 25% MDD PRS had 1.8-fold higher odds of being assigned to group 2, compared to those with the lowest MDD PRS (OR: 1.79, 95%CI: 1.06-3.02, P = 0.03). Those with the top 25% BMI PRS had 2.7-fold higher odds of being assigned to group 3 compared to those with the lowest BMI PRS distribution (OR: 2.68, 95%CI: 1.19-6.04, P = 0.017). Similarly, those with the top 25% BMI PRS had 1.7-fold higher odds of being assigned to group 2 compared to those with the lowest BMI PRS (OR: 1.71, 95%CI: 1.01-2.86, P = 0.04). All other findings were non-significant (P > 0.05).
Mediation analyses investigating whether IMID mediates the association between BMI PRS and depressive symptoms trajectories found no significant indirect effect of BMI PRS on depressive symptoms trajectories via IMID (Table S6). To investigate the possibility of an interaction between IMID vs. non-IMID and each PRS (Fig. S3), we performed a regression analysis stratified by the presence of an IMID or no IMID ( Table 2). The most notable finding was the increased BMI PRS effect estimate in the IMID sample (Group 3 vs. 1 = OR: 2.31, 95%CI: 1.44-3.71, P < 0.001). In contrast, BMI PRS was no longer significantly associated in the non-IMID sample (Group 3 vs. 1 = OR: 0.95, 95%CI: 0.50-1.66, P = 0.856, Table 2). Upon separately adjusting for genetic and non-genetic factors in the stratified analyses, similar findings were noted for MDD and BMI PRS compared to the full cohort, with again identifying an increased effect for BMI PRS in the IMID-only model (Table S5B).

Discussion
We aimed to understand the evolution of depressive symptoms in those with an IMID, and to determine if patient factors, including polygenic risk scores, were associated with the identified trajectories. We identified three distinct depressive symptoms trajectories, which were differentially associated with cumulative genetic burden for MDD, depressive symptoms, and BMI. Stratified regression analysis further identified that the BMI and MDD genetic burden were differentially associated with the IMID vs. the non-IMID/psychiatric sample.
Using a semi-parametric trajectory-based analysis, we stratified participants into 3 distinct trajectories based on longitudinal PHQ-9 scores. This method has been used to model longitudinal depressive symptoms in primary care [45], type 2 diabetes [56], and in the older population [31,46]. A systematic review of depressive symptom trajectories found that most studies identified 3-4 trajectories, where each varied in terms of stability and severity [30]. Similar to other studies [30,57,58], we found the largest trajectory group containing >50% of the participants belonged to a trajectory characterized by stable, minimal depressive symptoms (PHQ-9 score < 4 [59]). Other frequently reported trajectories that were also identified in the current study were that of persistently high depressive symptoms and containing <10% of the study sample [30,60], and of persistently moderate symptoms and containing 16-61% of the total sample [30,61]. We note many similarities between the trajectories identified here and the studies above, even though different depressive symptom scales were employed. Previous studies used largely either the Centre for Epidemiologic Studies Depression Scale (CES-D) or Beck Depression Inventory (BDI) [30]. Nonetheless, the PHQ-9 has been shown to be a valid screening measure for depression in the three IMIDs included here [12][13][14].
Given the lack of trajectory-based analysis using depressive symptoms specifically in those with an IMID, we looked to the general literature for comparisons. Non-genetic factors associated with trajectory group membership included smoking, and lower income or education in the current study; findings which are also seen in the broader literature [30,60,62]. Rates of lifetime depression in our study increased in a dose-dependent manner across the trajectories, with the lowest rate in the minimal depressive symptoms group (32.2%, Group 1) to >80% of the high depressive symptoms trajectory group (Group 3). Previous studies have also found that those in the higher depressive symptoms trajectory are more likely to have a formal diagnosis of psychiatric disorders, including depression [30,60].
All three polygenic risk scores investigated were associated with membership in the high depressive symptom trajectory group, when adjusting only for genetic ancestry. As expected, the polygenic risk for depressive symptoms and MDD were both associated with increasing depressive symptoms. The effect size of the PRS for all three traits were similar (OR = 1.34-1.66), with the largest effect for BMI. In addition, PRS for BMI was the only trait which remained statistically significant following the adjustment for additional potential confounders, such as age and sex. The link between BMI and depression in the general population is well-established, with a Mendelian randomization study suggesting that higher BMI partly causes depression [23,24,43]. Our stratified regression analysis further identified that the BMI PRS was significantly associated with the higher depressive symptom group in the IMID sample but was no longer significant in the non-IMID/ psychiatric sample. This may highlight a potential phenotypic difference in depressive symptoms and may point to a metabolic subtype of depressive symptoms in those with an IMID. Previous Mendelian randomization analyses found BMI is causally associated with different depressive symptoms from the PHQ-9, including tiredness, anhedonia, changes in appetite, and feelings of inadequacy [63]. Investigating subtypes of depression in those with multiple sclerosis has shown a higher risk for the somatic atypical depression symptom of weight gain [64]. IMID is unlikely to mediate the relationship between BMI and MDD given evidence from Mendelian randomization studies demonstrating IMID is not causally associated with MDD [24,27,28]. We acknowledge that differences exist in the genetics of depressive symptoms and MDD [29], and thus we undertook a mediation analysis for which also identified that IMID is unlikely to mediate the relationship between BMI and depressive symptoms trajectories. Further work in this area is needed to assess the intersection of BMI, depression, and IMIDs, and evaluate whether different therapeutic approaches for depression may be warranted in the IMID population.
This study had several strengths including the assessment of genetic and non-genetic factors in relation to longitudinal depressive symptoms in those with an IMID. The study participants had very high completion rates of the PHQ-9 for each study visit (>90%) over a 4-year period, which meant more robust modelling. The PRS were computed in our cohort using summary statistics of previous large-scale genome-wide association studies for each trait, with each study including on average, >500,000 participants. A limitation of our study is that we only utilized European genetic ancestry participants, as it is currently unadvisable to generate PRS in non-Europeans when they were originally developed in Europeans. Several studies have identified a link between racial groups and higher depressive symptoms [65,66], thus, as non-European genetic studies become more readily available, these should be considered for future genetic studies of depressive symptoms trajectories. In this cohort, we did not have access to information on the treatment of depression, which may have influenced persistence of depressive symptoms. We also did not know whether individuals with a lifetime diagnosis of MDD could be further classified as atypical MDD; a known MDD subtype associated with significant weight gain or increase in appetite and may have elucidated our finding of increased BMI PRS in the IMID sample and not in the non-IMID sample. A limitation of group-based trajectory modelling is that while comparisons can be made between trajectories, one cannot make any within-group comparisons [67].
Individuals in the persistently high depressive symptoms trajectory with an IMID are at high likelihood of experiencing clinical depression. Evidence from studies in the general population demonstrates that those experiencing a high level of persistent depressive symptoms are an important group given their association with poor outcomes [30]. Our study also identified that those with an IMID also experience chronically high depressive symptoms and may warrant greater clinical attention. Clinical relevance of our findings should be interpreted with caution as even though three stable latent trajectories were identified in this analysis, these results should not be interpreted as indicating that substantial positive and/or negative changes in mood status do not occur. Clinicians should remain vigilant for the occurrence of major mood episodes among their patients.
Future studies should assess the intersection of the three IMIDs included here along with their relationship with depressive symptoms and BMI, and include further subtyping of depression. In addition, as new GWAS' are published with larger sample sizes and different ancestries, this may improve the utility of PRS for depression in contexts whereby depression is a comorbidity, such as in those with an IMID.
In summary, the identified significant association between polygenicity for depressive symptoms, MDD, and BMI with PHQ-9 trajectory groups supports a role for genetic inheritance in the variability in depressive symptoms in those with an IMID. We have also highlighted a potential phenotypic difference in depressive symptoms in those with an IMID.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.  Characteristics of the participants included in the PHQ-9 group-based trajectories.