Extended exome sequencing identifies BACH2 as a novel major risk locus for Addison'

14 Pages·2016·1.06 MB·English

Original Article
doi:10.1111/joim.12569
BACH2
Extended exome sequencing identifies as a novel
major risk locus for Addison’s disease
D.Eriksson1,2,*,M.Bianchi3,*,N.Landegren1,4,J.Nordin3,F.Dalin1,4,A.Mathioudaki3,G.N.Eriksson5,
L.Hultin-Rosenberg3,J.Dahlqvist3,H.Zetterqvist3,4,(cid:1)A.Karlsson3,(cid:1)A.Hallgren1,4,F.H.G.Farias3,E.Mur(cid:3)en3,K.M.Ahlgren4,
A.Lobell4,G.Andersson6,K.Tandre4,S.R.Dahlqvist7,P.S€oderkvist8,L.R€onnblom4,A.-L.Hulting5,J.Wahlberg9,
O.Ekwall10,11,P.Dahlqvist7,J.R.S.Meadows3,S.Bensing2,5,K.Lindblad-Toh3,12,O.K€ampe1,2,4,†&G.R.Pielberg3,†
Fromthe1DepartmentofMedicine(Solna),CenterforMolecularMedicine,KarolinskaInstitutet;2DepartmentofEndocrinology,Metabolism
andDiabetesKarolinskaUniversityHospital,Stockholm;3ScienceforLifeLaboratory,DepartmentofMedicalBiochemistryandMicrobiology;
4ScienceforLifeLaboratory,DepartmentofMedicalSciences,UppsalaUniversity,Uppsala;5DepartmentofMolecularMedicineandSurgery,
KarolinskaInstitutet,Stockholm;6DepartmentofAnimalBreedingandGenetics,SwedishUniversityofAgriculturalSciences,Uppsala;
7DepartmentofPublicHealthandClinicalMedicine,Ume(cid:1)aUniversity,Ume(cid:1)a;8DepartmentofClinicalandExperimentalMedicine;
9DepartmentofEndocrinology,DepartmentofMedicalandHealthSciences,DepartmentofClinicalandExperimentalMedicine,Linko€ping
University,Linko€ping;10DepartmentofPediatrics,InstituteofClinicalSciences,SahlgrenskaAcademy;11DepartmentofRheumatologyand
InflammationResearch,InstituteofMedicine,SahlgrenskaAcademy,UniversityofGothenburg,Gothenburg,Sweden;and12BroadInstitute
ofMITandHarvard,Cambridge,MA,USA
Abstract. Eriksson D, Bianchi M, Landegren N, Methods. To understand the genetic background of
Nordin J, Dalin F, Mathioudaki A, Eriksson GN, Addison’sdisease,weutilizedtheextensivelychar-
Hultin-Rosenberg L, Dahlqvist J, Zetterqvist H, acterizedpatientsoftheSwedishAddisonRegistry.
Karlsson (cid:1)A, Hallgren (cid:1)A, Farias FHG, Mure(cid:3)n E, We developed an extended exome capture array
Ahlgren KM, Lobell A, Andersson G, Tandre K, comprising a selected set of 1853 genes and their
Dahlqvist SR, So€derkvist P, Ro€nnblom L, Hulting potential regulatory elements, for the purpose of
A-L, Wahlberg J, Ekwall O, Dahlqvist P, Meadows sequencing 479 patients with Addison’s disease
JRS, Bensing S, Lindblad-Toh K, Ka€mpe O, and 1394controls.
Pielberg GR (Karolinska Institutet; Metabolism
and Diabetes Karolinska University Hospital, Results. We identified BACH2 (rs62408233-A,
Stockholm; Uppsala University; Swedish Univer- OR = 2.01 (1.71–2.37), P = 1.66 9 10(cid:1)15, MAF
sity of Agricultural Sciences, Uppsala; Ume(cid:1)a 0.46/0.29 in cases/controls) as a novel gene
University, Ume(cid:1)a; Linko€ping University, Linko€ping; associated with Addison’s disease development.
University of Gothenburg, Gothenburg; Sweden; We also confirmed the previously known associa-
Broad Institute of MIT and Harvard, Cambridge, tions withtheHLAcomplex.
MA, USA). Extended exome sequencing iden-
Conclusion.WhilstBACH2hasbeenpreviouslyreported
tifies BACH2 as a novel major risk locus for
Addison’s disease. J Intern Med 2016; 280: 595– toassociatewithorgan-specificautoimmunediseases
co-inherited with Addison’s disease, we have identi-
608.
fiedBACH2asamajorrisklocusinAddison’sdisease,
independent of concomitant autoimmune diseases.
Background. Autoimmune disease is one of the lead- Our results may enable future research towards
ingcausesofmorbidityandmortalityworldwide.In preventivediseasetreatment.
Addison’s disease, the adrenalglands are targeted
by destructive autoimmunity. Despite being the Keywords: AddisonDisease,BACH2protein,genetic,
most common cause of primary adrenal failure, genetic association studies, genetic predisposition
littleis knownaboutits aetiology. todisease, human, polymorphism.
insufficiency, and it is referred to as autoimmune
Introduction
Addison’s disease (AAD) [1]. The ensuing loss of
Autoimmune destruction of the adrenal cortex is adrenal hormones is ultimately lethal if cortisol is
the predominating cause of primary adrenal not replaced [2]. The prevalence of AAD amongst
Caucasiansis87–144permillion[3–9],andforthe
last decades, an increasing incidence of AAD has
*Theseauthorscontributedequallytothiswork.
†Theseauthorsjointlysupervisedthiswork. been reported [3, 6, 7, 10]. Typically, the disease
ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine. 595
ThisisanopenaccessarticleunderthetermsoftheCreativeCommonsAttribution-NonCommercial-NoDerivsLicense,
whichpermitsuseanddistributioninanymedium,providedtheoriginalworkisproperlycited,theuseisnon-commercialand
nomodificationsoradaptationsaremade.
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
onset occurs in the fourth or fifth decade of life, regulatory elements, as well as their splice sites.
althoughindividualsofallagescanbeaffected[8], Using a well-defined AAD phenotype, based on all
withapreponderanceforwomen(60%–65%)[8,11]. available diagnostic modalities ranging from clini-
cal signs to autoantibody markers, we identified
Despite its rarity, AAD displays high concordance BACH2asanovelmajorrisklocusinAAD.
rates in twins and aggregation in families,
suggesting a strong genetic contribution [12–18]. Materials and methods
A recent twin concordance study (Skov J. et al.,
Studysubjects
submitted) reports a heritability higher than
that of, for instance, type 1 diabetes mellitus All participants in this study were recruited in
(0.69–0.88) [19–23] and Graves’ disease (0.79) Sweden.Intotal,1501controlswereretrievedfrom
[24] . directedsamplingofhealthyindividuals(n = 653),
or from blood donors (n = 848) [29], not known to
The success of genetic case–control disease stud- sufferfromAAD.Allthe700AADsamplesincluded
ies, including autoimmune diseases, depends to a were obtained from the Swedish Addison Registry
large extent on the phenotypic homogeneity of the (SAR). SAR incorporates detailed clinical informa-
case group. Diagnosis of AAD can be confirmed tion and serological data compiled by the respon-
with a serological marker, namely 21-hydroxylase sible physician. All cases fulfilled diagnostic
autoantibodies in patient sera [2, 25]. 21-hydro- criteria of primary adrenal insufficiency, that is
xylase autoantibodies are present in about 86%of low morning serum cortisol and elevated plasma
patients presenting with idiopathic primary adre- adrenocorticotropic hormone (ACTH) [2, 30]. In
nalinsufficiency,withaspecificityofabout95%in equivocal cases, insufficient response to corti-
experienced labs [8, 26]. Therefore, patients with cotropin stimulation confirmed the diagnosis.
AAD are likely to constitute one of the most Blood and serum samples were stored at (cid:1)80 °C
homogenous autoimmune case groups available untiluse.
for geneticstudies.
Tostrictlydefinethecasegroup,weexcludedcases
The Swedish Addison Registry (SAR) is one of the where other causes of adrenal failure were sus-
world’s largest AAD biobanks. At the commence- pected:caseswithfragmentaryclinicalinformation,
mentofthisstudy,theregistryincluded700cases, adrenoleukodystrophy/adrenomyeloneuropathy,
representing approximately half of the estimated autoimmune polyendocrine syndrome type 1
patients with AAD currently living in Sweden [3]. (APS1), adrenal failure secondary to corticosteroid
The majority of these patients have concomitant treatment, congenital adrenal hyperplasia (CAH),
tissue-specific autoimmune diseases such as type adrenal metastases, isolated ACTH deficiency,
1 diabetes, autoimmune thyroid disease or perni- hypogonadotropic hypogonadism with adrenal
ciousanaemia.Thehighrateofcomorbiditymakes insufficiencyduetoDAX1mutationsortuberculosis
it possible to study the shared genetic basis of andcasesthathadundergoneadrenalectomy.Fur-
autoreactivity. thermore, AAD cases with autoantibodies against
interferon-a, interferon-x or interleukin-22, and
Similar to most other autoimmune diseases, AAD cases lacking 21-hydroxylase autoantibodies were
has a complex inheritance. Apart from the well- also excluded (Figure S1 in the supplementary
characterized association with the HLA complex information).AADcaseswithanticytokineantibod-
[27],littleisknownaboutthegeneticvariantsthat ies were excluded because of potentially undiag-
contribute todisease development [27].Untilnow, nosedAPS1.Toprioritizespecificityoversensitivity,
candidate gene studies have linked genes impli- we required positive results from two 21-hydroxy-
catedinco-occurringautoimmunediseasestoAAD laseassays,performedattwodifferentlaboratories.
[27, 28]. The rarity of AAD has however, made
extensive unbiased genomewide association stud- All study subjects gave their informed consent.
iesunfeasible. The study was performed in accordance with the
Declaration of Helsinki and was approved by the
Here, we performed exome sequencing of a com- local ethics committees (2008/296-31/2, M177-
prehensive selection of 1853 genes, extended to 07, 2012/52-31, 2009/013, 08-109M § 137/08,
include their promoters and other potential 2010/182-32).
596 ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine.
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
to 400-bp fragments (Covaris E220). 30 and 50
GenomicDNA samples
overhangswererepairedusingT4DNApolymerase
GenomicDNAwasisolatedfromwholebloodEDTA (EP0062; Thermo Fisher Scientific), Klenow frag-
usingQIAmpBloodMidiKit(51185;Qiagen,Venlo, ment (EP0052; Thermo Fisher Scientific) and T4
Netherlands)followingthemanufacturer’sprotocol, polynucleotide kinase (EK0032; Thermo Fisher
or using a LGC Genomics extraction service. The Scientific).AdenylationwithKlenowfragmentexo(cid:1)
concentration of the DNA was measured using (EP0422; Thermo Fisher Scientific) and ligation of
NanoDrop ND-1000 Spectrophotometer or Qubit NEXTflex-96TM DNA barcode adapters (514106;
2.0Fluorometer(ThermoFisherWaltham,MA,US) Bioo Scientific Austin, TX, US) using T4 DNA
withQubitdsDNAHSAssayKit(Q32851,Thermo- Ligase (EL0011; Thermo Fisher Scientific) were
Fisher). The rough assessment of high molecular performed according to the manufacturer’s proto-
weightDNAwasperformedbyelectrophoreticsep- col. Samples with a final fragment library concen-
arationof150 ngofgenomicDNAon1%agarosegel. tration of less than 130 ng lL(cid:1)1 (n = 118 control
samples) were amplified 4–7 cycles using Phusion
High-Fidelity PCR Master Mix (F531S; Thermo
Sequence capture array
Fisher Scientific) prior to hybridization. Pools of
A custom SeqCap EZ Choice XL Library eight samples were hybridized according to Roche
(06266517001; Roche NimbleGen, Basel, Switzer- NimbleGen protocol and sequenced with 100-bp
land) was designed targeting 1853 genes. Genes paired-end reads using Illumina HiSeq 2500 ver-
were selected based on their known role in basic sion 3or version4 chemistry.
immunefunctions,inflammationandautoimmune
disease.Asystematicinclusionofgenesassociated
Bioinformatic analysis
with immunological or autoimmune disease, for
instance type 1 diabetes and coeliac disease, was Rawreadsfromindividualsamplesweremappedto
madebasedoncandidatelociidentifiedingenome- the hg19 human reference genome with the Bur-
wideassociationstudies.Bysearchingthepathway rows-Wheeler aligner 0.7.4 [32], and duplicate
databasesKEGG(www.genome.jp/kegg)andInge- reads were marked by Picard tools 1.92 (http://
(cid:2)
nuity (QIAGEN, www.qiagen.com/ingenuity) with broadinstitute.github.io/picard). GATK 2.7.2 was
keywords such as autoimmunity, cytokines and T applied for realignment around indels [33–35].
cells,additionalgenesinvolvedinimmunepathways Alignment quality was evaluated using Samtools’
were retrieved and selected. Genes causing mono- flagstat [36] and Picard tools CalculateHsMetrics.
genicadrenaldiseasewerealsotargeted.Thegene Samples with mean target coverage less than 109
intervalsforallalternativetranscriptswereretrieved were excluded. Variants were called on individual
from human assembly NCBI36/hg18 using UCSC samples by GATK 3.2.2 HaplotypeCaller in gVCF
GenomeTableBrowser(https://genome.ucsc.edu/ mode, according to the procedures recommended
cgi-bin/hgTables) known gene and microRNA for single samples all-sites calling [34, 35]. Geno-
tracks. Specific regions ofinterest werethefollow- typingwasmadejointlyforthecompletecasegroup
ing:50and30UTRs,codingexons,potentialpromoter with GATK 3.2.2 GenotypeGVCFs, and thereafter,
regions ((cid:3)2 kb from transcription start sites of all only single nucleotide polymorphisms (SNPs) were
alternative transcripts) and splice sites (20 bp of considered in further analyses. Genotype filters
intronicsequencesadjacenttoexons).Inaddition, were applied with VCFtools [37] requiring a mini-
we targeted a selection of potential regulatory ele- mum read depth of 8 and a minimum genotype
mentsfromtheextendedgeneintervals((cid:3)100 kbof Phredqualityscoreof20[38].SNPswereassigned
thegene).Theseelementsweredefinedbasedon29 probability scores using the GATK 3.2.2 Vari-
mammals conservation, and the elements with antRecalibrator and filtered at tranche level 99.0.
SiPhyLOD-scorehigherthan7wereincluded[31]. Furthermore, heterozygous calls were kept if their
Thetotalsizeofthetargetedregionsofinterestwas variant allele balance across all samples was
>36 Mb,whereasthedesignedarrayprobescovered between 0.2 and 0.8. Positions deviating from
89.6%ofit(97.7%includingthe100-bpoffset). Hardy–Weinberg equilibrium (P < 10(cid:1)6) and
monomorphic sites were excluded. SNPs were
annotated with identifiers from dbSNP 137 [39]
Sequence capture experiment
andhadtheirimpactpredictedwithSnpEff4.1[40].
Thesequencinglibrarywaspreparedbyultrasoni- Annotations for noncoding variants were retrieved
ficationofupto1 lgofhighmolecularweightDNA using HaploReg v4.1 [41] and RegulomeDB [42].
ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine. 597
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
Furthermore, we used the IGV [43] and UCSC and single-variant association analyses. Lastly,
genome browser [44] for sequence alignments and variants were pruned out if their call rate in cases
functionalinvestigation,respectively. versus controls differed more than ten percentage
points (Figure S6).
Sample ancestry and relatedness estimation
Association analyses – common variants
We determined the study subjects’ ancestry based
on off-target genetic differences using the LASER Logistic regression was applied for association
software, University of Michigan [45] (Figure S2). analyses of common single variants with a minor
Wesetthedimensionsofthereferencespace(k)to allele frequency above 1% (MAF > 0.01) [47]. The
4, and we used 20 principal components (K) for logisticregressionanalysiswasperformed withan
calculation and projection of the study samples additive model, which is the most common model
after 10 iterations (r). Study subjects were super- for unknown modes of inheritance [48, 49]. Prin-
imposedontheHumanGenomeDiversityProjects cipalcomponentanalysis(PCA)wascomparedwith
referencepanelandpredictednon-Europeanswere multidimensionalscaling(MDS) incorrectingpop-
excluded if they did not fall within the European ulationstratificationinourdataset,withamethod
reference coordinates defined by Wang et al. [45]. adapted from Wang et al. [50] (Figures S7–S9,
Intersample relatedness coefficients were inferred Table S1). In our data set, MDS performed better
pairwise with VCFtools [37] that implements the in decreasing the inflation factor (Figure S10) and
algorithm by Manichaikul et al. [46]. The coeffi- the three first MDS-generated dimensions were
cients representing first-degree relatedness included as regression covariates to account for
(Φ > 0.177) were used to identify closely related population substructure. Alternative covariates
ij
pairs. The sample with the lowest mean target such as sex and amplification status were also
coverage in case–case and control–control pairs tested consecutively. Lambda values were calcu-
was excluded, as well as both samples in case– latedwithandwithoutthetop-associatedloci,and
control pairs. Furthermore, seeming twins, not P values were adjusted by means of genomic
bornon thesame day,were discarded. control [51] to account for cryptic relatedness
inflating the test statistic (Tables S2–S3). The
MDS and single-variant association analysis were
Quality control
performedinPLINK[52].PCAwascalculatedusing
Samples with discordance between reported sex EIGENSOFT[53].LDstructuresweremeasured(r2)
and sex inferred from sequencing data were and visualized withHaploview[54].
excluded. Furthermore, samples with a rate of
missing data, heterozygosity ratio, or transition–
Rarevariant gene-based test
transversion ratio deviating more than five stan-
dard deviations from the mean, were discarded. Togainstatisticalpowertoassociaterare variants
Finally, samples with singleton counts, calculated with disease, variants have to be evaluated collec-
onvariants called inthemajority of samples (0.95 tively in predefined groups. We aggregated genetic
call rate), deviating by five standard deviations or variants by the genes in which they were located
more, were removed. With this adjustment, we and calculated the risk associated with each gene
were able to pinpoint both low number of single- and intergenic space using the SKAT statistic
tons revealing duplicates or contamination and implemented in the SCORE-Seq software [55].
high singleton counts indicating low-quality DNA Aggregates including less than two variants were
or suggestingforeign ancestry (FiguresS3–S5). disregarded. Using gene boundaries defined by
SnpEff-based annotation, comprised of Ensembl
After removing all low-quality samples, variants GRCh37 genes and intergenic regions, we tested
with more than 20% missing data across all theeffectofrarevariants(MAF ≤ 0.01)[47]onAAD
subjects were excluded as well as subjects with susceptibility. We made two analyses with rare
more than 35% missing data (Figure S6). We set variants. First, we tested the whole set of rare
these filter thresholds for call rates to reach final variants, and in a second analysis, a subset
variantandindividualcallrates≥95%.Bylimiting comprised of only moderate- and high-impact
the number of variants and individuals discarded, variants. In both tests, we included three covari-
whilst keeping adequate average call rates, we ates of MDS dimensions but no genomic control
aimed at increasing power in both aggregate tests [56]. We used a Bonferroni-corrected threshold
598 ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine.
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
(3.3 9 10(cid:1)6, a = 0.05 corrected for the total num- 1%, across all study participants. To account for
ber of aggregates in the two analyses; 15 089) to multipletesting,weusedawell-establishedstatis-
evaluatethestatisticalsignificanceofourfindings. ticalsignificancethresholdforgenomewideassoci-
ation studies, P = 5 9 10(cid:1)8 [48], which is
conservative for our extended exome study. MDS
InsilicoHLA typing
plots revealed population substructures, separat-
HLA types were inferred from de novo assembly of ing study subjects sampled in northern versus
sequencing reads using MOSAIK aligner and the central/southernSweden(FigureS7).Theoverdis-
ATHLATES software [57, 58]. HLA reference defi- persionoftheteststatisticwaspartiallyreducedby
nitions were updated to the IMGT/HLA release includingthefirstthreeMDS-generateddimensions
3.22.0 with sequences aligned 10 October 2015. as covariates, andresidual inflationwas corrected
AllelesatlociA,BandCofclassI,andDPA,DPB, withgenomiccontroltoaccountforcrypticrelated-
DQA, DQB and DRB of class II were deduced. ness(k = 1.02)(TablesS3,FiguresS11–S12).
GC
Haplotype frequencies were compared both per
haplotype, and per haplotypic pair. P values were The single-variant association test identified two
computed from v2-test in comparison with the peaks of association on chromosome 6 (Fig. 1). In
mostequallydistributedhaplotype.Theoddsratio consecutive analyses, using no or alternative
95% confidence intervals were calculated accord- covariates, the peaks remained statistically signif-
ing toSzumilas et al.[59]. icant. With sequential conditioning of the lead
SNPs, we identified three statistically significant
independent signals (Table 1). We identified the
Results
geneBACH2(BTBdomainandCNChomolog2)asa
Our study initially included 700 primary adrenal novel AAD risk locus (rs62408233-A, OR = 2.01
insufficiency case samples and 1501 control sam- (1.71–2.37), P = 1.66 9 10(cid:1)15, MAF 0.46/0.29 in
ples.Wesetourfocusoncharacterizingthegenetic cases/controls). Minor allele frequencies from the
background of complex inherited autoimmune 1000 Genomes Project [60] and call rates for the
Addison’s disease (AAD). Therefore, we stringently associated SNPs in BACH2 are presented in
excluded173caseswhereothercausesofprimary Table S4. Furthermore, we found two signals of
adrenal failure were suspected, including all the association in the HLA region (rs41315836-G,
119 cases that were negative for 21-hydroxylase OR = 0.20 (0.13–0.29), P = 5.73 9 10(cid:1)16;
autoantibodies (Materials and methods), thereby rs17221059-A, OR = 2.29 (1.84–2.84),
ensuring a case group with homogeneous disease P = 9.59 9 10(cid:1)13). A comprehensive summary of
pathoaetiology. We also excluded 25 case and 28 top-associated SNPs is presented in Table 1 and
control samples estimated to be either of non- conditionaltestsinFigureS13.
European ancestry or highly related. After the
application of all quality control measures, 479
Separating AAD fromconcomitant autoimmunediseases
case and 1394 control samples remained for fur-
ther analyses. Of the 479 cases, 298 (62.2%) were BACH2 had previously been reported to have an
female, and the median age of onset was 34 years association withothertissue-specific autoimmune
(range 3–78). Figures S1–S5 depict in detail the diseases such as type 1 diabetes [61], autoim-
sample quality control process. The average target mune thyroid disease [62], coeliac disease [63],
coverage of the remaining individuals was inflammatory bowel disease [64] and vitiligo [65].
37.29 (SD12.9). Because AAD cases have a high degree of autoim-
mune comorbidity, the potential risk of confound-
Thefinaldatasetincluded103 120commonSNPs ing in an assorted AAD case group is substantial.
(MAF > 1%) and 304 338 rare SNPs (MAF ≤ 1%). It was also known that BACH2 had been associ-
The mean variant call rate (96%) and the mean atedwiththepresenceofthyroidperoxidase(TPO)
individual call rate (95%) indicated a data set of autoantibodies[62].BecausetheSwedishAddison
high quality. Registry (SAR) includes detailed clinical informa-
tion and SAR samples have been assayed for TPO
autoantibodies, we were able to identify cases
BACH2 isanovel AAD susceptibility locus
diagnosed with any, present or previous, autoim-
We conducted a single-SNP association analysis mune diseases besides AAD, and furthermore, to
comparingcommonvariantswithMAFhigherthan identify TPO autoantibody-positive cases [66].
ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine. 599
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
Fig.1 Associationbetweensingle-nucleotidepolymorphisms(SNPs)andautoimmuneAddison’sdisease(AAD),exhibited
asPvaluesandgenomiclocation.TheassociationbetweenSNPsandAAD,ascalculatedbyanadditivelogisticregression
model,is plottedas the negativelogof theP valueagainstchromosomallocation.ThreeMDS generated dimensionsare
included as covariates and P values are corrected with genomic control. The horizontal red line denotes the statistical
significancelevelof5910(cid:1)8.Bandsinshadesofbluedepictthechromosomeboundaries.
Table1 LeadvariantsindependentlyassociatedwithautoimmuneAddison’sdiseaseinsingle-variantassociationanalysis
Minorallele
Minor/Maj. frequency Association
Leadvariant Chromosome Position Locus allele Cases Controls OR(95%CI)a Pvalueb
rs41315836 6 31975151 HLA G/T 0.04 0.17 0.20(0.13–0.29) 5.73910(cid:1)16
rs17221059 6 33048348 HLA A/G 0.20 0.09 2.29(1.84–2.84) 9.59910(cid:1)13
rs62408233 6 90976609 BACH2 A/G 0.46 0.29 2.01(1.71–2.37) 1.66910(cid:1)15
aThe odds ratio (OR) represents the effect of the minor allele on the odds of observing AAD, under an additive logistic
regressionmodel,suchthatanORgreaterthanoneimpliesthattheminoralleleincreasesodds.bThepresentedPvalues
arecalculatedwithpopulationsubstructurecovariatesandcorrectedwithgenomiccontrol.
Recalculatingthelogisticregressionusingonlythe level (P < 5 9 10(cid:1)8), 26 of which were located in
cases with isolated AAD and without TPO autoan- BACH2. Conditioning on the top SNP in BACH2
tibodies (119 cases), the signal in BACH2 deflated all association signals at this locus,
remained statistically significant (P = 3.9 9 10(cid:1)8, underlining the linkage between associated mark-
Figure S14). TPO autoantibody-positive and ers in the region (Figure S13). The associated
autoantibody-negative cases had the same mean haplotype block (HSA6:90,815,190-91,029,762)
age, lowering the possibility that negative cases covers the 50 noncoding exons and introns of
had a genetic risk that was not yet penetrant BACH2 (Fig. 2). In general, it is complex to
(Figure S15, Table S5). decipher which single genetic variants have an
effect on disease development. The presence of
linkage disequilibrium and the difficulty of detect-
FromBACH2 variants todisease
ing the effects of single variants confound the
From the total of 103 120 common variants, 40 dissection of which genes have an important
exceeded our stringent statistical significance impact on disease predisposition. Therefore, we
600 ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine.
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
(a)
(b)
(c)
(d)
Fig.2 Association between single-nucleotide polymorphisms (SNPs) and autoimmune Addison’s disease (AAD), the
locationofnearbygenes,conservation,andthesurroundinglinkagedisequilibrium(LD)structure,inaregionofchromosome
6.ShownisaregionaroundBACH2anditsneighboringgenesonchromosome6q15.Inpanel(a),theassociationbetween
SNPsandAAD,ascalculatedbyanadditivelogisticregressionmodel,isplottedasthenegativelogofthePvalueagainst
chromosomal position. The horizontal red line denotes the level considered statistically significant, 5910(cid:1)8. SNPs are
coloredaccordingtotheirlevelofgenotypecorrelation(r2)withthemostsignificantSNP(redsquare,rs62408233),which
hasaPvalueof1.7910(cid:1)15.Panel(b)depictsthegenomiclocationsofthegenespresentintheregion,and(c)atrackwith
basewiseconservationbetween100vertebrates,bothfromthegenomebrowseratUCSC,assemblyhg19.Arrowsindicate
the direction of transcription (right to left for BACH2). The main BACH2 transcript, shown here, has three 50 non-coding
exons.Panel(d)illustratesthelinkagedisequilibrium(LD),asmeasuredbypairwisecorrelation(r2)ofSNPs,andplottedby
Haploview.DarkershadingmarkshighergenotypecorrelationandstrongerLD.
compiled a table with available epigenetic data associated SNPs (Table S6). Of the 26 associated
from both HaploReg v4.1 [41] and RegulomeDB SNPs in this region, all except two overlap
[42] to highlight relative functional importance of enhancer histone marks in blood (B and T
ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine. 601
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
lymphocytes) and four are associated with BACH2
Rarevariant contribution toAAD
eQTLs.
We included all detected variants with MAF lower
thanorequalto1%inagene-basedaggregatetest.
HLA typing confirmedknown AADrisk alleles
First, we performed association analysis with the
Using ATHLATES for HLA typing from sequencing whole set of variants. Thereafter, we selected only
reads [57], we recapitulate previously established high- and moderate-impact variants as enriching
Scandinavian associations in the HLA class II forharmfulallelescouldincreasestatisticalpower
region (Table 2). For instance, we could confirm in aggregate tests. The inflation factors (k) were
the risk associated with the haplotypes 1.07and1.02,respectively(TableS3).Neitherwith
DQB1*02:01-DQA1*05:01-DRB1*03:01, corre- the inclusion of all variants nor with only high- or
sponding to the serotype DR3-DQ2, and DQB1 moderate-impact variants didany geneexceed the
*03:02-DQA1*03:01-DRB1*04:04 corresponding statistical significance threshold (P = 3.3 9 10(cid:1)6)
to the serotype DR4-DQ8 [8, 67]. We could also (Figure S16). Hence, no gene aggregate, powered
confirm the protective effect of DQB1*06:03- withonlyrarevariants,wasassociatedwithAADin
DQA1*01:03-DRB1*13:01. Moreover, we could this study.
replicate the association with the risk allele
DRB1*04:03, previously identified in Russian
Discussion
patients with AAD [68]. This allelic association
was also present on the haplotype DQB1*03:02- This is the first large-scale study aiming at identi-
DQA1*03:01-DRB1*04:03. fying the genetic background of autoimmune
Addison’s disease (AAD). For this purpose, we
By further analysing haplotypes in the pairs in developed an array to capture a comprehensive
which they occur, we could confirm that homozy- selection of potential causative genes and sur-
gosityforriskhaplotypesislinkedtoanincreased rounding noncoding regions. Our Swedish collec-
susceptibility to AAD (6.4 < OR < 8.8) (Table 3). tion of AAD samples is one of the largest available
The pair of haplotypes associated withthe highest for investigation, and the detailed phenotypic
riskwasthecompoundheterozygousDQB1*02:01- characterization of all cases underpins the foun-
DQA1*05:01-DRB1*03:01 and DQB1*03:02-DQ dation forthis study.
A1*03:01-DRB1*04:04 (OR = 18 (7.2–46), P =
4.9 9 10(cid:1)12), also concordant with current knowl- Weidentifiedanovelrisklocus,BACH2,associated
edge[8,67,69]. with AAD and confirmed a well-established
Table2 MostfrequentHLAclassIIhaplotypesandtheirassociationswithautoimmuneAddison’sdisease
Cases Controls
DQBDQADRB n(%)a n(%)a Oddsratio(95%CI) Pvalueb
02:0105:0103:01 358(37.4) 341(12.2) 2.7(1.8–4.1) 3.7910(cid:1)6***
03:0203:0104:04 174(18.2) 149(5.3) 3(1.9–4.6) 1.9910(cid:1)6***
06:0201:0215:01 93(9.7) 430(15.4) 0.55(0.35–0.86) 1.2910(cid:1)2*
03:0203:0104:01 66(6.9) 172(6.2) 0.98(0.6–1.6) 1.0
04:0204:0108:01 35(3.7) 89(3.2) 1(ref)
03:0203:0104:03 18(1.9) 12(0.4) 3.8(1.7–8.7) 2.1910(cid:1)3**
06:0301:0313:01 15(1.6) 230(8.2) 0.17(0.086–0.32) 1.2910(cid:1)8***
Otherhaplotypes 119(20.8) 1365(49.0)
Total 958(100.2) 2788(99.9)
aThe table is restricted to include haplotypes with at least 10 cases and 10 controls. The remaining haplotypes are
summedin‘Otherhaplotypes’.bPvaluesfromthev2distributionarecalculatedwiththecountsfromthemostequally
distributedHLAhaplotypeasreference(ref).*Pvaluesmallerthan0.05,**Pvaluesmallerthan0.01,***Pvaluesmaller
than0.001.
602 ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine.
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
Table3 MostfrequentpairsofHLAclassIIhaplotypesandtheirassociationswithautoimmuneAddison’sdisease
Cases Controls Oddsratio
DQBDQADRB/DQBDQADRB n(%)a n(%)a (95%CI) Pvalueb
02:0105:0103:01/03:0203:0104:04 108(22.5) 20(1.4) 18(7.2–46) 4.9910(cid:1)12***
02:0105:0103:01homozygous 32(6.7) 17(1.2) 6.4(2.4–17) 3.0910(cid:1)4***
02:0105:0103:01/03:0105:0511:01 17(3.5) 7(0.5) 8.2(2.5–27) 6.9910(cid:1)4***
02:0105:0103:01/03:0105:0512:01 13(2.7) 5(0.4) 8.8(2.4–32) 1.5910(cid:1)3**
03:0203:0104:04homozygous 13(2.7) 5(0.4) 8.8(2.4–32) 1.5910(cid:1)3**
02:0105:0103:01/03:0103:0304:01 16(3.3) 14(1.0) 3.9(1.3–11) 2.3910(cid:1)2*
02:0105:0103:01/03:0203:0104:01 28(5.8) 31(2.2) 3.0(1.2–7.8) 3.1910(cid:1)2*
02:0105:0103:01/06:0201:0215:01 38(7.9) 52(3.7) 2.5(1–6) 7.0910(cid:1)2
02:0105:0103:01/04:0204:0108:01 12(2.5) 14(1.0) 2.9(0.96–8.7) 0.10
04:0204:0108:01/06:0201:0215:01 11(2.3) 13(0.9) 2.9(0.93–8.8) 0.12
03:0103:0304:01/03:0203:0104:01 6(1.3) 6(0.4) 3.4(0.85–13) 0.16
03:0203:0104:01/03:0203:0104:04 9(1.9) 12(0.9) 2.5(0.79–8.2) 0.2
06:0201:0215:01/03:0105:0512:01 6(1.3) 9(0.6) 2.2(0.61–8.3) 0.37
02:0105:0103:01/06:0301:0313:01 5(1.0) 33(2.4) 0.51(0.15–1.7) 0.44
03:0203:0104:04/06:0301:0313:01 5(1.0) 8(0.6) 2.1(0.54–8.3) 0.47
03:0203:0104:04/06:0201:0215:01 8(1.7) 27(2.1) 1(ref)
Otherhaplotypepairs 152(31.7) 1121(80.4)
Total 479(99.8) 1394(99.7)
aThetableisrestrictedtoincludehaplotypeswithatleast5casesand5controls.Theremaininghaplotypesaresummed
in ‘Other haplotype pairs’. bP values from the v2 distribution are calculated with the counts from the most equally
distributedHLAhaplotypeasreference(ref).*Pvaluesmallerthan0.05,**Pvaluesmallerthan0.01,***Pvaluesmaller
than0.001.
contribution of HLA risk alleles to the risk of characterization, BACH2 emerges as a major gene
developingthisdisease.Furthermore,wewereable candidate in organ-specific autoimmunity. Here,
to show that neither autoimmune comorbidity nor wepresenta215-kbBACH2regionassociatedwith
the presence of TPO autoantibodies was required AAD. Even though the new association does not
for the strong association signal in BACH2. From directly infer a causal relationship with AAD, it
the rare variant association analysis, we revealed providesdirectionforfurtherinvestigations.The50
noadditional candidateloci. noncoding region of BACH2, including both
upstream and intronic sequence and potentially
BACH2 encodes a basic leucine zipper transcrip- harbouring regulatory variants, poses a candidate
tion repressor with important functions in acti- locusfor functional studies.
vated B cells. Acting as a heterodimer with Maf
proteins, BACH2 is an important regulator of the The newly identified risk locus includes 26 previ-
immunoglobulinheavychaintranscription,aswell ously known common variants. Due to high link-
as class switching [70, 71]. Although initially age disequilibrium, itis complicated todissect the
linked to B-cell function, BACH2 has also been effect of single variants in the 50 noncoding region
foundtobecrucialforT-celldifferentiation.Asthe of BACH2. The vast majority of these variants are
most prominent T-cell super-enhancer, a lack of either located in enhancer/promoter regions in
BACH2 shifts the T-cell balance towards effector relevantcelltypesorhavebeenassociatedwiththe
cells (CD8+), at the expense of helper and regula- expressionofBACH2inlarge-scaleeQTL-mapping
tory T cells [72]. In the light of the BACH2 locus studies [73]. Therefore, considering both genetic
being associated with the development of other and epigenetic data, we are currently not able to
autoimmune diseases, and its functional define a single most important AAD-associated
ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine. 603
JournalofInternalMedicine,2016,280;595–608
D. Eriksson et al. BACH2 is a risk locus for Addison’s disease
BACH2 variant. Instead, we postulate that multi- are suggested tobe majorcontributors to complex
ple variants may contribute to the risk of develop- diseases [76].
ing AAD and could be selected for functional
studies.
Conclusion
Using in silico typing from sequencing reads, we WehaveidentifiedBACH2asamajorrisklocusin
could confirm the risk conferred by certain HLA AAD. Knowledge gained from this study may aid
haplotype pairs [8, 67–69]. The strongest risk thefuture development ofpreventive treatment.
factor is by far the combination of HLA alleles
DRB1*03:01 and DRB1*04:04, as previously
Author contributions
shown in other populations [8, 67–69]. Despite
the fact that we selected only patients with AAD The study was conceived by KL-T, OK and GRP.
reacting against the same autoantigen, 21-hydro- The study design and optimization was performed
xylase, we detect a great variability at the HLA byNL,JN,AM,EM,KMA,AL,GA,LR,JRSM,KL-T,
locus. Interestingly, the majority of AAD cases did OK and GRP. Acquisition of study material and
not carry any of the high-risk HLA alleles. This data was implemented by MB, NL, JN, FD, AM,
(cid:1) (cid:1)
suggests that it is not necessary to have complete GNE, LHR, JD, HZ, AK, AH, FHGF, EM, KT, SRD,
homogeneity at the HLA locus in order to elicit PS,LR,JRSM,A-LH,JW,OE,PD,SB,OKandGRP.
immune reactivity against 21-hydroxylase, in the Analysisandinterpretationofdatawascarriedout
context ofAAD. by DE, MB and GRP. The manuscript was drafted
by DE, MB, NL, OK and GRP, with input from all
Apart from HLA, other already published loci co-authors.
found to be associated with AAD in candidate
gene studies did not exceed our significance
Conflict ofinterest statement
threshold [27, 28]. This could be due to several
factors. First, a strict significance threshold is Dr. Ro€nnblom reports grants from AstraZeneca,
needed in our study to avoid false-positive asso- grants from the Swedish Research Council,
ciations when testing multiple hypotheses. Sec- grants from the KAW Foundation, grants from
ondly, associated alleles could be population- the Swedish Reumatism Foundation, grants from
specific and therefore may not be reproducible in the King Gustaf V0s 80-year Foundation, during
our study population. Lastly, directed candidate the conduct of the study, grants from UCB,
gene studies do not have the possibility to outside the submitted work. Dr. Olle Ka€mpe
account for population substructure derived from reports grants from the Swedish Research
genetic data and are therefore more prone to Council, the Novo Nordisk Foundation and the
spurious associations. Torsten and Ragnar So€derberg’s Foundations.
He is board member of Olink Biosciences. The
Even though correction for population stratifica- remaining authors had no conflict of interests to
tion is extensively enabled in our study, rare disclose.
variant association analysis is limited by the cur-
rentsamplesize[74].Becausenotallgenesandnot
Acknowledgements
all noncoding regions could be selected in our
study,whole-genomesequencingwouldbeadesir- We would like to thank all patients and control
able nextstep. Thatwould alsoenable theassess- subjects participating in this study for making
mentofthepotentialroleofstructuralvariationin this research possible. Sequencing was performed
thediseasedevelopment[75].Functionalstudiesof bythe SNP&SEQTechnology Platform inUppsala,
potentially regulatory variants and replication in supported by the Swedish Research Council and
otherpopulationsshouldbeperformedtobuildon the Knut and Alice Wallenberg Foundation. The
theassociation withBACH2. facility is part of the National Genomics Infras-
tructure (NGI) Sweden and Science for Life Labo-
The current study adopts the benefits of high- ratory. Computational resources were provided by
throughputsequencingtotackletheshortcomings the Swedish National Infrastructure for Comput-
ofcandidategenestudies.Byextendingtheexome ing (SNIC) at Uppsala Multidisciplinary Center for
with noncoding conserved elements, we gain the AdvancedComputationalScience(UPPMAX),clus-
opportunity to identify regulatory variants, which ter and storage project b2014026. We would
604 ª2016TheAuthors.JournalofInternalMedicinepublishedbyJohnWiley&SonsLtdonbehalfofAssociationforPublicationofTheJournalofInternalMedicine.
JournalofInternalMedicine,2016,280;595–608