Identification of NF-κB and PLCL2 as new susceptibility genes and highlights on a potential role of IRF8 through interferon signature modulation in systemic sclerosis

Introduction Systemic sclerosis (SSc) and primary biliary cirrhosis (PBC) are rare polygenic autoimmune diseases (AIDs) characterized by fibroblast dysfunction. Furthermore, both diseases share some genetic bases with other AIDs, as evidenced by autoimmune gene pleiotropism. The present study was undertaken to investigate whether single-nucleotide polymorphisms (SNPs) identified by a large genome-wide association study (GWAS) in PBC might contribute to SSc susceptibility. Methods Sixteen PBC susceptibility SNPs were genotyped in a total of 1,616 patients with SSc and 3,621 healthy controls from two European populations (France and Italy). Results We observed an association between PLCL2 rs1372072 (odds ratio (OR) = 1.22, 95% confidence interval (CI) 1.12 to 1.33, Padj = 7.22 × 10−5), nuclear factor-kappa-B (NF-κB) rs7665090 (OR = 1.15, 95% CI 1.06 to 1.25, Padj = 0.01), and IRF8 rs11117432 (OR = 0.75, 95% CI 0.67 to 0.86, Padj = 2.49 × 10−4) with SSc susceptibility. Furthermore, phenotype stratification showed an association between rs1372072 and rs11117432 with the limited cutaneous subgroup (lcSSc) (Padj = 4.45 × 10−4 and Padj = 0.001), whereas rs7665090 was associated with the diffuse cutaneous subtype (dcSSc) (Padj = 0.003). Genotype-mRNA expression correlation analysis revealed that the IRF8 protective allele was associated with increased interferon-gamma (IFN-γ) expression (P = 0.03) in patients with SSc but decreased type I IFN (IFIT1) expression in patients and controls (P = 0.02). In addition, we found an epistatic interaction between NF-κB and IRF8 (OR = 0.56, 95% CI 0.00 to 0.74, P = 4 × 10−4) which in turn revealed that the IRF8 protective effect is dependent on the presence of the NF-κB susceptibility allele. Conclusions An analysis of pleiotropic genes identified two new susceptibility genes for SSc (NF-κB and PLCL2) and confirmed the IRF8 locus. Furthermore, the IRF8 variant influenced the IFN signature, and we found an interaction between IRF8 and NF-κB gene variants that might play a role in SSc susceptibility. Electronic supplementary material The online version of this article (doi:10.1186/s13075-015-0572-y) contains supplementary material, which is available to authorized users.


Introduction
Systemic sclerosis (SSc) is an autoimmune disease (AID) characterized by vascular damage, autoantibody production, and fibrotic events. SSc is an orphan disease that is considered the most severe connective tissue disorder, representing an important medical challenge because of its debilitating progressive nature [1]. The precise aetiology of the disease remains unclear. However, as in other AIDs, interactions between environment and genetic factors are thought to play a key role in disease susceptibility [2]. Indeed, several reports of well-powered candidate gene studies, together with genome-wide association studies (GWASs), have established numerous SSc susceptibility genes, including STAT4, IRF5, TNFSF4, CD226, CD247, and BANK1 [3,4]. The vast majority of these susceptibility loci belong to pathways involved in immune responses or inflammation and were identified in other AIDs [3]. Previous work from our team demonstrated that polyautoimmunity affects up to a quarter of patients with SSc, highlighting the concept of shared autoimmunity in SSc [5]. It is well established that AIDs encompass a broad range of phenotypic manifestations and severity, indicating that the effects of these loci are not of equal magnitude and may not be associated in the same direction (risk or protection) among different diseases [3]. Fibroblast dysfunction, the SSc hallmark, can also be observed in primary biliary cirrhosis (PBC). PBC is a disease marked by progressive destruction of the liver and can co-occur with SSc, suggesting that these two AIDs may share common pathways and may be a result of common-variant loci of weak effect [6]. A recent GWAS carried out in patients with PBC identified 12 new loci involved in disease susceptibility and highlighted the importance of type I interferon (IFN), nuclear factor-kappa-B (NF-κB), and Toll-like receptor (TLR) signaling [7] in its pathogenesis. Given the evidence of shared autoimmune genes between SSc and other AIDs, we investigated the association of 16 single-nucleotide polymorphisms (SNPs), recently identified as susceptibility factors for PBC, with SSc and its subphenotypes.

