Interferon signaling and hypercytokinemia-related gene expression in the blood of antidepressant non-responders

Only 50% of patients with depression respond to the first antidepressant drug administered. Thus, biomarkers for prediction of antidepressant responses are needed, as predicting which patients will not respond to antidepressants can optimize selection of alternative therapies. We aimed to identify biomarkers that could predict antidepressant responsiveness using a novel data-driven approach based on statistical pattern recognition. We retrospectively divided patients with major depressive disorder into antidepressant responder and non-responder groups. Comprehensive gene expression analysis was performed using peripheral blood without narrowing the genes. We designed a classifier according to our own discrete Bayes decision rule that can handle categorical data. Nineteen genes showed differential expression in the antidepressant non-responder group (n = 15) compared to the antidepressant responder group (n = 15). In the training sample of 30 individuals, eight candidate genes had significantly altered expression according to quantitative real-time polymerase chain reaction. The expression of these genes was examined in an independent test sample of antidepressant responders (n = 22) and non-responders (n = 12). Using the discrete Bayes classifier with the HERC5, IFI6, and IFI44 genes identified in the training set yielded 85% discrimination accuracy for antidepressant responsiveness in the 34 test samples. Pathway analysis of the RNA sequencing data for antidepressant responsiveness identified that hypercytokinemia- and interferon-related genes were increased in non-responders. Disease and biofunction analysis identified changes in genes related to inflammatory and infectious diseases, including coronavirus disease. These results strongly suggest an association between antidepressant responsiveness and inflammation, which may be useful for future treatment strategies for depression.


Introduction
Depression is associated with repeated relapses and remissions. The 12-month prevalence of major depressive disorder (MDD) is approximately 6%, and the lifetime risk of depression is 15-18% [1]. In one-third of patients, depression is chronic and refractory to treatment, resulting in a significant reduction in quality of life [1,2]. Selective serotonin-reuptake inhibitors (SSRIs), which are widely used as first-line antidepressants, have a response rate of 39.6-68.0% [2] and remission rate of 23.5% [3]. The STAR*D study, a large-scale practical clinical trial, reported a response rate of 48.6% and remission rate of 36.8% in patients with depression who received an antidepressant as their first treatment step [4]. Given the low response rates, biological indicators that can predict the antidepressant response are desirable. If patients with depression who will not respond to antidepressants could be identified in advance, alternative therapies such as electroconvulsive therapy, repetitive transcranial magnetic stimulation, or cognitive behavioral therapy could be implemented.
A recent large-scale genome-wide association study (GWAS) identified three single-nucleotide polymorphisms related to antidepressant responsiveness, and a gene-enrichment analysis identified immune response genes [5]. However, other GWASs have reported that no antidepressant treatment response-related genes meeting genome-wide significance could be identified; thus, no consensus was established [6,7].
The expression of proinflammatory genes, including interleukin 1 beta (IL1B) and tumor necrosis factor alpha (TNFA), is related to treatment response [8]. A recent review describing antidepressant response candidate gene expression studies in peripheral blood revealed that mRNA levels of IL1B, IL11, TNFA, and FK506-binding protein 5 (FKBP5) may be potential predictive and mediator biomarkers [9]. In a comprehensive analysis without a prior hypothesis, the expression of genes related to inflammatory responses was increased in citalopram non-responders [10]. In another report, changes in the expression of genes involved in the immune response and inflammation were associated with antidepressant responses [11]. A study of plasma proteins related to treatment prognosis for depression indicated that cytokines, such as interleukins, are involved in treatment responsiveness [12]. However, the results of individual gene expression levels are not consistent among studies [13][14][15][16]. Therefore, further data accumulation is needed. In addition, to the best of our knowledge, only one study has reported the results of a comprehensive gene analysis validated by another method, such as quantitative real-time polymerase chain reaction (q-PCR), and in independent samples with further pathway analysis [11].
Thus, the present study aimed to comprehensively analyze gene expression in peripheral blood and identify biomarkers that could predict antidepressant treatment responsiveness using a novel data-driven approach based on statistical pattern recognition. We also aimed to examine the pathophysiology of the treatment response via pathway analysis of the identified genes.

