Systems Approach to Identify Common Genes and Pathways Associated with Response to Selective Serotonin Reuptake Inhibitors and Major Depression Risk

Despite numerous studies on major depressive disorder (MDD) susceptibility, the precise underlying molecular mechanism has not been elucidated which restricts the development of etiology-based disease-modifying drug. Major depressive disorder treatment is still symptomatic and is the leading cause of (~30%) failure of the current antidepressant therapy. Here we comprehended the probable genes and pathways commonly associated with antidepressant response and MDD. A systematic review was conducted, and candidate genes/pathways associated with antidepressant response and MDD were identified using an integrative genetics approach. Initially, single nucleotide polymorphisms (SNPs)/genes found to be significantly associated with antidepressant response were systematically reviewed and retrieved from the candidate studies and genome-wide association studies (GWAS). Also, significant variations concerning MDD susceptibility were extracted from GWAS only. We found 245 (Set A) and 800 (Set B) significantly associated genes with antidepressant response and MDD, respectively. Further, gene set enrichment analysis revealed the top five co-occurring molecular pathways (p ≤ 0.05) among the two sets of genes: Cushing syndrome, Axon guidance, cAMP signaling pathway, Insulin secretion, and Glutamatergic synapse, wherein all show a very close relation to synaptic plasticity. Integrative analyses of candidate gene and genome-wide association studies would enable us to investigate the putative targets for the development of disease etiology-based antidepressant that might be more promising than current ones.


