The role of genetic polymorphisms of interleukin-1 (IL-1R1 and IL-1RN) in primary knee osteoarthritis in Indonesia

This study aimed to evaluate the association of SNPs of the IL-1 family with the clinical severity of knee OA. This case‒control study was performed among 100 healthy knees and 130 osteoarthritis (OA) knees of people aged ≥ 50 years with a BMI ≥ 25 kg/m2. The possible correlations among clinical findings, radiographic evaluations, serum levels of IL-1R1 and IL-1Ra, and genotype analyses were evaluated. Three SNPs of IL-1R1, rs871659, rs3771202, and rs3917238, were associated with primary knee OA. Females with IL-1R1 SNP rs871659 allele A had a higher prevalence of primary knee OA. No correlation was found between SNPs of IL-1R1 and IL-1RN and clinical or radiologic severity or serum concentrations of IL-1R1 and IL-1Ra (p > 0.05). BMI and IL-1R1 rs3917238 genotype C/C were correlated with moderate-severe VAS scores. A correlation was also found between the EQ-5D-3L self-care dimension and obesity and between the EQ-5D-3L pain and usual activity dimensions and age ≥ 60 and obesity (p < 0.05). Radiologic severity was only associated with age ≥ 60 years (p < 0.05). We found the IL-1R1 SNPs rs871659, rs3771202, and rs3917238 to be predisposing factors for primary knee osteoarthritis. The clinical findings, radiographic severity, and serum concentrations of IL-1R1 and IL-1Ra were not correlated with these gene polymorphisms.


Materials and methods
This was a case-control study of primary knee osteoarthritis patients in Mongoloid who sought medical treatment at our orthopaedic outpatient clinic from October 2019 to August 2021. The control group was selected from patients' relative and friends who has no knee complaint who fulfilled the inclusion and exclusion criteria. The inclusion criteria for controls were healthy knees, age over 50 years with a BMI over 25 kg/m 2 , no knee pain (VAS score = 0), no functional knee limitations (Knee Injury and Osteoarthritis Outcome Score (KOOS) = 100), and normal knee radiology. For the case group, the inclusion criteria were patient age over 50 years with a BMI over 25 kg/m 2 who were diagnosed with primary knee osteoarthritis radiologically and had knee pain or functional limitations. All subjects in both groups agreed to follow all research procedures, including structured interview sessions and complete physical examinations by the principal investigator or appointed research assistant in this study. The exclusion criteria for the case and control group were patients with trauma or fracture around the knee, previous traumatic knee injury or any history of trauma, arthritis, ankylosing spondylitis, septic arthritis, other arthritis or any other systemic inflammatory or autoimmune disorder, and any cancer. Based on a sample size calculation with a power of 80%, β = 0.2 and α = 0.05 using binomial proportions, the sample for the control group was 100 patients, and the sample for the case group was 130 patients, who were divided into grade 2 with mild OA (65 patients) and grades 3-4 with moderate-severe OA (65 patients). We evaluated the clinical findings (using the VAS, EQ-5D-3L, and KOOS), radiological scores (Kellgren-Lawrence), laboratory findings (protein IL-1R1 and IL-1Ra) in plasma, and genetic polymorphisms of IL-1R1 and IL-1RN with 4 SNPs (IL-1R1_rs228704, IL-1RN_rs419598, IL-1RN_rs9005, and IL-1RN_rs315943) based on previous studies performed by Attur et al. 10 , Wu et al. 7 , and Ahmed et al. 8 . To screen other SNPs that might correlate with primary knee OA, we performed next-generation sequencing (NGS). This research was approved by the Ethics Review Committee of the Faculty of Medicine, Public Health and Nursing, Universitas Gadjah Mada.
Blood samples. Blood samples were collected after informed consent was obtained from patients with normal knees and those with knee osteoarthritis diagnosed by radiograph examination. Patients and families were informed that samples would be taken to establish solid diagnoses and for research purposes. For DNA isolation (genotype analysis), we used EDTA tubes to keep whole blood, while for ELISA (phenotyping analysis), we used blood serum to exclude any inflammatory disease and cancer.

