Abstract

Introduction

Systemic Sclerosis (SSc) is a multi-system disorder of connective tissue characterized by extensive cutaneous and visceral fibrosis. Although the etiopathogenesis of SSc is still unclear, both genetic susceptibility and environmental triggers are widely believed as two major contributors. Genetic association studies have reported a number of SSc-associated genes including HLA genes (e.g. HLA-DQB1, -DPB1 and -DRB1), STAT4, IRF5, CD247, TBX21, BANK1, C8orf13, PTPN22, TNFSF4, FAS, TNFAIP3, CD226, IRAK1, MECP2, MIF, ITGAM, PLD4, TLR-2, CAV1, IL2RA, NLRP1, OPN, PXK, JAZF, KIAA0319, IL-6, IL-21, CXCL8, CSK, PSD3, NFKB1, XRCC1, XRCC4, CCR6, IRF8, GRB10, SOX5, NOTCH4, TNIP1, PSOR1C1, RHOB, DNASEIL3, SCHIP1/IL12A and ATG5 (1–27). On the other hand, studies of environmental hazards contributing to SSc were mainly suggested by occupations with high incidence of SSc. Silica particles have been implicated as an environmental trigger of SSc in several epidemiological studies. Particularly, stonemasons and gold miners who were often exposed to silica particles were more likely to develop SSc (28, 29).

Fibrosis in connective tissues is the most prominent feature of SSc. A major component of connective tissues is fibroblasts. Fibroblasts synthesize, secrete, and maintain all major extracellular matrix (ECM) components that are significantly increased in fibrotic tissues, which directly contribute to fibrosis. Perturbed expression of some major ECM genes have been observed in fibroblasts of SSc patients (30, 31). In particular, COL1A2 and COL3A1 are two major genes encoding structure proteins collagen type I and III, respectively. These collagens are frequently over-expressed in fibroblasts, and are accumulated in fibrotic tissue of SSc patients. Connective tissue growth factor (CTGF) is a profibrotic cytokine that is also over-expressed in SSc patients, and it induces ECM production. Matrix metalloproteinases (MMP) are antifibrotic molecules that induce the degradation of the collagens and other ECM components. Unbalanced expressions of MMP and tissue inhibitor of metalloproteinases (TIMP) are involved in fibrotic process in SSc (30, 31). Our hypothesis is that inherited genetic variations modulate cellular responses to the environment, and in some situations in a pathologic manner. The purpose of this paper is to use a systems-based approach in a fibroblast model that integrates genome-wide genetic variants and the expression of some fibrosis-associated ECM genes to investigate their interactive responses to environmental perturbation.

Methods

Results

Initial analysis

After stringent quality control measures, 183 subjects including 85 cases and 98 controls were included in the final association analyses; see Supplemental Methods for details. Demographic information of the study subjects is shown Table 1. We first analyzed all Caucasians among the subjects (50 SSc cases and 65 controls) as a discovery cohort (Table S1). Using LMM model assuming an additive genetic model, we found several genetic loci with SNPs significantly associated with expression of the ECM genes in this cohort (Table 2). Among them, three previously reported SSc loci including multiple SNPs of TNFAIP3 (6q23) (10, 37), IL2RA (10p15-p14) (15) and ITGAM (16p11.2) (12, 38) showed strong associations with the expression of MMP1 and/or MMP3 of the fibroblasts in responses to dosage-dependent and/or time-course stimulations of silica particles (Table 2). SNPs of CASC9 (8q21.11), intergenic region between LINC00284 and SMIM2-AS1 (13q14.11), FGFR1OP (6q27) and THADA (2p21) also showed associations with the expression of specific ECM genes of the fibroblasts (Table 2), but these genes have not been reported as SSc-associated.

Table 1

Distribution of demographics among SSc cases and controls

Variable	SSc (n = 85) n (%)	Controls (n = 98) n (%)	P-value (χ2)
Age group
≤30	7 (8.24)	25 (25.51)
31–40	14 (16.5)	23 (23.47)
41–50	24 (28.24)	22 (22.45)
>50	40 (47.06)	28 (28.57)	0.0034
Racea
Caucasians	50 (58.82)	65 (66.33)
Hispanics	22 (25.88)	10 (10.20)
African Americans	13 (15.29)	23 (23.47)	0.015
Sex
Female	67 (78.82)	60 (61.22)
Male	18 (21.18)	38 (38.78)	0.016

