Generalised Anxiety Disorder – A Twin Study of Genetic Architecture, Genome-Wide Association and Differential Gene Expression

Generalised Anxiety Disorder (GAD) is a common anxiety-related diagnosis, affecting approximately 5% of the adult population. One characteristic of GAD is a high degree of anxiety sensitivity (AS), a personality trait which describes the fear of arousal-related sensations. Here we present a genome-wide association study of AS using a cohort of 730 MZ and DZ female twins. The GWAS showed a significant association for a variant within the RBFOX1 gene. A heritability analysis of the same cohort also confirmed a significant genetic component with h2 of 0.42. Additionally, a subset of the cohort (25 MZ twins discordant for AS) was studied for evidence of differential expression using RNA-seq data. Significant differential expression of two exons with the ITM2B gene within the discordant MZ subset was observed, a finding that was replicated in an independent cohort. While previous research has shown that anxiety has a high comorbidity with a variety of psychiatric and neurodegenerative disorders, our analysis suggests a novel etiology specific to AS.


Introduction
Generalised Anxiety Disorder (GAD) is a chronic disorder characterised by persistent worrying, anxiety symptoms, and tension independent of recent stressful events. GAD is one of the most prevalent anxiety disorders affecting approximately 5% of the adult population and is twice as common in women as in men [1,2]. Anxiety symptoms are frequently combined with an assortment of somatic and psychological complaints, such as sleep disruption, irritability, autonomic arousal, restlessness, fatigue and difficulties in concentrating. GAD often presents with comorbid disorders such as Major Depression Disorder, Bipolar Disorder and substance abuse [3] but it is an independent disorder which is described within the DSM-5 [4]. A strong aspect of GAD is a high level of anxiety sensitivity (AS), which is defined as the fear of arousal related sensations, arising from beliefs that these anxiety-related sensations have harmful consequences [5][6][7]. While it is unclear whether AS represents a "lower order facet of trait anxiety" or rather a distinct construct [8], current research suggests that anxiety sensitivity may usefully be conceptualised as a variable risk factor for anxiety-related problems [9]. Previous research has suggested that heightened AS may contribute to mechanisms underlying GAD, particularly the expression and processing of emotional cognitions [10,11].
A commonly used clinical metric to measure AS is the Anxiety Sensitivity Index (ASI) [12], which quantifies overall psychological and somatic symptoms related to anxiety [13][14][15]. The ASI is based upon a 16-item self-report questionnaire assessing a person's beliefs about the social and somatic consequences of anxiety symptoms. The questions are designed to assess the subject's fear of physical symptoms, fear of publicly observable anxious symptoms and fear of cognitive dyscontrol [16]. Previous studies using the ASI have repeatedly demonstrated that patients with GAD have a higher ASI score compared to healthy controls [17] and that the physical, social and cognitive facets of ASI are all indicators of clinical GAD [16,17]. The ASI has been shown to have strong internal consistency (a = 0.81-0.94), a good degree of test/retest reliability (r = 0.71-0.75), and a high degree of inter-item relatedness [13,14,18]. While MZ twins arise from the same single cell and therefore share a near identical genome, they show phenotypic discordance for many complex traits and disorders including autism (58 to 60%) [19] and schizophrenia (58%) [20]. This observation suggests that many complex traits, have a strong non-genetic component. The discordant MZ model is a powerful tool in order to identify non-genetic contributions to a phenotype including variation in gene expression not mediated by cis-or trans-eQTL effects [21]. Here we present a twin study design to identify genetic associations with AS, as well as a differential gene expression analysis of MZ twins who are discordant for AS.

Sample Collection
AS data was collected from participants previously recruited from the TwinsUK registry [22,23]. The 730 participants comprised 143 MZ and 222 DZ twin pairs with an age range of 38-84 years. The TwinsUK cohort has previously been shown to be comparable to population singletons in terms of disease-related and a range of lifestyle characteristics [24,25]. The participants completed a detailed health questionnaire which included the ASI questionnaire (see S1 Table). Data on potential confounding factors such as age, BMI, smoking and alcohol consumption were also collected as part of the study. Peripheral blood samples were collected, and LCLs were generated by EBV transformation of the B-lymphocyte component by the European Collection of Cell Cultures agency [26]. The St. Thomas' Research Ethics Committee (REC) approved on 20 September 2007 the protocol for the dissemination of data, including DNA, with REC reference number RE04/015. On 12 March 2008, the St Thomas' REC confirmed that this approval extended to expression data. Volunteers gave informed consent and signed an approved consent form before the biopsy procedure. Volunteers were supplied with an appropriate detailed information sheet regarding the research project.

