Abstract
Rheumatoid arthritis (RA) has an intricate etiology that includes environmental factors as well as genetics. Organophosphate esters (OPEs) are frequently used as chemical additives in many personal care products and household items. However, there has been limited research on their potential effects on rheumatoid arthritis (RA). The specific associations between OPEs and RA remain largely unexplored. This study investigates any potential associations between adult rheumatoid arthritis risk and exposure to OPEs. We investigated data from the National Health and Nutrition Examination Survey (NHANES) 2011–2018 among participants over 20 years old. In two models, multivariable logistic regression was utilized to investigate the relationship between exposure to OPEs and RA. Furthermore, subgroup analyses stratified by age, gender, and dose exposure response were evaluated. Generalized additive models and smooth curve fits were used to characterize the nonlinear relationship between RA and OPEs. In conclusion, 5490 individuals (RA: 319, Non-RA: 5171) were analyzed. Higher quantiles (Q4) of DPHP and DBUP showed a higher prevalence of RA than the lowest quantiles. Our findings show that adult RA prevalence is higher in those who have been exposed to OPEs (DPHP, DBUP). Interestingly, these correlations seem to be stronger among women, the elderly, those with higher BMIs, and those who have diabetes. The dose-response curve for DPHP and DBUP demonstrated an upward-sloping trend. In contrast, BCEP and BCPP showed a U-shaped relationship and an inverted U-shaped relationship with the probability of RA.BDCPP demonstrated a complex relationship with a peak at lower concentrations followed by a decrease. Finally, our study also concludes that exposure to OPEs plays a crucial role in the pathogenesis of RA.
Highlights
Urinary metabolites of OPEs (DPHP, DBUP) are associated with the increased prevalence of RA.
BDCPP appears to have a protective effect, reducing the risk of RA.
Associations are more prominent in females, participants over 60 years, and those with higher BMI.
Exposure to OPEs plays a critical role in the pathogenesis of rheumatoid arthritis (RA).
1. Introduction
Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disorder that can cause irreversible joint damage and significant disability [1]. It is characterized by inflammatory changes of the synovial tissues of joints, cartilage, and bone [2]. Although the exact pathogenesis of RA is still unknown, it has become evident that the incidence and development of RA is a multi-factorial interaction and is closely associated with major genetic and environmental factors or stresses [3]. According to previous studies, genetic and epigenetic influences accounted for approximately 50–60 % of the risk for developing RA, while the remainder can be explained by environmental exposures such as dust exposure, smoking, and especially the microbiome, which also represents an “internal” environment [4]. There appears to be an imperative reciprocation between the adaptive and innate immune system components. The association between environmental exposure and prevalence of RA is well established [5, 6, 7, 8]. Recent studies confirmed the augmenting evidence and substantiated that environmental pollutants like polycyclic aromatic hydrocarbons, phthalates, heavy metals, and volatile organic compounds play a quintessential role in the prevalence of RA[9]. Notably, exposures to pesticides, insecticides, and traffic-related pollutants have also been implicated in the heightened risk of RA, accentuating the intricate link between environmental pollution and autoimmune disease susceptibility [3]. While RA is a lifelong condition, it can be controlled and prevented with some lifestyle modification, including regular exercise and a nutritious diet [10]. Polybrominated diphenyl ethers (PBDEs) were banned in 2004 by the European Commission due to their environmental persistence, bio-accumulation, and associated toxicity. Organophosphate esters were introduced as an environmentally safer alternative to PBDEs [11]. OPEs are organic phosphoric acid esters containing ether alkyl chains or aryl groups. They are either halogenated or nonhalogenated. Halogenated organophosphates are used as flame retardants, while nonhalogenated organophosphates are mostly used as plasticizers in consumer products, textiles, and construction materials [12]. OPEs are not chemically bound in the products they are added to, and due to this, they are released into the environment via volatilization, leaching, and/or abrasion [11]. Human exposure to OPEs stems from various sources, including inhalation of dust, unintentional ingestion of contaminated foods, and dermal contact. Indoor environments containing these chemicals serve as the most common source of exposure [13]. Epidemiological studies revealed the adverse effects of OPEs on human health, such as endocrine disturbance, neurotoxicity, and carcinogenicity [14]. Numerous adverse effects of exposure to OPEs have been identified in both experimental and human observational studies [15]. Decreased total-body bone mineral density (BMD) and lumbar spine BMD in zebrafish and in vitro models, as well as exposure disrupted thyroid and sex hormones. Watanabe et al. showed that exposure to OPEs can increase the level of IFN-γ in the bodily fluid of mice [16]. OPEs exposure activated interferon signaling pathways in rats that can affect at both at the level of biogenesis and bioactivity of immune ligands. Recent experimental studies provided significant insights into the linkage between OPEs exposure causing immunotoxicity and immunomodulatory effects on immune cells and humoral mediators [17]. The human body metabolizes OPEs rapidly and excretes them mainly through the urine, where diesters and monoesters are representative biomarkers [18]. OPEs can damage the renal structure of mice through the NF-κB p65/NLRP3-mediated inflammatory response [19]. Although experimental studies demonstrate exposure to OPEs may alter the production of RA biomarkers cytokines IL-1β, IL-6, and TNF-α in macrophages. However, little is known about the relationship between OPEs exposure and RA in humans.
The present research conducts a cross-sectional study based on U.S. national population data from the National Health and Nutrition Examination Survey (NHANES), 2011 to 2018. The impact of urinary metabolites on OPEs was analyzed using multivariate logistic regression. Furthermore, weighted generalized additive models and restricted cubic splines were utilized to determine the non-linear relation between urinary metabolites and RA and identify the dose-response curve of each metabolite. Our results provide novel epidemiological evidence on the associations of urinary metabolites of OPEs exposures with RA risk and contribute to identifying the hazardous factors of RA.
2. Materials and Methods
2.1. Study design and population
The National Health and Nutrition Examination Survey (NHANES) is a cross-sectional survey conducted by the National Centre for Health Statistics (NCHS) program of the Center for Disease Control (CDC), which is designed to assess the nutritional and health status of the non-institutionalized US population. A complex, stratified multistage sampling design is applied to make the sample more representative. The survey is conducted every two years, and information covering demographic, dietary, examination, laboratory, and questionnaire data was collected. In our study, urinary metabolites of OPEs were analyzed from survey cycles (2011-2012, 2013-2014, 2015-2016, 2017-2018). All the participants enrolled were provided informed consent at the time of recruitment for the survey. Detailed information about the participant’s demographic, socio-economic, and dietary factors, as well as medical and health status, is collected through in-home interviews. A standardized medical and physical examination with bio-specimen (such as serum and urine) collection is conducted by highly trained medical personnel in specially equipped mobile examination centers (MECs). In our analysis, participants aged 20 years and older with complete data about urinary flame-retardant metabolites and rheumatoid arthritis were included for analysis. Finally, we have 5490 participants: 319 with rheumatoid arthritis and 5171 without rheumatoid arthritis. (Fig. 1)
Flow chart of participant selection
2.2. Assessment of RA
The diagnosis of RA status was assessed by self-report in a face-to-face interview during a questionnaire examination regarding the health condition. Arthritis-related questions in the medication condition questionnaire were completed before a physical examination using the computer-assisted personal interviewing (CAPI) system. At first, participants were asked, “Has a doctor or other health professional ever told you that you have arthritis?” stored under MCQ160A, and If the answer was “yes,” a follow-up question was asked: “Which type of arthritis was it? “Stored as MCQ195. Participants were classified as RA or Non-RA according to their answers to the second question.
2.3. Detection of OPEs and creatinine in urine
The urine specimens of enrolled participants were collected, aliquoted, and shipped to the National Center for Environmental Health (NCEH) according to the analytical guidelines. Due to the target analytes being diverse in different NHANES cycles, we selected five consecutively detected flame retardants for statistical analysis: diphenyl phosphate (DPhP), dibutyl phosphate (DBuP), bis(1,3-dichloro-2-propyl) phosphate (BDCIPP), bis(1-chloro-2-propyl) phosphate (BCIPP), bis(2-chloroethyl) phosphate (BCEP), and tetrabro-mobisphenol A (TBBA). Briefly, the analytical method is based on enzymatic hydrolysis of urinary conjugates of target analytes, automated offline solid phase extraction, and isotope dilution high-performance liquid chromatography-tandem mass spectrometry detection. The lower detection limits (LLOD) for urinary BCIPP, BDCPP, DBuP, and DPhP were 0.1 µg/L, and the LLOD for urinary TBBA was 0.05 µg/L.The creatinine in the urine is determined using an Enzymatic Roche coba 6000 analyzer (Roche Di-agnostics, Indianapolis, IN, USA). Urinary creatinine concentrations were measured to account for dilution-dependent sample variation in biomarker concentrations.Table 2 In our research, we selected chemicals for inclusion based on the metabolite’s detectable frequencies being greater than 50%. The distributions of flame retardant metabolite concentrations were described using detection frequency, GM (95% CI), and selected percentiles.
2.4. Assessment of Covariates
Directed Acyclic Graphs (DAGs) are developed based on existing literature and assumptions regarding plausible causal relationships among the covariates. We used the Daggity web application to establish the minimally sufficient adjustment sets of variables to estimate the association and total effects of OPEs on Rheumatoid Arthritis. Covariates included Age(20-60 years, Greater than 60 years), Gender(Female and Male), Ethnicity(Mexican American, Non-Hispanic White, Non-Hispanic Black, Non-Hispanic Asian, and other race), Education(Under high school, Above high school, High school or equivalent,), Marital status(Previously married, Married, and Never Married), the ratio of family income to poverty (PIR) was binned into four levels(Less than 1, 1-1.99, 2-3.99, Greater than or equal to 4), Alcohol status(Yes, No), Smoking(Every day, Some days, Not at all), body mass index (BMI) (kg/m2), Urinary Creatinine(mg/dL), Activity level(Vigorous Activity, Moderate Activity, None), Citizenship(Citizen, Non-citizen), Diabetes(Yes, No). Missing data for the covariates were coded as a missing indicator category for categorical variables (citizenship (n = 1), education levels (n = 7), marital status (n = 1), smoking (n = 3196), and alcohol status (n = 422)) and were imputed with the median for the continuous variables (BMI (n = 48), family income (n = 504) and urinary creatinine (n = 1).
2.5. Statistical Analysis
Continuous variables in covariates were described as means ± standard deviation and analyzed using the Mann-Whitney U test. Univariate analyses were used to assess exposure and covariate distributions. Categorical variables were expressed as numbers (percent) and compared by the Chi-square test (X2) or Fisher exact test. All OPEs exposure variables were natural-log-transformed due to a long right tail before running a multivariate logistics regression. Urinary OPEs metabolites were stratified based on quantiles (Q1: < 25th percentile, Q2: ≥ 25th to 50th percentile, Q3: ≥ 50th to 75th percentile, Q4: ≥ 75th percentile). For multivariate logistic regression, we used sampling weights to account for the complex NHANES survey design. Firstly, we fitted a crude model that assessed the association between each urinary metabolite of OPEs and RA status by comparing the second, third, and fourth quantiles to the first metabolite concentration using multivariate logistic regression. After that, two separate models were conducted for logistic regression: model 1 was adjusted for age and gender, and model 2 was adjusted for ethnicity, education, marital status, citizenship, PIR, alcohol status, smoking, activity level, BMI, diabetes, and urinary creatinine. Stratification analysis was also performed to evaluate the relationship between RA and OPEs in the US population by subgrouping the age and gender respectively based on the fact that incidence of RA majorly varies in different age group and gender. Additionally, Pearson correlations were performed on urinary OPEs to calculate the degree of correlation among them. We also assessed the nonlinear relationship between RA and OPEs with smooth curve fittings using generalized additive models (GAM). GAM is a flexible and smooth technique that captures the non-linearities in the data and helps us to fit nonlinear models. We implemented GAMs in R using the “gam” package. They are generalized versions of the linear model in which the predictors Xi depend either linearly or non-linearly on smooth non-linear functions like splines, polynomials, or step functions. Mathematical representation for Generalized additive models the regression equation:
where the functions f1, f2, f3 … fp are different non-linear functions on variables Xp. Variance inflation factors (VIF) and Spearman’s correlation analysis were employed to examine collinearity among variables for each model, respectively. A VIF value below 5 was considered indicative of the absence of multicollinearity [20]. Although GAM models are capable of offering a flexible and detailed exploratory analysis of the data by helping to identify complex and nonlinear relationships between OPEs exposures and RA risk without mandating strict assumptions about the shape of the dose-response curve. Restricted cubic spline (RCS) allows for precise modeling of the dose-response relationship, providing a detailed visualization of how variations in OPEs exposure levels correlate with changes in RA risk. This approach facilitates the identification of possible thresholds and critical exposure levels. We selected three knots for each OPE as the default and used the Wald test to identify nonlinear dose-response relationships. The logistic model with restricted cubic splines can be written as:
Where h1(V) and h2(V) are the spline basis functions defined by the chosen knots. Adjusted factors were the same as fully adjusted models. RCS analysis was performed by rms R package (version 6.7)
3. Results
3.1. Characteristics of the study population
This study included a total of 5490 adults, distinguished on the basis of RA status. Participants with rheumatoid arthritis (RA) (n = 319) and those without RA (Non-RA) (n = 5171). The mean age of participants was 48.35 years, with 50.9% being female. Participants with RA were of older age, predominantly female, non-Hispanic Black, and with higher BMI as compared to non-RA. We observed a higher prevalence of diabetes and less physical activity (P < 0.05). No significant differences were observed in urinary creatinine and alcohol consumption (P > 0.05).
3.2. Measurement of urinary metabolites of OPEs
The detection frequency, mean concentration, geometric mean concentration, and distribution of OPES are shown in Table 2.TBBA was excluded in subsequent analyses because 99.78% samples were below the LLOD. The detection frequency of the remaining flame retardants was above 50%, including DPHP (95.88%), BDCPP (93.75%), BCEP (82.17%), DBUP (51.07%) and BCPP (51.93%). OPEs exposure patterns varied largely across population subgroups. Higher levels of DPHP and DBUP were observed in females (P<0.05), whereas BDCPP and BCEP levels were higher in males (P<0.05). We observed a significant positive correlation among each pair of OPEs (P< 0.001), with the highest correlation between DPHP and BDCPP (coefficients = 0.53).
3.3. Association between urinary metabolites of OPEs and RA
Based on a sample size of 5,490 observations, the correlation analysis of five urinary metabolites, DPHP, BDCPP, BCPP, BCEP, and DBUP, revealed multiple significant relationships. BDCPP (r = 0.25), BCPP (r = 0.13), BCEP (r = 0.16), and DBUP (r = 0.10) showed positive associations with DPHP. Positive correlations were found between BDCPP and BCPP (r = 0.15), BCEP (r = 0.20), and DBUP (r = 0.02). The BCPP and BCEP showed the strongest positive association (r = 0.28), indicating a moderate link. Except for the correlation between DBUP and BDCPP, which was not significant after correction (p = 0.07), all correlations were statistically significant (p < 0.05). These results suggest that there may be shared sources or comparable metabolic pathways between higher levels of one metabolite and higher levels of another. In model 1, adjusted for age and gender, the relationship between different quartiles of the flame-retardant metabolite DPHP and the risk of RA was examined. Compared to Q1, participants with the increasing quantile of DPHP (Q3, OR: 1.56, 95% CI: 0.95-2.56), DBUP (Q2, OR: 1.52, 95% CI: 1.03-2.24; Q3, OR: 1.45, 95% CI: 0.96-2.21; Q4, OR: 1.60, 95% CI: 1.00-2.57) showed significantly higher risk of RA compared to Q1 and BDCPP (Q3, OR: 0.47, 95% CI: 0.29-0.76) showed lower risk of RA as compared to Q1. In model 2, which was adjusted for age category, gender, BMI, ethnicity, citizenship, education, marital status, Poverty income ratio, smoking, diabetes, alcohol, activity level, and urinary creatinine (log-transformed), BDCPP showed a lower risk of RA compared to (Q4, OR: 0.55, 95% CI: 0.33-0.91) (all p for trend<0.05) compared to Q1. Our methodology’s simultaneous use of GAMs and RCS ensures an extensive examination of the association between OPEs exposure and RA. The initial exploratory observations were provided by GAMs, which highlighted plausible nonlinear patterns. The dose-response relationship was investigated in detail and with specificity by RCS, significantly strengthened the findings’ accuracy and interoperability.
We utilized the generalized additive model (GAM) to explore the nonlinear relationship between RA and OPEs. The smooth function for BDCPP shows a slight non-linear pattern with a p-value of 0.711 estimated edf; 1.5, with a minor upward trend at higher levels. The smooth function for BCEP shows a notable non-linear pattern, with an initial decrease followed by an increase having a p-value of 0.53 and edf; 2.44. Variance inflation factors were calculated for gam utilizing vif.gam() function from mgcv.helper package to examine collinearity among variables. The plot from the RCS revealed distinct characteristics between OPEs and RA. DPHP and DBUP showed a positive association with increasing log odds of RA at higher concentrations. In contrast, BCEP showed a U-shaped relationship, suggesting a higher incidence of RA at both lower and higher concentrations with a huge dip in the center.
3.3.1. Subgroup analysis
The covariates adjusted in the subgroup analysis were age.cat, gender, BMI, ethnicity, citizenship, education, marital status, poverty income ratio, smoking, diabetes, alcohol, activity level, urinary creatinine (log-transformed), and the results from the subgroup analysis showed the association between flame retardants quantiles and RA was mainly present in females and in participants aged more than 60 years. Female had a growing risk of RA in the population with increasing quantiles of DPHP (Q2, OR: 2.52, 95% CI: 1.32-4.81; Q3, OR: 3.07, 95% CI: 1.37-6.88; Q4, OR: 4.04, 95% CI: 1.42-11.46), BCEP (Q3, OR: 2.67, 95% CI: 1.03-6.95), DBUP (Q4, OR: 2.92, 95% CI: 0.94-9.07) compared to Q1. An increased risk of RA was observed in participants aged 20-60 years and greater than 60 years. For DPHP (Q4, OR: 3.50, 95% CI: 0.96-12.70), BCPP (Q2, OR: 2.61, 95% CI: 1.02-6.68; Q3, OR: 3.48, 95% CI: 0.96-12.67; Q4, OR: 11.48, 95% CI: 1.38-5.69), DBUP (Q4, OR: 5.89, 95% CI: 1.03-3.68). The growing risk of RA in males was insignificant for most metabolites and participants over 60 years.
4. Discussion
To the best of our knowledge, this is the first large-sample study, a population-based study to investigate the association between exposure to OPEs and rheumatoid arthritis. From the association results, we observed that exposure to OPEs (DPHP and DBUP) was associated with a 1.45 – 1.60-fold increased prevalence of RA, whereas BDCPP appeared to be protective, reducing the risk of RA among participants between 20-60 years. High levels of OPEs were significantly associated with the elevated prevalence in females and participants between the ages of 20-60 years. No significant association was observed between BCPP and BCEP. We also fitted restricted cubic splines and observed a positive correlation between DPHP and DBUP with RA. Rheumatoid arthritis is a chronic autoimmune disease characterized by inflammation, pain, and joint stiffness [21]. Several studies have suggested a potential association between exposure to organophosphate esters and the development or exacerbation of rheumatoid arthritis. One study conducted on a population of agricultural workers found that those with occupational exposure to OPEs had a significantly higher risk of developing rheumatoid arthritis compared to those without exposure [22]. Another study conducted on a population of Gulf War veterans found a positive correlation between exposure to OPEs during the war and an increased risk of developing rheumatoid arthritis later in life [23]. These findings suggest that exposure to OPEs may play a role in the pathogenesis of rheumatoid arthritis [24]. The mechanism by which OPEs may contribute to the development of rheumatoid arthritis is still unclear. Some hypotheses suggest that OPEs may trigger an immune response, producing autoantibodies and subsequent joint inflammation [25]. However, it is important to note that these studies have limitations, and further research is needed to fully understand the association between OPEs exposure and rheumatoid arthritis. Additionally, there are challenges in studying the relationship between exposure to OPEs and rheumatoid arthritis due to the complex nature of both the disease and the physiochemical activity of OPEs.
Dose-response relationships between urinary OPEs metabolites concentrations and Rheumatoid Arthritis. Adjusted for age.cat, gender, BMI, ethnicity, citizenship, education, marital status, poverty income ratio, smoking, diabetes, alcohol, activity level, and urinary creatinine (log-transformed). The red rugs represent the distribution of concentration.
Associations between exposure to OPEs and RA. Blue bands represent the 95% confidence interval from the fit. The covariates used in the modeling are age.cat, gender, BMI, ethnicity, citizenship, education, marital status, Poverty income ratio, smoking, diabetes, alcohol, activity level, urinary creatinine (log-transformed).
We observed the curves for DPHP and DPHP, which suggest a positive association between higher exposure levels and an increased risk of RA. BCEP had a U-shaped relationship, with an increased risk of RA at lower and higher exposure; however, the associations were insignificant. The curves for BCPP and BDCPP show an initial increase in RA risk with exposure to OPEs and then a decrease at higher levels. The subgroup gender stratification suggests that exposure to OPEs is higher in females than males. The observed association of OPEs and RA herein is biologically plausible as nail polish is considered a significant source of short-term DPHP exposure and chronic exposure for frequent users or those occupationally exposed [26]. Studies have found that exposure to OPEs can disrupt both innate and adaptive immune systems, causing delayed or inappropriate immune responses to pathogens [15]. The proper equilibrium of extracellular immune ligands, which are essential for starting and controlling immune responses, is disrupted by OPEs.
Epidemiological studies have confirmed that exposure to certain OPEs reduces the adhesion capacity of macrophages, impairing their ability to respond to inflammatory signals and pathogens. According to a study on human cell lines, flame retardant exposure increased the release of pro-inflammatory cytokines like IL-1β, IL-6, and IL-17 and increased the expression of the membrane protein E-cadherin, which impairs LPS-induced inflammatory responses in macrophages [27]. Alturaiki, in a recent study, found higher levels of rheumatoid factor (RF), C-reactive protein (CRP), and erythrocyte sedimentation rate (ESR) in RA patients [28]. He also found serum levels of IL-1β, IL-6, IL-8, and CCL5, but not TNF-α, significantly increased in those with RA[29]. It is plausible that the production of these disease-relevant pro-inflammatory cytokines is intricately linked to the development of rheumatoid arthritis (RA) [30]. RA is well-documented for causing progressive destruction of synovial joints, primarily through irreversible damage to articular cartilage and bone. During RA inflammation, cytokines such as IL-1β, IL-6, and TNF-α play pivotal roles by activating the JAK/STAT pathway, which subsequently increases the expression of receptor activator of nuclear factor kappa-β ligand (RANKL) [31]. The activated JAK/STAT3 pathway also suppresses osteoprotegerin expression, mediated by the Wnt signaling pathway. RANKL, by binding to its receptor RANK on the surface of osteoclast precursors, activates the NF-κB pathway, thereby promoting osteoclast differentiation and bone resorption [15]. This cascade of molecular events not only leads to enhanced bone degradation but also perpetuates the inflammatory environment within the joint. Additionally, the chronic inflammation and cytokine milieu in RA can induce synovial fibroblasts to produce matrix metalloproteinases (MMPs), which further degrade the extracellular matrix of the cartilage. This creates a vicious cycle of inflammation and tissue destruction, exacerbating joint damage. Understanding these pathways highlights potential therapeutic targets for RA, such as inhibitors of the JAK/STAT path-way, RANKL-RANK interactions, and MMP activity, aiming to mitigate inflammation and prevent joint destruction [31, 32]. Furthermore, some studies have shown that exposure to flame retardants will potently induce the production of RA-relevant cytokines TNF-α, IL-6, and IL-1β [33].
5. Conclusion
In conclusion, our results indicate that exposure to flame retardants (DPHP and DBUP) is associated with an increased prevalence of RA, however, inverse association between urinary BDCPP concentration was observed. These associations appear to be more pronounced in female and older people. Given the cross-sectional design of the current study, these results should be interpreted with caution, and further investigations are needed to confirm our findings.
Appendix A. Sample Appendix Section
Not Applicable.
Data Availability
All data produced in the present study are available upon reasonable request to the authors.