Introduction

Frontal fibrosing alopecia (FFA) is a recently reported lichenoid and scarring inflammatory skin disorder associated with widespread cutaneous inflammation and irreversible hair loss, which occurs predominantly in women of postmenopausal age (Fig. 1)1,2. Since FFA was first identified by Kossard in 1994, there has been rapid increase in reported incidence culminating in intense clinical and public interest in the condition. FFA is often referred to as a dermatological epidemic with possible environmental trigger(s) implicated3. Nevertheless, the pathogenesis of FFA also includes a genetic component, as evidenced by frequent familial segregation4,5,6,7,8,9,10,11.

Fig. 1
figure 1

Clinical features of frontal fibrosing alopecia. Scalp with frontal hairline recession (a) involving the temporal areas bilaterally (b), as well as eyebrows (c). Histopathology (d) shows two hair follicles with focal interface changes, and a moderately dense perifollicular lymphoid cell infiltrate with perifollicular fibrosis, characteristic of FFA (×200)

Full size image

FFA is considered to be a clinical sub-variant of lichen planus, a more common inflammatory skin condition of unresolved aetiology, while also representing a variant of the prototypic primary lymphocytic cicatricial (or scarring) alopecia lichen planopilaris (LPP). A key molecular event in the pathology of scarring hair loss has been postulated to be the immune privilege collapse at the level of the immunologically shielded hair follicle bulge, which is home to epithelial hair follicle stem cells (eHFSC): T-cell mediated inflammatory presence culminates in stem cell apoptosis and irreversible alopecia12. Dissection of the genetic basis of FFA and its interplay with environmental risk factors, therefore, could provide insight into the molecular profile of lichenoid inflammation, scarring and mechanisms of immune privilege collapse. Furthermore, the identification of environmental triggers that interact with FFA genetic susceptibility loci could ultimately contribute to disease prevention by avoidance of exposures in genetically predisposed individuals (Supplementary Note 1).

To date there have been no systematic investigations into the molecular genetic basis of FFA or any other lichenoid inflammatory disorder. We hypothesise that common genetic variation contributes to FFA susceptibility and undertake a genome-wide association study and meta-analysis of two independent European cohorts of females with FFA and controls and investigate transcriptomic and metabolomic involvement in the disease.

Results

Genome-wide association study

We undertook a genome-wide association analysis across 8,405,903 common variants in a UK cohort of 844 FFA female cases and 3760 female controls. Inspection of the quantile-quantile plot indicated adequate control of confounding bias (λGC(MAF>0.05) = 1.03; Supplementary Table 2; Supplementary Figure 3). We observed genetic variants with genome-wide significant association (P < 5.0 × 10−8) with FFA at three genomic loci; 6p21.1, 2p22.2 and 15q26.1 (Table 1). We estimate the genome-wide SNP heritability for FFA as 46.66% (SE = 3.00%).

Table 1 Genome-wide significant loci for UK, Spain and meta-analysis
Full size table

In an attempt to replicate the observed associations at each of the three loci and identify additional FFA susceptibility loci, we performed a genome-wide association study across 7,964,651 common variants in our independent Spanish cohort comprising 172 affected females and 385 controls. We observed allelic associations with FFA at each of the three loci that had been implicated in FFA susceptibility in the UK cohort. The direction and magnitude of the effect of these associations was consistent between the UK and Spanish cohorts (Table 1). To identify further loci harbouring variation contributing to FFA risk we performed a statistical meta-analysis of the association summary statistics from the UK and Spanish cohorts. This revealed a single additional risk locus at 8q24.22, again with a consistent direction and magnitude of effect in both studies (Fig. 2 and Table 1) and a number of loci at which there is suggestive evidence of association (P < 5 × 10−5; Supplementary Table 3).

Fig. 2
figure 2

