A 3-factor epistatic model predicts digital ulcers in Italian scleroderma patients

European Journal of Internal Medicine 21 (2010) 347–353

Contents lists available at ScienceDirect

European Journal of Internal Medicine j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / e j i m

Original article

A 3-factor epistatic model predicts digital ulcers in Italian scleroderma patients Lorenzo Beretta a,⁎, Alessandro Santaniello a, Michael Mayo c, Francesca Cappiello a, Maurizio Marchini a,b, Raffaella Scorza a,b a b c

Referral Center for Systemic Autoimmune Diseases, Fondazione, "Ca' granda" Ospedale Maggiore Policlinico Milano, Italy University of Milan, Italy Department of Computer Science, University of Waikato, Hamilton, New Zealand

a r t i c l e

i n f o

Article history: Received 19 March 2010 Received in revised form 23 April 2010 Accepted 19 May 2010 Available online 23 June 2010 Keywords: Cytokine HLA Petri nets Epistasis Systemic sclerosis Digital ulcers

a b s t r a c t Background: The genetic background may predispose systemic sclerosis (SSc) patients to the development of digital ulcers (DUs). Methods: Twenty-two functional cytokine single nucleotide polymorphisms (SNPs) and 3 HLA class I and II antigens were typed at the genomic level by polymerase chain reaction in 200 Italian SSc patients. Associations with DUs were sought by parametric models and with the Multifactor Dimensionality Reduction (MDR) algorithm to depict the presence of epistasis. Biological models consistent with MDR results were built by means of Petri nets to describe the metabolic significance of our findings. Results: On the exploratory analysis, the diffuse cutaneous subset (dcSSc) was the only single factor statistically associated with DUs (p = 0.045, ns after Bonferroni correction). Gene–gene analysis showed that a 3-factor model comprising the IL-6 C-174G, the IL-2 G-330T SNPs and the HLA-B*3501 allele was predictive for the occurrence of DUs in our population (testing accuracy = 66.9%; p b 0.0001, permutation testing). Conclusion: Biological interpretation via Petri net showed that IL-6 is a key factor in determining DUs occurrence and that this cytokines may synergise with HLA-B*3501 to determine DUs onset. Owing to the limited number of patients included in the study, future research are needed to replicate our statistical findings as well as to better determine their functional meaning. © 2010 European Federation of Internal Medicine. Published by Elsevier B.V. All rights reserved.

1. Introduction Digital ulcers (DUs) are a frequent complication of systemic sclerosis (SSc) that may occur in up to 50% of patients during their disease history [1–4]. The occurrence of DUs has been associated with several clinical features of the disease in different populations, such as a younger age at onset of the disease [2–6], the male gender [2], the diffuse cutaneous subset (dcSSc) [2,6,7], the anti-topoisomerase I antibody [2,6], the presence of pulmonary hypertension [6], the presence of systemic inflammation [6,7] or a delay in introducing vasodilating therapy with calcium channel-blockers or iloprost (4). Whatever the association with other clinical features, ischemic DUs are thought to be the consequence of endothelial injury and small-vessel vasculopathy [8]. Yet, autoimmunity with involvement of T lymphocytes and mcrophages and the local release of cytokines with pro-inflammatory, as well as chemoattractant, properties are also relevant in the early onset of endothelial injury and thus in the development of DUs. Indeed, it has been shown that interleukin-6 (IL-6) plasma levels are increased in patients with SSc ⁎ Corresponding author. Referral Center for Systemic Autoimmune Diseases, Fondazione IRCCS Ospedale Maggiore Policlinico, Mangiagalli e Regina Elena, Via Pace 9, 20122 Milano, Italy. Tel.: +39 02 55035282; fax: +39 02 55035289. E-mail address: [email protected] (L. Beretta).

and DUs [7], and that this group of patients has increased plasma soluble CD40 ligand concentrations [9], reflecting both an inflammatory state and lymphocyte T activation. Also, endothelin-1 (ET-1) blockade by bosentan was shown effective in preventing the recurrence of new DUs in SSc, confirming, in this context, a pathophysiological role for ET-1 which is actively released by the damaged endothelium [10,11]. The secretion and/or bioavailability of cytokines with immunomodulatory, pro-fibrotic or pro-inflammatory function is regulated at the genetic level [12–16], hence it could be hypothesized that single nucleotide polymorphisms (SNPs) in or near cytokine genes may be relevant in determining vasculopathy and thus the appraisal of DUs in SSc patients. This hypothesis was verified in the present study as part of a project aimed at finding a possible role for genetic variants with well-known regulatory functions on cytokine production [17] and SSc or its clinical expressions [18–22]. Besides the possible association with DU and 22 cytokine gene SNPs, we also sought a possible association with the human leukocyte antigen system (HLA) allele, HLA-B*3501, that was previously shown to regulate in vitro ET-1 secretion [23] and endothelial cell apoptosis [24], and that was also associated with other vascular complications of SSc, such as isolated pulmonary hypertension [25]. In light of the multifactorial pathogenesis of SSc, we searched for high-order non-linear gene–gene interactions (epistasis) and built models of biochemical and

0953-6205/$ – see front matter © 2010 European Federation of Internal Medicine. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.ejim.2010.05.010

348

L. Beretta et al. / European Journal of Internal Medicine 21 (2010) 347–353

