Steroid hormone-related polymorphisms associate with the development of bone erosions in rheumatoid arthritis and help to predict disease progression: Results from the REPAIR consortium

Here, we assessed whether 41 SNPs within steroid hormone genes associated with erosive disease. The most relevant finding was the rheumatoid factor (RF)-specific effect of the CYP1B1, CYP2C9, ESR2, FcγR3A, and SHBG SNPs to modulate the risk of bone erosions (P = 0.004, 0.0007, 0.0002, 0.013 and 0.015) that was confirmed through meta-analysis of our data with those from the DREAM registry (P = 0.000081, 0.0022, 0.00074, 0.0067 and 0.0087, respectively). Mechanistically, we also found a gender-specific correlation of the CYP2C9rs1799853T/T genotype with serum vitamin D3 levels (P = 0.00085) and a modest effect on IL1β levels after stimulation of PBMCs or blood with LPS and PHA (P = 0.0057 and P = 0.0058). An overall haplotype analysis also showed an association of 3 ESR1 haplotypes with a reduced risk of erosive arthritis (P = 0.009, P = 0.002, and P = 0.002). Furthermore, we observed that the ESR2, ESR1 and FcγR3A SNPs influenced the immune response after stimulation of PBMCs or macrophages with LPS or Pam3Cys (P = 0.002, 0.0008, 0.0011 and 1.97•10−7). Finally, we found that a model built with steroid hormone-related SNPs significantly improved the prediction of erosive disease in seropositive patients (PRF+ = 2.46•10−8) whereas no prediction was detected in seronegative patients (PRF− = 0.36). Although the predictive ability of the model was substantially lower in the replication population (PRF+ = 0.014), we could confirm that CYP1B1 and CYP2C9 SNPs help to predict erosive disease in seropositive patients. These results are the first to suggest a RF-specific association of steroid hormone-related polymorphisms with erosive disease.