Manhattan plot showing the P values for the meta-analysis genome-wide association study. Each SNP was tested for association by logistic regression using an additive regression model; the interrupted line indicates the threshold for genome-wide significance (P = 5 × 10−8); the y axis has been collapsed for better illustration of all genomic signals; the continuous line represents the threshold for suggestive significance (P = 1 × 10−5); N = 5161 biologically independent subjects (Ncases = 1044 and Ncontrols = 4145)

Full size image

We sought to further investigate the allelic basis for the observed FFA association at 6p21.1, which is located within the MHC region. We undertook imputation of classical MHC Class I alleles and evaluated the association of each allele with FFA. The strongest evidence of association was observed for the HLA-B allele HLA-B*07:02 (PMeta = 9.44 × 10−117, OR = 5.22 (4.53–6.01); Supplementary Table 4), indicating that this is the most likely classical HLA allele to be underpinning the observed SNP associations in this region. Although full characterization of the role of HLA genes in FFA is complicated by the complex linkage disequilibrium structure across the region, a sequential conditional analysis with both SNPs and imputed HLA alleles indicates that there may be at least a further two independent HLA-alleles that contribute independently to disease risk (Supplementary Figure 3 and Supplementary Tables 4A and 4B).

To further investigate causal genes and alleles at the three remaining FFA susceptibility loci, we performed Bayesian fine-mapping of the association signals. This process identified a single putative causal variant with a posterior probability >0.5 of being the causal variant underlying the association signal at two of the three loci (Fig. 3 and Supplementary Table 5). At 2p21.2, rs1800440, is likely to be the causal variant underlying the association signal with a posterior probability of 0.98 (Fig. 4). The FFA protective allele is a missense allele (c.1358A>G p.Asn453Ser) in the CYP1B1 gene, which introduces a serine residue in the haem binding domain of the enzyme. In silico pathogenicity prediction tools predict that this allele has a deleterious effect on the function of the protein (SIFT = 0.015; CADD = 32)13,14, which is corroborated by published functional investigations of the p.Asn453Ser substitution15. At the 8q24.22 locus, rs760327 is the most likely causal allele (posterior probability = 0.68). The variant is located within intron 1 of the ST3GAL1 gene encoding the homonymous galactoside sialyltransferase enzyme, which has been studied in the context of T cell homeostasis16,17. At 15q.26.1, statistical fine mapping was unable to clearly resolve the causal variant at this locus (Supplementary Table 5), though the most likely causal variant rs34560261 (posterior probability = 0.4) is located within intron 1 of SEMA4B. Co-localization with skin eQTLs (Supplementary Figure 5) provides evidence that the same variant(s) underlying the observed FFA association at this locus are also associated with variation in the expression of SEMA4B in the skin (Pcoloc = 0.99), providing support that SEMA4B may be the causal gene at this locus.

Fig. 3
figure 3

Regional plots of lead signals at loci 15q26.1 and 8q24.22. a Regional plot at locus 15q26.1 (rs34560261); and b Regional plot at locus 8q24.22 (rs760327). The blue lines shows the fine-scale recombination rates (right y axis) estimated from individuals in the 1000 Genomes population; genes are highlighted with horizontal lines; the x axis shows the chromosomal position in Mb; variants within the 95% credible set are highlighted in red

Full size image
Fig. 4
figure 4

