SNP Variants in Major Histocompatibility Complex Are Associated with Sarcoidosis Susceptibility—A Joint Analysis in Four European Populations

Sarcoidosis is a multiorgan inflammatory disorder with heritability estimates up to 66%. Previous studies have shown the major histocompatibility complex (MHC) region to be associated with sarcoidosis, suggesting a functional role for antigen-presenting molecules and immune mediators in the disease pathogenesis. To detect variants predisposing to sarcoidosis and to identify genetic differences between patient subgroups, we studied four genes in the MHC Class III region (LTA, TNF, AGER, BTNL2) and HLA-DRA with tag-SNPs and their relation to HLA-DRB1 alleles. We present results from a joint analysis of four study populations (Finnish, Swedish, Dutch, and Czech). Patients with sarcoidosis (n = 805) were further subdivided based on the disease activity and the presence of Löfgren’s syndrome. In a joint analysis, seven SNPs were associated with non-Löfgren sarcoidosis (NL; the strongest association with rs3177928, P = 1.79E−07, OR = 1.9) and eight with Löfgren’s syndrome [Löfgren syndrome (LS); the strongest association with rs3129843, P = 3.44E−12, OR = 3.4] when compared with healthy controls (n = 870). Five SNPs were associated with sarcoidosis disease course (the strongest association with rs3177928, P = 0.003, OR = 1.9). The high linkage disequilibrium (LD) between SNPs and an HLA-DRB1 challenged the result interpretation. When the SNPs and HLA-DRB1 alleles were analyzed together, independent association was observed for four SNPs in the HLA-DRA/BTNL2 region: rs3135365 (NL; P = 0.015), rs3177928 (NL; P < 0.001), rs6937545 (LS; P = 0.012), and rs5007259 (disease activity; P = 0.002). These SNPs act as expression quantitative trait loci (eQTL) for HLA-DRB1 and/or HLA-DRB5. In conclusion, we found novel SNPs in BTNL2 and HLA-DRA regions associating with sarcoidosis. Our finding further establishes that polymorphisms in the HLA-DRA and BTNL2 have a role in sarcoidosis susceptibility. This multi-population study demonstrates that at least a part of these associations are HLA-DRB1 independent (e.g., not due to LD) and shared across ancestral origins. The variants that were independent of HLA-DRB1 associations acted as eQTL for HLA-DRB1 and/or -DRB5, suggesting a role in regulating gene expression.

Introduction

Sarcoidosis (MIM 609464) is a multiorgan inflammatory disorder of unknown etiology. The majority of the sarcoidosis patients (of European descent) have a favorable prognosis if involvement is solely pulmonary, but the clinical picture and prognosis vary (1). Current understanding views sarcoidosis primarily as a multifactorial disorder with heritability estimates up to 66%, with possible environmental triggers in addition to the susceptibility gene(s) (2–4).

Sarcoidosis is manifested by accumulation of activated CD4-positive T lymphocytes and macrophages at disease sites suggesting a functional role for antigen-presenting molecules and immune mediator genes (4). The search for genetic components, indicated on the basis of family and multi-ancestral studies (2, 5), has shown a strong role of the major histocompatibility complex (MHC) region at chromosome 6p21.3 in susceptibility to sarcoidosis. Recently, several genome-wide association studies have indicated other regions of interest as well (6–8), but none as influential as the MHC.

Clinical manifestations in sarcoidosis range from asymptomatic disease to severe loss-of-function, including an acute disease (Löfgren’s syndrome), which usually resolves spontaneously, a subacute disease, which also may resolve spontaneously or with treatment, and a chronic/progressive disease. Associations vary from protective to predisposing markers, or markers influencing clinical outcomes. Thus, previous studies have pointed out the importance to characterize patient groups according to clinical phenotypes to avoid inconsistency between studies (6, 9–12). Associations with sarcoidosis also vary between different ancestral groups. Most classical examples being strong association found between HLA-DRB1*03 and Löfgren’s syndrome in caucasian populations, whereas among the Japanese, this association is absent due the lack of HLA-DRB1*03 allele in the population (11).

The most probable pathophysiology of sarcoidosis, the dysregulation of the immune response (13), strongly suggests benefits from a better understanding of the role of the immune-mediating genes in sarcoidosis susceptibility. Based on the robust evidence of the MHC class III gene association to sarcoidosis, this particular region warrants further investigation. Probably, the most investigated variant in the MHC class III, the splice-site SNP in BTNL2 (OMIM 606000) gene (rs2076530 G>A), has so far shown contradictory results with different sarcoidosis phenotypes and populations (7, 14–16).

Here, we present results from a Finnish case-control discovery sample as well as three independent replication studies from the Swedish, Dutch, and Czech populations. Our aim was to investigate genetic variance and phenotype specific variants in five functional candidate genes: lymphotoxin alfa (LTA; OMIM 153440), tumor necrosis factor (TNF; OMIM 191160), advanced glycosylation end product-specific receptor (AGER; OMIM 600214), butyrophilin-like 2 (BTNL2; OMIM 606000), and HLA-DR alpha (HLA-DRA; OMIM 142860) within the MHC classes III and II with tag-SNP genotyping approach. We also address the question whether the tag-SNP associations were due to the strong linkage disequilibrium (LD) with HLA-DRB1. Replication of gene variants predisposing to disease and disease phenotypes in several European populations would indicate their important role in disease development.

Patients and Methods

Patients and Control Subjects

