Objective. Gene-gene interaction, or epistasis, is considered a ubiquitous component of complex human diseases such as systemic sclerosis (SSc). Epistasis is difficult to model by traditional parametric approaches; therefore, nonparametric computational algorithms, such as multifactor dimensionality reduction (MDR), have been developed. Methods. A total of 242 consecutive unrelated Italian SSc patients and an equal number of well-matched healthy controls were genotyped for 22 cytokine single-nucleotide polymorphisms (SNPs; 13 cytokine genes). The distribution of the SNPs between controls and SSc patients, controls and limited cutaneous SSc (lcSSc) patients, and controls and diffuse cutaneous SSc (dcSSc) patients was tested by the MDR constructive induction algorithm and by focused interaction testing framework (FITF), a logistic regression-based approach. Results. None of the studied SNPs had main independent effects on SSc or disease subset susceptibility, therefore no epistatic interaction was detectable by FITF. The MDR analysis showed a significant epistatic interaction among the interleukin-2 (IL-2) G-330T, IL-6 C-174G, and interferon-␥ AUTR5644T SNPs and the IL-1 receptor Cpst1970T, IL-6 Ant565G, and IL-10 C-819T SNPs in lcSSc and dcSSc susceptibility, respectively. The relevance of the single multilocus attributes constructed by the MDR inductive algorithm was then confirmed by the parametric approach (P < 0.001 for both controls versus lcSSc patients and controls versus dcSSc patients). Conclusion. We provide evidence for gene-gene interaction among cytokine SNPs in the context of SSc. The interaction among cytokine SNPs with a profibrotic or a regulatory function on profibrotic interleukins is relevant to the susceptibility to SSc subsets and it appears to be more important than the contribution of any single cytokine SNP.