Association of psychosocial factors and biological pathways identified from rare-variant analysis with longitudinal trajectories of treatment response in major depressive disorder

Abstract

Background

Antidepressant efficacy is influenced by a multitude of factors, yet predicting treatment outcomes remains challenging. This difficulty is partly due to the commonly employed dichotomous classifications of treatment response that rely on a single primary endpoint.

Methods

The study enrolled 972 patients diagnosed with depression, including both first-episode and recurrent cases. All patients received treatment with a single class of antidepressant medication over an eight-week period. Treatment response trajectories were identified through cluster analysis using normalized score change ratios from the 17-item Hamilton Rating Scale for Depression (HAMD-17) at baseline and weeks 2, 4, 6, and 8. The impact of psychosocial factors—including childhood trauma experience, social support, and family environment—on these response patterns was evaluated using ANOVA and Tukey’s HSD tests. Additionally, targeted exome sequencing was conducted to perform rare-variant burden and enrichment analyses to investigate genetic influences on antidepressant response.

Results

Three patterns of antidepressant treatment response were identified: gradual response (C1 cluster), early response (C2 cluster), and fluctuating response (C3 cluster). Notably, patients in the C3 cluster exhibited higher levels of suicidal ideation, alexithymia, and anhedonia after the treatment period, along with the highest baseline levels of family control (a subscale of the family environment). Our rare-variant analysis revealed genes associated with response efficiency between C1 and C2 clusters to be significantly enriched in the neurotrophin signaling pathway (odds ratio = 23.94; p-adjusted = 6.96e-05). In addition, genes linked to response volatility between C1 and C3 clusters were enriched in the regulation of inflammatory mediators of transient receptor potential (TRP) channels (odds ratio = 31.5; p-adjusted = 1.83e-07).

Conclusions

Our findings suggest that patients exhibiting a fluctuating response to antidepressant treatment may endure more severe clinical symptoms throughout the treatment course. The involvement of the neurotrophin signaling pathway and TRP channels in these response patterns highlights their potential as novel targets for therapeutic intervention in depression. This underscores the importance of personalized treatment strategies that consider the underlying genetic and psychological factors influencing antidepressant efficacy.

Background

Major depressive disorder (MDD) is a prevalent mental illness affecting over 300 million individuals worldwide [1]. Although antidepressant medications are the primary treatment for MDD, the significant variability in treatment responses is a major challenge. Approximately 70% of patients do not respond to a single trial of antidepressants, and 20–30% remain unresponsive even after multiple interventions [2]. The prevailing'trial and error'strategy employed by clinicians to identify the most suitable treatment underscores the pressing necessity for more precise and individualized therapeutic interventions.

Existing literature recognizes the involvement of psychosocial, clinical, and genetic factors in influencing antidepressant efficacy [3]. Psychosocial factors as predictors of treatment response are frequently inconsistent or insufficient for clinical application. For instance, while certain studies have demonstrated that higher levels of social support at baseline predict better antidepressant response, others have failed to establish a significant association between antidepressant response and social support [4, 5]. Additionally, specific adverse clinical symptoms, such as suicidal ideation, anhedonia, and alexithymia, have demonstrated correlations with the outcomes of antidepressant treatment [6,7,8].

Genetic variation is a promising predictor of antidepressant efficacy. Several genome-wide association studies (GWASs) focusing on common variants have been conducted. Nine of the published articles have revealed SNP signals associated with antidepressant efficacy of genome-wide significance [9,10,11,12,13,14,15,16,17]. Except for rs116692768, which was weakly replicated after adding the exome genotype to previously available genome-wide data, there were no consistent findings, even for biological pathways [17]. Notably, the largest GWAS integrating genetic data from a substantial MDD patient cohort failed to identify significant variants associated with remission or percentage improvement [18]. Rare variants are of particular interest, as they may be directly associated with functional disruptions and have a higher likelihood of being causal compared to common variants [19]. Investigating rare variants could provide valuable insights into biological mechanisms that are less detectable through studies of common variants, thereby enhancing the clinical applicability of genetic findings in precision psychiatry. Limited studies have focused on the impact of rare variants on antidepressant efficacy. Kang and his colleagues found that genetic indicators for non-remission in MDD patients were rare variants [20]. Two exome analyses of treatment-resistant depression (TRD) patients suggested some biological mechanisms associated with TRD [21, 22]. In addition, a small rare variant analysis identified 35 genes associated with antidepressant response [23]. The absence of robust findings, coupled with the limited number of studies conducted in this field underscores the need for a more comprehensive investigation.