Total of 805 sarcoidosis patients and 870 controls were included in the study. The discovery sample set consisted of 188 Finnish sarcoidosis patients and 150 Finnish healthy controls. Swedish (cases = 219, controls = 360), Dutch (cases = 180, controls = 180), and Czech (cases = 218, controls = 180) data sets were used for the replication and the joint analysis. In quality control, 40 sarcoidosis patients (1 Finnish, 29 Swedish, 10 Czech) and 11 controls (2 Swedish, 2 Czech, 7 Dutch) were excluded due to low quality of genotyping results.

Table 1 shows the sample characteristics. Although we had no age matched controls to the patients, the controls were unrelated healthy individuals recruited from the same geographic and ancestral background as the cases. More detailed sample characteristics and DNA extraction have been previously reported (10, 12, 17–19). Briefly, in the Finnish discovery sample, the patients were recruited from 17 pulmonary units throughout the country. The control population consisted of voluntary Finnish subjects attending a health survey, representing the Finnish population. The Swedish patients were recruited from single center, Karolinska University Hospital, Solna, Sweden. The control group consisted of consecutively collected Scandinavian blood donors. The Czech cases were gathered from the Olomouc University Hospital, Olomouc, Czech Republic, and controls were Czech nationality blood donors from the same region. The Dutch patients with sarcoidosis were from the St. Antonius Hospital, Nieuwegein, Netherlands. The control subjects comprised of healthy, Dutch Caucasian employees of the St. Antonius Hospital, and blood donors from Sanquin blood bank in the Netherlands. In all four study populations, the DNA was extracted from the buffy coat fraction of whole blood samples.

TABLE 1

Table 1. Sample characteristics.

Disease activity was determined using the generally accepted WASOG (World Association of Sarcoidosis and Other Granulomatous diseases) criteria. The sarcoidosis patients had been followed for at least 2 years and further subdivided into those with a resolved (normalized radiography and pulmonary lung function; no signs of active extrapulmonary disease) or a non-resolved disease (persistence of chest radiographic changes with clinical signs of disease activity). The resolution threshold in the Finnish, Swedish, and Czech populations was 2 years but 4 years in the Dutch. The patients were grouped into those with Löfgren’s syndrome [Löfgren syndrome (LS); n = 136] and non-Löfgren’s syndrome patients (NL; n = 629). Dutch cohort did not contain LS patients.

NL patients were grouped into those with a non-active disease [NL resolved (NLR); n = 249] and those with persisting activity at that time point [NL persistent (NLP); n = 337]. Forty-three patients (5.6%) had no subgroup information available, and we excluded them from the subgroup analysis. Due to small sample size, the LS patients were not subgrouped based on the disease activity.

All the subjects were of European descent and gave their written informed consent to participate in the genetic association study. The local Ethics Committee of the Department of Internal Medicine, Hospital District of Helsinki and Uusimaa, Finland approved the study protocol.

Genotyping

We selected MHC gene regions of HLA-DRA, LTA, TNF, AGER, and BTNL2 for SNP analysis and used the SNP tagging approach to identify a SNP or a set of SNPs associating with the trait. The tag-SNPs reduce the number of SNPs needed to cover the selected region and are chosen from the HapMap database. Information about the validation status, tagging quality, minor allele frequency (MAF) (>0.01), and gene structure was used for selecting the SNPs. In addition to HapMap database, the public dbSNP¹ database was used to select additional SNPs. Eventually, 89 SNPs tagging two 100 and 92 SNP variants were genotyped in the Finnish discovery sample [r² > 0.8; HapMap3 (20)]. The detailed information for SNP genotypes is provided in Table S1 in Supplementary Material. In the replication stage, 16 of the 89 SNPs (Table 2) were genotyped in the Swedish, Dutch, and Czech samples (Table S2 in Supplementary Material). Thirteen of the 16 SNPs were selected based on unadjusted P-values of P < 0.05 found in the discovery sample and only one SNP from each LD block (r² < 0.8) were chosen. In addition, three SNPs (rs1800624, rs3130349, rs3135365) were included based on the published studies (14) showing sarcoidosis association.

TABLE 2

Table 2. Disease associated SNPs from the major histocompatibility complex region and their allele frequencies in the Finnish sarcoidosis discovery sample set.

We did assays designing with AssayDesign software and performed the multiplex PCR and the iPLEX reaction using 9–10 ng of DNA as a template. The SNP allele separation based on the differences of the single base extension products used the Sequenom MassArray iPLEX system (Sequenom, San Diego, CA, USA). The HLA-DRB1 genotypes for the samples were previously published (10, 12, 18) including genotypes for 497 cases and 517 controls.

Statistical Analysis

We compared allele frequencies of 89 SNPs of the discovery sample and of 16 SNPs of the replication samples between different groups (NL vs. controls, LS vs. controls, NL resolved vs. NL persistent) by a case-control association analysis (Chi-square χ² test, PLINK software) (21) (Table S2 in Supplementary Material). We applied the following quality control filters: minimum call rate per sample of 90%, SNP MAF >0.01, and Hardy–Weinberg equilibrium >0.001. Total success rate for accepted SNP arrays was 99% in the discovery sample and 98% in the replication samples together.

We combined the results from SNP analyses of the discovery and replication samples for a joint analysis using the random effects model (PLINK software). We used the Benjamini–Hochberg False Discovery Rate method for correction for multiple testing. SNP associations and HLA-DRB1 alleles were adjusted for sex, and all the HLA-DRB1 alleles by a logistic regression [PAWStatistics 18.0, PLINK software (21)].

