Transcriptome-Wide Association Study Provides Insights Into the Genetic Component of Gene Expression in Anxiety

Anxiety disorders are common mental disorders that often result in disability. Recently, large-scale genome-wide association studies (GWASs) have identified several novel risk variants and loci for anxiety disorders (or anxiety traits). Nevertheless, how the reported risk variants confer risk of anxiety remains unknown. To identify genes whose cis-regulated expression levels are associated with risk of anxiety traits, we conducted a transcriptome-wide association study (TWAS) by integrating genome-wide associations from a large-scale GWAS (N = 175,163) (which evaluated anxiety traits based on Generalized Anxiety Disorder 2-item scale (GAD-2) score) and brain expression quantitative trait loci (eQTL) data (from the PsychENCODE and GTEx). We identified 19 and 17 transcriptome-wide significant (TWS) genes in the PsychENCODE and GTEx, respectively. Intriguingly, 10 genes showed significant associations with anxiety in both datasets, strongly suggesting that genetic risk variants may confer risk of anxiety traits by regulating the expression of these genes. Top TWS genes included RNF123, KANSL1-AS1, GLYCTK, CRHR1, DND1P1, MAPT and ARHGAP27. Of note, 25 TWS genes were not implicated in the original GWAS. Our TWAS identified 26 risk genes whose cis-regulated expression were significantly associated with anxiety, providing important insights into the genetic component of gene expression in anxiety disorders/traits and new clues for future drug development.