physiological networks that are consistent with high-dimensional models of disease susceptibility to better explain these associations. 2. Methods 2.1. Patients selection The patients used in this study consisted of two-hundred unrelated Italian SSc patients who were referred to our outpatient clinic. They provided written consent to have their DNA collected as well as their clinical data recorded and utilised for medical research. The vast majority of the patients (n = 181, 91%) fulfilled the classification criteria proposed by the American College of Rheumatology [26], yet we also considered patients (n = 19, 9%) with definite SSc who do not fulfil these criteria [1]. Patients were categorised as having Early SSc, the “sine scleroderma,” the limited cutaneous (lcSSc) or the dcSSc subset [1,27]. History of DUs was ascertained by reviewing the patients' medical records. Disease onset was determined by the patient's recall of the first non-Raynaud symptom clearly attributable to scleroderma [28] and were categorised into 3 groups according to the age at onset: ≤ 30 years, 31–50 years, N50 years. Restrictive lung disease was defined as a forced vital capacity (FVC) ≤70% of predicted, with a normal Tiffenau's index [1]; diffusing capacity for carbon monoxide (DLco) impairment was considered for values ≤70% of predicted [1]. Pulmonary hypertension was defined as a right-ventricular systolic pressure (RVSP) ≥ 40 mm Hg and was confirmed by right-heart catheterisation. Antinuclear antibodies (ANA) were determined by indirect immunofluorescence on Hep2 cells (Kallestad, Chaska, MN) using a standardised technique [29]; extractable nuclear antigens (ENAs) were determined by a commercial enzyme-linked immunoassay (ELISA) (Diamedix, Miami, FL). 2.2. Cytokine SNP determination Blood samples were collected in citrate and DNA was extracted using the DNA Isolation Kit For Mammalian Blood (Roche Diagnostics, Indianapolis, IN). Genotyping was performed by polymerase chain reaction (PCR) with sequence-specific primers (PCR-SSP) according to the 13th International Histocompatibility Workshop recommendations, and by using the commercial Cytokine kit provided by the University Clinic of Heidelberg (CTS-PCR-SSP TRAY, from the Institute of Immunology, Department of Transplantation Immunology, University of Heidelberg, Heidelberg, Germany) as previously described [21]. Overall, the following 22 interleukin (IL) SNPs were analysed: IL1 alpha (IL-1α) C-889T, IL-1 beta (IL-1β) C-511T and IL-1β C+3962T, IL-1 receptor (IL-1R) Cpst1970T, IL-1 receptor antagonist (IL-1Ra) Cmspal11100T, IL-2 G-330T, IL-2 G+160T, IL-4 G-1098T, IL-4 C-590T, IL-4 C-33T, IL-4 receptor antagonist alpha (IL-4Rα) A+1902G, IL-6 C174G, IL-6 Ant565G, IL-10 A-1082G, IL-10 C-819T, IL-10 A-590C, IL-12 A-1188C, transforming growth factor beta (TGF-β1) T/C codon 10, TGF-β1 G/C codon 25, interferon gamma (IFNγ) AUTR5644T, tumor necrosis factor alpha (TNFα) A-308G, TNFα A-238G. Typing results for the IL-4 G-1098T, IL-4 C-590T, IL-4 C-33T did not meet the quality requirements for interpretation (e.g. because of unequal or weak amplification results) and were therefore excluded from analysis. 2.3. HLA genotyping The following HLA class I and II antigens were considered for analysis: HLA*B3501, that was previously associated with an increased production of ET-1 or cell apoptosis in vitro and ex vivo [23,24], HLA-DR*11 and HLA-DR*07¸which were reported to be associated with SSc in Italian subjects [30].

HLA class I and class II antigens were typed at genomic level by polymerase chain reaction–sequence-specific oligonucleotide probes as previously described elsewhere [31,32]. Products of the PCR reaction were purified using Microcom columns (Amicon, Beverly, MA, USA), while products of the sequencing reaction were purified using CentriSep Spin Columns (Applied Biosystems, Monza, Italy). The sequencing of exons was done with BigDye Primers (Applied Biosystems). The sequencing reaction products were resuspended in 6 ml of formamide buffer (5:1), from which 1.75 ml was taken to be electrophoresed on 5% polyacrylamide (Applied Biosystems) 6 M urea gels with an ABI PRISM 377 DNA sequencer. DRB1 locus was typed by BigDye Terminator SBT Typing Kit (Applied Biosystems).

2.4. Statistical analysis The distribution of genotypes was tested for Hardy–Weinberg equilibrium (HWE) with the goodness-of-fit χ2 test at a significance level of 0.05. The distribution of clinical features, cytokine SNP variations, and of HLA alleles in SSc patients with or without DU was tested by the χ2 test or Fisher's exact test when necessary. Variation in a particular SNP was considered to be associated with the endpoint at a significance level of 0.05 after Bonferroni correction.

2.5. MDR The evaluation of gene–gene interactions was performed using the Multifactor Dimensionality Reduction (MDR) algorithm [33,34]. Using the MDR algorithm we constructed a series of combinations of twoto-four variables and then used a naïve Bayes classifier in the context of 10-fold cross-validation to estimate the testing accuracy (TA) of each best two-to-four factor model. A single best model was selected that maximized the TA. This is the model that is most likely to generalize to independent datasets. Statistical significance of the final model was evaluated using a 1000-fold permutation test to compare observed testing accuracies with those expected under the null hypothesis of no association. Permutation testing corrects for multiple testing by repeating the entire analysis on 10,000 datasets that are