Study population
We performed a case-control association study with a European Caucasian cohort consisting of 1,616 patients with SSc (1,022 French and 594 Italians) and 3,621 controls (2,384 French and 1,237 Italians). Only individuals with European ancestry (defined as all four grandparents being of European Caucasian ancestry) were included in the study. This cohort has already been used in several previous genetic studies, and high homogeneity between French and Italian samples has been demonstrated [4,8]. Control subjects were age-and sex-matched to patients with SSc. The characteristics of the patients with SSc are shown in Table 1. All patients with SSc were classified by LeRoy's cutaneous subtype and were phenotypically screened for anti-nuclear antibodies and pulmonary fibrosis (PF). The latter was defined as any of the following changes on computed tomography scan, such as in our previous studies: reticulations, honeycombing, and traction. In line with all of our previous projects, patients with SSc were also evaluated for the presence of other potential AIDs. To avoid any bias due to a possible excess of the risk alleles attributable to these patients, 23 SSc patients showing PBC co-occurrence were excluded from the main study. This study was approved by the Comité Consultatif de Protection des Personnes dans la Recherche Biomédicale (CCPPRB) Paris Cochin (dossier 2068). All participants enrolled in this study provided written informed consent.

Single-nucleotide polymorphism selection
The SNP selection was based on a previous PBC GWAS, which identified these SNPs as susceptibility genes to PBC [7]. Among these candidates, two SNPs, STAT4 rs10931468 and IRF5 rs12531711, located in loci reproducibly strongly associated with SSc in several studies [3,4], were excluded from the study. The list of the selected SNPs, as well as information about their localization and putative function, is shown in Additional file 1: Table S1.