To determine the genetic relationships among individuals belonging to four different European populations (Finnish, Swedish, Dutch, and Czech), we performed principal coordinates analyses (PCoA), a multivariate technique, using pairwise individual-by-individual linear genetic distance matrix with covariance standardized method in GenAlEx 6.5 tool. This linear genetic distance approach correlates genetic and geographic distances (Mantel test) within the dataset of 16 SNPs in healthy control subjects and sarcoidosis patients across the studied European populations. The method identified the major variation axis within our multidimensional genotype data set. As each successive axis explained proportionately less of the total genetic variation, the first two axes were used to reveal the major separation among individuals. Heterogeneity between the studies was evaluated using the I² metric with a heterogeneity threshold of I² < 25% (low heterogeneity) (22). The haplotypes were constructed and pairwise LD (r²) detected using SNP Haploview software (23). The LD structure in control samples from all study populations was compared to the LD structure of population genotype data originating from Phase 3 of the 1,000 Genomes Project (24). We estimated HLA-DRB1-SNP haplotype frequencies from allele data using the Bayesian method with PHASE v. 2.1.1 (23). The HLA-DRB1 low-resolution data (HLA-DRB1*01, *03, *04, *07 *08, *10, *11, *12, *13, *14, *15, *16) was available for the Finnish (n = 187 cases, n = 150 controls), and for a part of the Swedish (n = 184 cases, n = 187 controls) and Czech (n = 90 cases, n = 180 controls) samples. The HapMap3 proxy SNPs (r² > 0.8) were detected using SNAP software (25, 26). We followed the STREIS principles in immunogenomic data analysis (27).

Analysis of Expression Quantitative Trait Loci (eQTL) Data and Pathway Connectivity Analysis

To further investigate the possibly functional effects of the significant SNPs, we used GENe Expression Variation (GeneVar) (28) database to study the eQTL.

Results

A total of 805 sarcoidosis patients and 870 controls were included in this study. The discovery sample set consisted of Finnish sarcoidosis patients and controls. Swedish, Dutch, and Czech data sets were the replication samples and used for the joint analysis. We analyzed 89 SNPs for the discovery samples and selected 16 of these SNPs for the replication and the joint analysis.

Thirteen SNPs Showed Association with Sarcoidosis in the Discovery Sample

Seven SNPs showed association (uncorrected P < 0.05) with NL in the discovery sample (Table 2). The strongest association with NL was observed for rs3177928 (P = 0.001, OR = 2.17) located within downstream of HLA-DRA. Four SNPs in BTNL2 showed significant association with NL, of these, the most significant rs28362677 (P = 0.002, OR = 1.92) was a missense variant (Ser360Gly). Two SNPs (rs5007259 in BTNL2 promoter and rs6937545 in downstream of HLA-DRA) associated with LS. Six SNPs [rs9268528 (BTNL2 promoter), rs3135351, and rs3129843 (between genes HLA-DRA and BTNL2), rs9268644 and rs3129877 (HLA-DRA intronic), rs6937545 (downstream of HLA-DRA)] were associated with the disease course of sarcoidosis (Table 2). Neither TNF nor LTA variants were significant in the discovery sample or showed strong LD with other MHC variants, hence, not included in the replication stage. For discovery stage, detailed results of the allele frequencies from 89 SNPs are found in Table S1 in Supplementary Material.

Shared and Population-Specific SNP Associations with Sarcoidosis in the Replication Samples

Heterogeneity at association level was observed among the samples and BTNL2 SNPs showed population-specific differences (Table S2 in Supplementary Material). The sarcoidosis-associated BTNL2 missense SNP rs2076530 (16) and BTNL2 promoter SNP rs5007259 were associated in Swedish (rs2076530 with borderline association) and Czech samples (NL vs. controls), but were not replicated in the Dutch sample. In Dutch samples, the most significant exonic BTNL2 SNP rs28362677 (NL 93% vs. controls 83%; P = 0.00016, OR = 2.5) was not replicated among Swedish and Czech samples.

PCoA for Population Stratification

Due to different European origins, the samples were checked for ancestral homogeneity using the PCoA. In all study populations, the healthy control subjects seemed to be more diverse than the sarcoidosis patients. The Finnish and Czech healthy controls were most diverged, while Czech and Swedish sarcoidosis patients were among least diverged. Overall, a homogenous clustering with absence of any major cluster(s) among the individuals and low percent of variations was determined (Figures S1 and S2 in Supplementary Material).

Joint Analysis Combining Results from Discovery and Replication Samples

The allelic associations found in the discovery and replication samples (after filtering, 765 sarcoidosis patients and 859 controls) were used for the joint analysis. Four of the SNPs were replicated in all four populations: rs3763313 (NL vs. C; BTNL2 upstream-variant), rs3177928 (NL vs. C; HLA-DRA downstream variant), rs5007259 (L vs. C; BTNL2 upstream-variant), and rs6937545 (L vs. C; HLA-DRA downstream variant).