Regional plot at locus 2p22.2 and CYP1B1 enzyme structure. a Regional plot of lead signal at locus 2p22.2 demonstrating the lead causal variant (rs1800440); the blue lines shows the fine-scale recombination rates (right y axis) estimated from individuals in the 1000 Genomes population; genes are highlighted with horizontal lines; the x axis shows the chromosomal position in Mb; variants within the 95% credible set are highlighted in red (N = 1, posterior probability = 0.98); b CYP1B1 enzyme structure drawn via the SWISS-MODEL repository with the site affected by the residue change from Asparagine (N) to Serine (S) at position 453 magnified (https://swissmodel.expasy.org/repository)

Full size image

Plasma metabolomic analysis

At the 2p22.2 locus statistical fine mapping of the causal allele to a functional missense variant in CYP1B1 implicates variation in xenobiotic and endogenous hormone metabolism18,19,20,21,22 as a potential mechanism influencing disease susceptibility. To evaluate if there are systematic differences in metabolomic profiles between FFA cases and controls we compared levels of plasma metabolites in 52 treatment-naïve FFA cases and 35 ethnicity-, gender-, age- and BMI-matched healthy controls (Supplementary Figure 6). We did not observe differences in the distribution of individual metabolite levels between cases and controls of the magnitude detectable by this experiment following multiple-testing correction (Supplementary Table 6), nor did we observe any enrichment of xenobiotic or endogenous hormone metabolites in the extremes of the distribution of observed mean differences.

Transcriptomic analysis

To further investigate genes and biological pathways involved in FFA pathobiology we performed transcriptome sequencing in lesional scalp skin from seven treatment-naïve postmenopausal FFA cases and compared to transcriptome profiles from scalp skin in seven healthy matched controls. Differences in transcript abundance between cases and controls in these bulk tissue samples was observed for 117 genes with a greater representation of transcripts from 80 genes in affected tissue and 37 genes with reduced representation of their respective transcripts (Supplementary Tables 7–9; Supplementary Figure 7). Of these, only C2, within the MHC at 6p21.33, is located within 1 Mb of any of the FFA associated loci and only two of the 117 genes (CCL19 and EPSTI1) are located within 1 Mb of any variant with moderate evidence of association (P < 5 × 10−5) with FFA. Investigation of the enrichment of gene sets and pathways indicate that immune genes are over-represented amongst the differentially expressed genes (DEGs) (Supplementary Tables 8–12). Notably, four of the 10 most extreme DEGs are genes that have an established role in the interferon gamma (IFNγ) pathway.

Discussion

We have identified four genomic loci, at which genetic variation is robustly associated with the lichenoid inflammatory and scarring dermatosis FFA in two independent cohorts of European ancestry.

The strongest effect on FFA susceptibility is observed at 6p21.1 which is located within the MHC region. Through imputation of classical HLA alleles we implicate the Class I allele HLA-B*07:02 as conferring a five-fold increase in risk of FFA. The highly polymorphic HLA genes and their encoded proteins play a key role in self and non-self immune recognition and are known to determine susceptibility to numerous infectious and auto-immune disorders23. HLA-B*07:02 itself has previously been reported to be associated with HIV progression but has not been implicated in the susceptibility to human disease24. The hair follicle bulge region and the outer root sheath express low levels of HLA-A, HLA-B and HLA-C and these are key to rendering immune privilege25,26. HLA-B*07:02 may facilitate the process of hair follicular autoantigen presentation culminating in the auto-inflammatory lymphocytic destruction of the hair follicle bulge and its resident epithelial hair follicle stem cells. Investigation of differentially expressed genes between affected and unaffected scalp tissue further highlights the importance of genes encoding the components of innate and adaptive immunity and, notably, the IFNγ pathway, which is an important regulator of antigen presentation. Also relevant to a putative role of T cell dysfunction in FFA, the lead variant at the 8q24.22 locus is located in intron 1 of ST3GAL1, which encodes a membrane bound sialyltransferase. Changes in cell surface glycan structures have been implicated in human T cell lymphocyte activation and maturation27 and ST3GAL1 itself has been implicated in immunity by homeostatically controlling CD8+T cells16,17.

At 2p22.2 we observe strong evidence that the causal variant underlying the association at this locus is a missense variant in CYP1B1, a ubiquitously expressed gene encoding the Cytochrome P450 1B1 microsomal enzyme, also known as xenobiotic monooxygenase and aryl hydrocarbon hydroxylase. This enzyme contributes to the oxidative metabolism of oestradiol and oestrone to their corresponding hydroxylated catechol oestrogens28,29,30. Functional investigation of allelic variation in CYP1B1 has shown that the FFA protective p.Asn453Ser allele increases the rate of CYP1B1 degradation leading to reduced intracellular CYP1B1 levels15. Given the established role of CYP1B1 in sex hormone metabolism, alongside the female preponderance of FFA and its rapid and recent increase in incidence, it is plausible to speculate that an increase in exposure of females to a CYP1B1 substrate, whether endogenous or exogenous, may contribute to the development of FFA. The temporal relationship between the introduction of oral contraceptives in the 1960s and the appearance of FFA in the published literature in the 1990s should be fully explored with a well-designed gene-environment interaction study. Nevertheless, no striking differences in such substrates nor any other metabolites were identified in our metabolomic study although this may reflect the limited power of this initial investigation to observe more subtle disruption of this metabolic pathway. It should be noted that CYP1B1 has also been implicated in human immune cell regulation31,32 and the potential that FFA susceptibility at the 2p22.2 is mediated through immune pathways cannot be excluded.

In summary, in this exploration of the molecular genetics of FFA susceptibility we identify common alleles at four genomic loci that contribute to disease risk. The putative biological impact of this genetic variation indicates that the disease is a complex immuno-inflammatory trait underpinned by risk alleles in MHC Class I molecule-mediated antigen processing and T cell homeostasis and function. The insight into the pathobiology of FFA from the genetic susceptibility loci combined with the observation that there is an increase in transcripts encoding components of the IFNγ pathway in affected scalp tissue suggest that drugs such as JAK inhibitors, some of which have already proved effective in alopecia areata33 and trialled in lichen planopilaris34, may prove to be efficacious for FFA.

Methods

Clinical resource

Ethical approval was granted by the Northampton NRES Committee, UK (REC 15/EM/0273) and the study was conducted in accordance with the Declaration of Helsinki (https://www.wma.net/policies-post/wma-declaration-of-helsinki-ethical-principles-for-medical-research-involving-human-subjects/). We ascertained two independent cohorts of highly clinically consistent female cases of classic FFA, diagnosed by specialist dermatology clinics in the UK and Spain. All recruited cases were of European ancestry and diagnosed with FFA on the basis of the following clinical and histopathological features (recently proposed as diagnostic criteria)35: (1) cicatricial alopecic involvement of the frontal, temporal/parietal hair margin; (2) bilateral eyebrow loss; (3) clinical, trichoscopic (or histological) evidence of lichenoid perifollicular inflammatory presence; (4) facial or body hair loss; (5) absence of multifocal scalp involvement and other signs suggestive of classic LPP or its Graham-Little-Piccardi-Lasseur subvariant.

All research participants provided written informed consent for participation in the study. The individual depicted in Fig. 1 provided informed consent for publication of her clinical images.

Genotyping and genome-wide association analysis

For the UK cohort, genome-wide genotyping of cases was undertaken using Infinium OmniExpressExome BeadChip array (Illumina) and an unselected female control cohort from the English Longitudinal Study of Aging (ELSA) project (http://www.elsa-project.ac.uk), genotyped on the Infinium Omni2.5M BeadChip array (Illumina). We retained variants that were assayed with the same probe design on both genotyping arrays and excluded variants with a call rate of <99% or which deviated from Hardy–Weinberg Equilibrium (P < 10−4). Individuals with a call rate of <99% or extensive heterozygosity were also excluded. A subset of 46,789 variants in linkage equilibrium (r2 < 0.2 between each pair) was used to evaluate relatedness between individuals using the KING software package (KING; version 2.1.1). Individuals were thus removed from the study such that no two individuals had estimated relatedness closer than third degree (Kinship coefficient > 0.0442). Principal component analysis of the same set of 46,789 variants was performed and individuals outlying the main cluster (implying non-European ancestry) were also excluded from further analysis.

In the Spanish cohort, FFA case genotyping was performed using the Infinium OmniExpressExome BeadChip array (Illumina). Genotypic data for unaffected controls were obtained from 1061 individuals from the INfancia y Medio Ambiente (INMA) project (Valencia, Sabadell and Menorca, Spain http://www.proyectoinma.org) genotyped on the Omni1-Quad BeadChip (Illumina). Genotype calling, quality control and imputation were undertaken following the same protocol as described for the UK cohort across all variants that are assayed on both genotyping arrays with the same probe design.

Following quality control, genome-wide imputation was performed for both cohorts using the Michigan Imputation Server, with the Haplotype Reference Consortium (HRC) reference haplotype panel (www.haplotype-reference-consortium.org). All variants with an imputation info score <0.7 or a minor allele frequency of <0.005 were excluded from downstream analysis. This process of data generation and QC resulted in a combined total of 7,039,930 variants successfully genotyped or imputed in a combined total of 1016 cases and 4145 controls.

Genome-wide association testing was performed on all variants with MAF > 0.005 using a logistic Wald association test (EPACTS), including the first five principal components as covariates. Association testing was performed based on 844 affected females and 3760 female controls from the UK cohort and separately for 172 affected females and 385 female controls from the Spanish cohort. Association summary statistics were subsequently combined for 7,039,930 variants across the UK and Spanish cohorts via a standard error-weighted meta-analysis using METAL36. To evaluate potential confounding bias due to population stratification or residual cryptic relatedness we calculated the genomic inflation factor (λGC) for variants with MAF > 0.05 and the LD score regression intercept for each cohort and the combined meta-analysis37.

Causal variant identification and evaluation

For the UK and Spanish case-control cohorts, imputation of classical HLA alleles to two- and four-digit resolution was performed with the SNP2HLA tool, based on the genotypes of 1297 MHC region single nucleotide polymorphisms (SNPs) genotyped in both phases38. In the UK cohort, dosage-based association testing was performed in PLINK v1.9 for all 208 alleles that were well imputed (r2 > 0.9), using a logistic regression framework that included the same covariates as the full GWAS39. Replication of specific variants of interest was undertaken in the Spanish replication cohort in the same way, and standard-error weighted meta-analysis was performed using the meta package in R40. To test for multiple independent association signals, stepwise conditional analysis was performed: at each round of testing, the dosage of the HLA allele achieving the lowest discovery phase association p-value in the previous round was added to the list of covariates. This process was iterated until no allele achieved genome-wide significance (PMeta < 5 × 10−8). To ascertain which specific variants underlie the observed allelic associations, a similar stepwise analysis was performed using imputed HLA-region SNPs in place of HLA alleles. To verify that no additional genome-wide significant independent signals remained after the final iteration, we also tested for association of the HLA-region SNPs in the full GWAS dataset after conditioning on all independently-associated HLA alleles.

At 2p22.2, 8q24.22 and 15q26.1, fine-mapping was undertaken to identify putative causal variant(s) underlying the observed association signal41. For each locus, we constructed a credible set of variants considered most likely to be causal based on evidence for association as quantified by their Bayes factor42.

In order to explore the correlation between genetic variation and tissue expression, eQTL colocalization analysis was performed between the observed FFA association signals and sun-exposed skin cis-eQTL data from the Gene-Tissue Expression (GTEx) project database43. Candidate skin eQTLs were identified by looking into whether any variant within the FFA risk loci was associated with varied expression of nearby genes. Bayesian testing for colocalization between the FFA association signal and the skin eQTL signal was undertaken using a set of overlapping variants for the two datasets, employing the R package coloc tool44, with a defined prior probability of colocalization of P = 10−5.

Heritability estimation

FFA heritability explained by genome-wide SNPs (MAF > 0.01) was estimated using the genomic-relatedness-based restricted maximum-likelihood (GREML) approach, implemented in the software tool package GCTA45. Heritability estimates were expressed on the liability scale using an estimated prevalence of FFA of one in 5000.

Transcriptomic analysis

We performed transcriptome profiling with RNA-sequencing of scalp skin from seven cases of European ancestry and seven matched controls (Supplementary Table 1). All seven cases were clinically evaluated prior to obtaining skin biopsy from actively inflamed parietal scalp skin for histologic confirmation of FFA. Samples from cases were only subjected to downstream processing if they were confirmed to be actively inflamed upon histological evaluation. Macroscopically unremarkable parietal scalp skin was also harvested from healthy controls undergoing plastic facial surgery and all control tissue specimens were also examined microscopically and confirmed to be histologically normal (Supplementary Figure 1). Total RNA was isolated from each tissue sample using the RNeasy Plus Universal kit (Qiagen, Valencia, USA) as per the manufacturer’s protocol and instructions. Samples with RNA Integrity Number (RIN) < 8 were rejected from further processing.

Whole transcriptome RNAseq libraries were prepared using the Agilent SureSelect strand-specific RNA library preparation kit (Agilent, Santa Clara, USA) and multiplexed sequencing was performed on the HiSeq 2500 platform (Illumina, San Diego, USA).

Processing of the raw transcriptomic data files was conducted using an established analytical pipeline (Supplementary Figure 2). EdgeR software package in the (R-based) Bioconductor was utilized to undertake differential expression analysis46, following the trimmed mean of M-values (TMM) normalization method47. Genes with very low expression (defined as genes with counts per million (CPM) < 1 in at least seven samples) were discarded and not considered for further differential expression analysis. Transcript abundance was estimated and compared between the two groups with an exact (Robinson & Smyth) negative binomial (NB) test in EdgeR (Supplementary Figure 2). The P value distribution was obtained and Benjamini-Hochberg (BH) adjusted P values were estimated: genes with a false discovery rate FDR ≤ 0.05 and log fold change LFC ≥ 1 or ≤ −1 were considered significantly differentially expressed and extracted. Such differentially expressed genes were used to generate a heat map using the R package gplots2 via the START interface48.

Gene set enrichment and pathway analysis (GSEA) was performed as implemented in GAGE (R package) using C2 (pathways) and C5 (gene ontology) gene set collections from the Molecular Signatures Database (MSigDB)49,50. The p values for the gene set enrichment were calculated using default Stouffer test in GAGE and FDRs were generated using the BH procedure. Differentially expressed HLA genes were excluded from the analysis because of the challenges of accurately quantifying the expression levels of these highly polymorphic genes (due to the difficulty of correctly mapping divergent reads to a single reference genome).

Plasma metabolomic analysis

Fifty-two treatment-naïve post-menopausal female FFA cases of European ancestry (median age 64; mean BMI 24.7) and 35 matched controls (median age 58; mean BMI 24.5) were recruited. Peripheral venous blood was collected and centrifuged (at 1300 g for 15 min) to separate plasma, which was aliquoted and stored at −80 oC until required for further analysis. Metabolomic profiling of samples was undertaken by Metabolon (Durham, NC, USA) by subjecting plasma samples to methanol extraction prior to splitting into aliquots for analysis by ultra-high pressure liquid chromatography/mass spectrometry (UHPLC/MS)51. Metabolites were identified by automated comparison of ion features to a reference library of chemical standards followed by visual inspection for quality control52. Missing values were presumed to be below the limits of detection and were therefore imputed to the compound minimum.

Metabolomic data analysis was undertaken, having accounted for medicinal drug-related by-products and discarded unnamed molecules. We constructed a heat map illustrating group differences at individual and group average level using MetaboAnalyst 4.053. Univariate comparison of abundance of 947 named small-molecule (<500Da) metabolites between cases and controls was performed using the Mann-Whitney U test.