INTRODUCTION
Anxiety disorders are common mental disorders (Huang et al., 2019) and a leading cause of disability (GBD 2017 Disease andInjury Incidence andPrevalence Collaborators, 2018). As an adaptive behavioral response, anxiety is usually a normal emotional experience in daily life. However, individuals affected by anxiety disorders frequently have feelings of excessive, intense, and persistent worry and fear about the actual or anticipated event. If the degree of anxiety exceeds a certain range and the abnormal feeling of fear and worry lasts for a long time, normal anxiety experiences may turn into anxiety disorders. Anxiety disorders are highly prevalent mental disorders (with a lifetime prevalence rate of over 20% (Kessler et al., 2012)) that bring huge economic burdens to the world (Baxter et al., 2014). The etiology of anxiety disorders remains poorly understood. However, it is found that genetic factors play pivotal roles. The heritability of anxiety disorders were estimated to range from 20 to 60% (Polderman et al., 2015;Craske et al., 2017), indicating substantial genetic components of anxiety disorders. To uncover the genetic basis of anxiety disorders/traits, several genetic studies have been performed (Sipilä et al., 2010;Stein et al., 2011;Nivard et al., 2014). Nevertheless, large-scale GWASs of anxiety disorders/traits have not emerged until recently (Otowa et al., 2016;Meier et al., 2019;Levey et al., 2020;Purves et al., 2020). For example, by using the UK biobank sample, Purves et al. identified five genome-wide significant (GWS) risk loci for anxiety disorder (Purves et al., 2020). In addition, Meier et al. reported that PDE4B is associated with anxiety and stress-related diagnoses in the Danish population (Meier et al., 2019). Otowa et al. also found that a noncoding RNA region located in 3q12.3 is associated with anxiety disorder (Otowa et al., 2016). Notably, a recent large-scale GWAS (N 199,611) conducted by Levey et al. (i.e., the Million Veteran Program, MVP) reported six GWS loci for Generalized Anxiety Disorder 2-item scale (GAD-2) GWAS (Levey et al., 2020). There are also other GWASs for specific anxiety disorders/ traits, such as Generalised Anxiety Disorder (Davies et al., 2015;Dunn et al., 2017) , panic disorder (Otowa et al., 2009;Erhardt et al., 2011;Otowa et al., 2012), phobic anxiety (Walter et al., 2013), and social anxiety . Although these GWASs have identified several risk variants and loci for anxiety disorders/traits, the number of reported anxiety disorders/traits risk loci remains small compared with other psychiatric disorders (including schizophrenia (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014; Pardiñas et al., 2018;Lam et al., 2019), depression (Wray et al., 2018), and bipolar disorder (Stahl et al., 2019)).
To date, most of the risk variants reported by anxiety disorders/traits GWASs are located in noncoding regions (Otowa et al., 2016;Meier et al., 2019;Levey et al., 2020;Purves et al., 2020), suggesting that these genetic variations confer risk of anxiety disorders/traits by modulating gene expression. However, we currently know little about the target genes regulated by these reported risk variants. In addition, the joint effects of the reported anxiety disorders risk loci on gene expression across human brain tissues remain elusive. The emergence of TWAS (Gusev et al., 2016) provides an opportunity to infer genes whose cis-regulated expression may be associated with diseases and traits. TWAS first leverages genotype and expression data from an external reference panel to identify the associations between genetic variations and gene expression. Then it identifies genes whose cis-regulated expression are associated with diseases or phenotypes by integrating gene expression data (i.e., external reference panel) with genome-wide associations from large-scale GWASs. TWASs have been widely used to identify risk genes for psychiatric disorders and neurodegenerative disease, including schizophrenia (Gandal et al., 2018), depression (Dall'Aglio et al., 2020, Parkinson's disease , and attention deficit hyperactivity disorder (Liao et al., 2019). Moreover, TWASs have identified multiple novel risk genes that were not implicated by the original GWASs (Liao et al., 2019). In this study, we performed a TWAS (using FUSION software (Gusev et al., 2016)) on anxiety traits by integrating genome-wide associations from a recent large-scale GAD-2 score GWAS of anxiety disorder (Levey et al., 2020) (i.e., the Million Veteran Program) and gene expression data from multiple brain tissues. We identified 26 TWS genes, of which 19 were from the PsychENCODE gene expression reference and 17 were from GTEx. It is noteworthy that 10 genes (KANSL1-AS1, CRHR1, CRHR1-IT1, SPPL2C, RP11-707O23.5, RP11-259G18.1, MAPT-AS1, LRRC37A4P, PLEKHM1, and DND1P1) showed TWS in both datasets (i.e., PsychENCODE and GTEx), implying these genes are promising candidate genes for anxiety. Moreover, we noticed that some novel risk genes identified by our TWAS were not implicated in the original GWAS. Our study identifies genes whose expression change may confer risk of anxiety disorders/ traits, shading lights on the genetic component of gene expression in anxiety disorders/traits, and providing a start point for mechanistic investigation and future drug development.

Genome-Wide Association Studies Summary Statistics
The GAD-2 anxiety GWAS summary statistics were from a recent large-scale GWAS of anxiety disorder (Levey et al., 2020) and were downloaded from dbGAP (https://dbgap.ncbi. nlm.nih.gov/, phs001672) (Levey et al., 2020). This GWAS was performed by the Million Veteran Program (MVP) which is a large biobank that includes information about genetic, environmental, and medical records (Gaziano et al., 2016). GAD-2 (Kroenke et al., 2010) was used to measure the degree of anxiety. European Americans (N 175,163) and African Americans (N 24,448) were included in the MVP. Considering that gene expression weights data were based on European populations, only GWAS summary statistics from European Americans were used in this study. Samples were genotyped with customized MVP Affymetrix genotyping array. Details about sample information, phenotyping, data processing, and statistical analyses were provided in the original paper (Levey et al., 2020).

Gene Expression References Panels
Two sets of gene expression weights (SNP-gene expression associations) were used to perform the TWAS analysis. The first gene expression reference panel was from the PsychENCODE (Gandal et al., 2018). Briefly, gene expression and genotypes of 1,321 individuals were downloaded from the PsychENCODE website (http://resource.psychencode.org) and gene expression weights (SNP-gene expression relations) were generated by using the R script (FUSION.compute_weights.R) provided by the FUSION software. Genotype data from the PsychENCODE were used as the LD (linkage disequilibrium) reference panel (a file (PLINK format) (Purcell et al., 2007) which contains the genotype of a reference population that matches the GWAS sample ancestry. TWAS can use a personalized reference panel or the one provided by FUSION). The second expression panel was from GTEx (Battle et al., 2017) (V7) and only expression weights from brain tissues (a total of 13) were used. For TWAS using the GTEx dataset, the LD reference panel (from European ancestry) was downloaded from the FUSION website (https://data.broadinstitute.org/alkesgroup/ FUSION/LDREF.tar.bz2) provided by Gusev et al. (2016) More detailed information about FUSION software, PsychENCODE, and GTEx were provided in the original papers (Gusev et al., 2016).

Transcriptome-Wide Association Study
To identify genes whose genetically-regulated expression are associated with risk of anxiety, we performed TWAS analyses by integrating GWAS summary statistics of anxiety and gene expression weights from the PsychENCODE and GTEx by using FUSION software (http://gusevlab.org/projects/fusion/) (Gusev et al., 2016) (with the using of default settings). The basic process were as follows: Firstly, FUSION computed TWAS expression weights (i.e., SNP-gene expression correlations) using five linear models, including BLUP, BSLMM, LASSO, Elastic Net, and top SNPs from the reference expression panels (i.e., PsychENCODE and GTEx). When performing transcriptomic imputation, FUSION calculated an out-sample R 2 using a fivefold crossvalidation of each model to determine the best performing prediction model for a gene. Then the imputed gene expression was used to test the association with anxiety. Data from the PsychENCODE (http://resource.psychencode.org/ #Derived) (Gandal et al., 2018) and Gusev et al. (https://data. broadinstitute.org/alkesgroup/FUSION/LDREF.tar.bz2) (Gusev et al., 2016) were used as LD reference panels to account for LD structure. Transcriptome-wide significant genes were corrected by the Bonferroni correction approach (p 3.52 × 10 -06 , 0.05/14,223 (the number of genes in the weight files) for PsychENCODE. Variable P thresholds were applied to GTEx as the gene number included in different brain weight files were different. GTEx P thresholds 0.05/(the number of genes in each GTEx weight file)/13 (13 was the number of GTEx reference panels included in this study).

Joint/Conditional Analysis
We performed Joint/Conditional analysis on transcriptomewide significant (TWS) loci by utilizing the script FUSION.post_process.R implemented in FUSION (Gusev et al., 2016). The Joint/Conditional analysis aimed at exploring to what extent the GWAS signals remain significant after removing TWAS significant signals (i.e., test if the GWAS signals are still significant after removing the expression weights from TWAS). This analysis requires the top significant genes from TWAS analysis, LD reference panel, and GWAS summary statistics as input. Each SNP association from anxiety GWAS was conditioned on the joint model (one SNP at a time) and permutation test (100,000) was used to determine the p-value of Joint/Condition analysis results. The genes that passed the permutation test represent promising driven genes (i.e., these genes are unlikely to be co-localized due to chance). Our TWAS analysis pipeline is in agreement with the original FUSION paper. Default settings and parameters (recommended by FUSION software) were used in Joint/ Conditional analysis.

Tissue and Cell-Type Enrichment Analysis
To explore if the genetic associations identified by anxiety GAD-2 GWAS were enriched in specific tissues and cell types, we utilized MAGMA  (implemented in FUMA (Watanabe et al., 2017)) (https://fuma.ctglab.nl/) to conduct tissue and cell-type-specific enrichment analysis (with the using of default parameters and settings). FUMA detected tissue-specific enrichment of genome-wide associations using expression data from 53 GTEx tissues (Battle et al., 2017). Briefly, FUMA first defined the differentially expressed genes (DEGs) for each tissue (by comparing the expression level of a gene in a specific tissue to all other tissues). The defined DEGs for each tissue contained genes with the greatest expression discrepancy in the specific tissue compared with other tissues. These DEGs were then used as input in enrichment analysis with MAGMA. Details, principles, and procedures have been described in the original paper (Watanabe et al., 2017) and the FUMA website.
Cell-type enrichment analysis was also performed by MAGMA . The single-cell gene expression data of the mouse central nervous system was obtained from Zeisel et al. (2018). Data was processed following the instructions by Bryois et al. (2020). Gene expression data of 160,769 single cells were analyzed and genes that were not expressed were removed. We calculated gene expression specificity as described by Bryois et al. (2020), and the top 10% genes ranked by gene expression specificity of each cell was remained for MAGMA gene set analysis. Mouse gene ids were mapped to their corresponding human genes (one-to-one orthologous) based on MGI annotations (http://www.informatics.jax.org/homology.shtml). FDR (false discovery rate) was used for multiple correction.

Pathway and Gene Ontology Enrichment Analysis
We utilized gene network v2.0 (https://genenetwork.nl/) (Deelen et al., 2019) to analyze the co-expression pattern and to investigate if TWS genes are enriched in specific pathways/GO terms. GeneNetwork contains a gene expression matrix from public RNA-seq datasets with 31,499 samples and 56,435 genes. Firstly, a PCA (principal component analysis) was performed on this gene expression matrix, and 1,588 principal components were selected. Secondly, for each PC, a t-test was performed for genes that belonged to a GO term/pathway and other genes. Thirdly, a Mann-Whitney U test was performed to determine the Z-score of the annotated genes and genes that were not in the network. Gene sets including REACTOME, GO, and KEGG pathways were used. For details on RNA-seq data processing, PCA, co-regulation analysis, and other statistical analysis, please refer to the original paper of network v2.0 39

Transcriptome-Wide Association Study Identified 26 Risk Genes for Anxiety
We conducted TWAS by integrating GAD2-GWAS summary statistics of anxiety and two sets of brain gene expression reference panels (i.e., eQTL data) from PsychENCODE and GTEx. We identified 19 TWS genes (from four genomic loci, 2p23.3, 3p21.31, 3p21.1 and 17q21.31) when PsychENCODE brain eQTL data were used ( Figure 1A; Table 1;  Supplementary Table S1). Of note, genes from these four genomic loci did not show genome-wide associations with anxiety in the original GWAS (Levey et al., 2020).
Though the sample size of the PsychENCODE is large, most of the brain tissues used for eQTL analysis were from the dorsolateral prefrontal cortex. We thus further carried out TWAS using brain SNP-gene expression weights (SNP-gene expression relations) from 13 GTEx brain tissues (Battle et al., 2017). A total of 17 TWS genes spanning three genomic loci (3p21.31, 17q21.31, 20q13.33) were identified in this analysis ( Figure 1B, Table 1, Supplementary Figure S1-S12, Supplementary Table S2). Intriguingly, genes from two genomic loci did not show genome-wide significant associations with anxiety in the original GWAS. We compared the results from these two analyses and identified 10 overlapping genes that showed transcriptome-wide significant associations with anxiety in both datasets ( Table 1). These overlapping genes include KANSL1-AS1, CRHR1, CRHR1-IT1, SPPL2C, RP11-707O23.5, RP11-259G18.1, MAPT-AS1, LRRC37A4P, PLEKHM1, and DND1P1. These findings identified genes whose genetically regulated expression may confer risk of anxiety, suggesting that genetic variants may confer anxiety risk by regulating the expression of these genes.

Conditional Analysis of Transcriptome-Wide Association Study Significant Loci
As several TWS genes were identified in each risk locus, we next sought to explore which gene drove the TWAS signal (i.e., to test if the transcriptome-wide association signal is conditionally independent or not). We therefore performed Joint/ Conditional tests for TWAS significant regions on 3p21.31, and 17q21.31, two regions supported by TWAS of two different brain eQTL datasets (PsychENCODE and GTEx). We identified several independent transcriptome-wide significant genes from both of the brain eQTL datasets. For example, we found that RNF123 explains most of the variance (0.725) at its locus in PsychENCODE dataset (rs34484573 lead SNP P GWAS 2.7 × 10 -06 , conditioned on RNF123 lead SNP P GWAS 0.0091) ( Figure 2A). In addition, two genes on 17q21.31, including MAPT−AS1and KANSL1−AS1 ( Figure 2B), are also jointly significant in PsychENCODE dataset. These two genes jointly explain most (0.894) a large proportion of the variance at this locus (rs2696689 lead SNP P GWAS 2.5 × 10 -07 , conditioned on MAPT-AS1 and KANSL1−AS1 lead SNP P GWAS 0.60). In addition, the 17q21.31 region also showed TWS in GTEx frontal cortex (BA9) region ( Figure 1B), the TWS gene DND1P1 explained the most of the variance (0.951) in this loci (rs17631676 lead SNP P GWAS 7.2 × 10 -07 , conditioned on DND1P1 lead SNP P GWAS 0.27) (Supplementary Figure S13).

Anxiety Associations Were Enriched in Brain Tissues and Neurons
To explore if the genome-wide associations of anxiety were enriched in specific tissues, we performed tissue-specific enrichment analysis using MAGMA  (implemented in FUMA (Watanabe et al., 2017) (Figure 3, Supplementary Table S3). The MVP GAD-2 score anxiety GWAS summary statistics were used as input for MAGMA. Overall, we found significant enrichment of anxiety associations in brain tissues (Figure 3, Supplementary   Table S3), including the cortex (p 0.0026), frontal cortex (p 0.0031), anterior cingulate cortex (p 0.013), etc. Interestingly, we noticed that anxiety associations were also enriched in the testis (p 0.0021). We then investigated the enrichments of anxiety associations in different brain cell types and found significant enrichments in neurons from the hindbrain (p 0.00018), inhibitory (p 0.0012), and excitatory (p 0.014) neurons, as well as cholinergic and monoaminergic neurons (p 0.026). These data prioritized the possible tissues and cell types that were affected by anxiety risk genes (Figure 3, Supplementary Table S4).

Potential Biological Relevance of Transcriptome-Wide Association Study Significant Genes
To explore the potential biological implications of TWAS findings, we utilized the gene network (V2) (Deelen et al., 2019) (https://genenetwork.nl/gene-list/) to carry out gene network analysis. The genes were clustered based on the coexpression relationship based on public RNA-seq data (n 31,499). Our gene network analysis showed that the 10 TWS genes (in both PsychENCODE and GTEx datasets) ( Table 1) formed two clusters (or networks) based on their co-expression  Figure S14). The first cluster includes CRHR1, LRRC37A4P, KANSL1-AS1, SPPL2C, DND1P1, MAPT-AS1. And the other cluster contains CRHR1-IT1, RP11-259G18.1, PLEKHM1, and RP11-707O23.5. These results suggest that anxiety risk genes tend to form network to exert their biological effects.
We next performed pathway and GO Enrichment analysis (based on REACTOME, GO, KEGG) to test if the TWS genes were enriched in specific biological pathways ( Table 2). REACTOME enrichment analysis showed that the TWS genes are enriched in myogenesis processes, including CDO in myogenesis (p 2.8 × 10 -4 ) and myogenesis (p 2.8 × 10 -4 ). Of note, GO analysis showed that the TWS genes were significantly enriched in the GO biological process Wnt signaling pathway, calcium modulating pathway (p 1.9 × 10 -03 ), indicating that risk variants may confer risk of anxiety by affecting these biological processes. is all genes that located in the loci (usually gray), the genes with marginally TWAS association were marked in blue, genes that are jointly significant are in green. The bottom panel is the Manhattan plot of the original GWAS summary statistics data before (gray) and after (blue) conditioning on the green genes. The arrows indicate the association result of lead SNPs before (black arrow) and after (red arrow) conditional/joint analysis.

DISCUSSION
In this study, we reported the first TWAS of anxiety to identify genes whose genetically regulated expression may confer risk of anxiety. Our study identified 26 TWS genes, of which 25 were not nominated in the original GWAS, indicating that TWAS is a powerful approach to nominate novel risk genes that are not implicated by GWAS. And 10 TWS genes were identified in both PsychENCODE and GTEx datasets. Among the 10 overlapping genes, we found that LRRC37A4P is widely expressed in several neuronal cell types. None of these genes showed cell-type specific expression (Polioudakis et al., 2019) Figure  S15-S20). Furthermore, we also found that genes located in 3p21.31 and 17q21.31 are significantly associated with anxiety. We noticed that these two loci also showed significant associations with other psychiatric disorders. For example, FIGURE 3 | Tissue and cell-type enrichment analysis of anxiety disorder summary statistics (MAGMA). Tissues and cell types that showed significant enrichment (p < 0.05) were marked in red.

(Supplementary
Frontiers in Genetics | www.frontiersin.org September 2021 | Volume 12 | Article 740134 7 MAPT gene (located in 17q21.31) is a well-known candidate risk gene for Parkinson's disease (Zabetian et al., 2007;Simón-Sánchez et al., 2009;Edwards et al., 2010). In addition, the 3p21.31 and 17q21.31 are also significantly associated with Post-traumatic stress disorder (PTSD) in a recent GWAS study by the MVP (Gelernter et al., 2019). Intriguingly, a recent TWAS on PTSD also indicated that genes located in these two regions are related to PTSD (Girgenti et al., 2021). Notably, RNF123 was one of the top genes that were significantly associated with PTSD at the transcriptome-wide significance level (Girgenti et al., 2021). The shared TWAS findings may be due to the high genetic correlation (approximately 50-60%) between PTSD and anxiety . In addition, RNF123 was also reported to be dysregulated in depression and may serve as a clinical biomarker for depression (Glahn et al., 2012;Teyssier et al., 2013). These results collectively indicate that RNF123 may play an important role in stress-related disorders, including anxiety, PTSD, and depression. Though anxiety disorder, PTSD, and depression are diagnosed as different mental disorders (as they represent potentially fundamentally different dimensional endophenotypes (Insel and Cuthbert, 2009)), the genetic risk factors that these diseases shared may represent general genetic risk factors of these mental diseases (Smoller, 2016). Taken together, these lines of evidence indicate that these two genomic regions may harbor authentic risk genes for anxiety and genetic variants likely confer risk of anxiety by modulating the expression of genes located in 3p21.31 and 17q21.31.
RNF123 (Ring Finger Protein 123) encodes E3 ubiquitinprotein ligase which is a subunit of the Kip1 ubiquitinationpromoting complex (KPC). Previous studies revealed the important function of RNF123, including regulation of cell cycle process (Kamura et al., 2004) and innate antiviral signaling (Wang et al., 2016). Besides, several studies also showed that RNF123 may have a role in cancer, including aggressive glioblastoma and melanoma (Iida et al., 2017;Wang et al., 2020). However, the function of RNF123 in the human central nervous system and the corresponding mechanisms about how it confers the risk of mental disorders remain to be investigated.
Two major reasons might account for the observations that some of the loci were significant in original GWAS but not in TWAS. Firstly, TWAS was designed to detect the gene-trait association by performing transcriptome imputation of GWAS data with a reference trained gene expression predictive model (constructed from a dataset with both genotype and gene expression) (Gusev et al., 2016;Wainberg et al., 2019). If the genetic variations in the GWS loci are not functionally as cisregulation variations (which affect gene expression), TWAS is not able to nominate candidate genes in these loci. Secondly, the power of the TWAS may be limited by the sample size of the training set we used, e.g., the psychENCODE dataset we used in this study can pinpoint more TWS genes than a single GTEx panel. We also identified loci that were not significant in original GWAS but significant in TWAS. In TWAS analysis, we can identify TWS loci that are not nominated by GWS, which is the power and advantages of TWAS (Gusev et al., 2016;Liao et al., 2019). For TWS loci not reported by GWS, a possible reason is that the statistic power is not enough in the current GWAS. For example, the leading GWAS association signal in TWS locus (17q21.31) is rs142925250 (5.95E-08, Table 1), which is close to GWS level (5.00E-08). With the increase of sample size, rs142925250 may reach the genome-wide significance level. Finally, we also identified loci that reached GWS and TWS simultaneously. The GWAS result is not easy to interpret because of complex LD (linkage disequilibrium) (several linked genes may show associations with a GWAS trait). TWAS leverages a reference expression panel (samples with gene expression and genotype data simultaneously) to perform transcriptome imputation on a target GWAS dataset, and the predicted gene expression of the GWAS dataset was used to identify the gene-trait relationship. So TWAS can help to nominate gene-trait relationships by integrating eQTL data and GWAS data, which may help to explain the functions of the GWS loci (Wainberg et al., 2019).
As the predicted expression change of the TWS genes may contribute to disease risk, we also explored the expression of TWS genes in individuals with anxiety disorder and controls, using the published differential gene expression study by Wingo et al. (Wingo and Gibson, 2015) Gene expression in the blood of 157 anxiety disorder cases and 179 controls (GSE61672) were measured by Wingo et al (Wingo and Gibson, 2015). We performed differential expression analysis using NCBI GEO2R  (Barrett et al., 2013), and found no overlap between the differentially expressed genes (DEGs) and our TWAS findings. We hypothesized that this might be due to the tissue difference between DEG analysis (blood) and TWAS (brain).
We performed MAGMA analysis for two purposes. First, the original GAD-2 GWAS papers did not perform MAGMA geneset enrichment analysis to interpret the GWAS results (Levey et al., 2020). We therefore performed gene-set enrichment analysis by MAGMA (de Leeuw et al., 2015) to explore if genome-wide associations were enriched in specific tissues or cell types (using expression data from GTEx and single-cell RNA sequencing dataset). Second, we aimed at finding out the mostrelevant tissues for TWAS analysis when using GTEx reference panels. Our MAGMA tissue enrichment analysis indicated that the brain tissues in GTEx were suitable for performing TWAS. As expected, we found that the anxiety associations were enriched in brain tissues and neurons. Of note, significant enrichment in the testis was also observed. A possible reason for this observation is that only limited genome-wide associations were reported by the original GWAS. Though the sample size of the MVP is relatively large (N 175,163), we noticed that only five genome-wide significant (GWS) loci have been identified in the original GWAS of anxiety. The number of GWS risk loci for anxiety is quite small compared with other psychiatric disorders such as schizophrenia and depression. In addition to sample size, we found that the heritability were different among anxiety disorders (20-60%) (Polderman et al., 2015;Craske et al., 2017), SCZ (approximately 80%) (Sullivan et al., 2003), and depression (30-40%) (Sullivan et al., 2000). The heritability estimated from the GWAS is lower than the true heritability as GWAS can only capture the common variations (Manolio et al., 2009). The heritability of anxiety traits estimated from GAD-2 based GWAS is about 5.58% (n 175,163) (Levey et al., 2020), the heritability of schizophrenia explained by GWAS is approximately 23-24% (n 135,236) (Lam et al., 2019), and the heritability of depression estimated from GWAS is about 8.7% (n 480,359) (Wray et al., 2018). These results suggest that the polygenetic nature varies among these diseases. Second, gene expression patterns are highly similar between human brain and testis (Guo et al., 2005), which may result in the observation of significant enrichment of anxiety associations in testis. Finally, it is possible that the anxiety risk genes also play a role in testis. More work is needed to explain this interesting observation.
Our study highlighted the importance of TWAS (which integrates gene expression data with GWAS summary) in nominating risk genes for anxiety. However, genetics can only explain a small proportion of anxiety disorders/traits, which indicates that environmental factors are also involved in anxiety disorders/traits. Among the various environmental factors, psychosocial stress plays a vital role in anxiety (Bartlett et al., 2017). For example, early life adversity (Nemeroff CB, 2004;Maccari et al., 2014;Wiegand et al., 2021) is an important environmental risk factor for social anxiety disorder. The stress environmental factors like early life adversity would trigger the DNA methylation change (Provençal and Binder, 2015), and epigenome-wide association study (EWAS) have reported many genes that showed epigenetic changes in anxiety disorders, including CFAP46 (Ziegler et al., 2019), SLC43A2, and TNXB (Wiegand et al., 2021). Cognitivebehavioral therapy (CBT) (Barlow et al., 2000) is an useful psychosocial intervention for anxiety disorders and previous study has reported that the methylation level of IL1R1 was changed in subjects of panic disorder (one of the anxiety disorders) treated with CBT (Ziegler et al., 2019). These findings indicate that epigenetic regulation also plays an import role in anxiety disorders/traits.
Our study has several limitations. Firstly, we only included the GAD-2 score based GWAS summary statistics from the MVP project in our current study. We compared the phenotyping methods of this GWAS and other published anxiety GWASs (Otowa et al., 2009;Erhardt et al., 2011;Otowa et al., 2012;Walter et al., 2013;Davies et al., 2015;Otowa et al., 2016;Dunn et al., 2017;Stein et al., 2017;Meier et al., 2019;Levey et al., 2020;Purves et al., 2020), and found that the phenotyping methods were not identical across studies. We believed that it would be best to perform an independent replication study based on the same phenotyping method. Secondly, though several TWS genes are identified, more efforts are needed to uncover the roles of these identified risk genes in anxiety etiology. Thirdly, although TWAS could identify candidate risk genes for anxiety, the causal variants in the GWAS/TWAS significant loci could not be pinpointed by TWAS. More efforts are needed to reveal the genetic mechanisms and pathogenesis of the genes in the GWAS/TWAS significant loci in anxiety.
In summary, our study uncovered the gene-trait relationships of anxiety traits for the first time. We identified 26 TWS genes for anxiety. Our results provide novel insight into anxiety disease etiology and facilitate further mechanism studies and future drug development.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
XJL and JL conceived, designed, and supervised the whole study. XS, WL, LL, XL, JL, and JY conducted all the analyses, including SNP-gene expression weights, TWAS, conditional analyses, tissue, and cell-type enrichment analyses, and pathway analysis. JL drafted the manuscript, XJL oversaw the project and finalized the manuscript. All authors revised the manuscript critically and approved the final version.

FUNDING
This study was equally supported by the Distinguished Young Scientists grant of the Yunnan Province (202001AV070006) and Innovative Research Team of Science and Technology department of Yunnan Province (2019HC004) to XJL. Also was supported by the National Nature Science Foundation of China (31970561 to XJL), the CAS "Light of West China" Program to JL, and Yunnan Fundamental Research Projects (202001AT070099) to JL. One of the brain eQTL datasets used in this study were generated as part of the CommonMind Consortium supported by funding from Takeda Pharmaceuticals Company Limited, F. Hoffman-La Roche Ltd. and NIH grants R01MH085542, R01MH093725, P50MH066392, P50MH080405, R01MH097276, RO1-MH-075916, P50M096891, P50MH084053S1, R37MH057881 and R37MH057881S1, HHSN271201300031C, AG02219, AG05138 and MH06692. Brain tissue for the study was obtained from the following brain bank collections: the Mount Sinai NIH Brain and Tissue Repository, the University of Pennsylvania Alzheimer's Disease Core Center, the University of Pittsburgh NeuroBioBank and Brain and Tissue Repositories and the NIMH Human Brain Collection Core. CMC Leadership: Pamela Sklar, Joseph Buxbaum (Icahn School of Medicine at Mount Sinai), Bernie Devlin, David Lewis (University of Pittsburgh), Raquel Gur, Chang-Gyu Hahn (University of Pennsylvania), Keisuke Hirai, HiroyoshiToyoshiba (Takeda Pharmaceuticals Company Limited), Enrico Domenici, Laurent Essioux (F. Hoffman-La Roche Ltd.), Lara Mangravite, Mette Peters (Sage Bionetworks), Thomas Lehner, Barbara Lipska (NIMH). The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS.