# The association of ARRB1 polymorphisms with response to antidepressant treatment in depressed patients

Kenneth Chappell, Abd El Kader Ait Tayeb, Romain Colle, Jérôme Bouligand, Khalil El-Asmar, Florence Gressier, Séverine Trabado, Denis Joseph David, Bruno Feve, Laurent Becquemont, et al.

► To cite this version:

Kenneth Chappell, Abd El Kader Ait Tayeb, Romain Colle, Jérôme Bouligand, Khalil El-Asmar, et al. The association of ARRB1 polymorphisms with response to antidepressant treatment in depressed patients. Frontiers in Pharmacology, 2022, 13, pp.974570. 10.3389/fphar.2022.974570. hal-04021229

HAL Id: hal-04021229

https://hal.science/hal-04021229

Submitted on 6 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY

Renan Pedra de Souza,
Universidade Federal de Minas Gerais,
Brazil

REVIEWED BY

Gualberto Ruaño,
Hartford Hospital, United States
Hussein Kadhem Al-Hakeim,
University of Kufa, Iraq

*CORRESPONDENCE

Kenneth Chappell,
kennyrchappell@gmail.com

¹These authors have contributed equally to this work

SPECIALTY SECTION

This article was submitted to Pharmacogenetics and Pharmacogenomics, a section of the journal Frontiers in Pharmacology

RECEIVED 21 June 2022

ACCEPTED 11 October 2022

PUBLISHED 26 October 2022

CITATION

Chappell K, Ait Tayeb AEK, Colle R, Bouligand J, El-Asmar K, Gressier F, Trabado S, David DJ, Feve B, Becquemont L, Corruble E and Verstuyft C (2022), The association of ARRB1 polymorphisms with response to antidepressant treatment in depressed patients.
*Front. Pharmacol.* **13**:974570.
doi: **10.3389/fphar.2022.974570**

COPYRIGHT

© 2022 Chappell, Ait Tayeb, Colle, Bouligand, El-Asmar, Gressier, Trabado, David, Feve, Becquemont, Corruble and Verstuyft. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The association of *ARRB1* polymorphisms with response to antidepressant treatment in depressed patients

Kenneth Chappell¹*, Abd El Kader Ait Tayeb¹,², Romain Colle¹,², Jérôme Bouligand³,⁴, Khalil El-Asmar¹,⁵, Florence Gressier¹,², Séverine Trabado³,⁴, Denis Joseph David¹,⁶, Bruno Feve⁷, Laurent Becquemont¹,⁸, Emmanuelle Corruble²,⁶† and Céline Verstuyft¹,⁴,⁹†

¹Université Paris-Saclay, UMR 1018, CESP-Inserm, Team MOODS, Faculté de Pharmacie, Bâtiment Henri MOISSAN, Orsay, France, ²Service Hospitalo-Universitaire de Psychiatrie de Bicêtre, Hôpitaux Universitaires Paris-Saclay, Assistance Publique-Hôpitaux de Paris, Hôpital de Bicêtre, Le Kremlin Bicêtre, France, ³INSERM UMR-S U1185, Faculté de Médecine, University Paris-Saclay, Le Kremlin Bicêtre, France, ⁴Service de Génétique Moléculaire, Pharmacogénétique et Hormonologie, Hôpitaux Universitaires Paris-Saclay, Assistance Publique-Hôpitaux de Paris, Hôpital de Bicêtre, Le Kremlin Bicêtre, France, ⁵Department of Epidemiology and Population Health, Faculty of Health Sciences, American University of Beirut, Beirut, Lebanon, ⁶CESP, MOODS Team, INSERM UMR 1018, Faculté de Médecine, University Paris-Saclay, Le Kremlin Bicêtre, France, ⁷Sorbonne Université-INSERM, Centre de Recherche Saint-Antoine, UMR S938, Institut Hospitalo-Universitaire ICAN, Service d'Endocrinologie, Hôpital Saint-Antoine, Assistance Publique-Hôpitaux de Paris, Centre de Référence des Maladies Rares de l'Insulino-Sécrétion et de l'Insulino-Sensibilité, Paris, France, ⁸Centre de Recherche Clinique Paris-Saclay, Hôpitaux Universitaires Paris-Saclay, Assistance Publique-Hôpitaux de Paris, Hôpital de Bicêtre, Le Kremlin Bicêtre, France, ⁹Centre de Ressources Biologiques Paris-Saclay, Hôpitaux Universitaires Paris-Saclay, Assistance Publique-Hôpitaux de Paris, Hôpital de Bicêtre, Le Kremlin Bicêtre, France

**Introduction:** β-arrestin 1, a protein encoded by *ARRB1* involved in receptor signaling, is a potential biomarker for the response to antidepressant drug (ATD) treatment in depression. We examined *ARRB1* genetic variants for their association with response following ATD treatment in METADAP, a cohort of 6-month ATD-treated depressed patients.

**Methods:** Patients (n = 388) were assessed at baseline (M0) and after 1 (M1), 3 (M3), and 6 months (M6) of treatment for Hamilton Depression Rating Scale (HDRS) changes, response, and remission. Whole-gene *ARRB1* variants identified from high-throughput sequencing were separated by a minor allele frequency (MAF)≥5%. Frequent variants (i.e., MAF≥5%) annotated by RegulomeDB as likely affecting transcription factor binding were analyzed using mixed-effects models. Rare variants (i.e., MAF<5%) were analyzed using a variant set analysis.

**Results:** The variant set analysis of rare variants was significant in explaining HDRS score changes ($T = 878.9; p = 0.0033$) and remission ($T = -1974.1; p = 0.034$). Rare variant counts were significant in explaining response ($p = 0.016$), remission ($p = 0.022$), and HDRS scores at M1 ($p = 0.0021$) and M3 ($p<0.001$). rs553664 and rs536852 were significantly associated with the HDRS score (rs553664: $p = 0.0055$ | rs536852: $p = 0.046$) and remission (rs553664: $p = 0.026$ | rs536852: $p = 0.012$) through their interactions with time. At M6,
significantly higher HDRS scores were observed in rs553664 AA homozygotes (13.98 ± 1.06) compared to AG heterozygotes (10.59 ± 0.86; *p* = 0.014) and in rs536852 GG homozygotes (14.88 ± 1.10) compared to AG heterozygotes (11.26 ± 0.95; *p* = 0.0061). Significantly lower remitter rates were observed in rs536852 GG homozygotes (8%, *n* = 56) compared to AG heterozygotes (42%, *n* = 105) at M6 (*p* = 0.0018).

**Conclusion:** Our results suggest *ARRB1* variants may influence the response to ATD treatment in depressed patients. Further analysis of functional *ARRB1* variants and rare variant burden in other populations would help corroborate our exploratory analysis. β-arrestin 1 and genetic variants of *ARRB1* may be useful clinical biomarkers for clinical improvement following ATD treatment in depressed individuals.

**Clinical Trial Registration:** clinicaltrials.gov; identifier NCT00526383

KEYWORDS

high-throughput sequencing, pharmacogenetics, major depressive disorder, β
Arrestin, ARRB1

# 1 Introduction

Major Depressive Disorder (MDD) is the leading contributor to global disability (World Health Organization, 2017). Its main treatment option, antidepressant drugs (ATD), is only modestly effective, with more than 65% of ATD-treated patients failing to achieve remission (Trivedi et al., 2008). Given the rising global burden of MDD and the lack of effective treatments, it is necessary to identify robust biomarkers to predict treatment response.

Two potential candidates, β-arrestin 1 and 2, are ubiquitously expressed proteins encoded by the *ARRB1* and *ARRB2* genes, respectively (Parruti et al., 1993; Lefkowitz and Shenoy, 2005). Both have functions in receptor desensitization, especially with many G protein-coupled receptor (GPCR) families, though roles in protein scaffolding, clathrin-mediated endocytosis, and G-protein-independent signaling are known (Shukla et al., 2011; Bond et al., 2019). Indeed, the β-arrestins serve as scaffolds for various proteins, including the phosphoinositide-3, extracellular signal-regulated 1 and 2 (ERK1/2), and several mitogen-activated protein (MAP) kinases, whose signaling pathways are implicated in the treatment of mood disorders (Golan et al., 2009; Shukla et al., 2011; Bond et al., 2019). Additionally, β-arrestin 1 has been observed to function as a cytoplasm-to-nucleus messenger and nuclear protein scaffold (Golan et al., 2009). Interactions with GPCRs include serotonin, dopamine, adrenergic, and melatonin receptors, which are associated with different brain pathways and implicated in the pathophysiology of MDD (Kovacs et al., 2009; Petit et al., 2018).

Furthermore, biased agonists can preferentially activate downstream pathways. With respect to GPCRs, this includes classical G-protein- and/or β-arrestin-mediated pathways. Importantly, biased signaling may lend itself to the develop of

novel therapeutic strategies in diverse pharmacological domains (Violin and Lefkowitz, 2007; Shukla et al., 2011; Bond et al., 2019). An example of biased agonism with respect to antidepressants was demonstrated with the tricyclic antidepressant, desipramine, a norepinephrine reuptake inhibitor and α2-adrenergic receptor ligand that was observed to lead to β-arrestin-mediated α2-adrenergic receptor internalization and downregulation in mouse embryonic fibroblasts (Cottingham et al., 2011).