In the joint analysis, all the other SNPs (P-value for random effects < 0.05), except rs9268528, were associated with at least one of the disease phenotypes (I² < 25%) (Table S2 in Supplementary Material). Seven SNPs in AGER, BTNL2, and HLA-DRA associated with NL, two most significant being downstream variant of HLA-DRA (rs3177928, P = 0.0000002, OR = 1.9) and BTNL2 promoter SNP (rs3763313, P = 0.000001, OR = 1.6). Eight SNPs associated with LS (Table S2 in Supplementary Material). Three of these SNPs had P < 10⁻⁸: rs3130349 in AGER (OR = 0.39), rs3129843 between genes BTNL2 and HLA-DRA (OR = 3.4), and rs6937545 in HLA-DRA (OR = 2.3). Fifty-seven percent (78/136) of LS patients were Swedish. Five SNPs were associated with the disease course of sarcoidosis (BTNL2 exonic SNP rs28362677, BTNL2 promoter SNP rs5007259, SNPs rs3135351 and rs3129843 between genes BTNL2 and HLA-DRA, and HLA-DRA downstream SNP rs3129877) (Table S2 in Supplementary Material). All the associations, except rs5007259 in NLR vs. NLP, remained significant after correcting for multiple comparisons (Table S2 in Supplementary Material).

Haplotype Analysis Showed Population-Specific Structures

Haplotype analysis was concordant with the single allele findings showing that population-specific structures occur (Figure 1). Haplotype 1 with HLA-DRB1*01 was more common in the Finnish sample (cases 6%, controls 14%) than in the Swedish (cases 3%, controls 8%) or the Czech sample (cases 3%, controls 5%). The haplotype 3 with HLA-DRB1*03 was the most frequent haplotype in the Swedish sample (cases 19%, controls 12%) probably due to the high number of LS patients that commonly share the allele. Haplotype 28 with HLA-DRB1*16 was missing or rare in the Finnish and Swedish sample.

FIGURE 1

Figure 1. HLA-DRB1-SNP haplotype frequencies. HLA-DRB1 data were available for Finnish, Swedish, and Czech samples. Subjects with missing SNP genotype were excluded from the analysis. SNP order is the same as in the Table 2. The analysis shows population-specific differences in haplotype structures. Haplotype frequencies (2n) > 0.01 are presented.

HLA-DRB1-Independent Association in Four SNPs

Addressing the question of whether the SNP associations were secondary to the HLA-DRB1 associations found previously (10, 12, 18, 29), we performed logistic regression adjusting for the low-resolution HLA-DRB1 alleles (Table 3). In the comparisons of “NL vs. controls (recessive model),” SNPs rs3135365 and rs3177928, HLA-DRB1*15 and *16 were suggested as independent variants (P = 0.015, OR = 2.15, CI 95% = 1.16–3.99; P < 0.001, OR = 2.19, CI 95% = 1.54–3.15; P < 0.001, OR = 3.49, CI 95% = 1.87–6.53, P = 0.044, OR = 0.43, CI 95% = 0.19–0.98; respectively). For LS using a dominant model, the downstream SNP of HLA-DRA (rs6937545; P = 0.012, OR = 3.49, CI 95% = 1.32–9.25) and HLA-DRB1*03 (P < 0.001, OR = 4.67, CI 95% = 2.18–10.12), *13 (P = 0.034, OR = 2.41, CI 95% = 1.07–5.42), and *14 (P = 0.007, OR = 5.38, CI 95% = 1.59–18.2) alleles were suggested as significant predisposing markers. In the analysis of NLR vs. NLP (recessive model), no independent associations were detected in studied SNPs. None of the AGER SNPs remained significant after adjusting for HLA-DRB1 alleles showing strong LD with HLA-DRB1.

TABLE 3

Table 3. Disease associated SNPs in the major histocompatibility complex region in the joint analysis (Finnish, Swedish, Dutch, and Czech) and associations after adjusting for HLA-DRB1 low-resolution alleles.

LD Structure Confirmed SNP Selection Strategy

Figures 2–4 visualize genotyped SNPs (P < 0.05 and I² < 25%) and the LD structure (r²) based on HapMap data for both cases and controls. Estimated recombination rates are plotted to show the LD structure around the associated SNPs. Gene annotations were adapted from the University of California at Santa Cruz Genome Browser.² As none of the genotyped SNPs had strong (r² > 0.8) pairwise LD in this study, it confirmed the selection strategy where only one SNP from each LD block was genotyped. LD structure of control populations corresponded to the 1,000 Genomes population data (data not shown).

FIGURE 2

Figure 2. The meta-analysis of non-Löfgren sarcoidosis (NL) vs. controls with Finnish, Swedish, Dutch, and Czech samples. The linkage disequilibrium (LD) information (r²) is shown by color-coded dots and the recombination rates where the peaks are the hotspots are based on the HapMap data. The LD structure (r²) in the current study for the top significant SNPs (I² < 25) are shown on the right. After adjusting for HLA-DRB1 alleles, the SNPs that remained significant are pointed with an arrow.

FIGURE 3

Figure 3. The meta-analysis of Löfgren (LS) vs. controls with Finnish, Swedish, and Czech samples. The linkage disequilibrium (LD) information (r²) is shown by color-coded dots and the recombination rates where the peaks are the hotspots are based on the HapMap data. The LD structure (r²) in the current study for the top significant SNPs (I² < 25) are shown on the right. After adjusting for HLA-DRB1 alleles, the SNPs that remained significant are pointed with an arrow.

FIGURE 4

Figure 4. The meta-analysis of sarcoidosis subgroups (resolved sarcoidosis, NLR vs. persistent sarcoidosis, NLP) with Finnish, Swedish, Dutch, and Czech samples. The linkage disequilibrium (LD) information (r²) is shown by color-coded dots and the recombination rates where the peaks are the hotspots are based on the HapMap data. The LD structure (r²) in the current study for the top significant SNPs (I² < 25) are shown on the right. After adjusting for HLA-DRB1 alleles, the SNPs that remained significant are pointed with an arrow.

