HLA-B27 may modulate the interaction between ERAP1 polymorphisms and smoking in ankylosing spondylitis patients.

BACKGROUND: Ankylosing spondylitis (AS) is an autoimmune disease that affects the enthesis and synovial membrane of the spine, the sacroiliac vertebrae and peripheral joints. Genetic susceptibility to AS is mainly due to the presence of the HLA-B*27 (B27) allele, and endoplasmic reticulum aminopeptidase-1 (ERAP-1) plays a key role in antigen processing and presentation to HLA class I molecules. Tobacco consumption is one of the main environmental factors involved in the pathogenesis of various diseases, including AS. The objective of the present study was to evaluate the association and the interactive effects of variants of the ERAP1 gene with smoking in modulating the risk of AS. METHODS AND
RESULTS: A case-control study in the Mexican population. The association of two functional variants of the ERAP1 gene (rs30187 and rs27044) in patients with AS was analyzed by the allelic discrimination method using TaqMan probes. B27 was typified by PCR-SSP. The interaction between the variants of ERAP1 and B27 and smoking was assessed using the multifactorial dimensionality reduction (MDR) method. There was no significant association of the two variants of ERAP1 in the cases compared with the controls (P > 0.05); however, a strong interaction between the variants and smoking could be demonstrated, with entropy values of 4.97% for rs30187 and 5.13% for rs27044. In addition, these interaction effects were increased in patients carrying the B27 allele.
CONCLUSIONS: The rs30187 and rs27044 variants of the ERAP1 gene appear to potentiate the effect of smoking in patients with AS carrying the B27 allele.

Introduction

Ankylosing spondylitis (AS) is an inflammatory autoimmune disease that is part of the group of spondyloarthropathies (SpAs), which share clinical, genetic and radiographic characteristics. Seronegative SpAs are those that lack rheumatoid factor or anti-nuclear antibodies in serum, including AS, psoriatic arthritis, reactive arthritis and spondylitis associated with nonspecific inflammatory bowel diseases [1-3]. AS is characterized by inflammation of the peripheral and sacroiliac joints of the axial skeleton, and in some cases, it can cause enthesitis, pain in the area where a tendon, muscle or ligament inserts into bone. Later in the course of the disease, pain and stiffness appear in the lumbar back and hips. Within the natural history of the disease, fusion of vertebrae (ankylosis) is common, causing loss of flexibility and forward bending of the spine. Epidemiological data indicate that AS occurs in one per 200 individuals, and both the incidence and prevalence of AS in the general population depend on the frequency of occurrence of the B27 allele. It is estimated that between 0.5 and 1.0% of carriers of this marker in any population have the disease; in Mexico, it is estimated that the prevalence of AS is 0.09%. This disease affects men more than women at a 5:1 ratio, and the first symptoms appear between the second and third decades of life [4-6].

Although the etiology of AS has not been fully clarified, it has been proposed that it develops as a result of complex interactions between genetic components and environmental factors. Although 90% of patients with AS are carriers of the HLA-B*27 (B27) allele, only 1% to 5% of B27-positive people will develop AS, which indicates that genes other than B27 may play important roles in its etiology [7]. In this sense, recent genome-wide association studies (GWAS) have identified approximately 40 genes involved in AS, the most consistent of which is endoplasmic reticulum aminopeptidase 1 (ERAP1), which can contribute to up to 26% of the genetic susceptibility to AS. In fact, ERAP1 is considered the second most important gene after B27 in the pathogenesis of AS, although its seems to occur only in B27-positive patients [8-12].

ERAP-1 is a zinc (Zn)-dependent endoplasmic reticulum (ER) aminopeptidase that is not part of the main histocompatibility complex (MHC), but its biological function consists of antigenic processing and cutting long peptides to the lengths required for presentation by MHC Class I to T lymphocytes [13, 14]. The enzymatic activity of some cytosolic aminopeptidases, such as ERAP-1, can be modified due to the presence of heavy metals such as cadmium (Cd), which is found in various environmental sources such as fertilizers, batteries, some foods (mollusks and crustaceans) and tobacco smoke [15-18]. Additionally, the presence of single nucleotide polymorphisms (SNPs) in ERAP1 can potentially alter the enzymatic activity of its product by modifying the protein structure [19, 20].

Based on this information, we suggest that there may be a gene–environment interaction between ERAP1 and tobacco consumption. Therefore, the present study was designed to evaluate the association and interactive effects of functional SNPs of the ERAP1 gene and exposure to tobacco smoke as they relate to the modulation of the risk of AS.