Lower *ARRB1* mRNA and β-arrestin 1 protein levels have been observed in murine models of depression, as well as in human mononuclear leukocytes of depressed patients, compared to healthy controls (Avissar et al., 2004; Matuzany-Ruban et al., 2005; Golan et al., 2013; Lipina et al., 2013; Mendez-David et al., 2013). These levels returned to levels similar to those of healthy controls in both murine models and depressed patients following 4 weeks of antidepressant treatment, a change that preceded clinical improvement (Golan et al., 2013). A timeline of studies examining β-arrestin 1 and *ARRB1* in the context of depression and its treatment with ATDs is summarized in Figure 1. These associations suggest that β-arrestin 1 may be a promising candidate biomarker for ATD treatment response in MDD patients (Mendez-David et al., 2013). Factors that may alter its expression or function, such as genetic variants, are thus of interest.

*ARRB1*, located on chromosome 11q13 and comprising 16 exons, encodes for two splicing isoforms of β-arrestin 1 (Parruti et al., 1993; Calabrese et al., 1994; O'Leary et al., 2016) (see Figure 2). Single nucleotide polymorphisms (SNP) of *ARRB1* have been studied principally in the neuropsychiatric contexts of nicotine dependence and depression (Sun et al., 2008; Chang et al., 2015). In a European cohort, a haplotype constructed from 6 SNPs within *ARRB1* was significantly associated with nicotine dependence (Sun et al., 2008). In a meta-analysis of 3 cohorts of depressed individuals of European
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

ancestry, the intronic rs17133921 SNP was associated with MDD, specifically in the subgroup of females with recurrent MDD and an onset age of MDD prior to the age of 31 (Shyn et al., 2011). Furthermore, the rs12274033 SNP upstream of *ARRB1*, and a haplotype including it, were both associated with the response to mirtazapine treatment in a Korean population of MDD (Chang et al., 2015). These findings suggest that genetic polymorphisms of *ARRB1* may have an influence in different psychiatric contexts and could thus serve as potential candidate biomarkers for the response to ATD treatment in the context of MDD.

Advances in genetic analysis thanks to genome-wide association studies (GWAS) and high-throughput sequencing (HTS) technology have greatly increased the identification of genetic variants associated with disease. Furthermore, the Encyclopedia of DNA elements (ENCODE) project has mapped the vast majority of the human genome to functional elements involved in protein coding, transcription factor (TF) binding, and chromatin structure (The ENCODE Project Consortium, 2012). Importantly, these data can be used to annotate genetic
variants with regulatory information and offer insight into
their potential functional consequences in various disease
contexts. RegulomeDB is one such database that uses these
and other predictive data to rank genetic variants on their
likelihood to lie in a functional element and thus have a
functional consequence (Boyle et al., 2012).

However, most variants detected by HTS have a minor
allele frequency (MAF) less than the 5% threshold principally
used to delineate common variants (Momozawa et al., 2016).
The analysis of these rarer variants is thus generally more
difficult and underpowered. One approach to this issue is
grouping variants into sets—defined by genes or genomic
regions—and analyzing the association between this
cumulative variable and a trait of interest (Chen et al.,
2019). Recently, the optimal unified sequence kernel
association test (SKAT-O) variant set method was used to
identify significant associations between rarer variants in
ADH1C and GSTO1 and mRNA expression from a set of
303 pharmacokinetic-related genes (Klein et al., 2019).
Importantly, SKAT-O considers both burden and SKAT
tests, which groups variants with respect to the mean or
the variance, respectively. It thus optimizes test selection
depending on whether most variants act in the same
direction (burden) or in opposing directions, i.e., both
protectively and deleteriously (SKAT) (Lee et al., 2012).

Given that both β-arrestin 1 levels and *ARRB1* genetic
variants have been associated with MDD and ATD treatment
response, the objectives of this study were to analyze the
association of clinical response following ATD treatment in
a cohort of 6-month ATD-treated MDD patients, with 1) the
set of rare variants and 2) frequent variants with likely
functional consequences, identified from exonic, intronic,
and 5'- and 3'-untranslated regions of *ARRB1*.

2 Materials and methods

## 2.1 Study design

Do Antidepressants Induce Metabolic Syndromes (METADAP) is a 6-month prospective, multicentric, and observational cohort study carried out in a psychiatric setting (Corruble et al., 2015). Patients with a current major depressive episode (MDE) in the context of MDD were treated in naturalistic conditions and assessed before and after beginning a new ATD treatment. This study was registered by the French National Agency for Medicine and Health Products Safety (ANSM) and the Commission Nationale de l'Informatique et des Libertés (CNIL). It was approved by the Ethics Committee of Paris-Boulogne (France) and conformed to international ethical standards (ClinicalTrials.gov identifier: NCT00526383).

2.2 Patient population

Males and females 18-65 years of age without a serious
medical condition were recruited. Inclusion criteria were
presentation with a current MDE in the context of MDD
(DSM-IVTR), as assessed by the Mini International
Neuropsychiatric Interview (MINI) and a score ≥18 on the
17-item Hamilton Depression Rating Scale (HDRS)
(Hamilton, 1960), as well as need for a new ATD
treatment. Patients with psychotic symptoms or other
mental disorders—such as psychotic disorder, bipolar
disorder, alcohol or drug dependence, or an eating
disorder—or those who were pregnant or who had organic
brain syndromes or serious medical conditions, were
excluded. Patients using antipsychotics or mood stabilizers
before inclusion and/or for at least 4 months during the year
before inclusion were also excluded. Concurrent use with
antipsychotics, mood stabilizers, and/or stimulants was not
allowed during the study. Benzodiazepines, at the minimum
effective dose and for the minimum time period, and
psychotherapies were allowed (Corruble et al., 2015).
Measures and samples were obtained prior to beginning
ATD treatment (M0), and after 1 month (M1), 3 months
(M3), and 6 months (M6) of ATD treatment. Data were
collected from November 2009 to March 2013 in
6 university hospital (CHU) psychiatry departments across
France in Paris (CHU Fernand-Widal and CHU Saint-
Antoine), Le Kremlin-Bicêtre (CHU Bicêtre), Grenoble
(CHU Grenoble Alps), Lille (CHU Michel Fontan 1), and
Besançon (CHU of Besançon) (Corruble et al., 2015).