. Demographic and clinical characteristics of RA patients included in the discovery and replication cohorts. Data are means ± standard deviation or n (%). Abbreviations: RF, rheumatoid factor; ACPA: anticitrullinated protein antibodies; DAS28, disease activity score; DMARDs, disease-modifying anti-rheumatic drugs. P < 0.05 in bold. *RF and anti-CCP data were available in 809 and 673 RA patients in the discovery population, respectively. *RF and anti-CCP data were available in 564 and 477 RA patients with erosive disease in the discovery population, respectively. *RF and anti-CCP data were available in 245 and 196 RA patients without erosive disease in the discovery population, respectively *RF and anti-CCP data were available in 422 and 151 RA patients in the replication population, respectively. *RF and anti-CCP data were available in 300 and 109 RA patients with erosive disease in the replication population, respectively. *RF and anti-CCP data were available in 122 and 42 RA patients without erosive disease in the replication population, respectively. ∂ Information about methotrexate treatment was available in 761 and 416 patients in the discovery and replication populations, respectively. ∂ Information about methotrexate treatment was available in 524 and 293 patients with erosive disease in the discovery and replication populations, respectively. ∂ Information about methotrexate treatment was available in 237 and 123 patients without erosive disease in the discovery and replication populations, respectively.
anti-CCP-stratified analyses were also carried out and we included RF as interaction term in the overall logistic regression analysis to evaluate whether there was any statistically significant effect modification by these factors. Haplotype analysis using the same variables for adjustment was conducted using the R package Haplo.stats 45 . In order to facilitate eventual meta-analyses, the major allele was set as reference allele. All tests were conducted using the statistical software STATA (v.12) and R (http://www.r-project.org). In order to account for multiple testing, we set a P value of 0.00074 as significance study-wide threshold. The P value was calculated considering the number of independent polymorphisms analyzed (n = 34, MeffLi method) 46 but also the number of inheritance models tested (dominant and recessive).
Linkage disequilibrium (LD) and haplotype analysis. We performed haplotype frequency estimation and haplotype association analysis adjusted for age, sex and country of origin using the haplo.stats 45 . Haplotype frequencies were determined using the Expectation-maximization (EM) algorithm and haplotypes were reconstructed using SNPtools 47 and Haploview 48 . Block structures were determined according to the method of Gabriel et al. 49 ( Supplementary Fig. 1).  and RF-specific associations between SNPs and the risk of developing bone erosions, we genotyped the most interesting markers in a replication population from the DREAM registry consisting of 436 RA patients (307 RA patients with bone erosions and 129 patients without erosive disease). Demographic and clinical parameters of this population are also included in Table 1. We performed a meta-analysis of the data obtained in the discovery population with those from the DREAM registry and we pooled the Odds Ratios (ORs) for the most interesting polymorphisms using a fixed-effect model. Coefficients with a P-value ≤ 0.05 were considered significant. I 2 statistic was used to assess heterogeneity between studies.
functional analysis of the estrogen-related variants. Cytokine stimulation experiments were conducted in the 500 Functional Genomics (500FG) cohort from the Human Functional Genomics Project (HFGP; http://www.humanfunctionalgenomics.org/), which was designed to determine the influence of genomic variation on the variability of immune responses. The HFGP study was approved by the Arnhem-Nijmegen Ethical Committee (no. 42561.091.12) and biological specimens were collected after informed consent was obtained. We investigate whether any of the 41 estrogen-related SNPs correlated with cytokine levels (IFNγ, IL1β, IL6, TNFα, IL17, and IL22) after the stimulation of peripheral blood mononuclear cells (PBMCs), macrophages or whole blood from 408 healthy subjects with LPS (1 or 100 ng/ml), PHA (10 μg/ml), and Pam3Cys (10 μg/ml). After log transformation, linear regression analyses adjusted for age and sex were used to determine the correlation of selected SNPs with cytokine expression quantitative trait loci (cQTLs). All analyses were performed using R software (http://www.r-project.org/). In order to account for multiple comparisons, we used a significant threshold of 0.00025 (0.05/34 independent SNPs × 6 cytokines). Details on PBMCs isolation, macrophage differentiation and stimulation assays have been reported elsewhere [50][51][52] . Briefly, PBMCs were washed twice in saline and suspended in medium (RPMI 1640) supplemented with gentamicin (10 mg/mL), L-glutamine (10 mM) and pyruvate (10 mM). PBMC stimulations were performed with 5 × 10 5 cells/well in round-bottom 96-wells plates (Greiner) for 24 hours in the presence of 10% human pool serum at 37 °C and 5% CO 2 . Supernatants were collected and stored in −20 °C until used for ELISA. LPS (100 ng/ ml), PHA (10 μg/ml) and Pam3Cys (10 μg/ml) were used as stimulators for 24 or 48 hours. Whole blood stimulation experiments were conducted using 100 μl of heparin blood that was added to a 48 well plate and subsequently stimulated with 400 μl of LPS and PHA (final volume 500 ul) for 48 hours at 37 °C and 5% CO 2 . Supernatants were collected and stored in −20 °C until used for ELISA. Concentrations of human IFNγ, IL1β, IL6, TNFα, IL17, and IL22 were determined using specific commercial ELISA kits (PeliKine Compact, Amsterdam, or R&D Systems), in accordance with the manufacturer's instructions.
Once we assessed the correlation of estrogen-related SNPs with cytokine levels, we used the HaploReg SNP annotation tool (http://www.broadinstitute.org/mammals/haploreg/haploreg.php) to further investigate the functional consequences of each specific variant. We also assessed whether any of the potentially interesting markers correlated with mRNA expression levels of their respective genes using data from public eQTL browsers (GTex portal; www.gtexportal.org/home/ and Blood eQTL browser; https://genenetwork.nl/bloodeqtlbrowser/). correlation between steroid hormone levels and hormone-related Snps. We also measured serum levels of seven steroid hormones (androstenedione, cortisol, 11-deoxy-cortisol, 17-hydroxy progesterone, progesterone, testosterone and 25 hydroxy vitamin D3) in the 500FG cohort, which includes 531 healthy subjects. Complete protocol details of steroid hormone measurements have been reported elsewhere 52 . Hormone levels and genotyping data were available for a total of 406 subjects.
After log-transform, correlation between steroid hormone levels and steroid hormone-related SNPs was evaluated by linear regression analysis adjusted for age and sex (or for age when men and women were analysed separately). In order to avoid a possible bias, we excluded from the analysis those subjects that were using oral contraceptives or those subjects in which this information was not available. A total of 279 healthy subjects (107 women and 272 men) were finally available for analysis. Significance threshold was set to 0.00021 considering the number of independent SNPs tested (n = 34) and the number of hormones determined (n = 7). predictive models and discriminative accuracy. The value of steroid hormone-related variants for prediction of prognosis and disease progression in seropositive and seronegative RA patients was assessed using stepwise logistic regression. Models were built including demographic variables (age and sex) and genetic polymorphisms that showed significant associations with erosive disease in the single-SNP analysis (P < 0.05). The genetic model was then compared with the reference model including demographic variables. The area under the curve (AUC) of a receiver operating characteristic (ROC) curve analysis and −2 log likelihood ratio (LR) tests were used to assess whether the genetic models fitted significantly better the data compared to their respective reference models. Finally, we run randomization tests to confirm whether the improved predictive ability of each genetic model was consistent after 50.000 iterations. All tests were conducted using R software (http:// www.r-project.org/). ethics approval. The

Results
Erosive RA patients had a similar age than those patients without bone erosions (59.66 ± 12.47 vs. 58.95 ± 14.30, P = 0.50) and had a significantly higher female to male ratio (462/105 = 4.4 vs. 182/67 = 2.71, P = 0.007; Table 1). Overall, the percentage of RA patients with positive RF and anti-CCP was 70.58% and 72.80% respectively, and these percentages were slightly higher among those patients with bone erosions (72.52% and 74.21%) than in those without erosive disease (66.12% and 69.39%). The mean of disease follow-up was 18.30 years whereas the mean of DAS28 was 5.63. Six hundred and three patients received methotrexate (79.24%) and 632 patients (77.45%) were treated with biologic therapies. With the exception of gender, none of the demographical or clinical variables significantly differ between patient with and without erosive disease ( Table 1).
Association of steroid hormone-related polymorphisms with the risk of having bone erosions. Selected polymorphisms did not show deviation from HWE in the control population (patients without erosive disease; P > 0.001). Logistic regression analysis adjusted for age, gender and country of origin revealed that carriers of the ESR2 rs1271572T/T genotype tended to have a decreased risk of developing erosive disease than those subjects carrying the G allele (OR = 0.55, P = 0.004: Table 3). Although the association of the ESR2 rs1271572T/T genotype with a decreased risk of having bone erosions remained only marginally significant after correction for multiple testing, we found a significant RF-specific effect of this SNP to modulate the risk of having erosive disease. Seropositive patients carrying the ESR2 rs1271572T/T genotype had a significantly reduced risk of developing erosive disease (OR = 0.38, P = 0.0002) whereas an opposite but not significant effect was found in seronegative patients (OR = 1.08, P = 0.83; P Int = 0.018; Table 3). Importantly, the association of the ESR2 rs1271572 SNP with a reduced risk of developing bone erosions in seropositive patients remained significant after correction for multiple testing (P < 0.00074). Although the association was not replicated in the DREAM cohort, the meta-analysis of our data and those from the DREAM registry, including 1252 RA patients, confirmed that seropositive patients carrying the ESR2 rs1271572T/T genotype had a decreased risk of developing erosive disease (OR RF+ = 0.52, P = 0.00074) whereas a totally opposite but not statistically significant effect was found in seronegative patients (OR RF− = 1.28, P = 0.42; P Het = 0.33; Table 4). No significant anti-CCP effect modification was found for this SNP to modulate the risk of developing bone erosions (P Int = 0.11; Supplementary Table 1), which suggest that ESR2 locus might play a relevant role in determining disease progression in a RF-dependent manner. In agreement with these results, we found a RF-specific effect of the ESR2 CGTA haplotype to determine the risk of developing erosive disease. Seropositive RA patients carrying the ESR2 CGTA haplotype (and, therefore, not harboring the ESR2 rs1271572 protective allele) showed an increased risk of developing bone erosions (OR RF+ = 1.63, P = 0.0051) whereas an opposite but not significant effect was detected in seronegative patients (OR RF− = 0.93, P = 0.99; Supplementary Table 2). An overall haplotype analysis including 1252 RA patients from the discovery and replication populations confirmed the RF-specific association of the ESR2 CGTA haplotype with an increased risk of developing bone erosions (OR RF+ = 1.44, 95%CI 1.13-1.84; P = 0.0036 and OR RF− = 0.89, 95%CI 0.62-1.26; P = 0.51). According to publicly available gene expression datasets (GTex portal and Haploreg), the ESR2 rs1271572 variant strongly correlate with ESR2 mRNA expression levels in whole peripheral blood (P = 3.1•10 −9 ) but also in primary B cells, lymphoblastoid cell lines (from P = 1.98•10 −6 to P = 3.47•10 −10 ) and several tissues (ranging from P = 2.60•10 −5 to P = 8.33•10 −23 ; Supplementary Table 3). Intriguingly, a similar level of correlation with gene expression was found for other variants belonging to the ESR2 CGTA haplotype (Supplementary Table 3), which strongly suggested that the ESR2 rs1271572 SNP or ESR2 CGTA haplotype might represent an eQTL for ESR2.
Similarly, we also found a RF-specific effect of the CYP2C9 rs1799853 and CYP1B1 rs10012 SNPs to determine the risk of developing bone erosions. Seropositive patients carrying the CYP2C9 rs1799853T/T or CYP1B1 rs10012G/G genotypes had a significantly reduced chance of developing bone erosions (OR = 0.16, P = 0.0007 and OR = 0.42, P = 0.0040) whereas an opposite but not significant effect was observed in seronegative patients (OR = 2.71, P = 0.23; P Int = 0.003 and OR = 1.70, P = 0.35; P Int = 0.031; Table 3). The effect of the CYP2C9 rs1799853 polymorphism on the risk of developing erosive disease in seropositive patients remained statistically significant after correction for multiple testing (P < 0.00074), which suggested a role of the CYP2C9 gene in modulating disease progression in RA. In accordance with these findings, we found that seropositive patients carrying the CYP2C9 AT haplotype had a significantly decreased risk of developing erosive disease (OR RF+ = 0.61, P = 0.0075) whereas no effect was observed in seronegative patients (OR RF− = 0.87, P = 0.57; Supplementary Table 2). No significant anti-CCP effect modification was found for CYP2C9 and CYP1B1 variants to determine the appearance of bone erosions (P Int = 0.88 and P Int = 0.27) underlying the importance of considering RF when evaluating the impact of the CYP2C9 and CYP1B1 loci on the risk of developing erosive disease. Importantly, when we attempted to validate the RF-specific association of the CYP1B1 rs10012G/G genotype with a decreased risk of having erosive disease in the replication population, we found that seropositive patients carrying the CYP1B1 rs10012G/G genotype had a significantly decreased risk of developing bone erosions (OR RF+ = 0.30, P = 0.0051) whereas an opposite but not significant effect was found in seronegative patients (OR RF− = 5.97, P = 0.10; P Int = 0.012; Table 4). The meta-analysis of both populations confirmed the strong RF-specific effect of this SNP to determine the risk of developing bone erosions (OR RF+ = 0.38, P = 0.000081; P Het = 0.52 vs. OR RF− = 2.22, P = 0.11; P Het = 0.31). Although we attempted to validate the association of the CYP2C9 rs1799853T/T genotype with a decreased risk of having erosive disease, the relatively small size of the replication population did not allow us to perform the association analysis according to a recessive model of inheritance. However, we found a RF-specific effect on the risk of having erosive disease for a neighboring SNP (rs1057910), which suggested a RF-dependent effect of the CYP2C9 locus to modulate the risk of erosive disease (OR RF+ = 2.75, P = 0.027 vs. OR RF− = 0.54, P = 0.47; Table 4). The meta-analysis of both cohorts confirmed the RF-specific effect of the CYP2C9 rs1057910 SNP on the risk of developing bone erosions (OR RF+ = 2.68, P = 0.0022 vs. OR RF− = 1.08, P = 0.83; P Het = 0.34; Table 4).
In line with these findings, we also observed an additional RF effect modification of the FcγR3A rs396991 and SHBG rs6259 SNPs to determine the risk of having erosions. Seronegative patients carrying the FcγR3A rs396991C (2019) 9:14812 | https://doi.org/10.1038/s41598-019-51255-0 www.nature.com/scientificreports www.nature.com/scientificreports/ allele had a significantly reduced chance of developing bone erosions (OR = 0.45, P = 0.013) whereas an opposite but not significant effect was observed in seropositive patients (OR = 1.18, P = 0.46; P Int = 0.028; Table 3). Furthermore, seropositive subjects carrying the SHBG rs6259A allele showed an increased risk of developing bone erosions (OR RF+ = 1.87, P = 0.015) whereas an opposite but not significant effect was detected in seronegative patients (OR RF− = 0.66, P = 0.19). Interestingly, although the effect was stronger in seronegative patients, we could validate the RF-specific effect of the SHBG rs6259 SNP on the risk of developing erosions in the replication population (OR RF+ = 1.17, P = 0.63 vs. OR RF− = 0.22, P = 0.009; P Int = 0.013; Table 4) and the meta-analysis of both cohorts confirmed that the effect of this marker was dependent on the RF status (OR RF+ = 1.55, P = 0.033 vs. OR RF− = 0.48, P = 0.0087; P Het = 0.14; Table 4). Although we could not validate the RF-specific association of the FcγR3A rs396991 SNP with bone erosions in the DREAM registry, the meta-analysis of both cohorts confirmed   www.nature.com/scientificreports www.nature.com/scientificreports/ the RF-specific effect of this SNP to modulate the risk of developing erosive disease (OR RF− = 0.47, P = 0.0067 vs. OR RF+ = 1.02, P = 0.93). None of these two SNPs showed a significant effect modification by anti-CCP (P Int = 0.85) suggesting again that RF, rather than anti-CCP, is a driver factor influencing the impact of the steroid hormone-related loci on disease progression in RA.
Finally, an overall association analysis revealed that carriers of the ESR1 rs1801132G allele showed a decreased risk of developing bone erosions (OR = 0.71, P = 0.034). Although we could not validate this association in the replication population, we found that this SNP showed a significant RF-specific effect to modulate the risk of developing bone erosions but according to a recessive model of inheritance. Thus, seropositive carriers of the ESR1 rs1801132G/G genotype showed a decreased risk of developing bone erosions (OR RF+ = 0.39, P = 0.004) whereas an opposite but not statistically significant effect was observed in seronegative subjects (OR RF− = 1.43, P = 0.57; Table 4). Furthermore, we found a similar RF-specific effect for the ESR1 rs9340799 SNP that was not detected in the discovery population (OR RF+ = 0.42, P = 0.009 vs OR RF− = 8.33, P = 0.011). Considering that none of these associations survived after correction for multiple testing and that the effect of ESR1 SNPs on the risk of developing erosive disease seemed to depend on the inheritance model applied, these results suggested a complex relationship between the ESR1 locus and bone erosion probably mediated by more than one SNP. In support of this notion, we found that 3 large ESR1 haplotypes (ESR1 CTATTTTCTA , ESR1 CTATTCTTCA , and ESR1 GTATTCTCTA ) were significantly associated with a decreased risk of having erosive disease (P = 0.0094, P = 0.0021 and P = 0.0023, respectively; Supplementary Table 2). correlation of selected Snps with steroid hormone levels. Besides the strong genetic association with erosive disease identified for the CYP2C9 rs1799853 SNP and its known role in controlling the metabolism of a wide range of drugs (with the T allele acting as poor metabolizer), we found that this coding variant strongly correlated with serum vitamin D3 levels in women (P = 0.00085 and P = 0.0019; Fig. 1) whereas no effect was detected in men. Although the association of the CYP2C9 rs1799853 SNP with reduced levels of vitamin D3 in women remained borderline significant, this finding suggested that this marker might have a role in the modulation of bone homeostasis and vitamin D3-mediated immune responses.
On the other hand, we found that the ESR1 rs2881766G/G genotype or G allele weakly correlated with progesterone levels (P = 0.0076 and P = 0.0071) and that the ESR1 rs851984, ESR1 rs2077647 , ESR1 rs2071454 , ESR1 rs3798577 and ESR1 rs910416 variants mapped among histone marks in several cell types including osteoblasts and a wide variety of immune cells.
correlation of steroid hormone Snps with cytokine levels. Interestingly, we also found that carriers of the ESR2 rs4986938T allele had reduced levels of TNFα after the stimulation of PBMCs with Pam3Cys for 24 h (P = 0.0022; Fig. 2A). These results along with those reporting that the ESR2 rs4986938 and ESR2 rs1271572 SNPs map among histone marks in multiple cell types including osteoblasts and different subsets of immune cells, suggest a possible functional role of the ESR2 variants in modulating the risk of developing bone erosions likely through the modulation of ESR2 expression. In addition, we found that the presence of the CYP2C9 rs1799853T allele correlated with an increased production of IL1β after the stimulation of PBMCs with LPS or PHA for 24 h or 48 h (P = 0.0057 and P = 0.0058; Fig. 2B,C), which also pointed to a functional role of this marker in determining the presence of bone erosions.
In addition, we found that the ESR1 rs3798577 variant correlated with TNFα and IL6 levels after the stimulation of human macrophages with LPS for 24 h (P = 0.00083 and 0.0011; Fig. 2D,E). Finally, we found that carriers of the FcγR3A rs396991C allele showed a significantly increased production of TNFα after stimulation of macrophages with LPS for 24 h (P = 1.97•10 −7 ; Fig. 2F). Of note, the association of the FcγR3A rs396991C allele with increased levels of TNFα in macrophages survived after correction for multiple testing, which strongly suggested a functional effect of this variant to modulate macrophage-mediated immune responses, a key factor influencing the risk of developing erosive disease. On the contrary, although it was tempting to speculate that ESR1, ESR2, CYP2C9 SNPs might also exert their effect on the risk of developing erosive disease through the modulation of steroid hormones or steroid hormone-mediated immune responses, it is important to mention that none of the associations between ESR1, ESR2 or CYP2C9 SNPs and cytokine levels survived after correction for multiple testing, which suggested a modest functional impact of these polymorphisms on the risk of developing bone erosions.
Usefulness of steroid hormone-related Snps to predict erosive disease. As a whole, our data suggest that the attributable effect of the CYP1B1, CYP2C9, ESR1, ESR2, SHBG, and FcγR3A loci to modulate the risk of developing bone erosions in RA patients might be dependent on the presence of either missense or intronic polymorphisms that affect the immune responses to a greater or lesser extent. Considering the strength of the RF-specific associations found for SNPs within CYP1B1, CYP2C9, SHBG, ESR1, ESR2, FcγR3A and CYP3A4 loci in the discovery and/or replication populations, we decided to assess whether SNPs within these loci could be useful to differentially predict disease progression in seropositive and seronegative patients. Our results showed that the addition of 5 steroid hormone-related SNPs within the CYP1B1, CYP2C9, CYP3A4, ESR2 and SHBG loci to a model including demographic variables significantly improved the ability to predict the appearance of bone erosions in seropositive patients (AUC Genetic = 0.73 vs. AUC Demographic = 0.63; P = 2.46•10 −8 ; Table 5) whereas no significant predictive value was found for these SNPs in seronegative patients (AUC Genetic = 0.61 vs. AUC Demographic = 0.59; P = 0.36; Fig. 3). The consistency of this predictive analysis was confirmed through a permutation test that showed that none of the 50.000 permuted models for each group showed a better prediction capacity than the genetic model (Average sorted AUC = 0.644, Z-score = 6.79 and P Z_score (50.000perm) = 5.67•10 −12 ). Even though the lack of patients carrying the CYP2C9 rs1799853T/T genotype and the relatively small size of the replication population hampered the validation of this 5-SNP model to predict erosive disease, we attempted

Discussion
The present study reports, for the first time, both overall and RF-specific associations of steroid hormone-related polymorphisms with the risk of developing erosive RA. The most relevant effect was found for SNPs within CYP1B1, CYP2C9, and ESR2 genes. We observed that seropositive RA patients carrying the CYP1B1 rs10012G/G , CYP2C9 rs1799853T/T , and ESR2 rs1271572T/T , genotypes had a significantly reduced risk of developing bone erosions during the course of the disease whereas an opposite but not significant effect was found in seronegative patients. Although the relatively small size of the replication population hampered the validation of these associations according to a recessive model of inheritance, we could validate the RF-specific association of the CYP1B1 rs10012 SNP with the risk of developing erosive disease in the replication population and the meta-analysis of the discovery and replication cohorts confirmed the strong RF effect modification of this SNP to determine the risk of bone erosions. In addition, although we could not validate the RF-specific association of the CYP2C9 rs1799853 variant in the replication population due to the lack of patients carrying the CYP2C9 rs1799853T/T genotype, we found a RF-specific effect on the risk of having erosive disease for a neighbouring SNP within the CYP2C9 locus (rs1057910) that was further confirmed through meta-analysis. Although this SNP was not in linkage disequilibrium (LD) with the rs1799853 and, therefore, does not represent the same association signal, these results support the idea that the CYP2C9 locus might influence the risk of developing bone erosions in a RF dependent manner and likely through the modulation of the hormone metabolism and hormone-dependent immune responses.
Whilst the CYP1B1 locus is located on chromosome 2p21-22, CYP2C9 belongs to the CYP2C family, a gene cluster (CYP2C19-CYP2C9-CYP2C8) located on chromosome 10q23.3. Together with CYP1A2 and CYP3A4,  Table 4. Replication of the most interesting associations between estrogen-related polymorphisms and risk of developing erosive disease (DREAM registry) and meta-analysis. Abbreviations: SNP, single nucleotide polymorphism; OR, odds ratio; CI, confidence interval. Data on RF was available in 422 RA patients in the replication population. Estimates were adjusted for age and sex. P ≤ 0.05 in bold. The RF-specific effect modification of steroid hormone SNPs was determined by logistic regression using RF as interaction term. Meta-analysis was conducted following a fixed effect model. † Estimates calculated according to a dominant model of inheritance. § Estimates calculated according to a recessive model of inheritance.
Scientific RepoRtS | (2019) 9:14812 | https://doi.org/10.1038/s41598-019-51255-0 www.nature.com/scientificreports www.nature.com/scientificreports/ CYP1B1 and CYP2C9 catalyze the conversion of estrogens to genotoxic catechol estrogens (estradiol 4-and 2-hydroxylation, respectively) 53 , which are key processes that allow the binding of catechol estrogens to ESR1 and ESR2. At low concentrations, CYP2C9 is also implicated in the 17beta-hydroxy dehydrogenation of estradiol creating estrone, which is one of the 3 natural estrogens with multiple immunomodulatory actions. Given that both the CYP1B1 rs10012 and CYP2C9 rs1799853 SNPs are coding variants that alters their respective protein amino acid sequences (R48G and Arg144Cys) and appear to decrease the activity of the enzyme but also proper folding and stability, it seems to be plausible to hypothesize that the presence of these SNPs could determinate estrogen-dependent immune responses and thereby modulate the risk of developing bone erosions. Our functional studies also demonstrated that the CYP2C9 rs1799853 SNP correlated with serum vitamin D3 levels, which suggest that the CYP2C9 rs1799853 SNP might also affect disease progression through the regulation of vitamin D3-mediated immune responses. However, we need to interpret these results with caution as we only found a significant correlation with vitamin D3 in women whereas no effect was seen in men. These results, together with those reporting that carriers of the CYP2C9 rs1799853T allele have an increased production of IL1β after the stimulation of PBMCs or whole blood with LPS or PHA, suggest that the protective effect attributed to this coding variant might not only depend on vitamin D3 but other factors such as RF or even other stimuli or substrates. In this regard, our team has reported that the CYP2C9 rs1799853 polymorphism is strongly associated with poor response to anti-TNF drugs in RA 54 , suggesting that this missense variant might modulate the strength of immune responses through the regulation of the metabolism of endogenous compounds but also compounds administered exogenously.
Although we also attempted to validate the RF-specific association of the ESR2 rs1271572T/T genotype with a decreased risk of having erosion in the replication population, we only found a modest and not significant effect of this variant to determine erosive disease. However, the direction of the effect in seropositive and seronegative patients was similar to the one observed in the discovery population and the meta-analysis of both cohorts confirmed that the effect of this SNP on the risk of developing bone erosions was modified by RF. In support of a RF-specific effect of this variant to influence the risk of erosive RA, we found that seropositive carriers of the ESR2 CGTA haplotype had a decreased risk to develop erosive RA whereas no effect was detected in seronegative patients. Interestingly, an overall haplotype analysis also revealed a significant association of 3 common The ESR2 and ESR1 genes (14q23.2 and 6q25.1 respectively) encode the estrogen receptor beta (ESRβ) and alpha (ESRα) that are highly expressed in synovial cells 55 and bone 56 but also in most of the immune cells 57 . Although a number of experimental studies have shown that female RA patients have worse prognosis and higher disease activity and health assessment questionnaire scores in comparison with male patients, it is also well established that steroid hormones have both pro-and anti-inflammatory roles in RA. Although it has been reported, for instance, that the activation of ESRs by estradiol (E 2 ) often leads to joint protection and the maintenance of bone density (by inhibiting bone resorption) 58 and that the withdrawal of estrogens drastically increases the severity of the disease (by promoting joint destruction, bone erosions and physical disability) 59 , it has been also reported that RA patients have high levels of estrone in the synovial fluid compared to healthy individuals and that estrogens can also induce pro-inflammatory responses through the activation of different mechanisms involving humoral immunity 60 , multiple transcription factors (such as c/EBPβ, STAT-1, NFκB) and oxidative stress pathways (especially those involving iNOs) 61,62 . Furthermore, it has been reported that estrogens are able to promote pro-inflammatory pathways including B-and T-cell proliferation 63 , thymocyte maturation 64 , cell trafficking 65 and the expression of specific adhesion molecules 63 . Although the existing paradox with respect to the immunomodulating role of steroid hormones in RA remain unresolved, it seems to be reasonable to hypothesize that the presence of ESR2 polymorphisms that correlate either with ESR2 mRNA expression levels may influence on the risk of developing bone erosions in RA likely through the modulation of ESR2-dependent tolerogenic immune responses. In addition, although the association of the ESR2 and ESR1 polymorphisms with serum hormone levels or TNF and IL6 levels in stimulated macrophages did not remain significant after correction for multiple testing, our functional findings were in agreement with the genetic results suggesting a protective effect of ESR2 polymorphisms and specific ESR1 haplotypes on the risk of developing erosive RA. In addition, our genetic and functional results were also concordant with data of previous studies reporting that the presence of certain SNPs, microsatellites or even specific haplotypes within estrogen receptor genes is associated with bone mineral density and influences the risk of developing bone erosions 66,67 affecting RA patients 68 but also subjects diagnosed with other chronic inflammatory diseases 69 and bone degenerative diseases 70,71 .
Finally, this study also showed a weak but still interesting RF-specific effect of the SHBG rs6259 and FcγR3A rs396991 SNPs to determine the risk of having erosions. Importantly, we could validate the RF-specific effect of the SHBG rs6259 SNP on the risk of developing erosions in the replication population and the meta-analysis of both cohorts confirmed that the effect of this marker was strongly dependent on the RF status. On the other hand, although the RF-specific association of the FcγR3A rs396991 SNP with bone erosions was not statistically significant in the DREAM registry, the meta-analysis of both cohorts confirmed that the effect of this SNP on the risk of bone erosions was modified by the RF.
Whereas little is known about the role of the SHBG locus (17p13) in determining RA progression, a number of experimental studies have shown that the FcγR3A locus (1q23) is involved in the recognition of IgG1 and IgG3 by NK cells and macrophages and that the activation of this receptor by IgG and IgG-RF immunocomplexes might lead to the initiation of a range of sustained and harmful inflammation events that, if chronified, may cause joint and bone destruction and promote the onset of RA [72][73][74] . In this context and considering the number of studies reporting association of the FcγR3A rs396991 SNP with autoimmune diseases [75][76][77][78][79][80][81][82] , the response to a wide range of biologic drugs [83][84][85][86][87][88][89] and an exacerbated production of TNFα after stimulation of macrophages with LPS for 24 h but also the reported association of the SHBG rs6259 with serum SHBG levels 90 , we hypothesize that these SNPs might also play a relevant role in determining bone erosions and disease progression.
Considering the noticeable RF-specific impact of the CYP1B1, CYP2C9, ESR2, FcγR3A and SHBG SNPs on the risk of developing erosive disease but also their functional implication in modulating hormone levels and/ or immune responses, we decided to assess if the presence of steroid hormone-related SNPs could be useful to reliably predict the appearance of bone erosions in seropositive and seronegative patients separately. To do that, we built a genetic model including demographic variables and those SNPs that were consistently associated with the risk of developing bone erosions in seropositive patients. After removing the SNPs that were not significantly associated with erosive disease in the model, we obtained a model including 5 SNPs within the CYP1B1, CYP2C9, CYP3A4, ESR2, and SHBG loci that significantly improved the ability to predict the risk of developing erosive disease when compared with a reference model including demographic variables. The predictive capacity of these SNPs was restricted to seropositive patients since the addition of the same SNPs (or any other genetic marker) to a model built with demographic variables in seronegative patients did not show any predictive value. The predictive   Table 5. Discriminative value AUC for the model including estrogen-related variants in the discovery and replication populations. a Including age and gender as variables never dropped from models and when are compared with a baseline model with AUROC = 0.5. P ≤ 0.10 in bold (stepwise threshold). *All SNPs showing a significant association with erosive disease (P < 0.10) were initially added to the model in the discovery population. ‡ A sort analysis in the discovery population revealed that this model showed an AUC value systematically higher than those observed in 50.000 randomized models: Average AUC of null distribution (50.000 models) = 0.644 Z score = 6.79, P Z_score-value_(50.000perm) = 5.67•10 −12 . Ϯ All SNPs were forced to be included in the replication population with the exception of the CYP2C9 rs1057910 that was included due to the impossibility to calculate association estimates for the CYP2C9 rs1799853 SNP. (2019) 9:14812 | https://doi.org/10.1038/s41598-019-51255-0 www.nature.com/scientificreports www.nature.com/scientificreports/ ability of the genetic model in seropositive patients was consistent as no similar models were found after performing a 50.000 permutation test. When we attempted to confirm the utility of this model in the DREAM registry, we found that the CYP1B1 and CYP2C9 SNPs in seropositive patients showed a consistent predictive value for the development of bone erosions. These results suggest that CYP1B1 and CYP2C9 SNPs alone or in combination with other clinical and genetic markers might help to improve the ability to predict the appearance of bone erosions in seropositive patients (~70% of RA patients). Additional studies including these and other genetic and clinical markers are urgently needed to improve our ability to predict disease progression in RA.
This study has strengths and weaknesses. The strengths of this study include a relatively large and well-characterized population and the meta-analyses conducted considering results from the DREAM registry. In the discovery population, we had 80% of power to detect an odds ratio of 1.68 (α = 0.00074) for a SNP with a frequency of 0.25, which underlined the feasibility of the study design. Another important strength of this study is the development of cytokine stimulation experiments and the measurement of seven serum steroid hormones in a large cohort of healthy subjects, which allowed us to investigate the functional role of the most interesting markers in modulating immune responses but also to test their impact on determining steroid hormone levels. A drawback is the multicenter nature of this study that placed inevitable limitations such as the impossibility of using available scores to better define bone erosions (Sharp van der Heijde, Genant, SENSE, and Ratingen scores). Given the cross-sectional approach of the study, we had also intrinsic limitations such as a possible bias due to variations in treatments and follow-up time among study participants. Finally, it is important to mention that the selection of SNPs for this study was influenced by the limited research funds and that the relatively small size of ROC curves summarize the accuracy of prediction for genetic and demographic models in seropositive and seronegative patients. The genetic models (marked in blue) included SNPs that were significantly associated with erosive disease in seropositive patients (either in the single-SNP or haplotype analyses) whereas the demographic models included demographic variables (age and gender as covariates; marked in green) for seropositive and seronegative patients.