Study participants

Participants were drawn from eight independent genome-wide association studies early-onset MPB samples: 23andMe (9,009 cases; 8,491 controls); Bonn (581 cases; 416 controls); CoLaus (622 cases; 655 controls); Nijmegen Biomedical Study (145 cases; 247 controls); QIMRB1 (216 cases; 1,162 controls); QIMRB2 (59 cases; 498 controls); THISEAS (52 cases; 150 controls); and TwinsUK (163 cases; 210 controls). A detailed description of the studies and their phenotype definitions is provided in Supplementary Note 1. All eight studies were approved by the respective institutional ethics review committees (specified in Supplementary Note 1), and written informed consent was obtained from all participants prior to inclusion.

Genome-wide association analyses

A summary of the study specific genotyping platforms, imputation methods and GWAS is provided in Supplementary Table 1. For each of the eight studies fulfillment of the following two criteria was required: (i) a minor allele frequency (MAF) of >1%; and (ii) a call rate of >98%, a variance ratio of ≥0.3 (MACH) or a proper info statistic of ≥0.4 (IMPUTE2). Imputed data were analysed using logistic regression and the dosage data options of either PLINK or SNPTEST.

Meta-analysis

The present meta-analysis was performed in accordance with the GWAS meta-analysis standards outlined in de Bakker et al.61 A fixed effects model was used to combine the logistic regression effect estimates of individual studies into a joint estimate, as implemented in METAL62. Results were cross-validated using the respective implementations in METAINTER and YAMAS. The Cochran’s Q statistics and the I2 measurement were used to test for cross-study heterogeneity. P-values for the test of heterogeneity were calculated according to Higgins et al.63 The meta-analysis included a total of N=8,004,650 SNPs that were available in the 23andMe cohort and (i) at least four additional studies for non-X-chromosomal SNPs, and (ii) three additional studies for X-chromosomal SNPs. The QQ-plot of the meta-analysis and a plot of the across-study homogeneity are shown in Supplementary Figs 2 and 3.

Identification of independent risk loci

SNP associations with a P-value of <10−6 (METAL) were extracted from the meta-analysis and filtered for a heterogeneity P-value (P Het ) of <0.01. SNPs separated by a distances of ≤100 kb on the autosomes and ≤500 kb on the X-chromosome were assigned to the same genomic region. Loci were considered to be significantly associated with MPB if at least one SNP within the defined region showed an association to MPB at P<5 × 10−8 (METAL) and (i) had a MAF of ≥0.05; or (ii) was supported by at least one SNP in linkage disequilibrium (LD; r2>0.5), with a MAF of ≥0.05 and an association with MPB at P<10−6 (METAL). MPB-risk loci were considered to be independent if the lead SNPs of the regions showed an LD of r2<0.2.

Test for polygenicity

For the meta-analysis overall, an inflation factor of λ=1.303 was observed. After removal of LD-SNPs surrounding the two major loci (chr20, chrX), the inflation decreased to λ=1.20. After exclusion of all 63 associated regions, the inflation decreased to λ=1.16. This indicates that a relevant proportion of the observed inflation is driven by genome-wide significant associations. Furthermore, it is reasonable to assume from this finding that the remaining inflation is driven primarily by true genetic association. This is consistent with previous reports of a polygenetic contribution to MPB, and with observations from meta-analyses of other complex genetic traits in which large numbers of common risk factors were identified64,65. To determine whether the inflation was attributable to polygenicity or to confounding factors, the LD score regression method66 was applied using LD Score v1.0.0. The analysis confirmed that most of the observed inflation is because of polygenicity (63.5%). The residual genomic inflation is λ GC =1.09, indicating that the observed association findings are not due to population stratification.

To estimate the proportion of potentially-false positive findings, association P-values (logistic regression) and effect directions for the 63 lead SNPs were compared between the 23andMe cohort (23andMe) and a meta-analysis of the remaining seven cohorts (MAwo23andMe) (Supplementary Table 8). For 62 of the 63 loci (98%) the analysis revealed a consistency of effect directions between 23andMe and MAwo23andMe. Of these loci, 76% (47/62) showed at least nominal significant association to MPB in the MAwo23andMe study (METAL), with 72% (34/47) also achieving significant association at P<0.05/63=7.9 × 10−4 (METAL P-value Bonferroni corrected for number of risk loci). Lastly, in contrast to the self-report based phenotyping used in the 23andMe cohort, phenotyping in the Bonn sample (BN) was based on clinical assessment. Thus for true positive association findings, stronger effect sizes would be expected on average in the BN sample. This was indeed the case: for 52 out of 61 lead SNPs available in the BN sample (85%) which would be highly unlikely if these loci included a relevant proportion of false positives. On the basis of these results, no correction was made for the residual inflation. The level of significance of the 63 MPB lead SNPs after double GC correction with the residual inflation factors from the LD score regression analyses is however indicated in Table 1.