Ethical approval
This study was conducted in accordance with the principles embodied in the latest version of the Declaration of Helsinki. The Ethics Committee of Hiroshima University (approval number H-35), the Institutional Ethics Committee of the University of Tokushima Graduate School (approval number R3-24), and the Institutional Review Board of the Yamaguchi University Hospital (approval number H30-172) approved this study. All participants provided informed consent before participation.

Participants
Participant set 1: Thirty participants were recruited from the Hiroshima University Hospital and collaborating clinics between 2012 and 2018 for the Depression Biomarker Project approved by the Ethics Committee of Hiroshima University. Patients with MDD were diagnosed by trained psychiatrists following the Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition (DSM-IV) criteria using unstructured interviews, information from medical records, and use of the Mini-International Neuropsychiatric Interview (M.I.N.I.) [17] by a research psychiatrist. All participants could provide written informed consent because the severity of depressive symptoms was moderate. No patient had taken any antidepressants for ≥1 month preceding their inclusion in the study. Participant set 2: Sixteen unmedicated patients with MDD were recruited from Tokushima University Hospital between 2004 and 2012 for the Research on Analysis of Genes for Mental Disorders study, which has been approved by the Institutional Ethics Committee of the University of Tokushima Graduate School. MDD was diagnosed based on the DSM-IV criteria by at least two trained psychiatrists. None of the patients had any other psychiatric disorders (axis I or II). The exclusion criteria were the use of non-steroidal antiinflammatory agents, steroids, antidepressants, or anticonvulsants at the time of blood collection. Participant set 3: Eighteen patients with MDD were recruited at Yamaguchi University Hospital between 2019 and 2020 using community posters for the Depression Stratification Project. This project was approved by the Institutional Review Board of Yamaguchi University Hospital. The diagnosis of MDD was based on the DSM-5 criteria by trained psychiatrists. Participants were screened using the Japanese version of the M.I.N.I. and/or by clinical interview. The exclusion criteria included current or previous substance abuse/dependence, other psychotic illnesses, any neurological disease, family history of hereditary neurological disorder, endocrine disease, head trauma, or other severe medical conditions (e.g., liver failure).
The depressive symptoms of all participants were assessed using the Structured Interview Guide for the Hamilton Depression Rating Scale (SIGH-D) [18,19] before and after treatment in accordance with the planned evaluation period (set 1: 6-8 weeks, set 2: 8 weeks, set 3: 4-12 weeks). Participants whose scores improved by ≥ 50% compared with their pre-treatment scores were defined as the responder group, and those whose scores improved by <50% were defined as the non-responder group. The demographic information of all participants is shown in Table

RNA isolation from human blood samples
Participant set 1: Total RNA was isolated from blood samples following our previous study [20]. In brief, for each participant, 7 mL of venous blood sample was collected in EDTA-Na 2 tubes and then stored in a deep freezer at − 80 • C for several years. Once thawed, the blood samples were immediately mixed with the lysis buffer of the NucleoSpin RNA blood kit (Takara Bio, Kusatsu, Japan). RNA was isolated following the protocol established by the manufacturer.
Participant sets 2 and 3: At the Tokushima University Hospital and Yamaguchi University Hospital, venous blood samples were collected from participants into PAXgene tubes (Becton, Dickinson and Company, Franklin Lakes, NJ, USA). RNA isolation from the filled PAXgene tubes was performed using the PAXgene Blood RNA Kit (Qiagen, Venlo, Netherlands) following the protocol established by the manufacturer.
The quality and quantity of RNA were measured with a NanoDrop One spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The RNA integrity numbers were measured with an Agilent Bioanalyzer with the Agilent RNA 6000 nano or pico kit (Agilent Technologies, Santa Clara, CA, USA) following the protocol established by the manufacturer.