Table 1 Clinical and demographic characteristics. Demographic and clinical characteristics of 200 Italian systemic sclerosis (SSc) patients with or without a personal history of digital ulcers (DU). lcSSc, limited SSc; dcSSc, diffuse SSc; ANA, antinuclear antibodies; ACA, anticentromere antibodies; Scl70, anti-topoisomerase I antibody; FVC, forced vital capacity; DLco, diffusing capacity for carbon monoxide; PAH, pulmonary hypertension (right-heart catheterisation). Variable Subset, n (%) Early SSc Sine SSc lcSSc dcSSc⁎ Autoantibodies, n (%) ANA ACA Scl70 Females, n (%) Age at onset, n (%) ≤30 years 31–50 years N50 years Disease duration, years FVC ≤ 70% predicted, n (%) DLco ≤ 70% predicted, n (%) PAH, n (%)

SSc (n = 200) 14 (7) 4 (2) 137 (68.5) 45 (22.5) 195 82 81 179

(97.5) (42) (40) (89.5)

17 (8.5) 70 (35) 97 (48.5) 14 ± 3.1 49 (24.5) 79 (39.5) 17 (8.5)

DU (n = 111) 5 (4.5) 2 (1.8) 73 (65.8) 31 (27.9)

No DU (n = 89) 9 (10.1) 2 (2.3) 64 (71.9) 14 (15.7)

107 (96) 44 (40) 44 (40) 97 (87)

88 40 37 82

(99) (45) (42) (92)

15 (13.5) 41 (36.9) 49 (49.6) 14.8 ± 2.7 29 (26) 46 (41) 11 (9.9)

12 (13.5) 29 (32.6) 48 (53.9) 13.5 ± 3.6 20 (22) 33 (37) 6 (6.7)

⁎ p = 0.045 vs the other disease subtypes (not significant after Bonferroni correction).

L. Beretta et al. / European Journal of Internal Medicine 21 (2010) 347–353

consistent with the null hypothesis. Models were considered significant at the 0.05 level. To perform MDR analysis the dataset was pre-processed by the Tuned ReliefF (TurfF) algorithm: the most interesting attributes were selected by the TurfF algorithm [35], removing those variables with TurfF scores close or below to zero, that is those attributes deemed irrelevant or noisy to class prediction. Only the unfiltered attributes were then used to build the interaction models [21]. To perform MDR analysis, HLAs were considered at the allelic level and heterozygous and homozygous individuals were grouped together to reduce the possibility of including multilocus cells with

349

few or no cases in high-order epistatic models, thus increasing the chance of false-positive findings [8]. MDR analysis was both carried out to build gene–gene interaction models as well as to test for interactions with clinical features; in the latter cases all the studied genetic variants and the patients' clinical characteristics were considered in the searching strategy. Finally, we used measures of interaction information to provide a statistical interpretation of gene–gene interactions models and dendograms. To conduct the analyses the open-source MDR software package (v.2.0) and the MDR permutation tool available from www.epistasis.org, were used. The interaction information

Table 2 Cytokine SNPs genotype and allele distribution. Allele and genotype frequency of 19 cytokine single nucelotide polymorphisms (SNPs) in 200 Italian systemic sclerosis (SSc) patients with or without digital ulcers (DU). None of the SNPs was associated with DU (chi-square test) with Bonferroni adjustment for multiple comparisons. Genetic variant

Position

IL-1α

−889

IL-1β

SSc (n = 400)

DU (n = 222)

No DU (n = 178)

n (%)

n (%)

n (%)

C T

310 (77.5) 90 (22.5)

180 (81.2) 42 (18.5)

130 (72.5) 49 (27.5)

−511

C T

266 (66) 134 (34)

149 (67.1) 73 (32.9)

117 (64.6) 61 (35.4)

+ 3962

C T

320 (80) 80 (20)

179 (80.6) 43 (19.4)

141 (79.29) 37 (20.8)

IL-1R

pst11970

C T

283 (70.8) 117 (29.2)

157 (70.7) 65 (29.3)

126 (70.8) 52 (29.2)

IL-1RA

Mspa111100

C T

85 (21.3) 315 (78.8)

52 (23.4) 170 (76.6)

33 (18.5) 145 (81.5)

IL-4RA

+ 1902

A G

335 (83.8) 65 (16.2)

182 (82) 40 (18)

153 (86) 25 (14)

IL-12

−1188

A C

293 (73.3) 107 (26.7)

159 (71.6) 63 (28.4)

134 (75.3) 44 (24.7)

IFNγ

UTR5644

A T

219 (54.8) 181 (45.2)

118 (53.2) 104 (46.8)

101 (56.7) 77 (43.3)

TGF-β1

Codon 10

C T

202 (50.5) 198 (49.5)

118 (53.2) 104 (46.8)

84 (47.2) 94 (52.8)

Codon 25

C G

41 (10.3) 359 (89.8)

25 (11.3) 197 (88.7)

16 (9) 162 (91)

−308

A G

33 (8.3) 267 (91.8)

21 (9.5) 201 (90.5)