Four SNPs Act As Local eQTL for HLA-DRB1/-DRB5

Using the GeneVar database (28), we analyzed the available eQTL-SNP-gene data to identify eQTL genes surrounding independent SNPs [lymphoblastoid cell lines, CEU (30)]. The analysis suggested SNPs rs3135365, rs3177928, and rs6937545 act as cis-acting eQTL for HLA-DRB1, and rs3135365, rs6937545, and rs5007259 act as eQTL for HLA-DRB5. Figure S3 in Supplementary Material presents the SNP-probe association plots.

Discussion

Major histocompatibility complex region contains immune mediator genes that have showed either predisposing or protective effects to sarcoidosis (12, 15). Here, we report a set of 16 SNPs located in the MHCs that are associated with sarcoidosis and shared between four European populations with distinct ancestral origins (Finnish, Swedish, Dutch, and Czech). Only two of these SNPs have been previously associated with sarcoidosis. After conditioning with HLA-DRB1 alleles, three of the 16 SNPs remained independently associated with sarcoidosis.

We identified two novel HLA-DRA downstream variants that were independent of HLA-DRB1 alleles: rs3177928 associated with NL and rs6937545 associated with LS. In a high-density genetic mapping study of extended MHC region in four European populations and one black African descent population, NL sarcoidosis showed similar association pattern: main associations were found in variants located in the MHC class II region (31). MHC II region associations were found in LS patients as well. Recent non-sarcoidosis studies showed that rs3177928 was associated with lipoprotein metabolism (32) and connected with inflammatory mechanisms (33–35). HLA-DRA encodes the α-subunit of the HLA-DR cell surface receptor and expresses on antigen-presenting cells. HLA-DRA shows strong LD structure and its variants regulate the expression of other MHC genes (36). Both rs3177928 and rs6937545 are also cis-acting eQTL affecting the expression of HLA-DRB1 and are not in LD with other HLA-DRA gene variants previously associated with sarcoidosis (6–8). Interestingly, HLA-DRA has shown to have SNP–SNP (8) interaction with ANXA11, another sarcoidosis associated gene. ANXA11 has an essential role in cell division and apoptosis (6, 7).

The increasing line of evidence suggests that BTNL2, a member of the immunoglobulin gene superfamily, is a key element for sarcoidosis predisposition (6, 15, 16, 37). BTNL2 regulates T cells and is expressed on dendritic cells (15, 38). Here, the novel variant in BTNL2 5’ region associated independently with distinct sarcoidosis phenotype (rs3135365 with NL) supporting the importance of the BTNL2 promoter region in sarcoidosis disease development. Several other BTNL2 SNPs also reached the level of significance. However, they showed either heterogeneity between populations (e.g., exonic SNP rs28362677 and 5’ upstream SNP rs5007259) or the association weakened after adjustments with HLA-DRB1, as occurred with the splice-site variant of BTNL2 (rs2076530) that was previously associated with chronic sarcoidosis in subjects of European descent (6, 16, 30, 37, 39–41). The strong LD between BTNL2 and HLA-DRB1 complicates the interpretation of the results and the inconsistency with previous studies may be due to differences in the study design (heterogeneity in diagnosis) and population origin (16, 42, 43). Given that, we suggest that district variants of BTNL2 are present only in certain populations, like exonic SNP rs28362677 in the Dutch sample that may explain the lack of replication of rs2076530 (11, 40). Fine-mapping studies of BTNL2 in large multi-population sample would shed light on the allelic heterogeneity of BTNL2.

HLA allele variation in Europe follows the North to Southeast axis and as expected, different HLA haplotype profiles were observed in this study (44). We noticed that the use of different population data was advantageous for the HLA-DRB1 adjustments, because the populations have different chromosomal regions where LD breaks down as well as distinct allele frequencies (4). HLA-DRB1*03 has been associated with sarcoidosis with favorable prognosis and Löfgren’s syndrome in several populations with European descent, especially in Scandinavian countries (10, 39, 45–47). In Japan, the HLA-DRB1*03 allele is rare, explaining the non-association of the variant with sarcoidosis (48). We showed that LS-associated SNPs were inherited as a conserved block with the HLA-DRB1*03, especially in the Swedish sample (10), and after conditioning with HLA-DRB1 alleles, only one SNP remained significant (rs6937545 in HLA-DRA downstream). Given the strong LD within MHC genes, we highlight the importance to use of HLA genotypes in addition to SNP data when aiming to pinpoint the causal variants within MHC region.

There are several strengths in our study. The sample size in the joint analysis of the discovery and replication sets was relatively large for case-control study in the context of sarcoidosis and its incidence. Our fine-mapping strategy of MHC candidate genes in different ancestral origins aimed to detect variants associated with different disease phenotypes shared in different populations. Previous studies have pointed out that the precise clinical characterization of the patients is essential, because sarcoidosis is a highly heterogeneous disease. To overcome the disease heterogeneity, we subdivided the patients according to clinical phenotypes (49, 50).