Phenotyping analysis.
To evaluate serum IL-1Ra, human IL-1Ra standard was reconstituted with distilled water to make a standard stock solution with a concentration of 20,000 pg/mL. The standard stock solution was then diluted at a 1:1 ratio with Calibrator Diluent RD5-33 (2X dilution per level) to achieve the following concentrations (Calibrator Diluent RD5-33 was also used as standard 0 pg/mL). A total of 100 µL of Assay Diluent RD15 was added to each well. Standards and samples were put into the wells according to the following procedure. Plates were incubated using Lovibond TC 135 S (Tintometer GmbH, Dortmund, Germany) for 2 h at room temperature (± 20 °C). The plates were washed with wash buffer 4 times using an automated washer, Biochrom Anthos Fluido 2 Microplate Washer (Biochrom Ltd., Cambridge, UK). A total of 200 µL of human IL-1Ra conjugate was added to each well. Plates were incubated for 2 h at room temperature (± 20 °C). The plates were washed with wash buffer 4 times using an automated washer. Substrate solution was prepared by mixing Color Reagent A and Color Reagent B at a ratio of 1:1. A total of 200 µL of substrate solution was added to each well. Plates were incubated for 30 min at room temperature. After the addition of substrate and incubation for 30 min, the liquid in the well appeared blue. A total of 50 µL of Stop Solution was added to each well. After the addition of Stop Solution, the liquid in the well changed colour from blue to yellow. Absorbance was measured with a microplate reader Bio-Rad model 680 (Bio-Rad Laboratories Inc., CA, USA) and Microplate Manager ver. 5.2.1 software (Bio-Rad Laboratories Inc., CA, USA) at a wavelength of 450 nm and wavelength correction at 540 nm. Standard curves were made based on optical density (OD) and standard concentrations of IL-1Ra. The concentration in the sample was determined based on the standard curve equation.
To evaluate serum IL-1R1, the capture antibody stock solution was diluted with PBS. A total of 100 µL of the diluted capture antibody was added to each well. The plate was covered with a plate sealer and then incubated overnight at room temperature (± 20 °C). The solution in the well was discarded, and then the plate was washed with wash buffer 3 times. A total of 300 µL of reagent diluent was added to each well (blocking stage). Plates were incubated for a minimum of 1 h at room temperature. The plates were washed with wash buffer 3 times. The plates were then ready to be used for adding samples. Human IL-1R1 standard was reconstituted with distilled water to make a standard stock solution with a concentration of 8000 pg/mL. The standard stock solution was then diluted at a 1:1 ratio (2X dilution) or 1:2 ratio (3X dilution) with Reagent Diluent to achieve the following concentrations (Reagent Diluent was also used as standard 0 pg/mL). Standards and samples were put into the wells in the following procedure. A total of 100 µL of standard 1 to standard 9 was added to wells A1 to D3. A total of 100 µL of the sample was added to wells E3 to H12. Plates were incubated for 2 h at room temperature (± 20 °C). The plates were washed with wash buffer 3 times using an automated washer. A total of 100 µL of detection antibody was added to each well. Plates were incubated for 2 h at room temperature. The plates were washed with wash buffer 3 times using an automated washer. A total of 100 µL of streptavidin-HRP was added to each well. Plates were incubated for 20 min at room temperature. The plates were washed with wash buffer 3 times using an automated washer. Substrate solution was prepared by mixing Color Reagent A and Color Reagent B at a ratio of 1:1. A total of 100 µL of substrate solution was added to each well. Plates were incubated for 20 min Genotyping analysis. All genotyping was performed in a clinical genetics laboratory (Eijkman Institute for Molecular Biology). For genotyping analysis, we first designed primers for the targeted sequences of the IL-1R1 and IL-1RN genes. DNA was extracted from blood samples using standard protocols. Before polymerase chain reaction (PCR), DNA was diluted to adjust concentrations to within a range compatible with multiplex PCR conditions. Then, we performed sequence amplification using PCR and next-generation sequencing (NGS), followed by analysing the data obtained. We used NGS to screen other polymorphisms that might correlate with primary knee OA. This protocol used the Nextera XT DNA Library Prep Kit (Illumina) according to the manufacturer's instructions. Purified DNA from the same sample was pooled together. Pooled DNA was fragmented using an enzyme or sonicator and tagged with indices. All samples were amplified again using a thermal cycler to increase the amount of DNA in the library. Excess reagent was cleaned from the amplified library. The randomness of library amplification could lead to different amounts of amplified library for each sample. Therefore, it was necessary to normalize the libraries to even the amounts of each amplicon. All cleaned DNA libraries were pooled together for cluster generation. The index in each DNA sample acted as a unique marker to differentiate each sample for data analysis. We used GRCh38(hg38) as a reference for genotype analysis. The sequence of each DNA strand was read at the same time when new DNA strands were being synthesized in the flow cell. Cluster generation and sequencing were performed in a sequencer (Illumina MiSeq Next Generation Sequencer). DNA sequencing was analysed with bioinformatics. The primer sequences for this research are described in Table 1.