12 (6.7) 166 (93.3)

−238

A G

26 (6.5) 374 (93.5)

17 (7.7) 205 (92.3)

9 (5.1) 169 (94.9)

−330

G T

137 (34.3) 263 (65.7)

79 (35.6) 143 (64.4)

58 (32.6) 120 (67.4)

+ 160

G T

254 (63.5) 146 (36.5)

143 (64.4) 79 (35.6)

111 (62.4) 67 (37.6)

−174

C G

123 (30.8) 277 (69.2)

60 (27) 162 (73)

63 (35.4) 115 (34.6)

nt565

A G

105 (26.3) 295 (73.7)

53 (23.9) 169 (76.1)

52 (29.2) 126 (70.8)

−1082

A G

235 (58.8) 165 (41.3)

124 (55.9) 98 (44.1)

111 (62.4) 67 (37.6)

−819

C T

319 (79.8) 81 (20.2)

186 (83.8) 36 (16.2)

133 (74.7) 45 (25.3)

−590

A C

81 (20.3) 319 (79.7)

36 (16.2) 186 (83.8)

45 (25.3) 133 (74.7)

TNFα

IL-2

IL-6

IL-10

Allele

Genotype

CC CT TT CC CT TT CC CT TT CC CT TT CC CT TT AA AG GG AA AC CC AA AT TT CC CT TT CC CG GG AA AG GG AA AG GG GG GT TT GG GT TT CC CG GG AA AG GG AA AG GG CC CT TT AA AC CC

SSc (n = 200)

DU (n = 111)

No DU (n = 89)

n (%)

n (%)

n (%)

120 70 10 83 98 19 131 58 11 100 83 17 5 75 120 144 47 9 103 87 10 63 93 44 47 108 45 2 37 161 0 33 167 0 26 174 17 103 80 105 44 51 17 89 94 10 85 105 6 103 31 127 65 8 8 65 127

72 37 2 46 57 8 74 31 6 54 49 8 4 44 63 77 28 6 54 51 6 33 52 26 30 58 23 2 21 88 0 21 90 0 17 94 8 63 40 57 23 28 8 44 59 5 43 63 34 56 21 79 28 4 4 28 79

48 (53.9) 33 (37.1) 8 (9) 37 (41.6) 41 (46.1) 11 (12.4) 57 (64) 27 (30.3) 5 (5.6) 46 (51.7) 34 (38.2) 9 (10.1) 1 (1.1) 31 (34.8) 57 (64) 67 (75.3) 19 (21.3) 3 (3.4) 49 (55.1) 36 (40.4) 4 (4.5) 30 (33.7) 41 (46.1) 18 (20.2) 17 (19.1) 50 (56.2) 22 (24.7) 0 (0) 16 (18) 73 (82) 0 (0) 12 (13.5) 77 (86.5) 0 (0) 9 (10.1) 80 (89.9) 9 (10.1) 40 (44.9) 40 (44.9) 45 (50.6) 21 (23.6) 23 (25.8) 9 (10.1) 45 (50.6) 35 (39.3) 5 (5.6) 42 (47.2) 42 (47.2) 32 (36) 47 (52.8) 10 (11.2) 48 (53.9) 37 (41.6) 4 (4.5) 4 (4.5) 37 (41.6) 48 (53.9)

(60) (35) (5) (41.5) (49) (9.5) (65.5) (29) (5.5) (50) (41.5) (8.5) (2.5) (37.5) (60) (72) (23.5) (4.5) (51.5) (43.5) (5) (31.5) (46.5) (22) (23.5) (54) (22.5) (1) (18.5) (80.5) (0) (16.5) (83.5) (0) (13) (87) (8.5) (51.5) (40) (52.5) (22) (25.5) (8.5) (44.5) (47) (5) (42.5) (52.5) (3) (51.5) (15.5) (63.5) (32.5) (4) (4) (32.5) (63.5)

(64.9) (33.3) (1.8) (41.4) (51.4) (7.2) (66.7) (27.9) (5.4) (48.6) (44.1) (7.2) (3.6) (39.6) (56.8) (69.4) (25.2) (5.4) (48.6) (45.9) (5.4) (29.7) (46.8) (23.4) (27) (52.3) (20.7) (1.8) (18.9) (79.3) (0) (18.9) (81.1) (0) (15.3) (84.7) (7.2) (56.8) (36) (51.1) (20.7) (25.2) (7.2) (39.6) (53.2) (4.5) (38.7) (56.8) (30.6) (50.5) (18.9) (71.2) (25.2) (3.6) (3.6) (25.2) (71.2)

350

L. Beretta et al. / European Journal of Internal Medicine 21 (2010) 347–353