Definition of credible sets of SNPs

For each risk locus, a credible set of SNPs was defined, as reported in Ripke et al.65 (Supplementary Table 9). In brief, here the connection between the deviance and the single-variant test statistics is used to calculate the probability that a specific variant is associated with MPB for the given set of data. These probabilities are used to determine sets of variants within the risk loci—that is, credible sets of SNPs—which contain the causal variant 99% of the time.

Estimation of heritability

The Bonn cohort was used to estimate the heritability of MPB, since individual genotypes were available. All imputed variants with an imputation info score of >0.6 were extracted, and the corresponding genotypes were called using a probability threshold of 0.9. The resulting dataset comprised 997 male individuals and 8,271,786 SNPs. Variants with a MAF of >0.01 and a genotype SNP-call rate of >95% were then filtered. Only SNPs with a P-value of >5% for the corresponding test of differences in missingness between cases and controls were retained. After outlier detection with EIGENSTRAT based on autosomal data and the exclusion of individuals based on a genetic similarity of >0.05, calculated by GCTA, we obtained the Genetic Relationship Matrix (GRM) for the autosomes and the X-chromosome with respect to 991 male individuals. With these two GRMs and 10 eigenvectors from the EIGENSTRAT step above, the heritability components were estimated using GCTA. The estimated heritability on the observed scale was 0.666 (0.234) for the autosomes and 0.277 (0.054) for the X-chromosome (standard errors in brackets). To obtain the estimated heritability on the liability scale, as described in Lee et al.67, a straightforward extension of this transformation was used to account for the extreme case and extreme control sampling. Table 2 in Hamilton68 was used to divide the population into four bins. On the basis of age structure and Hamilton-Norwood (HN) grade, the transformation was calculated under the assumption that controls were selected from the first bin and cases from the last, thus reflecting the extreme sampling scheme. The resulting factor was approximately 0.509. This generated an estimated heritability for autosomal chromosomes of 0.339 (0.119), and 0.141 (0.027) for the X-chromosome.

Estimation of explained phenotypic variance

The amount of phenotypic variance explained by the 63 genome-wide significant risk loci was estimated as the correlation coefficient r2 in a multivariate linear model, in which the significant index SNPs were predictors and the outcome was modelled as 1=case, 0=control.

Genotype-risk score analysis

Using the lead SNPs from the 63 genome-wide significant loci, a genotype-risk score was constructed based on the weighted number of susceptibility alleles in an independent replication sample from the Heinz–Nixdorff Recall (HNR) cohort (N=1,201). Individuals with HN grade>II (N=1,108) were defined as cases, and individuals with HN-grade I or II at age <65 years were defined as controls (N=93). The weights were established using the beta coefficient from the meta-analysis. The resulting genotype score was divided into four quartiles. The risk for MPB was then tested in each quartile using a logistic regression model and the lowest quartile as a reference. Two models were tested: (I) genotype-risk score; and (II) age-adjusted risk score. The results of this analysis are shown in (Table 2).

DEPICT analysis

DEPICT was used to identify plausible candidate genes at each of the 63 risk loci. As recommended in Pers et al.69, all independent SNPs with a P-value of <5 × 10−8 were included in the DEPICT analysis. Independent SNP sets were generated by retaining the most significant SNP from each set of SNPs with a pairwise LD of r2>0.1 and a physical distance of <500 kb. Pairwise LD coefficients were computed based on the imputing panel used in the eight GWASs (1000 Genomes Project Phase I CEU, GBR and TSI genotype data). The results of these analyses are shown in Supplementary Tables 4–6.

Identification of hair follicle eQTLs

DNA from peripheral blood and RNA from occipital scalp hair follicle samples were obtained from 125 volunteer healthy male donors of German descent (mean age 27.9 years). Genome-wide genotyping and imputation of blood DNA samples were performed using Illumina’s Human OmniExpress-12v1.0 bead array and IMPUTE2 (1000 Genomes, Phase I, June 2014). Whole-transcriptome profiling was performed on Illumina’s HT-12v4 bead arrays after amplification and biotinylation of the hair follicle-RNA using the TotalPrep-96 RNA Amplification Kit (Illumina, San Diego, CA, USA). Expression data were quantile normalized, and only probes with a detection P-value of <0.01 (Illumina GenomeStudio Software) in at least 5% of the samples were taken into account. The selected expression probes were subsequently filtered for: a unique alignment to the transcriptome; a perfect or good probe quality, as reported in the R package illuminaHumanv4.db; and mapping to an ENTREZ gene ID. After quality control and filtering, data for 14,687 expression probes and 6,593,881 SNPs were included in the eQTL analysis. Associations between gene expression levels and SNP genotypes in cis (distance between SNP and expression probe ≤1 Mb) was tested in MatrixEQTL using an additive linear regression model. Association tests were corrected for the top five principle components. All eQTL findings at a false discovery rate (FDR) of <0.001 were considered significant.