Open in a separate window

^aSelf-reported race and corrected by STRUCTURE analysis

Table 2

Associations between expression of the ECM genes and genetic backgrounds in Caucasians and mixed populations by meta-analysis.

Peak SNP	Chr	Positiona	Type	Gene(s)b	RA	Pf	P-FE	P-RE	HetISq (%)	HetP	Directionc	Estimates(W|B|H)e	StdErr(W|AA|H)e	MAF(W|AA|H)e	responseType
rs58417815*	8	76159689	ncRNA_intronic	CASC9	C	5.10E-08	-g	-	-	-	+??	1.77|NA|NA	0.29|NA|NA	0.074|NA|NA	COL3A1_dose
rs7823944	8	76182788	ncRNA_intronic	CASC9	T	3.47E-08	2.00E-07	7.10E-02	51.20	0.13	+−+	1.54|−0.03|1.00	0.25|0.73|1.05	0.052|0.033|0.021	COL3A1_dose
rs79411652	6	138138945	intergenic	OLIG3|TNFAIP3d	C	3.87E-11	-	5.69E-01	88.80	0.00	+−?	2.29|−0.80|NA	0.32|0.97|NA	0.039|0.033|NA	MMP1_dose
rs58905141	6	138132123	intergenic	OLIG3|TNFAIP3**	G	3.97E-11	3.46E-11	0.402918	43.1	0.17	++−	2.29|62.41|−2.18	0.32|32.09|65.44	0.039|0.014|0.016	MMP1_dose
rs41290329	10	6054083	UTR3	IL2RA	C	1.51E-08	-	-	-	-	+??	3.43|NA|NA	0.57|NA|NA	0.013|NA|NA	MMP1_dose
rs79365263	13	44618508	intergenic	LINC00284|SMIM2-AS1	C	3.13E-09	-	-	-	-	+??	2.22|NA|NA	0.35|NA|NA	0.024|NA|NA	MMP1_dose
rs79411652	6	138138945	intergenic	OLIG3|TNFAIP3	C	5.55E-18	-	5.12E-01	93.30	0.00	+−?	3.69|−0.91|NA	0.41|0.99|NA	0.039|0.033|NA	MMP3_dose
rs58905141	6	138132123	intergenic	OLIG3|TNFAIP3**	G	6.43E-18	6.00E-18	1.32E-17	0	0.81	+++	3.69|28.75|22.87	0.41|34.61|128.32	0.039|0.014|0.016	MMP3_dose
rs41290329	10	6054083	UTR3	IL2RA	C	1.02E-11	-	-	-	-	+??	5.29|NA|NA	0.74|NA|NA	0.013|NA|NA	MMP3_dose
rs79365263	13	44618508	intergenic	LINC00284|SMIM2-AS1	C	5.33E-08	-	-	-	-	+??	2.80|NA|NA	0.49|NA|NA	0.024|NA|NA	MMP3_dose
rs12926702	16	31240971	intergenic	TRIM72|ITGAM	C	1.33E-10	-	3.42E-01	87.90	0.00	+−?	5.09|−0.35|NA	0.76|1.52|NA	0.014|0.016|NA	MMP3_dose
rs78409037	2	43740411	intronic	THADA	C	4.51E-05	-	9.14E-02	94.10	0.00	+?+	0.97|NA|3.63	0.23|NA|0.55	0.069|NA|0.041	MMP3_timecourse
rs77533229	2	43479638	intronic	THADA**	G	2.68E-08	-	-	-	-	+??	3.35|NA|NA	0.60|NA|NA	0.010|NA|NA	MMP3_timecourse
rs79411652	6	138138945	intergenic	OLIG3|TNFAIP3	C	3.52E-10	-	7.42E-01	93.80	0.00	+−?	1.86|−1.02|NA	0.29|0.64|NA	0.039|0.033|NA	MMP3_timecourse
rs58905141	6	138132123	intergenic	OLIG3|TNFAIP3**	G	3.73E-10	3.23E-10	0.446757	29.1	0.24	++−	1.85|38.02|−11.59	0.29|21.66|53.94	0.039|0.014|0.016	MMP3_timecourse
rs75701002	6	167443918	intronic	FGFR1OP	C	2.18E-08	5.68E-08	5.17E-02	47.50	0.17	+?+	2.39|NA|0.32	0.42|NA|1.32	0.026|NA|0.021	MMP3_timecourse
rs41290329	10	6054083	UTR3	IL2RA	C	1.03E-10	-	-	-	-	+??	3.25|NA|NA	0.47|NA|NA	0.013|NA|NA	MMP3_timecourse