measures that we used are all implemented in this MDR software package. 2.6. Petri nets Whilst MDR is a powerful strategy to detect statistical epistasis, it is difficult to dissect the true biological significance of high-order gene–gene interactions found at the population level. The use of discrete dynamic system tools called Petri nets (PNs) has thus been advocated as a strategy for developing plausible biochemical or metabolic systems that are consistent with genetic models of disease susceptibility [36,37]. Briefly, to generate deterministic PNs consistent with MDR models, an extension of the multi-start random Hill-Climbing approach described by Mayo [37] was used. This algorithm takes as input a set of genotype combinations. Each genotype also has an associated output, being the presence of absence of some disease (in this case, DUs). The algorithm then attempts to find the simplest Petri net that “explains” the presence or absence of the disease for each genotype combination, and a Petri net that explains all of the observed genotype/disease correlations is referred to as a “perfect” net. The following rules were applied in the search strategy: a) given a “gene unit” as described in [38], the weight of the arc leading to the “gene-product place” is genotype-dependent. All the other components of the net, such as the number of places and arc weights are determined by the search strategy; b) the weight of genotypedependent arcs is randomly determined during the searching strategy such that, given a quantity n assigned to the wild-type allele we have the following genotype-dependent assignments: wild-type = 2 ⁎ n, heterozygous = n + (x ⁎ n) and homozygous = 2 ⁎ (x ⁎ n), where x = 2, 3 or 4. These numbers roughly correspond to cases were the mutant allele is associated with a 2-to-4-fold increase in the gene product and

are consistent with several biological observations [13,14]; c) the maximum capacity of a place is set to 10 tokens, and otherwise enabled transitions that may increase the place capacity over the maximum allowed value are inhibited from firing; d) given an initial genotype/marking, the network is allowed to “run” until no further transitions are enabled; at this point, the number of tokens at a special place representing the “toxic reaction” (i.e. the accumulation of substances favouring the onset of digital ulcers) is measured, and if this quantity exceeds or equals 9 token (being 90% of the maximum capacity of a place) the presence of the disease is assumed; e) Petri nets that may require more than 100 transitions to fire for any genotype are automatically excluded from the search; f) when multiple transitions are enabled to fire, an arbitrary delay ordering of transitions is assumed, that is, each transition Tm has a unique index m, such that transition Tm will always fire before transition Tn whenever m b n; g) no infinite source or infinite sinks of tokens are allowed; and h) when multiple perfect models are found, the most parsimonious one, i.e. that is the one with the lowest number of arcs and firing transitions, is chosen. 3. Results The clinical and the demographic characteristics of the patients are reported in Table 1. The presence of digital ulcers was statistically associated with the dcSSc subset on the exploratory analysis (p = 0.045 vs the other subtypes). However, this association did not stand after correction for multiple comparisons. Sixty-three (31.5%) patients carried the HLA-B*3501 allele, 130 (65%) patients the HLA-DR*11 allele, and 28 (14%) patients the HLADR*07 allele. The prevalence of HLA class I and II alleles was evenly distributed in SSc patents with or without DU (HLA-B*3501: 34% vs 30%, HLA-DR*11: 69% vs 62% and HLA-DR*07: 12% vs 15%,

Fig. 1. Multifactor Dimensionality Reduction (MDR) model. Top panel. Multidimensional cells resulting from the interaction of the HLA-B*3501 allele, the IL-6 C-174G and IL-2 G330T single nucleotide polymorphisms (SNPs). Left bars represent the number of patients with digital ulcers (DUs), right bars the number of patients without DUs. High-risk genotype combinations are depicted as dark-shaded cells. The scattered distribution of high-risk cells across the multidimensional space is typical of non-linear epistatic interactions. Bottom panel. Dendrogram resulting from information theory analysis. The IL-6 and the IL-2 SNPs synergistically interact to determine the occurrence of DUs, this bi-factorial component further interacts in a non-linear manner with HLA-B*3501 to explain the presence of DUs.

L. Beretta et al. / European Journal of Internal Medicine 21 (2010) 347–353

respectively). The allele and genotype frequencies for the studied cytokine SNPs in SSc patients with or without DU are reported in Table 2. No statistical differences in allele or genotype distribution were observed between the two groups of patients either. The best MDR model (Fig. 1, top panel) included: the IL-2 G-330T, the IL-6 C-174G SNPs, and the HLA-B*3501 allele; overall, this model had a TA = 66.9% (sensitivity = 69.4%, specificity = 64.4%) and was highly significant after 10,000-fold permutation test (p b 0.001). The intermediate models generated before reaching the best 3-way model were the following: 1-way, IL-6 C-174G, TA = 52.6%; 2-way, IL-6 C174G and IL-2 G-330T, TA = 58.7%. Dendograms plotted in Fig. 1, bottom panel, better illustrate the relationship among the genetic variations included in the best MDR model. As it can be observed, the IL-2 G-330T and the IL-6 C-174G SNPs show a strong synergic effect that further interacts in a multiplicative manner with HLA-B*3501 to determine the disease status. When both the genetic and the clinical variables were considered together no better model was found. Only the 3-way genetic interaction model was thus deemed to be associated with DU in SSc patients. The Petri model consistent with MDR analysis results is shown in Fig. 2a. This net yielded 100% accuracy in assigning cases to high- or low-risk cells on the basis of the IL-2 G-330T, the IL-6 C-174G genotypes and the HLA-B*3501 allele. As it can be observed, the IL-6 C-173G SNP the HLA-B*3501 allele converge into a biochemical reaction/pathway to determine the occurrence of DUs. Yet, as it can be appreciated in Fig. 2b, IL-6 appears to be a key factor within this network. Indeed, the activation of its gene and the production of this cytokine is, in most cases, sufficient to activate HLA-B*3501 responses. Finally, an IL-6 amplifying loop may further sustain the phenomena that lead to DUs occurrence; this loop can be further modulated by IL2. All these aspects are consistent with information theory (dendograms) analysis and would confirm the correctness of the Petri net model.

