Sample collections

We conducted a GWAS across nine discovery data sets (the majority overlapping with Genetic Consortium for AN, as part of the Wellcome Trust Case Control Consortium 3 (WTCCC3/WTCCC3 samples), resulting in a total of 2158 cases and 15 485 ancestrally matched controls (Table 1 and Figure 1). All AN cases were female. AN diagnosis was made via semistructured or structured interview, or population assessment strategy using Diagnostic and Statistical Manual of Mental Disorders (DSM)-IV criteria for AN. The amenorrhea criterion was not applied, as this has been shown not to be diagnostically relevant9 and has since been dropped from DSM-5.10, 11 All cases met criteria for lifetime AN.

Table 1 Final numbers of cases and controls after QC Full size table

Figure 1 Geographical distribution of samples across Europe. (a) Distribution of cases across Europe; 375 USA cases are not shown in this diagram. (b) Distribution of controls across Europe; 873 USA controls are not shown in this diagram. PowerPoint slide Full size image

Exclusion criteria included confounding medical diagnoses, for example, psychotic conditions, developmental delay or medical or neurological conditions causing weight loss.

Ancestry-matched controls were selected for each AN case set. Both male and female controls were used (Table 1). These were obtained either from existing collaborations, or through genotyping repository (dbGaP) access. Each site obtained ethical approval from the local ethics committee, and all participants provided written informed consent in accordance with the Declaration of Helsinki.

Population prevalence of AN in these populations ranged from 0.4 to 3% (refs 12, 13, 14, 15, 16, 17, 18; Table 1).

Genotyping

Cases were genotyped on either the ‘Infinium HumanCoreExome-12 BeadChip Kit (Illumina, San Diego, CA, USA),19 or the ‘Infinium HumanCoreExome-24 BeadChip Kit (Illumina),20 at the Wellcome Trust Sanger Institute. Where possible, controls were selected from existing studies with matching genotyping platforms to cases. Three control cohorts had been genotyped on the ‘Infinium HumanExome-12 BeadChip Kit’ (Table 1). To ameliorate potential confounding due to chip effects,21 chip-type quality control (QC) was carried out, and ~14 000 single-nucleotide polymorphisms (SNPs) removed.

Quality control

Genotypes were called using the GenCall22 and Zcall23 algorithms. At each of these genotype-calling stages, QC was performed for each population and for cases and controls separately (Supplementary Table 1). The final number of SNPs included in the analyses is given in Table 2.

Table 2 Final number of SNPs per population Full size table

Controlling for population stratification

In order to account for population stratification, a principal components analysis was carried out for each cohort separately using the smartpca software.24

Population outliers were identified by merging each population with central European 1000 Genomes data.25

Variance explained by each PC was plotted for each population. In order to be both conservative and consistent across populations, the first 10 principal components were included as covariates in the association testing.

Association testing

Unbalanced case–control ratios can lead to anticonservative P-value estimates.26 This study includes a number of unbalanced strata (Table 1). The likelihood ratio test has been shown to have low type-I error rate across both balanced and unbalanced cohorts,26 and was chosen as the association test for this study.

A lower cutoff of minor allele count of 5 and MAF of 0.1% was used. Association testing was performed for each cohort separately using SNPtest.27 In the cohorts with mixed sex controls (all except Italy and Norway), sex was also included as a covariate.

The standard genome-wide significance threshold of P⩽5 × 10−8 was applied.

Meta-analysis

Summary statistics across cohort were meta-analyzed using an inverse variance-based test in METAL.28 In order to test the heterogeneity of the results, Cochran’s Q and the I2 statistic were computed.

Assigning variants to genes

Variants identified associated at P⩽1 × 10−4 were assigned to genes using Ensembl (release 83; Ensembl Genome Browser).29, 30 For each variant, all predicted consequences (for example, missense, non-synonymous, and so on) and associated gene transcripts were downloaded and compared. Each variant was associated with only one predicted consequence and one Ensembl gene ID (Ensembl Genome Browser).29

Cluster plot checking

Cluster plots were created for all SNPs reaching P⩽1 × 10−4 in any analysis (cohort-specific or meta-analysis) using ScatterShot.31 SNPs were visually inspected for each cohort, and for cases and controls separately. In instances where multiple cohorts were merged (for example, UK cases), cluster plots were checked separately for each original cohort.

Burden testing

The potential aggregation of rare variants in cases compared with controls was investigated using a gene-based approach. Burden tests were carried out using the Zeggini–Morris burden test32 as implemented in rvtests (Rvtests - Genome Analysis Wiki).

All SNPs with MAF between 0.1 and 5% were included; similar to the single-point analysis, a lower bound of minor allele count=5 was used. A list of genes and locations was obtained from the UCSC genome browser (Table Browser: www.genome.ucsc.edu). All genes with at least two qualifying variants in at least two populations were used, resulting in a total of 9083 genes.

Burden tests were carried out for each population individually, and the results meta-analyzed using Stouffer’s method, weighted according to effective sample size.33

The genome-wide significance threshold for burden testing is computed in a similar manner to that for single-point analysis, using Bonferroni correction for the number of genes tested. This results in a genome-wide significance threshold of 5.5 × 10−6.

Pathway analysis

One of the key motivations of studying complex psychiatric disorders such as AN is the desire to unearth biological pathways underlying disease development. Pathway analysis was performed using summary statistics from the meta-analysis for the full data set.

Four pathway databases were used: the Kyoto Encyclopedia of Genes and Genomes (KEGG),34, 35 the Reactome pathway database (REACTOME),36 PANTHER pathway (PANTHER)37, 38 and the Gene Ontology database (GO).39, 40 These were curated to remove redundancy, resulting in a total set of 1836 pathways.

The analysis was run once on a merged set of 235 KEGG,34, 35 REACTOME36 and PANTHER37, 38 pathways, and once for the 1601 GO pathways.39, 40

Pathway analysis was carried out using MAGMA.41 MAGMA was selected for its ability to deal robustly with linkage disequilibrium (LD) between markers, correct for gene length and deal accurately with rare variants. To our knowledge, MAGMA was first used to annotate SNPs to genes. This analysis was repeated twice. In the first analysis, variants were assigned only to the gene they were in, resulting in 68.73% of the variants being assigned to 13 400 genes. In the second analysis, variants were assigned allowing a 20 kb window in both directions from the gene. This procedure included 75.44% of variants across 18 118 genes.

SNP P-values were used to create gene scores. The European panel of the 1000 Genomes project was used as a reference set to estimate LD between SNPs. The analysis also requires the sample size of the study to be specified; because of the unbalanced nature of the study, the effective sample sizes were given here.

Gene P-values were calculated using MAGMA.41 The top 10% of SNPs per gene were used. Significance was defined using a false discovery rate of 5%.42

There is a risk when assigning SNPs to genes using MAGMA that some highly associated SNP might be assigned to multiple overlapping genes, and thus distort pathway results. SNP–gene assignments were checked for all pathways that reached false discovery rate-corrected significance. No instances of SNPs being assigned to multiple genes were found across these pathways.

Replication

SNPs reaching P<1 × 10−4 in the discovery stage were prioritized for replication. In total, 16 SNPs were selected.

Replication was carried out using two data sets: one existing in silico data set and one set for de novo genotyping. The in silico data set came from an existing GWAS of AN,7 genotyped on the Illumina HumanHap610 platform. This data set included 1033 cases and 3733 controls. All cases included in this study were female. Controls were both male and female. The de novo replication cohort consisted of 266 self-volunteered female UK cases, collected through the charity Charlotte’s Helix (www.charlotteshelix.net). All participants were adults and had been diagnosed with AN by their clinician. In addition, all participants completed an online questionnaire based on the Structured Clinical Interview43 for the Diagnostic and Statistical Manual of Mental Disorders-IV Module H. The Structured Clinical Interview has been used extensively in epidemiological investigations. The Structured Clinical Interview eating disorder module was modified to capture information on lifetime history of eating disorders including AN, and includes questions on body mass index, age of onset, and experience of eating disorders. DNA from the saliva samples was extracted using standard protocols and was quantified using pico-green. Samples were genotyped on the Infinium HumanExome 12 Beadchip, genotypes were called using GenCall and Zcall algorithms and stringent QC was performed pre- and post-call. In all, 1500 ancestry-matched controls (55% female) were obtained from the UK Household Longitudinal Study.

De novo genotyping was performed using the iPLEX Assay and the MassARRAY System (Agena Bioscience, San Diego, CA, USA) (formerly Sequenom). Sample and SNP QC were carried out within each replication data set, using an 80% sample call rate and a 90% SNP call rate threshold, and a Hardy–Weinberg equilibrium threshold of 10−4. Five samples and one SNP were removed using these criteria.

Post-QC, 15 SNPs and 261 de novo cases remained. The de novo replication analysis therefore included 15 SNPs, 261 cases and 1500 controls. Genotypes for 12/16 SNPs were available in the in silico replication cohort, across 1033 in silico cases and 3733 controls.

Expression analysis

Gene expression data were obtained from the Genotype-Tissue Expression (GTex project) web portal, data release version 6 (dbGap Accession phs000424.v6.p1).44, 45, 46

Power

The sample sizes used in this study are small in the context of other psychiatric phenotypes. Power to identify genome-wide significant signals was calculated using Quanto.47, 48 This study is adequately powered to detect low-frequency alleles with large effect sizes and common alleles with substantial effect sizes (80% power to detect common alleles with odds ratio (OR)>1.5; low-frequency alleles with OR>2, Supplementary Figure 1).

Data availability

Genotypes of European cases included in this study are publicly available through the European Genome-Phenome Archive (EGA), under accession number EGAS00001000913, data set EGAD00010001043, with the exception of German and Dutch genotypes. Genotypes for cases from the United States of America may be obtained through dbGaP. Summary statistics are available for download from the PGC website (https://www.med.unc.edu/pgc/results-and-downloads).