In conclusion, we found novel SNPs in BTNL2 and HLA-DRA regions associating with sarcoidosis. Our finding further establishes that polymorphisms in the HLA-DRA and BTNL2 have a role in sarcoidosis susceptibility. This multi-population study demonstrates that at least a part of these associations are HLA-DRB1 independent (e.g., not due to LD), and shared across ancestral origins. The variants that were independent of HLA-DRB1 associations, acted as eQTL for HLA-DRB1 and/or -DRB5, suggesting a role in regulating gene expression. Future functional studies with larger sample are required to reveal the causal regulatory variation at this locus and the immunogenetic basis related to sarcoidosis.

Ethics Statement

The study was performed with approval of institutional ethical committees at respective centers (Ethics Commitee of the Department of Internal Medicine, Hospital District of Helsinki and Uusimaa, Finland; Ethics Committee of the University Hospital and Medical Faculty of Palacky University, Olomouc, Czech Republic; Ethics Committee of the Karolinska University Hospital, Solna, Sweden; Ethics Committee the University Medical Center Utrecht, Netherlands).

Author Contributions

M-LL, OS, AP, MP, AW, and EL conceived and designed the work. JG, CM, MP, AE, JG, VK, FM, LP, AP, MR, MS, M-LL, and OS contributed to data acquisition. AW and EL gathered SNP information and did all the analyses and interpretation. AK and VA assisted with data analyzing. AW and EL drafted the manuscript after its revision for important intellectual context by OS, AP, MP, VA, FM, AK, M-LL, CM, AE, JG, and MR. AW and EL finalized the article. All authors have read and approved the final manuscript and agreed to be accountable for all aspects of the work.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank Leena Saraste for assisting with language issues.

Funding

This study was supported by Nummela Foundation for Medical Research, Hengityssairauksien tutkimussäätiö (HES), the Helsinki Biomedical Graduate School, the League of European Research Universities (LERU), IGA_LF_2016_009/2017_014 and LO1304, the Swedish Heart-Lung Foundation, the Swedish Research Council, through the regional agreement on medical training and clinical research (ALF) between Stockholm County Council and Karolinska Institutet, and Karolinska Institutet. ESF (Europe) COST of Action BM0803.

Supplementary Material

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fimmu.2017.00422/full#supplementary-material.

Footnotes

1. ^http://www.ncbi.nlm.nih.gov/projects/SNP.
2. ^http://genome.ucsc.edu/.

References

1. Iannuzzi MC, Rybicki BA, Teirstein AS. Sarcoidosis. N Engl J Med (2007) 357:2153–65. doi:10.1056/NEJMra071714

CrossRef Full Text | Google Scholar

2. Sverrild A, Backer V, Kyvik KO, Kaprio J, Milman N, Svendsen CB, et al. Heredity in sarcoidosis: a registry-based twin study. Thorax (2008) 63:894–6. doi:10.1136/thx.2007.094060

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Newman LS, Rose CS, Bresnitz EA, Rossman MD, Barnard J, Frederick M, et al. A case control etiologic study of sarcoidosis: environmental and occupational risk factors. Am J Respir Crit Care Med (2004) 170:1324–30. doi:10.1164/rccm.200402-249OC

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Spagnolo P, Grunewald J. Recent advances in the genetics of sarcoidosis. J Med Genet (2013) 50:290–7. doi:10.1136/jmedgenet-2013-101532

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Rybicki BA, Iannuzzi MC, Frederick MM, Thompson BW, Rossman MD, Bresnitz EA, et al. Familial aggregation of sarcoidosis. A case-control etiologic study of sarcoidosis (ACCESS). Am J Respir Crit Care Med (2001) 164:2085–91. doi:10.1164/ajrccm.164.11.2106001

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Adrianto I, Lin CP, Hale JJ, Levin AM, Datta I, Parker R, et al. Genome-wide association study of African and European Americans implicates multiple shared and ethnic specific loci in sarcoidosis susceptibility. PLoS One (2012) 7:e43907. doi:10.1371/journal.pone.0043907

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hofmann S, Franke A, Fischer A, Jacobs G, Nothnagel M, Gaede KI, et al. Genome-wide association study identifies ANXA11 as a new susceptibility locus for sarcoidosis. Nat Genet (2008) 40:1103–6. doi:10.1038/ng.198

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Levin AM, Iannuzzi MC, Montgomery CG, Trudeau S, Datta I, McKeigue P, et al. Association of ANXA11 genetic variation with sarcoidosis in African Americans and European Americans. Genes Immun (2013) 14:13–8. doi:10.1038/gene.2012.48

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Rossman MD, Thompson B, Frederick M, Maliarik M, Iannuzzi MC, Rybicki BA, et al. HLA-DRB1*1101: a significant risk factor for sarcoidosis in blacks and whites. Am J Hum Genet (2003) 73:720–35. doi:10.1086/378097

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Grunewald J, Brynedal B, Darlington P, Nisell M, Cederlund K, Hillert J, et al. Different HLA-DRB1 allele distributions in distinct clinical subgroups of sarcoidosis patients. Respir Res (2010) 11:25. doi:10.1186/1465-9921-11-25

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sato H, Woodhead FA, Ahmad T, Grutters JC, Spagnolo P, van den Bosch JM, et al. Sarcoidosis HLA class II genotyping distinguishes differences of clinical phenotype across ethnic groups. Hum Mol Genet (2010) 19:4100–11. doi:10.1093/hmg/ddq325

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Wennerström A, Pietinalho A, Vauhkonen H, Lahtela L, Palikhe A, Hedman J, et al. HLA-DRB1 allele frequencies and C4 copy number variation in Finnish sarcoidosis patients and associations with disease prognosis. Hum Immunol (2012) 73:93–100. doi:10.1016/j.humimm.2011.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Spagnolo P, Richeldi L, du Bois RM. Environmental triggers and susceptibility factors in idiopathic granulomatous diseases. Semin Respir Crit Care Med (2008) 29:610–9. doi:10.1055/s-0028-1101271

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Campo I, Morbini P, Zorzetto M, Tinelli C, Brunetta E, Villa C, et al. Expression of receptor for advanced glycation end products in sarcoid granulomas. Am J Respir Crit Care Med (2007) 175:498–506. doi:10.1164/rccm.200601-136OC

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Nguyen T, Liu XK, Zhang Y, Dong C. BTNL2, a butyrophilin-like molecule that functions to inhibit T cell activation. J Immunol (2006) 176:7354–60. doi:10.4049/jimmunol.176.12.7354

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Valentonyte R, Hampe J, Huse K, Rosenstiel P, Albrecht M, Stenzel A, et al. Sarcoidosis is associated with a truncating splice site mutation in BTNL2. Nat Genet (2005) 37:357–64. doi:10.1038/ng1519

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Heron M, van Moorsel CH, Grutters JC, Huizinga TW, van der Helm-van Mil AH, Nagtegaal MM, et al. Genetic variation in GREM1 is a risk factor for fibrosis in pulmonary sarcoidosis. Tissue Antigens (2011) 77:112–7. doi:10.1111/j.1399-0039.2010.01590.x

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Foley PJ, McGrath DS, Puscinska E, Petrek M, Kolek V, Drabek J, et al. Human leukocyte antigen-DRB1 position 11 residues are a common protective marker for sarcoidosis. Am J Respir Cell Mol Biol (2001) 25:272–7. doi:10.1165/ajrcmb.25.3.4261

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Veltkamp M, van Moorsel CH, Rijkers GT, Ruven HJ, Grutters JC. Genetic variation in the toll-like receptor gene cluster (TLR10-TLR1-TLR6) influences disease course in sarcoidosis. Tissue Antigens (2012) 79:25–32. doi:10.1111/j.1399-0039.2011.01808.x

PubMed Abstract | CrossRef Full Text | Google Scholar

20. International HapMap 3 Consortium, Altshuler DM, Gibbs RA, Peltonen L, Altshuler DM, Gibbs RA, et al. Integrating common and rare genetic variation in diverse human populations. Nature (2010) 467:52–8. doi:10.1038/nature09298

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet (2007) 81:559–75. doi:10.1086/519795

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Higgins JP, Thompson SG. Quantifying heterogeneity in a meta-analysis. Stat Med (2002) 21:1539–58. doi:10.1002/sim.1186

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet (2001) 68:978–89. doi:10.1086/319501

CrossRef Full Text | Google Scholar

24. Machiela MJ, Chanock SJ. LDlink: a web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants. Bioinformatics (2015) 31(21):3555–7. doi:10.1093/bioinformatics/btv402

CrossRef Full Text | Google Scholar

25. Johnson AD, Handsaker RE, Pulit SL, Nizzari MM, O’Donnell CJ, de Bakker PI. SNAP: a web-based tool for identification and annotation of proxy SNPs using HapMap. Bioinformatics (2008) 24:2938–9. doi:10.1093/bioinformatics/btn564

PubMed Abstract | CrossRef Full Text | Google Scholar

26. International HapMap Consortium. The international HapMap project. Nature (2003) 426:789–96. doi:10.1038/nature02168

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Hollenbach JA, Mack SJ, Gourraud PA, Single RM, Maiers M, Middleton D, et al. A community standard for immunogenomic data reporting and analysis: proposal for a STrengthening the REporting of immunogenomic studies statement. Tissue Antigens (2011) 78:333–44. doi:10.1111/j.1399-0039.2011.01777.x

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yang TP, Beazley C, Montgomery SB, Dimas AS, Gutierrez-Arcelus M, Stranger BE, et al. Genevar: a database and java application for the analysis and visualization of SNP-gene associations in eQTL studies. Bioinformatics (2010) 26:2474–6. doi:10.1093/bioinformatics/btq452

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Voorter CE, Drent M, van den Berg-Loonen EM. Severe pulmonary sarcoidosis is strongly associated with the haplotype HLA-DQB1*0602-DRB1*150101. Hum Immunol (2005) 66:826–35. doi:10.1016/j.humimm.2005.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Stranger BE, Nica AC, Forrest MS, Dimas A, Bird CP, Beazley C, et al. Population genomics of human gene expression. Nat Genet (2007) 39:1217–24. doi:10.1038/ng2142

CrossRef Full Text | Google Scholar