Heritability calculation
A heritability analysis of the cohort of 730 MZ and DZ female twins was conducted with the structural equation modeling package Mx [27]. The goodness of fit of the genetic models was evaluated by comparing them with an unconstrained saturated model, which estimates the maximum number of parameters without partitioning variance into genetic and environmental components. This gives rise to a likelihood ratio chi-square tests(-2LL; significant tests indicate significant deterioration in fit) and Akaike's Information Criterion (AIC; a parsimony fit index, with lower values indicating the more suitable model) [28]. The best-fitting models were selected on the basis of parsimony.

Genotyping and GWAS
Genotyping of the cohort of 730 MZ and DZ female twins was carried out with a combination of Illumina arrays (HumanHap300, HumanHap610Q, 1M-Duo, and 1.2MDuo 1M). Intensity data for each of the three arrays were pooled separately (with 1M-Duo and 1.2MDuo 1M pooled together), and genotypes were called with the Illuminus calling algorithm with the use of a threshold on a maximum posterior probability of 0.95. Samples were imputed into the 1000 Genomes Phase 1 reference panel (data freeze, 10/11/2010) [29] using IMPUTE2 [30] and filtered (MAF<0.05, IMPUTE info value < 0.8) [26]. A GWAS analysis on the full cohort of twins using AS as the phenotype was implemented using the program GEMMA [31]. The analysis used age as a covariate and also incorporated the family structure of the twins (see Fig 1A and 1B, Table 1).

RNA processing
Samples were prepared for sequencing with the IlluminaTruSeq sample preparation kit (Illumina, San Diego, CA) according to the manufacturer's instructions and were sequenced on a HiSeq2000 machine. An average of 28 million exonic reads per sample was generated. Afterwards, the 49-bp sequenced paired-end reads were mapped to the GRCh37 reference genome [32] with BWA v0.5.9 [33] allowing 2 mismatches in the seed (30 first bases of the read).
Gene expression quantification and the MZ discordant model. Genes defined as protein coding in the GENCODE 10 annotation [34] were selected for analysis. RPKM values calculated from the RNA-seq data were root mean transformed and PEER software [35] was used to remove 50 latent factors, after which the data were quantile normalised. Normalised expression data from the LCL tissues were examined in a subset of the samples containing MZ twins discordant for AS with one twin with an AS <10 (control) with a sibling with an ASE score >15 points higher (case). 25 discordant MZ pairs matching the criteria could be identified. A linear mixed effect model was fitted for the scaled exon values of LCL tissue using the R package lmer [36]. The model incorporated ASI, age, BMI, smoking and alcohol consumption of the twins as fixed effects predictors and family structure as a random effect.