Introduction
Major depressive disorder (MDD) is the third largest cause of burden of disease and it is responsible for almost 80% of psychiatric hospitalizations. According to recently conducted world mental health surveys, MDD is experienced by 10-15% people in their lifetime [1] and it can lead to high incidence of suicide. Over 800,000 lives are lost yearly due to suicide, which translates to 3000 suicide deaths every day [2]. It is almost more than half a century now since the first antidepressant drug was discovered, starting from non-specific monoamine oxidase inhibitors (MAO-I) and tricyclic antidepressants (TCA) to target specific selective serotonin reuptake inhibitors (SSRIs). Among them, SSRIs have proven to be the most effective drugs to date, yet approximately 30-40% of depressive patients do not or partially respond to the therapy whereas 60-75% fail to achieve complete remission [3]. This may be attributed to the poor understanding of MDD pathophysiology and lack of etiology-based drugs. Candidate gene-based studies and GWAS have shown that the clinical heterogeneity in therapeutic outcome is also influenced by a variety of genetic (single nucleotide polymorphisms, SNP; copy number variations, CNV; insertions, I; deletions, D etc.), pathophysiological and environmental factors [4,5].
With the advent of high throughput technologies it became possible to generate numerous genomic datasets to identify genetic markers associated with complex disorders and predict drug response. However, these datasets are not consistent enough to turn these findings into clinical practice. The inconsistency may be attributed to the difference in ethnicity of the studied population, underpowered study design, and various confounding factors like age of onset, the severity of disease, gender, etc. Therefore, there is an impending need for integrating these datasets so that we can have a more advance understanding of disease etiology as well as drug response. In pursuit of a better understanding of the problem, several integrated approaches have been put forward for studying the interactions between disease-associated genes and proteins [6]. Network and pathway analysis of candidate genes involved in MDD has provided important information about gene interaction and regulation in MDD [7]. A very recent study has explored the common link for pathogenesis between MDD and glioblastoma using the transcriptomics convergence method, thus providing a new approach to analyze the available huge genetic datasets [8].
With the genetic data on SSRI response and MDD available across diverse populations, a systematic review will help us to better interpret (or curate) the findings of these studies. Therefore, in this study, we first systematically reviewed the literature concerning the genetics of SSRI response (studies on responder versus non-responder patients on SSRIs) and MDD (studies on MDD cases versus healthy controls). Further, we used an integrative genetics approach and retrieved common genes and pathways involved in SSRI response and MDD manifestation to find out potential drug targets for etiology-specific antidepressants development. We identified 29 overlapping genes from a systematic literature search concerning SSRI response and MDD development. Finally, we performed functional enrichment analysis using the WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) [9] (http://webgestalt.org/) to identify the top pathways that are inter-linked between antidepressants response outcome as well as MDD, on the molecular basis. Hence, by integrating genes associated with disease and drug response studies we found a list of genes and pathways which can be validated and might be novel molecular targets for etiology-based antidepressant development.

Results
A systematic review of candidate gene studies for SSRI response and systematic literature search of GWAS for SSRI response and MDD both revealed 245 and 800 genes to be significantly associated with SSRI response and MDD pathophysiology, respectively. Using these sets of genes, we retrieved

Data Extraction from Candidate Gene Studies of SSRI Response
All the 186 significantly associated SNPs (p-value ≤ 0.05) from the 90 selected articles are summarized in Table 1 (complete table in Table S1 online). Sample size along with their phenotype and drugs studied in each of the articles are also detailed in the same. In some research studies, authors had performed two types of analysis, one, responder versus non-responder and second, remitter versus non-remitter; here, we have tabulated all the significant SNPs reported in either kind of analysis.

Data Extraction from Candidate Gene Studies of SSRI Response
All the 186 significantly associated SNPs (p-value ≤ 0.05) from the 90 selected articles are summarized in Table 1 (complete table in Table S1 online). Sample size along with their phenotype and drugs studied in each of the articles are also detailed in the same. In some research studies, authors had performed two types of analysis, one, responder versus non-responder and second, remitter versus non-remitter; here, we have tabulated all the significant SNPs reported in either kind of analysis. N.A., data not available; M, male; F, female; FP, Follow-up period in weeks; OR, odds ratio; CI, confidence interval; Score, cumulative score for methodological quality assessment (see Tables S1 and S2 for detailed scoring) * Study name shows remission data. OR calculated using reported frequencies from the respective article. All the studies represent significant polymorphisms (p ≤ 0.05) with their corresponding genes.

Data Processing and Genetic Co-Occurrence of SSRI Response and MDD
A total of 186 significant SNPs from a systematic review of SSRI candidate gene studies and 423 SNPs from GWAS of SSRI response were merged to form the first data, Set A, representing SSRI response. Additionally, 3884 SNPs from GWAS of MDD formed the second data, Set B, relating to MDD pathophysiology. Beside this, both sets of SNPs were enriched with probable genetic variations found to be in LD (r 2 = 1) using the HaploReg v4.1 online tool. Now enriched Set A and Set B comprise 718 and 4769 SNPs, respectively. Enriched Set A and B thus obtained were mapped to genes using the SCAN and Ensembl VEP online tools. After excluding unmapped genes, pseudogenes, orfs, miRNAs, and uncharacterized genes, these annotated genes were assigned their respective HGNC IDs. This obtained a pool of 245 and 800 genes associated with SSRI response and MDD etiology, respectively (Table S5 online). Furthermore, we identified the overlap between SSRI-response genes (Set A) with that of MDD genes (Set B) and found 29 common genes. Since all the MDD associated genes were extracted from GWAS articles, i.e., an unbiased source, therefore, all of the 29 overlapped genes can be considered as unbiased molecular players co-occurring between SSRI response and MDD susceptibility ( Figure 2a). All the GWAS were performed on MDD patients vs. healthy controls.

Data Processing and Genetic Co-Occurrence of SSRI Response and MDD
A total of 186 significant SNPs from a systematic review of SSRI candidate gene studies and 423 SNPs from GWAS of SSRI response were merged to form the first data, Set A, representing SSRI response. Additionally, 3884 SNPs from GWAS of MDD formed the second data, Set B, relating to MDD pathophysiology. Beside this, both sets of SNPs were enriched with probable genetic variations found to be in LD (r² = 1) using the HaploReg v4.1 online tool. Now enriched Set A and Set B comprise 718 and 4769 SNPs, respectively. Enriched Set A and B thus obtained were mapped to genes using the SCAN and Ensembl VEP online tools. After excluding unmapped genes, pseudogenes, orfs, miRNAs, and uncharacterized genes, these annotated genes were assigned their respective HGNC IDs. This obtained a pool of 245 and 800 genes associated with SSRI response and MDD etiology, respectively (Table S5 online). Furthermore, we identified the overlap between SSRI-response genes (Set A) with that of MDD genes (Set B) and found 29 common genes. Since all the MDD associated genes were extracted from GWAS articles, i.e., an unbiased source, therefore, all of the 29 overlapped genes can be considered as unbiased molecular players co-occurring between SSRI response and MDD susceptibility (Figure 2a).

Functional Enrichment of SSRI Response and MDD Gene Pool
The WebGestalt web-tool was used to further enrich the genes obtained from systematic literature mining. The GO term and pathway enrichment analysis were performed keeping default GO slim classification, 0.05 significance level, and the minimum number of genes in a category was set at 3. The significance value was adjusted by the false discovery rate (FDR) analysis using the Benjamini-Hochberg (BH) procedure (Table S7 and S8 online). Set A of the SSRI-response genes were significantly enriched in 298 GO terms, involving 228 GO biological process terms, 15 GO molecular function terms, and 55 GO cellular component terms along with 24 KEGG pathway terms. Set B of the MDD etiology genes was significantly enriched in 94 GO terms, involving 7 GO biological process terms, 1 GO molecular function terms, and 86 GO cellular component terms along with 32 KEGG

Functional Enrichment of SSRI Response and MDD Gene Pool
The WebGestalt web-tool was used to further enrich the genes obtained from systematic literature mining. The GO term and pathway enrichment analysis were performed keeping default GO slim classification, 0.05 significance level, and the minimum number of genes in a category was set at 3. The significance value was adjusted by the false discovery rate (FDR) analysis using the Benjamini-Hochberg (BH) procedure (Tables S7 and S8 online). Set A of the SSRI-response genes were significantly enriched in 298 GO terms, involving 228 GO biological process terms, 15 GO molecular function terms, and 55 GO cellular component terms along with 24 KEGG pathway terms. Set B of the MDD etiology genes was significantly enriched in 94 GO terms, involving 7 GO biological process terms, 1 GO molecular function terms, and 86 GO cellular component terms along with 32 KEGG pathway terms. On overlapping, these results showed nine common KEGG pathways (Table S6 online) and 43 common GO terms as shown in Figure 2b,c. This overlap directs a common molecular ground between the mechanism behind antidepressant response and MDD pathophysiology. Since the biased data source (candidate study available literature) in Set A was overlapped with Set B which has an unbiased source (GWAS), common enriched biological pathways must be unbiased and can be used as putative future novel drug targets for MDD. Moreover, functional enrichment of overlapping gene set, i.e., Set C pinpoints the significant pathways involved in both antidepressant response and MDD etiology. pathway terms. On overlapping, these results showed nine common KEGG pathways (Table S6 online) and 43 common GO terms as shown in Figure 2b, 2c. This overlap directs a common molecular ground between the mechanism behind antidepressant response and MDD pathophysiology. Since the biased data source (candidate study available literature) in

Discussion
Major depressive disorder is a progressive brain disease with one of the leading cause of disability-adjusted life years affecting approximately 10-15% of the population worldwide [82]. Various treatment regimens are present for MDD, but SSRI pharmacotherapy is being commonly

Discussion
Major depressive disorder is a progressive brain disease with one of the leading cause of disability-adjusted life years affecting approximately 10-15% of the population worldwide [82]. Various treatment regimens are present for MDD, but SSRI pharmacotherapy is being commonly used and often recommended as a first-line treatment option in moderate-to-severe depression because of their higher efficacy and lower side effects compared to other antidepressants. Majorly, all the antidepressants available up to date are aimed at symptomatic management rather than complete cure owing to poorly understood MDD disease etiology. Since the response to antidepressant treatment varies markedly between individuals due to considerable clinical heterogeneity, therefore, the role of genetic predictors of antidepressant response and MDD disease per se are of utmost importance for the development of better treatment regimen which would further improve clinical management of MDD.
In this study, we used an integrative genetics approach to unbiasedly elucidate the possible association between antidepressant therapy and MDD etiology. Initially, an extensive literature search was done for identification of genetic variants involved in antidepressant response and disease etiology, followed by functional gene set enrichment, and thereby distinguishing the GO terms and molecular pathways involved in both antidepressant response and MDD susceptibility. For the first time in this uniquely designed study, to the best of our knowledge, has incorporated the amalgam of two exhaustive datasets of disease, concerning its pathophysiology and drug response genes, in search of precise drug targets to accelerate future drug development with minimal side effects.
For quality assessment of the studies included in this systematic review, based on the modified criteria used by Guin D et al. [83], the present study has incorporated the most comprehensive quality assessment scoring of research articles for screening SSRI response articles. Our systematic review has retrieved a total of 55 genes from the 90 research articles which were found to be involved in various synaptic transmission and neuronal development pathways. For instance, rs1083801 (GRM7, a glutamate receptor gene), was found to be most significantly associated with the early response with SSRIs [54]. The TPH2 gene and its variant, rs4760815, has been reported as probable risk factors for the development of MDD and also associated with SSRI response [84]. In addition, other genetic variants from TPH2, (rs11179027 and rs17110532); glutamate receptor ionotropic, kainate 2 gene, GRIK2 (rs543196), glutamic acid decarboxylase gene, GAD1 (rs3828275) and SLC6A4 (rs2066713) were also reported to be strongly associated with SSRI response [15]. SERTPR often regarded as 5-HTTLPR is a serotonin transporter gene promoter polymorphism which is been exhaustively studied in both remission rate and response rate of antidepressants [85]. A review and meta-analysis study by Kato and Serretti [86] has extensively studied and reported a significant association between the 5-HTTLPR variant and better response to antidepressants. Zill et al. [87] have identified functional polymorphism in the β1 adrenergic receptor 1165G>C (rs1801253), which was found to be involved in conferring the faster response towards antidepressant treatment, but did not influences the depressive phenotype. Furthermore, HTR2A, a serotonergic receptor gene variant rs7997012 is found to be associated with SSRI response as well as remission [56,59,88,89]. However, Illi A et al. [90] and Kishi T et al. [28] have reported a negative association of rs7997012 with SSRI response or remission. Thus, to overcome such inconsistencies across several research studies, we need an innovative approach to eliminate the biases and limitation of candidate gene studies [91]. Hence in the present study, we have opted for "candidate gene studies and GWAS overlapping" approach to overcome the lack-of-harmony among candidate gene studies. As a result, the systematic review of candidate gene studies and GWAS was merged to ascertain the unbiased genes and molecular pathways involved in both antidepressant response and MDD etiology.
After merging findings from the systematic review and GWAS studies, all the genetic variants were grouped into Set A and Set B followed by their LD enrichment, gene annotation, and functional enrichment analysis, to reveal multiple relevant pathways among antidepressant response and MDD etiology ( Figure 1). Furthermore, 56 KEGG pathways (24 from Set A and 32 from Set B) were identified with 9 commonly enriched pathways, namely Cushing syndrome, Retrograde endocannabinoid signaling, Axon guidance, cAMP signaling pathway, Glutamatergic synapse, Insulin secretion, Gap junction, Gastric acid secretion, and Salivary secretion, from Set A and Set B. In the view of MDD etiology, the most interesting molecular pathways are Axon guidance, Glutamatergic synapse, and Cushing syndrome. These 9 co-occurring pathways were further overlapped with the molecular pathways enriched from the 29 common genes implicated in antidepressant response and MDD etiology. This results in five overlapping molecular pathways, namely, Cushing syndrome, Axon guidance, cAMP signaling pathway, Insulin secretion, and Glutamatergic synapse. As MDD is a complex disorder, multiple pathways are expected to be involved in its etiology and our result is in concordance with the same showing the interplay of stress, neurodevelopmental, synaptic plasticity, and metabolic physiology in MDD etiology and antidepressant response.
Cushing syndrome, as being the top pathway of the analysis, is a metabolic disorder caused by overproduction of cortisol (glucocorticoid) produced by the adrenal cortex of the adrenal gland in response to glucocorticoids medication or in stressful condition. The adrenal gland is a part of stress-responsive hypothalamic pituitary adrenal (HPA) axis which consists of stimulating and feedback mechanism involving the hypothalamus, pituitary, and adrenal gland, and thus regulates the production of glucocorticoids. The HPA axis dysregulation has been implicated in the pathophysiology of many neuropsychiatry traits including MDD [92][93][94]. Genes associated with the homeostatic response to environmental stressors particularly lies in the HPA axis [95].
Corticotrophin-releasing hormone (CRH), released by the hypothalamus, stimulates the release of adrenocorticotropic hormone (ACTH) from pituitary which in turns regulate the levels of cortisol secreted from the adrenal gland. Hence CRH and its downstream effects are of prime importance which is controlled by cortisol level through feedback mechanism and two CRH receptors i.e., CRHR1 and CRHR2 [96]. Alteration in CRHR1 and CRHR2 activity may lead to HPA dysregulation and can be a major risk factor for depressive symptoms. Bradley et al. [97] demonstrated the role of genetic variants in CRHR1 as moderators of the effects of child abuse on adult depressive symptoms in two independent populations which was further confirmed by a replication study which reported the association of TAT (rs7209436, rs110402, and rs242924) haplotype in CRHR1 in predicting the adult depression [98,99]. Woody et al. [100] also supported the hypothesis of protective CRHR1 haplotype (TAT) as their result suggested the development of brooding among children without the protective CRHR1 haplotype. Moreover, animal studies also indicate the prenatal glucocorticoid exposure can cause epigenetic instability in CRHR1 promoter which further increases the risk for the affective disorder in offspring's across two generations [101]. Epistatic interaction between AVPR1b, CRHR1, and BDNF genes has also been reported to be involved in susceptibility to MDD [102,103]. Ressler et al. [104] reported the involvement of interaction of 5-HTTLPR S allele with CRHR1 haplotype in predicting adult depression in individuals with child abuse. The SNPs rs110402, rs242924 rs3779250, rs7209436, and rs173365 from CRHR1 and CRHR2 genes were reported to be positively associated with MDD in the Japanese population [105]. However, a CRHR2-based 10 SNP study showed no significant allelic or genotypic differences among unipolar patients and matched healthy controls [106].
In addition, there are evidences for an antidepressant response via glucocorticoid receptors (GRs) where it has been demonstrated that sertraline increases human hippocampal neurogenesis via a GR-dependent mechanism [107]. A recent study has shown the significant association of rs41423247 polymorphism with fluoxetine response in depressed patients [10]. Another study of haplotype-tag SNPs (rs1876828, rs242939 and rs242941) in CRHR1 shown the significantly improved response to antidepressants among highly anxious patients homozygous for the GAG haplotype, suggesting the possible role of CRHR1 and other stress-inflammatory pathway genes in variable antidepressant response [108,109].
Increasing evidence suggests the dysfunction of glutamatergic neurotransmission impairs neuroplasticity in the brain and this leads to major depressive disorder [110,111]. Many novel targets in relation to this pathway are evidenced for causing the depression phenotype. In the past decade, NMDA has gained the utmost attention with respect to the biology of depression and also serves as a potential target for drug development and treatment. As under pathological condition, elevated levels of glutamate resulted in the impairment of synaptic plasticity and even excitotoxicity. Alternate depression hypothesis includes exposure to stress inculpates the release of the high level of stress hormone i.e., cortisol from the adrenal gland. Correspondingly, adequate literature demonstrates that stress and glucocorticoids are responsible for alteration in the expression and activity of vesicular proteins of the neurotransmission of glutamate [112]. Emerging evidences from post-mortem studies has reported the dysregulation of genes and increased glutamate levels, thus highlighting the role of altered glutamate signaling in MDD patients [113]. Decrease expression of presynaptic genes in MDD patients such as SYN3, SNAP25, essential for vesicular release of neurotransmitters have been reported. Likewise, a significant downregulation of postsynaptic genes was also reported in the dentate gyrus (DG) and CA1 of MDD patients such as AMPA receptors, specifically GLU1 and GLU3 [114]. In the same note, several NMDA gene variations have been reported to have an association with an abnormality in NMDA receptors. One such study had shown significant association of GRIN1 (rs4880213) with depression. Further studies revealed that variations in the GRIN2B were associated with schizophrenia, psychiatric disorders, and brain plasticity [115]. Similarly, downregulation of NMDA receptor subunits GRIN1A and GRIN2B, as well as PSD-95 have been demonstrated in the anterior prefrontal cortex of MDD subjects [116]. Likewise, metabotropic glutamate receptor 7 (GRM7) encodes the protein mGluR7 mediates the glutamate neurotransmission, found to be involved in the development of the major depressive disorder. Li W et al. [117] studied the association of genetic variation of the GRM7 gene with MDD and schizophrenia and reported the significant association of GRM7 gene variation (rs779706) with MDD and (rs2229902 and rs9870680) with schizophrenia in the Han Chinese population. Thus, it is hypothesized that dysfunction of ionotropic and metabotropic receptors are associated with the depressed phenotype. However, in a resequencing study, none of the GRM7 variants had shown the significance level with MDD in the Dutch cohort [118].
Also, a number of studies support the association between the antidepressant response and GRIK4 in MDD patients [89,119,120]. Genes encoded ionotropic glutamate receptors were studied with respect to citalopram treatment and demonstrated a significant association of two SNPs (rs4825476 and rs2518224) located within GRIA3 and GRIK2, respectively, with the treatment-emergent suicidal ideation in MDD patients [121]. Moreover, in order to explore the potential targets for treatment-resistant depressive patients, emerging evidence from clinical trials supported the use of glutamate receptor modulators for the treatment of depression and these include non-competitive NMDA receptor antagonists such as ketamine, subunit (NR2B)-specific NMDA receptor antagonists, NMDA receptor glycine-site partial agonists and metabotropic glutamate receptor (mGluR) modulators [122]. Therefore, glutamate pathways and its associated receptors are important and further insights and detailed understanding could help us to target the accurate site for future drug development.
Neuronal circuit formation involves a molecular cascade of events such as axon guidance, where axons move to their target cells in a complex, constantly changing the environment. It has been speculated that change in this gene-environment interaction may lead to alteration in axon guidance followed by its implication in neuropsychiatric disorders. Furthermore, a meta-analysis by S Jovanova O et al. [123], has reported three methylated sites associated with depressive symptoms and were also found to be involved in axon guidance as a pathway in major depression. Interestingly, research articles have also shown that miRNAs implicated in axon guidance are involved in differential antidepressant response where miRNAs namely miR-146a-5p, miR-146b-5p, miR-221-3p, miR-24-3p, miR-26a-5p are known to be involved in axon guidance and were also associated with antidepressant response in MDD patients [124,125]. Our current finding of Slit guidance ligand 3 gene, SLIT3 in antidepressant and MDD gene set is in consensus with findings from Glessner JT et al. [126], where they have performed genome-wide copy number variation scan of large cohort of MDD patients and controls and has observed 5q35.1 as the most significant locus harboring the SLIT3 gene which is integral to repulsive axon guidance. Thus synaptic plasticity mediated via axon guidance is a topic of new research and can further be studied for better antidepressant development.
Since the majority of intracellular messenger cascades are regulated either by G protein-coupled receptors (GPCRs) or protein tyrosine kinases (PKAs) and are major activator of various cellular and molecular signaling pathways such as cortisol secretion in the adrenal cortex contributing to Cushing's syndrome and neurotransmitters such as serotonin and norepinephrine, which are known to mediate the effects of antidepressant treatments modulates secondary messenger cascades via interaction with GPCRs or PKAs [127,128]. Hence, receptor activation induced by ligands such as neurotransmitters confers to cAMP generation via stimulating adenylate cyclase (AC), and their binding to the G-protein subtype leads to activation of PKA, which is an important factor for driving several biological functions either by phosphorylation or dephosphorylation of specific target proteins, undermines the antidepressant actions. Furthermore, the cAMP response element binding protein (CREB), a transcription factor that mediates the actions of the cAMP cascade, is a substrate for PKA, is involved in regulating gene expression, and has the capability to modulate their transcriptional activity, which is important for cellular adaptions during antidepressant administration.
Dowlatshahi et al. [129] have reported low CREB levels in post-mortem temporal cortex of naive major depressive disorder patients as compared to MDD patients treated with antidepressants. Odagaki et al. [130] have even demonstrated the increased immunoreactivities of phosphorylated CREB as well as total CREB levels in the prefrontal cortex of depressed suicide victims and specifically in antidepressant drug-free subjects. Also, with a slight trend for increased levels of PKA-Cα in depressed suicide victims only as compared to healthy controls. As a downstream consequence, the expression of various target genes critical to the organization of neuronal networks and synaptic plasticity, like neurotrophin, brain-derived neurotrophic factor (BDNF), and neuropeptide Y (NPY) is also increased contributing towards antidepressant-mediated changes in structural remodeling, neuronal plasticity and synaptic restructuring [131,132]. Henceforth, evidence from such studies confer the ability of HTR receptors in either stimulating or inhibiting the AC-cAMP-PKA signaling transduction pathway but further invokes solicitation and functional validation of pathway genes in undermining the action of antidepressants, as this transduction pathway is highly regulated by several other factors such as stress, apoptosis, inflammation, and others.
Hence pathway enrichment analysis of "biased and unbiased" merged data is assured to facilitate our understanding of the underlying molecular mechanism of the complex trait of anti-depressant response and major depression as a disease, to circumvent the symptomatic respite and design a definite therapy. Application of this elaborated analytical tactics in translational research concerning complex disorders like MDD would be beneficial after in-vitro and in-vivo validation of the top promising pathways and genes involved in antidepressant response and MDD etiology.

Methodology
The complete workflow of the study is represented in Figure 4. We initiated with integrating significantly associated SNPs from candidate genes studies and GWAS concerning SSRI response and from GWAS of MDD etiology. Therefore all the SNPs from candidate gene studies and GWAS of SSRI response were merged in one group and another group consist significant genetic variations from MDD GWAS. Further, they were annotated into genes followed by functional and pathway enrichment analysis.

A systematic Review of Antidepressant Response Candidate Gene Studies
Among currently available antidepressants, selective serotonin reuptake inhibitors (SSRIs) like escitalopram, sertraline, and fluoxetine are the most commonly prescribed drugs [133] considering their higher efficacy and lower side-effects. Hence in the present study, we have considered "SSRIs" and "antidepressants" interchangeably in this manuscript. A systematic literature search was performed in accordance with PRISMA guidelines [134]. The MEDLINE and Web of Science databases were searched using Medical Subject heading (MeSH) terms "selective serotonin reuptake inhibitor", "SSRI", "pharmacogenetics", "pharmacogenomics", "response", "treatment outcome", "SNP", "variant", "polymorphism" with AND/OR Boolean operators to extract all the studies evaluating association of SNPs with SSRI treatment outcome in patients with MDD. The search and study selection was carried out independently by three authors (AS, HG, and PS) covering the articles published till 28th February 2018.

A systematic Review of Antidepressant Response Candidate Gene Studies
Among currently available antidepressants, selective serotonin reuptake inhibitors (SSRIs) like escitalopram, sertraline, and fluoxetine are the most commonly prescribed drugs [133] considering their higher efficacy and lower side-effects. Hence in the present study, we have considered "SSRIs" and "antidepressants" interchangeably in this manuscript. A systematic literature search was performed in accordance with PRISMA guidelines [134]. The MEDLINE and Web of Science databases were searched using Medical Subject heading (MeSH) terms "selective serotonin reuptake inhibitor", "SSRI", "pharmacogenetics", "pharmacogenomics", "response", "treatment outcome", "SNP", "variant", "polymorphism" with AND/OR Boolean operators to extract all the studies evaluating association of SNPs with SSRI treatment outcome in patients with MDD. The search and study selection was carried out independently by three authors (AS, HG, and PS) covering the articles published till 28th February 2018.
The searches were confined to human and English language studies. All the articles that were review, meta-analysis, commentary, editorial, clinical trial, letter, randomized control trial, technical report were excluded. Articles were sorted for their relevance at two stages, first using the title and second using their abstract. At the first stage of title screening, duplicates, reviews, systematic The searches were confined to human and English language studies. All the articles that were review, meta-analysis, commentary, editorial, clinical trial, letter, randomized control trial, technical report were excluded. Articles were sorted for their relevance at two stages, first using the title and second using their abstract. At the first stage of title screening, duplicates, reviews, systematic reviews, meta-analysis, non-human, and co-morbid studies were removed. Secondly, the abstracts of all remaining articles were retrieved and screened based on the inclusion and exclusion criteria of the study. Inclusion criteria required genetic association of SSRI response in MDD patients, with age range of 18-75 years, and that response should be inferred based on gold standard severity rating scales like Hamilton Depression Rating Scale (HAM-D) or Montgomery Asberg Depression Rating Scale (MADRS). Whereas studies involving MDD as co-morbidity or MDD patients with other severe medical illness, psychiatric disorder or substance abuse were excluded. Detailed exclusion criteria are given in Figure 1.

Data Extraction and Quality Assessment of Antidepressant Response Candidate Gene Articles
Data, from selected full-text articles, were extracted by AS, HG, and PS and checked by DG and RK. Ethnicity, used response criteria, sample size, genes with information of genetic variant, related genotypic or allelic frequency and respective reported p-value, drug, dose, and follow-up period were included in the data collection table. Ethnicity was classified as reported in the respective article, else the country in which the study was conducted, assumed to be the individual's ethnicity. A cut-off p-value ≤ 0.05 was used for extracting SNPs from the selected articles. Methodological quality of each article was assessed by two independent reviewers (AS and HG) using modified criteria for quality assessment, as used by Guin D et al. [83]. The quality assessment was scored on 8 parameters (Table  S2 online), with a positive score awarded for each detail present in the study, the lack of detail was described as 0. Conflicting scores were reached to a consensus upon discussing with DG and RK. If the score obtained was 7 or higher, the study was considered as high quality.

A Systematic Literature Search of GWAS of SSRI Response and MDD Susceptibility
Data were also extracted from GWAS, using keywords "Depression" or "Depressive Disorder", correlating genetic variability with SSRI response in patients with MDD, and distinguishing MDD patients from healthy controls. To maintain homogeneity in studies, only those articles were selected where researchers have opted peripheral blood as their source of DNA extraction. In the case of disease vs. control GWAS, the diagnosis has been made by a psychiatrist/clinical psychologist, otherwise, a study has been excluded from the pool. Similarly, in GWAS correlating SNPs with SSRI response, if the response were adjudicated using standard depression severity rating scales the study was included, else excluded. Studies considering meta-analysis were also excluded. A cut-off p-value ≤ 0.0001, was used for extracting SNPs reported in respective GWAS.

Data Processing of Candidate Gene Studies for SSRI Response and GWAS of SSRI Response and MDD Susceptibility
SNPs from candidate and GWAS studies were retrieved and categorized into two data sets, Set A containing SNPs associated with SSRI response (data from both systematic review and GWAS of SSRI response), and Set B comprising SNPs associated with MDD susceptibility (data only from GWAS of MDD disease). Here, we were investigating the targets for etiology based antidepressant development using genomic integration of disease and drug response, one dataset needed to be unbiased and hence we decided not to include candidate gene studies for MDD. Thus, Set A can be considered to be a biased data set as it contains SNPs from candidate gene studies as well as from GWAS, whereas Set B is genetically unbiased data set as it contains data from GWAS only. As interpatient genetic variability can modulate SSRI clinical response and MDD etiology, therefore, in order to widen the genetic region to find out the biological relevance of the associated SNPs, we have incorporated all the probable genetic variations found to be in linkage disequilibrium (LD) with previously extracted SNPs of Set A and Set B. Furthermore, pairwise LD among the extracted SNPs of set A and B were performed using HaploReg v4.1 which uses mammalian conservation algorithm from GERP and SiPhy-omega with LD threshold (r 2 ) = 1, corresponds to 100% LD [135]. SNPs that are found to be in 100% LD in all the four major population i.e., African, American, Asian, and European were included in the respective SNP sets. In addition, SCAN (SNP and CNV Annotation Database) [136] (http://www.scandb.org/newinterface/about.html) and VEP (Variant Effect Predictor -Ensembl) [137] (https://asia.ensembl.org/info/docs/tools/vep/index.html) online tools were employed for gene annotation for all the identified SNPs, which were extracted and enriched from LD. The SCAN utilizes two ways for SNP annotation, i.e., relative position based and eQTLs (expression quantitative trait loci) method. Relative position-based method identifies a gene based on SNP position (intronic, inter-genic, etc.) or an intergenic variant can be annotated to a gene if it is in LD with any other variant present in the gene. An eQTLs-based method annotates a gene whose quantitative expression is altered by input SNP. Whereas VEP annotates an SNP to a gene if it has a functional effect on that gene, transcript, protein sequence or regulatory regions. The SNPs, which remained unmapped, were excluded from further analysis. Genes thus obtained were assigned HGNC (HUGO Gene Nomenclature Committee) (https://www.genenames.org/) IDs manually [138]. The pseudogenes, hypothetical loci, non-coding RNAs, non-protein coding genes, open reading frames (orfs), microRNA (miRNA), and uncharacterized genes were excluded from each dataset. Further, a third data set, Set C, was also framed which contained all the overlapping genes from Set A and Set B. Set C was used to distinguish the most significant pathways from the commonly enriched pathways among Set A and Set B.

Functional Enrichment Analysis
For each of the three datasets, functional enrichment analysis was conducted using a WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) [9]. Gene Ontology (GO) term enrichment analysis (http: //www.geneontology.org/) and pathway enrichment analysis with Kyoto Encyclopedia of Genes and Genomes (KEGG) engine (https://www.genome.jp/kegg/) was performed keeping default GO slim classification, 0.05 significance level and a minimum number of genes in a category was set at 3. For enrichment analysis of each dataset, the significance value was adjusted by the false discovery rate (FDR) analysis using the Benjamini-Hochberg (BH) procedure. The GO terms and molecular pathways overlapping between Set A and Set B were extracted, which further overlaid with the enriched pathways from Set C to identify pathophysiological grounds which can be targeted to develop etiology based antidepressants.

Conclusions
Burgeoning literature evidence has so far pointed out the clinical relevance of antidepressant response and MDD etiology individually, but such studies are not consistent enough to manifest these findings into clinical practice. In the present study, we have focused on highlighting the possible confounding factors responsible for antidepressant response and MDD pathophysiology altogether. Thus we have shortlisted significant genes and pathways implicated in the same which can be further utilized as novel molecular targets for the development of more efficacious antidepressant drugs. Application of this elaborated analytical approach in translational research concerning complex disorders like MDD would be beneficial after in-vitro and in-vivo validation.
Author Contributions: A.S. conceived and wrote the manuscript. A.S. and R.K. contributed to the study design. A.S., H.G., and P.S. extracted data from the literature. A.S. did the gene set enrichment analysis shown in the manuscript. In case of disagreement in systematic review or cumulative quality assessment, D.G. and H.K. helped in resolving. A.S., D.G., N.K., and P.S. contributed in writing manuscript. R.K., D.V., H.K., L.S., and J.Y. contributed for critical evaluation of the study and improving the manuscript. M.S. and R.K.C. provided valuable insights for selection of relevant articles from clinical perspectives. R.K. conceived, interpreted, and supervised the study design. All the authors read and approved the final manuscript.