351

The biological interpretation of statistical models of epistasis is a challenging task as what is found at the population level does not necessarily reflect what ultimately happens at the biological level [36]. To solve this issue we applied a discrete dynamic system tool, specifically Petri nets, to provide a plausible biological interpretation of the MDR model. The most interesting aspect that emerge from the model we describe (Fig. 2) is the major effect played by IL-6 in DUs susceptibility. Notably, to date only a few studies have sought for DUs biological markers in scleroderma subjects [7,9] and among these, IL-6 has emerged as relevant in Italian patients. The effect of IL-6 appears to be partially sustained and amplified by an autocrine loop, which has already been described in autoimmune diseases [41]. It has also been observed that in autoimmune disorders other cytokines may act

4. Discussion DUs are a frequent complication of SSc severely affecting wellbeing and contributing significant morbidity of affected patients [39,40]. Several authors have sought for factors that may predispose or be associated with the onset of DU in scleroderma patients [2–8]. Many pathophysiological processes that underlie their occurrence have been depicted [8,10,11]. Yet, patients otherwise similar as far as clinical and demographical characteristics as well as exposure to environmental factors, may have a completely different burden of DUs and may or may not experience the appraisal of DUs during their disease history. In the present study, we demonstrated that the individual genetic background may be — at least to some extent — responsible for this clinical variability regardless from the patient's clinical characteristics. Indeed, the genetic association model we describe is not influenced by clinical factors otherwise linked the occurrence of DUs, such as the dcSSc subset [2,5,7], autoantibodies [2,5], the age at onset of the disease [2,4,6] or the presence of pulmonary hypertension [5]. Complex human diseases as well as their variable phenotypic expression are likely to arise from a plethora of interacting factors. Thus it is not surprising that among the genetic variables we studied, none was singularly associated with DUs in Italian SSc patients, but rather it was their non-linear (epistatic) interaction to be relevant to DUs occurrence. The finding that our interaction model, which results from the study of a limited number of genetic variants retains a relatively high predictive accuracy is not unexpected either. Yet, whilst the candidate-gene approach we used may ignore a number of potentially important SNPs, on the other hand it greatly reduces the chances of reporting false-positive findings and of describing a spurious result [21].

Fig. 2. Perfect Petri net model. a. Perfect Petri net for the Multifactor Dimensionality Reduction (MDR) model depicted in Fig. 1. The maximum capacity of each place is 10 tokens. From left to right the IL-6 C-174G , theIL-2 G-330T and the HLA-B*3501 gene units are shown. Respectively, the weights of gene-dependent arcs (*) are: wildtype = 4, 2, 2; heterozygous = 6, 5, 6; and homozygous = 8, 8, 6. When multiple transitions are enabled the one with the lowest index will fire first. See text for further details about Petri nets. b. Simplification of the Petri net model. One or more triggering events activate the network; the IL-6 gene is preferentially activated (bold line), followed by the IL-2 gene (solid line) or HLA-B*3501 (dotted line). Following activation a series of self-enhancing reactions will follow to eventually cause digital ulcers in relation to the individual genetic background.

352

L. Beretta et al. / European Journal of Internal Medicine 21 (2010) 347–353

as co-enhancing factors for IL-6-driven autocrine processes [42]. In our model, this co-enhancing effect is exerted by IL-2, a cytokine that may indeed promote IL-6 production in monocytes [43] or synergize with IL-6 to promote B lymphocyte activation [44], as well as autoantibody production [45], a marker of DU occurrence in other SSc populations [2,5]. Lastly, IL-6 would promote the expression of HLA-B*3501, a necessary co-factor in DUs susceptibility according to our model. Whilst the precise nature of IL-6 and HLA-B*3501 interaction is at the moment elusive and matter of further researches, the inducing effect of IL-6 on MHC class I expression has long been described [46]. Summarising, in the present paper, we provide evidence that the genetic background may predispose Italian SSc patients to DU occurrences. The causative factors we found are likely to interplay in a system where IL-6 exerts a key function and whose construct is biologically plausible and supported by several lines of evidence, but that constitutes, at the same time, a starting point for further translational research. 5. Learning points • Gene–gene interaction may be more relevant than single-genotype associations in determining clinical complications of systemic sclerosis. • The IL-6 C-174G, the IL-2 G-330T SNPs and the HLA-B*3501 allele epistatically interact to predict digital ulcers occurrence in Italian scleroderma patients. • System biology may provide biological model to explain statistical models of epistasis.