Open in a separate window

RA: risk allele; P: P-value in Caucasians; P-FE: P-value by fixed-effect meta-analysis; P-RE: P-value in random-effect meta-analysis; HetISq: heterogeneity measure I²; HetP: Cochran's heterogeneity test p-valie; Estimates/StdErr/MAF: per risk allele effect estimates/standard error/minor allele frequency in different ethnic groups; W|AA|H: Caucasians|African Americans|Hispanics; responseType: the ECM gene expression in response to silica stimulation in elevated-dosage or time-course.

^aGRCh37/hg19 assembly;

^bGene in which the peak signal is located in, or genes flanking peak signal if intergenic;

^cDirection shown in W|AA|H order, where “?” represents not applicable for that specific ethnic group;

^dBold-face genes are known SSc risk loci;

^e“NA”: not applicable for that specific ethnic group;

^fP-value based on standard normal approximation;

^g”−”: meta-analysis was not performed when only one ethnicity was present.

^*Imputed, otherwise Peak SNP was genotyped;

^**The peak SNP in Caucasians was different from that in fixed-effect meta-analysis with both SNPs listed.

Rs79411652 in the upstream of TNFAIP3 was associated with MMP1 dose-response, MMP3 dose-response and time-course response. The minor allele frequency (MAF) of rs79411652 was 3.9% and 3.3%, respectively, in the Caucasian and African American samples. However, the effect size of this SNP was quite large as shown in Figure 4A, B and C, with 3.6 to 12.7 fold change of gene expression between heterozygous alleles and homozygous major allele. Thus, we had high enough power (over 80%) to detect this association. In addition, regional plots, as shown in Figure 3A, B and C, indicate that there was a cluster of SNPs in high LD with rs79411652 showing association signals, which strengthen the findings about the potential pleiotropic effects of TNFAIP3 on MMP1 and MMP3.

Open in a separate window

Figure 3

Regional plots of association results between previously identified SSc risk loci and gene expression of fibroblasts in the Caucasian cohort

500kb upstream and downstream of the peak SNPs are plotted. Panels A to C: genetic locus at 6q23.3 containing SNPs of TNFAIP3 in association with expression of (A) MMP1 in dose-response assays, (B) MMP3 in dose-response assays, and (C) MMP3 in time-course assay. Panels D to F: genetic locus at 10p15.1 containing SNPs of IL2RA in association with expression of (D) MMP1 in dose-response assays, (E) MMP3 in dose-response assays, and (F) MMP3 in time-course assays. Round dots represent genotyped SNPs while square dots represent imputed SNPs. Purple dots correspond to the peak SNPs. Linkage disequilibrium (r²) with the peak SNP (rs79411652 of TNFAIP3 and rs41290329 of IL2RA, respectively) is shown by different colors.

Open in a separate window

Figure 4

Dose-response and time-course expression trajectories of the ECM genes by genotypes of the peak SNPs in TNFAIP3 and IL2RA in the Caucasian cohort

Panels A-F correspond to panels A-F in Figure 3. Red solid line: sample average trajectory for the heterozygous alleles; blue solid line: sample average trajectory for the homozygous major alleles; gray solid lines: individual trajectories; dashed horizontal grey line at 0 of the y-axis: the reference line representing constant gene expression across time/doses. Y-axis corresponds to the absolute quantity change (AQchange); x-axis corresponds to silica dosage (in µg) or time course (in days).