Of the 643 patients included in METADAP, 19 had major
protocol deviations and were excluded. Of the 624 available
for analysis, 519 patients provided samples for genetic studies.
From these, 400 have undergone HTS. Two samples were
removed due to technical difficulties. As our objective was to
examine clinical response following ATD treatment,
10 patients without follow-up measures (i.e., at M1, M3,
and M6) were removed. Thus, 388 patients were analyzed.
Genetic analyses may have included fewer than 388 patients if
overall call quality was high, but individual calls were poor
and removed (discussed further below). The reasons for
dropout included loss to follow-up (n = 87, 48.6%), ATD
changes (n = 75, 41.9%), use of unauthorized drugs, including
antipsychotics, mood stabilizers, and stimulants (n = 8, 4.5%),
the presence of an exclusion criterion (i.e., unstable medical
condition, psychiatric disorder, substance abuse, or
pregnancy) during follow-up (n = 8, 4.5%), or death (n = 1,
0.5%). Missing data at M1 (n = 21) included 8 dropouts
(38.1%; 4 lost to follow-up and 4 due to ATD changes),
while missing data at M3 (n = 121) included 110 dropouts
(90.9%; 62 lost to follow-up, 36 due to ATD changes, 6 due to
unauthorized drug use, 5 due to the presence of an exclusion
criterion, and 1 due to death. Information about ancestry was
self-reported. Caucasian ethnicity was defined as having
Caucasian parents, African as having Sub-Saharan African
and/or Afro-Caribbean parents, and Asian as having East
Asian, Central Asian, and/or South Asian parents (Morales
et al., 2018). Three Asian patients were present within the
analyzed population and were combined with the seven
patients of mixed ethnicity. Education was classified into
3 levels: primary [e.g., elementary (≤5 years)]; secondary
[e.g., middle school and high school (>5 and ≤12 years)];
tertiary [e.g., any college- or university-level education
(>12 years)]. All patients provided written informed
consent for study participation and for genetic analyses.

## 2.3 Antidepressant treatment

ATD monotherapies were prescribed by a psychiatrist in a “real
world” psychiatric treatment setting as previously described (Corruble
et al., 2015). ATDs belonged to 1 of 4 classes: selective serotonin
reuptake inhibitors (SSRI), serotonin norepinephrine reuptake
inhibitors (SNRI), tricyclic antidepressants (TCA), or other ATD
treatments. In this ancillary study of 388 patients, 40% were
prescribed an SSRI, 41% an SNRI, 6% a TCA, and 9% another
ATD treatment; 4% received electroconvulsive therapy (ECT) instead
of an ATD monotherapy. If a change in treatment was required
during follow-up, the patient was dropped from the study.

## 2.4 Assessment of antidepressant treatment response

The HDRS was used to assess depression severity at M0
(i.e., baseline) and response to treatment at M1, M3, and M6.
Responders to treatment were classified by an improved HDRS
score of ≥50% relative to baseline, remitters by a HDRS
score ≤7 after at least 4 weeks of treatment, as recommended
by the American College of Neuropsychopharmacology (ACNP)
Task Force (Rush et al., 2006). Clinical assessments were
performed blind to genotyping results. Each interview and
diagnostic assignment was reviewed by a senior psychiatrist.
For each patient, all visits were reviewed by the same psychiatrist.

## 2.5 High-throughput sequencing and sequence alignment

Patient DNA samples were sequenced using a targeted panel
of genes involved in mood disorders and ATD metabolism.
Sequencing of the whole *ARRB1* gene was performed,
including exons, introns, and 5'- and 3'-untranslated regions.
The HTS protocol and variant calling methods are as previously
described (Bouali et al., 2017; Chappell et al., 2021).

## 2.6 Variant filtration and selection pipeline

Variant Call Format data were loaded into R (v4.1.0) (R Core
Team, 2020) using the vcfR package (v1.9.0) (Knaus and
Grünwald, 2017). Variant calls were annotated for call quality.
Specifically, variant calls with a sequencing depth (DP) < 20,
SNPs with a quality score (QUAL) < 275, insertions/deletions
(indels) with a QUAL<770, heterozygous calls with an allele
balance (AB) < 0.34 or >0.79, and homozygous calls with an
AB<0.96, were annotated as poor-quality calls. Variants with a
call rate (# poor-quality calls/# of calls) < 95% were removed
(Anderson et al., 2010). Frequent and rare variants were defined
by a MAF≥5% or <5% (corresponding to both low-frequency and
rare variants (Momozawa et al., 2016), but defined here as "rare"
for simplicity), respectively. Frequent variants with a
RegulomeDB category ranking of 1 or 2, corresponding to
variants likely to directly affect TF binding—and in the case
of rank 1, linked to expression of a gene target (Boyle et al.,
2012)—were prioritized for individual analysis in association
with clinical measures.

## 2.7 Variant set analysis

A variant set analysis was performed on all rare variants
passing quality control (QC) using the SMMAT function from
the GMMAT package (v1.3.2), which accounts for repeated
measures in longitudinal data and can analyze both continuous
and binary variables (Chen et al., 2019). The variables to be
explained were the HDRS score, response rate, and remission
rate. A null model was constructed including age, sex, ATD class,
and visit (i.e., time) as covariables. Gaussian distributions with an
identity link function were used to model the HDRS score, while
binomial distributions with a logit link function were used to
model response and remission rates. The genomic relationship
matrix required by the SMMAT function was constructed using
the snpgdsGRM function of the SNPRelate package (v1.26.0)
(Zheng et al., 2012). The SMMAT function was run using
default settings, except that the “use.minor.allele” argument was
set to True—corresponding to the use of the minor allele as the
coding allele rather than the alternative allele—and the
“MAF_range” was set to 0.00-0.50. The SKAT-O test,
corresponding to a linear combination of the burden and
SKAT statistics, was used (Lee et al., 2012). As a single
genomic region was analyzed, a threshold of *p* < 0.05 was
considered significant.

## 2.8 Haplotype analysis and functional annotation of prioritized variants

Hardy-Weinberg equilibrium (HWE) and linkage
disequilibrium (LD) were assessed using Haploview (v4.2) in
the subgroup of Caucasian patients comprising 91% of the patient population (Barrett et al., 2005). Variants were considered to be in LD and assigned to a haplotype block as defined by an r² ≥ 0.8. The genetic variant with the highest MAF from each haplotype block was selected as a proxy for further analysis. Genetic variants were analyzed according to an additive model. Data from RegulomeDB (v2.0.3) and HaploReg (v4.1), including LD information, chromatin state, bound proteins, and altered TF motifs obtained from sources including the 1000 Genomes Project, ENCODE, and JASPAR, were used to annotate variants (Boyle et al., 2012; Ward and Kellis, 2016). LD between analyzed variants and undetected variants or variants outside of our sequencing range were further assessed using LDLink (Machiela and Chanock, 2015).

## 2.9 Data analysis

Statistical analyses were performed in R (v4.1.0) (R Core Team, 2020). Quantitative variables were analyzed using nonparametric Kruskal-Wallis rank-sum tests, while qualitative variables were analyzed using Fisher Exact tests. Mixed-effects models were used to analyze these longitudinal data. Importantly, missing data (e.g., due to dropout) does not result in a complete loss of information as these models allow for an average estimate to be calculated from the non-missing data (Mallinckrodt et al., 2003). Linear mixed-effects models and generalized linear mixed-effects logistic models were constructed with the lme4 package (v1.1-27.1) (Bates et al., 2015). The main variable to explain was the HDRS score. Other variables to explain were response and remission rates. Age, sex, ATD class, and visit were included a priori as fixed-effect covariables in all models. Demographic variables that significantly differed between genotype groups (i.e., *p* < 0.05) were also added as fixed-effects covariables. Individual was added as the random effect to account for repeated measures. Rare variant counts, variant genotype, and its interaction with visit, were the main explanatory variables examined. Rare variant counts were calculated by summing the number of MAF<5% variants present (i.e., heterozygous or homozygous genotype) in each individual. Given the exploratory nature of our study, we considered *p* < 0.05 as significant (Bender and Lange, 2001). Significance of fixed effects was assessed with the Satterthwaite method in linear mixed-effects models and with the Wald chi-square test in generalized linear mixed-effects models. In the event of a significant main or interaction effect, post hoc comparisons were carried out by performing linear and logistic regressions at individual timepoints (i.e., M1, M3, and M6) or by using the emmeans package (Lenth, 2022). Briefly, models were input into the emmeans function (with the “type” argument set to response) to estimate emmeans (for the HDRS) and probabilities (for response and/or remission) according to variant genotypes

<table><caption>TABLE 1 Sociodemographic characteristics. Sociodemographic characteristics of METADAP patients are shown for the whole sample. Kruskal-Wallis tests were used to compare age, tobacco consumption, onset age of MDE, and baseline HDRS scores (presented as mean ± standard deviation (m±sd)). Fisher Exact tests were used to compare sex, socio-education status, ethnicity, smoking status at baseline, MDE recurrence, prescribed ATD drug, and dropout rates across study time (presented as the number of patients and percentage). \*: *p* < 0.05; \*: \*: *p* < 0.01; \*: \*: *p* < 0.001. ATD, antidepressant drug; ECT, electroconvulsive therapy; HDRS, 17-item Hamilton Depression Rating Scale; m, mean; M1, after 1 month of treatment; M3, after 3 months of treatment; M6, after 6 months of treatment; MDE, major depressive episode; n, number of patients; *p*: *p*-value; SNRI, serotonin norepinephrine reuptake inhibitor; sd, standard deviation; SSRI, selective serotonin reuptake inhibitor; TCA, tricyclic antidepressant.</caption><thead><tr><th></th><th></th><th>Total sample n = 388</th></tr></thead><tbody><tr><td>Age (in years) (m±sd)</td><td></td><td>45.4 ± 13.3</td></tr><tr><td>Female [n (%)]</td><td></td><td>263 (68)</td></tr><tr><td rowspan="3">Education [n (%)]</td><td>Primary</td><td>36 (9)</td></tr><tr><td>Secondary</td><td>168 (43)</td></tr><tr><td>Tertiary</td><td>183 (47)</td></tr><tr><td rowspan="3">Ethnicity [n (%)]</td><td>Caucasian</td><td>353 (91)</td></tr><tr><td>African</td><td>24 (6)</td></tr><tr><td>Mixed</td><td>10 (3)</td></tr><tr><td>Current smoker [n (%)]</td><td></td><td>149 (38)</td></tr><tr><td>Pack years (m±sd)</td><td></td><td>15.2 ± 14.9</td></tr><tr><td>Recurrent MDE [n (%)]</td><td></td><td>284 (73)</td></tr><tr><td>Onset age MDE (m±sd)</td><td></td><td>35.3 ± 14.5</td></tr><tr><td>Baseline HDRS (m±sd)</td><td></td><td>24.8 ± 4.9</td></tr><tr><td rowspan="5">Prescribed ATD [n (%)]</td><td>SSRI</td><td>157 (40)</td></tr><tr><td>SNRI</td><td>158 (41)</td></tr><tr><td>TCA</td><td>24 (6)</td></tr><tr><td>Other</td><td>34 (9)</td></tr><tr><td>ECT</td><td>15 (4)</td></tr><tr><td rowspan="3">Missing at follow-up [n (%)]</td><td>M1</td><td>21 (5)</td></tr><tr><td>M3</td><td>121 (31)</td></tr><tr><td>M6</td><td>179 (46)</td></tr></tbody></table>

after controlling for other covariables. Comparisons were performed with the contrast function (with the “method” argument set to pairwise and the “infer” argument set to True) using the emmeans output.

# 3 Results

## 3.1 Patient demographics

Sociodemographic characteristics for the whole cohort are shown in Table 1. The mean age was 45.4, 68% of subjects were women, and 91% of subjects were Caucasian. 38% were current smokers at baseline and 73% had previously experienced an MDE.
<figure><img src="image_5.png" /></figure>

## 3.2 Genetic variant selection and linkage analysis

Genetic variant selection is summarized in Figure 3. A total of 953 unique genetic variants were identified from extracted Variant Call data, of which 832 were SNPs and 121 were indels (see Supplementary Table S1). 195 variants with call rates <95%

(119 SNPs and 76 indels) were removed. Thus, 758 variants were available for analysis. 643 variants were rare (i.e., MAF<5%) (602 SNPs and 41 indels), of which 610 had a MAF<1% and 415 were identified in 1 individual (see Supplementary Table S1). Among the 115 frequent variants (i.e., MAF≥5%) (111 SNPs and 4 indels), 12 SNPs had a RegulomeDB category rank of 1 or 2 [i.e., linked to expression of a gene target (rank 1) and likely to directly affect TF binding (ranks 1 and 2)] and were prioritized for individual analysis (see Figure 3). Information for each prioritized SNP is shown in Table 2. Prioritized SNPs did not significantly deviate from HWE. Ten of the twelve prioritized SNPs were intronic, eight of which were located in intron 1.

Among the 12 prioritized SNPs, 3 pairs were observed in LD. rs501372 and rs567807 comprised haplotype block 1 (ht1), rs553664 and rs506233 (ht2), and rs569796 and rs504683 (ht3) (see Table 2 and Supplementary Figure S1). rs501372, rs553664, and rs504683 were selected as proxy SNPs for further analysis from each haplotype block. Nine SNPs were thus analyzed.

## 3.3 Variant set and accumulation analyses of rare variants

SMMAT was used to assess the association of the 643 rare (i.e., MAF<5%) variants with clinical measures over time. The 643 rare variants were significant in explaining changes in the HDRS score over time in burden ($T = 878.9; \Phi = 89194.6; p = 0.0033$) and SKAT-O tests ($p = 0.0069$) and remission rates in the burden test ($T = -1974.1; \Phi = 869281.1; p = 0.034$), but not response rates, though a statistical trend was observed in the burden test ($T = -1773.9; \Phi = 937912.7; p = 0.07$) (see Table 3). To further

<table><caption>TABLE 2 Genomic and other information about the 12 prioritized frequent SNPs. For the prioritized SNPs (i.e., MAF≥5% and RegulomeDB category rank of 1 or 2), the rs#, genomic position relative to the human genome assembly, hg19, and genomic region are shown. RegulomeDB rankings correspond to likely to affect binding and linked to expression of a gene target (1b, 1f) and likely to affect binding (2b, 2c). The major and minor alleles and minor allele frequency (MAF (%)) with respect to the Caucasian subpopulation are also shown. Hardy-Weinberg equilibrium p-values (HWE), haplotype blocks (Haplotype block), and linkage disequilibrium r² values (LD (r²)) relative to the Caucasian subpopulation are also given. HWE, Hardy-Weinberg equilibrium; MAF, minor allele frequency; LD, linkage disequilibrium.</caption><thead><tr><th>rs #</th><th>Position</th><th>Region</th><th>RegulomeDB rank</th><th>Major</th><th>Minor</th><th>MAF (%)</th><th>HWE</th><th>Haplotype block</th><th>LD (r²)</th></tr></thead><tbody><tr><td>rs2279130</td><td>74977185</td><td>3'UTR</td><td>2b</td><td>C</td><td>T</td><td>7.2</td><td>0.5</td><td></td><td></td></tr><tr><td>rs501372</td><td>74993593</td><td>intron 5</td><td>2b</td><td>C</td><td>A</td><td>38.2</td><td>0.24</td><td>ht1</td><td>0.96</td></tr><tr><td>rs877711</td><td>74994352</td><td>exon5</td><td>2b</td><td>G</td><td>A</td><td>11.8</td><td>0.18</td><td></td><td></td></tr><tr><td>rs567807</td><td>75000189</td><td>intron 2</td><td>2b</td><td>A</td><td>G</td><td>37.2</td><td>0.13</td><td>ht1</td><td></td></tr><tr><td>rs553664</td><td>75005209</td><td>intron 1</td><td>2b</td><td>G</td><td>A</td><td>44.4</td><td>0.58</td><td>ht2</td><td>0.98</td></tr><tr><td>rs506233</td><td>75012535</td><td>intron 1</td><td>1f</td><td>A</td><td>G</td><td>43.7</td><td>0.38</td><td>ht2</td><td></td></tr><tr><td>rs536852</td><td>75017436</td><td>intron 1</td><td>2b</td><td>G</td><td>A</td><td>49.1</td><td>0.98</td><td></td><td></td></tr><tr><td>rs1676887</td><td>75019296</td><td>intron 1</td><td>1f</td><td>G</td><td>A</td><td>10.4</td><td>0.38</td><td></td><td></td></tr><tr><td>rs113636971</td><td>75033020</td><td>intron 1</td><td>2b</td><td>A</td><td>G</td><td>6.4</td><td>0.88</td><td></td><td></td></tr><tr><td>rs569796</td><td>75050842</td><td>intron 1</td><td>1b</td><td>G</td><td>C</td><td>29.9</td><td>0.48</td><td>ht3</td><td>0.98</td></tr><tr><td>rs504683</td><td>75051233</td><td>intron 1</td><td>2c</td><td>A</td><td>G</td><td>30.2</td><td>0.39</td><td>ht3</td><td></td></tr><tr><td>rs561923</td><td>75058004</td><td>intron 1</td><td>2b</td><td>C</td><td>A</td><td>25.4</td><td>0.8</td><td></td><td></td></tr></tbody></table>
<table><caption>TABLE 3 SMMAT results. Data and results for variant set analyses of the association between rare (i.e., MAF<5%) variants passing quality control and HDRS scores, response rates, and remission rates. All models included age, sex, antidepressant class, and visit as covariables. Shown are the mean allele frequency, the burden test score statistic and variance, and the associated *p*-values of the burden, SKAT, and SKAT-O tests. bold text*: *p* < 0.05. AF, allele frequency; *p*: *p*-value.</caption><thead><tr><th>Clinical measure</th><th>Minimum frequency</th><th>Mean frequency</th><th>Maximum frequency</th><th>Burden stat</th><th>Burden variance</th><th>Burden p</th><th>SKAT p</th><th>SKAT-O p</th></tr></thead><tbody><tr><td>HDRS</td><td>0.00129</td><td>0.00442</td><td>0.5</td><td>878.9</td><td>89194.6</td><td><strong>0.0033*</strong></td><td>0.41</td><td><strong>0.0069*</strong></td></tr><tr><td>Response</td><td></td><td></td><td></td><td>-1773.9</td><td>937912.7</td><td>0.07</td><td>0.47</td><td>0.12</td></tr><tr><td>Remission</td><td></td><td></td><td></td><td>-1974.1</td><td>869281.1</td><td><strong>0.034*</strong></td><td>0.37</td><td>0.06</td></tr></tbody></table>

<table><caption>TABLE 4 Mixed-effects model results. Associations for each of the clinical measures (HDRS, response, and remission) with 1) the rare (i.e., MAF<5%) variant count and its interaction with visit (i.e., time) and 2) the 9 prioritized frequent (i.e., MAF≥5% and RegulomeDB rank of 1 or 2) SNPs analyzed and their interactions with visit. Significant (i.e., *p* < 0.05) associations are highlighted in dark gray and bolded. Significance of fixed effects was assessed with the Satterthwaite method in mixed-effects models of the HDRS score and with the Wald chi-square test in mixed-effects models of response and remission. All models were adjusted for age, sex, antidepressant class, and any sociodemographic characteristics that significantly (bold text and *: *p* < 0.05) differed between genotypes (see Supplementary Table S2).</caption><thead><tr><th rowspan="2"></th><th colspan="6">HDRS</th><th colspan="3">Response</th><th colspan="3">Remission</th></tr><tr><th>SS</th><th>MS</th><th>df (num)</th><th>df (den)</th><th>F</th><th>p</th><th>χ<sup>2</sup></th><th>df</th><th>p</th><th>χ<sup>2</sup></th><th>df</th><th>p</th></tr></thead><tbody><tr><td colspan="13"><strong>MAF&lt;5% count</strong></td></tr><tr><td>MAF&lt;5% count</td><td>804.264</td><td>804.264</td><td>1</td><td>438.056</td><td>26.87</td><td><strong>&lt;0.001*</strong></td><td>5.81</td><td>1</td><td><strong>0.016*</strong></td><td>4.56</td><td>1</td><td><strong>0.033</strong></td></tr><tr><td>Visit:MAF&lt;5% count</td><td>348.786</td><td>116.262</td><td>3</td><td>918.602</td><td>3.88</td><td><strong>0.0090*</strong></td><td>0.53</td><td>2</td><td>0.77</td><td>1.99</td><td>2</td><td>0.37</td></tr><tr><td colspan="13"><strong>Prioritized variants (MAF≥5% and RegulomeDB rank ≤ 2)</strong></td></tr><tr><td>rs2279130</td><td>59.196</td><td>29.598</td><td>2</td><td>401.121</td><td>0.98</td><td>0.38</td><td>1.48</td><td>2</td><td>0.48</td><td>4.11</td><td>2</td><td>0.13</td></tr><tr><td>Visit:rs2279130</td><td>122.998</td><td>20.500</td><td>6</td><td>914.026</td><td>0.68</td><td>0.67</td><td>3.03</td><td>4</td><td>0.55</td><td>1.94</td><td>4</td><td>0.75</td></tr><tr><td>rs501372</td><td>5.277</td><td>2.638</td><td>2</td><td>371.393</td><td>0.09</td><td>0.92</td><td>1.16</td><td>2</td><td>0.56</td><td>1.50</td><td>2</td><td>0.47</td></tr><tr><td>Visit:rs501372</td><td>168.779</td><td>28.130</td><td>6</td><td>882.364</td><td>0.93</td><td>0.47</td><td>4.64</td><td>4</td><td>0.33</td><td>2.00</td><td>4</td><td>0.73</td></tr><tr><td>rs877711</td><td>29.797</td><td>14.898</td><td>2</td><td>411.229</td><td>0.50</td><td>0.61</td><td>1.17</td><td>2</td><td>0.56</td><td>1.10</td><td>2</td><td>0.58</td></tr><tr><td>Visit:rs877711</td><td>195.206</td><td>32.534</td><td>6</td><td>942.398</td><td>1.09</td><td>0.36</td><td>0.63</td><td>4</td><td>0.96</td><td>4.59</td><td>4</td><td>0.33</td></tr><tr><td>rs553664</td><td>147.286</td><td>73.643</td><td>2</td><td>399.038</td><td>2.51</td><td>0.08</td><td>0.49</td><td>2</td><td>0.78</td><td>3.76</td><td>2</td><td>0.15</td></tr><tr><td>Visit:rs553664</td><td>541.248</td><td>90.208</td><td>6</td><td>919.440</td><td>3.08</td><td><strong>0.0055*</strong></td><td>8.86</td><td>4</td><td>0.06</td><td>11.02</td><td>4</td><td><strong>0.026*</strong></td></tr><tr><td>rs536852</td><td>213.512</td><td>106.756</td><td>2</td><td>385.710</td><td>3.53</td><td><strong>0.030*</strong></td><td>0.85</td><td>2</td><td>0.66</td><td>3.85</td><td>2</td><td>0.15</td></tr><tr><td>Visit:rs536852</td><td>390.451</td><td>65.075</td><td>6</td><td>886.704</td><td>2.15</td><td><strong>0.046*</strong></td><td>7.11</td><td>4</td><td>0.13</td><td>12.84</td><td>4</td><td><strong>0.012*</strong></td></tr><tr><td>rs1676887</td><td>50.854</td><td>25.427</td><td>2</td><td>372.885</td><td>0.85</td><td>0.43</td><td>1.36</td><td>2</td><td>0.51</td><td>1.78</td><td>2</td><td>0.41</td></tr><tr><td>Visit:rs1676887</td><td>194.069</td><td>32.345</td><td>6</td><td>887.639</td><td>1.08</td><td>0.37</td><td>1.71</td><td>4</td><td>0.79</td><td>5.43</td><td>4</td><td>0.25</td></tr><tr><td>rs113636971</td><td>10.977</td><td>5.489</td><td>2</td><td>379.137</td><td>0.18</td><td>0.83</td><td>0.23</td><td>2</td><td>0.89</td><td>0.14</td><td>2</td><td>0.93</td></tr><tr><td>Visit:rs113636971</td><td>94.369</td><td>15.728</td><td>6</td><td>896.752</td><td>0.52</td><td>0.79</td><td>0.57</td><td>4</td><td>0.97</td><td>0.74</td><td>4</td><td>0.95</td></tr><tr><td>rs504683</td><td>64.300</td><td>32.150</td><td>2</td><td>373.296</td><td>1.05</td><td>0.35</td><td>2.47</td><td>2</td><td>0.29</td><td>2.47</td><td>2</td><td>0.29</td></tr><tr><td>Visit:rs504683</td><td>74.684</td><td>12.447</td><td>6</td><td>892.312</td><td>0.41</td><td>0.87</td><td>0.79</td><td>4</td><td>0.94</td><td>0.79</td><td>4</td><td>0.94</td></tr><tr><td>rs561923</td><td>17.422</td><td>8.711</td><td>2</td><td>366.270</td><td>0.28</td><td>0.75</td><td>1.03</td><td>2</td><td>0.60</td><td>0.31</td><td>2</td><td>0.86</td></tr><tr><td>Visit:rs561923</td><td>37.649</td><td>6.275</td><td>6</td><td>883.377</td><td>0.20</td><td>0.98</td><td>2.06</td><td>4</td><td>0.72</td><td>2.62</td><td>4</td><td>0.62</td></tr></tbody></table>

