|Home | About | Journals | Submit | Contact Us | Français|
Anorexia nervosa and bulimia nervosa are common and severe eating disorders (EDs) of unknown etiology. Although genetic factors have been implicated in the psychopathology of EDs, a clear biological pathway has not been delineated. DNA from two large families affected by EDs was collected, and mutations segregating with illness were identified by whole-genome sequencing following linkage mapping or by whole-exome sequencing. In the first family, analysis of twenty members across three generations identified a rare missense mutation in the estrogen-related receptor α (ESRRA) gene that segregated with illness. In the second family, analysis of eight members across four generations identified a missense mutation in the histone deacetylase 4 (HDAC4) gene that segregated with illness. ESRRA and HDAC4 were determined to interact both in vitro in HeLa cells and in vivo in mouse cortex. Transcriptional analysis revealed that HDAC4 potently represses the expression of known ESRRA-induced target genes. Biochemical analysis of candidate mutations revealed that the identified ESRRA mutation decreased its transcriptional activity, while the HDAC4 mutation increased transcriptional repression of ESRRA. Our findings suggest that mutations that result in decreased ESRRA activity increase the risk of developing EDs.
Anorexia nervosa (AN) and bulimia nervosa (BN) are severe, debilitating eating disorders (EDs) characterized by marked disturbances in body image and abnormal eating patterns. In a meta-analysis of studies from Europe and the United States, AN has a 1-year prevalence of 370 cases per 100,000 young females, while BN affects 1,000 per 100,000 females each year (1). Although current treatment strategies focused on refeeding and family therapy do benefit patients, EDs still have one of the highest rates of mortality of all mental illnesses, with an annual crude mortality rate (the total number of deaths in the study population over 1 year) of 0.51% for AN and 0.17% for BN, with men being affected at approximately one-tenth of this rate (1). EDs are thought to occur as a result of a complex interaction between genetic predisposition and environmental risk factors. While genetic factors are estimated to contribute 50%–80% of the risk of developing an ED (2), to date, several studies using both genome-wide analysis (3, 4) and candidate gene (5) approaches have failed to identify specific genes that predispose to the development of an ED. One alternative approach to large-scale linkage and association studies is to characterize rare single-gene mutations in large families severely affected by mental illness.
Here, we have implemented this approach to identify signaling pathways that contribute to the development of EDs through the analysis of two large families with multiple members suffering from AN or BN. We report that mutations affecting the functional interaction of the transcription factor estrogen-related receptor α (ESRRA) and the transcriptional repressor histone deacetylase 4 (HDAC4) are associated with the development of EDs. ESRRA is an orphan nuclear receptor with sequence homology to estrogen receptors α and β and does not require estrogen binding for transcriptional activity (6). ESRRA has a demonstrated role in energy balance and metabolism (7), is upregulated by exercise (8) and calorie restriction (9) in peripheral tissues, and is a transcriptional target of the estrogen receptor (10). By contrast, HDAC4 is a transcriptional repressor that has been implicated in numerous processes relevant to EDs, including locomotor activity, body weight homeostasis, and neuronal plasticity (11–13).
Genetic analysis of twenty members of the AN-1 family (Figure (Figure1A)1A) was conducted on all affected and unaffected family members using a combination of microarray linkage and whole-genome sequencing. Linkage analysis identified a region on chromosome 11 from 44.1 to 64.3 cM with a lod score greater than 3.5 (Figure (Figure1B).1B). No other region of the genome showed linkage approaching this level of significance. Therefore, we performed whole-genome sequencing on DNA from two affected family members (AN 10 and 13) in order to identify mutations in the linkage region with a frequency lower than 1.0% in the population shared by both sisters. A missense mutation of G to A resulting in an arginine-to-glutamine substitution at position 188 of the gene encoding ESRRA was the only mutation predicted to affect the coding sequence identified in this region of chromosome 11 (Table (Table1).1). Subsequent genotyping of all family members for the G-to-A mutation showed that all ten affected family members inherited this mutation from AN 16, while nine of ten unaffected family members proved to be homozygous for the reference G allele at this position (Table (Table22).
Position 188 of ESRRA occurs in the “hinge” region between the DNA-binding domain and the ligand-binding domain of ESRRA, within a dibasic motif that is conserved across species (Figure (Figure1C).1C). We next used two computer algorithms, PolyPhen-2 and SIFT, to predict whether the mutation in this position would be deleterious (14). The R188Q mutation is predicted to be tolerated by SIFT (0.192), but possibly damaging by PolyPhen-2 (0.838). Therefore, the R188Q mutation in ESRRA was identified as the most likely candidate mutation in the AN-1 family because of the segregation with illness, the known functions of the gene, and a possible functional effect of the mutation.
Eight members of the AN-2 family (Figure (Figure2A)2A) were analyzed in parallel with the AN-1 family using whole-exome sequencing of five affected family members (AN 74, 80, 81, 82, and 83). All mutations that segregated with a phenotype and had an allele frequency of less than 1.0% in the 1000 Genomes database were identified, and six mutations on the distal end of chromosome 2 were found to cosegregate with affected status (Table (Table3).3). This observation suggests that the affected family members inherited this segment of chromosome 2. Subsequent linkage analysis assuming a dominant mode of inheritance using the six variants yielded a lod score of 2.4 at the region.
Of the six mutations, one was silent (BOK), and three occurred in the 3′ untranslated region (UTR) (ANO7, Septin2, and ING5). Of the two missense mutations that changed the amino acid sequence, the mutation in ATG4B was conserved (valine to isoleucine) and was not predicted to affect protein function (PolyPhen-2 = 0.02 and SIFT = 0.4). By contrast, a C-to-T mutation resulting in a change of alanine to threonine at position 786 of histone deacetylase 4 (HDAC4) was predicted to be deleterious to protein function (PolyPhen-2 = 0.965 and SIFT = 0.01). In addition, the alanine at position 786 of HDAC4 occurs near the end of α helix 8 of the deacetylase domain (15). This residue is completely conserved among class IIa HDAC members as well as among homologs across species, including zebrafish and fruit flies (Figure (Figure2B).2B). All six affected family members inherited this mutation, while both unaffected family members were homozygous for the reference C allele at this position (Table (Table4).4). Based on this combined approach, we assigned HDAC4 as the top candidate to emerge from analysis of this second family.
Several lines of evidence suggest that ESRRA and HDAC4 may have a functional interaction. ESRRA is a transcription factor implicated in mitochondrial biogenesis and energy partitioning (6), while HDAC4 is primarily recognized as a repressor of multiple transcription factors including MEF, RUNX, and STAT1 (16). Importantly, HDAC4 binds to and represses the ligand-independent transcriptional activity of the structurally related estrogen receptor α (17). Also, ESRRA activity is repressed by nuclear receptor corepressor 1 (NCOR1), which forms a functional complex with HDAC4 (18, 19). Finally, ESRRA and HDAC4 activity have both been implicated in several common physiologic processes, including the conversion of glycolytic to oxidative muscle fibers (11, 18), bone formation (20, 21), and cardiac function (22, 23). Given our findings, we hypothesized that ESRRA and HDAC4 may also function in a common biological pathway related to the development of EDs.
To test our hypothesis, we first expressed wild-type versions of the genes in HeLa cells and performed immunoprecipitation studies using antibodies against ESRRA and HDAC4 to probe for potential interaction. Using an anti-ESRRA antibody to purify all proteins bound to ESRRA, followed by Western blot detection of HDAC4, we demonstrated a specific interaction between the two proteins (Figure (Figure3A).3A). In order to confirm the specificity of the interaction, we then reversed the direction of immunoprecipitation and demonstrated by Western blotting that ESRRA is in the group of proteins immunoprecipitated with an HDAC4 antibody (Figure (Figure3A).3A). Next, we used the same approach to test whether the mutations we identified in each family affected this ESRRA-HDAC4 interaction. Each point mutation was introduced into wild-type plasmids, and pairwise combinations were expressed in HeLa cells. In all cases, we detected that the ESRRA-HDAC4 interaction was unaffected by the mutations (Figure (Figure3A).3A). After demonstrating that the two proteins interact in an in vitro overexpression system, we tested whether the interaction occurs under more physiologically relevant conditions by performing immunoprecipitation with anti-ESRRA antibodies from lysate obtained from mouse cortex, followed by detection of HDAC4 binding by Western blotting (Figure (Figure3B).3B). (The full, uncut gels are shown in the Supplemental Material.)
HDAC4 is a repressor of multiple transcription factors (16) including estrogen receptor α (17), which is structurally related to ESRRA. Therefore, we next tested our hypothesis that HDAC4 might affect the transcriptional activity of ESRRA using an in vitro system to measure ESRRA-induced changes in mRNA levels of known target genes. We first transfected ESRRA and its requisite coactivator peroxisome proliferator–activated receptor γ coactivator 1 α (PGC-1α) (24) into HeLa cells, with increasing concentrations of ESRRA and measured the relative abundance of mRNAs for α-ketoglutarate dehydrogenase (OGDH), aspartate aminotransferase cytoplasmic (GOT1), and medium-chain acyl-coenzyme A dehydrogenase (MCAD), genes known to be targets of the PGC-1α/ESRRA pathway (25, 26). As expected, transfection of PGC-1α-ESRRA increased the expression of all three genes (Figure (Figure4A).4A). We next determined whether HDAC4 expression affects ESRRA activity. While holding the concentration of PGC-1α-ESRRA constant, increasing amounts of HDAC4 were transfected, and the mRNA levels of the same genes were measured. Increasing HDAC4 levels induced a dose-dependent suppression of PGC-1α-ESRRA transcriptional activity, indicating that HDAC4 represses ESRRA activity (Figure (Figure44B).
Finally, we tested whether the mutations we identified in the AN-1 and AN-2 families carried a functional consequence for ESRRA-HDAC4 activity by repeating transcriptional activity measurements using mutant forms of the protein. The ESRRA-R188Q mutation significantly reduced the induction of target genes compared with wild-type ESRRA, indicating that the mutation decreases transcriptional activity of the protein (Figure (Figure4C).4C). We then measured the transcriptional repression activity of HDAC4-A786T in the PGC-1α-ESRRA assay. Compared with wild-type HDAC4, the A768T mutation significantly increased the transcriptional repression of ESRRA target genes (Figure (Figure4D).4D). Hence, both mutations decrease the expression of ESRRA target genes.
In the present study, we analyzed two large families with multiple family members affected by EDs. Such family studies are often complicated by environmental and psychological factors that encourage disordered eating patterns, complicating the assignment of affected and unaffected status. In particular, many family members commented during diagnostic interviews that they felt pressured by other family members with EDs to be thin. However, several features of these families suggest genetic involvement, including the large percentage of affected members, young age of onset, extremely low body weight, extended course of illness, and, in one family, multiple affected males. In order to minimize the miscategorization of participants, only a full diagnosis of AN or BN was included as “affected,” while the less severe diagnoses, such as “eating disorder, not otherwise specified,” were included as “subsyndromal.” Even by these stricter criteria, ten affected and ten unaffected family members were identified from the AN-1 family, and six affected and two unaffected members were identified in the AN-2 family.
While both families are severely affected by EDs, several interesting differences were noted. The AN-1 family has a high comorbidity of obsessive-compulsive disorder (OCD), with eight of ten affected family members also having a comorbid Axis I diagnosis of OCD, featuring obsessiveness, anxiety, perfectionism, and calorie restriction (27). By contrast, the AN-2 family has significantly lower rates of comorbid anxiety disorders and OCD, but higher rates of BN. Future studies will need to determine whether these phenotypic differences are related to the different biological effects of the ESRRA and HDAC4 mutations identified in the respective families.
Several important questions arise from the current study. First, it will be important to determine whether these mutations are found in sporadic cases of EDs. Both mutations are found in less than 1% of the population, indicating that a sample size of several thousands of cases will probably be required to make this determination. It will also be important to determine whether genetic variations in other genes in this pathway are more common in patients with an ED. One association study evaluating 182 candidate genes did not include either ESRRA or HDAC4, but did include PGC-1α, the coactivator of ESRRA (5). While, overall, no significant differences were reported, the rs11934028 SNP in PGC-1α was one of the top five SNPs identified in patients with AN without binge eating, with a trend P value of 0.0037. In addition to the numerous targets of ESRRA, several gene products have been reported to affect ESRRA activity, including RaRF (28), RIP140 (29), miR-137 (29), DPF2 (30), and PROX1 (31). Additionally, two other members of the estrogen receptor superfamily, estrogen-related receptor β and estrogen-related receptor γ, have similar DNA-binding motifs for ESRRA and have also been implicated in metabolic processes, suggesting a potential redundancy with certain aspects of ESRRA function (6, 32). Therefore, analysis of the estrogen-related receptors and their signaling pathways may also yield insights into the genetic basis of EDs.
While HDAC4 is expressed in the brain (33) and has a known role in synaptic plasticity (12, 34), much less is known about the function of ESRRA in the brain. Several lines of evidence, however, support our finding that decreases in ESRRA activity may be relevant to neuronal dysfunction in ED patients. First, ESRRA is highly expressed in the brain, including cortical regions implicated in EDs (35), such as the anterior cingulate cortex and the insular cortex (36). Additionally, ESRRA has a well-established role in mitochondrial biogenesis, and several papers have identified a role for mitochondrial functioning in neuronal plasticity. For instance, mitochondria are abundant in the synapse and affect several key processes, including neurotransmitter synthesis and release, calcium buffering, and ATP synthesis (37, 38), and the expression of mitochondrial genes is induced by synaptic activity in rodents (39). Furthermore, disruption of metabolic pathways in the mitochondria has been identified in synaptic vesicles prepared from the anterior cingulate cortex of mice with high anxiety (40), a trait frequently observed in ED patients (41). Finally, ESRRA induces the expression of both monoamine oxidase-A and -B, suggesting a potential role in the metabolism of monoamine neurotransmitters such as serotonin and dopamine (42).
Another important question is the role of estrogen signaling in HDAC4-ESRRA activity. Given the female predominance of EDs, it is possible that estrogen signaling may mediate this risk by altering ESRRA-HDAC4 activity. Both ESRRA and HDAC4 have predicted estrogen receptor α (ERα) binding sites in their promoters, (43) and estrogen has previously been shown to affect ESRRA expression (10). Conversely, ESRRA and ERα bind to similar DNA elements, with ERα binding as a dimer to the inverted repeat sequence of AGGTCA, while ESRRA binds as a monomer to the consensus extended half-site TnAAGGTCA (6). While initial studies suggested crosstalk between ESRRA and ERα at the level of DNA binding (6), these findings have yet to be supported by genome-wide approaches to identifying common binding targets of ERα and ESRRA. One interesting clue may come from the families themselves, in which the AN-1 family has two male members affected, while males in the AN-2 family appear to be unaffected carriers. Future studies will need to fully define the interaction of these two systems.
Clinicians at the University of Texas (UT) Southwestern Medical Center identified patients with multiple family members affected with EDs. Family members that consented to participation in the study had blood drawn for the collection of DNA and received a diagnostic interview (Structured Clinical Interview for DSM-IV-TR Axis I Disorders) administered by a trained research team member. Consensus diagnosis was obtained by agreement of five clinicians specializing in EDs, including one child psychiatrist, two adult psychiatrists, one PhD psychologist, and one Master’s psychologist. Affected status was assigned to anyone with a lifetime history of AN or BN for two reasons. First, patients with AN or BN will frequently convert to the other diagnosis within their lifetime, suggesting a shared biological risk (44). Second, BN is more frequently found in relatives of patients with AN than in the general population, indicating that both disorders may share a common genetic risk (45). Individuals with no history of ED symptoms were designated as unaffected. Individuals with ED symptoms that did not meet the criteria for AN or BN, such as night eating or self-induced vomiting without a history of binging or being underweight, were designated as subsyndromal, since they could not be clearly identified as either affected or unaffected and were not included in linkage determinations.
Twenty-three family members agreed to participate in the interview, ten of whom were found to be affected, ten unaffected, and no consensus could be reached on three individuals. The affected and unaffected family members were genotyped using the Illumina HumanLinkage-24 array at the Microarray Core Facility of UT Southwestern Medical Center. A genome-wide, model-based multipoint linkage analysis was performed using the MORGAN package (46) under the assumption of a dominant genetic model with penetrance of 0.95, a sporadic rate of 0.001, and a disease–predisposing allele frequency of 0.001. Whole-genome sequencing was performed on two affected individuals (10 and 13) by Complete Genomics Inc. using a sequencing-by-ligation method (47). The linkage interval was effectively interrogated with average unique scores of 67.5 and 63.5, respectively, with greater than 98% of the regions in the interval covered by greater than 20.
For analysis of the AN-2 family, five samples (AN 74, 80, 81, 82, 83) from affected family members were subjected to whole-exome sequencing performed by the UT Southwestern McDermott Sequencing Core facility. Because of the smaller number of affected family members, we selected whole-exome sequencing because it allowed us to perform a comprehensive analysis at a substantially reduced cost. Indexed DNA libraries were constructed by the sequencing core using the TruSeq DNA Sample Preparation Kit, according to the manufacturer’s protocol. Target enrichment was then performed by pooling four indexed samples and using the human TruSeq Exome Enrichment Kit for each pool, also according to the manufacturer’s protocol. This kit targets 62 Mb of exomic sequences in the human genome, including 5′UTR, 3′UTR, microRNA, and other noncoding RNA regions. The resulting libraries were sequenced using the HiSeq 2000, producing 100-bp paired-end reads.
Reads were demultiplexed using CASAVA, version 1.8.2 (Illumina), allowing up to one mismatch in the index sequence, and clusters not passing filter were excluded from the resulting FASTQ files (Illumina). The reads were then quality checked using FastQC, version 0.8.0 (Babraham Bioinformatics). Each sample’s reads were aligned to the hg19 build of the human genome reference sequence using Burrows-Wheeler Aligner (BWA), version 0.6.2, and duplicates were removed using Picard MarkDuplicates, version 1.41. Sample level realignment around indels and base quality recalibration were performed using GATK methods, according to the Broad Institute’s best practices guidelines. Coverage was also calculated using GATK, and an average coverage of 20× was achieved across targeted regions. SNPs and indels in the targeted regions were called using GATK’s UnifiedGenotyper, followed by hard filtering (indels) and VariantRecalibration (SNPs). During recalibration, data from the OMNI and HapMap SNP chips were used as training sets. Those variants passing filters were annotated using successive ANNOVAR runs (2012-05-25), wherein dbSNP IDs, 1000 Genomes frequencies, polyPhen, and SIFT scores were added. ANNOVAR also provided the genomic location, and when variants were found to be exonic, it annotated the gene name, transcript ID, functional change, nucleotide, and amino acid change of each variant.
Coverage of 20× was achieved for greater than 84% of all baits in the five samples. The annotated exome data were then analyzed based on the following assumptions: (a) autosomal dominant transmission; (b) 90% penetrance; and (c) less than 1% frequency in the 1000 Genomes database. Variations that segregated with illness were confirmed by Sanger sequencing. Six variants on chromosome 2 were found to segregate with affected status. This linkage interval was effectively interrogated by exome sequencing, with greater than 95% of all coding regions in the interval covered by five or more sequencing reads in all samples.
Human ESRRA and mouse PGC-1α cDNAs were provided by Carol Mendelson of UT Southwestern Medical Center (48). Human HDAC4 cDNA was provided by Eric Nestler of Mount Sinai School of Medicine (New York, New York, USA) (49). Wild-type FLAG-tagged ESRRA and HA-tagged HDAC4 were subcloned into pcDNA3.3 (K8300-01; Invitrogen), and point mutations corresponding to the human mutations were introduced into the sequences using the GeneArt Site-Directed Mutagenesis System (A13282; Invitrogen). HeLa cells (ATCC) were plated at approximately 25% confluency in either 6-cm (gene expression studies) or 10-cm (protein expression studies) dishes 24 hours prior to transfection and transfected with the stated amount of DNA using X-tremeGENE HP DNA Transfection Reagent (06366236001; Roche). Empty vector (EV) was used to equalize all amounts of transfected DNA. For protein studies, cells were harvested, washed in ice-cold PBS, resuspended in NP-40 buffer (50 mM Tris-HCL, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, and 1% NP-40 at pH 7.4), passed through a 25-gauge needle three times and then subjected to three rapid freeze-thaw cycles in liquid nitrogen. Mouse frontal cortex and hippocampus were rapidly dissected from 12-week-old female C57BL6 mice, rinsed in ice-cold PBS, and placed in NP-40 buffer. Tissue was crudely separated using a tissue pestle (12-141-364; Fischer Scientific) and then sonicated by three brief pulses with an ultrasonic processor (40%; GEX 30PB; Cole-Parmer). Lysates were cleared by ultracentrifugation, protein concentration was measured by the BCA method (23228; Thermo Fisher Scientific), and samples were rotated overnight at 4°C with the following antibodies: anti-FLAG M2 magnetic beads (M8823; Sigma-Aldrich); anti-HA magnetic beads (88836; Thermo Scientific); and anti-ESRRA (04-1134; Millipore) conjugated with protein G–coupled Dynabeads (100.07D; Invitrogen). Beads were separated using a Dyanl System magnetic bar (Invitrogen) and washed four times in wash buffer. Protein was eluted by boiling in Laemmli Sample Buffer (161-0737; Bio-Rad) supplemented with 5% β-mercaptoethanol and separated by SDS-PAGE, probed using the specified antibodies anti-HDAC4 (ab79521; Abcam) or anti-ESRRA (04-1134; Millipore), and visualized by enhanced chemiluminescence.
mRNA was purified using the RNeasy Plus Mini Kit (74134; QIAGEN) and measured using a NanoDrop 2000 (Thermo Fisher Scientific). Complementary DNA (cDNA) was synthesized using SuperScript VILO (11755-050; Invitrogen). Quantitative PCR was performed on an Applied Biosystems Model 7900HT with an annealing temperature of 50°C using 50 ng cDNA and iQ SYBR Green Supermix (170-8882; Bio-Rad) with 200 nM of the following primers: 36B4: CACTGGTCTAGGACCCGAGAAG/GGTGCCTCTGAAGATTTTCG; GOT1: GTGAGGAAGGTCGAACAGAAG/CCCTAAGAAGTCAGCTCCAATC; MCAD: GATGCCATCACCCTCGTGTAAC/AAGCCCTTTTCCCCTGAAG; and OGDH: GACCGGAAGGAGCAGAATAAG/CTTGGCCCTAGCAAGATGAA. Experiments were performed in triplicate.
A two-tailed Student’s t test (Prism 6; GraphPad Software) was used to compare differences in transcriptional activity levels of wild-type and mutant ESRRA and HDAC4 proteins. P < 0.05 was considered significant.
Human studies were reviewed and approved by the IRB of UT Southwestern Medical Center. Written informed consent was obtained from participants or their guardians. All mouse studies were approved by the IACUC of the University of Iowa.
The authors would like to thank Helen Hobbs for providing access to the AN-1 family and for her assistance with experimental procedures. The authors would also like to thank Alex Bassuk and Elizabeth Campbell for their assistance analyzing the whole-exome data, and Barbara Gilbert and Marianne Nash for their assistance with collecting blood samples. The Disease Oriented Scholars Program and a High-Impact/High-Risk grant from UT Southwestern Medical Center supported this work.
Conflict of interest: The authors have declared that no conflict of interest exists.
Citation for this article: J Clin Invest. 2013;123(11):4706–4713. doi:10.1172/JCI71400.