Materials and methods

Study population

One hundred twenty-three individuals of Mexican descent older than 18 years of age who were geographically matched were included in this case–control study. Fifty-eight of them were men and women diagnosed with axial AS who were seen at the outpatient clinic of the Rheumatology Service of the Luis Guillermo Ibarra Ibarra National Institute of Rehabilitation (INR-LGII). The diagnosis of AS was based on the modified New York criteria [21]. Patients with other autoimmune diseases were not included in the study. Sixty-five men and women who declared that they were in good health, without symptoms suggestive of disease and without a family history of AS, were selected as healthy controls; they were patients’ companions, blood bank donors and personnel from the Human Resources and Human Communication Departments of the INR-LGII. This study was conducted under the criteria established in the Declaration of Helsinki and was derived from a protocol with registration number INR-51/19, which was approved by the Ethics and Research Committee of the INR-LGII.

Smoking parameters

All participants were given a validated questionnaire on tobacco consumption that collected the number of cigarettes smoked per day and their time of consumption.

DNA extraction

A venous blood sample obtained by venipuncture was collected into a tube with EDTA-K2 anticoagulant. Genomic DNA was isolated from 200 µl of blood using a commercial kit (QIAmp 250 DNA Blood Kit, Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The DNA obtained was quantified by spectrophotometry using a Nanodrop 2000 (Thermo Scientific), adjusted to 50 ng/μl with molecular-grade water and kept at − 80 °C until analysis.

B27 typing

B27 typing was performed by PCR-SSP in both study groups using the commercial kit HLA-Fluo-Gene B27 (Inno-Train, Diagnostik GmbH, Germany) and genomic DNA. PCR was performed with a Veriti thermal cycler (Applied BioSystems) following the manufacturer’s instructions. The fluorescence readings were interpreted using FluoVista equipment (Inno-Train, Germany) with Fluogene v1.5.5 software.

ERAP1 polymorphism genotyping

The nonsynonymous functional variants of the ERAP1 gene rs27044 (Gln730Glu) and rs30187 (Lys528Arg) were selected for genotyping in this study (cat. number 4351379 ThermoFisher Scientific, C___3056870_10 and C___3056885_10, respectively). The selection criteria were based on the fact that both variants had been previously explored in other populations; additionally, these two variants have been previously reported by our study group [12, 13, 22, 23] and that the minor allele frequency (MAF) in the Mexican population was greater than 1% (http://www.ensembl.org). Genotyping was performed using the allelic discrimination method with TaqMan probes (Applied Biosystems, Foster City, CA, USA). For PCR, a mixture was prepared with 5 μl of TaqMan Universal PCR Master Mix (Applied BioSystems, Warrington, United Kingdom), 0.25 μl of TaqMan probe at 20×, 0.25 μl of water and 4.5 μl of genomic DNA. All samples were amplified using StepOne Plus real-time PCR equipment (Applied Biosystems) following the manufacturer’s protocol. Allelic discrimination analysis was performed with StepOne v2.3 software (Applied BioSystems).

Statistical analysis

Data were analyzed using software SPSS 21.0 for Windows® (SPSS Inc., Chicago, IL, USA). Continuous variables are expressed as the mean ± standard deviation (SD), while categorical variables are expressed as frequencies and percentages. The χ2 test and Student’s t‑test were used to compare clinical characteristics between the two groups, a P-value < 0.05 was considered statistically significant. The χ2 test was also used to determine the Hardy‑Weinberg equilibrium (HWE) in the control group. Differences between the non-continuous variables and genotype distribution frequency were assessed using the χ2 test. The odds ratio (OR) and 95% confidence interval (CI) were evaluated using binary logistic regression analysis. The Bonferroni’s test was used to determine the significance level to correct multiple test errors, in which taking into account the two selected SNPs, an adjusted P-value < 0.025 (α/number of loci) was considered statistically significant. The statistical power was calculated using the Piface application freely available online (Java Applets for power and sample size, http://www.stat.uiowa.edu/~rlenth/Power, accessed on 30 march 2022).

The multifactor dimensionality reduction (MDR) method was applied to assess interactive effects of ERAP1 variants and tobacco exposure in modulating AS risk. The analysis was conducted by MDR software 3.0.2 (Computational Genetics Laboratory, University of Pennsylvania, USA; available for free at https://sourceforge.net/projects/mdr/). MDR is a data mining strategy for detecting and characterizing nonlinear interactions among discrete attributes (e.g. SNPs, smoking, gender) that are predictive of a discrete outcome (e.g. case–control status, outcome, etc.). It combines attribute selection, attribute construction and classification with cross-validation to provide a powerful approach to modeling interactions. The best n-factor interaction model in predicting AS risk was identified with the maximum cross-validation consistency (CVC) and the optimal testing accuracy. An interaction map was prepared to show the interaction of individual factors in the best predictive model through information gain values (entropy percentage).

Results

Characteristics of the study population

At the time of the present study, 58 patients with AS and 65 healthy controls were analyzed with rigorous inclusion criteria since 2019 year. The characteristics of the study population are shown in Table 1. The cases were younger than the controls, 46.1 ± 13.5 years vs. 53.9 ± 5.7 years, respectively (P < 0.001); 38.0% of the patients were women, and 62.0% were men. Regarding B27 typing, 70.7% of patients were positive and 29.3% were negative, while 98.5% of the controls were negative and 1.5% were positive (P < 0.001). The average age at the time of AS diagnosis was 36.1 ± 13.9 years. The values on the Bath Ankylosing Spondylitis Functional Index (BASFI) and the Bath Ankylosing Spondylitis Disease Activity Index (BASDAI) scales were 4.9 and 3.73, respectively. There were no significant differences in tobacco consumption between the two study groups (P = 0.228), and only there were 11 smoking subjects (patients and controls) who are B27 positive. Considering to the B27 as an exposure factor, the statistical power reaches 100%; however; for cases and controls who are B27 positive, but who are exposed to smoking, the statistical power was 76.2%.

Table 1

Demographic and clinical characteristics of AS patients and healthy controls

	AS patients (New York criteria)(n = 58)	Controls (n = 65)	P-value
Age (years)	46.1 ± 13.5	53.9 ± 5.7	< 0.001*
Gender
Female (%)	22 (38.0)	60 (92.3)	< 0.001**
Male (%)	36 (62.0)	5 (7.7)
HLA-B27 (%)
Negative	17 (29.3)	64 (98.5)	< 0.001***
Positive	41 (70.7)	1 (1.5)
Age at diagnosis (years)	36.1 ± 13.9	NA
BASFI	4.9	NA
BASDAI	3.73	NA
Smoking (%)
No	42 (72.4)	53 (81.5)	0.228**
Yes	16 (27.6)	12 (18.5)

The variables are expressed as the mean ± standard deviation (SD). *P-values were estimated using t-test, α = 0.05; **P-values were estimated using χ2 test; ***P-values were estimated using Fisher´s exact test, α = 0.05; significant P-values are in bold

AS ankylosing spondylitis; HLA human leukocyte antigen; BASFI bath ankylosing spondylitis functional index; BASDAI bath ankylosing spondylitis disease activity index; NA not applicable

Distribution of the gene and allelic frequencies of the studied SNPs

The distribution of the rs30187 and 27044 genotypes of the ERAP1 gene in the control group was consistent with the HWE (P > 0.05), and the MAF values in the Mexican population were (T) = 0.39 and (G) = 0.34, respectively. Table 2 shows the distribution of the genotypes of the patients with AS and healthy controls. For the two SNPs of ERAP1 studied, no statistically significant associations with AS were demonstrated under any of the three inheritance models that were analyzed (P > 0.025).

Table 2

Analyses of the association of two SNPs with ankylosing spondylitis

SNP	AS cases N (%)	Control N (%)	OR*	95% CI	P-value
ERAP1 rs30187
C/C	20 (34.4)	23 (35.4)	1.00	(Reference)
C/T	28 (48.3)	32 (49.2)	2.09	(0.13–1.73)	0.262
T/T	10 (17.2)	10 (15.4)	1.60	(0.18–2.09)	0.444
Dominant model
C/C	20 (34.4)	23 (35.4)	1.00	(Reference)
C/T + T/T	48 (65.5)	42 (64.6)	1.31	(0.63–2.72)	0.461
Recessive model
C/C + C/T	30 (82.7)	55 (84.6)	1.00	(Reference)
T/T	10 (17.2)	10 (15.4)	1.83	(0.68–4.89)	0.223
HWE		0.836
ERAP1 rs27044
C/C	23 (39.6)	27 (41.5)	1.00	(Reference)
C/G	26 (44.8)	29 (44.6)	1.96	(0.14–1.83)	0.302
G/G	9 (15.5)	9 (13.8)	2.12	(0.13–1.70)	0.254
Dominant model
C/C	23 (39.6)	27 (41.5)	1.00	(Reference)
C/G + G/G	34 (60.3)	38 (58.4)	1.05	(0.50–2.16)	0.864
Recessive model
C/C + C/G	49 (84.4)	56 (86.1)	1.00	(Reference)
G/G	9 (15.5)	9 (13.8)	1.14	(0.42–3.10)	0.793
HWE		0.786

AS ankylosing spondylitis; ERAP1 endoplasmic reticulum aminopeptidase 1; OR odds ratio; CI confidence interval; SNP single nucleotide polymorphism; HWE Hardy–Weinberg equilibrium

*Adjusted for age and gender

Evaluation of gene–gene and gene–environment interactions: MDR results

An exhaustive MDR analysis revealed that the best interaction model for predicting the development of AS was the interaction of B27, smoking and the two polymorphisms of ERAP1 (rs30187 and rs27044). This model had a maximum balanced test precision of 0.8513, a CVC of 10/10 and a significant P-value of 0.0016 (Table 3). Figure 1 shows the interaction map of all the factors based on the measurements of entropy between the individual variables. For this case, a strong interaction effect was observed between smoking and the two variants, rs27044 and rs30187, with information gains of 5.13% and 4.97%, respectively.

Table 3

Results of MDR analysis

	Number of the risk factors	Testing balanced accuracy	CVC	P-value*
1	B27	0.8381	9/10	0.0101
2	B27, smoke	0.8458	9/10	0.0107
3	B27, rs30187, rs27044	0.8509	10/10	0.0130
4	B27, smoke, rs30187, rs27044	0.8513a	10/10	0.0016
5	Smoke, rs30187, rs27044	0.3983	10/10	0.0210
6	Smoke, rs30187	0.4447	8/10	0.6773
7	Smoke	0.5456	7/10	0.7032

The model with the maximum testing balnced accuracy and maximum CVC was considered as the best model

MDR multifactor dimensionality reduction, CVC cross validation consistency

*P-values were based on 1000 permutations

aThe best interaction model in MDR analysis

Fig. 1

Interaction map for ankylosing spondylitis risk. The interaction model describes the percentage of the entropy (information gain) that is explained by each factor or 2-way interaction. Values inside nodes indicate information gain of individual attributes or main effects, whereas values between nodes show information gain of pairwise combinations of attributes or interaction effects. Positive entropy (plotted in red or orange) indicates interaction, which can be interpreted as a synergistic or nonadditive relationship; while negative entropy (plotted in yellow-green or green) indicates independence or additivity (redundancy). (Color figure online)

Through a conditioned analysis, we showed that the presence of B27 increases the interaction between smoking and the two variants, rs30187 and rs27044. In the absence of B27, the entropy values between smoking and these two polymorphisms were − 0.18% and 1.11%, respectively; while the individual effects of each factor were 0.81% and 0.37%, respectively (Fig. 2a). However, the presence of B27 increases both the percentage of entropy of the two polymorphisms (to 2.24% each) and the individual effects (to 6.37% each) (Fig. 2b).

Fig. 2

Conditioned analysis for B27 negative (−) or B27 positive (+). The presence of B27 increases the interactive effects between smoking and ERAP1 polymorphisms

Discussion

Chronic exposure to tobacco smoke through cigarette smoking is a serious public health problem that can lead to the development of serious diseases such as lung cancer, cardiovascular disease and chronic obstructive pulmonary disease (COPD) [24, 25]. At the musculoskeletal level, there is consistent evidence that smoking is an important risk factor in the development of rheumatoid arthritis (RA), mainly in men, through its interaction with the citrullination process or with anti-citrullinated peptide antibodies (ACPAs) [26, 27]. In the case of osteoporosis, smoking affects bone calcium homeostasis, and in osteoarthritis (OA), there is inflammation with a notable increase in oxidative stress [17]. With respect to AS, clinical studies indicate that tobacco use is associated with erectile dysfunction, cardiovascular risk, significant changes in acute phase reactants, systemic inflammation and decreased physical activity [28-30]. In addition, patients with AS who smoke experience increased pain and fatigue, leading to a poor quality of life [31].

It has been shown that the interaction between genetic and environmental factors, such as tobacco smoke, is crucial for the development of RA and OA [32, 33]. In the present study, we analyzed both the association of polymorphisms of the ERAP1 gene with AS and the effects of their interaction with smoking on the risk of developing AS. Our results showed that, despite the great representativeness of the MAF of the two variants of ERAP1 that were analyzed in our population, there was no significant association with the genetic susceptibility to AS. Several studies have been conducted to examine the association between SNPs of the ERAP1 gene and AS, and the results are inconsistent. Some argue that the association with ERAP1 occurs only in B27-positive patients, which suggests that there could be an epistatic gene–gene interaction [11, 22, 34]. Recently, it was discovered that ERAP1 variants interact genetically with the B27 and B*4001 alleles in patients with AS, indicating that ERAP1 risk variants can be present in both B27-positive and B27-negative/B*4001-positive patients [35]. However, other reports show that there is no association with ERAP1 in either B27-positive or B27-negative patients [36]. A study conducted in patients of Australian origin with AS, whether B27 positive or B27 negative, in the presence of risk or protective ERAP1 variants did not show significant differences in the expression of ER stress markers or proinflammatory cytokines [37]. Apparently, the association with B27 in patients with AS occurs mainly in populations of European Caucasian origin and is present at a very low or almost null frequency in populations of Asian origin [38]. It is possible that epistatic mechanisms are closely related to the genetic structure of populations in such a way that the interaction between ERAP1 and B27 may differ among them and that they may therefore manifest very particular phenotypes.

The ERAP1 gene is located on chromosome 5q15, and its product, the ERAP-1 protein, is a metallopeptidase whose enzymatic activity is strongly affected by the presence of the rs30187 and rs27044 variants; the first appears to decrease it, and the second increases it, but its effects on the ability to process antigenic peptides are not yet conclusive [39, 40]. Another important factor that alters the enzymatic activity of ERAP-1 is the presence of heavy metals, such as Cd, in tobacco smoke. A possible mechanism for the involvement of Cd in ERAP-1 enzymatic activity is that because it is an enzyme that depends on Zn2+ and Cd2+ having the same valence, its biological actions are described as ionic mimicry of Zn2+, and it may even compete with Zn2+ for the same target site [41] since they use the same transport mechanisms for entry and for the regulation of cellular homeostasis [42].

We used MDR to evaluate the interaction between the rs30187 and rs27044 variants of ERAP1 and smoking in patients with AS. Our results revealed interactions with a high degree of synergy between smoking and the ERAP1 variants, as well as a discrete epistatic interaction between rs30187 and rs27044. Age and sex play important roles in the development of AS, and through the MDR method, we were able to corroborate this finding due to the strong interaction observed between these two variables (Fig. 1). Overall, the best interaction model we obtained was for the interaction between B27, smoking, rs30187 and rs27044 (Table 3). Subsequently, we performed a conditioned analysis to determine the effect of B27 on the interaction of smoking with the two variants of ERAP1. Interestingly, the presence of B27 increases the entropy (or interaction) between smoking and the two variants (Fig. 2). Each of these variables exerts an individual effect, and as they are added, the risk of developing AS increases. B27 is the most important factor in genetic susceptibility to AS, but we confirmed that other factors, such as the ERAP1 gene and smoking, act synergistically, enhancing the risk.

Several studies of SNP associations under the SNP × SNP pair scheme have been proposed to evaluate this interaction. We decided to perform an analysis of the interaction between all of the variables to highlight each of the individual effects; however, we are aware of the limitations of our work. First, the sample size was small because during the recruitment period, there was limited access of the population to medical services at our institution due to the COVID-19 pandemic, which affects the statistical power, so the results must be taken with caution. Second, our study focused on a single population sampled from a single hospital, and therefore, it would be worthwhile to replicate the study with a multicenter sample to corroborate our findings. Third, there are numerous SNPs within the ERAP1 gene and other genes (such as those in the IL-23 pathway, KIRs and TLR4, among others) that were not considered but impact the modulation of AS. Finally, it is possible that, in addition to Cd, the lead, arsenic and chromium present in tobacco smoke could favor the development of AS; however, these elements were not considered.

Conclusions

In summary, the rs30187 and rs27044 variants of ERAP1 were not associated with genetic susceptibility to AS in our population; however, we suggest that they play an important role in enhancing the effects of smoking in patients carrying B27. More studies are needed to explore in depth the mechanisms of interaction between these polymorphisms, as well as those of other genes and other environmental factors, in the development of AS.