Figures 3, ​,44 and S3 also show regional plots and genotype-specific expression trajectories for other previously identified SSc loci. SNP rs41290329 located in the 3’ untranslated region (UTR) of IL2RA was significantly associated with MMP1 dose-response, MMP3 dose-response and time-course expression profiles (Figure 3D, E and F) and with 9.5 to 39 fold change of gene expression of MMP1/MMP3 between heterozygous alleles and homozygous major allele (Figure 4D, E and F). Of note, rs41290329 is a low frequency SNP with an MAF of 1.3% in Caucasians and monomorphic in African Americans and Hispanics in our samples as well as in the 1000 Genomes Project. In addition, SNP rs12926702 was significantly associated with MMP3 dose-response expression profile (Figure S3 H) with 34 fold change of gene expression of MMP3 between heterozygous alleles and homozygous major alleles. Rs12926702 is located at the intergenic region of TRIM72 and ITGAM. It is also a low frequency SNP with a MAF of 1.4% and 1.6% in Caucasians and African Americans, respectively, and monomorphic in our Hispanics samples. Given the low MAF of the peak SNPs and the current moderate sample size, we caution that future studies with larger sample size are warranted to confirm the signals identified here and obtain more reliable effect size estimates.

To evaluate whether the association signals of TNFAIP3 and IL2RA could be better explained by multiple associated SNPs in LD or single peak SNPs, we performed haplotype-based analysis of the chromosomal regions defined by 200 SNPs upstream and downstream of the peak SNPs, rs79411652 of TNFAIP3 and rs41290329 of IL2RA. We defined haplotypes based on 3, 5 or 7 consecutive SNPs in sliding windows. The peak signals of haplotype-based analyses agreed well with those of single SNP-based association analyses, although the former yielded no as strong signals as the latter (Figure S2).

In addition to the signals within or near SSc-associated genes, we also identified several genetic loci that have not been reported to be associated with SSc (Table 2), including rs7823944 of CASC9 associated with COL3A1 dose-response expression (Figure S3 A), rs79365263 in the intergenic region of LINC00284 and SMIM2-AS1 associated with MMP1 and MMP3 dose-response expression (Figure S3 D and G), rs78409037 of THADA and rs75701002 of FGFR1OP associated with MMP3 time-course expression (Figure S3 I and K). In particular, rs7823944 of CASC9 was supported by a cluster of imputed SNPs in high LD (r² around 0.80) (Table 2 and Figure S3 A), while rs79365263 of LINC00284|SMIM2-AS1 was clustered with a number of genotyped SNPs in high LD (Figure S3 D and G).

HLA class II genes are well-documented in association with SSc. The nominal association signals (5 × 10⁻⁸ < P < 5 × 10⁻²) of HLA class II genes were observed in the Caucasian cohort. Other SSc-loci did not show significant association with gene expression of the fibroblasts.

As elaborated in Supplemental Methods, the standard normal distribution-based p-values for the SNP effect in the LMM can be too liberal, i.e., not controlling the Type I error at the nominal level. We therefore calculated the more accurate t-distribution-based p-values using the computationally demanding Satterthwaite approximation for the peak SNPs in Table 2. As shown in Table S2, although all of the p-values became larger as expected, the peak SNPs in the SSc-associated loci of TNFAIP3, IL2RA and ITGAM, remained significant.

Meta-analysis of Caucasian, African American and Hispanic cohorts

We conducted meta-analysis of the Caucasian (50 SSc patients and 65 controls), African American (13 SSc patients and 23 controls) and Hispanics (22 SSc patients and 10 controls) cohorts via the fixed-effect inverse-variance weighting method as well as random-effect meta-analysis if determined necessary; see the Supplemental Methods for details. Of note, in the presence of potential heterogeneous genetic effects across the three populations, fixed-effect meta-analysis may not be appropriate. To investigate this, we calculated the heterogeneity measure I² and Cochran’s Q test p-value (Table 2), which indicated that most of genome-wide significant SNPs in the Caucasian cohort, might have heterogeneous genetic effects across different populations (I² > 56%) (39). We, therefore, also performed random-effect meta-analysis for these SNPs. Noticeably, SNP rs58905141 of TNFAIP3 remained genome-wide significantly associated with MMP3 dose-response expression by both fixed-effect and random-effect meta-analysis, even though I² and Cochran’s test results suggested that the random-effect model was not necessary for this SNP. Other SNPs did not show significance with random-effect meta-analysis. The Manhattan plots for fixed-effect and random-effect meta-analysis are, respectively, shown in Figure 1 and Figure S1. Meta-analysis results for all the genome-wide significant SNPs in the discovery analysis of the Caucasian samples are shown in Table S3. The non-significant random-effect meta-analysis results for most of the peak SNPs in Caucasians may be attributed to the low MAF of the SNPs, smaller sample sizes in the African American and Hispanic cohorts and disparate LD patterns across populations, leading to diminished statistical power.