Emerging evidence also suggests that immune and inflammatory processes play a significant role in the pathophysiology of depression and treatment response [24]. Elevated levels of pro-inflammatory cytokines, such as interleukin-6 (IL-6) and tumor necrosis factor-alpha (TNF-α), have been associated with poor response to conventional antidepressants [25]. A subset of depressed patients demonstrates heightened inflammatory activation, which may contribute to resistance to standard treatment approaches [26].

The poor achievement in predicting antidepressant efficacy may be attributed to the inappropriate classification of outcomes. In MDD trials, results are often defined using a predetermined cutoff score measured at a single primary endpoint [27]. The definition ignores rich information in the process of symptom change over time and reduces the potential statistical power of measurements at multiple time points [28]. Identifying distinct trajectories of symptom change can provide a foundation for improving the design of antidepressant clinical trials [29]. We attempted to use the change rate of the HAMD-17 scores to directly reflect the individual treatment response to explore alternative categorization of outcomes [30].

In this study, we examined the patterns of treatment response in MDD patients using multi-point longitudinal data. We focused on the role of rare variants in antidepressant efficacy by targeted exome sequencing. In addition, we also examined differences in psychosocial factors and symptom characteristics between different patterns.

Methods

Study design overview

The overall study design flow is demonstrated in Fig. 1. This study of antidepressant efficacy was carried out in a naturalistic clinical setting. It predominantly involved untreated adult inpatients having current episodes of nonpsychotic unipolar depression.

Fig. 1

Study design flow chart. HAMD-17, 17-item Hamilton depression rating scale; CTQ-SF, childhood trauma questionnaire-short form; FES-CV, family environment scale-Chinese version; SSRS, social support rating scale; BSI-CV, Beck scale for suicide ideation-Chinese version; SHAPS: Snaith-Hamilton pleasure scale; TAS-20: 20-item Toronto alexithymia scale; ANOVA, Analysis of Variance; Tukey HSD, Tukey honestly significant difference

Participants

This study was conducted from 2010. Participants enrolled in this study were Chinese Han patients recruited from the Zhongda Hospital MDD inpatient database. Patients were diagnosed by two independent senior psychiatrists, and if diagnoses differed, the final diagnosis was confirmed by a third psychiatrist who was blinded to previous evaluations. Clinical assessments, including psychiatric evaluations, administration of the HAMD-17, and psychosocial scales, were conducted by trained psychiatrists at Zhongda Hospital. The exclusion criteria included other documented diagnostic history on Axis 1 of the DSM-IV, independent diagnosis of intellectual disability, pregnancy, lactation, personality disorder, primary organic disease(e.g., major neurocognitive disorders such as Alzheimer’s disease, severe traumatic brain injury, epilepsy, multiple sclerosis, or other progressive neurological diseases), and other medical conditions that might significantly impair psychiatric assessment or confound the diagnosis of MDD. Additionally, a history of electroconvulsive therapy within the last six months, contraindication to repetitive transcranial magnetic stimulation (rTMS), and experiencing a manic episode within 12 months post-entry led to exclusion. The study was approved by the hospital ethical committee (2016ZDSYLL100-P01), and all participants signed written informed consent. In total, 972 MDD patients remained for the downstream analysis.

Antidepressant treatment