χ², chi-square value; den, denominator; df, degrees of freedom; F, F-value; HDRS, 17-item Hamilton Rating Depression scale; MAF, minor allele frequency; MS, mean sum of squares; num, numerator; *p*: *p*-value; SS, sum of squares; ":" designates an interaction between model terms.

analyze the associations of rare variant accumulation with clinical measures, rare variant counts were calculated and analyzed. After controlling for age, sex, ATD class, and visit, a significant main effect of the rare variant count was observed in the mixed-effects models of response [χ² (df = 1, *n* = 388) = 5.81, *p* = 0.016] and remission [χ²

(df = 1, *n* = 388) = 4.56, *p* = 0.033], as was a significant interaction with visit in the mixed-effects model of the HDRS score (F<sub>3,918</sub> = 3.88, *p* = 0.0090) (see Table 4). In simplified mixed-effects models (i.e., without interaction) controlling for the same factors, the rare variant count was indeed significant in models of response (odds
<figure><img src="image_6.png" /></figure>

**FIGURE 4**
Association of MAF < 5% variant count with clinical measures. Model estimates adjusted for the mean age (45.36), sex (male), and antidepressant class (SSRI) of (A) the probability of response, (B) the probability of remission, (C) the HDRS score at M1, and (D) the HDRS score at M3 (y-axis) according to the MAF < 5% variant count (x-axis). HDRS, 17-item Hamilton Depression Rating Scale; M1, 1 month after beginning antidepressant treatment; M3, 3 months after beginning antidepressant treatment; MAF, minor allele frequency.

ratio (OR) = 0.95; 95% confidence interval (95%CI) [0.91–0.99]; *p* = 0.016) and remission (OR = 0.92; 95%CI [0.86–0.99]; *p* = 0.022) (see Figures 4A,B). Following Bonferroni correction, it was also a significant factor explaining the HDRS score in multiple linear regressions at M1 (coefficient = 0.26; 95%CI [0.11–0.41]; *p* = 0.0021) and M3 (coefficient = 0.46; 95%CI [0.28–0.64]; *p* < 0.001) (see Figures 4C,D), but not at M6 (coefficient = 0.25; 95%CI [−0.01–0.51]; *p* = 0.19).