Validation Data
Validation for the GWAS analysis was taken from a current meta-analysis of anxiety disorders. This study, conducted in over 18,000 subjects of European descent, uses two phenotypic strategies to test association with common variants genome-wide that index genetic susceptibility shared across the primary anxiety disorders. The first compared cases with unaffected controls; the second used factor analysis to estimate a quantitative factor score that indexed a common latent risk phenotype [37].
Expression analysis was replicated in an independent dataset of U219 Affy gene expression arrays for 586 twins from the Netherlands Study of Depression and Anxiety (NESDA) and Netherlands Twin Register (NTR) cohorts [38,39] The phenotype used for the NESDA Samples is the Beck Anxiety Inventory (BAI) [40], which is a 27-item questionnaire ranging from 0 to 63 also with high internal consistency (Cronbach's α = 0.92) [41]. A score of 0-9 refers to normal severity, whereas a score of 10-18 refers to mild severity, a score of 18-29 refers to moderate severity, and a score higher than 29 refers to severe anxiety symptoms. 22 MZ twins highly discordant for BAI were identified for the replication analysis.

Ethics Statement
The study was approved by the St Thomas' Hospital Research Ethics Committee (REC reference: EC04/015) in accordance with the St Thomas' Hospital Local Ethics Committee. All participants in the study provided written informed consent.

Heritability of ASI
Heritability analysis of the ASI score using Mx estimated h 2 as being 0.4445 (CI: 0.3162-0.5548). This is comparable with previous estimate of the heritability of anxiety related traits (see Discussion).

GWAS study
In the AS GWAS analysis, a single variant of genomewide significance (assuming genomewide significance to have an unadjusted p-value < 5x10 -8 [42]) occurred within the coding region of the RBFOX1 gene (rs13334105,p value: 4.39 x10 -8 , MAF = 5.164, T/G). Additionally, several SNPs of suggestive significance also occured within the coding region of the RBFOX1 gene(see Fig 2A and  2B). The significant association in RBFOX1 was compared with the anxiety meta-analysis data (see Data Validation). Although the specific variant rs13334105 was not included in the analysis, both the case control and factor scores approach each identified a SNP within the RBFOX1 intronic region with unadjusted p-values <0.05. The variants were rs17664315 (p = 0.00037) and rs7203983 (p = 0.00250) for case control and factor scores respectively (see Fig 3A and 3B).

Discordant MZ twin model
An alternative approach to the study of complex traits is the use the twin MZ discordant model [43,44] whereby differences in gene expression can be identified by simultaneously controlling for genotypic variants. Here, we implemented a linear mixed effect model in order to identify differential expression amongst twins discordant for ASI. After correcting for the number of exons interrogated (116,528), an adjusted p-value < 0.05 was observed for two out of six exons (Exon 5 and Exon 6) within the gene ITM2B with a further two exons (Exon 3 and Exon 4) approaching significance. The exons of greatest significance within the model are shown in Table 2. A plot of the raw case control expression data is plotted in Fig 4. 22 MZ twin pairs discordant for BAI were selected from the U219 Affy gene expression arrays (see Validation Data) In the U219 array, multiple 25 bp probes targeting exon 4 (chr13:48832933-48833083), exon 5 (chr13: 48835273-48836232) and exon 6 (chr13: 48832262-48832372)of ITM2B (see Fig 5) were used to assess gene expression levels. Expression values were residualised by age, sex, BMI, cell counts and technical covariates and the difference between residualized gene expression were computed between pairs. T-tests for this difference on each of the three exons showed p-values of 0.03, 0.3 and 0.002 for exons 4, 5 and 6 respectively. As was the cases with the RNA-seq data, twins with higher anxiety scores in the replication set showed higher expression levels of ITM2B (see Fig 6). The fixed effects incorporated showed no significant association with AS.

Discussion
Although blood-derived samples may have limited application in the study of neurological function due to the limitations of tissue specificity, previous studies have suggested they can be used as a surrogate tissue for the evaluation of genetic factors in psychiatric disorders [45,46]. Our heritability analysis of AS (h 2 = 0.445) confirmed previous work suggesting an important genetic contribution to anxiety-related disorders. A study calculating the heritability of AS as a unifactorial construct produced a comparable estimate of h 2 = 0.45 [47], while in a separate study the heritability of AS was estimated as being 0.48 [48]. The result is also comparable to a previous measure of anxiety using the Anxiety-Related Behaviours Questionnaire (ARBQ), which assesses four anxiety traits [49] and recorded a heritability of 0.52.
Our GWAS analysis also suggested a genetic basis to ASI through the identification of a SNP of genomewide significance and a further eight SNPs of near significance occurring within the coding region of the binding protein, fox-1 homolog (RBFOX1) gene (see Fig 2A and 2B, Table 1). A further GWAS incorporating rare variants (MAF 1-5%) identified another variant of genomewide significance with the RBFOX1 gene. RBFOX1 (also known as ataxin 2-binding protein 1 (A2BP1) or hexaribonucleotide-binding protein 1 (HRNBP1) codes for the Fox1 protein and regulates alterative splicing which controls the gene expression coordinating neuronal brain activity [50,51]. Abnormalities within the limbic system are the main factor in eliciting anxiety-like symptoms [52,53]. The limbic system is characterised by a significant density of monoaminergic and GABAergic neurons. Patients with GAD have abnormal GABAergic activity resulting from the down-regulation of Gammaamino butyric acid (GABA-A) receptor, an ionotropic receptor and ligand-gated ion channel [54]. Positive modulators of GABA receptors generally possess anxiolytic (anti-anxiety) activity Genome-Wide Association of Generalised Anxiety Disorder while negative modulators produce anxiogenic-like effects. The action of anxiolytic drugs work by suppressing this overactive transmission by increasing the GABA-A signal [52]. A probable mechanism by which RBFOX1 influences GABAergic neuronal function is shown in Fig 7 [55]. Immunostaining techniques have shown that there is an increased expression of the active  Fox1 isoform in the nucleus following neuronal depolarisation [56]. Fox1 may regulate the action of Gabrg2, a gene which encodes the gamma 2 subunit of GABA(A) receptor, by alternative splicing [57]. Gabrg2 regulation is controlled by the binding of Fox1 near the splice sites of the downstream intron using (U)GCAUG targeting. This sequence inserts 8 amino acids on the intracellular loop between the M3 and M4 transmembrane domain site of the GABA receptor [58][59][60], creating a putative serine phosphorylation site for protein kinase C. [61,62] Phosphorylation of serine has shown in to alter the tertiary structure of the γ 2 subunit of the GABA (A) receptor, which has control of the GABA-A receptor function [61][62][63][64].
In clinical settings, variations within the RBFOX1 gene have been shown to be associated with autism spectrum disorder (ASD), schizophrenia, epilepsy and attention deficit hyperactivity disorder (ADHD) [65]. This pathology has been clinically evaluated in a case study  involving 16-year-old male who has a heterozygous deletion spanning 100.3 kb in the 5' region of the RBFOX1. The patient showed speech behavioural abnormalities including anxiety like behaviour towards certain stimuli. This supports a previous study on patients that had copy number variation in the region of RBFOX1, which was shown to be associated with a variety of neurodevelopmental and neuropsychiatric disorders [65].

Differential gene expression between MZ twins discordant for ASI
The MZ twin discordance analysis showed statistically higher significant expression of exons 5 and 6 in the ITM2B (BRI2) gene for the anxiety cohort (see Table 3). ITM2B gene encodes a transmembrane protein which is processed at the C-terminus by furin or furin-like proteases to produce a small secreted peptide which inhibits the deposition of beta-amyloid [66]. A point mutation at the stop coding of ITM2B causes a build up of beta-amyloid plaques which can led to neuronal cell death. ITM2B is therefore strongly associated with the neurodegenerative diseases, such as British and Danish Familial Dementia [67]. The conditions are characterised by significant neuronal loss, a loss of cognitive and memory functions, disruptions to muscle movement and the induction of special rigidity [66,67].
The accumulation of beta-amyloids within the amygdala has been shown to drive fear and anxiety-like behavior in transgenic mice models [68]. These observed behaviours are distinct from anxiety symptoms associated with ageing [69]. This model suggests a mechanism for the neuronal connectivity loss that is characteristic of amygdala in GAD [70] and is consistent with studies that demonstrate anxiety-like symptoms being present in neurodegenerative disorders [71]. That said, increased anxiety may also be due to memory loss and physical impairment that are symptomatic of neurodegeneration.

Conclusion
Our analyses suggest an involvement of the RBFOX1 gene in the development of anxiety-related conditions such as GAD. The GWAS results were supported by an independent study which also identified variants in the RBFOX1 gene which were linked to AS. We further identified two exons of ITM2B as having significantly higher expression levels in individuals with higher levels of AS than their MZ twin, an effect we saw replicated in an independent sample. This also suggests that both genetic and non-genetic variation is crucial to the development of the disorder. Our results may indicate novel etiological mechanisms and paths in the disease development or maintenance of the disease. Further investigation would be necessary to determine whether these variants are specific to AS or whether they are common to more general anxiety disorders.