All MDD patients received treatment with a single antidepressant medication initiated at the minimum effective doses recommended by standard clinical guidelines. Among these patients, 61.2% received selective serotonin reuptake inhibitors (SSRIs), 31.3% serotonin and noradrenaline reuptake inhibitors (SNRIs), and 7.5% other types of antidepressants. Dosages were adjusted as clinically indicated based on patient response and tolerability throughout the eight-week treatment period. A subset of patients also received concurrent repetitive transcranial magnetic stimulation (rTMS) during the first two weeks (14 consecutive days). Concomitant psychotropic medications were prohibited except for low-dose benzodiazepine anxiolytics (e.g., alprazolam 0.4–0.8 mg/day or estazolam 1–2 mg/day) for insomnia when clinically necessary. Patients who altered their antidepressant medication or demonstrated significant non-compliance were excluded from the study.

Assessment of treatment outcome, psychosocial factors, and clinical symptoms

The severity of depressive symptoms was measured by the HAMD-17 at baseline and at scheduled follow-up time points at weeks 2, 4, 6, and 8. We evaluated several psychosocial factors that might have influenced the treatment response at baseline, including childhood trauma experience, family environment, and social support. Childhood trauma was measured using the Childhood Trauma Questionnaire-Short Form (CTQ-SF), which covered five areas: sexual abuse, physical abuse, physical neglect, emotional abuse, and emotional neglect [31, 32]. We used all the ten dimensions of family environment in the family environment scale-Chinese version (FES-CV), including cohesion, expressiveness, conflict, independence, achievement orientation, intellectual cultural orientation, active recreational orientation, moral-religious emphasis, organization, and controls [33, 34]. Social support was analyzed through the Social Support Rating Scale (SSRS), which encompassed objective and subjective support scores, utility of support, and a total social support score [35, 36].

We also evaluated several clinical symptoms, including suicidal ideation, alexithymia, and anhedonia at baseline and week 8 post-treatment. Suicidal ideation was measured using the Beck scale for suicide ideation-Chinese version (BSI-CV), where higher scores represented more severe ideation [37, 38]. Alexithymia, defined by impairments in understanding and managing emotions [39], was measured using the 20-item Toronto Alexithymia Scale (TAS-20). The scale's total scores ranged from 20 to 100, where higher values signify more severe alexithymia [40, 41]. Anhedonia was evaluated using the Snaith-Hamilton pleasure scale (SHAPS). It comprised 14 items rated on a 4-point Likert scale, with higher scores reflecting higher levels of anhedonia [42, 43].

Targeted exome sequencing

To study the rare variants associated with antidepressant response, we conducted targeted exome sequencing for 1309 genes (Supplementary Table S1). The targeted genes were chosen based on pathways reported to be associated with MDD and/or with regulating antidepressant efficacy. They were sequenced using the Illumina MiSeq high-throughput sequencing platform (Illumina, San Diego, CA, USA). The data were processed with Illumina's Bcl2fastqv2.16.0.10 software for demultiplexing into individual Fastq files, followed by quality trimming of low-quality reads and bases (Q < 15) using the FastX tool. High-quality reads were aligned to the NCBI human reference genome (hg19) using Burrows-Wheeler Aligner (BWA). Post-alignment processing included duplicate marking with Picard's MarkDuplicates, and realignment and recalibration of base quality scores using SAMtools and GATK. To identify single nucleotide variants (SNVs), criteria of minimal depth coverage > 20 × and a quality score > 30 in over 80% of sequenced subjects were applied. SNVs with a minor allele frequency (MAF) of < 1% were considered rare variants. After thorough in-house quality control of clinical and sequencing data, 56,552 rare SNVs within the 1,309 targeted genes were successfully annotated.

Statistical analysis

Treatment outcomes evaluation