Gene expression in peripheral blood mononuclear cells
To assess a possible genotype-phenotype correlation, we selected the three associated genes-PLCL2, NFKB1, and IRF8-to perform gene expression assays of their mRNA. Moreover, because of compelling data on IRF8-associated pathways, we also investigated IFN-γ and type I IFN (IFIT1) mRNA gene expression. We randomly selected peripheral blood nuclear cell (PBMC) samples from 39 patients (who did not receive immunosuppressive drugs) and 24 healthy controls. The total RNA was obtained by using an RNA extraction kit (RNeasy Mini Kit; Qiagen, Hilden, Germany) followed by a c-DNA reverse transcription step (SuperScript II Reverse Transcriptase; Invitrogen To test SNP × SNP epistasis, the allele frequencies were compared in patients and control subjects. The PLCL2*T alleles, NF-kB*G alleles, and IRF8*A alleles were considered binary random variables, and the dependency of disease on pairwise combinations of these covariates was tested with a logistic regression model. The fitted model included each one of two tested factors as first-order components and an interaction term reflecting the extent of the allelic by allelic epistasis. We applied a logistic regression assuming a double-dominant model. Exact calculations of ORs and P values were made with Logxact 8 software (Cytel).
We applied the Bonferroni correction although it is very conservative; 16 tag SNPs were considered to test the association between each SNP and SSc (P values were therefore multiplied by 16 to get the adjusted ones). For disease phenotypes correction, we took into account seven subtypes (P values were multiplied by 7 to get the adjusted ones). An unpaired two-sided t test was used for mRNA expression levels analysis. An adjusted P value of less than 0.05 was considered statistically significant.

Single-marker analysis
Sixteen SNPs were genotyped in two European populations (French and Italian). Among the SNPs of interest, minor allele frequencies (MAFs) range from 11% to 44%; therefore, we provide a power calculation for these two MAFs. For an 11% MAF, our power to detect an association is more than 99% for an OR of 1.5 and is 81% to detect an OR of 1.3. For a 44% MAF, our power is more than 99% to detect an OR of 1.3 or higher. Genotype frequencies were in Hardy-Weinberg equilibrium in the control population for all of the SNPs investigated.
When the 23 SSc-PBC patients were analysed separately, no significant association was observed when they were compared with the 3,621 controls for any SNPs; furthermore, their addition to the whole sample (23 SSc-PBC + 1,628 SSc patients) did not modify the results (significant association for PLCL2, IRF8, and NF-κB, with unchanged magnitude of the effects). In regard to the group of SSc patients having another AID (n = 224), single-marker analyses revealed significant association with PLCL2 rs1372072 only (P = 0.001; OR = 1.375, 95% CI 1.13 to 1.67).

Follow-up of the raised hypotheses for single markers
SNPs belonging to the genes highlighted herein have been previously reported as associated or with a trend for association in previous works. For some genes highlighted in our study, some variants not primarily herein studied have been reported as associated with SSc or as having a trend for association with SSc. Indeed, Martin et al. [10] reported a suggestive association of NF-κB rs1598859, and a robust association has been reported for IRF8 rs11642873 [11]. Furthermore, a recent work did suggest some association in two loci for which we herein identified a nominal association but not reaching significance after multi-test correction as we observe for DDX6 rs7130875 and IL12A rs77583790 [12].
In our sample, we therefore investigated these markers and compared the association with the one observed with our newly identified SNPs. To that end, we used a subset of our patients for whom we have performed genome-wide typing in the past (sample of 564 patients with SSc and 488 healthy controls) [13]. In regard to NF-κB, we first observed in HapMap that rs1598859 and rs7665090 are not in complete linkage disequilibrium (LD) (D′ = 0.77 and R 2 = 0.31; SNPs are on the same block according to D′ value but are not in phase according to R 2 ) and then found from our GWAS data that the MAFs were 0.355 in patients with SSc and 0.382 in controls (uncorrected P value = 0.2). For IRF8, an established locus, we first observed that rs11642873 and rs11117432 are not in complete LD (D′ = 0.707 and R 2 = 0.456). Within our subset with GWAS data, we identified that the rs11642873 MAFs were 0.14 in the patients with SSc and 0.174 in the controls (P value = 0.034 without correction; OR = 0.772, 95% CI 0.61 to 0.98). The P value is less significant than the one obtained for rs11117432 in the same subset (MAF = 0.156 versus 0.213; P = 0.0008), and the OR also shows a weaker effect (OR = 0.68, 95% CI 0.55 to 0.85). We also performed haplotype analyses to estimate the effects of the two SNPs and that relates to conditional analyses; we observed that the haplotype combination including the single rs11642873 risk variant (CG haplotype) was not associated with SSc risk (frequency of 0.108 versus 0.117; P = 0.499) but that the haplotype combination including the single rs11117432 risk variant (AA haplotype) was significantly associated with SSc risk (0.124 versus 0.156; P = 0.032). We proceeded similarly for DDX6 rs7130875 and observed that it was not in LD with rs6421571 (D′ = 0.839 and R 2 = 0.084). DDX6 rs7130875 was not included in our array, but we used the perfect proxy rs4499035 and identified that the MAFs were 0.261 in patients with SSc and 0.252 in the controls (uncorrected P value = 0.69). IL12A rs77583790 SNP was not included in the array used in our GWAS, and no accurate imputation could be done. Therefore, we have performed a dedicated genotyping by using the TaqMan method in a subset of the patients and controls for whom DNA was still available and who did participate in the present study; we found that rs77583790 minor A allele was present in 2.33% of 1,585 patients with SSc and in 1.36% of 1,842 healthy controls showing a nominal association but not supporting any correction for multiple testing (uncorrected P = 0.04). We have thereafter performed some meta-analysis by using the data included in the cited works [10][11][12] together with ours, thanks to Pr Javier Martin, who provided the genotypic data from his original studies (the sample size was not always the same as in the original study which could include several cohorts from various origins). As shown in Additional file 3: Table S3, neither the P value nor the OR was improved for NF-κB, DDX6, and IL12A genes. This suggests that adding our data does not strengthen the previous findings and does not support a  significant association for NF-κB rs1598859 and DDX6 rs7130875. In regard to IL12A rs77583790, our results do not provide an independent replication; however, the meta-analysis shows a more significant P value when the French population is added, suggesting that this locus is further confirmed by our data. For IRF8 rs11642873, our data decrease the P value, and although the OR is unchanged, the CI is reduced, suggesting that the association is both stronger and of higher magnitude (Additional file 3: Table S3).  Figure 1 shows the ORs for SSc patients with one or two and three to six risk alleles. Our results demonstrate that the risk of SSc increases proportionally as the number of alleles increases, with an OR of 1.70 when six risk alleles are considered (Figure 1).

Genetic interaction between associated single-nucleotide polymorphisms
To investigate possible gene-gene interactions between the associated SNPs, we applied a logistic regression (multivariate analysis) on PLCL2 rs1372072*T, NF-κB rs7665090*G, and IRF8 rs1117432*A alleles. Of the most interest, we found an epistatic interaction between IRF8 and NF-κB (OR = 0.56, 95% CI 0.00 to 0.74, P = 4 × 10 −4 ), which in turn demonstrated that the protective effect of the IRF8 A allele is revealed when the NF-κB G allele is present (Figure 2A and B) Figure 2C).