mRNA and miRNA expression analysis

Whole-transcriptome profiling of hair follicle-RNA samples (N=125) was performed as described above. MiRNA profiling of hair follicle samples (N=25) was performed on the AffymetrixGeneChip miRNA 4.0 (Affymetrix, Santa Clara, CA, USA) using a total of 250 ng of hair follicle miRNA. Poly(A) tailing and biotinylation were performed with the AffymetrixGeneChip Hybridization, Wash and Stain Kit, in accordance with the manufacturer’s instructions. Data were analysed using the Affymetrix Expression Console software (v.1.4) (Affymetrix). MiRNAs were defined as ‘present’ (N=1,169) or ‘absent’ (N=1,409) according to the implemented RMA (robust multichip average) and DABG (detected above background) methods.

Identification of miRNA target genes at MPB-risk loci

MiRNA target genes at MPB-risk loci (±500 kb from MPB lead SNP) were identified using the miRWalk2.0 algorithm (last accessed 24 March 2016). Only validated target genes and genes that were predicted by the miRWalk algorithm and at least three additional implemented databases were taken into account.

Overlap between eQTLs and MPB-risk variants

To identify regulatory effects at known and novel MPB-risk loci, the present meta-analysis data were compared with the eQTL data from: (i) the Blood eQTL Browser (http://genenetwork.nl/bloodeqtlbrowser; last accessed on 21 December 2015); (ii) the GTEx Browser (http://www.gtexportal.org/home); and (iii) an unpublished eQTL dataset from human hair follicle. MPB-risk variants were considered to coincide with an eQTL if the MPB lead SNP itself, or any SNP in r2>0.5, showed an eQTL effect. The complete list of overlapping eQTL findings is provided in Supplementary Table 2.

Enhancer enrichment analysis

Enhancers (including promoters) were defined by chromatin immunoprecipitation-sequencing (ChIP-seq) peaks of histone 3 lysine 27 acetylation (H3K27ac). A total of 140 H3K27ac data sets were downloaded from the roadmap epigenomics project (http://www.roadmapepigenomics.org), and reprocessed using DFilter. Only peaks with P<10−10 (DFilter) were retained. Four in house data sets from a human balding and a non-balding dermal papilla cell line with and without 10 nM of DHT treatment were included70. In these lines, a H3K27ac antibody (Abcam, Cat. ab4729; 5 μg antibody per 100 μg chromatin) was used for ChIP. The MPB credible SNPs were overlapped with the enhancer (H3K27ac) sets of each of the 144 cell lineages and tissues. For tissue (or cell lineage) A, the remaining 143 enhancer sets were pooled together to define a tissue agnostic superset. The superset was classified in feature categories based on: (i) distance to transcription start site; (ii) GC content; (iii) length of the enhancer site. To obtain for the enhancer set of tissue A an appropriate null distribution, we randomly drew for every enhancer of tissue A an enhancer from a matched feature category of the superset. The number of credible SNP overlaps across this matched enhancer set was computed, and a million permutations of feature matched enhancer sets were conducted and overlapped with the credible SNP set to define a distribution for tissue A. This procedure was conducted for each of the 144 tissue specific enhancer sets. For Fig. 2, enhancers of cell lineages/tissues that belonged to the same tissue type were merged, resulting in 23 grouped enhancer sets. The above framework was repeated for these grouped enhancer sets. The results of these analyses are provided in Supplementary Data 1.

Ingenuity pathway enrichment analysis

To test for an enrichment of MPB candidate genes in canonical pathways, the Ingenuity Pathway Analysis was used (Qiagen, Hilden, Germany). IPA considers 658 pathways, and calculates enrichment based on the right-tailed Fisher’s exact test. All genes within 500 kb of the MPB lead SNPs were included in the analysis. Only pathways with ≥3 annotated genes were taken into account. The complete list of nominally significant pathways is provided in Supplementary Table 7.

NHGRI GWAS catalogue

Previously reported GWAS associations that mapped to the associated regions and showed an r2≥0.3 and/or D′>0.8 (1000 Genomes Project, Phase 3) with the respective MPB index SNP were extracted from the NHGRI GWAS catalogue (http://www.ebi.ac.uk/gwas/; last accessed on 14 January 2016). This resulted in the identification of 124 overlapping associations, as shown in Supplementary Data 3.

Data availability

All data that support the findings of this study are available from the corresponding author upon reasonable request.