Data analysis.
Genotyping performance was assessed as the total genotyping success rate, genotype concordance rate between duplicates, and Hardy-Weinberg equilibrium (HWE). SNPs with minor allele frequencies < 0.01 were excluded. We analysed the data using crosstabulation with chi-square tests for nominal data. The independent t test or Mann-Whitney U test was used for ratio/interval data, depending on the data distribution. We also evaluated the correlation of genetic polymorphisms with clinical findings using coefficient contingency for nominal data. Logistic regression analysis was used to examine the association between SNP variants and sex, age, BMI, clinical findings, radiologic severity, and laboratory findings. We used the smallest frequencies of alleles and genotypes as a reference, as we did not have any genotype mapping information for the Indonesian population. All statistical analyses were performed using IBM Statistical Package for the Social Sciences version 24.0 (SPSS Inc., Chicago, IL, USA).
Ethical approval. The

Results
The total number of patients in the control group was 100, with 72 (72%) females and 28 (28%) males and an average age of 57.35 ± 6.34 years. In the case group, there were 130 patients, divided into 65 patients in the mild subgroup and 65 patients in the moderate-severe subgroup. In the mild subgroup, 48 (73.85%) patients were female, and 17 (26.15%) were male, with an average age of 58.91 ± 6.85 years. The moderate-severe subgroup consisted of 51 females (78.46%) and 14 males (21.54%), with an average age of 62.00 ± 7.27 years, as shown in Table 2.
We did not find any difference in the frequency of 4 SNPs (IL-1R1_rs228704, IL-1RN_rs419598, IL-1RN_ rs9005, and IL-1RN_rs315943) that we targeted. However, we found another SNP that showed differences between groups (IL-1R1_rs871659 and IL-1R1_rs3771202) in terms of genotype (Table 3) or (IL-1R1_rs871659, IL-1R1_rs3771202, and IL-1R1_rs3917238) allele (Table 4). Therefore, the SNPs that we analysed further were IL-1R1_rs871659, IL-1R1_rs3771202, and IL-1R1_rs3917238. The gene variant rs871659 is located on chromosome 2:102,155,395 (GRCh38.p13), rs3771202 is located on chromosome 2:102,156,209 (GRCh38.p13) and rs3917238 is located on chromosome 2:102,156,623 (GRCh38.p13). All three gene variants are located in the IL-1R1 gene and are intron variants. In genotype analysis, we found that rs871659 and rs3771202 were associated with the incidence of mild OA compared to controls, and we found that rs871659 A/A and rs3771202 G/G had  -GCA AGG TGT TCC GTC TCT T-3'  5'-GCC AAG GAC ATA GTC AAC AAT-3'   2  IL-1RN 2  5'-TGT ATG AAA GAT GGC TGT GC-3'  5'-CAA GTC CCC AGT CTT CTC AAA-3'   3  IL-1R1  5'-ATG TGT TCT TCC TTC CCC AGG-3' 5'-ACA GAG ACG AGT ATG AGA AAT GAC-3 www.nature.com/scientificreports/ the highest prevalence among mild OA cases. In allele analysis, we found that rs871659 allele A, rs3771202 allele G, and rs3917238 allele T were associated with the incidence of mild OA compared to controls. We found a higher prevalence of the IL-1R1 rs3917238 genotypes C/T (OR = 1.7, 95% CI: 1.4-2.1) and T/T (OR = 1.9, 95% CI: 1.5-2.5) among those 3 SNPs. A higher prevalence was also found among people aged ≥ 60 years (OR = 2.1, 95% CI: 1.8-2.4) and those with obesity (OR = 3.6, 95% CI: 1.9-6.8). When we performed multivariate analysis by genotype, we found that only older age (≥ 60 years) with an OR of 2.4 (95% CI: 1.3-4.3) and obesity with an OR of 4.4 (95% CI: 1.4-14.1) were associated with primary knee OA, but not for the SNP genotype (Table 5). However, the result was different when stratified by sex. After stratification by sex, those variables remained associated with knee OA only among females, among whom those who were older, had a higher BMI and had SNP rs871659 allele A showed a higher risk of knee OA. We also found that rs3771202 allele G was associated with a higher risk of knee OA among males (Table 6).
We analysed the association between sex, age, BMI, and SNPs of IL-1R1 and IL-1RN and clinical findings. We found an association between moderate-severe VAS scores and obesity with an OR of 2.8 (95% CI: 1.0-7.5). We also found a higher prevalence of moderate-severe VAS scores in IL-1R1_rs3917238 (C/C) (Supplementary Table 1).
We analysed the association between sex, age, BMI, and genetic polymorphisms of IL-1R1 and IL-1RN with EQ-5D-3L. None of the variables were associated with problems in the EQ-5D-3L mobility or anxiety dimension. Regarding the self-care dimension in the EQ-5D-3L, we found that older age was associated with self-care problems, with an OR of 2.3 (95% CI: 1.1-4.9). None of the SNPs were associated with self-care problems (Supplementary Table 2). When we analysed the association of the EQ-5D-3L usual activity dimension with sex, 3) were associated with problems in the EQ-5D-3L pain dimension, while none of the SNPs were associated with pain problems (Supplementary Table 4). We found an association between radiologic severity and age. Severe OA was associated with older age (≥ 60 years), with an OR of 2.7 (95% CI: 2.1-3.5). We also found a higher prevalence of severe OA in  www.nature.com/scientificreports/ IL-1R1_rs871659 (G/A), IL-1R1_rs3771202 (C/G) and IL-1R1_rs3917238 (C/C); however, after adjustment in multivariate analysis, only older age was associated with severe OA (Supplementary Table 5). We analysed the association between genetic polymorphisms of IL-1R1 and IL-1RN and IL-1R1 and IL-1Ra protein concentrations in plasma in primary knee OA cases in Indonesia. After stratification by the covariate, we found an association between the gene variant rs871659 genotype A/A, rs3771202 genotype G/G and rs3917238 genotype T/T and the plasma concentration of IL-1R1 protein (Supplementary Table 6). We also found an association between the gene variant rs3917238 genotype T/T and the plasma concentration of IL-1Ra protein (Supplementary Table 7).
We evaluated the association between the haplotypes of the three SNPs and the OA group in this study. Each existing SNP had 2 allele variants, A and G for rs871659, G and C for rs2771202, and T and C for rs3917238. We used a combination of alleles from these three SNPs to perform haplotype analysis. From all samples analysed, only 3 haplotypes from these three SNPs were found, namely, AGT, AGC, and GCC. Major haplotypes were used as a reference to calculate ORs and p values. Based on calculations, GCC haplotypes were associated with OA between the control group and the mild group (Supplementary Table 8) but not between the control and severe groups (Supplementary Table 9). The GCC haplotype was also associated with knee OA (Supplementary Table 10). Regarding HWE, we found that in the control versus mild group, rs871659 and rs3771202 did not follow the HWE (p value < 0.05) (Supplementary Table 11).