Acknowledgments The present work was supported by a grant from the Fondazione CARIPLO (CAssa di RIsparmio delle Province LOmbarde) and by a grant from GILS, “Gruppo Italiano per la Lotta alla Sclerodermia.” References [1] Ferri C, Valentini G, Cozzi F, Sebastiani M, Michelassi C, La Montagna G. Systemic Sclerosis Study Group of the Italian Society of Rheumatology (SIR-GSSSc). Systemic sclerosis: demographic, clinical, and serologic features and survival in 1,012 Italian patients. Medicine (Baltimore) 2002;81:139–53. [2] Tiev KP, Diot E, Clerson P, Dupuis-Siméon F, Hachulla E, Hatron PY, et al. Clinical features of scleroderma patients with or without prior or current ischemic digital ulcers: post-hoc analysis of a nationwide multicenter cohort (ItinerAIR-Sclerodermie). J Rheumatol 2009;36:1470–6. [3] Walker UA, Tyndall A, Czirják L, Denton C, Farge-Bancel D, Kowal-Bielecka O, et al. Clinical risk assessment of organ manifestations in systemic sclerosis: a report from the EULAR Scleroderma Trials And Research group database. Ann Rheum Dis 2007;66:754–63. [4] Hachulla E, Clerson P, Launay D, Lambert M, Morell-Dubois S, Queyrel V, et al. Natural history of ischemic digital ulcers in systemic sclerosis: single-center retrospective longitudinal study. J Rheumatol 2007;34:2423–30. [5] Sunderkötter C, Herrgott I, Brückner C, Moinzadeh P, Pfeiffer C, Gerss J, et al. Comparison of patients with and without digital ulcers in systemic sclerosis: detection of possible risk factors. Br J Dermatol 2009;160:835–43. [6] Caramaschi P, Martinelli N, Volpe A, Pieropan S, Tinazzi I, Patuzzo G, et al. A score of risk factors associated with ischemic digital ulcers in patients affected by systemic sclerosis treated with iloprost. Clin Rheumatol 2009;28:807–13. [7] Alivernini S, De Santis M, Tolusso B, Mannocci A, Bosello SL, Peluso G, et al. Skin ulcers in systemic sclerosis: determinants of presence and predictive factors of healing. J Am Acad Dermatol 2009;60:426–35. [8] Chung L, Fiorentino D. Digital ulcers in patients with systemic sclerosis. Autoimmun Rev 2008;5:125–8. [9] Allanore Y, Borderie D, Meune C, Lemaréchal H, Weber S, Ekindjian OG, et al. Increased plasma soluble CD40 ligand concentrations in systemic sclerosis and association with pulmonary arterial hypertension and digital ulcers. Ann Rheum Dis 2005;64:481–3. [10] Matucci-Cerinic M, Seibold JR. Digital ulcers and outcomes assessment in scleroderma. Rheumatology (Oxford) 2008;47:v46–7. [11] Heresi GA, Minai OA. Bosentan in systemic sclerosis. Drugs Today (Barc) 2008;44: 415–28.