To assess treatment outcomes, we employed two approaches utilizing multi-time point longitudinal data measured by HAMD-17. Initially, we adhered to traditional measures of antidepressant efficacy. Initially, we adhered to traditional definitions of antidepressant efficacy, categorizing outcomes at each assessment time point (weeks 2, 4, 6, and 8 post-treatment) into remission (HAMD-17 score < 7) and non-remission (HAMD-17 score ≥ 7), as well as response (≥ 50% reduction in HAMD-17 score from baseline) and non-response (< 50% reduction from baseline). To enrich our analysis with temporal data, we constructed decision trees to classify patients based on their ability to maintain remission or response across four time points. To enrich our analysis with longitudinal response patterns, we constructed decision trees to classify patients according to their ability to maintain treatment efficacy across the entire follow-up period (weeks 2, 4, 6, and 8). Specifically, patients were categorized as follows: Stable remission (maintained HAMD-17 scores consistently below 7 across all subsequent time points); Unstable remission (achieved remission at one time point but failed to consistently maintain remission); Consistently not-remitted (never achieved remission at any assessed time points); Stable response (maintained response consistently after initial achievement); Unstable response (or fluctuating response) (achieved response at one time point but failed to consistently maintain response); Consistently non-responsive (never achieved response criteria throughout the entire assessment period). Furthermore, we implemented a time course analysis using the “timeclust” function from the TCseq R package (TCseq: Time course sequencing data analysis. R package version 1.20.0), which was designed for the analysis of longitudinal data, and offered advanced tools for differential analysis across various time points as well as for the identification and visualization of temporal patterns in complex data sets [44]. We calculated the change ratio of HAMD-17 scores at weeks 2, 4, 6, and 8 by subtracting the score at each time point from the baseline score, then dividing by the baseline score. These change ratios, normalized to fit a normal distribution, served as inputs for the “timeclust” function. For the timeclust function, we utilized the default parameters, which included the k-means clustering method with a maximum of 10 random starts and a limit of 100 iterations for convergence. To control for the potential confounding effects of concurrent rTMS treatment, we performed a separate analysis excluding these patients, ensuring a more accurate assessment of antidepressant response trajectories.

Psychosocial factors and clinical symptoms comparison among the clusters

To examine the influence of Psychosocial factors on treatment clusters and explore variations in clinical symptoms among them, we utilized Analysis of Variance (ANOVA) to identify group differences. For analyses with significant ANOVA results (p < 0.05), we applied a False Discovery Rate (FDR) correction to account for multiple testing and control the false positive rate. Post hoc pairwise comparisons were subsequently conducted using Tukey's Honestly Significant Difference (HSD) test only for variables that remained significant after FDR correction (adjusted p < 0.05). To evaluate symptom improvement following treatment, we calculated change scores by subtracting baseline scores from those at week 8.

Rare-variants burden analysis and pathway enrichment analysis

To identify the genes that show different rare variants burden among different response clusters, rare-variants burden analysis was conducted using rare variant test software for next-generation sequencing data (Rvtests). Rvtests was developed as a comprehensive tool to support genetic association analysis and meta-analysis [45]. We applied the combined multivariate and collapsing (CMC) method within Rvtests to conduct association tests for rare variants [46]. The CMC method is a type of direct burden test that collapses variants into a single score. Our assumption is that most rare variants in the target regions contribute to the phenotype in a consistent direction. While methods like SKAT allow for varied effect directions, CMC provides greater power under the expectation of a high proportion of uniformly causal variants, making it more suitable for our analysis framework. Prior to analysis, covariates such as sex, age, and genetic principal components were adjusted to reduce confounding influences. We used genes as the unit to group variants based on the annotation of the human genome version 19. Different burden analyses for different groups of clusters were calculated for each tested gene. The FDR approach was applied to correct for multiple comparisons and consider the genes with FDR < 0.05 as response-associated genes. If there are no genes that survive the multiple testing, we also consider genes with p-value < 0.05.

For the response-associated genes, we then performed enrichment analysis using Enrichr to get the functional pathway based on annotated gene sets representing prior biological knowledge. The analysis employed the default settings provided by the platform, using a p-value cutoff of 0.05 for statistical significance and applying the Benjamini–Hochberg method for multiple testing correction. To capture a wide range of biological pathways, we utilized extensive gene set libraries available in Enrichr, including KEGG, Reactome, and others. Using the response-associated genes as an input set and the tested genes as a background, we checked whether response-associated genes significantly overlapped with annotated gene sets.

Results

Classification of antidepressant treatment outcomes: traditional definitions