mRNA expression
We then investigated the potential influence of the three associated gene variants PLCL2, NF-κB, and IRF8 on their mRNA expression levels in 38 patients with SSc and 24 controls. After stratification according to the various genotypes, we observed no genotype-phenotype correlation. However, we further analyzed a possible effect of the variant IRF8 rs11117432 on IFIT1 and IFN-γ mRNA relative expression in the same 38 patients with SSc (nine of them with the IRF8 rs11117432 AA or GA genotypes and 29 with the GG genotype) and 24 controls (12 with the IRF8 rs11117432 AA or GA genotypes and 12 with the GG genotype) according to the published data on IRF8 effects on IFN signature. Using a multiple linear regression test, we found out that the IRF8 rs11117432*A allele was associated with IFIT1 mRNA expression levels independently of the status of the included individuals (P = 0.022, data not shown). When all individuals were considered, our results demonstrated that the IRF8 rs11117432*A variant was associated with decreased levels of IFIT1 mRNA expression (P = 0.02). When patients and controls were considered separately, a trend was observed in the comparison of the IRF8 rs11117432*A carriers with the non-carriers (P = 0.10 for patients, P = 0.06) ( Figure 3A). We also observed increased levels of IFN-γ mRNA expression in SSc patient carriers of IRF8 rs11117432*A when compared with the non-carriers (P = 0.03) with a trend for difference for the whole group (P = 0.11) but no difference in controls (P = 0.35) ( Figure 3B). There was no effect of the epistatic interaction on IFIT1 or IFN-γ expression (P = 0.59 for IFIT1 and P = 0.50 for IFN-γ).