[12] Suarez A, Castro P, Alonso R, Mozo L, Gutierrez C. Interindividual variations in constitutive interleukin-10 messenger RNA and protein levels and their association with genetic polymorphisms. Transplantation 2003;75:711–7. [13] Pociot F, Molvig J, Wogensen L, Worsaae H, Nerup J. A TaqI polymorphism in the human interleukin-1β (IL-1β) gene correlates with IL-1β secretion in vitro. Eur J Clin Invest 1995;22:396–402. [14] Fishman D, Faulds G, Jeffery R, Mohamed-Ali V, Yudkin JS, Humphries S, et al. The effect of novel polymorphisms in the interleukin-6 (IL-6) gene on IL-6 transcription and plasma IL-6 levels, and an association with systemic-onset juvenile chronic arthritis. J Clin Invest 1998;102:1369–76. [15] Hoffmann SC, Stanley EM, Darrin Cox E, Craighead N, DiMercurio BS, Koziol DE, et al. Association of cytokine polymorphic inheritance and in vitro cytokine production in anti-CD3/CD28-stimulated peripheral blood lymphocytes. Transplantation 2001;72:1444–50. [16] Hollegaard MV, Bidwell JL. Cytokine gene polymorphism in human disease: online databases. Genes Immun 2006;7:269–76. [17] Uboldi de Capei MU, Dametto E, Fasano ME, Rendine S, Curtoni ES. Genotyping for cytokine polymorphisms: allele frequencies in the Italian population. Eur J Immunogenet 2003;30:5–10. [18] Beretta L, Bertolotti F, Cappiello F, Barili M, Masciocchi M, Toussoun K, et al. Interleukin-1 gene complex polymorphisms in systemic sclerosis patients with severe restrictive lung physiology. Hum Immunol 2007;68:603–9. [19] Beretta L, Cappiello F, Barili M, Scorza R. Proximal interleukin-10 gene polymorphisms in Italian patients with systemic sclerosis. Tissue Antigens 2007;69:305–12. [20] Beretta L, Santaniello A, Cappiello F, Barili M, Scorza R. No evidence for a role of the proximal IL-6 G/C-174 single nucleotide polymorphism in Italian patients with systemic sclerosis. J Cell Mol Med 2007;11:896–8. [21] Beretta L, Cappiello F, Moore JH, Barili M, Greene CS, Scorza R. Ability of epistatic interactions of cytokine single-nucleotide polymorphisms to predict susceptibility to disease subsets in systemic sclerosis patients. Arthritis Rheum 2008;59:974–83. [22] Beretta L, Cappiello F, Moore JH, Scorza R. Interleukin-1 gene complex single nucleotide polymorphisms in systemic sclerosis: a further step ahead. Hum Immunol 2008;69:187–92. [23] Santaniello A, Salazar G, Lenna S, Antonioli R, Colombo G, Beretta L, et al. HLA-B35 upregulates the production of endothelin-1 in HLA-transfected cells: a possible pathogenetic role in pulmonary hypertension. Tissue Antigens 2006;68:239–44. [24] Salazar G, Colombo G, Lenna S, Antonioli R, Beretta L, Santaniello A, et al. HLA-B35 influences the apoptosis rate in human peripheral blood mononucleated cells and HLA-transfected cells. Hum Immunol 2007;68:181–91. [25] Grigolo B, Mazzetti I, Meliconi R, Bazzi S, Scorza R, Candela M. Anti-topoisomerase II alpha autoantibodies in systemic sclerosis-association with pulmonary hypertension and HLA-B35. Clin Exp Immunol 2000;121:539–43. [26] Anon. Subcommittee for Scleroderma Criteria of the American Rheumatism Association Diagnostic and Therapeutics Criteria Committee. Preliminary criteria for the classification of systemic sclerosis (scleroderma). Arthritis Rheum 1980;23:581–90. [27] LeRoy EC, Medsger Jr TA. Criteria for the classification of early systemic sclerosis. J Rheumatol 2001;28:1573–6. [28] White B, Bauer EA, Goldsmith LA, Hochberg MC, Katz LM, Korn JH. Guidelines for clinical trials in systemic sclerosis (scleroderma). I. Disease-modifying interventions. The American College of Rheumatlogy Committee on Design and Outcomes in Clinical Trials in Systemic Sclerosis. Sem Arthritis Rheum 1995;38:351–60. [29] Bayer PM, Bauerfeind S, Bienvenu J, Fabien N, Frei PC, Gilburd B, et al. Multicenter evaluation study on a new HEp2 ANA screening enzyme immune assay. J Autoimmun 1999;13:89–93. [30] Frezza D, Giambra V, Tolusso B, De Santis M, Bosello S, Vettori S, et al. Polymorphism of immunoglobulin enhancer element HS1, 2A: allele *2 associates with systemic sclerosis. Comparison with HLA-DR and DQ allele frequency. Ann Rheum Dis 2007;66:1210–5. [31] Pozzi S, Longo A, Ferrara GB. HLA-B locus sequence based typing. Tissue Antigens 1999;53:275–81. [32] Marsh SGE. HLA class II region sequences. Tissue Antigens 1998;51:467–507. [33] Ritchie MD, Hahn LW, Roodi N, Bailey LR, Dupont WD, Parl FF, et al. Multifactor dimensionality reduction reveals high-order interactions among estrogen metabolism genes in sporadic breast cancer. Am J Hum Genet 2001;69: 138–47. [34] Hahn LW, Ritchie MD, Moore JH. Multifactor dimensionality reduction software for detecting gene–gene and gene–enviroment interactions. Bioinformatics 2003;19:376–82. [35] Moore JH, White BC. Tuning ReliefF for genome-wide genetic analysis. Lect Notes Comput Sci 2007;4447:166–75. [36] Moore JH, Boczkob EM, Summar ML. Connecting the dots between genes, biochemistry, and disease susceptibility: systems biology modeling in human genetics. Mol Genet Metab 2005;84:104–11. [37] Mayo M. Learning Petri net models of non-linear gene interactions. ByoSystems 2005;82:74–82. [38] Mayo M, Beretta L. Evolving concurrent Petri net models of epistasis. LNAI 2010;5991:166–75. [39] Mouthon L, Mestre-Stanislas C, Bérezné A, Rannou F, Guilpain P, Revel M, Pagnoux C, Guillevin L, Fermanian J, Poiraudeau S. Impact of digital ulcers on disability and health-related quality of life in systemic sclerosis. Ann Rheum Dis 2010;69:214–7. [40] Nihtyanova SI, Brough GM, Black CM, Denton CP. Clinical burden of digital vasculopathy in limited and diffuse cutaneous systemic sclerosis. Ann Rheum Dis 2008;67:120–3.

L. Beretta et al. / European Journal of Internal Medicine 21 (2010) 347–353 [41] Kitani A, Hara M, Hirose T, Suzuki K, Kawakami M, Kawaguchi Y, et al. Autostimulatory effects of IL-6 on excessive B cell differentiation in patients with systemic lupus erythematosus: analysis of IL-6 production and IL-6R expression. Clin Exp Immunol 1992;88:75–83. [42] Ogura H, Murakami M, Okuyama Y, Tsuruoka M, Kitabayashi C, Kanamoto M, et al. Interleukin-17 promotes autoimmunity by triggering a positive-feedback loop via interleukin-6 induction. Immunity 2008;29:628–36. [43] Musso T, Espinoza-Delgado I, Pulkki K, Gusella GL, Longo DL, Varesio L. IL-2 induces IL-6 production in human monocytes. J Immunol Feb. 1992;148:795–800.

353

[44] Lipsky PE, Davis LS, Cush JJ, Oppenheimer-Marks N. The role of cytokines in the pathogenesis of rheumatoid arthritis. Springer Semin Immunopathol 1989;11: 123–62. [45] Splawski JB, McAnally LM, Lipsky PE. IL-2 dependence of the promotion of human B cell differentiation by IL-6 (BSF-2). J Immunol 1990;144:562–9. [46] Kuwana M, Medsger Jr TA, Wright TM. Analysis of soluble and cell surface factors regulating anti-DNA topoisomerase I autoantibody production demonstrates synergy between Th1 and Th2 autoreactive T cells. J Immunol 2000;164:6138–46.