## 3.4 Clinical response according to prioritized variants

Sociodemographic characteristics and significant differences between prioritized SNP genotypes are shown in Supplementary Table S2. Mixed-effects model associations of clinical measures with genotype and genotype × visit interactions for the

9 prioritized SNPs are summarized in Table 4. In mixed-effects models of the HDRS score controlling for age, sex, ATD class, ethnicity, and smoking status—and baseline HDRS scores for rs553664—the genotype × visit interaction was a significant explanatory factor for rs553664 (*F*<sub>6,919</sub> = 3.08, *p* = 0.0055) and rs536852 (*F*<sub>6,886</sub> = 2.15, *p* = 0.046) (see Table 4). In mixed-effects models of remission controlling for these same variables, the genotype × visit interaction was also significant for rs553664 [*χ*<sup>2</sup> (df = 4, *n* = 388) = 11.02, *p* = 0.026] and rs536852 [*χ*<sup>2</sup> (df = 4, *n* = 388) = 12.84, *p* = 0.012]. Following Bonferroni correction, a significantly higher HDRS score was observed in rs553664 AA homozygotes (13.98 ± 1.06) compared to AG heterozygotes (10.59 ± 0.86) at M6 (coefficient = 3.39, 95%CI [1.29–5.49], *p* = 0.014) (see Figure 5A) and in rs536852 GG homozygotes (14.88 ± 1.10) compared to AG heterozygotes (11.26 ± 0.95) at M6 (coefficient = 3.62, 95%CI [1.54–5.71], *p* = 0.0061) (see Figure 5B). A significantly lower remitter rate
<figure><img src="image_7.png" /><figcaption>**FIGURE 5**
Average HDRS scores according to rs553664 and rs536852 genotypes over the course of antidepressant treatment. Average HDRS scores (y-axis) according to rs553664 (A) and rs536852 (B) genotypes (see legend) are shown across time (x-axis) after controlling for age, sex, antidepressant class, ethnicity, and significantly different demographic factors (see Supplementary Table S3). Error bars correspond to the standard error estimates from mixed-effects models. \*: *p* < 0.0056 (0.05/9) \*: \*: *p* < 0.0011 (0.01/9). HDRS, 17-item Hamilton Depression Rating Scale; M0, baseline, prior to beginning antidepressant treatment; M1, 1 month after beginning antidepressant treatment; M3, 3 months after beginning antidepressant treatment; M6, 6 months after beginning antidepressant treatment.</figcaption></figure>

was also observed in rs536852 GG homozygotes (8%, *n* = 56) compared to AG heterozygotes (42%, *n* = 105) at M6 (OR = 0.12, 95%CI [0.39–0.37], *p* = 0.0018) (see Figure 6). No significant differences in response rates, or in remitter rates between rs553664 genotypes, were observed.

## 3.5 Functional analysis of ARRB1 polymorphisms

Annotations from RegulomeDB (v2.0.3) and HaploReg (v4.1) for rs553664 and rs536852—associated with clinical

<figure><img src="image_8.png" /></figure>

measures—and any variants in LD with them are shown in Table 5. rs553664 was observed to be in LD with rs506233 and an intronic deletion, rs35731045. rs536852 was not observed to be in LD with any variants within ±500,000 base pairs of its genomic position in non-Finnish European populations. It was most closely linked to rs553664 and rs506233, though these associations (*r*² = 0.67 and *r*² = 0.68, respectively) were weaker than in our sample (see Supplementary Figures S1, S2). rs553664, rs506233, and rs35731045 were annotated as binding FOS, NFIC, USF1, and IKZF1. rs536852 was annotated as binding POLR2A, FOS, and JUN. Altered TF binding motifs included the CCCTC-binding factor (CTCF). Each variant had at least one annotation with histone marks/chromatin state (see Table 5).

## 4 Discussion