Discussion
Osteoarthritis is a degenerative disease that is influenced by multiple factors, such as age, BMI, sex, mechanical abnormalities, previous injury, and genetics 11 . Recently, more studies have focused on the relationship between genetic polymorphisms and osteoarthritis. This is the first case-control study to be performed in Indonesia that focuses on the relationship between genetic polymorphisms and knee osteoarthritis and its correlation to clinical findings.
We evaluated the incidence of primary OA, and we compared the healthy knee (control) with the knee OA group (mild, moderate-severe groups) and matched them based on age, sex, and BMI. We still found differences in age and BMI, where patients aged ≥ 60 years had a higher incidence of knee OA than those aged < 60 years. The incidence of primary knee OA was also higher according to BMI, specifically for overweight and obese patients. We also found that age ≥ 60 years was correlated with primary knee OA. Obesity also had a correlation with primary knee OA. These results were also seen in previous studies that showed age and obesity to be associated with knee OA [12][13][14][15][16] . However, we could not find any association between sex and primary knee OA. This finding is different from that of a previous study 12,17 that found female sex to be a factor associated with knee OA. www.nature.com/scientificreports/ In our study, we performed genotype analysis on IL-1RN based on multiple previous studies in East Asia and showed that several polymorphisms were related to primary knee OA and its severity. A study by Attur in 2009 showed the role of IL-1RN variants in radiographic severity. The study revealed that carriers of the CTA haplotype (rs419598/rs315952/rs9005) had a decreased risk for severe knee OA. Another polymorphism at rs419598 was also reported to be associated with radiographic severity 10 . Another study by Wu and colleagues in 2012 showed that an IL-1RN gene variant (rs419598/rs9005/rs315943) was associated with increased progression of knee OA 7 . A study by Ahmed in 2016 showed that genetic variants in IL-1RN (rs9005 and rs315943) were associated with the risk of OA development 8 .
Another study conducted in different geographical areas showed genetic variants of IL-1A (rs1800587), IL-1B (rs1143634), IL-1B (rs16944), and IL-1RN (VNTR). Kaarvatn and his colleagues found a clear association that only affected women. These genetic variants were associated with a sixfold lower risk of knee OA. This finding implicated a novel sex-specific issue to consider regarding OA susceptibility 18 . These findings showed that geographical differences might influence the development of certain diseases and that the precision of ethnicity in sample selection should prevent confounding variables that might interfere with the results of the study. www.nature.com/scientificreports/ However, our study showed no association between IL-1RN gene variants rs419598, rs9005, and rs315943 and IL-1R1 rs2287047. Conflicting results were also found in previous systematic reviews and meta-analyses regarding the effect of genetic polymorphisms on the incidence of primary knee OA and clinical, radiologic, and laboratory findings 19 . Studies on the relationship of polymorphisms of the IL-1 family with the severity of knee OA are still scarce.
To provide better gene variant mapping in Indonesian knee OA cases, we expanded primer selection to all variants located in IL-1RN and IL-1R1. This action resulted in different polymorphisms that were related to primary knee OA cases in Indonesia. We found 3 SNPs located in the IL-1R1 gene to be correlated with primary knee OA. These three variants (rs871659, rs3771202, and rs3917238) were found to be associated with predisposing factors for primary knee OA. We found a higher prevalence of knee OA in rs3917238 C/T and T/T. However, when we performed a multivariate analysis with age, sex, and BMI, the impact was diminished. This is probably because the impact of this SNP is less than the impact of age and BMI in primary knee OA. The other reason is probably that this SNP was detected in small amounts; therefore, we need a larger sample size to evaluate whether this SNP impacts the incidence of primary knee OA. However, together with age and BMI, this SNP affected the incidence of primary knee OA. We also found a higher prevalence of severe VAS scores in rs3917238 genotype C/C. We did not find any SNP that correlated with the EQ-5D-3L mobility, EQ-5D-3L self-care, EQ-5D-3L usual activity, EQ-5D-3L pain, or EQ-5D-3L anxiety. We also did not find any genetic factor in the incidence of primary knee OA that correlated with serum levels of IL-1R1 and IL-1Ra.
Osteoarthritis has been considered a genetic disease for the past decade. Some studies have shown a correlation between genetics and osteoarthritis. Almost 100 DNA polymorphisms have been found to correlate with polymorphic DNA 20 . The DNA sequence is copied into an RNA sequence. This messenger RNA (mRNA) is used to produce the protein. A study performed by Shorter showed that miRNA can bind to mRNA and silence gene expression 21 . This results in a decrease or change in the resulting protein. Several studies have shown the important role of miRNA in regulating the homeostasis of joint tissues, especially in OA [22][23][24] . miRNA itself can suppress the expression of MMP13, IL-6, and IL-8 25 or can dampen IL-1β expression [26][27][28] . In this study, we found 3 SNPs that might correlate with OA. However, we did not find any difference in the serum levels of IL-1R1 or IL-1Ra. This finding can be due to several reasons. First, the sample source was serum. Second, miRNA itself can silence the genetic expression of SNPs, which can cause no difference in the protein levels of IL-1R1 and IL-1Ra. Third, the effect of miRNA in causing knee OA is stronger than the effect of the SNP itself. The correlation between miRNA and SNP expression needs to be evaluated.
The three gene variants analysed in our study were located in the intron region of the IL-1R1 gene, with no known function to date, and have not been linked to knee OA in previous studies. A study by Wu in 2016 found a correlation of the rs3771202 gene variant in IL-1R1 with severe preeclampsia 29 , and another study of rs3771202 conducted by Golozar in 2014 showed its correlation to oesophageal squamous cell carcinoma 30 . The other two variants (rs871659 and rs3917238) had never been found to be correlated with any condition. This finding showed the novelty of our study and could lead to further and larger genetic analysis studies in the general population. This finding can also be added to genomic data to be used for further research on knee OA worldwide.
The Hardy-Weinberg equilibrium results showed that these 2 new SNPs did not follow the HWE rule between the control group and the mild OA group. HWE states that genetic variation in a population will always be constant from one generation to the next if there are no external factors (such as mutations, recombination during sexual reproduction, genetic drift, gene migration/gene flow, and natural selection). This does not directly mean that there were external factors that influenced the appearance of these SNPs in the population (patients with mild OA) but may be due to the small number of samples used in this study. The sample size that we used in this study was only 165 (100 controls and 65 individuals with mild OA). To determine whether these two SNPs indeed deviate from the HWE and are not simply due to bias in the data, analysis with a larger sample size is needed.
This study has several limitations. First, in this study, we only focused on a small region of chromosome 2 depicting the IL-1R1 and IL-1RN genes and found four variants that were not correlated but found other 3 variants that were correlated. Only very small regions of the genome can be investigated at a time, meaning that important genes may be overlooked using this method. Second, we only examined IL-1Ra and IL-1R1 from serum and not from synovial fluid. We had difficulty in performing aspiration of the normal knee joint; therefore, we focused on blood examination. One study evaluated the concentration of IL-1β in serum and synovial fluid, and the results showed that the serum level of IL-1β was below the detection limit for RA and OA but not in synovial fluid 31 . Another study evaluated IL-1Ra in serum compared to synovial fluid in the temporomandibular joint. The study showed that the concentration of IL-1Ra in plasma was lower than that in synovial fluid 32 . This may be the reason we could not find any difference in the concentrations of plasma IL-1Ra and IL-1R1 within groups. Third, to examine the genetic effect, a larger sample size is needed. However, in this study with this number of samples, the post hoc power analysis was 84.3%, which showed no type II error in this study. There is no study on the SNP we evaluated and found in Indonesia. Therefore, we cannot evaluate the real homozygote that can be used as a reference. However, even if we decide which homozygote to use as a reference, it will not change the results, only the direction of the correlation. Fourth, we did not elaborate on the role of epigenetics of the interleukin-1 gene. The studies by Genemaras et al. 33 and Shen et al. 34 showed the implication of epigenetics in IL-1β and have been implicated in the progression of OA. A separate thorough study of epigenetics in interleukin-1 in the Indonesian population should be performed and analysed together with our results to provide a better understanding of why the gene variant found in this study was not correlated with clinical, radiographic or serum biomarkers.

Conclusions
In this study, we did not find the correlation between IL-1RN gene variants rs419598, rs9005, and rs315943 and IL-1R1 rs2287047 and knee OA. But we found three SNPs (IL-1R1_rs871659, IL-1R1_rs3771202, and IL-1R1_ rs3917238) that might be predisposing factors for primary knee OA. The SNP rs871659 allele A was a predisposing factor for knee OA among females, while the SNP rs3771202 allele G was a predisposing factor for knee OA among males. The evaluation of clinical findings (VAS, KOOS, and EQ-5D-3L), radiographic severity, and laboratory markers of serum IL-1Ra and IL-1R1 protein showed no association with the genetic polymorphisms of the IL-1R1 and IL-1RN genes. Studies on predisposing factors for primary knee OA are still inconclusive. Therefore, further studies with a larger sample size are needed, especially case-control studies.

Data availability
The data supporting the findings of this study are available from the corresponding author upon reasonable request.