During the eight-week follow-up period, 27 patients were lost to follow-up due to inability to contact or refusal to participate in follow-up visits, resulting in incomplete HAMD-17 data at the endpoint. We conducted a comparative analysis that included demographic variables (gender, age, educational level) and baseline HAMD-17 scores. There were no significant differences between the two groups in gender, age, or educational level (p > 0.05). However, patients who dropped out had significantly higher baseline HAMD-17 scores than those who completed the study (p = 0.00046) (Supplementary Table S2).

Using decision-tree analysis based on longitudinal HAMD-17 scores, we identified distinct groups of antidepressant treatment outcomes. Regarding remission, 269 patients showed consistent non-remission, never meeting remission criteria throughout the follow-up period; 609 patients achieved sustained remission, reaching remission criteria at a certain point and maintaining it until week 8; and 94 patients exhibited non-sustained remission, initially achieving remission but failing to maintain it until week 8 (Supplementary Fig S1). Similarly, in terms of treatment response, 44 patients demonstrated consistent non-response, never reaching the ≥ 50% improvement threshold; 859 patients showed sustained response, achieving and maintaining ≥ 50% improvement until week 8; and 69 patients showed non-sustained response, initially meeting response criteria but failing to sustain it until the end of week 8 (Supplementary Fig S2). Notably, 294 of the 972 patients fell into different classifications depending on whether remission or response criteria were applied. Specifically, 187 patients met response criteria (≥ 50% improvement) without ever reaching remission (HAMD-17 < 7), primarily because individuals with high baseline severity could demonstrate substantial improvement but still fail to meet absolute remission criteria. Conversely, patients with lower baseline severity might achieve remission even if their overall percentage improvement was relatively modest (Fig. 2A). This discrepancy underscores the limitations of using either remission or response criteria in isolation, as they may not fully capture the complexity of individual treatment trajectories. To further illustrate this issue, representative HAMD-17 trajectories from individual patients across different categories are shown in Fig. 2B, demonstrating substantial heterogeneity in symptom progression even among patients classified similarly under single-dimension outcome measures.

Fig. 2

A: Assignment of patients based on classifications for both response and remission derived from decision-tree analysis. B: Representative individual examples illustrating distinct antidepressant response trajectories identified within patient clusters. Each solid line represents one patient's change in HAMD-17 scores over the treatment period. The three dashed horizontal lines indicate clinically relevant cutoff points on the HAMD-17 scale. Red line: HAMD-17 score of 24, differentiating severe depressive states; Brown line: HAMD-17 score of 17, distinguishing moderate depressive severity; Green line: HAMD-17 score of 7, indicating clinical remission

Clusters of treatment response trajectories

To represent patients’ response trajectories over time, the patients were grouped into three clusters based on their HAMD-17 score change ratio over time (Fig. 3). The first cluster (C1 cluster) was a slow response cluster, which showed a continuing increase in their change ratio in the HAMD-17 score. We defined it as the gradual response group. There were 402 patients in the C1 cluster, which represented 41.35% of the patients. The second cluster (C2 cluster) was an early response cluster, which showed an increase in the HAMD-17 score change ratio in the first four weeks and stayed well for the rest of the four weeks. We defined it as the early response group. There were 449 patients in the C2 cluster, which represented 46.19% of the patients. The third cluster (C3 cluster) was diverse, showing a fluctuating course of alternating improvement and worsening and ultimately poor treatment response. We defined it as the fluctuating response group. There were 121 patients in the C3 cluster, which represented 12.56% of the patients. Table 1 summarizes the demographic data of the patients classified into three antidepressant response clusters. The distribution of antidepressants among the three clusters was assessed using a Chi-square test. The results indicated no significant differences in medication use between the clusters (χ² = 5.456, p = 0.204) (Supplementary Table S3). Patients who only had pharmacotherapy (n = 579) showed similar response trajectories and patterns (Supplementary Fig S3). To assess the potential impact of missing data on our results, we applied multiple imputation to generate five different datasets and conducted clustering analysis on each imputed dataset. The results showed that unsupervised clustering consistently grouped patients into three clusters across all imputed datasets, indicating that the missing data had a relatively small impact on the clustering results and did not alter the main findings of our study.