The present study is an ancillary investigation of the METADAP cohort that aimed to analyze the association between genetic variants of ARRB1 and response following ATD treatment. Overall, 758 genetic variants passed quality control measures and were further filtered based on allelic
<table><caption>TABLE 5 HaploReg and RegulomeDB annotations. For rs553664 and rs536852, the rs#’s in linkage disequilibrium (i.e., r² > 0.8) with the given SNP, the associated histone marks/chromatin state in brain tissues, the bound proteins, and known and discovered altered motifs, are given.</caption><thead><tr><th>rs #</th><th>rs #’s in LD (r² > 0.8)<sup>a</sup></th><th>Histone marks<sup>a</sup>/ chromatin state<sup>b</sup> (brain)</th><th>Bound proteins</th><th>Altered TF motifs</th></tr></thead><tbody><tr><td>rs553664</td><td>rs506233, rs35731045</td><td>Strong transcription<sup>b</sup></td><td>NFYA<sup>a</sup>, USF1<sup>a</sup>, USF1<sup>b</sup>, USF2<sup>b</sup>, RAD51<sup>b</sup>, NFYC<sup>b</sup>, FOS<sup>b</sup>, EMSY<sup>b</sup>, CBFA2T2<sup>b</sup>, ZNF316<sup>b</sup></td><td>Brachyury_1<sup>a</sup>, HP1-site-factor<sup>a</sup>, IRX3<sup>b</sup>, TBXT<sup>b</sup></td></tr><tr><td>rs506233</td><td>rs553664, rs35731045</td><td>Enhancer<sup>a</sup>, promoter<sup>a</sup>, enhancers<sup>b</sup>, strong transcription<sup>b</sup></td><td>NFIC<sup>b</sup></td><td>Myc Known8<sup>a</sup></td></tr><tr><td>rs35731045</td><td>rs553664, rs506233</td><td>Strong transcription<sup>b</sup></td><td>IKZF1<sup>b</sup></td><td>ERalpha-a_disc1<sup>a</sup>, ERalpha-a_disc2<sup>a</sup>, ERalpha-a known1<sup>a</sup>, ERalpha-a known4<sup>a</sup>, Esr2, RAR<sup>a</sup>, RXRA_disc1<sup>a</sup>, RXRAKnown5<sup>a</sup>, SF1<sup>a</sup>, STAT_disc3<sup>a</sup></td></tr><tr><td>rs536852</td><td></td><td>Enhancer<sup>a</sup>, promoter<sup>a</sup>, enhancers<sup>b</sup>, strong transcription<sup>b</sup></td><td>CFOS<sup>a</sup>, CJUN<sup>a</sup>, FOS<sup>b</sup>, JUN<sup>b</sup>, POLR2A<sup>b</sup></td><td>CTCF_disc9<sup>a</sup>, Irf_disc4<sup>a</sup>, Pax-4_3<sup>a</sup>, SP1Known2<sup>a</sup>, TATA_disc4<sup>a</sup>, TR4_disc3<sup>a</sup>, PAX4<sup>b</sup></td></tr></tbody></table>

frequency and their functional annotations as likely to affect TF binding. Analyses were carried out to assess the association of clinical response following ATD treatment with 1) the cumulative effect of rare (i.e., MAF<5%) variants and 2) frequent (i.e., MAF≥5%) variants with likely functional consequences. Twelve SNPs were prioritized, of which three pairs were observed to be in LD. Thus, 9 SNPs (3 proxies and 6 individual) were analyzed.

The variant set analysis using SMMAT (Chen et al., 2019) and leveraging the SKAT-O framework (Lee et al., 2012), suggests that the accumulation of rare ARRB1 variants is associated with HDRS score changes and remission rates over the course of ATD treatment, as the burden test results were significant. As such, increased variant accumulation in ARRB1 may have an impact on clinical response following ATD treatment. Furthermore, as reflected by the significant finding in the burden test, but not the SKAT test, a large proportion of these variants likely have causal effects on ATD response that, overall, act in the same direction (Lee et al., 2012). This was further highlighted by the significant association of greater rare variant counts with higher HDRS scores—notably at M1 and M3—and overall lower response and remission rates during treatment. Of note, certain rare variants in our cohort, including those in LD with the rs12274033 polymorphism that was previously associated with remission rates following mirtazapine treatment in a Korean population (Chang et al., 2015), may contribute to the observed associations with HDRS score changes and remission rates.

The 9 prioritized SNPs we analyzed have not been previously examined. They were prioritized for analysis from over 900 variants identified from HTS data spanning the ARRB1 gene, including exons, introns, and 5'- and 3'-untranslated regions. Given that ARRB1 mRNA and β-arrestin 1 protein levels have been previously associated with depression and

ATD treatment (Avissar et al., 2004; Matuzany-Ruban et al., 2005; Golan et al., 2013), variants annotated as likely to lie in a functional element—and thus be more likely to impact ARRB1 expression or β-arrestin function—were prioritized using RegulomeDB (Boyle et al., 2012). RegulomeDB has previously been used to select and study genetic variants in various disease backgrounds (Lee et al., 2015; Hong et al., 2018; Huo et al., 2019). Specifically, variants with a category ranking of 1 or 2, corresponding to variants likely to directly affect TF binding, were prioritized for further analysis. To our knowledge, this is the first study to leverage functional annotations to prioritize variants for analysis in the context of MDD and its treatment with ATDs.

Protein binding annotations from HaploReg and RegulomeDB, utilizing ChIP-seq data from ENCODE, suggest that rs553664 (and its proxies, including rs506233) and rs536852, associated with the HDRS score and remission in METADAP, lie within functional elements that bind various proteins including USF1, NFIC, FOS, JUN, and POLR2A. Genetic risk factors that overlap with TF-binding sites for USF1, NFIC, and POLR2A have been identified in relation to schizophrenia and MDD (Huo et al., 2019; Li et al., 2020).

FOS and JUN can dimerize to form a member of the Activator Protein-1 (AP-1) TF family, a leucine-zipper TF shown to regulate proliferation, differentiation, and apoptosis (Garces de los Fayos Alonso et al., 2018). Both FOS and JUN are expressed in the brain where they have context-dependent influences on neurodegeneration and/or neuroprotection (Herdegen and Waetzig, 2001). Recently, AP-1 was observed to have a role in fluoxetine response in a murine model of heightened anxiety (Chottekalapanda et al., 2020), and was also found to be activated by other factors such as BDNF, which is itself a factor influencing the response to different ATD treatments (Colle et al., 2015). It is known that ATDs can

<footnote><sup>a</sup>Annotation from HaploReg v4.1.</footnote>

<footnote><sup>b</sup>Annotation from RegulomeDB, v2.0.3. LD, linkage disequilibrium; TF, transcription factor.</footnote>
modify the activation of different genes through their modulation of serotonergic pathways (Lesch and Heils, 2000). Indeed, several TFs have been associated with ATD response or have been shown to be altered in depression (Pittenger and Duman, 2008; Sasayama et al., 2013). Moreover, the beta-arrestins act as scaffolds for kinases, including c-Jun N-terminal kinases of the JNK axis, which are implicated in neuronal plasticity, memory, and neuronal maturation (Herdegen and Waetzig, 2001; Shukla et al., 2011). As such, the potential impact of rs53664 and/or rs536852 may have on AP-1 and/or the binding of other TFs could alter β-arrestin 1 expression and, consequently, any downstream functions of β-arrestin 1, including the transmission of receptor-mediated signals during ATD treatment or the scaffolding of proteins.

A putative binding motif for CTCF was also associated with rs536852. CTCF is a conserved zinc-finger TF with important roles in transcriptional regulation and 3D genome organization (Li et al., 2020). In a study of risk alleles for schizophrenia, SNPs lying in CTCF binding sites were enriched among those analyzed (Huo et al., 2019). Furthermore, in a study of MDD genetic risk factors, 11 of the 34 TF binding-disrupting SNPs analyzed were specific to CTCF binding sites, suggesting that disruption of CTCF binding may be a shared mechanism among genetic risk variants of MDD (Li et al., 2020). ATD-induced transcriptional alterations in the brain may also be influenced by CTCF (Piechota et al., 2015). Of course, whether CTCF binding is disrupted and, if so, the biological and clinical consequences of this, should be verified and further investigated.

Our findings demonstrate an association between genetic variants of *ARRB1* and clinical outcomes of response following ATD treatment. These findings suggest that these genetic factors might influence various mechanisms involved in the response to ATD therapy. Given that biased agonists can preferentially activate or bypass β-arrestin 1-mediated processes (Violin and Lefkowitz, 2007; Cottingham et al., 2011; Shukla et al., 2011; Bond et al., 2019), patients carrying either an excess of rare *ARRB1* variants or the rs553664 and rs536852 polymorphisms may benefit from the use of ATD treatments that 1) favor G-protein-mediated pathways and 2) circumvent β-arrestin-mediated pathways and the desensitization of GPCRs. Alternatively, these patients could also benefit from the addition of other therapeutical strategies, such as psychotherapies (i.e., talk therapies) or brain stimulation therapies, including ECT and repetitive transcranial magnetic stimulation.

Our study has several limitations. First, the proportion of missing data was relatively high, though the proportion at 6 months (46.1%) was comparable to the proportions in the STAR*D cohort at 9 (43%) and 12 weeks (61%) (Trivedi et al., 2006; Muthén et al., 2011). Additionally, mixed-effects models are a robust method to control for the bias imposed by these

