Introduction
Anorexia nervosa (AN) is a severe eating disorder with a lifetime prevalence of 0.9–4%, characterized by extreme low body weight, fear of gaining weight, and compensatory weight-loss behaviors such as dietary restrictions, purging, and excessive exercise (Zipfel, Giel, Bulik, Hay, & Schmidt, Reference Zipfel, Giel, Bulik, Hay and Schmidt2015). Despite having one of the highest mortality rate of any psychiatric disorder and an increased risk of suicide (Bulik, Flatt, Abbaspour, & Carroll, Reference Bulik, Flatt, Abbaspour and Carroll2019; Chesney, Goodwin, & Fazel, Reference Chesney, Goodwin and Fazel2014; Mitchell & Peterson, Reference Mitchell and Peterson2020), relatively little is known about the biological mechanisms underlying AN, and effective, evidence-based treatments are scant, especially for adults (Watson & Bulik, Reference Watson and Bulik2013). Twin studies have established AN-heritability between 50% and 60%, indicating a considerable contribution of genetic factors to disease liability (Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019; Yilmaz, Hardaway, & Bulik, Reference Yilmaz, Hardaway and Bulik2015).
Genetic studies of AN provide evidence of both psychiatric and metabolic etiology. The largest AN genome-wide association study (GWAS) to date (N Cases = 16 992) uncovered eight AN-associated loci (Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019), and determined single-nucleotide polymorphism (SNP)-based heritability (h 2SNP) of 11–17%, similar to other psychiatric disorders. In addition, genetic correlations demonstrate significant shared genetic variation between AN and other psychiatric disorders, including major depressive disorder (MDD), obsessive-compulsive disorder, anxiety disorders (Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019), schizophrenia (Duncan et al., Reference Duncan, Yilmaz, Gaspar, Walters, Goldstein, Anttila and Bulik2017; Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019), alcohol use disorder (Munn-Chernoff et al., Reference Munn-Chernoff, Johnson, Chou, Coleman, Thornton, Walters and Agrawal2021), as well as with physical activity (Hübel et al., Reference Hübel, Gaspar, Coleman, Finucane, Purves, Hanscombe and Breen2019a; Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019). Significant negative genetic correlation between AN and anthropometric and metabolic traits [including body mass index (BMI), fat mass, obesity, type-2 diabetes, leptin, and insulin-related traits; Duncan et al., Reference Duncan, Yilmaz, Gaspar, Walters, Goldstein, Anttila and Bulik2017; Hübel et al., Reference Hübel, Gaspar, Coleman, Finucane, Purves, Hanscombe and Breen2019a; Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019] have also been observed, further indicating metabolic components to AN disease risk. Studies of AN polygenic risk scores indicate additional genetic associations of AN risk variants with anthropometric, behavioral, and psychiatric traits (Hübel et al., Reference Hübel, Abdulkadir, Herle, Loos, Breen, Bulik and Micali2021), and with weight trajectories even among healthy adults (Xu et al., Reference Xu, Johnson, Birgegård, Jordan, Kennedy, Landén and Huckins2021).
Although GWASs may provide powerful insights into genetic associations with AN, they cannot pinpoint gene- or tissue-level associations. Therefore, in this study we apply transcriptomic imputation (TI) to identify tissue-specific gene associations with AN. TI approaches leverage gene expression predictor models derived from large, well-curated gene expression datasets [e.g. the Genotype-Tissue Expression project (GTEx) (Lonsdale et al., Reference Lonsdale, Thomas, Salvatore, Phillips, Lo, Shad and Moore2013), CommonMind Consortium (CMC) (Fromer et al., Reference Fromer, Roussos, Sieberts, Johnson, Kavanagh, Perumal and Sklar2016; Huckins et al., Reference Huckins, Dobbyn, Ruderfer, Hoffman, Wang, Pardiñas and Stahl2019), and Depression Genes and Networks (DGN) (Battle et al., Reference Battle, Mostafavi, Zhu, Potash, Weissman, McCormick and Koller2014)]. These models may be applied to predict genetically regulated gene expression (GReX) in large genotyped cohorts, without the need to collect tissue samples (Huckins et al., Reference Huckins, Dobbyn, Ruderfer, Hoffman, Wang, Pardiñas and Stahl2019).
Observational genetic studies are powerful approaches to detect disease-associated variants, but cannot address subthreshold or prodromal disease states, and do not speak to clinical relevance. Therefore, we leverage a phenome-wide association study (pheWAS) design to probe clinical outcomes with AN-associations. pheWAS effectively query the full electronic health record (EHR) to identify diagnoses and traits associated with a gene or variant (Denny et al., Reference Denny, Ritchie, Basford, Pulley, Bastarache, Brown-Gentry and Crawford2010; Pendergrass et al., Reference Pendergrass, Brown-Gentry, Dudek, Torstenson, Ambite, Avery and Ritchie2011; Smoller, Reference Smoller2018), and have been used to replicate associations and identify new pleiotropic consequences of GWAS variants, including psychiatric disorders (Denny et al., Reference Denny, Bastarache, Ritchie, Carroll, Zink, Mosley and Roden2013; Leppert et al., Reference Leppert, Millard, Riglin, Davey Smith, Thapar, Tilling and Stergiakouli2020; Zheutlin et al., Reference Zheutlin, Dennis, Karlsson Linnér, Moscati, Restrepo, Straub and Smoller2019). EHRs contain vast amount of longitudinal health data, including diagnoses, medications, laboratory tests, vital signs, and family medical history, and, coupled with genetic data through hospital-based biobanks, can provide an understanding of the clinical spectrum of disease and disease progression across the lifetime of the patient. Exploring the associations of AN-genetics with clinical phenotypes has the potential to clarify how some of these GWAS variants functionally contribute to AN disease risk, symptomatology, and clinical presentation.
Here, we explored the associations between AN-GReX and the clinical phenome. We performed TI using S-PrediXcan on the most recent PGC-ED AN GWAS (Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019) to first find GReX associated with AN, and then tested for clinical associations of these genes with structured EHR-encoded phenotypes using pheWAS. We further investigated the effects of BMI and sex on the GReX-phenotype associations through stratification of biobank individuals. Understanding the clinical associations of aberrant gene expression across the phenome may clarify the biological mechanisms of relevant AN GWAS risk variants.
Methods
Transcriptomic imputation
We performed TI using S-PrediXcan (Barbeira et al., Reference Barbeira, Dickinson, Bonazzola, Zheng, Wheeler, Torres and Im2018) on the largest available summary statistics of AN (16 992 cases/55 525 controls) (Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019). We tested for the association of GReX using available GTEx (Lonsdale et al., Reference Lonsdale, Thomas, Salvatore, Phillips, Lo, Shad and Moore2013), DGN (Battle et al., Reference Battle, Mostafavi, Zhu, Potash, Weissman, McCormick and Koller2014), and CMC (Fromer et al., Reference Fromer, Roussos, Sieberts, Johnson, Kavanagh, Perumal and Sklar2016; Hoffman et al., Reference Hoffman, Bendl, Voloudakis, Montgomery, Sloofman, Wang and Roussos2019) predictor models (Barbeira et al., Reference Barbeira, Dickinson, Bonazzola, Zheng, Wheeler, Torres and Im2018; Gamazon et al., Reference Gamazon, Wheeler, Shah, Mozaffari, Aquino-Michaels, Carroll and Im2015; Huckins et al., Reference Huckins, Dobbyn, Ruderfer, Hoffman, Wang, Pardiñas and Stahl2019) across 50 tissues with AN case–control status. We established two thresholds for significance; first, Bonferroni correction for all genes tested within each tissue (p tissue, online Supplementary Table S1), and second, correcting for all tissues and genes tested (p Experiment = 3.75 × 10−8). We performed a gene-set analysis for 53 significant AN-genes using FUMA (v1.3.6) GENE2FUNC (Watanabe, Taskesen, van Bochoven, & Posthuma, Reference Watanabe, Taskesen, van Bochoven and Posthuma2017), including all S-PrediXcan genes (N = 28 454) as the background genes. We defined significant gene sets using p FDR < 0.05.
BioMe™
The Mount Sinai BioMe™ Biobank includes genotype and EHR data from 31 704 individuals (online Supplementary Table S2). Individuals were genotyped on the Illumina Global Screening Array; quality control (QC), and imputation of the genotyping data for BioMe™ is described elsewhere (Zheutlin et al., Reference Zheutlin, Dennis, Karlsson Linnér, Moscati, Restrepo, Straub and Smoller2019). Ancestry QC of BioMe™ individuals using principal component analysis in PLINK (Chang et al., Reference Chang, Chow, Tellier, Vattikuti, Purcell and Lee2015; Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender and Sham2007) is described in the online Supplemental Methods. After QC, a total of 31 585 individuals were available for analysis.
pheWAS
We calculated GReX for all S-PrediXcan p tissue significant genes (N = 53) across 45 tissues from GTEx, DGN, and CMC in BioMe™ using PrediXcan-2 (Barbeira et al., Reference Barbeira, Dickinson, Bonazzola, Zheng, Wheeler, Torres and Im2018; Gamazon et al., Reference Gamazon, Wheeler, Shah, Mozaffari, Aquino-Michaels, Carroll and Im2015). We excluded sex-specific tissues (ovary, uterus, vagina, prostate, and testis) from this analysis, and genes with GReX variance (gVAR) less than 0.002 (online Supplementary Table S3C).
Logistic regressions between AN-GReX and categorical phenotypes were performed using the pheWAS R package (Carroll, Bastarache, & Denny, Reference Carroll, Bastarache and Denny2014), adjusting for age, sex, and the first five genotype-derived principal components. pheWAS were run individually for each BioMe™ ancestry cohort, and results were meta-analyzed using an inverse-variance based approach in METAL (Willer, Li, & Abecasis, Reference Willer, Li and Abecasis2010). We required at least 10 occurrences of each phenotype in each population group for inclusion in the analysis, and an overall effective sample size N eff > 100 (Equation (1)). We excluded associations with significant Cochran's Q-test for heterogeneity scores across cohorts (pHet > 0.001).
Effective sample size (N eff):
‘Encounter Diagnoses’ were recorded at each patient visit using the International Classification of Disease (ICD) coding system. Phecodes were assigned from Encounter Diagnoses by grouping ICD-9 and ICD-10 diagnostic codes (Wu et al., Reference Wu, Gifford, Meng, Li, Campbell, Varley and Wei2019). We defined cases as individuals with at least two counts of a code. Those with zero counts were considered ‘controls’, and those with only one count were set to missing. After QC, our dataset included 2178 unique Encounter Diagnosis codes and 1093 unique phecodes. Due to the high correlation between Encounter Diagnosis and phecode files, we combined all results from both pheWAS and performed an FDR correction in R to determine significance (p FDR < 0.05).
In addition to diagnostic codes, BioMe™ EHR data include allergies, vital signs (weight, height, blood pressure, pulse, pulse oximetry, respirations, and temperature), lab results, family history, personal history, medication use, obstetrics and gynecology outcomes (OB/GYN), and other social and behavioral traits. Summaries, descriptive statistics, and QC descriptions for these phenotypes can be found in online Supplementary Methods. We tested for AN-GReX-pheWAS associations within each category, establishing significance using a Bonferroni correction (Equation (2)).
Bonferroni-correction for EHR data at tissue- and experiment-significance thresholds:
Stratification by BMI
Given the nature of AN symptomology, we tested whether our significant AN-GReX associations varied based on BMI. We stratified BioMe™ individuals into three BMI categories: low (<first BMI quartile, N = 6976), mid (first—third BMI quartiles, N = 13 972) and high (>third BMI quartile, N = 9350; online Supplementary Fig. S1A). Quartiles were defined based on ancestry- and sex-specific BMI distributions within our data (online Supplementary Fig. S1B; Table S2; Supplemental Methods). We repeated our pheWAS within these categories, for all AN-GReX associations reaching significance in our overall analysis. BMI-stratified pheWAS were adjusted for age, sex, and the first five genotype-derived PCs.
Stratification by sex
To determine whether our associations varied based on sex, we stratified the BioMe™ cohort by females (N = 17 968) and males (N = 12 417), with sex determined by genotype, and repeated our AN-GReX pheWAS analyses within each group. Here, we included sex-specific tissues (ovary, uterus, vagina, prostate, and testis) for the respective sex. Sex-stratified pheWAS were adjusted for age and the first five genotype-derived PCs.
Testing for hidden case contamination
It is possible that some of the clinical associations within our study stem from undiagnosed ED and AN cases, rather than a direct effect of gene expression on phenomic expression. We term this ‘diagnostic contamination’. To test whether this effect may drive the associations we observe, we simulated the effects of diagnostic contamination occurring at rate p within our biobank sample, with effect β on GReX.
Full derivations and simulations are shown in our online Supplementary Material; briefly, we derive the expression among controls, and cases contaminated at rate p; the difference in expression between the two groups; the expected variance among cases, and pooled across all samples, and the expected statistical significance (T-score, and p value). Next, we simulated gene expression distributions for (i) 1000 cases and 1000 controls; (ii) 1000 cases and 10 000 controls; (iii) 1000 cases and 30 000 controls; each for a range of p (0.1, 0.5, 1, 2, 5, 10, 20, 30, 40, 50%) and β (1/10, 1/5, 1/4, 1/3, 1/2, 1, 2, 3, 4, 5, 10) values. We repeated our simulations 10 000 times at each case–control proportion- p- β combination, and demonstrated that our formulae accurately estimate the desired values (online Supplementary Fig. S2).
In order to test whether diagnostic contamination may account for the associations observed within our pheWAS, we calculated the expected impact of diagnostic contamination for two genes (NCKIPSD-Aorta; SEMA3F-Spinal Cord) with the largest effect sizes β observed in our S-PrediXcan analysis, under two scenarios:
(1) Diagnostic contamination occurs among our cases only, assuming the normal population rate for AN: 0.9% among women; 0.3% among men.
(2) Contamination occurs at significantly higher levels than might be expected in the population (0.6, 1.2, 3, 6%), and all of these samples fall into our case category.
For each scenario, we repeated our calculation using 500 and 1000 cases, matched with 1000, 5000, and 30 000 controls.
Results
S-PrediXcan identifies gene–tissue associations with AN
We applied S-PrediXcan to PGC-AN GWAS summary statistics, and identified 218 gene–tissue associations, including 53 unique genes across 12 loci (Fig. 1). Two loci on chromosomes 2 and 3 reached experiment-wide significance (p < 3.75 × 10−8), while 10 loci reached tissue-specific significance (Fig. 2, Table 1; online Supplementary Table S3A). Our S-PrediXcan AN-genes are significantly enriched in gene sets for traits including inflammatory bowel disease (IBD: Crohn's and ulcerative colitis) (p adj = 1.3 × 10−25), sleep duration (p adj = 6.9 × 10−17), blood protein levels (p adj = 3.2 × 10−16), intelligence (p adj < 1.5 × 10−6), regular attendance at a religious group, gym, or sports club (p adj<1.46 × 10−10), mood instability (p adj = 3.3 × 10−9), AN (p adj = 7.4 × 10−8), and body fat distribution (P adj < 5.0 × 10−7) (online Supplementary Table S4A and B).
Graphical depiction of S-PrediXcan and PrediXcan TI and pheWAS analyses. (a) We used S-PrediXcan predictor models for 50 different tissue types to impute GReX in Watson et al. (Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019) AN GWAS and found 53 genes whose GReX was associated with AN. We then (b) imputed GReX for those 53 AN-genes in our BioMe™ cohort and (c) performed a pheWAS across available EHR phenotypes. The pheWAS analyses were run within each ancestral population and then meta-analyzed using an inverse-variance approach in METAL. Secondary analyses included stratifying individuals in BioMe™ by BMI and sex, and running the pheWAS analyses within each stratification group.
S-PrediXcan results for PGC-ED AN GWAS (Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019) (N Cases = 16 992, N Controls = 55 525). S-PrediXcan TI of the PGC-ED AN GWAS summary statistics to determine GReX, with tests for association of GReX with AN disease status. Manhattan plot of S-PrediXcan gene–tissue associations with AN for 50 tissues. Each point represents a different gene–tissue association result; i.e. the same gene may have multiple points within a peak. Experiment-wide significance threshold of p < 3.75 × 10−8 (solid line); tissue-specific significance threshold of p < 2.45 × 10−5 (dotted line).
S-PrediXcan loci results
Summary of S-PrediXcan of PGC-ED AN GWAS results. Genes associated with AN from each locus are shown, as well as the top finding Z score, p value, gene and tissue for each locus. p values shaded in bold indicate experiment-wide significant loci (p < 3.75 × 10−8).
AN-GReX is associated with autoimmune and autoinflammatory diseases
We performed a pheWAS to identify clinical associations with dysregulation of our AN-genes. We identified 16 FDR-significant gene–tissue associations with four Phecode- and Encounter Diagnosis phenotypes (FDR < 0.05): type 1 diabetes (N Cases = 408), celiac disease (N Cases = 63), peptic ulcer (N Cases = 254), and unspecified immunodeficiency (N Cases = 55) (Fig. 3; online Supplementary Table S5). These associations are driven in part by the MHC-gene CLIC1; predicted upregulation of CLIC1 in spleen was associated with celiac disease and type 1 diabetes (p = 1.30 × 10−11, p = 5.08 × 10−20, respectively) in the overall cohort. We did not observe any significant effect of patient BMI group on these associations (online Supplementary Fig. S3A and C). Upregulation of CLIC1 in females was significantly associated with type 1 diabetes (N Cases = 230) and celiac disease (N Cases = 53) (Females-Spleen-CLIC1, p < 1.49 × 10−9) (online Supplementary Fig. S3B and D) .
AN-GReX associations with BioMe™ Diagnosis codes. GReX-Tissue-Phenotype associations for (a) encounter diagnosis ICD-10 codes (N = 2178) and (b) phecodes (N = 1093) for the BioMe™ cohort (N = 30 585). Diagnosis codes are plotted along the x-axis and grouped by category, with the −log(10) p value associations along the y-axis. FDR-significant diagnosis codes are labeled (FDR-adjusted p < 0.05). FDR-significant p value threshold p = 3.4 × 10−7 (blue dashed line).
pheWAS AN-GReX tissue associations with cholesterol phenotypes
pheWAS associations with highest, lowest, and mean total cholesterol, LDL cholesterol and HDL cholesterol measures for all stratification groups: overall (N = 23 661), high-BMI (N = 6910), mid-BMI (N = 11 243), low-BMI (N = 5324), females (N = 13 907), and males (N = 9669). AN-GReX association with phenotypes of (A) highest recorded total cholesterol (mg/dl), (B) highest recorded LDL cholesterol (mg/dl), (C) highest recorded HDL cholesterol (mg/dl), (D) mean total cholesterol (mg/dl), (E) mean LDL cholesterol (mg/dl), (F) mean HDL cholesterol (mg/dl), (G) lowest recorded total cholesterol (mg/dl), (H) lowest recorded LDL cholesterol (mg/dl), and (I) lowest recorded HDL cholesterol (mg/dl). Bonferroni experiment-wide significant associations are marked in bold (p Experiment = 0.05/(9 phenotypes × 45 tissues) = 1.23 × 10−4). Bonferroni tissue-specific value threshold p < 0.0056 (0.05/9 phenotypes).
BMI moderates the effect of MGMT-GReX on cholesterol levels
Next, we looked at the association of AN-GReX and measures of highest recorded, lowest recorded and mean total cholesterol (mg/dl), high-density lipoprotein (HDL) cholesterol (mg/dl), and low-density lipoprotein (LDL) cholesterol (mg/dl) (Table 2). In the full cohort, RPS26 and SUOX GReX were significantly associated with the highest HDL measures, and MGMT GReX was associated with the highest cholesterol at within-tissue significance (N = 23 357) (p < 5.1 × 10−4) (Fig. 4C). Associations of AN-GReX with lipid measures varied across BMI-stratified groups; predicted downregulation of MGMT in the stomach, liver, esophagus, and cells was associated with the highest total cholesterol and LDL levels, and in dorsolateral prefrontal cortex and hippocampus with mean cholesterol and LDL among high-BMI individuals (N High BMI = 6910) (High-BMI-Stomach-MGMT, p = 1.55 × 10−7) (Fig. 4C; online Supplementary Table S6).
Context-specific pheWAS associations. (a) Concordance of context-specific experiment-wide significant AN-GReX clinical associations with AN direction of effect. We compared the direction of effect (DoE) for each experiment-wide gene–tissue-pheWAS association with the DoE of that gene–tissue pair for AN from our S-PrediXcan results (online Supplemental Methods). For those phenotypes concordant with AN, this may indicate that genetic regulation of those AN-genes is more similar to individuals with AN in individuals with those clinical outcomes. (b) Schematic of the proportion of concordance of AN-GReX pheWAS associations with AN S-PrediXcan associations. Associations with similar direction of effect to AN (green) identified as ‘concordant’, associations with opposite direction of effect (purple) identified as ‘discordant’. (c) Context-specific associations of AN-GReX with lipid phenotypes of highest recorded, lowest recorded, and mean measures of total cholesterol (mg/dl), HDL cholesterol (mg/dl), and LDL cholesterol (mg/dl). Experiment-wide significant threshold set at p = 0.05(9 phenotypes × 45 tissues) = 1.2 × 10−4. Tissue-specific threshold set at 0.05/(9 phenotypes) = 0.0056. Context-specific associations of AN-GReX with pain location for (d) experiment-wide significant associations (p = 0.05/(99 phenotypes × 45 tissues) = 1.1 × 10−5), (e) all context-specific associations (experiment-wide and tissue-specific) with generalized pain, and (f) with foot pain. Tissue-specific threshold for pain location set at 0.05/(99 phenotypes) = 5.0 × 10−4.
AN-GReX is associated with tobacco use
Multiple genes were associated with categorical and continuous measures of tobacco use in the overall cohort at within-tissue significance (online Supplementary Fig. S4, Table S7, Supplementary Results). We found moderating effects of patient BMI on tobacco use: among High-BMI groups, different genes were associated with the number of years of tobacco use at experiment-wide significance (N High BMI = 8129) (High-BMI-Pituitary-CCDC36, p = 4.2 × 10−6; High-BMI-Colon, Transverse-P4HTM, p = 3.4 × 10−5; High-BMI-Stomach-P4HTM, p = 2.2 × 10−5) (online Supplementary Fig. S4, Table S7).
AN-GReX is associated with measures of pain score and pain location
We tested for association between AN-GReX and (1) pain location, (2) pain score overall, and (3) pain score by body location. AN-GReX was associated with multiple measures of pain and location across all stratification groups (Fig. 4D–F; online Supplementary Table S8). Our results revealed a large number of AN-GReX-pain associations specific to the Low-BMI category (Fig. 4D–F; online Supplementary Table S8); in particular, we note 60 gene–tissue associations with presence of foot pain (N Cases = 433), one of which passed experiment-wide significance (Low-BMI-Small intestine-GPX1, p = 7.6 × 10−6), which may indicate a propensity to excess exercise, or exercise-related injuries among these individuals.
Upregulation of CLIC1 is associated with glucagon medication
Upregulation of CLIC1 was associated with glucagon, a hormone used to treat severe hypoglycemia, in the overall cohort (N Cases = 129) (Overall-Spleen-CLIC1, p = 4.20 × 10−9) and in females (N Cases = 78) (Females-Spleen-CLIC1, p = 1.4 × 10−8) (online Supplementary Fig. S5, Table S9). We additionally find significant associations between the upregulation of CLIC1 and insulin medications in females (Females-Subcutaneous adipose-CLIC1, p < 9.6 × 10−7).
Identifying potential case contamination effects
Our BioMe™ patients have not all been explicitly assessed for eating disorders, and information regarding eating disorder diagnoses and assessments earlier in life may be omitted from the records due to incomplete clinical history assessments. Therefore, it is possible that diagnostic contamination within some of our sample is responsible for the associations observed within our data. We tested whether such contamination may drive the associations observed in our study, assuming the following different contamination scenarios (online Supplementary Table S10), for two genes with large S-PrediXcan effect sizes: (i) that contamination occurs among cases, assuming the highest common estimate for AN prevalence (0.9% among women, 0.3% among men); (ii) that contamination occurs within our biobank, assuming the highest common estimate for AN prevalence (0.6%, 190 cases), and that all of these samples fall into our ‘case’ category.
Our model indicates that diagnostic contamination among cases is unlikely to result in significant pheWAS associations in our data. Contamination occurring at population ED prevalence (0.6%) did not result in significant associations for any of the case–control scenarios in our model (online Supplementary Table S10); nor did contamination at 1.2% or 3%. Assuming 6% contamination resulted in potentially significant associations (estimated p < 2.8 × 10−8; see Methods and online Supplementary Methods for more details).
Discussion
Our analysis identified novel AN-genes associated with metabolic, anthropometric, autoimmune, and psychiatric phenotypes (online Supplementary Table S4B). In particular, the experiment-wide significant locus on chromosome 3 overlaps with a known GWAS peak for IBD (de Lange et al., Reference de Lange, Moutsianas, Lee, Lamb, Luo, Kennedy and Barrett2017), and includes genes associated with Crohn's and ulcerative colitis (online Supplementary Table S4B) (Hedman et al., Reference Hedman, Breithaupt, Hübel, Thornton, Tillander, Norring and Bulik2019; Mårild et al., Reference Mårild, Størdal, Bulik, Rewers, Ekbom, Liu and Ludvigsson2017; Raevuori et al., Reference Raevuori, Haukka, Vaarala, Suvisaari, Gissler, Grainger and Suokas2014; Zerwas et al., Reference Zerwas, Larsen, Petersen, Thornton, Quaranta, Koch and Bulik2017), in line with earlier AN GWAS showing genetic overlap with autoimmune disease (Gibson & Mehler, Reference Gibson and Mehler2019). Variants in one of our significant S-PrediXcan genes, GPR75, have been recently shown to have a protective effect against obesity, and have been associated with lower body weight overall (Akbari et al., Reference Akbari, Gilani, Sosina, Kosmicki, Khrimian, Fang and Lotta2021).
To our knowledge, our results include the first AN-association within the MHC locus (CLIC1, chloride intracellular 1), a region that has been associated with many other psychiatric (Ripke et al., Reference Ripke, Neale, Corvin, Walters, Farh, Holmans and O'Donovan2014; Stahl et al., Reference Stahl, Breen, Forstner, McQuillin, Ripke, Trubetskoy and Sklar2019) and autoimmune disorders. In particular, CLIC1 has previously been identified as associated with schizophrenia (The Autism Spectrum Disorders Working Group of The Psychiatric Genomics Consortium, 2017), autism (The Autism Spectrum Disorders Working Group of The Psychiatric Genomics Consortium, 2017), MDD (Zhu et al., Reference Zhu, Zhu, Liu, Shi, Shen, Yang and Liang2019), post-traumatic stress disorder (Marchese et al., Reference Marchese, Cancelmo, Diab, Cahn, Aaronson, Daskalakis and Feder2021), neuroticism (Baselmans et al., Reference Baselmans, Jansen, Ip, van Dongen, Abdellaoui, van de Weijer and Bartels2019) and depressive phenotypes (Baselmans et al., Reference Baselmans, Jansen, Ip, van Dongen, Abdellaoui, van de Weijer and Bartels2019). Importantly, CLIC1 variants have also been associated with complement component C4 and C3 protein levels in the blood (Yang et al., Reference Yang, Sun, Gao, Tan, Zhang, Hu and Mo2012), which through the complement cascade (Sekar et al., Reference Sekar, Bialas, de Rivera, Davis, Hammond, Kamitaki and McCarroll2016) are involved in immunological functions of pathogen clearance and in synaptic pruning and neuronal connectivity (Stephan, Barres, & Stevens, Reference Stephan, Barres and Stevens2012). CLIC1 encodes a chloride ion channel protein involved in many necessary cellular functions including the regulation of cell membrane potential, and the proliferation and differentiation of cells (Li et al., Reference Li, Mao, Wang, Chen, Wang, Zhai and Chen2018), including a role in axonal outgrowth of neurons (Carlini et al., Reference Carlini, Verduci, Cianci, Cannavale, Fenoglio, Galimberti and Mazzanti2020).
We also identified multiple genes that have been previously associated with decreased sleep and insomnia phenotypes (including a core circadian clock gene ARNTL), and with a range of psychiatric disorders (Alloy, Ng, Titone, & Boland, Reference Alloy, Ng, Titone and Boland2017; Boivin, Reference Boivin2000; Bunney & Bunney, Reference Bunney and Bunney2013; Bunney et al., Reference Bunney, Li, Walsh, Stein, Vawter, Cartagena and Bunney2015; Etain et al., Reference Etain, Jamain, Milhiet, Lajnef, Boudebesse, Dumaine and Bellivier2014; Huckins et al., Reference Huckins, Dobbyn, McFadden, Wang, Ruderfer, Hoffman and Stahl2017; Karatsoreos, Reference Karatsoreos2014; Kripke, Mullaney, Atkinson, & Wolf, Reference Kripke, Mullaney, Atkinson and Wolf1978; Mansour et al., Reference Mansour, Talkowski, Wood, Chowdari, McClain, Prasad and Nimgaonkar2009; Melo, Abreu, Linhares Neto, de Bruin, & de Bruin, Reference Melo, Abreu, Linhares Neto, de Bruin and de Bruin2017; Meyrer, Demling, Kornhuber, & Nowak, Reference Meyrer, Demling, Kornhuber and Nowak2009; Murray, Allen, Trinder, & Burgess, Reference Murray, Allen, Trinder and Burgess2002), including eating disorders (Allison, Spaeth, & Hopkins, Reference Allison, Spaeth and Hopkins2016), as well as with satiety and hunger (Arble, Bass, Laposky, Vitaterna, & Turek, Reference Arble, Bass, Laposky, Vitaterna and Turek2009; Herpertz et al., Reference Herpertz, Albers, Wagner, Pelz, Köpp, Mann and Hebebrand2000).
In a patient population, expression of AN genes is associated with metabolic and autoimmune phenotypes
Unlike GWAS, which include carefully constructed case–control cohorts, pheWAS encompass all individuals within a healthcare system, including patients with subthreshold or partial presentations of a disorder, or individuals with commonly comorbid or co-diagnosed conditions. Importantly, because individuals are not ascertained for any specific disorder, they represent a more comprehensive clinical picture of the comorbidities and symptomatology associated with AN gene expression.
Examining the consequences of aberrant predicted gene expression among these patients (i.e. testing for GReX associations) may reveal clinical and biological consequences of these genes; for example, studying whether AN-associated genes have anthropometric and metabolic consequences among adults with no evidence for previous AN- or ED-diagnoses may disentangle whether certain endophenotypes present as a cause or consequence of AN. Food avoidance and restriction may arise due to gastrointestinal (GI) complaints and distress that provoke these behaviors and precede development of AN (Zucker & Bulik, Reference Zucker and Bulik2020). Similarly, autoimmune disorders of the GI tract, such as celiac disease and Crohn's disease, show a bidirectional relationship with AN, with previous diagnosis of a GI-associated autoimmune disorders increasing the risk of AN and vice versa (Hedman et al., Reference Hedman, Breithaupt, Hübel, Thornton, Tillander, Norring and Bulik2019). Our pheWAS results of AN-GReX associations with GI symptoms such as abdominal pain, ascites, and peptic ulcer, as well as GI-related autoimmune disorders such as celiac disease, suggest AN-GReX may contribute directly to these diseases and phenotypes, and that food aversive behaviors and gastric distress may be genetically regulated in these individuals, rather than occurring as a consequence of AN. We found a sex-specific association of CLIC1-GReX with celiac disease in females (online Supplementary Fig. S3B), indicating sex-specific regulation of AN-genes may be contributing to the disparity in autoimmune diagnoses and symptomatology between the sexes (Fuchs et al., Reference Fuchs, Kurppa, Huhtala, Collin, Mäki and Kaukinen2014; Ludvigsson & Murray, Reference Ludvigsson and Murray2019). Our main pheWAS diagnosis results of celiac disease and type 1 diabetes further show concordance with AN in the direction of effect of the tissue GReX, indicating a similarity in the genetic regulation of AN-genes between patients with these disorders and individuals with AN (online Supplementary Fig. 4A; Supplemental Results). Although these associations could be the consequence of undiagnosed AN-individuals in our biobank, they more likely reflect real biological associations of expression of those particular genes with the phenotype.
Our results further confirm the contribution of metabolic factors to AN etiology; we see very robust association of AN-GReX with type 1 diabetes, the hyperglycemic hormone glucagon and various forms of insulin. The association of AN with metabolic traits and abnormalities has been fairly well established, including with insulin resistance (r g = −0.29), fasting insulin (r g = −0.24), leptin (r g = −0.26) and type 2 diabetes (r g = −0.22), and a positive genetic correlation with HDL cholesterol (r g = 0.21) (Adams, Reay, Geaghan, & Cairns, Reference Adams, Reay, Geaghan and Cairns2021; Ilyas et al., Reference Ilyas, Hübel, Stahl, Stadler, Ismail, Breen and Kan2019; Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019). Notably, our results point to a similar role of aberrant glycemic regulation in the etiology of AN. Future analyses including EHR-derived lab results (LabWAS) studies may further elucidate AN-genes associated with abnormal metabolic regulation and clinical features.
Context-specific associations reveal a role for BMI in the regulation of AN-gene expression
We stratified our pheWAS analyses by BMI in order to observe whether the effects of predicted AN-gene expression on clinical outcomes were moderated by BMI status. Understanding how BMI influences associations between AN-GReX and the clinical phenome may give us further insight into the biological mechanisms leading to those outcomes and the conferral of AN disease risk. The context-specific associations we see between BMI status and the association of AN-GReX with multiple phenotypes indicate the potential differences in AN-gene regulation in an environment of higher or lower BMI.
Among our BMI-stratified pheWAS analyses, we identified a number of associations between AN-GReX and foot pain, exclusive to the Low-BMI group. We hypothesize that these associations may stem from exercise-related injuries among this group. Multiple genes in our study, many of which were associated with foot pain in Low-BMI individuals, have been previously associated with regular gym attendance (online Supplementary Table S4B, p adj = 3.4 × 10−10), or measures of physical activity. Excessive and compulsive exercise is a behavior often seen in individuals with AN (Bulik, Reference Bulik and Thompson2004; Shroff et al., Reference Shroff, Reba, Thornton, Tozzi, Klump, Berrettini and Bulik2006), and evidence is suggestive that hyperactivity increases risk of chronic AN (Achamrah, Coëffier, & Déchelotte, Reference Achamrah, Coëffier and Déchelotte2016; Strober, Freeman, & Morrell, Reference Strober, Freeman and Morrell1997). Similarly, studies have shown a genetic correlation between AN and physical activity (r g = 0.17) (Hübel et al., Reference Hübel, Gaspar, Coleman, Hanscombe, Purves, Prokopenko and Breen2019b; Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019). These genes may reflect a genetic liability to compulsive physical activity. Given the known association of AN disease with low bone mineral density and propensity for bone fracture (Solmi et al., Reference Solmi, Veronese, Correll, Favaro, Santonastaso, Caregaro and Stubbs2016), our results may reflect the result of genes associated with compulsive exercise and bone mineral density contributing to increased osteoarticular pain.
We saw distinct modulation by BMI of MGMT-GReX-cholesterol associations in our cohort. Variants in the MGMT gene have been shown to be associated with AN (Watson et al., Reference Watson, Yilmaz, Thornton, Hübel, Coleman, Gaspar and Bulik2019), as well as anthropometric phenotypes such as waist-to-hip ratio and body height (Kichaev et al., Reference Kichaev, Bhatia, Loh, Gazal, Burch, Freund and Price2019). CTNNB1, which was associated with both highest and lowest weight measures in our cohort, distinctly with measures of highest weight in High-BMI individuals, has been previously shown to be associated with lean body mass (Hübel et al., Reference Hübel, Gaspar, Coleman, Finucane, Purves, Hanscombe and Breen2019a). We saw distinct associations of PROS1-GReX with lowest weight in Low-BMI individuals, which is a gene involved in many biological processes, including inflammation (Suleiman, Négrier, & Boukerche, Reference Suleiman, Négrier and Boukerche2013). Our results appear to be due to the effect of BMI on AN-GReX, rather than the direct effect of BMI/weight on the phenotypes (online Supplementary Fig. S4B). Future studies are needed to assess the specific role of BMI on AN-gene regulation.
Although the hypothesis-free, phenome-wide design of our study allows for powerful detection of clinical and biological associations of AN-risk genes, the same design also bears some notable limitations and caveats. One key caveat is the lack of diagnostic detail available for the patients in our study. EHR analyses leverage large, existing data sets to rapidly amass cohorts for analysis, to yield insights into whole phenome consequences of genotype and GReX associations; however, the scale and scope of these studies precludes deep phenotyping or performing clinical interviews. This lack of diagnostic precision may arise from a number of factors.
First, we make use of ICD codes within the medical record in order to infer diagnoses; since these are primarily used by clinicians for billing purposes, they likely provide an imperfect proxy for true disease state. In order to mitigate spurious results stemming from imperfectly assigned codes, we require ‘cases’ to have at least two instances of an ICD code; individuals with only one code were set to missing. Furthermore, for each code we restricted our meta-analysis to include only phenotypes with at least 10 cases and 10 controls, and required a total effective sample size >100 for inclusion in our final analysis.
Second, it is possible that our patients regularly receive treatment at multiple different hospital systems; as such, we may be capturing only partial data for each of our patients. In order to mitigate this, we restricted our analysis to patients with multiple data points within our EHR. Future analyses that seek to harmonize or meta-analyze patient data across EHR (e.g. NYC-CDRN, PsycheMERGE) are ongoing, and will be vital to disentangling this effect further.
Importantly, our patients have not all been explicitly assessed for EDs, and information regarding ED diagnoses and assessments earlier in life may be omitted from the records due to incomplete clinical history assessments. Therefore, it is possible that case-contamination within some of our sample is responsible for the associations observed within our data. In order to address this, we performed a simple experiment to simulate the effect sizes expected if undiagnosed AN contamination drives our result. We estimated the expected association statistics observed due to possible levels of contamination among pheWAS cases. Our model indicates that diagnostic contamination at 0.06, 1.2, or 3% is insufficient to account for the effect sizes observed within our pheWAS analysis. Moreover, since these genes are selected due to their large S-PrediXcan effect sizes, we expect that the contamination effects observed for these two genes will be greater than any others in our study; as such, we do not believe that our findings are attributable to hidden case contamination.
Our results illustrate the clinical outcomes associated with differences in AN-gene expression. Characterization of the phenotypic associations with AN-gene expression in a clinical setting can give us more insight into the biological mechanisms underlying AN and, consequently, how to diagnose and treat the disorder. By understanding the associations of AN-gene expression with symptomatology, prodromal or subthreshold disease states, we may gain insights into the biology of the disease, and perhaps identify therapeutic targets and opportunities for clinical interventions. For example, if GI complaints are truly the consequence of aberrant AN-gene expression, and contribute to disordered eating due to GI distress, treatment of those symptoms may help alleviate other AN symptoms or prevent development of later AN (Riedlinger et al., Reference Riedlinger, Schmidt, Weiland, Stengel, Giel, Zipfel and Mack2020; Wiklund et al., Reference Wiklund, Kuja-Halkola, Thornton, Hübel, Leppä and Bulik2019). An understanding of the clinical associations of AN-gene expression can further augment the definition of AN, and could allow clinicians to more broadly identify individuals at greater risk of AN, or those who present with symptom constellations that do not yet meet the established diagnostic threshold.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291721004554.
Acknowledgements
We are deeply grateful for the mentorship of Pamela Sklar, whose guidance and wisdom we miss daily. We strive to continue her legacy of thoughtful, innovative, and collaborative science.
Financial support
JJ and LMH were supported by funding from the Klarman Family Foundation (2019 Eating Disorders Research Grants Program) and the NIMH (R01MH118278). CMB is supported by NIMH (R01MH120170; R01MH124871; R01MH119084; R01MH118278; and R01MH124871); Brain and Behavior Research Foundation Distinguished Investigator Grant; Swedish Research Council (Vetenskapsrådet, award: 538-2013-8864); Lundbeck Foundation (Grant no. R276-2018-4581). This study was supported in part through the resources and staff expertise provided by the Charles Bronfman Institute for Personalized Medicine and The BioMe TM Biobank Program at the Icahn School of Medicine at Mount Sinai. Research reported in this paper was supported by the Office of Research Infrastructure of the National Institutes of Health under award numbers S10OD018522 and S10OD026880. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.
Conflict of interest
CMB reports: Shire (grant recipient, Scientific Advisory Board member); Idorsia (consultant); Pearson (author, royalty recipient); Equip Health Inc. (clinical advisory board). ML declares that over the past 36 months, he has received lecture honoraria from Lundbeck Pharmaceutical (no other equity ownership, profit-sharing agreements, royalties, or patent). The remaining authors declare no competing interests.
Open access