Fig. 3

HAMD-17 score change rate trajectories and resulting cluster shape characteristics for all patients. X-axis: observation time in weeks; Y-axis: normalized HAMD-17 score change ratio; membership: each sample has a membership value ranging from 0 to 1, which indicates its degree of belonging to the cluster

Table 1 Demographic data among different treatment response clusters

Differences in the family environment and social support among the clusters

We found significant differences in family control and social support but not in childhood trauma experience among the three clusters. The differences in family environment factor-control between the three clusters were significant [1–4], C2 = 3 [1–4], C3 = 3 [1.75–5.25]; p = 0.0489) (Fig. 4a). The control factor reflected the extent to which the family is organized hierarchically and the rigidity of family rules and procedures [33]. Patients in the C3 cluster had significantly higher levels of control compared with those in the C1 (p = 0.0096) and C2 (p = 0.031). For social support, there were significant differences in support utility ([4–8], C2 = 5 [4–8], C3 = 5 [4–7]; p = 0.0339) among the three clusters (Fig. 4b). Patients in the C1 cluster had significantly higher levels of support utility than those in C3 (p = 0.0075). There were no significant differences between C1 and C2 (p = 0.17), or between C2 and C3 (p = 0.16) in support utility.

Fig. 4

Violin plots of clinical psychological scale scores for the three clusters of patients. ANOVA analysis and Tukey HSD test of clinical psychological scale scores among three clusters were conducted. Figures above each two violin plots are p-values by Tukey’s HSD test. The numbers 1, 2, and 3 on the X-axis represent the C1, C2, and C3 cluster, respectively. a. Family control: Distribution of the control subscale score from the Family Environment Scale among clusters; b. Support utility: Distribution of support utility scores among clusters; c. BSI at baseline: Distribution of baseline BSI scores across clusters; d. BSI at week 8: Distribution of Beck Scale for Suicide Ideation scores at week 8 among clusters. e. SHAPS at week 8: Distribution of Snaith-Hamilton Pleasure Scale scores at week 8; f. Change in SHAPS (week 8 - baseline): Change scores for SHAPS from baseline to week 8; g. TAS-20 total at week 8: Distribution of Toronto Alexithymia Scale total scores at week 8; h. Change in TAS-20 (week 8 - baseline): Change scores for TAS-20 total scores from baseline to week 8

Differences in clinical symptoms among the clusters

There were significant differences in suicidal ideation among the three clusters at baseline (median [IQR]: C1 = 2 [0–9], C2 = 3 [0–11], C3 = 7 [0–13]; p = 0.00161) and at week 8 (median [IQR]: C1 = 0 [0–0], C2 = 0 [0–0], C3 = 0 [0–0]; p = 1.88e-06) (Fig. 4c, d). Patients in the C1 cluster had significantly lower scores at baseline (C1 vs. C2, p = 0.029; C1 vs. C3, p = 0.00038). Patients in the C3 cluster had significantly higher scores than those in C1 cluster (p = 2.3e-06) and C2 cluster (p = 3.0e-07) at week 8. For SHAPS, there were significant differences in [10–13], C2 = 11 [9–13], C3 = 9 [5–11]; p < 2.2e-16) and change scores from baseline to week 8 (median [IQR]: C1 = 5 [2–8], C2 = 5 [2–9], C3 = 1 [−1.25–4]; p = 1.85e-15) (Fig. 4e, f). Patients in the C3 cluster showed significantly higher levels of anhedonia compared with those in the C1 cluster (p < 2.2e-16) and C2 cluster (p < 2.2e-16) and fewer change scores of anhedonia from baseline to week 8 (C1 vs. C3, p = 5.3e-14; C2 vs. C3, p = 4.7e-14). For TAS, there were significant differences in total TAS scores at week 8 (median [IQR]: C1 = 50 [47–56], C2 = 50 [47–56], C3 = 55 [51.5–60]; p = 5.8e-11) and change scores from baseline to week 8(median [IQR]: C1 = −9 [−14– −3], C2 = −9 [−15– −4], C3 = −4 [−11–2]; p = 6.65e-06) among the three clusters (Fig. 4g, h). Patients in the C3 cluster had a significantly higher score in TAS at week 8 (C1 vs. C3, p = 1.3e-09; C2 vs. C3, p = 4.8e-12) and fewer change scores from baseline to week 8 (C1 vs. C3, p = 5.7e-05; C2 vs. C3, p = 6.0e-07).

Rare variance burden difference between clusters

To estimate the rare variance difference between different treatment response clusters, we performed rare variance burden tests between the C1 and C2 clusters, as well as between the C1 and C3 clusters. The C1 cluster was used as the reference. No gene was found with FDR < 0.05. We found 50 genes (p < 0.05) with rare variance burden differences between patients in the C1 and C2 clusters (Supplementary Table S4), and 64 genes (p < 0.05) with rare variance burden differences between patients in the C1 and C3 clusters (Supplementary Table S5). The top ten genes identified in the burden tests are presented in Table 2. The genes with burden differences between the C1 and C2 clusters represent response efficiency-related genes, while those with burden differences between the C1 and C3 clusters represent response volatility-related genes.

Table 2 Top ten genes in the burden tests

Functional enrichment of genes with different rare variants between clusters

To illustrate the functional information of the rare variance burden genes, we performed functional enrichment of the treatment response-related rare variant genes and found the response efficiency-related (C1 vs. C2) genes were enriched in neuronal-related pathways. For example, the neurotrophin signaling pathway (Fig. 5A). The response volatility-related (C1 vs. C3) genes were enriched in inflammatory-related pathways, for example, inflammatory mediator regulation of transient receptor potential (TRP) channels (Fig. 5B). P-values and the genes involved are listed in Supplementary Table S6-7.

Fig. 5

The top ten enriched KEGG pathways for the treatment response-related rare variant genes. A: The top ten enriched KEGG pathways of response efficiency-related genes. B: The top ten enriched KEGG pathways of response volatility-related genes

Discussion

Previous studies for predicting antidepressant efficacy have been largely unsuccessful due to the limited definition of treatment response using a single time point. In this study, we used longitudinal data from multi-point measurements and performed a cluster-based method to determine treatment response patterns. These treatment response clusters showed differences in clinical symptoms, psychosocial factors, and relevant biological pathways. The results of functional enrichment of rare variants shed light on the underlying biological pathways related to different patterns of responses.

Group-based trajectory methods can incorporate more information to capture heterogeneity in efficacy that traditional approaches would overlook, and improve the analysis of clinical studies [47]. A study using latent-class trajectory analysis provided more consistent predictors of antidepressant response than traditional endpoint analyses [48]. Our categorization of outcome truly reflected the situation of clinical treatment, so further analysis based on the classification can be more advantageous in determining the components that influence antidepressant efficacy.

The results of functional enrichment revealed some meaningful findings. The response efficiency-related (C1 vs. C2) genes were mainly concentrated in the neurotrophin signaling pathway, which suggested a connection with the rapid onset of antidepressant action. Brain-derived neurotrophic factor (BDNF), a major neurotrophic factor in the central nervous system (CNS), exerts its neurotrophic effects mainly through the activation of tropomyosin-related kinase receptor B (TrkB) [49]. Increased BDNF expression and signal transduction may play a key role in the rapid activation mechanism of antidepressants [50, 51]. Currently, the direct evidence supporting these mechanisms primarily comes from animal experiments. A study in rodents indicated that impaired processing and release of BDNF significantly attenuate rapid antidepressant responses to treatments like scopolamine [50]. From a clinical perspective, we strengthen the evidence that neurotrophic factors and their signal transduction may be involved in the mechanism of the rapid onset of antidepressant action.