31. Rivera NV, Ronninger M, Shchetynsky K, Franke A, Nöthen MM, Müller-Quernheim J, et al. High-density genetic mapping identifies new susceptibility variants in sarcoidosis phenotypes and shows genomic-driven phenotypic differences. Am J Respir Crit Care Med (2016) 193(9):1008–22. doi:10.1164/rccm.201507-1372OC

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Teslovich TM, Musunuru K, Smith AV, Edmondson AC, Stylianou IM, Koseki M, et al. Biological, clinical and population relevance of 95 loci for blood lipids. Nature (2010) 466:707–13. doi:10.1038/nature09270

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Ivanišević J, Kotur-Stevuljević J, Stefanović A, Jelić-Ivanović Z, Spasić S, Videnović-Ivanov J, et al. Dyslipidemia and oxidative stress in sarcoidosis patients. Clin Biochem (2012) 45:677–82. doi:10.1016/j.clinbiochem.2012.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Salazar A, Mañá J, Pintó X, Argimón JM, Hurtado I, Pujol R. Corticosteroid therapy increases HDL-cholesterol concentrations in patients with active sarcoidosis and hypoalphalipoproteinemia. Clin Chim Acta (2002) 320:59–64. doi:10.1016/S0009-8981(02)00046-3

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Salazar A, Pinto X, Mana J. Serum amyloid A and high-density lipoprotein cholesterol: serum markers of inflammation in sarcoidosis and other systemic disorders. Eur J Clin Invest (2001) 31:1070–7. doi:10.1046/j.1365-2362.2001.00913.x

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Stranger BE, Montgomery SB, Dimas AS, Parts L, Stegle O, Ingle CE, et al. Patterns of cis regulatory variation in diverse human populations. PLoS Genet (2012) 8:e1002639. doi:10.1371/journal.pgen.1002639

CrossRef Full Text | Google Scholar

37. Wijnen PA, Voorter CE, Nelemans PJ, Verschakelen JA, Bekers O, Drent M. Butyrophilin-like 2 in pulmonary sarcoidosis: a factor for susceptibility and progression? Hum Immunol (2011) 72:342–7. doi:10.1016/j.humimm.2011.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Arnett HA, Escobar SS, Gonzalez-Suarez E, Budelsky AL, Steffen LA, Boiani N, et al. BTNL2, a butyrophilin/B7-like molecule, is a negative costimulatory molecule modulated in intestinal inflammation. J Immunol (2007) 178:1523–33. doi:10.4049/jimmunol.178.3.1523

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Spagnolo P, Sato H, Grutters JC, Renzoni EA, Marshall SE, Ruven HJ, et al. Analysis of BTNL2 genetic polymorphisms in British and Dutch patients with sarcoidosis. Tissue Antigens (2007) 70:219–27. doi:10.1111/j.1399-0039.2007.00879.x

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Rybicki BA, Walewski JL, Maliarik MJ, Kian H, Iannuzzi MC; ACCESS Research Group. The BTNL2 gene and sarcoidosis susceptibility in African Americans and Whites. Am J Hum Genet (2005) 77:491–9. doi:10.1086/444435

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Milman N, Svendsen CB, Nielsen FC, van Overeem Hansen T. The BTNL2 A allele variant is frequent in Danish patients with sarcoidosis. Clin Respir J (2011) 5:105–11. doi:10.1111/j.1752-699X.2010.00206.x

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Sato H, Spagnolo P, Silveira L, Welsh KI, du Bois RM, Newman LS, et al. BTNL2 allele associations with chronic beryllium disease in HLA-DPB1*Glu69-negative individuals. Tissue Antigens (2007) 70:480–6. doi:10.1111/j.1399-0039.2007.00944.x

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Coudurier M, Freymond N, Aissaoui S, Calender A, Pacheco Y, Devouassoux G. Homozygous variant rs2076530 of BTNL2 and familial sarcoidosis. Sarcoidosis Vasc Diffuse Lung Dis (2009) 26:162–6.

PubMed Abstract | Google Scholar

44. Riccio ME, Buhler S, Nunes JM, Vangenot C, Cuénod M, Currat M, et al. 16(th) IHIW: analysis of HLA population data, with updated results for 1996 to 2012 workshop data (AHPD project report). Int J Immunogenet (2013) 40:21–30. doi:10.1111/iji.12033

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Morais A, Alves H, Lima B, Delgado L, Gonçalves R, Tafulo S. HLA class I and II and TNF-alpha gene polymorphisms in sarcoidosis patients. Rev Port Pneumol (2008) 14:727–46. doi:10.1016/S0873-2159(15)30284-1

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Bogunia-Kubik K, Tomeczko J, Suchnicki K, Lange A. HLA-DRB1*03, DRB1*11 or DRB1*12 and their respective DRB3 specificities in clinical variants of sarcoidosis. Tissue Antigens (2001) 57:87–90. doi:10.1034/j.1399-0039.2001.057001087.x

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Martinetti M, Tinelli C, Kolek V, Cuccia M, Salvaneschi L, Pasturenzi L, et al. “The sarcoidosis map”: a joint survey of clinical and immunogenetic findings in two European countries. Am J Respir Crit Care Med (1995) 152:557–64. doi:10.1164/ajrccm.152.2.7633707

CrossRef Full Text | Google Scholar

48. Morimoto T, Azuma A, Abe S, Usuki J, Kudoh S, Sugisaki K, et al. Epidemiology of sarcoidosis in Japan. Eur Respir J (2008) 31:372–9. doi:10.1183/09031936.00075307

CrossRef Full Text | Google Scholar

49. Berlin M, Fogdell-Hahn A, Olerup O, Eklund A, Grunewald J. HLA-DR predicts the prognosis in Scandinavian patients with pulmonary sarcoidosis. Am J Respir Crit Care Med (1997) 156:1601–5. doi:10.1164/ajrccm.156.5.9704069

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Spagnolo P, Sato H, Grunewald J, Brynedal B, Hillert J, Mañá J, et al. A common haplotype of the C-C chemokine receptor 2 gene and HLA-DRB1*0301 are independent genetic risk factors for Lofgren’s syndrome. J Intern Med (2008) 264:433–41. doi:10.1111/j.1365-2796.2008.01984.x

CrossRef Full Text | Google Scholar