Open in a separate window

Figure 1

Manhattan plots of fixed-effect meta-analysis association p-values for dose-response and time-course ECM gene expressions in fibroblasts in response to silica stimulation

A. MMP1 gene expression in silica dose response showed significant signals at 6q23.3 (TNFAIP3), 10p15.1 (IL2RA) and 13q14.11 (LINC00284). B. MMP3 gene expression in silica dose response showed significant signals at 6q23.3 (TNFAIP3), 10p15.1 (IL2RA), 13q14.11 (LINC00284) and 16p11.2 (TRIM72|ITGAM). C. MMP3 gene expression in silica time course assays showed significant signals at 2p21(THADA), 6q23.3 (TNFAIP3) and 6q27(FGFR10P|CCR6). Red line corresponds to the p-value cutoff for genome-wide significance (5 × 10⁻⁸); blue line corresponds to the suggestive significance cutoff (1 × 10⁻⁵). P-values from Caucasian-only analysis are displayed if fixed-effect meta-analysis was not performed due to missing association results in African Americans and Hispanics.

Discussion

SSc is a fibrotic disease with complex genetic traits, to which environmental factors may trigger a pathological process toward systemic fibrosis. Functional association studies of genetic polymorphisms in connection with environmental triggers are challenge. It is unknown whether the SSc susceptibility loci contribute to specific fibrotic process in the presence of an environmental hazard. In comparison with conventional gene by environment (GxE) analysis based on observational studies (40, 41), the studies herein are novel in that we examined profibrotic effects of silica particles in a complex cell model of human fibroblasts that were genetically profiled with Immuno-chip SNPs.

A previously reported SSc locus at 6q23.3 with SNPs of TNFAIP3 (10, 37) appeared to be persistently in association with expression of MMP1 and/or MMP3 genes of the fibroblasts in responses to dosage-dependent and time-course stimulations of silica particles in both Caucasian only analysis and meta-analysis of diverse populations (Table 2). TNFAIP3 stands for tumor necrosis factor-α-induced-protein 3. It encodes A20, an ubiquitin-modifying enzyme that inhibits NF-kB activity and is a key regulator of TNF-mediated immune responses and inflammation (42, 43). Inflammation has been reported as a major consequence induced by silica particles in in vivo and in vitro studies (35, 44). Our previous studies indicated that silica might activate fibroblasts to overexpress the ECM genes through proinflammatory process (32, 35). Therefore, the results suggest that an interaction between genetic factors of TNFAIP3 and environmental factor of silica stimulation may control fibroblasts to express MMP1/MMP3 genes, which may be associated with development of SSc through an inflammatory and immune-control mechanism. Of note, intronic SNP rs5029939 in TNFAIP3 was reported to be associated with SSc in a previous case-control study (10). LD analysis by the SNAP tool (45) indicated that this SNP was in perfect LD (r² = 1 and D’=1) with our top hit SNP rs79411652 in the discovery analysis of Caucasian samples, and in high LD with our meta-analysis top hit SNP rs58905141 (r² = 0.80 and D’=1) associated with MMP3 dose-response and time-course expression as well as MMP1 dose-response expression. Regional plots (Figure 3A, C and E) also suggest that these SNPs were likely to tag the same functional locus/loci. We annotated the top SNPs of TNFAIP3 with RegulomeDB (46) according to the ENCODE database (47). While rs79411652 was annotated to category 7, i.e., no known biological function, rs58905141 was annotated to the second highest functional category “2b”, which is likely to affect binding of transcription factors FOXA1 and SMARCA4 by ChIP-seq experiments, binding of STAT1 by motif analysis, and is in DNase footprint by DNase-seq and in DNase peak by ChIP-seq experiments. This suggests that rs58905141 may play a functional role in regulating TNFAIP3 expression and directly or indirectly regulate MMP1/MMP3 expression in fibroblasts. Further replication and functional studies of the signals in TNFAIP3 are warranted to confirm this hypothesis.