missing data (Mallinckrodt et al., 2003). Second, hard filters were used to filter out poor-quality variants, which could have removed variants near the cutoff threshold. Importantly, these variants could have had significant associations with clinical measures. Third, it should be noted that the functions of β-arrestin 1 and 2 often overlap. Thus, any impact that variants of *ARRB1* may have on β-arrestin 1 expression or function may be compensated for by β-arrestin 2, at least within the cytoplasm (Shukla et al., 2011). The strengths of this study include its prospective and naturalistic design, which allows for the analysis of treatment response across a 6-month period and better reflects “real-world” clinical practice. Additionally, this study leveraged functional annotations to select and prioritize variants with likely functional consequences on *ARRB1* expression for analysis since reduced β-arrestin 1 levels have been previously associated with depression severity and the response to ATD treatment (Avissar et al., 2004; Matuzany-Ruban et al., 2005; Golan et al., 2013). Finally, to date, this is the largest cohort analysis of *ARRB1* variants in relation to response following ATD treatment in patients with MDD. The findings from this exploratory analysis should be replicated in independent cohorts of depressed individuals undergoing ATD treatment.

To conclude, *ARRB1* genetic variants were associated with clinical measures of response following ATD treatment in a cohort of ATD-treated depressed patients. Specifically, we observed associations between clinical measures and 1) rare variant accumulation and 2) two frequent variants with probable functional consequences, rs553664 and rs536852. These findings suggest that variants of *ARRB1*, and perhaps especially those within genomic regulatory regions, may contribute to clinical response following ATD treatment. Functional analyses would help confirm the potential effect(s) that these polymorphisms have on clinical improvement following antidepressant treatment. Further study in larger cohorts would also help corroborate the findings of our exploratory analysis.

## Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: data are under the protection of health data regulation set by the French National Commission on Informatics and Liberty (Commission Nationale de l'Informatique et des Libertés, CNIL), in line with European regulations and the Data Protection Act, and the Comité de protection des personnes (CPP, equivalent to the Research Ethics Committee). The data can be available upon reasonable request to the principal investigator of the study (emmanuelle.corruble@aphp.fr). The French law forbids us to provide free access to METADAP data. Please feel free to contact us should you have any additional questions.
## Ethics statement

The studies involving human participants were reviewed and approved by Ethics Committee of Paris-Boulogne. The patients/ participants provided their written informed consent to participate in this study.

## Author contributions

KC, CV, and JB designed the study. EC, BF, LB, RC, DD, ST, and FG acquired the data. KC, AA, and KE-A analyzed the data. KC and CV wrote the article. All authors contributed to manuscript revision and read and approved the submitted version.

## Funding

The METADAP study was funded by a national grant (PHRC, AOM06022) and sponsored by the Assistance Publique-Hôpitaux de Paris (AP-HP) (ClinicalTrials.gov Identifier: NCT00526383).

## Acknowledgments

We would like to thank the participants of the METADAP study, the staff of the Centre de Recherche en Epidemiologie et

## References

Anderson, C. A., Pettersson, F. H., Clarke, G. M., Cardon, L. R., Morris, A. P., and Zondervan, K. T. (2010). Data quality control in genetic case-control association studies. *Nat. Protoc.* 5 (9), 1564–1573. doi:10.1038/nprot.2010.116

Avissar, S., Matuzany-Ruban, A., Tzukert, K., and Schreiber, G. (2004). [[bgr]]-Arrestin-1 levels: Reduced in leukocytes of patients with depression and elevated by antidepressants in rat brain. *Am. J. Psychiatry* 161 (11), 2066–2072. doi:10.1176/appi.ajp.161.11.2066

Barrett, J. C., Fry, B., Maller, J., and Daly, M. J. (2005). Haploview: Analysis and visualization of LD and haplotype maps. *Bioinformatics* 21 (2), 263–265. doi:10.1093/bioinformatics/bth457

Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. *J. Stat. Softw.* 67 (11), 1–48. doi:10.18637/jss.v067.i01

Bender, R., and Lange, S. (2001). Adjusting for multiple testing—When and how? *J. Clin. Epidemiol.* 54 (4), 343–349. doi:10.1016/S0895-4356(00)00314-0

Bond, R. A., Lucero Garcia-Rojas, E. Y., Hegde, A., and Walker, J. K. L. (2019). Therapeutic potential of targeting β-arrestin. *Front. Pharmacol.* 10, 124. doi:10.3389/fphar.2019.00124

Bouali, N., Francou, B., Bouligand, J., Imanci, D., Dimassi, S., Tosca, L., et al. (2017). New MCM8 mutation associated with premature ovarian insufficiency and chromosomal instability in a highly consanguineous Tunisian family. *Fertil. Steril.* 108 (4), 694–702. doi:10.1016/j.fertnstert.2017.07.015

Boyle, A. P., Hong, E. L., Hariharan, M., Cheng, Y., Schaub, M. A., Kasowski, M., et al. (2012). Annotation of functional variation in personal genomes using RegulomeDB. *Genome Res.* 22 (9), 1790–1797. doi:10.1101/gr.137323.112

Calabrese, G., Sallese, M., Stornaiuolo, A., Palka, G., and De BlAsi, A. (1994). Assignment of the β-arrestin 1 gene (ARRB1) to human chromosome 11q13. *Genomics* 24 (1), 169–171. doi:10.1006/geno.1994.1594

Santé des Populations, the members of the Psychiatry Department at Bicêtre Hospital and the Centre de Ressources Biologiques Paris-Saclay (CRB Paris-Saclay) for the storage of METADAP biological samples (BRIF number: BB-0033-00089).

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.974570/full#supplementary-material

Chang, H. S., Won, E. S., Lee, H-Y., Ham, B. J., Kim, Y. G., and Lee, M. S. (2015). Association of ARRB1 polymorphisms with the risk of major depressive disorder and with treatment response to mirtazapine. *J. Psychopharmacol.* 29 (5), 615–622. doi:10.1177/0269881114554273

Chappell, K., Francou, B., Habib, C., Huby, T., Leoni, M., Cottin, A., et al. (2021). Galaxy is a suitable bioinformatics platform for the molecular diagnosis of human genetic disorders using high-throughput sequencing data analysis. Five years of experience in a clinical laboratory. *Clin. Chem.* 2021, 313–321. doi:10.1093/clinchem/hvab220

Chen, H., Huffman, J. E., Brody, J. A., Wang, C., Lee, S., Li, Z., et al. (2019). Efficient variant set mixed model association tests for continuous and binary traits in large-scale whole-genome sequencing studies. *Am. J. Hum. Genet.* 104 (2), 260–274. doi:10.1016/j.ajhg.2018.12.012

Chottekalapanda, R. U., Kalik, S., Gresack, J., Ayala, A., Gao, M., Wang, W., et al. (2020). AP-1 controls the p11-dependent antidepressant response. *Mol. Psychiatry* 25 (7), 1364–1381. doi:10.1038/s41380-020-0767-8

Colle, R., Gressier, F., Verstuyft, C., Deflesselle, E., Lepine, J. P., Ferreri, F., et al. (2015). Brain-derived neurotrophic factor Val66Met polymorphism and 6-month antidepressant remission in depressed Caucasian patients. *J. Affect. Disord.* 175, 233–240. doi:10.1016/j.jad.2015.01.013

Corruble, E., El Asmar, K., Trabado, S., Verstuyft, C., Falissard, B., Colle, R., et al. (2015). Treating major depressive episodes with antidepressants can induce or worsen metabolic syndrome: Results of the METADAP cohort. *World Psychiatry* 14 (3), 366–367. doi:10.1002/wps.20260

Cottingham, C., Chen, Y., Jiao, K., and Wang, Q. (2011). The antidepressant desipramine is an arrestin-biased ligand at the α2a-adrenergic receptor driving receptor down-regulation *in vitro* and *in vivo*. *J. Biol. Chem.* 286 (41), 36063–36075. doi:10.1074/jbc.M111.261578
Garces de los Fayos Alonso, I., Liang, H.-C., Turner, S. D., Lagger, S., Merkel, O., and Kenner, L. (2018). The role of activator protein-1 (AP-1) family members in CD30-positive lymphomas. *Cancers* 10 (4), 93. doi:10.3390/cancers10040093

Golan, M., Schreiber, G., and Avissar, S. (2013). Antidepressant-induced differential ubiquitination of β-arrestins 1 and 2 in mononuclear leucocytes of patients with depression. *Int. J. Neuropsychopharmacol.* 16 (8), 1745–1754. doi:10.1017/S1461145713000291

Golan, M., Schreiber, G., and Avissar, S. (2009). Antidepressants, β-arrestins and GRKs: From regulation of signal desensitization to intracellular multifunctional adaptor functions. *Curr. Pharm. Des.* 15 (14), 1699–1708. doi:10.2174/138161209788168038

Hamilton, M. (1960). A rating scale for depression. *J. Neurol. Neurosurg. Psychiatry* 23 (1), 56–62. doi:10.1136/jnnp.23.1.56

Herdegen, T., and Waetzig, V. (2001). AP-1 proteins in the adult brain: Facts and fiction about effectors of neuroprotection and neurodegeneration. *Oncogene* 20 (19), 2424–2437. doi:10.1038/sj.onc.1204387

Hong, M. J., Yoo, S. S., Choi, J. E., Kang, H. G., Do, S. K., Lee, J. H., et al. (2018). Functional intronic variant of SLC5A10 affects DRG2 expression and survival outcomes of early-stage non-small-cell lung cancer. *Cancer Sci.* 109 (12), 3902–3909. doi:10.1111/cas.13814

Huo, Y., Li, S., Liu, J., Li, X., and Luo, X. J. (2019). Functional genomics reveal gene regulatory mechanisms underlying schizophrenia risk. *Nat. Commun.* 10, 670. doi:10.1038/s41467-019-08666-4

Klein, K., Tremmel, R., Winter, S., Fehr, S., Battke, F., Scheurenbrand, T., et al. (2019). A new panel-based next-generation sequencing method for ADME genes reveals novel associations of common and rare variants with expression in a human liver cohort. *Front. Genet.* 10, 7. doi:10.3389/fgene.2019.00007