The enrichment analysis of response volatility-related genes (C1 vs. C3) revealed a connection with inflammation, particularly inflammatory mediator regulation of TRP channels. TRP channels are cation channels that influence the functional state of Ca²⁺ signaling in various immune cells and can thus modulate immune and inflammatory responses [52]. Elevated levels of proinflammatory cytokines have been observed in both the cerebrospinal fluid and peripheral blood of depressed patients, suggesting that inflammatory dysregulation may contribute to treatment resistance or inadequate therapeutic response [53]. Notably, TRP channels regulate not only the expression of peripheral proinflammatory genes but also the production of proinflammatory molecules in the CNS by microglia, further linking their activity to inflammation-related mechanisms in depression [54, 55]. Animal studies specifically indicate that activation of TRP channels such as TRPV1 and TRPA1 enhances central inflammation and depressive-like behaviors. Conversely, inhibition or genetic deletion of these TRP channels mitigates inflammation-induced depressive behaviors, highlighting their promise as novel antidepressant targets [56,57,58]. However, direct clinical validation of these preclinical mechanisms remains scarce. Our findings thus provide additional clinical-level support for TRP channels as promising therapeutic targets in depression.

The analysis of family environment and social support among the three groups yielded noteworthy results. Patients in the fluctuating response group (C3) showed the highest control levels, implying that higher control levels may be detrimental to antidepressant treatment. Previous research has established a strong association between family control and depression symptoms [59]. However, there is currently no direct evidence demonstrating that reducing family control through family therapy improves antidepressant treatment outcomes. Our findings suggest a possible relationship between high family control and a fluctuating response pattern, but further research is needed to determine the causal nature of this association. Therefore, investigating whether interventions targeting family dynamics, such as family therapy, could contribute to improved treatment outcomes in MDD warrants further exploration. The result of the comparison of social support between the gradual response group (C1) and the fluctuating group (C3) aligns with previous research that has found that greater baseline social support is linked to better antidepressant treatment outcomes [60]. Clinical symptoms were also examined among the three groups. We discovered that the recovery of alexithymia was less pronounced in the fluctuating response group. Previous studies have reported that alexithymia can worsen the severity of symptoms and negatively affect antidepressant treatment in MDD patients [39]. Psychosocial interventions such as psychological group therapy have been shown to alleviate alexithymia [61, 62]. Therefore, for patients with high alexithymia scores and large efficacy fluctuations, strengthening psychological treatment may help improve alexithymia and ultimately lead to better treatment response. In addition, patients in the fluctuating response group demonstrated higher suicidal ideation and less improvement in anhedonia, supporting a less successful therapy because anhedonia and suicidal ideation are core symptoms of depression [63].

Some limitations are worth noting. First, although our results are encouraging, they should be validated in an independent cohort with larger sample sizes. Second, this eight-week follow-up longitudinal study has provided valuable insight into the antidepressant trajectories; however, future studies should be expanded to include a longer follow-up period to examine the long-term trajectories of depressive symptoms in patients with MDD after clinical treatment. Third, our study was a naturalistic therapy study that included various antidepressant categories. Further research using mono-therapy with a larger sample size would benefit the underlying mechanism study. Additionally, it is important to note that patients who dropped out of the study had significantly higher baseline HAMD-17 scores than those who completed the study. This pattern of missing data suggests that data were not Missing Completely at Random (MCAR) but were instead conditionally missing in relation to baseline depression severity. Specifically, patients with higher baseline HAMD-17 scores exhibited a higher likelihood of missing follow-up assessments, indicating a Missing at Random (MAR) mechanism. This association may introduce bias into the interpretation of the primary outcome, and therefore, the results should be viewed with caution. Future studies should take into account potential dropout bias and consider methods such as imputation to address missing data.

Conclusions

Our results suggest that social support and family control may affect depression treatment response. Patients exhibiting a fluctuating response to antidepressant treatment may endure more severe clinical symptoms throughout the treatment course. We provide novel evidence for inflammatory mechanisms, and specifically the TRP pathway, as potential therapeutic targets for depression, as well as an association between neurotrophin signaling and mechanisms of rapid onset of antidepressant action.