Discussion
The present study identified NF-κB and PLCL2 as new SSc susceptibility genes and confirmed the IRF8 locus.
We also showed potential effects on the IFN signature in regard to the IRF8 gene variant and an interaction between IRF8 and NF-κB genes.
Our results corroborate previously reported associations in regard to the IRF8 gene and SSc susceptibility. Here, we demonstrate, for the first time, the specific association between IRF8 rs11117432 SNP and SSc susceptibility and provide data suggesting a stronger association for rs11117432 as compared with rs11642873, although our meta-analysis results strengthen the previous findings for this marker [11]. Indeed, our haplotype analyses further confirm a stronger association for rs11117432 for a subset of our patients, but a denser genotyping would be required to better narrow a potential causal variant. IRF8 is a nuclear protein, which upon stimulation by pathogenassociated molecular pattern molecules (PAMPs) moves into the cytoplasm and activates NF-κB and TLR signaling pathways, leading to cytokine production and thereby regulating inflammatory responses [14]. The consequences of unregulated inflammation are associated with the development of pathologic fibrosis such as that seen in SSc [15]. Here, we show that the variant IRF8 rs11117432 increases IFN-γ but decreases IFIT1 gene expression. Despite reaching significance, our results show subtle differences, and given the complex regulation of protein synthesis and activity, more work is required to confirm the functional impacts of the IRF8 variant. IFNs are a group of cytokines that play an important role in the regulation of inflammation [16]. Type I IFNs, such as IFIT1, seem to contribute to the upregulation of TLR-induced inflammation and transforming growth factor-beta (TGF-β)-responsive genes expressed by fibroblasts, whereas IFN-γ is a potent suppressor of the pro-inflammatory TGF-β signaling pathway, which in turn reduces the production of collagen I synthesis critical to sustain profibrotic processes in the extracellular matrix [17][18][19]. Therefore, our results are in agreement with the role played by those specific cytokines and support the protective role of the IRF8 rs11117432 variant. Thus, IRF8 polymorphisms might have an effect on pro-and anti-inflammatory modulation, along with regulation of the extracellular matrix, and collagen deposition in fibrotic diseases. A recent study of serologic and cytokine profiles in patients with systemic lupus erythematosus also demonstrated an association between an IRF8 variant (rs17445836) and decreased levels of IFIT1 [20]. Moreover, according to HapMap data, the SNPs rs11117432 and rs17445836 are located in a regulatory region, underlying potential functional effects of those SNPs. In SSc, IRF8 SNPs rs11642873 and rs2280381 were already reported as being associated with disease susceptibility. It must be pointed out that rs11117432 and rs2280381 belong to the same block of LD, suggesting a strong association between the IRF8 region and SSc [11,21]. Of interest, we also identified that the association was preferentially observed in the lcSSc subset and its usual associated anti-nuclear auto-antibody (anti-centromere), suggesting a potential role for the pathogenesis in this subgroup of patients.
Epidemiological studies indicate that lcSSc and anticentromere-positive patients are often associated with other AIDs, including secondary Sjögren thyroiditis and PBC, where IFIT1 is known to play a critical role [5,6].
As a complex polygenic AID, SSc is thought to have innumerous genes play a role in its aetiology and heterogeneity. However, only a handful of genes have been identified as SSc risk variants. Moreover, these variants confer relatively small increments in risk and might explain only a small proportion of the trait heritability, suggesting that other loci remain to be discovered. Here, we describe for the first time the association between NF-κB rs7665090 and SSc susceptibility. A previous work suggested that NF-κB rs1598859 had a suggestive association with SSc [10]; our results do support such an association after the analyses of a subset of our sample and by the performance of a meta-analysis combining our results together with the previously reported ones [10]. NF-κB is present in the cytoplasm of virtually all cell types in an active form associated with the inhibitory κB (IκB) and has a key role in primarily immune responses related to antigen recognition, being commonly associated with the development of autoimmune and chronic inflammatory diseases [22]. Of interest, the dysregulation of NF-κB activity has already been implicated in the pathogenesis of systemic lupus erythematosus [23]. In SSc, the functional implications of NF-κB polymorphisms are still unclear. The analysis of NF-κB rs7665090 variant demonstrated no effect on the levels of NF-κB mRNA expression; this is a false-negative result, which may have been due to the small sample size of PBMCs obtained from individuals with SSc. Further studies with a larger sample size might better elucidate the influence of NF-κB rs7665090 variant in gene expression patterns. It is worth mentioning that, in a previous work, our group identified an association between TNIP1 and TNFAIP3 with SSc, genes that are known negative regulators of the NF-κB signaling pathway, suggesting that NF-κB disturbances might have a role in SSc susceptibility [13].
Our single-marker analyses also revealed an association of PLCL2 rs1372072 (OR = 1.23, 95% CI 1.12 to 1.33, P adj = 7.22 × 10 −5 ) with SSc. Phospholipase C (PLC) proteins are important mediators of the calcium-protein kinase C signaling pathway known to be expressed in hematopoietic cells, which play a role in B-cell receptor signaling and other immune responses [24]. The PLCL2 rs4535211 variant has been reported in RA and psoriatic arthritis. The variant was associated with RA risk, whereas in psoriatic arthritis it was found to be protective [25,26]. According to HapMap data, rs1372072 and rs4535211 are in high LD (D′ = 0.90) in a block that contains the PLCL2 gene although not in phase (R 2 = 0.44). Therefore, our findings imply that this locus might contribute to several AIDs; this is an observation that must be confirmed by denser mapping and additional functional experiments to further increase our understanding of the role of PLCL2 in AIDs and SSc. We followed up two other genes (DDX6 and IL12A) for which signals of association with SSc were recently described by another study [12] and because we observed for these 2 genes a nominal association in the first step of our own study for variants harbored by the same genes even if the results of the first step were not confirmed in our combined analyses. For these additional SNPs that are a common variant of DDX6 and a rare variant of IL12A, we could not confirm the signal of association in our sample, ruling out independent replication. When a meta-analysis was performed, no association was observed for DDX6 rs7130875, but the original association observed for IL12A rs77583790 was strengthened, thus confirming this latter locus.
As compared with linkage analyses, association studies have much less power in cases of different mutations acting in different families; however, they are far more sensitive than family studies to study complex polygenic associations. In such diseases, the phenotype is associated with the joint affects of several weakly contributing variants across various loci. However, GWASs have received criticism and skepticism because of the structure of the knowledge it has produced, which contrasts with the highly penetrant Mendelian genetic discoveries. Indeed, interpreting and understanding the molecular mechanisms related to noncoding variants are challenging given the diversity of noncoding functions, the incomplete annotation of regulatory elements, and the probable still uncovered mechanisms of regulatory control. If some studies, through extensive experimental work, uncovered some molecular mechanisms responsible for disease association [27], far more disease-associated variants remain uncharacterized. However, expanding the catalog of non-coding variants associated with human diseases is still meaningful. Indeed, several promises have emerged; for example, the application to GWAS data of gene-set enrichment analyses has provided meaningful data regarding some regulatory mechanisms [28]. The new and complementary approaches, including integration of annotation, higher resolution, and functional maps, may provide the basis for deciphering the effects of non-coding variants in pathways or cell types taking into account the weak additive interactions discovered on GWAS and the potential epistatic interactions between the variants. Such complementary strategies might finally improve the interpretation of disease-associated variants and may help identify the most likely causal regulator. These post-genomic functional studies are eagerly awaited in order to try to elucidate the biological consequences of non-coding variants.
Herein, we identified a novel epistatic interaction between NF-κB and IRF8 genes, which may contribute to SSc susceptibility. This gene-gene interaction might explain the small OR obtained in our genetic association studies, which increase our understanding of missing heritability. Although the biological consequences of this interaction remain unknown, our results highlight the importance of combined analysis of genetic associations in candidate gene studies.
Taken together, our data support the shared genetic bases that are involved in SSc and PBC susceptibility. Of note, we highlight the participation of genes that are related to antigen recognition mediated by TLR and B cells, emphasizing their role in SSc susceptibility. In our mRNA expression analysis, we used resting PBMCs; hence, it is suggested that further studies use PBMC activators to enhance the expression signal.