Knaus, B. J., and Grünwald, N. J. (2017). vcfr: a package to manipulate and visualize variant call format data in R. *Mol. Ecol. Resour.* 17 (1), 44–53. doi:10.1111/1755-0998.12549

Kovacs, J. J., Hara, M. R., Davenport, C. L., Kim, J., and Lefkowitz, R. J. (2009). Arrestin development: Emerging roles for β-arrestins in developmental signaling pathways. *Dev. Cell.* 17 (4), 443–458. doi:10.1016/j.devcel.2009.09.011

Lee, S., Wu, M. C., and Lin, X. (2012). Optimal tests for rare variant effects in sequencing association studies. *Biostat. Oxf. Engl.* 13 (4), 762–775. doi:10.1093/biostatistics/kxs014

Lee, S. Y., Hong, M. J., Jeon, H.-S., Choi, Y. Y., Choi, J. E., Kang, H. G., et al. (2015). Functional intronic ERCC1 polymorphism from regulomeDB can predict survival in lung cancer after surgery. *Oncotarget* 6 (27), 24522–24532. doi:10.18632/oncotarget.4083

Lefkowitz, R. J., and Shenoy, S. K. (2005). Transduction of receptor signals by beta-arrestins. *Science* 308 (5721), 512–517. doi:10.1126/science.1109237

Lenth, R. (2022). emmeans: Estimated marginal means, aka least-squares means. Available at: https://CRAN.R-project.org/package=emmeans.

Lesch, K. P., and Heils, A. (2000). Serotonergic gene transcriptional control regions: Targets for Antidepressant drug development? *Int. J. Neuropsychopharmacol.* 3 (1), 67–79. doi:10.1017/S1461145700001747

Li, S., Li, Y., Li, X., Liu, J., Huo, Y., Wang, J., et al. (2020). Regulatory mechanisms of major depressive disorder risk variants. *Mol. Psychiatry* 25 (9), 1926–1945. doi:10.1038/s41380-020-0715-7

Lipina, T. V., Fletcher, P. J., Lee, F. H., Wong, A. H. C., and Roder, J. C. (2013). Disrupted-In-Schizophrenia-1 Gln31Leu polymorphism results in social anhedonia associated with monoaminergic imbalance and reduction of CREB and β-arrestin-1, 2 in the nucleus accumbens in a mouse model of depression. *Neuropsychopharmacology* 38 (3), 423–436. doi:10.1038/npp.2012.197

Machiela, M. J., and Chanock, S. J. (2015). LDlink: A web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants. *Bioinformatics* 31 (21), 3555–3557. doi:10.1093/bioinformatics/btv402

Mallinckrodt, C. H., Sanger, T. M., Dubé, S., DeBrota, D. J., Molenberghs, G., Carroll, R. J., et al. (2003). Assessing and interpreting treatment effects in longitudinal clinical trials with missing data. *Biol. Psychiatry* 53 (8), 754–760. doi:10.1016/S0006-3223(02)01867-X

Matuzany-Ruban, A., Avissar, S., and Schreiber, G. (2005). Dynamics of beta-arrestin1 protein and mRNA levels elevation by antidepressants in mononuclear leukocytes of patients with depression. *J. Affect. Disord.* 88 (3), 307–312. doi:10.1016/j.jad.2005.08.007

Mendez-David, I., El-Ali, Z., Hen, R., Falissard, B., Corruble, E., Gardier, A. M., et al. (2013). A method for biomarker measurements in peripheral blood mononuclear cells isolated from anxious and depressed mice: β-Arrestin 1 protein levels in depression and treatment. *Front. Pharmacol.* 4, 124. doi:10.3389/fphar.2013.00124

Momozawa, Y., Akiyama, M., Kamatani, Y., Arakawa, S., Yasuda, M., Yoshida, S., et al. (2016). Low-frequency coding variants in CETP and CFB are associated with susceptibility of exudative age-related macular degeneration in the Japanese population. *Hum. Mol. Genet.* 25 (22), 5027–5034. doi:10.1093/hmg/ddw335

Morales, J., Welter, D., Bowler, E. H., Cerezo, M., Harris, L. W., McMahon, A. C., et al. (2018). A standardized framework for representation of ancestry data in genomics studies, with application to the NHGRI-EBI GWAS Catalog. *Genome Biol.* 19 (1), 21. doi:10.1186/s13059-018-1396-2

Muthén, B., Asparouhov, T., Hunter, A., and Leuchter, A. F. (2011). Growth modeling with non-ignorable dropout: Alternative analyses of the STAR*D antidepressant trial. *Psychol. Methods* 16 (1), 17–33. doi:10.1037/a0022634

O'Leary, N. A., Wright, M. W., Brister, J. R., Ciufo, S., Haddad, D., McVeigh, R., et al. (2016). Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation. *Nucleic Acids Res.* 44, D733–D745. doi:10.1093/nar/gkv1189

Parruti, G., Peracchia, F., Sallese, M., Ambrosini, G., MasiniM.Rotilio, D., et al. (1993). Molecular analysis of human beta-arrestin-1: Cloning, tissue distribution, and regulation of expression. Identification of two isoforms generated by alternative splicing. *J. Biol. Chem.* 268 (13), 9753–9761. doi:10.1016/S0021-9258(18)98412-7

Petit, A. C., El Asmar, K., David, D. J., Gardier, A. M., Becquemont, L., Feve, B., et al. (2018). The association of β-arrestin2 polymorphisms with response to antidepressant treatment in depressed patients. *Prog. Neuropsychopharmacol. Biol. Psychiatry* 81, 74–79. doi:10.1016/j.pnpbp.2017.10.006

Piechota, M., Golda, S., Ficek, J., Jantas, D., Przewlocki, R., and Korostynski, M. (2015). Regulation of alternative gene transcription in the striatum in response to antidepressant drugs. *Neuropharmacology* 99, 328–336. doi:10.1016/j.neuropharm. 2015.08.006

Pittenger, C., and Duman, R. S. (2008). Stress, depression, and neuroplicity: A convergence of mechanisms. *Neuropsychopharmacology* 33 (1), 88–109. doi:10.1038/snj.npp.1301574

R Core Team (2020). *R: A language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing. Available at: https://www. R-project.org/.

Rush, A. J., Kraemer, H. C., Sackeim, H. A., Fava, M., Trivedi, M. H., Frank, E., et al. (2006). Report by the ACNP Task Force on response and remission in major depressive disorder. *Neuropsychopharmacology* 31 (9), 91841–91853. doi:10.1038/ sj.npp.1301131

Sasayama, D., Hiraishi, A., Tatsumi, M., Kamijima, K., Ikeda, M., Umene- Nakano, W., et al. (2013). Possible association of CUX1 gene polymorphisms with antidepressant response in major depressive disorder. *Pharmacogenomics J.* 13 (4), 354–358. doi:10.1038/tpj.2012.18

Shukla, A. K., Xiao, K., and Lefkowitz, R. J. (2011). Emerging paradigms of β-arrestin-dependent seven transmembrane receptor signaling. *Trends biochem. Sci.* 36 (9), 457–469. doi:10.1016/j.tibs.2011.06.003

Shyn, S., Shi, J., Kraft, J., Potash, J. B., Knowles, J. A., Weissman, M. M., et al. (2011).Novel loci for major depression identified by genome-wide association study of Sequestrated Treatment Alternatives to Relieve Depression and meta-analysis of three studies/article-title. *Mol. Psychiatry* 16 (2), 202–215. doi:10.1038/mp.2009.125

Sun, D., Ma, J. Z., Payne, T. J., and Li, M. D. (2008). Beta-arrestins 1 and 2 are associated with nicotine dependence in European American smokers. *Mol. Psychiatry* 13 (4), 398–406. doi:10.1038/sj.mp.4002036

The ENCODE Project Consortium (2012). An integrated Encyclopedia of DNA elements in the human genome. *Nature* 489 (7414), 57–74. doi:10.1038/nature11247

Trivedi, M. H., Hollander, E., Nutt, D., and Blier, P. (2008). Clinical evidence and potential neurobiological underpinnings of unresolved symptoms of depression. *J. Clin. Psychiatry* 69 (2), 246–258. doi:10.4088/jcp.v69n0211

Trivedi, M. H., Rush, A. J., Wisniewski, S. R., Nierenberg, A. A., Warden, D., Ritz, L., et al. (2006). Evaluation of outcomes with citalopram for depression using measurement-based care in STAR*D: Implications for clinical practice. *Am. J. Psychiatry* 163 (1), 28–40. doi:10.1176/appi.ajp.163.1.28

Violin, J. D., and Lefkowitz, R. J. (2007). Beta-arrestin-biased ligands at seven-transmembrane receptors. *Trends Pharmacol. Sci.* 28 (8), 416–422. doi:10.1016/j. tips.2007.06.006

Ward, L. D., and Kellis, M. (2016). HaploReg v4: Systematic mining of putative causal variants, cell types, regulators and target genes for human complex traits and disease. *Nucleic Acids Res.* 44, D877–D881. doi:10.1093/nar/gkv1340

World Health Organization (2017). *Depression and other common mental disorders: Global health estimates*. Geneva: World Health Organization.

Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., and Weir, B. S. (2012). A high-performance computing toolset for relatedness and principal component analysis of SNP data. *Bioinformatics* 28 (24), 3326–3328. doi:10.1093/ bioinformatics/bts606