RNA sequencing
For participant set 1, sequencing libraries were constructed employing the TruSeq Stranded Total RNA with Ribo-Zero Gold LT sample prep kit (Illumina, San Diego, CA, USA) following the protocol established by the manufacturer. After ribosomal RNA was removed, reverse transcription was performed to prepare the libraries. The libraries were pooled after quantification by bioanalyzer analysis and fluorometry with the Qubit dsDNA HS assay kit and a Qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MO, USA). Sequencing of paired-end fragments (75 bp × 2) was performed on a NextSeq 500 sequencing platform (Illumina; BGI, Beijing, China) to a depth of 13-35 (average: 21) million fragments.
For participant set 2, sequencing libraries were constructed with the TruSeq RNA Library Preparation Kit v2 (Illumina), and sequencing was performed using ṯhe HiSeq 4000 System (Illumina; BGI, Beijing, China). Data for each sample were separated to generate FASTQ files. Next-generation sequencing data were cleaned with cutadapt (version 1.8.3) [21] and cmpfastq_pe.pl (http://compbio.brc.iop.kcl.ac.uk/software/cmpfastq_pe.php). After a step of quality control, the filtered short reads were mapped to the reference genome (hg38) with STAR (version 2.5.1b) [22]. Strand-specific or non-strand-specific transcripts per million were calculated for each sample with RSEM (version 1.3.3.) [23] and were then normalized using the trimmed mean of M value method. P-values were calculated using a likelihood ratio test. Gene functional analysis was performed using an ingenuity pathway analysis (Qiagen, Hilden, Germany).

q-PCR
cDNA synthesis and q-PCR were conducted following our previous study [20]. In brief, cDNA was synthesized with a PrimeScript RT reagent kit (Takara Bio, Kusatsu, Japan) with oligo (dT) primers and 1 μg of total RNA. It was then mixed with SYBR Premix Ex TaqII (Takara Bio, Kusatsu, Japan) and specific primers. Amplification was conducted for 50 cycles with a StepOnePlus real-time PCR system (Thermo Fisher Scientific, Waltham, MO, USA). Each cycle included amplification for 15 s at 95 • C and 1 min at 60 • C. The primers we used are presented in Supplementary Table 2. All measurements were conducted in duplicate. To quantify the levels of candidate gene expression and those of Beta-2-microglobulin (B2M; the internal control), a calibration curve was created with 10, 5, 2.5, 1.25, 0.625, and 0.3125 fM synthetic DNA, which included the target sequence (Integrated DNA Technologies, Coralville, IA, USA). Fold change was calculated by the ratio of each gene to the B2M expression level of the internal control, with the mean value of the responders as 1. The synthetic DNA sequences we used are shown in Supplementary Table 3.

Statistical analysis
The distribution of participant age, SIGH-D score before starting treatment, antidepressant dose, and mRNA expression data from q-PCR were investigated with Student's t-tests, whereas the distribution of participant sex was investigated with Fisher's exact test with EZR v1.54 (http://www.jichi.ac.jp/saitama-sct/SaitamaHP.files/statmedEN.html) [24]. Holm's multi-testing correction was performed for the q-PCR results (Supplementary Table 4). The significance level was set to 0.05.

Data analysis 2.7.1. Classifier design by discrete bayes decision rule
The discrete Bayes classifier handles discrete data. x 1 and x 2 are two markers and their range is exclusively divided into divisions.
Suppose that x 1 has two divisions and x 2 has three divisions. Moreover, the discretized data of a patient belong to the first division x 1 (1) in marker x 1 and the third division . Then, P ( are defined as follows: and where n 1 1(1) denotes the number of training samples for ω 1 belonging to the division x 1(1) and n 1 2(3) denotes the number of the training samples for ω 1 belonging to the division x 2(3) Then, the class-conditional probability P ( for fx1 is given by Using n 2 1(1) and n 2 2(3) for ω 2 , we similarly get and for ω 2 is also given by The posterior probabilities of classes ω i are given as follows: In this study, we assumed that in the a priori probability, P(ω i ),P(ω 1 ) = P(ω 2 ). The patient was then assigned to the class with the maximum posterior probability. For additional details, refer to previous studies [25,26].

Marker selection
The following eight candidate genes were used to predict antidepressant response: ISG15, RSAD2, HERC5, IFIT1, IFI6, IFI44, IFI44L, and IFIT3. Among them, the combination of two markers was selected under conditions and examined. There were 30 training samples from participant set 1 and 34 test samples from participant sets 2 and 3. Note that marker selection was performed using only training samples.
We used the leave-one-out method [27] to identify an optimal combination of markers. According to the leave-one-out method, one training sample was selected as a sub-test sample from the 30 training samples and the remaining 29 training samples were assigned as sub-training samples (Supplementary Figure 1).
In pattern recognition fields, markers cannot be selected based on their individual effectiveness. Therefore, the combination of markers should be carefully selected. For this purpose, we first selected one combination of two markers. The discrete Bayes classifier was designed using 29 sub-training samples, and the re-substitution estimate for the candidate marker combination was obtained by reclassifying the 29 sub-training samples.
Next, feature criteria were calculated. Note that one sub-test sample was not used. This process was then repeated until all candidate marker combinations ( 8 C 2 = 28) had been evaluated. Among the 28 two-marker combinations, one marker combination with either maximal sensitivity subject to a specificity of ≥50% or a maximal F1 measure was selected to predict antidepressant response. Using the leave-one-out method, the above process was repeated 30 times (i.e., until each training sample had been selected once only as a sub-test sample). As a result, the most frequent combination of two markers identified in the leave-one-out loop was considered the optimal combination of markers. When the number of markers was three, marker selection was conducted again according to the same procedure.

Results
To identify genetic markers predictive of treatment response, we first conducted a genome-wide gene expression analysis without an a priori hypothesis. Among the untreated patients with MDD of participant set 1, we retrospectively assessed the treatment response after 6-8 weeks and divided them into the treatment responder (n = 15) and treatment non-responder groups (n = 15). We then performed RNA-Seq analysis using RNA obtained before starting treatment. Patient demographic data are shown in Table 1. There were no differences in sex, age, SIGH-D score before starting treatment, or antidepressant dose at the time of treatment response determination between groups. The SIGH-D score after treatment was significantly lower in the responder group than in the nonresponder group (P = 3.22E-08).
When we compared responders and non-responders, we identified 19 genes (P < 0.01 and a fold change of >1.5; Table 2). Among the 19 genes, the top 10 genes with an average count of 50 or more were validated using q-PCR. We found significant expression changes in eight genes (P < 0.05) (Fig. 1); however, OAS3 did not show significant changes according to q-PCR (P = 0.08), and Table 2 Genes with altered expression in non-responders of participant set 1 compared with those in responders based on RNA-Seq.   Table 4. Abbreviations: RES, responder; NRES, non-responder; q-PCR, qualitative real-time polymerase chain reaction.
SIGLEC1 was not successfully amplified using q-PCR. Thus, the remaining eight genes were selected as candidate predictive markers of antidepressant responsiveness. We then examined the untreated patients with MDD from an independent institution (participant set 2), whom we retrospectively divided into a treatment response group (n = 11) and treatment-resistance group (n = 5) after 8 weeks and performed comprehensive expression analysis without an a priori hypothesis using RNA-Seq data obtained before starting the treatment. The demographic data are presented in Table 1. There were no differences in sex, age, SIGH-D score before treatment initiation, or antidepressant dose at the time of treatment response determination after treatment between groups. The SIGH-D score after treatment was significantly lower in the responder group than in the non-responder group (P = 1.64E-05). The P-values and fold changes of eight candidate genes from participant set 1 are shown in Supplementary Table 5. We also identified 87 genes that met the criteria of P-values of <0.01 and fold changes of >1.5 in the responders compared with those in the non-responders according to the genome-wide gene expression analysis without an a priori hypothesis (Supplementary Table 6).
Furthermore, among the patients with MDD from another independent institute (participant set 3), we retrospectively examined the treatment response after 4-12 weeks and divided the patients into the responder (n = 11) and non-responder (n = 7) groups. The demographic data are provided in Table 1. There were no differences in sex, age, SIGH-D score before treatment, or antidepressant dose between groups. q-PCR was performed for these eight candidate genes identified from participant set 1 (Supplementary Figure 2).
Marker selection was then performed for eight candidate genes. In marker selection, feature criteria, such as accuracy, sensitivity, and F1 measure, were estimated using the discrete Bayes classifier designed with training samples from participant set 1. Varying the number of markers from two to three, all combinations of these eight genes were evaluated in the training samples. According to the feature criteria, two combinations were selected ( Table 3). The two combinations were evaluated using this discrete Bayes classifier in the independent test samples from participant sets 2 and 3 ( Table 3). The combination of HERC5 and IFI6 demonstrated high accuracy (79%). The combination of HERC5, IFI6, and IFI44 exhibited the highest accuracy (sensitivity of 91%, specificity of 75%, and accuracy of 85%).
Finally, for functional analysis, we performed pathway analysis using 19 candidate genes identified in the comprehensive analysis of participant set 1. In the ingenuity canonical pathway, interferon signaling (P = 1.38E-08) and the role of hypercytokinemia/ hyperchemokinemia in the pathogenesis of influenza (P = 4.79E-07) were significantly enhanced in non-responders relative to responders (Supplementary Table 7). Similar results were obtained using the 87 genes identified in participant set 2, an independent sample (Supplementary Table 8). These results suggested that markedly enhanced cytokine levels and interferon signaling were present in the antidepressant non-responder group.

Discussion
Immune responses and inflammation have been implicated in the pathogenesis of MDD. We performed a comprehensive gene expression analysis without an a priori hypothesis using blood samples from patients with MDD who were not yet taking antidepressants associated with treatment response. The expression levels of the identified genes were validated using q-PCR. Gene expression was further validated using independent samples by RNA sequencing (participant set 2) or q-PCR (participant set 3). Although we could not compare gene expression levels directly among samples from different facilities owing to different protocols for RNA purification and analysis, we could separate responders from non-responders with a discrimination accuracy of 85% using a combination of three genes, HERC5, IFI6, and IFI44, for independent test samples. Most importantly, we could classify the treatment response in completely independent test samples obtained from different facilities using genes identified in training samples via the data-driven approach. Our results suggested that the combination of these genes could be used as a predictive marker of treatment response. In recent depression treatment guidelines, antidepressant medication is only one option for first-line treatment, as is psychological therapy [28,29]. Prediction of potential antidepressant non-responders can be useful for a shared decision-making process. During the treatment selection phase, psychological therapy alone or in combination of psychological treatment and pharmacotherapy may be chosen, or electroconvulsive therapy may be recommended for cases of more severe depression.
We performed pathway analysis using identified gene sets and found that gene expression related to hypercytokinemia and interferon signaling was increased in non-responders. These results are in line with those of previous reports that inflammatory signals are involved in antidepressant responsiveness [8][9][10][11]. Abbreviations: RES, responder; NRES, non-responder.
The gene clusters associated with hypercytokinemia and interferon signaling were significantly increased in treatment nonresponders in three independent sample groups, and this provides a basis for speculating on the pathological mechanisms underlying antidepressant responsiveness. Interferon signaling induces depression and may be associated with increases in cytokine levels. For example, in patients with hepatitis C who were depressed after receiving interferon-alpha, increased spinal fluid IL-6 and monocyte chemotactic protein 1 levels were noted. Spinal fluid IL-6 was also negatively correlated with 5-hydroxyindole acetic acid, a serotonin metabolite [30]. Interferon-alpha has also been reported to cause reduced serum brain-derived neurotrophic factor levels, which may contribute to the mechanism underlying depression [31]. Peripheral interferon increased the expression of interferon-stimulated genes such as Rsad2, Isg15, and Ifit3 in the cortical microglia of mice and induced microglial activation [32]. In humans, microglia are considered to play crucial roles in various psychiatric disorders, including depression [33].
Regarding inflammatory cytokines and treatment responsiveness, as mentioned in the Introduction, several studies have demonstrated increased inflammatory cytokine levels in treatment-resistant depression. Therefore, the administration of anti-inflammatory agents has been proposed as a potential novel treatment for MDD. For example, depressive-like behavior and impaired neurogenesis in the hippocampus induced by interferon-alpha can be improved by minocycline administration in mice [34]. In humans, a large meta-analysis of randomized controlled trials (RCTs) of anti-inflammatory treatments reported the antidepressant effects of non-steroidal anti-inflammatory drugs (NSAIDs) and anti-inflammatory drugs [35][36][37]. Moreover, the tolerability of NSAIDs was reported to be superior to that of a placebo [37]. In contrast, large RCTs have yielded negative results on the preventive effect of small doses of aspirin on depression [38]. Thus, anti-inflammatory drugs are considered to have promising antidepressant effect; however, the selection of target patients is also important. For example, if the gene expression changes identified in this study are established as biomarkers and anti-inflammatory drugs are introduced after the identification of patients who are likely to be resistant to treatment, more effective therapeutic interventions are possible.
Notably, COVID-19 was identified in the analysis of diseases and biofunctions. It has been reported that 23.0% of COVID-19 patients had pre-existing depression [39]. Patients with psychiatric disorders, including depression, have increased mortality, hospitalization, and intensive care unit admission after severe acute respiratory coronavirus 2 infection [40]. It has also been reported that the levels of inflammatory markers, including IL-6, are increased in patients with COVID-19 [41]. The involvement of immune dysfunction in the development of depression has also been proposed [42]. Recently, OAS3, a gene identified in this study, was noted to be associated withCOVID-19 severity according to a GWAS analysis [43]. COVID-19-induced depression may be a factor affecting antidepressant resistance. In contrast, SSRIs may prevent the development of COVID-19 [44]. It will be important to monitor whether COVID-19 patients with depression respond well to antidepressants.
This study has several limitations. Firstly, the number of cases was small. However, depression is not a homogeneous disease but a heterogeneous syndrome. We previously reported that DNA methylation and gene expression profiles differ according to age of onset and sex [45][46][47][48]. Rather than merely increasing the number of participants, it is also important to collect more distinctive and uniform cases in future studies. In this study, reproducible changes in gene expression were identified in patients from three independent institutes focusing on antidepressant responsiveness. The second limitation was that this was an observational study rather than a double-blind interventional study with a placebo; hence, the effects of a placebo and nocebo on treatment response cannot be ignored. Thirdly, it was difficult to standardize the classes and dosages of antidepressants and match the timing of assessment because this study was conducted in a real-world clinical practice setting to reduce the burden on patients. However, the fact that we were able to identify reproducible treatment response markers among patients may be a strength of this study. To address these limitations, multicenter, double-blinded, large-scale interventional studies need to be conducted. The fourth limitation is that the methods of RNA purification and expression analysis differed among institutions. In our previous study, when RNA-Seq was conducted with RNA purified using different methods, significant differences in housekeeping gene expressions were identified [20]. It is necessary to standardize RNA purification methods and establish discrimination criteria at each institution for clinical application as predictive markers of treatment responses. The final limitation is that only RNA analysis was performed. A more robust biomarker could be identified by setting up a complex biomarker with plasma or DNA methylation.
Despite these limitations, the results of this study form the basis for the development of predictive markers of the therapeutic response to depression and may provide insights into the development of novel therapeutic strategies for depression, such as antiinflammatory drugs.