Two other reported SSc-associated genetic loci with SNPs of IL2RA (10p15-p14) and ITGAM (16p11.2) genes were significantly associated with the expression of MMP1 and/or MMP3 in silica response assays in Caucasians (Table 2). IL2RA stands for interleukin 2 receptor α. IL2RA is involved in various pathways that help control the differentiation of effector cells, T-cell proliferation, and immune tolerance (48). A meta-analysis study with a large Caucasian European cohort revealed an association between rs2104286 of IL2RA and ACA positive SSc patients (15). LD analysis using the SNAP software showed that MMP1/MMP3 expression associated rs41290329 was in complete LD with rs2104286 (D’ = 1), but in low LD as measured by r² (=0.06), due to the disparate MAFs of the two SNPs (0.013 and 0.225 in the 1000 Genomes CEU population, respectively). Functional analysis by RegulomeDB indicated that biological function for rs41290329 was unknown. As for ITGAM that encodes the α subunit of the αMβ2-integrin, it has been recently identified as an autoimmune disease risk gene (49). It is expressed on the surface of leukocytes, and it regulates adhesion of neutrophils and monocytes, cell activation, which is important for innate immunity (50). Rs1143679 of ITGAM was associated with SSc in European cohorts (12) and a large-scale meta-analysis (38). LD analysis suggested that MMP3 dose-response expression associated rs12926702 was in complete LD with rs1143679 (D’ = 1), but in low LD as measured by r² (=0.003), due to the low MAFs of both SNPs (0.033 and 0.083 in the 1000 Genomes CEU population, respectively). Considering the low MAFs of rs41290329 and rs12926702, in addition to our moderate sample size, especially for non-Caucasian samples, we caution that the findings on IL2RA and ITGAM need to be confirmed in future studies with larger sample size.

In addition, several other genetic loci and genes unrelated to SSc susceptibility were also strongly associated with expression of the ECM genes of silica-stimulated fibroblasts in Caucasian-only analysis and/or meta-analysis (Table 2). Whether these associations represent a general impact of specific genetic loci on immune-mediated diseases is worth further investigation.

It is worth noting that many of those known SSc-loci did not show significant impact on gene expression of the fibroblasts in this study. For instance, the HLA class II genes that confer the strongest susceptibility to SSc (1) showed only nominal association signals (see Figure 2), which suggest that potential genetic impact of these genes to profibrotic responses of silica-stimulated fibroblasts may be less significant in this study model. It is likely that genetic impact of specific SSc-loci may be determined by interactions of specific types of cells and environmental triggers, as well as specific responding genes.

Open in a separate window

Figure 2

Chromosome 6 Manhattan plots of fixed-effect meta-analysis association p-values for dose-response and time-course ECM gene expressions in fibroblasts in response to silica stimulation

Panels A, B and C correspond to Figure 1 A, B and C with a zoom in for chromosome 6 only. Green dots show loci on or near known SSc risk genes. Red line corresponds to the p-value cutoff for genome-wide significance (5 × 10⁻⁸); blue line corresponds to the suggestive significance cutoff (1 × 10⁻⁵). P-values from Caucasian-only analysis are displayed if fixed-effect meta-analysis was not performed due to missing association results in African Americans and Hispanics.

The current study employed the cost-effective Immunochip targeted at 186 gene regions which were identified in previous association studies of autoimmune diseases. Although we successfully performed 1000 Genomes-based imputation to boost the 130,000 directly genotyped SNPs to around 800,000 total SNPs, the coverage was likely to be very low beyond the 186 gene regions. As a result, we cannot exclude the possibility that other genetic variants that were not investigated in the current study might be associated with the ECM gene expression in the fibroblasts. In addition, heterogeneity of patients’ clinical status is another concern. This warrants further investigation in the future.

In summary, this is the first attempt to study interactions between genome-wide genetic variants and environmental factors of SSc in a complex fibroblast model. The results indicated that a previously identified SSc locus of TNFAIP3 was strongly and persistently associated with silica-induced profibrotic responses of the fibroblasts in both Caucasians and meta-analysis of mixed populations. Similar associations were observed in two other reported SSc loci of IL2RA and ITGAM in Caucasians. In addition, some non-SSc loci may also impact a general response of fibroblasts to silica particles. SSc subtype analysis of the associations was limited by small sample size. Therefore, this notion may be further verified in other cohorts with larger sample size and in functional studies beyond the ECM genes.

Supplementary Material

Acknowledgments

Footnotes

References