Functional link between sarcoidosis-associated gene variants and quantitative levels of bronchoalveolar lavage fluid cell types
Muntasir Abo Al Hayja1† Susanna Kullberg1† Anders Eklund1 
Leonid Padyukov2,3 Johan Grunewald1,3 
Natalia V. Rivera1,2,3*- 1Division of Respiratory Medicine, Department of Medicine Solna, Karolinska Institutet, Karolinska University Hospital, Stockholm, Sweden
- 2Division of Rheumatology, Department of Medicine Solna, Karolinska Institutet, Karolinska University Hospital, Stockholm, Sweden
- 3Center of Molecular Medicine (CMM), Karolinska Institutet, Stockholm, Sweden
Background: Sarcoidosis is an inflammatory disease that affects multiple organs. Cell analysis from bronchoalveolar lavage fluid (BALF) is a valuable tool in the diagnostic workup and differential diagnosis of sarcoidosis. Besides the expansion of lymphocyte expression-specific receptor segments (Vα2.3 and Vβ22) in some patients with certain HLA types, the relation between sarcoidosis susceptibility and BAL cell populations’ quantitative levels is not well-understood.
Methods: Quantitative levels defined by cell concentrations of BAL cells and CD4+/CD8+ ratio were evaluated together with genetic variants associated with sarcoidosis in 692 patients with extensive clinical data. Genetic variants associated with clinical phenotypes, Löfgren’s syndrome (LS) and non-Löfgren’s syndrome (non-LS), were examined separately. An association test via linear regression using an additive model adjusted for sex, age, and correlated cell type was applied. To infer the biological function of genetic associations, enrichment analysis of expression quantitative trait (eQTLs) across publicly available eQTL databases was conducted.
Results: Multiple genetic variants associated with sarcoidosis were significantly associated with quantitative levels of BAL cells. Specifically, LS genetic variants, mainly from the HLA locus, were associated with quantitative levels of BAL macrophages, lymphocytes, CD3+ cells, CD4+ cells, CD8+ cells, CD4+/CD8+ ratio, neutrophils, basophils, and eosinophils. Non-LS genetic variants were associated with quantitative levels of BAL macrophages, CD8+ cells, basophils, and eosinophils. eQTL enrichment revealed an influence of sarcoidosis-associated SNPs and regulation of gene expression in the lung, blood, and immune cells.
Conclusion: Genetic variants associated with sarcoidosis are likely to modulate quantitative levels of BAL cell types and may regulate gene expression in immune cell populations. Thus, the role of sarcoidosis-associated gene-variants may be to influence cellular phenotypes underlying the disease immunopathology.
1. Introduction
Sarcoidosis is a systemic inflammatory disease of unknown etiology. It is postulated that its cause might result from the interplay between genetic and environmental factors that contribute to disease pathogenesis and the formation of noncaseating granulomas (1).
Sarcoidosis is characterized by at least two distinct clinical phenotypes, Löfgren’s syndrome (LS) and non-Löfgren’s syndrome (non-LS) (2, 3). LS is characterized by acute onset, erythema nodosum and/or bilateral periarticular ankle swelling, and enlarged hilar lymph nodes with a high likelihood of spontaneous remission (4). Patients with non-LS, on the other hand, often present with a more gradual, insidious onset and a heterogeneous disease course (3).
Bronchoalveolar lavage (BAL) fluid has a valuable role in the diagnosis of sarcoidosis. BAL cells comprise several cell subtypes, including macrophages, lymphocytes, and neutrophils, as well as small quantities of basophils and eosinophils (5).
In the appropriate clinical and radiographic findings that are compatible with sarcoidosis, an increased total BAL cell counts, CD4+ T cells, and CD4+/CD8+ ratio (> 3.5) strongly supports the diagnosis of the disease (6).
Recent literature suggests sarcoidosis may be an autoimmune disease (7–9). Immunological studies revealed various types of immune cells implicated in sarcoidosis’s pathogenesis. For instance, it is suggested that a persistent T-helper 1 (Th1) immune response is mediated by antigen-presenting cells (APC), resulting in granuloma formation (10, 11). The role of macrophages and CD4+ T cell subtypes has been shown in the immunological profiling of sarcoidosis. Additionally, low levels of IFN-gamma and higher levels of T-bet and ROR-gamma T (12–15) have been reported in LS and non-LS. Nevertheless, a complete understanding of the implication of immune cell populations in the disease is missing.
Genetic investigations have suggested that the genetic architecture of sarcoidosis shares similarities with other autoimmune disorders (7, 16). In particular, the Human Leukocyte Antigen (HLA) region polymorphism plays a significant role in sarcoidosis, especially the HLA-DRB1 gene polymorphism. Moreover, recent evidence in the field showed that non-HLA gene polymorphism (as reviewed in (17)) also plays a significant role in defining sarcoidosis susceptibility. Our group has revealed that sarcoidosis clinical phenotypes LS and non-LS showed different genetic architectures but shared genomic loci in the major histocompatibility complex (MHC) class II region (18).
Although the biological function of genetic variants associated with sarcoidosis is yet to be elucidated, studies on gene expression suggest that genetic variants associated with the disease play a significant role in gene expression signatures via expression quantitative trait loci (eQTLs) (19). Studies on eQTL mapping showed that genetic variation could modify gene expression via cis-eQTL regulatory variants, resulting in specific transcriptomic signatures (19, 20). Genetic variants associated with the disease often act as regulatory switches that regulate gene expression in specific compartments leading to the diversification of disease phenotypes. It is postulated that cis-eQTL genetic variants explain a large proportion of the so-called missing heritability (21).
Since immune and inflammatory functions related to sarcoidosis development are presented by multiple immune cells, it is essential to evaluate regulatory mechanisms for controlling cell representation in relevant tissues. BALF is an essential proxy for the representation of immune cells in the lung and is often used in diagnostic procedures for sarcoidosis. The influence of sarcoidosis-associated genetic variants on quantitative levels of BAL cells has not been explored. Therefore, in the present study, we sought to investigate the relationship between sarcoidosis-associated genetic variants and quantitative levels of BALF cell populations in sarcoidosis clinical phenotypes LS and non-LS.
2. Materials and methods
2.1. Study design and participants
In this study, all included sarcoidosis patients were referred to the Department of Respiratory Medicine, Karolinska University Hospital, Stockholm, Sweden, for an investigational workup of clinically suspected pulmonary sarcoidosis (Table 1). Sarcoidosis diagnosis was based on clinical and radiographical presentation, the absence of an alternative diagnosis, and following the World Association of Sarcoidosis and other Granulomatous Disorders (WASOG) recommendations (22). BAL was performed at the time of the diagnosis. Scadding staging system and spirometry, including diffusion (DLCO), were conducted in most patients. Peripheral blood (PB) sampling for HLA-typing was performed for all patients. Patients signed informed consent, and the study was approved by the Regional Ethical Review Board (Stockholm, Sweden) according to the Declaration of Helsinki.
Table 1. Clinical and laboratory characteristics of studied Löfgren’s syndrome (LS) and non-Löfgren’s syndrome (non-LS) (n = 692).
2.2. Differential BAL cells data
As previously described, a bronchoalveolar lavage procedure was performed in all patients as part of the diagnosis protocol (5). Briefly, a flexible fiber-optic bronchoscope was passed trans-nasally or trans-orally after light sedation and local anesthesia. BAL was performed using sterile phosphate-buffered saline (PBS) solution at 37°C was instilled sequentially. Each aliquot was gradually sucked and collected in a siliconized plastic vessel kept on ice at 4°C. The BAL fluid was gently passed through a Dacron gauze (Millipore, Cork, Ireland) and centrifuged at 400 g for 10 min at 4°C. The pellet was resuspended in PBS. The cells were counted in a Bürker chamber. Cytospins of BAL cells were stained with May-Grünwald Giemsa, and differential cell counts were performed (5). Briefly, staining of centrifugal smears was performed according to May-Grunewald Giemsa. Smears were then fixed in methanol for (5 min), 18 ml Methanol, and 33 ml MAY-Grunewald (12 min). Smears were then rinsed with water 2.5 ml Giemsa solution and 48 ml Phosphate buffer 60 mM PH 6.8 (20 min). Smears were then rinsed with water and let dry. May-Grunewald (Mg-500) Sigma-Aldrich and Giemsa (GS500) Sigma-Aldrich were used. BALF cell analysis was defined by the percentage of cells and quantitative levels based on cell concentrations (× 106/L). BALF analysis was performed for macrophages, lymphocytes (CD3+, CD4+, and CD8+ cells), neutrophils, basophils, and eosinophils. The CD4+/CD8+ ratio was calculated for diagnostic purposes and included as a variable in the study.
2.3. Cohort description and BAL cell populations
A Swedish cohort of 692 sarcoidosis cases (288 LS and 404 non-LS) with BALF cell data was examined. Comparisons of distributions of BAL cell concentrations (× 106/L) for macrophages, lymphocytes, CD3+, CD4+, CD8+ cells, neutrophils, basophils, and eosinophils were conducted using the Wilcoxon Rank Sum test function in R software. The same test was applied for the CD4+/CD8+ ratio. Cohort characteristics are in Table 1.
2.4. Statistical analysis
Single nucleotide polymorphisms (SNPs) associated with LS (n = 1,933) and non-LS (n = 173) published in Rivera et al. (18) were assessed. SNPs were extracted from the ImmunoChip SNP array.
The distribution of BAL cell concentrations (× 106/L) was assessed using the ‘ggqplot’ function in the ggplot2 R package and two normality tests Shapiro–Wilk normality test function (23) and the Anderson-Darling normality test in the nortest R package was applied as secondary assessment. Normalization of BAL cell concentration levels (× 106/L) was applied using a log10-transform.
Outlier detection was conducted using the Grubb test function in the outliers R package. Outliers with a value of p of ≤ 0.05 were removed from the dataset. After removing outliers, normality tests were applied for a final assessment.
Association was assessed on normally distributed log10-transformed BAL cell populations and CD4+/CD8+ ratio in the LS and non-LS groups using an additive model and multivariate linear regression model. The regression coefficient defined by the parameter Beta (or genetic effect) signifies the degree of change in the outcome per unit of change in the quantitative level of the cell type tested. The standard error (SE) signifies the standard deviation of the error in the regression model. Association analysis was performed using PLINK version 1.9 beta (24). The log10-transformed BAL cell subtype was regressed on sex, age, and correlated BAL cell types.
Correlation matrix analysis was conducted using Spearman’s rank correlation method and applied to log10-transformed BAL cell differential counts using the ‘rcorr’ function as implemented in the Hmisc R package. The correlation matrix was visualized as a correlogram using the ‘corrplot’ function in the corrplot R package.
A nominal value of p of < 0.05 and a false discovery rate (FDR) < 0.05 (to account for multiple testing) were considered significant. FDR was calculated using the ‘p.adjust’ function in the R package. To account for high linkage disequilibrium (LD; r2 ≥ 0.8) among associated SNPs, LD pruning was conducted using PriorityPruner version 0.01.4 software1 with default parameters, which identified and prioritized tag-SNPs.
2.5. eQTL enrichment analysis
To explore the biology of genetic variants associated with BAL cell concentrations, a gene-based cross-tissue expression quantitative trait (eQTLs) enrichment was conducted using Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA GWAS) (25). FUMA eQTL databases included the BIOS QTL browser (26), DICE (27), eQTLcatalogue2, eQTLGen (28), GTExt version 83, and scRNA_eQTLs (29) publicly available sources. The selection of eQTL sets was focused on lung, whole blood, and immune cells eQTLs as relevant to the immunopathology of sarcoidosis.
3. Results
3.1. Clinical characteristics
Clinical characterization of the cohort, including age, sex, chest X-ray according to Scadding, BAL cell analysis, BAL differential cell counts in percentages and concentration levels, lung function, and HLA-DRB1 genotypes, was performed in LS and non-LS groups, respectively (Table 1).
Normality assessment on distributions of BAL cell concentrations (× 106/L), log10-transformed BAL cell concentrations, and outlier characterization in LS and non-LS groups are shown in Supplementary Figures S1, S2 in Supplementary Info.
Results from correlation matrix analysis using log10-transformed BAL cells quantitative levels and the Spearman correlation method in LS and non-LS groups showed lymphocytes, macrophages, and neutrophils were correlated (r2 ≥ 0.90 and p < 0.05). Similarly, CD3+, CD4+, and CD8+ cell types were correlated (r2 ≥ 0.99 and p < 0.05). Correlograms are shown in Supplementary Figures S3, S4. Linear regression models (illustrated in Supplementary Figure S5) facilitated the identification of several SNP associations with quantitative levels of BAL cells in LS (Table 2) and non-LS (Table 3).
Table 2. Association results of LS BAL cell types (log10-tranformed of quantitative levels).
Table 3. Association results of non-LS BAL cell types (log10-transformed of quantitative levels).
3.2. Löfgren’s syndrome (LS) group
Log10-transformed BAL macrophages associated with 261 LS-SNPs (denoted by 44 tag-SNPs after LD pruning). The lead SNP was rs9268861 (Beta = 0.09233, SE = 0.0243, p = 1.88 × 10−4, FDR = 2.5 × 10−2) located 17 kb of the 3′ of HLA-DRA (complete results see Supplementary Table S1). Eighty-five LS-SNPs (denoted by 27 tag-SNPs) were associated with log10-transformed BAL lymphocytes. The lead SNP was rs6457617 (Beta = 0.158, SE = 0.04782, p = 1.12E-03, FDR = 1.91E-02) located 29 kb 5′ of HLA-DQB1 (complete results, see Supplementary Table S2). Cell subgroups within lymphocytes were also identified. Specifically, 41 LS SNPs (denoted by 13 tag-SNPs) were associated with log10-transformed BAL CD4+. The lead SNP was rs1265061 (Beta = −0.01842, SE = 0.00776, p = 2.08E-02, FDR = 5.00E-02) located 741 bp 3′ of C6orf15 (complete results, see Supplementary Table S3). Sixteen LS SNPs (denoted by 10 tag-SNPs) were associated with log10-transformed BAL CD8+. The lead SNP was rs4947350 (Beta = −0.01842, SE = 0.00776, p = 2.08E-02, FDR = 5.00E-02) located 13 kb 3′ of HLA-DOB (complete results, see Supplementary Table S4). Twenty-seven LS SNPs (denoted by 9 tag-SNPs) were associated with log10-transformed BAL CD3+. The lead SNP was rs4947350 (Beta = 0.01792, SE = 0.007555, p = 2.09E-02, FDR = 4.18E-02) located 13 kb 3′ of HLA-DOB (complete results, Supplementary Table S5). Twenty LS SNPs (denoted by 5 tag-SNPs) were associated with log10-transformed BAL CD4+/CD8+. The lead SNP was rs2248462 (Beta = 0.1171, SE = 0.04368, p = 8.03E-03, FDR = 2.28E-02) located in AIF1 (complete results, see Supplementary Table S6). Ninety-three LS-SNPs (26 tag-SNPs) were associated with log10-transformed BAL neutrophils. The lead SNP was rs3130573 (Beta = −0.1218, SE = 0.04123, p = 3.48E-03, FDR = 4.91E-02) located in PSORS1C2 (complete results, see Supplementary Table S7). Seventy-seven LS SNPs (denoted by 22 tag-SNPs) were associated with log10-transformed BAL basophils. The lead SNP was rs2248462 (Beta = 0.6113, SE = 0.1958, p = 4.97E-03, FDR = 2.93E-02) located 6.6 kb 3′ of HCG26 (complete results, see Supplementary Table S8). Thirty-seven LS SNPs (denoted by 16 tag-SNPs) were associated with log10-transformed BAL eosinophils. The lead SNP rs3823417 (Beta = 0.2176, SE = 0.08061, p = 7.81E-03, FDR = 4.28E-02) located in PSORS1C1 (complete results, see Supplementary Table S9).
3.3. Non-Löfgren’s syndrome (non-LS) group
Sixty-nine non-LS SNPs (18 tag-SNPs) were associated with log10-transformed BAL macrophages. SNP-associated were located on chromosomes 6 and 10 (Table 3).
The lead SNP was rs2076536 (Beta = −0.06018, SE = 0.01919, p = 1.88 × 10−3, FDR = 2.03 × 10−2) located in C6orf10, (complete results, see Supplementary Table S10). Six non-LS SNPs (denoted by 3 tag-SNPs) were associated with log10-transformed BAL CD8+ cells. The lead SNP was rs4934167 (Beta = 0.05736, SE = 0.02209, p = 1.11E-02, FDR = 4.19E-02; complete results, see Supplementary Table S11). Thirty-one non-LS SNPs (7 tag-SNPs) were associated with log10-transformed BAL basophils. The lead SNP was rs2858332 (Beta = 0.2344, SE = 0.06674, p = 8.87E-04, FDR = 2.75E-02) located 28 kb 5′ of HLA-DQA2 (complete results, see Supplementary Table S12). 35 SNPs (7 tag-SNPs) were associated with log10-transformed BAL eosinophils. The lead SNP rs2857700 (Beta = 0.1616, SE = 0.06228, p = 1.01E-02, FDR = 4.66E-02) located 11 kb 5′ of AIF1 (complete results, see Supplementary Table S13). For complete results of associations, see Supplementary Tables S1-S13, available in figshare.4
As a secondary analysis, we evaluated the effect of the HLA-DBR1*03 allele on non-LS by excluding individuals with the HLA-DRB1*03 allele, given previous knowledge that carrying this allele can lead to a good prognosis, and re-ran the analysis. The findings from this assessment showed similar results for quantitative levels of macrophages, CD8, basophils, and eosinophils and associations with other cell types, including CD3, CD4, and neutrophils (for complete results, see Supplementary Table S14).
3.4. Cross-tissue and immune cell expression quantitative trait loci (eQTL) enrichment
Results from eQTL enrichment for LS and non-LS associated with quantitative levels of BALF cell types are summarized in Table 4. Briefly, 143 LS-SNPs associated with BAL cell types were eQTL SNPs (eSNP). Within this set, 46 tag-SNPs were found. eSNPs were correlated with the expression of various immune cells, peripheral blood, and lung using available resources. The most significant LS tag-SNPs for each BAL cell type coupled with eQTL data are shown in Table 5. Likewise, eQTL enrichment identified 35 non-LS SNPs associated with BALF cell types as eSNPs. Within the 35 eSNPs, eight were defined as tag-SNPs. The top non-LS tag SNPs coupled with eQTL data are shown in Table 6. Complete results of eQTL enrichment for LS and non-LS SNPs associated with quantitative levels of BAL cell types are shown in Supplementary Tables S15, S16, https://doi.org/10.6084/m9.figshare.20485038.
Table 4. Summary of eQTL mapping for LS and non-LS SNPs associated with BAL cell types.
Table 5. Top LS SNPs associated with BAL cell types identified as eQTL SNPs.
Table 6. Top non-LS tag-SNPs associated with BAL cell types identified as eQTL SNPs.
4. Discussion
In the present study, we assessed the relationship between LS- and non-LS-associated genetic variants (defined by SNPs) and quantitative levels of BALF cell populations in individuals with sarcoidosis. We also assessed significant findings with eQTL data across tissue and immune cell types relevant to sarcoidosis. Our results showed that sarcoidosis SNPs were significantly associated with the quantitative levels of BALF cells, suggesting a role of DNA variation on BAL cell phenotypes.
Genetic variants associated with LS (18) were significantly associated with quantitative levels of BALF macrophages, lymphocytes (CD3+, CD4+, CD8+, and CD4+/CD8+ ratio), neutrophils, basophils, and eosinophils, suggesting a wide range of cellular immuno-phenotypes implicated in the disease. LS-SNPs were primarily located on chromosome 6 from 25,844,710 to 33,032,975 base pairs, highlighting the role of the extended major histocompatibility complex (× MHC) region (30) in LS. The xMHC is a vast region with a complex linkage disequilibrium pattern. The ×MHC has been implicated in autoimmune infections and inflammatory diseases (31).
We also observed significant non-LS SNPs associated with the quantitative levels of BALF macrophages, CD8+ cells, basophils, and eosinophils. Non-LS SNPs were located on chromosome 6 from 31,449,327 to 33,073,322 base pairs, highlighting the role of the MHC class II subregion, a well-known and extensively studied locus in various autoimmune diseases. Besides this locus, we also observed non-LS SNPs on chromosome 10 from 82,142,176 to 82,206,348 associated with BALF cell quantitative levels. This locus is nearby the ANXA11 gene, a well-documented locus associated with sarcoidosis (32). ANXA11 has also been previously reported as a disease modifier and found to be associated with disease chronicity (33–35). A study by Hofmann et al., showed the downregulation of Annexin A11 upon activating CD8+ T cells and CD19+ B cells in sarcoidosis (36). A review by Kaneko et al., explained the orchestration of immune cells, such as cytotoxic CD4 and CD8 T cells and B cells and antibodies against Annexin A11 described as potential autoantigen could play key roles in the promotion of inflammation and fibrosis (37). The protein encoded by ANXA11 is a 56-kD antigen recognized by sera in patients with various autoimmune diseases, as reviewed in (38). Compiling these data, we offer further evidence that biological mechanisms underlying sarcoidosis may be autoimmune-mediated. Moreover, LS and non-LS SNPs are associated with other autoimmune diseases (39). Indeed, our findings showed how DNA variation associated with sarcoidosis risk modifies the expression of quantitative levels in BAL cells, which may have an implication on disease development, such as chronic phenotypes, as we have reported in a previous study (40).
The association between sarcoidosis (LS and non-LS) SNPs and quantitative levels of BAL macrophages is noteworthy. Macrophages are regarded as the main source of TNF-alpha, a primary cytokine for granuloma formation. The association between sarcoidosis SNPs and BAL CD8+ cells, particularly in non-LS patients, is an attractive finding, as higher levels of cytotoxic lymphocytes have been reported to be associated with chronic phenotypes (41).
Regarding basophils and eosinophils, there is growing evidence that these cell types may have a functional role in developing autoimmune and inflammatory diseases (42, 43). Eosinophils have strong cytotoxic properties and are modulators of other immune cells, which could apply in sarcoidosis. On the other hand, basophils have been reported to play a role in the pathogenesis of lupus nephritis (44) via the activation of autoreactive IgE, thymic stromal lymphopoietin, or Toll-like receptors (TLRs). TLR-mediated mechanisms and SNPs associated with different TLRs have been reported to be associated with disease progression in sarcoidosis (45). Interestingly, basophils have been suggested to be involved in B cell survival, activation, and differentiation (46). As new reports emerge, B cells may be key players in the immunopathogenesis of sarcoidosis (47–49). Thus, further investigations are needed to understand the role of basophils and eosinophils in sarcoidosis.
Cumulative evidence shows that eQTLs are strongly enriched for heritability across complex traits (50). eQTL mapping studies have shown that disease-associated loci revealed by genome-wide association studies are enriched in eQTLs and regulatory genes. However, eQTLs explain only a small proportion of the expression heritability as many of eQTLs may be in linkage-disequilibrium (LD) with disease-associated loci, as described in (51). Umans et al. depict the pearls and pitfalls of eQTLs in disease-associated loci and offer insights for future studies.
In this study, we took advantage of comprehensive eQTL sources, which identified several sarcoidosis-associated SNPs associated with quantitative levels of BAL cell types. Chiefly, sarcoidosis-associated SNPs were highly correlated with gene expression of relevant tissues and cell types, including immune cells, blood, and lung. Thus, suggesting a plausible functional mechanism between sarcoidosis gene-variants and regulation of gene expression in lung, blood, and immune cell types relevant to the disease.
The mapping of bulk and single-cell eQTLs on quantitative levels of BAL cells advances our knowledge previously reported in (18), where we showed that sarcoidosis-associated loci were enriched with cis-eQTLs and regulatory SNPs. Specifically, from bulk RNA of blood, we identified over 3,000 cis-eQTLs in LS and over 500 in non-LS, suggesting that regulation of gene expression is likely a dominant mediator for disease risk. Moreover, from bulk RNA in immune cells, we identified several cis-eQTLs significantly correlated with gene expression of B cells, monocytes, macrophages, neutrophils, natural killer (NK) cells, CD8 T cells, and CD4 T cells (i.e., naïve and memory regulatory T cells (Tregs), CD4 Th1 and Th1.17 cells, and CD4 Th2 cells). We also identified cis-eQTLs in the lung and blood derivates (platelets and PBMCs). Notably, although the specificity of cell type in bulk RNA limits the accuracy of these findings, we believe that besides CD4 Th1 cells, other immune cell types that have been previously neglected in the past, are involved in promoting granuloma formation, as highlighted in recent reviews (52, 53). Using single-cell RNA data, we identified 3 single-cell eQTLs in PBMCs and CD4 T cells in LS. sc-eQTLs allow to capture of dynamic processes in different cell-specific types and estimate the variability in gene expression in each cell (54).
Our work refines the immunobiology of sarcoidosis and suggests new incentives for further investigations using single-cell transcriptomics studies and studies for discovering response eQTLs, which are suitable for investigating immune responses under disease treatment conditions (51). Indeed, our study offers continuing knowledge by providing new evidence on the specificity of gene expression in immune cells and their cell type-specific features. It is also worth noting that most sarcoidosis variants lie in non-coding genomic regions, so disease mechanisms are likely driven by common regulatory variants that modify gene expression of target genes in cells and tissues relevant to the disease (55). Nonetheless, our findings should be interpreted with caution in a comparative study as additional evidence is needed to characterize causal variants and determine gene causality.
4.1. Limitations of the study
Although our study provides comprehensive information on the relationship between sarcoidosis susceptibility and quantitative levels of BAL cell types, there are some limitations. First, our study is limited by the small number of genetic variants associated with sarcoidosis (LS and non-LS). It is essential to highlight that these variants are derived from the ImmunoChip SNP-array, a customized array focused on genetic variants associated with immune-mediated disorders. The mapping of this SNP array covers a small fraction of genetic variants in the genome. Second, functional interpretation of sarcoidosis-associated genetic variants on the expression of BAL cells is missing as we do not have BAL transcriptomic data in our cohort. eQTL annotation for BALF is also not present in available sources.
5. Conclusion
Our finding indicates that sarcoidosis genetic variants associated with LS and non-LS are significantly associated with quantitative levels of BAL cell types. The identification of sarcoidosis-associated SNPs as regulators of gene expression via lung, immune cells, and blood eQTLs offers new insights into the immunobiology of sarcoidosis.
Data availability statement
Summary statistics for this study can be found in the GWAS catalog https://www.ebi.ac.uk/gwas/ under accession numbers GCST005540 and GCST005543. Further association results on SNPs and eQTL enrichment for LS and non-LS SNPs associated with BAL cell types quantitative levels can be found in Supplementary Tables S1–S13, available in Figshare https://doi.org/10.6084/m9.figshare.20485038.
Ethics statement
The studies involving human participants were reviewed and approved by Ethical Review Board in Stockholm. The patients/participants provided their written informed consent to participate in this study.
Author contributions
NR designed and executed data analysis. MA was involved in data analysis and drafted the first version of the manuscript. NR, LP, JG, SK, MA, and AE revised the manuscript and provided feedback. SK, JG, and AE phenotyped patients included in the study. All authors contributed to the article and approved the submitted version.
Funding
This project is funded by the Swedish Heart-Lung Foundation, awarded to NR (grant no. 20170664, 20200506), SK (grant no. 20200163), and JG (grant no. 20190478); Karolinska Institutet Foundation awarded to NR (grant no. FS-2018:0007); Swedish Research Council awarded to LP (grant no. 2018-02884) and JG (grant no. 2019-01034). MA was supported by grants from the Swedish Heart-Lung Foundation, The Swedish Research Council, King Gustaf V’s, and Queen Victoria’s Freemasons’ Foundation through the regional agreement on medical training and clinical research (ALF) between Stockholm County Council and Karolinska Hospital, and Karolinska Institutet. The computations and genomic data handling were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX), partially funded by the Swedish Research Council through grant agreement no. 2018-05973.
Acknowledgments
The authors thank all patients and the personnel involved in this study. We thank our research nurses, Margitha Dahl, Helene Blomqvist, and Susanne Schedin.
Conflict of interest
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.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmed.2023.1061654/full#supplementary-material
Footnotes
1. ^http://prioritypruner.sourceforge.net/
2. ^https://www.ebi.ac.uk/eqtl/Studies/
References
1. Grunewald, J, Grutters, JC, Arkema, EV, Saketkoo, LA, Moller, DR, and Muller-Quernheim, J. Sarcoidosis. Nat Rev Dis Primers. (2019) 5, 45. doi: 10.1038/s41572-019-0096-x
2. Costabel, U, and Hunninghake, GW. ATS/ERS/WASOG statement on sarcoidosis. Sarcoidosis statement committee. American Thoracic Society. European Respiratory Society. World Association for Sarcoidosis and Other Granulomatous Disorders. Eur Respir J. (1999) 14:735–7. doi: 10.1034/j.1399-3003.1999.14d02.x
3. Valeyre, D, Prasse, A, Nunes, H, Uzunhan, Y, Brillet, PY, and Muller-Quernheim, J. Sarcoidosis. Lancet. (2014) 383:1155–67. doi: 10.1016/S0140-6736(13)60680-7
4. Grunewald, J, and Eklund, A. Lofgren's syndrome: human leukocyte antigen strongly influences the disease course. Am J Respir Crit Care Med. (2009) 179:307–12. doi: 10.1164/rccm.200807-1082OC
5. Olsen, HH, Grunewald, J, Tornling, G, Skold, CM, and Eklund, A. Bronchoalveolar lavage results are independent of season, age, gender and collection site. PLoS One. (2012) 7:e43644. doi: 10.1371/journal.pone.0043644
6. Drent, M, Mansour, K, and Linssen, C. Bronchoalveolar lavage in sarcoidosis. Semin Respir Crit Care Med. (2007) 28:486–95. doi: 10.1055/s-2007-991521
7. Kaiser, Y, Eklund, A, and Grunewald, J. Moving target: shifting the focus to pulmonary sarcoidosis as an autoimmune spectrum disorder. Eur Respir J. (2019) 54:1802153. doi: 10.1183/13993003.021532018
8. Haggmark, A, Hamsten, C, Wiklundh, E, Lindskog, C, Mattsson, C, Andersson, E, et al. Proteomic profiling reveals autoimmune targets in sarcoidosis. Am J Respir Crit Care Med. (2015) 191:574–83. doi: 10.1164/rccm.201407-1341OC
9. Korsten, P, Tampe, B, Konig, MF, and Nikiphorou, E. Sarcoidosis and autoimmune diseases: differences, similarities and overlaps. Curr Opin Pulm Med. (2018) 24:504–12. doi: 10.1097/MCP.0000000000000500
10. Grunewald, J, Spagnolo, P, Wahlström, J, and Eklund, A. Immunogenetics of disease-causing inflammation in sarcoidosis. Clin Rev Allergy Immunol. (2015) 49:19–35. doi: 10.1007/s12016-015-8477-8
11. Moller, DR, Rybicki, BA, Hamzeh, NY, Montgomery, CG, Chen, ES, Drake, W, et al. Genetic, immunologic, and environmental basis of sarcoidosis. Ann Am Thorac Soc. (2017) 14:S429–36. doi: 10.1513/AnnalsATS.201707-565OT
12. Sakthivel, P, and Bruder, D. Mechanism of granuloma formation in sarcoidosis. Curr Opin Hematol. (2017) 24:59–65. doi: 10.1097/MOH.0000000000000301
13. Linke, M, Pham, HTT, Katholnig, K, Schnoller, T, Miller, A, Demel, F, et al. Chronic signaling via the metabolic checkpoint kinase mTORC1 induces macrophage granuloma formation and marks sarcoidosis progression. Nat Immunol. (2017) 18:293–302. doi: 10.1038/ni.3655
14. Grunewald, J, and Eklund, A. Role of CD4+ T cells in sarcoidosis. Proc Am Thorac Soc. (2007) 4:461–4. doi: 10.1513/pats.200606-130MS
15. Kaiser, Y, Lepzien, R, Kullberg, S, Eklund, A, Smed-Sorensen, A, and Grunewald, J. Expanded lung T-bet(+)ROR gamma T+ CD4(+) T-cells in sarcoidosis patients with a favourable disease phenotype. Eur Respir J. (2016) 48:484–94. doi: 10.1183/13993003.00092-2016
16. Lareau, CA, DeWeese, CF, Adrianto, I, Lessard, CJ, Gaffney, PM, Lannuzzi, MC, et al. Polygenic risk assessment reveals pleiotropy between sarcoidosis and inflammatory disorders in the context of genetic ancestry. Genes Immun. (2017) 18:88–94. doi: 10.1038/gene.2017.3
17. Calender, A, Weichhart, T, Valeyre, D, and Pacheco, Y. Current insights in genetics of sarcoidosis: functional and clinical impacts. J Clin Med. (2020) 9. doi: 10.3390/jcm9082633
18. Rivera, NV, Ronninger, M, Shchetynsky, K, Franke, A, Nothen, MM, Muller-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:1008–22. doi: 10.1164/rccm.201507-1372OC
19. Zhou, T, Zhang, W, Sweiss, NJ, Chen, ES, Moller, DR, Knox, KS, et al. Peripheral blood gene expression as a novel genomic biomarker in complicated sarcoidosis. PLoS One. (2012) 7:e44818. doi: 10.1371/journal.pone.0044818
20. Gonzalez-Garay, ML, Knox, KS, and Garcia, JGN. Identification of complicated sarcoidosis genes via GWAS SNP and eQTL analyses, C36. Clin. Stud. Sarcoidosis :A4821–1 MeetingAbstracts.A4821 https://www.atsjournals.org/doi/abs/10.1164/ajrccm-conference.2018.197.1_MeetingAbstracts.A4821.
21. Gamazon, ER, Segrè, AV, Van De Bunt, M, Wen, X, Xi, HS, Hormozdiari, F, et al. Using an atlas of gene regulation across 44 human tissues to inform complex disease- and trait-associated variation. Nat Genet. (2018) 50:956–67. doi: 10.1038/s41588-018-0154-4
22. Statement on Sarcoidosis. Joint statement of the American Thoracic Society (ATS), the European Respiratory Society (ERS) and the world Association of Sarcoidosis and Other Granulomatous Disorders (WASOG) adopted by the ATS Board of directors and by the ERS executive committee. Am J Respir Crit Care Med. (1999) 160:736–55.
23. Royston, JP. An extension of Shapiro and Wilk-W test for normality to large samples. J R Stat Soc C-Appl. (1982) 31:115–24.
24. Chang, CC, Chow, CC, Tellier, LC, Vattikuti, S, Purcell, SM, and Lee, JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience. (2015) 4:7. doi: 10.1186/s13742-015-0047-8
25. Watanabe, K, Taskesen, E, Van Bochoven, A, and Posthuma, D. Functional mapping and annotation of genetic associations with FUMA. Nat Commun. (2017) 8:1826. doi: 10.1038/s41467-017-01261-5
26. Zhernakova, DV, Deelen, P, Vermaat, M, Van Iterson, M, Van Galen, M, Arindrarto, W, et al. Identification of context-dependent expression quantitative trait loci in whole blood. Nat Genet. (2017) 49:139–45. doi: 10.1038/ng.3737
27. Schmiedel, BJ, Singh, D, Madrigal, A, Valdovino-Gonzalez, AG, White, BM, Zapardiel-Gonzalo, J, et al. Impact of genetic polymorphisms on human immune cell gene expression. Cells. (2018) 175:1701–1715.e16. doi: 10.1016/j.cell.2018.10.022
28. Võsa, U, Claringbould, A, Westra, H-J, Bonder, MJ, Deelen, P, Zeng, B, et al. Large-scale cis- and trans-eQTL analyses identify thousands of genetic loci and polygenic scores that regulate blood gene expression. Nat Genet. (2021) 53:1300–10. doi: 10.1038/s41588-018-0154-4
29. Van Der Wijst, MGP, Brugge, H, De Vries, DH, Deelen, P, Swertz, MA, and Franke, L. Single-cell RNA sequencing identifies celltype-specific cis-eQTLs and co-expression QTLs. Nat Genet. (2018) 50:493–7. doi: 10.1038/s41588-018-0089-9
30. Horton, R, Wilming, L, Rand, V, Lovering, RC, Bruford, EA, Khodiyar, VK, et al. Gene map of the extended human MHC. Nat Rev Genet. (2004) 5:889–99. doi: 10.1038/nrg1489
31. Matzaraki, V, Kumar, V, Wijmenga, C, and Zhernakova, A. The MHC locus and genetic susceptibility to autoimmune and infectious diseases. Genome Biol. (2017) 18:76. doi: 10.1186/s13059-017-1207-1
32. Mirsaeidi, M, Gidfar, S, Vu, A, and Schraufnagel, D. Annexins family: insights into their functions and potential role in pathogenesis of sarcoidosis. J Transl Med. (2016) 14:89. doi: 10.1186/s12967-016-0843-7
33. Li, Y, Pabst, S, Kubisch, C, Grohe, C, and Wollnik, B. First independent replication study confirms the strong genetic association of ANXA11 with sarcoidosis. Thorax. (2010) 65:939–40. doi: 10.1136/thx.2010.138743
34. Mrazek, F, Stahelova, A, Kriegova, E, Fillerova, R, Zurkova, M, Kolek, V, et al. Functional variant ANXA11 R230C: true marker of protection and candidate disease modifier in sarcoidosis. Genes Immun. (2011) 12:490–4. doi: 10.1038/gene.2011.27
35. 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
36. 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
37. Kaneko, N, Moriyama, M, Maehara, T, Chen, H, Miyahara, Y, and Nakamura, S. Orchestration of immune cells contributes to fibrosis in IgG4-related disease. Immuno. (2022) 2:170–84. doi: 10.3390/immuno2010013
38. Wang, J, Guo, C, Liu, S, Qi, H, Yin, Y, Liang, R, et al. Annexin A11 in disease. Clin Chim Acta. (2014) 431:164–8. doi: 10.1016/j.cca.2014.01.031
39. Polychronakos, C. Fine points in mapping autoimmunity. Nat Genet. (2011) 43:1173–4. doi: 10.1038/ng.1015
40. Al Hayja, MA, Wahlström, J, Kullberg, S, Darlington, P, Eklund, A, and Grunewald, J. Bronchoalveolar lavage fluid cell subsets associate with the disease course in Löfgren's and non-Löfgren's sarcoidosis patients. Respir Med. (2021) 186:106521. doi: 10.1038/s41588-018-0154-4
41. Parasa, VR, Forsslund, H, Enger, T, Lorenz, D, Kullberg, S, Eklund, A, et al. Enhanced CD8(+) cytolytic T cell responses in the peripheral circulation of patients with sarcoidosis and non-Lofgren's disease. Respir Med. (2018) 138:S38–44. doi: 10.1016/j.rmed.2017.10.006
42. Diny, NL, Rose, NR, and Čiháková, D. Eosinophils in autoimmune diseases. Front Immunol. (2017) 8:484. doi: 10.3389/fimmu.2017.00484
43. Sharma, M, and Bayry, J. Basophils in autoimmune and inflammatory diseases. Nat Rev Rheumatol. (2015) 11:129–31. doi: 10.1038/nrrheum.2014.199
44. Charles, N, Hardwick, D, Daugas, E, Illei, GG, and Rivera, J. Basophils and the T helper 2 environment can promote the development of lupus nephritis. Nat Med. (2010) 16:701–7. doi: 10.1038/nm.2159
45. Cooke, G, Kamal, I, Strengert, M, Hams, E, Mawhinney, L, Tynan, A, et al. Toll-like receptor 3 L412F polymorphism promotes a persistent clinical phenotype in pulmonary sarcoidosis. QJM Int J Med. (2018) 111:217–24. doi: 10.1093/qjmed/hcx243
46. Merluzzi, S, Betto, E, Ceccaroni, AA, Magris, R, Giunta, M, and Mion, F. Mast cells, basophils and B cell connection network. Mol Immunol. (2015) 63:94–103. doi: 10.1016/j.molimm.2014.02.016
47. Kudryavtsev, I, Serebriakova, M, Starshinova, A, Zinchenko, Y, Basantsova, N, Malkova, A, et al. Imbalance in B cell and T follicular helper cell subsets in pulmonary sarcoidosis. Sci Rep. (2020) 10:1059. doi: 10.1038/s41598-020-57741-0
48. Bauer, L, Müller, LJ, Volkers, SM, Heinrich, F, Mashreghi, M-F, Ruppert, C, et al. Follicular helper–like T cells in the lung highlight a novel role of B cells in sarcoidosis. Am J Respir Crit Care Med. (2021) 204:1403–17. doi: 10.1164/rccm.202012-4423OC
49. Ando, M, Goto, A, Takeno, Y, Yamasue, M, Komiya, K, Umeki, K, et al. Significant elevation of the levels of B-cell activating factor (BAFF) in patients with sarcoidosis. Clin Rheumatol. (2018) 37:2833–8. doi: 10.1007/s10067-018-4183-2
50. Hormozdiari, F, Gazal, S, Van De Geijn, B, Finucane, HK, Ju, CJT, Loh, P-R, et al. Leveraging molecular quantitative trait loci to understand the genetic architecture of diseases and complex traits. Nat Genet. (2018) 50:1041–7. doi: 10.1038/s41588-018-0148-2
51. Umans, BD, Battle, A, and Gilad, Y. Where are the disease-associated eQTLs? Trends Genet. (2021) 37:109–24. doi: 10.1016/j.tig.2020.08.009
52. Locke, LW, Schlesinger, LS, and Crouser, ED. Current sarcoidosis models and the importance of focusing on the granuloma. Front Immunol. (2020) 11: 1719. doi: 10.3389/fimmu.2020.01719
53. Zhang, H, Costabel, U, and Dai, H. The role of diverse immune cells in sarcoidosis. Front Immunol. (2021) 12:788502. doi: 10.3389/fimmu.2021.788502
54. Eling, N, Richard, AC, Richardson, S, Marioni, JC, and Vallejos, CA. Correcting the mean-variance dependency for differential variability testing using single-cell RNA sequencing data. Cell Syst. (2018) 7:284–294.e12. doi: 10.1016/j.cels.2018.06.011
Keywords: sarcoidosis, single nucleotide polymorphisms, immunogenetics, eQTLs, gene regulation, HLA, major histocompatibility complex, broncoalveolar lavage
Citation: Abo Al Hayja M, Kullberg S, Eklund A, Padyukov L, Grunewald J and Rivera NV (2023) Functional link between sarcoidosis-associated gene variants and quantitative levels of bronchoalveolar lavage fluid cell types. Front. Med. 10:1061654. doi: 10.3389/fmed.2023.1061654
Edited by:
Paolo Cameli, University of Siena, ItalyCopyright © 2023 Abo Al Hayja, Kullberg, Eklund, Padyukov, Grunewald and Rivera. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Natalia V. Rivera, ✉ natalia.rivera@ki.se
†These authors share first authorship