Association between Organophosphate Ester Exposure and Insulin Resistance with Glycometabolic Disorders among Older Chinese Adults 60-69 Years of Age: Evidence from the China BAPE Study

Abstract

Background: Organophosphate esters (OPEs) are common endocrine-disrupting chemicals, and OPE exposure may be associated with type 2 diabetes (T2D). However, greater knowledge regarding the biomolecular intermediators underlying the impact of OPEs on T2D in humans are needed to understand biological etiology. Objectives: We explored the associations between OPE exposure and glycometabolic markers among older Chinese adults 60-69 years of age to elucidate the underlying mechanisms using a multi-omics approach. Methods: This was a longitudinal panel study comprising 76 healthy participants 60-69 years of age who lived in Jinan city of northern China. The study was conducted once every month for 5 months, from September 2018 to January 2019. We measured a total of 17 OPEs in the blood, 11 OPE metabolites in urine, and 4 glycometabolic markers (fasting plasma glucose, glycated serum protein, fasting insulin, and homeostatic model assessment for insulin re-sistance). The blood transcriptome and serum/urine metabolome were also evaluated. The associations between individual OPEs and glycometabolic markers were explored. An adverse outcome pathway (AOP) was established to determine the biomolecules mediating the associations. Results: Exposure to five OPEs and OPE metabolites (trimethylolpropane phosphate, triphenyl phosphate, tri-iso-butyl phosphate, dibutyl phosphate, and diphenyl phosphate) was associated with increased levels of glycometabolic markers. The mixture effect analysis further indicated the adverse effect of OPE mixtures. Multi-omics analyses revealed that the endogenous changes in the transcriptional and metabolic levels were associated with OPE exposure. The putative AOPs model suggested that triggers of molecular initiation events (e.g., insulin receptor and glucose transporter type 4) with subsequent key events, including disruptions in signal transduction pathways (e.g., phosphatidylinositol 3-kinase/protein kinase B and insulin secretion signaling) and biological functions (glucose uptake and insulin secretion), may constitute the diabetogenic effects of OPEs. Discussion: OPEs are associated with the elevated risk of T2D among older Chinese adults 60-69 years of age. Implementing OPE exposure reduc-tion strategies may help reduce the T2D burden among these individuals, if the relationship is causal.

Headnote

Background: Organophosphate esters (OPEs) are common endocrine-disrupting chemicals, and OPE exposure may be associated with type 2 diabetes (T2D). However, greater knowledge regarding the biomolecular intermediators underlying the impact of OPEs on T2D in humans are needed to understand biological etiology.

Objectives: We explored the associations between OPE exposure and glycometabolic markers among older Chinese adults 60-69 years of age to elucidate the underlying mechanisms using a multi-omics approach.

Methods: This was a longitudinal panel study comprising 76 healthy participants 60-69 years of age who lived in Jinan city of northern China. The study was conducted once every month for 5 months, from September 2018 to January 2019. We measured a total of 17 OPEs in the blood, 11 OPE metabolites in urine, and 4 glycometabolic markers (fasting plasma glucose, glycated serum protein, fasting insulin, and homeostatic model assessment for insulin re-sistance). The blood transcriptome and serum/urine metabolome were also evaluated. The associations between individual OPEs and glycometabolic markers were explored. An adverse outcome pathway (AOP) was established to determine the biomolecules mediating the associations.

Results: Exposure to five OPEs and OPE metabolites (trimethylolpropane phosphate, triphenyl phosphate, tri-iso-butyl phosphate, dibutyl phosphate, and diphenyl phosphate) was associated with increased levels of glycometabolic markers. The mixture effect analysis further indicated the adverse effect of OPE mixtures. Multi-omics analyses revealed that the endogenous changes in the transcriptional and metabolic levels were associated with OPE exposure. The putative AOPs model suggested that triggers of molecular initiation events (e.g., insulin receptor and glucose transporter type 4) with subsequent key events, including disruptions in signal transduction pathways (e.g., phosphatidylinositol 3-kinase/protein kinase B and insulin secretion signaling) and biological functions (glucose uptake and insulin secretion), may constitute the diabetogenic effects of OPEs.

Discussion: OPEs are associated with the elevated risk of T2D among older Chinese adults 60-69 years of age. Implementing OPE exposure reduc-tion strategies may help reduce the T2D burden among these individuals, if the relationship is causal. https://doi.org/10.1289/EHP11896

Introduction

Diabetes mellitus, mainly type 2 diabetes (T2D), is one of the four most prevalent noncommunicable diseases.1,2 In 2021, 536.6 (10.5%) million adults 20-79 years of age worldwide were esti-mated to have diabetes, which resulted in 6.7 million deaths and health expenditures totaling USD $966 billion.3 Diabetes impacts aging societies, affecting an estimated 19.3% of adults 65-99 years of age worldwide.4 Most of the diabetes burden falls on low- and medium-income countries, such as China and India.5 In the past decade, the number of diabetic patients in China increased from 90 (age-adjusted prevalence: 8.8%) to 140 million (10.6%),3 and the number of older individuals (>65 years of age) with diabetes was 35.5 million, ranking first globally and accounting for 25% of the elderly diabetes patients worldwide.4

The risk of T2D is multifactorial and determined by an interplay of genetic and environmental factors.6 Family history of diabetes, older age, unhealthy diet, physical inactivity, smoking, and over-weight and obesity are associated with an increased risk of T2D; of these, overweight and obesity was found to be the strongest risk factor.7 Of note, the age-standardized prevalence of general obesity among Chinese adults, is 14.0%,8 considerably lower than that among American adults (43.3%).9 However, the estimated age-standardized prevalence of diabetes is almost the same among Chinese and American adults (12.4% vs. 14.3%).10,11 Genetic factors, a thrifty phe-notype,12-14 and unique dietary patterns (e.g., low intake of whole grains but high intake of refined grains)15 are suspected to be responsi-ble for the increased risk of T2D in the Chinese population. In addi-tion, since industrialization advancements and immense growth took place after the 1980s, chemical pollution has been the leading envi-ronmental cause of diseases and deaths in China.16 Emerging pollu-tants are ever-evolving and are widely used in the industry and for daily necessities worldwide, and their production is gradually shifting to China owing to the increasingly strict regulations in Europe and America. Extensive evidence suggests that disproportionate exposure to environmental pollutants may be an underappreciated contributor to disparities in the incidence of T2D.17

Long-term exposure to environmental pollutants combined with frailty and a decline of immune function makes the elderly

population particularly susceptible to T2D.18 Cohort studies have revealed that exposure to endocrine-disrupting chemicals (EDCs) disturbed the onset and progression of T2D among older Swedish and Korean adults.19-21 As a substitute for polybrominated di-phenyl ethers, organophosphate esters (OPEs) are commonly used as flame retardants and plasticizers,22 and the global annual con-sumption of OPEs exceeded 680,000 tons in 2015, with 30% of the usage was reported in China.23 As a result, OPEs are ubiquitously detected in Chinese individuals in the microgram to milliliter range.24 Such high OPE exposure levels adversely impact the reproductive, endocrine, cardiovascular, and nervous systems.25,26 In vivo and in vitro studies have reported that OPEs may cause dia-betogenic effects27-29; however, evidence confirming these find-ings at the population level is relatively scarce,30,31 and the mechanisms underlying the indicated association remain largely unknown. Further research and intervention are required to reduce the individual and societal burden of diabetes among the Chinese elderly population. Such studies should incorporate a comprehen-sive clarification of modifiable diabetogenic OPEs and the mecha-nistic standpoints in terms of the etiology of T2D. Furthermore, high-throughput multi-omics (e.g., transcriptomics and metabolo-mics) approaches can be used to fill in the knowledge gaps between exposure to pollutants and the pathophysiology of T2D.32-34

The present study was part of a well-characterized longitudinal multi-omics study, the China Biomarkers of Air Pollutant Exposure (BAPE) Study, conducted among Chinese individuals 60-69 years of age, with monthly longitudinal monitoring conducted over 5 months.35 The study aimed to a) explore the associations between internal OPE exposure and glycometabolic markers, b) identify the key OPEs associated with glycometabolic markers using mixture effect analysis, and c) elucidate the potential biological perturba-tions and mechanisms underlying the associations between key OPEs and glucose homeostasis using multi-omics approaches and an adverse outcome pathways (AOPs) model.

Methods

Study Design and Population

This investigation draws on the China BAPE longitudinal panel study that included a visit each month for 5 months to healthy older Chinese adults 60-69 years of age. In brief, 76 healthy adults 60-69 years of age (50% males, 100% Chinese Han nation-ality) were recruited from within the Jinan megacity, China. Each participant was assessed five times over a 5-month period from 10 September 2018 to 19 January 2019, with a 1-month interval between each assessment. To minimize the potential impact of different dietary patterns, all participants were freely provided three standardized and nutritionally balanced meals (e.g., rice, meat, vegetable, and fruit; in line with local dietary habits) per day for 5 continuous days before bio-sample collection and phys-ical examination. All of the participants were required to finish a questionnaire regarding family information, personal informa-tion, and time-activity patterns prior to the physical examination. During the physical examination, fasting venous blood samples and midstream urine samples were collected from the participants at 0700 hours and immediately stored at - 80 C until further processing. Of the 76 participants, 58 (76%) finished 5 visits; 12 (16%) finished 4 visits; 3 (4%) finished 3 visits; and 3 (4%) finished 2 visits. A total of 353 person-visits were included in our analyses. The study protocol was approved by the ethics commit-tee of the National Institute of Environmental Health (NIEH, Chinese Center for Disease Control and Prevention, No. 201816). Explicit written informed consent was acquired from all the participants. The detailed protocol of the China BAPE Study has been described in a previous publication.35

Characterization of OPE Exposure

Seventeen OPEs were measured in the blood samples (n = 352), and 11 OPE metabolites (m-OPEs) were evaluated in the urine samples (n = 353) (Table S1). Analytical procedures for the measurement of OPEs and the internal concentrations of OPEs were previously reported.36 Urine creatinine (Cr) was measured using a Flex reagent cartridge in a modified kinetic Jaffe assay (model RxL; Dade Behring). Cr-adjusted concentrations of urine m-OPEs were used in further analyses. Concentrations less than the limit of detection (LOD) were imputed using the LOD divided by 2 for each OPE bio-marker. The pairwise correlations of the OPE exposures were calcu-lated using the Spearman coefficient. All participants were self-reported nonsmoking participants; thus, plasma cotinine levels were measured using Hypersil GOLD C18 Selectivity HPLC Columns (ThermoFisher Scientific) interfaced with an LC-Q-Exactive Orbitrap Mass Spectrometer (ThermoFisher Scientific) to assess whether the participants were either actively or passively exposed to tobacco.

Assessment of Glycometabolic Markers

Glycometabolic marker levels, including fasting plasma glucose (FPG) and glycated serum protein [GSP, an indicator of short-to-medium-term (latest 2-3 wk) average blood glucose levels], were measured in all blood samples (n = 353) at Calibra Diagnostics Co. Ltd. using the Cobas 8000 c702 and Cobas 6000 c501 (Roche) modular analyzer series, respectively. Fasting insulin (FINS) levels were determined using a Millipex Human Metabolic Hormone Panel V3 (HMHEMAG-34K-07; Merck) with a flu-orescence detection system (Magpix; Luminex Corporation) and the xPonent 4.2 (Luminex Corporation) and Bio-Plex Manager (version 6.1; Bio-Rad) software, according to the manufacturer's protocol. Homeostatic model assessment for insulin resistance (HOMA-IR) was performed using the fol-lowing formula: [FPG (in millimoles/liter) FINS (in micro units/ milliliter) ] /22.5.37

Transcriptome Analysis

Total RNA was extracted from the leukocytes of venous blood samples (n = 346). In brief, leukocytes were isolated from 5 mL of whole blood. After centrifugation at 3,000 rpm for 5 minutes, the leukocytes aggregated in the medium were extracted and washed two times with phosphate-buffered saline. The leukocytes were then lysed in 1 mL TRIzol Reagent (Invitrogen Corp.). Next, the total RNA was assessed for quality and quantified using a NanoDrop ND-2000 (ThermoFisher Scientific) and Agilent 2100 Bioanalyzer. Library preparation with 100 ng of the extracted RNA was conducted through the TruSeq Stranded Total RNA Library Prep Kit (Illumina, Inc.). Last, the RNA was sequenced using the Illumina Hiseq X Ten System (Novogene). Genes were identified by using HISAT2 software and comparing the data with the human database.38 Gene expression quantification was conducted using the featureCounts read summarization program in the Subread software package.39

Nontargeted Metabolome Assessment

Nontargeted metabolomics analyses were performed in serum (n = 353) and urine (n = 346) samples. Sample processing, quality control and data extraction, metabolite identification, data curation, data quantification, and data normalization were conducted according to previously published protocols.40 In brief, the samples were precipitated with methanol for 2 min under severe vibration (GenoGrinder 2000), and then the metabolites were extracted by centrifugation. The supernatant was then separated into five fractions: Two fractions were analyzed by two different reverse phase/ultra-performance liquid chromatography-mass spectrum (RP/UPLC-MS/MS) methods using positive ion mode electrospray ionization (ESI). The remaining three fractions were analyzed by hydrophilic interaction liquid chromatography (HILIC)/UPLC-MS/MS with negative ion mode ESI, RP/UPLC-MS/MS with negative ion mode ESI, or kept for backup, respec-tively. The samples were put in TurboVap (Zymark) for a short time to remove the organic solvent and then kept at -80 C over-night until subsequent analysis.

The metabolomics analysis was conducted using Waters ACQUITY UPLC and Q-Exactive high resolution/accurate mass spectrometer (Thermo Fisher), which were connected to an Orbitrap mass analyzer running at 35,000 mass resolution and a heated ESI source. The first and second aliquots of sample extracts were analyzed following acidic positive ion conditions that were chromatographically optimized for more hydrophilic and hydro-phobic compounds, respectively. The third aliquot of sample extracts was gradient eluted from a separate dedicated C18 column and analyzed using basic negative ion-optimized conditions. The fourth aliquot of sample extracts was eluted from a HILIC column and analyzed via negative ionization. The mass spectrum analysis used dynamic exclusion that alternated between MS or data-dependent multistage mass spectrometry (MSn) scans with scan range covering 70-1,000 m/z.

After the raw mass spectrum data extraction and peak recog-nition, the metabolites were identified by comparison with an in-ternal library. The library was based on certified standards that included a) m/z, b) retention time, and c) chromatographic data of all molecules. In addition, metabolite identifications were based on the following three criteria: a) retention time within 100 ms tolerance, b) accurate mass match to the internal library (10-ppm mass tolerance), and b) MS/MS fractions based on the ions present in the experimental spectrum compared with the ions present in the library spectrum. To ensure data quality, all peaks were manually checked. Then the area under the curve was used to quantify and all data were normalized before further analysis.

For quality control (QC) purposes, three types of controls were analyzed in concert with the experimental samples: a) process blanks, b) mixed matrix samples (10 uL of each sample), and c) a cocktail of QC standards spiked into samples. All the experimental samples were randomly distributed on the platform and QC samples were also evenly distributed among the injections.

Statistical Analyses

All data were analyzed according to the following pipeline: a) a lin-ear mixed-effects model (LMM) was conducted to explore the associations between OPE exposures and the glycometabolic markers; b) quantile g-computation (qgcomp)41 was used to assess the effects of OPE mixtures on the glycometabolic markers and to identify the most important OPEs with relative positive weights >10% within the OPE mixtures; c) the associations between each key OPE and each biomolecule (transcripts and metabolites obtained using multi-omics profiling) were examined using LMM, and the biomolecules mediating the impact of each key OPE on a specific glycometabolic marker were determined using the causal inference test (CIT)42; and d) Integrative Pathway Analysis (IPA) of the biomolecular intermediators was conducted to investigate the underlying biological mechanisms43 (Figure 1). The detailed processes are described in the following sections.

LMM

The measurement values of certain OPE exposures and glycometabolic markers were logarithm (base 10)-, square root-, power one-third-, square-, or cubic-transformed to approach Gaussian distri-bution (Table S1). The maximum missing rate for certain exposure variables (blood OPEs) was 0.28% owing to the runout of one blood sample. After imputing the missing data for exposures using a chained equation (mice package with the predictive mean match-ing method),44 the exposures were standardized through z-score normalization to denote a change of 1 standard deviation (SD) in the glycometabolic marker values to facilitate their comparisons. A LMM with participant-specific intercepts and "unstructured" co-variance structure was used to independently assess the associations between exposure to each OPE and each glycometabolic marker. Each main model was adjusted for the set of predefined adjustment factors according to previous studies45,46: age (continu-ous), sex (female/male), body mass index [BMI; continuous; cal-culated as weight divided by height squared (in kilograms per meter squared)], education level (below primary school, primary school, junior school, senior high school, and university), financial income [continuous; annual household income (in 10,000 CNY)], plasma cotinine concentration (continuous), and other diet (the total frequency of extra food consumption in addition to the pro-vided standardized meals during the 3-d investigation). Among them, information on age, sex, education level, and financial income was assessed at the first visit for all of the participants, whereas other time-varying factors (BMI, plasma cotinine concentration, and other diet) were assessed at all five visits. Other diet information was recorded using a daily time-activity questionnaire at each visit, and the summarized results are shown in Table S2. Stratification analysis of the associations between OPE exposures and glycometabolic markers were conducted by sex. Multiple hypothesis testing-corrected p-values were obtained by calculating the false discovery rate (FDR) and the estimated proportion of false discoveries made vs. the number of total discoveries made at a given significance level (oc).47 In multiple testing corrections, FDR was statistically significant at <5%. b estimates and standard errors from the models were converted to percentage change values with the 95% confidence intervals (95% CIs) associated with 1-SD increases in individual OPE concentrations.

qgcomp

The OPEs that were found to be significantly associated with each glycometabolic marker after multiple testing using LMM were packed as a chemical mixture. The qgcomp method was performed to assess the effects of the OPE mixtures on the glycometabolic markers based on parameter inference via the qgcomp package in R. This approach combined the inferential simplicity of weighted quantile sum regression (WQS) with the flexibility of g-computation.41 The advantage of qgcomp was that exposures could interact with outcomes in any directions. Gaussian distributions were specified as link functions, and parameter q was set to four in the linear model. Five hundred bootstrap iterations were performed to calculate the 95% CI for each mixture. With the exception of predefined adjustment factors, visit time was included as an adjusting covariate for qgcomp analysis. OPEs with positive weights of >10% were identified as key components.

CIT

To identify the biomolecular intermediates, the causal relationship inference of OPEs-biomolecular intermediators-glycometabolic markers was assessed using the CIT.42 In brief, each OPE-metabo-lite/transcript-glycometabolic marker relationship was individu-ally analyzed to classify the components as consequential to/ mediated by/independent of expression of gene or metabolite. In this study, causal inference was defined by the following four criteria: a) the OPEs and glycometabolic markers were significantly correlated, b) the OPEs were significantly associated with the expression of biomolecules after adjusting for the glycometa-bolic markers, c) the expression of biomolecules was signifi-cantly associated with the glycometabolic markers after adjusting for OPEs, and d) the OPEs were independent of the glycometabolic markers after adjusting for the expression of biomolecules. In addition to predefined adjustment factors, visit time (month of sample collection) was included as an adjust-ment covariate for the CIT model. To summarize the p-values for the whole CIT, the intersection-union test framework was used as the maximum p-value of the four test components.42

IPA

IPA canonical pathway analysis of all the significant biomolecular intermediators (both genes and metabolites) of the causal relation-ship between each OPE and glycometabolic markers was per-formed using IPA software (version 68752261; Qiagen; https:// digitalinsights.qiagen.com/IPA). The corresponding biological pathways of these transcripts and metabolites were identified for each OPE independently. Canonical pathways with a p < 0:05 (Fischer's exact test) were regarded as statistically significant. IPA helps to uncover new biological insights and interpretations by combining an enormous schemata of existing knowledge from the literature with a massive collection of gene and metabolite expression measurements.

Sensitivity Analysis

Sensitivity analysis was performed to assess various models by controlling for one of the following covariates in the main model: month of sample collection, tea consumption (number of cups over 3 d), and frequency of alcohol consumption (over 3 d). A two-sided p < 0:05 in sensitivity analysis indicated sta-tistical significance. All statistical evaluations were conducted using R (R Development Core Team) with the lme4 and qgcomp packages.

Results

Overview of OPE Exposure and Glycometabolic Markers

Participant characteristics are presented in Table 1. Of all partici-pants, 38 (50%) were women, and the mean SD participant age was 64:5 4:5 y. The concentrations of the 17 blood OPEs and 11 urinary OPE metabolites over five longitudinal visits are summar-ized in Figure 2A and Table S1 (see also Excel Table S1). The con-centrations of tributyl phosphate (TBP), triphenyl phosphate (TPHP), and tri(1-chloro-2-propyl) phosphate (TCPP) were rela-tively high in the blood samples, whereas those of bis(2-chloroethyl) phosphate (BCEP), di(2-ethylhexyl) phosphate (DEHP), and bis(1,3-dichloro-2-propyl) phosphate (BDCPP) were relatively high in the urine samples. The pairwise Spearman corre-lations of the OPE exposures are shown in Figure S1. Meanwhile, four glycometabolic markers (FPG, GSP, FINS, and HOMA-IR) reflecting the statusof insulin resistance and glycometabolic home-ostasis are shownin Figure 2B (see also Excel Table S2).

Association between OPE Exposure and Glycometabolic Markers

The percentage change in estimates and 95% CIs of the glycome-tabolic markers for a 1-SD increase in the concentration of each OPE are shown in Figure 3A and Excel Table S3. After multiple testing corrections, 6, 7, 3, and 9 OPEs were found to be signifi-cantly associated with FPG [TPHP, trimethylolpropane phosphate (TMPP), 2-ethylhexyl diphenyl phosphate (EHDPP), diphenyl phosphate (DPHP), dibutyl phosphate (DBP), and BCEP], GSP [TPHP, TMPP, EHDPP, hydroxyphenyl 2-ethylhexyl diphenyl phosphate (EHDPP-OH), DPHP, DBP, and BCEP], FINS (TPHP, DPHP, and DBP), and HOMA-IR [TPHP, tri-n-butyl phosphate (TnBP), TMPP, tri-iso-butyl phosphate (TiBP), TBP, EHDPP, DPHP, DBP, and BCEP], respectively. Of these OPEs, blood TPHP and urinary DBP, and DPHP showed an undesirable posi-tive influence on all four glycometabolic markers. The results of the sensitivity analysis were consistent with the findings of the main model (Figure S2 and Excel Table S4). A stratification analy-sis of effect modification of sex was conducted (Figure S3 and Excel Table S5), and the correlation directions and significance between OPE exposure and glycometabolic marker changes in different sex were generally consistent.

Identification of the Key OPEs

After multiple testing corrections, OPEs that were associated with the four glycometabolic markers were grouped as four different mixtures for subsequent qgcomp analyses. Figure 3B shows the log-transformed concentrations of four OPE mixtures that were posi-tively associated with the z-scores of FPG (p< 0:001), GSP (p = 0:001), FINS (p= 0:002), and HOMA-IR (p= 0:005), respec-tively (see also Table S3). Each quartile increment in the mixture concentration was associated with elevated z-scores of 0.06 (95% CI: 0.03, 0.08), 0.04 (95% CI: 0.02, 0.05), 0.06 (95% CI: 0.02, 0.09), and 0.12 (95% CI: 0.04, 0.20) for FPG, GSP, FINS, and HOMA-IR, respectively. In Figure 3C and Table S4, the weight of each OPE reflects the contribution of the correlated compo-nents to the overall mixture effect. Within the mixtures, TMPP had the highest weight (54%) for elevated FPG, followed by DBP (24%) and TPHP (14%), whereas exposure to BCEP and DPHP contributed minimally (<10%) to the overall mixture positive effect. TMPP contributed to 46% of the overall positive mixture effect of elevated GSP, followed by DBP and TPHP (26% and 19%, respectively). TPHP and TiBP had the highest weights of contribution to FINS and HOMA-IR for positive weight, at 44% and 29%, respectively. For negative weight, EHDPP had the highest weights for FPG, GSP, and HOMA-IR with 100%, 97%, and 57%, respectively. Ultimately, five OPEs

(TMPP, TPHP, TiBP, DBP, and DPHP) with positive weights of >10% were considered the key OPEs and were suspected to be responsible for the increased risks of glycometabolic disorder and therefore were included in further analyses.

Longitudinal Correlations of OPEs-Metabolites

The numbers of endogenous metabolites, based on metabolo-mics classifications in the serum and urine samples, are shown in Figure 4A. The numbers and directions of the associations between the five key OPEs and the respective serum and urine metabolomes are presented in Figure 4B,C. A total of 1,354 associations were found (710 and 644 for the serum and urine metabolomes, respectively; FDR <0:05; Excel Tables S6 and S7). The network diagram shows the exposure-metabolite pair-wise associations with FDR <0:001 for the serum and urine metabolomes (Figure 4D,E; Excel Tables S6 and S7; raw serum and urine metabolome data in Excel Tables S8 and S9). Specifically, for the serum metabolome, the top three OPEs of high connective degree included TPHP, TMPP, and DBP. The top three metabolite classes with positive associations were peptides (20%), cofactors and vitamins (12%), and energy (12%), whereas those with negative associations were peptides (18%), carbohydrates (15%), and nucleotides (12%) (Figure 4F; Table S5). For the urine metabolome, DBP, DPHP, and TPHP were the three OPEs with the highest connective degree and were found to have mostly negative correlations. The top three metabolite classes with positive associations were carbohy-drates (10%), secondary metabolites (10%), and cofactors and vitamins (8%), whereas those with negative associations were nucleotides (17%), global metabolites (12%), and cofactors and vitamins (11%) (Figure 4G; Table S6).

Longitudinal Correlations of OPEs-Transcripts

The associations between five key OPEs and the blood transcriptome were also explored, and the numbers and directions of the associations are displayed using volcano plots (Figure 5A; Excel Table S10). A total of 16,585 associations (8,673 positive and 7,912 negative) were significant. The network diagram shows the OPE-transcript pairwise associations of the top 100 transcripts selected for absolute effect value, with an FDR of <0:001 for each OPE (Figure 5B; Excel Table S10). Specifically, the top three OPEs correlated with gene expressions of high connective degree were TPHP (3,875 positive and 3,476 negative associations), DBP (3,056 positive and 2,653 negative associations), and TMPP (1,741 positive and 1,783 negative associations) (Figure 5C).

Biomolecules Mediating the Effects of OPEs on Glycometabolic Markers

A tripartite network plot of OPEs-metabolites-glycometabolic markers constructed using all of the significant metabolites and transcripts as potential biomolecular intermediators is presented in Figure 6A (see also Excel Table S11). A total of 106 serum and 29 urinary metabolites were inferred as biomolecular intermediators, with an FDR of <0:2. Specifically, 87, 60, 31, and 10 metabolites were inferred as biomolecular intermediators of the effects of TMPP, TPHP, DPHP, and DBP on glycometabolic markers, respectively (the number of common and specific metabolites of each OPE are shown in Figure S4A). In addition, 1,651 transcripts were inferred as biomolecular intermediators, with an FDR of <0:1, and a network plot of OPEs-transcript- glycometabolic markers is presented in Figure 5D (top 150 genes selected for the absolute effect value for each OPE; see also Excel Table S12). Specifically, 1,248, 649, and 116 transcripts were inferred as biomolecular intermediators of the effects of TMPP, TPHP, and DBP on glycometabolic markers, respectively (the number of common and specific transcripts of each OPE are shown in Figure S4B).

IPA for Biomolecular Intermediators

Biomolecular intermediators of the diabetogenic effects of TPHP, TMPP, DBP, and DPHP were included in the IPA analysis. Representative canonical pathways (p< 0:05) shared by at least two OPEs are presented in Figure 6B and Excel Table S13. Among these, several pathways [e.g., peroxisome proliferator-activated receptor alpha (PPARa)/retinoid X receptor alpha (RXRa) activa-tion pathway, leptin signaling in obesity, G-protein alpha-s (Gas) signaling, myelocytomatosis viral oncogene (MYC)-mediated apo-ptosis signaling, and pyroptosis signaling] simultaneously appeared in the analysis results of at least three key OPEs. Core perturbed genes and metabolites within the abovementioned representative pathways were integrated to propose putative AOPs for impaired glucose homeostasis after OPE exposure (Figure 7A). These AOPs

began with possible molecular initiation events (MIEs): activation of first apoptosis signal ligand (FasL) and protease-activated receptor (PAR), as well asinhibitionof tumor necrosis factor (TNF), TNF receptor (TNFR), insulin receptor (IR), lysophosphatidic acid receptor (LPAR), and glucose transporter type 4 (GLUT4). Alterations of these ligands and membrane receptors caused a series of key events (KEs) at the molecular/cellular levels, such as KE-1, those involved in aberrant expression of kinases [phosphatidylinositol 3-kinase (PI3K), protein kinase B (AKT), and glycogen synthase kinase-3 (GSK3)], the apoptosis/autophagy regulator Bcl2-antagonist of cell death protein (BAD), the transcription factor nuclear factor kappa B (NF-jB), metabolic enzymes (e.g., citrate synthase and succinate dehydrogenase), and tricarboxylic acid (TCA) cycle-related metabolites (e.g., citrate and succinate). The downstream KEs at the organ/system levels (KE-2; e.g., those involved in glucose uptake, glycogen synthesis, insulin secretion, energy metabolism, oxidative stress, and inflammation) were subsequently altered and were responsible for perturbations in glucose homeostasis eventually leading to T2D. In general, the up-regulation or down-regulation trends of most MIEs and KEs were consistent among the key OPEs, except for certain molecules (e.g., PAR, PPARa, BAD, and GSK3) for urinary DBP and DPHP.

Discussion

OPEs are widely used as a substitute for polybrominated di-phenyl ethers and are ubiquitously detectable in environmental matrices and human bodies, which has caused considerable concern in terms of their adverse effects worldwide. The present study, for the first time, comprehensively profiled the associations between OPE exposure and insulin resistance and glycometabolic homeostasis among healthy older Chinese adults 60-69 years of age and systematically explored the underlying biological mecha-nisms using multi-omics profiling.

Identification of the Key OPEs

In this study, exposure to five key OPEs and m-OPEs (TMPP, TPHP, TiBP, DBP, and DPHP) was found to be associated with elevated levels of glycometabolic markers, independent of the traditional risk factors for T2D (e.g., age, sex, BMI, and diet). To the best of our knowledge, the observed associations were first reported at the population level; moreover, some associations had, to varying extents, a priori credibility based on the previous exper-imental literature (Table S7).28,29,48-62 Specifically, although TMPP, TPHP, and DPHP were widely detected in the Chinese population at concentration ranges of 0:24-1:68 ng=mL,24,63,64 to the bestof our knowledge, our studyisthefirst toreveal the adverse correlations between blood TPHP and its metabolite DPHPin urine and glucose homeostasis within humans. These findings are in line with those of previous experimental studies reporting the diabeto-genic effects of TPHP and DPHP in vivo and in vitro. For instance, exposure to TPHP results in a higher glucose level and HOMA-IR index and inhibits the level of adiponectin, an insulin-sensitizing hormone, in adult male mice49 and pubertal mice.27 Likewise, numerous in vitro studies have confirmed the diabetogenic effects of TPHP and DPHP on 3T3-L1 adipose and HepG2 cells.59-61 However, no previous studies have yet reported the relationships between exposure to TMPP, TiBP, and DBP and abnormal glucose metabolism. Thus, further health hazard assessments through epi-demiological studies or in vivo and in vitro experiments are war-ranted. In addition, it is noteworthy that the dominant pollutants differed between specific glycometabolic markers. For instance, TMPP and TPHP were the leading contributors of elevated FPG and FINS levels, respectively, possibly because of their different diabetogenic potential in triggering MIEs or downstream KEs. Therefore, subsequent CIT and multi-omics analyses were con-ducted to determine the biomolecular intermediators and obtain mechanistic insights into specific exposure-outcome associations.

Determination of Biomolecular Intermediators and Relevant Pathways

In this study, we comprehensively characterized the biomolecular changes in genes and metabolites induced by the key OPEs in healthy older adults for the first time. Some endogenous changes were first unveiled at the population level, including PI3K/AKT signaling [PI3K, AKT, and mitogen-activated protein kinase kinase kinase 8 (MAP3K8) inhibition], IR signaling [IR and growth factor receptor-bound protein 2 (GRB2) inhibition and IR substrate activation], PPARa=RXRa signaling (protein kinase A and PPARa activation, and RXRa and NF-jB inhibition), insulin secretion signaling [eukaryotic initiation factor 2 (EIF2) activa-tion], and TCA cycle processes (citrate synthase, succinate dehy-drogenase, and citrate inhibition). Our results were consistent with

those of a previous in vitro study that reported that PI3K, AKT, and PPARa mRNAs were significantly down-regulated after TPHP treatment of L02 cells.65 Dysregulation in PI3K/AKT signaling, especially AKT activity inhibition, increases insulin resistance among patients with T2D.66 Besides, an in vivo study conducted in mice, showed that TPHP reduces the abundance of citrate and suc-cinate within the TCA cycle52; moreover, the down-regulation of succinate dehydrogenase can impair glucose metabolism in islet cells in T2D.67 Similar to our findings, previous toxicological studies have reported that TPHP could induce cell apoptosis,68,69 autophagy,70 oxidative stress,71 and glucocorticoid regulation dis-ruption72-74 (Table S8).48,60,68,75-93 Alterations of some other pathways (e.g., leptin signaling,94 relaxin signaling,95 and Ga12=13 signaling96) are also closely related to insulin resistance and glucose metabolism, although they have not yet been reported to induce diabetogenic effects upon OPE exposure. In addition, to a certain extent, blood TMPP and TPHP share the most common bio-molecular intermediators and pathways. In contrast, urinary DBP and DPHP have relatively unique molecular characteristics, sug-gesting differential modes of action and diabetogenic potentials, which are perhaps caused by their different chemical structures.

Putative AOPs Linking OPE Exposure and the Outcomes of T2D

AOPs can address the existing gaps in biological knowledge and provide evidence-based mechanistic insights into the associations between environmental exposure and health outcomes. To the best of our knowledge, our study is the first to present a putative network of AOPs linking OPE exposure with possible MIEs to down-stream KEs that ultimately manifest as outcomes of T2D. The AOPs presented herein revealed that OPE exposure could trigger various MIEs, for example, aberrant alterations of ligands and membrane receptors (TNF/FasL, IR, TNFR, and GLUT4) that can subsequently disturb the downstream signal transduction and metabolic pathways [e.g., PI3K-AKT signaling, IR signaling, PPARa=RXRa activation, and TCA cycle (as KE-1)] and biologi-cal processes [e.g., glucose metabolism (as KE-2)]. These path-ways primarily communicate through the PI3K-AKT signaling hub, which is a well-recognized major effector of metabolic insulin action. Specifically, IRs stimulate a signaling cascade, leading to the phosphorylation/activation of IR substrates and the subsequent activation and phosphorylation of the PI3K cascade and AKT, respectively.97 This further regulates GSK3 activation98 and GLUT4 expression on the cell membrane,99 contributing to glyco-gen synthesis and glucose transportation, respectively. Moreover, we observed broad perturbations in the downstream events of AKT phosphorylation, with abnormal pathophysiological functions (e.g., apoptosis, autophagy, oxidative stress, and inflammation), resulting in the outcome of T2D.100 Our findings are also in line with those of previous in silico and in vitro analyses showing that TPHP exposure can trigger the PI3K/AKT pathway as a KE,65 highlighting its potential role as a promising preventive interven-tionor pharmacological target for OPE-associated T2D.

Strengths and Limitations

The main strength of this study is the multidimensional integration of the transcriptome and metabolome using two different biologi-cal markers (i.e., blood and urine), providing a reliable assessment of OPE exposure and endogenous biomolecular perturbations while considering the critical aging window.101 Besides, the five longitudinal repeated measurements helped avoid biases intro-duced by exposure measurement errors of one-time sampling and presented the molecular scenarios of the effects of differential exposure at individual levels, with the potential to explore causal associations. In addition, the availability of detailed and longitudi-nal data pertaining to participant demographics and time-activity helped us adjust for confounding factors in the statistical analyses and minimize the effects of other risk factors on our findings.

This study also has some limitations. First, the complete exclu-sion of potential uncontrolled confounders was impossible because reliance on uncontrolled longitudinal associations limits the ability to make causal determinations. Second, the population representa-tion and sample size of our study could be inadequate, warranting further evaluations among different ethnic groups and larger popu-lation sizes with a nested case-control design or validations in toxi-cological research in vivo. Third, as potential EDCs, OPEs may have sexually dimorphic effects on glucose metabolism. Given the relatively small sample size and complex mixture of OPE expo-sures, it was underpowered to assess sex-stratified analyses in our exploratory study statistically. Future studies with a larger sample size may evaluate the sex-specific effects of these OPEs on glucose metabolism. Last, owing to the observational design of this study and the high dimensional feature of multi-omics data, caution is warranted with regard to whether these associations are causal.

Conclusions and Implications

The present exploratory study suggests that several OPEs may be independently and jointly associated with abnormal glycometabo-lism. High production volumes and the broad applications of OPEs in China, combined with the persistence and mobility characteris-tics, have led to the ever-increasing bioaccumulation of OPEs, pos-sibly contributing to the high risk of T2D among the Chinese population. The putative AOPs linking OPE exposure to T2D, including the activation of a series of MIEs (e.g., IR, TNFR, and PAR) and the disruption of multiple downstream KEs (e.g., PI3K/ AKT, PPARoc/RXRoc activation, and insulin secretion), were also established (Figure 7B). Our findings advance the existing knowl-edge on the environmental bases of the etiology of T2D among older Chinese adults 60-69 years of age and shed light on the evi-dence linking chemical safety and noncommunicable disease in public health.

Acknowledgments

We thank all the participants of the China Biomarkers of Air Pollutant Exposure (BAPE) Study, the Dianliu Community, the Ankang Community Hospital, the Shandong Center for Disease Control and Prevention (CDC), the Jinan CDC, the China BAPE Study team, and the Metabolomics Facilities in Tsinghua University Protein Research Center, as well as Hangzhou Calibra Diagnostics Ltd., Dian Diagnostics, and Metabolon.

This study was financially supported by the National Research Program for Key Issues in Air Pollution Control of China (DQGG0401), the National Natural Science Foundation of China (82025030, 81941023, and 92043301), and the National Key Research and Development Program of China (2016YFC0206500 and 2022YFC3702700) toX.S.

References

References

1. Weil AR. 2015. The growing burden of noncommunicable diseases. Health Aff (Millwood) 34(9):1439, PMID: 26355043, https://doi.org/10.1377/hlthaff.2015.0974.

2. WHO (World Health Organization). 2017. Noncommunicable Diseases Progress Monitor2017. https://www.who.int/publications/i/item/9789241513029 [accessed May 14,2022].

3. Sun H, Saeedi P, Karuranga S, Pinkepank M, Ogurtsova K, Duncan BB, et al. 2022. IDF Diabetes Atlas: global, regional and country-level diabetes preva-lence estimates for 2021 and projections for 2045. Diabetes Res Clin Pract 183:109119, PMID: 34879977, https://doi.org/10.1016/j.diabres.2021.109119.

4. Sinclair A, Saeedi P, Kaundal A, Karuranga S, Malanda B, Williams R. 2020. Diabetes and global ageing among 65-99-year-old adults: findings from the International Diabetes Federation Diabetes Atlas, 9th edition. Diabetes Res Clin Pract 162:108078, PMID: 32068097, https://doi.org/10.1016/j.diabres.2020.108078.

5. Zhang P, Gregg E. 2017. Global economic burden of diabetes and its implications. Lancet Diabetes Endocrinol 5(6):404-405, PMID: 28456417, https://doi.org/ 10.1016/S2213-8587(17)30100-6.

6. Stanc kov A, Laakso M. 2016. Genetics of type 2 diabetes. Endocr Dev 31:203-220, PMID: 26824439, https://doi.org/10.1159/000439418.

7. Kahn SE, Hull RL, Utzschneider KM. 2006. Mechanisms linking obesityto insulin resistance and type 2 diabetes. Nature 444(7121):840-846, PMID: 17167471, https://doi.org/10.1038/nature05482.

8. Zhang X, Zhang M, Zhao Z, Huang Z, Deng Q, Li Y, et al. 2020. Geographic var-iation in prevalence of adult obesity in China: results from the 2013-2014 National Chronic Disease and Risk Factor Surveillance. Ann Intern Med 172(4):291-293, PMID: 31658469, https://doi.org/10.7326/M19-0477.

9. Liu B, Du Y, Wu Y, Snetselaar LG, Wallace RB, Bao W. 2021. Trends in obesity and adiposity measures by race or ethnicity among adults in the United States 2011-2018: population based study. BMJ 5(1):21 27, PMID: 33727242, https://doi.org/10.1136/bmj.n365.

10. Wang L, Peng W, Zhao Z, Zhang M, Shi Z, Song Z, et al. 2021. Prevalence and treatment of diabetes in China, 2013-2018. JAMA 326(24):2498-2506, PMID: 34962526, https://doi.org/10.1001/jama.2021.22208.

11. Wang L, Li X, Wang Z, Bancks MP, Carnethon MR, Greenland P, et al. 2021. Trends in prevalence of diabetes and control of risk factors in diabetes among US adults, 1999-2018. JAMA 326(8):1-13, PMID: 34170288, https://doi.org/10. 1001/jama.2021.9883.

12. Hales CN, Barker DJ. 1992. Type 2 (non-insulin-dependent) diabetes mellitus: the thrifty phenotype hypothesis. Diabetologia 35(7):595-601, PMID: 1644236, https://doi.org/10.1007/BF00400248.

13. Li Y, He Y, Qi L, Jaddoe VW, Feskens EJM, Yang X, et al. 2010. Exposure to the Chinese famine in early life and the risk of hyperglycemia and type 2 diabetes in adulthood. Diabetes 59(10):2400-2406, PMID: 20622161, https://doi.org/10. 2337/db10-0385.

14. Lincoff AM, Nicholls SJ, Riesmeyer JS, Barter PJ, Brewer HB, Fox KAA, et al. 2017. Evacetrapib and cardiovascular outcomes in high-risk vascular disease. N Engl J Med 376(20):1933-1942, PMID: 28514624, https://doi.org/10.1056/NEJMoa1609581.

15. Li Y, Wang DD, Ley SH, Vasanti M, Howard AG, He Y, et al. 2017. Time trends of die-tary and lifestyle factors and their potential impact on diabetes burden in China. Diabetes Care 40(12):1685-1694, PMID: 29046327, https://doi.org/10.2337/dc17-0571.

16. Landrigan PJ, Fuller R, Acosta NJR, Adeyi O, Arnold R, Basu NN, et al. 2018. The Lancet Commission on pollution and health. Lancet 391(10119):462-512, PMID: 29056410, https://doi.org/10.1016/S0140-6736(17)32345-0.

17. Ruiz D, Becerra M, Jagai JS, Ard K, Sargis RM. 2018. Disparities in environmen-tal exposures to endocrine-disrupting chemicals and diabetes risk in vulnerable populations. Diabetes Care 41(1):193-205, PMID: 29142003, https://doi.org/10. 2337/dc16-2765.

18. Yehuda AB, Zinger A, Durso S. 2014. The older patient with diabetes: a practical approach. Diabetes Metab Res Rev 30(2):88-95, PMID: 24123811, https://doi.org/ 10.1002/dmrr.2485.

19. Lind L, Zethelius B, Salihovic S, van Bavel B, Lind PM. 2014. Circulating lev-els of perfluoroalkyl substances and prevalent diabetes in the elderly. Diabetologia 57(3):473-479, PMID: 24337155, https://doi.org/10.1007/s00125-013-3126-3.

20. Lind PM, Zethelius B, Lind L. 2012. Circulating levels of phthalate metabolites are associated with prevalent diabetes in the elderly. Diabetes Care 35(7):1519- 1524, PMID: 22498808, https://doi.org/10.2337/dc11-2396.

21. Choi YH, Kim JH, Hong YC. 2015. Sex-dependent and body weight-dependent associations between environmental PAHs exposure and insulin resistance: Korean urban elderly panel. J Epidemiol Community Health 69(7):625-631, PMID: 25669219, https://doi.org/10.1136/jech-2014-204801.

22. Blum A, Behl M, Birnbaum L, Diamond ML, Phillips A, Singla V, et al. 2019. Organophosphate ester flame retardants: are they a regrettable substitution for polybrominated diphenyl ethers? Environ Sci Technol Lett 6(11):638-649, PMID: 32494578, https://doi.org/10.1021/acs.estlett.9b00582.

23. Hu Z, Yin L, Wen X, Jiang C, Long Y, Zhang J, et al. 2021. Organophosphate esters in China: fate, occurrence, and human exposure. Toxics 9(11):310, PMID: 34822701, https://doi.org/10.3390/toxics9110310.

24. Hu L, Yu M, Li Y, Liu L, Li X, Song L, et al. 2022. Association of exposure to organophosphate esters with increased blood pressure in children and ado-lescents. Environ Pollut 295:118685, PMID: 34923060, https://doi.org/10.1016/j. envpol.2021.118685.

25. Hu W, Gao P, Wang L, Hu J. 2023. Endocrine disrupting toxicity of aryl organo-phosphate esters and mode of action. Crit Rev Environ Sci Technol 53(1):1-18, https://doi.org/10.1080/10643389.2022.2050147.

26. Yang Y, Chen P, Ma S, Lu S, Yu Y, An T. 2022. A critical review of human internal exposure and the health risks of organophosphate ester flame retardants and their metabolites. Crit Rev Environ Sci Technol 52(9):1528-1560, https://doi.org/ 10.1080/10643389.2020.1859307.

27. Wang C, Le Y, Lu D, Zhao M, Dou X, Zhang Q, et al. 2020. Triphenyl phosphate causes a sexually dimorphic metabolism dysfunction associated with disor-dered adiponectin receptors in pubertal mice. J Hazard Mater 388:121732, PMID: 31796355, https://doi.org/10.1016/j.jhazmat.2019.121732.

28. Du Z, Zhang Y, Wang G, Peng J, Wang Z, Gao S. 2016. TPhP exposure disturbs carbohydrate metabolism, lipid metabolism, and the DNA damage repair system in zebrafish liver. Sci Rep 6:21827, PMID: 26898711, https://doi.org/10.1038/ srep21827.

29. Green AJ, Graham JL, Gonzalez EA, La Frano MR, Petropoulou SSE, Park JS, et al. 2017. Perinatal triphenyl phosphate exposure accelerates type 2 diabetes onset and increases adipose accumulation in UCD-type 2 diabetes melli-tus rats. Reprod Toxicol 68:119-129, PMID: 27421578, https://doi.org/10.1016/j. reprotox.2016.07.009.

30. Bo Y, Zhu Y. 2022. Organophosphate esters exposure in relation to glucose homeostasis and type 2 diabetes in adults: a national cross-sectional study from the National Health and Nutrition Survey. Chemosphere 301:134669, PMID: 35460677, https://doi.org/10.1016/j.chemosphere.2022.134669.

31. Luo K, Aimuzi R, Wang Y, Nian M, Zhang J. 2020. Urinary organophosphate esters metabolites, glucose homeostasis and prediabetes in adolescents. Environ Pollut 267:115607, PMID: 33254666, https://doi.org/10.1016/j.envpol.2020. 115607.

32. Wang S, Yong H, He XD. 2021. Multi-omics: opportunities for research on mechanism of type 2 diabetes mellitus. World J Diabetes 12(7):1070-1080, PMID: 34326955, https://doi.org/10.4239/wjd.v12.i7.1070.

33. Beulens JWJ, Pinho MGM, Abreu TC, den Braver NR, Lam TM, Huss A, et al. 2022. Environmental risk factors of type 2 diabetes-an exposome approach. Diabetologia 65(2):263-274, PMID: 34792619, https://doi.org/10.1007/s00125-021-05618-w.

34. Price EJ, Vitale CM, Miller GW, David A, Barouki R, Audouze K, et al. 2022. Merging the exposome into an integrated frameworkfor "omics" sciences. iScience 25(3):103976, PMID: 35310334, https://doi.org/10.1016/j.isci.2022. 103976.

35. Tang S, Li T, Fang J, Chen R, Cha Y, Wang Y, et al. 2021. The exposome in practice: an exploratory panel study of biomarkers of air pollutant exposure in Chinese people aged 60-69 years (China BAPE Study). Environ Int 157:106866, PMID: 34525388, https://doi.org/10.1016/j.envint.2021.106866.

36. Hou M, Fang J, Shi Y, Tang S, Dong H, Liu Y, et al. 2021. Exposure to organophosphate esters in elderly people: relationships of OPE body burdens with indoor air and dust concentrations and food consumption. Environ Int 157:106803, PMID: 34365320, https://doi.org/10.1016/j.envint.2021.106803.

37. Matthews DR, Hosker JP, Rudenski AS, Naylor BA, Treacher DF, Turner RC. 1985. Homeostasis model assessment: insulin resistance and beta-cell func-tion from fasting plasma glucose and insulin concentrations in man. Diabetologia 28(7):412-419, PMID: 3899825, https://doi.org/10.1007/BF00280883.

38. Kim D, Langmead B, Salzberg SL 2015. HISAT: a fastspliced alignerwith lowmem-ory requirements. Nat Methods 12(4):357-360, PMID: 25751142, https://doi.org/10. 1038/nmeth.3317.

39. Liao Y, Smyth GK, Shi W. 2014. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7):923-930, PMID: 24227677, https://doi.org/10.1093/bioinformatics/btt656.

40. Jia X, Jin Q, Fang J, Shi Y, Hou M, Dong H, et al. 2022. Emerging and legacy perand polyfluoroalkyl substances in an elderly population in Jinan, China: the exposure level, short-term variation, and intake assessment. Environ Sci Technol 56(12):7905-7916, PMID: 35584234, https://doi.org/10.1021/acs.est.2c00381.

41. Keil AP, Buckley JP, O'Brien KM, Ferguson KK, Zhao S, White AJ, et al. 2020. A quantile-based g-computation approach to addressing the effects of exposure mixtures. Environ Health Perspect 128(4):47004, PMID: 32255670, https://doi.org/10.1289/EHP5838.

42. Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genet 10:23, PMID: 19473544, https://doi.org/10.1186/1471-2156-10-23.

43. Kr mer A, Green J, Pollard J Jr, Tugendreich S. 2014. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 30(4):523-530, PMID: 24336805, https://doi.org/10.1093/bioinformatics/btt703.

44. White IR, Royston P, Wood AM. 2011. Multiple imputation using chained equations: issues and guidance for practice. Stat Med 30(4):377-399, PMID: 21225900, https://doi.org/10.1002/sim.4067.

45. Boyle M, Buckley JP, Quir s-Alcal L 2019. Associations between urinary organophosphate ester metabolites and measures of adiposity among U.S. children and adults: NHANES 2013-2014. Environ Int 127:754-763, PMID: 31003058, https://doi.org/10.1016/j.envint.2019.03.055.

46. Ospina M, Jayatilaka NK, Wong LY, Restrepo P, Calafat AM. 2018. Exposure to organophosphate flame retardant chemicals in the U.S. general population: data from the 2013-2014 National Health and Nutrition Examination Survey. Environ Int 110:32-41, PMID: 29102155, https://doi.org/10.1016/j.envint.2017.10.001.

47. Storey JD. 2002. A direct approach to false discovery rates. J R Stat Soc Series B Stat Methodol 64(3):479-498, https://doi.org/10.1111/1467-9868.00346.

48. Liu X, Zhao X, Wang Y, Hong J, Shi M, Pfaff D, et al. 2020. Triphenyl phosphate per-meates the blood brain barrier and induces neurotoxicity in mouse brain. Chemosphere 252:126470, PMID: 32443258, https://doi.org/10.1016/j.chemosphere. 2020.126470.

49. Wang D,Yan S,Yan J,Teng M, Meng Z, Li R, etal. 2019. Effects of triphenyl phosphate exposure during fetal development on obesity and metabolic dysfunctions in adult mice: impaired lipid metabolism and intestinal dysbiosis. Environ Pollut 246:630-638, PMID: 30605818, https://doi.org/10.1016/j.envpol.2018.12.053.

50. Wang L, Huang X, Laserna AKC, Li SFY. 2018. Untargeted metabolomics reveals transformation pathways and metabolic response of the earthworm Perionyx excavatus after exposure to triphenyl phosphate. Sci Rep 8(1):16440, PMID: 30401822, https://doi.org/10.1038/s41598-018-34814-9.

51. Philbrook NA, Restivo VE, Belanger CL, Winn LM. 2018. Gestational triphenyl phosphate exposure in C57Bl/6 mice perturbs expression of insulin-like growth factor signaling genes in maternal and fetal liver. Birth Defects Res 110(6):483-494, PMID: 29316351, https://doi.org/10.1002/bdr2.1185.

52. Wang D, Zhu W, Chen L, Yan J, Teng M, Zhou Z, et al. 2018. Neonatal triphenyl phosphate and its metabolite diphenyl phosphate exposure induce sex- and dose-dependent metabolic disruptions in adult mice. Environ Pollut 237:10-17, PMID: 29466770, https://doi.org/10.1016/j.envpol.2018.01.047.

53. Adams S, Wiersielis K, Yasrebi A, Conde K, Armstrong L, Guo GL, et al. 2020. Sex- and age-dependent effects of maternal organophosphate flame-retardant exposure on neonatal hypothalamic and hepatic gene expression. Reprod Toxicol 94:65-74, PMID: 32360330, https://doi.org/10.1016/j.reprotox.2020.04.001.

54. Walley SN, Krumm EA, Yasrebi A, Kwiecinski J, Wright V, Baker C, et al. 2021. Maternal organophosphate flame-retardant exposure alters offspring energy and glucose homeostasis in a sexually dimorphic manner in mice. J Appl Toxicol 41(4):572-586, PMID: 32969501, https://doi.org/10.1002/jat.4066.

55. Vail GM, Walley SN, Yasrebi A, Maeng A, Conde KM, Roepke TA. 2020. The interactions of diet-induced obesity and organophosphate flame retardant exposure on energy homeostasis in adult male and female mice. J Toxicol Environ Health A 83(11-12):438-455, PMID: 32546061, https://doi.org/10.1080/ 15287394.2020.1777235.

56. Krumm EA, Patel VJ, Tillery TS, Yasrebi A, Shen J, Guo GL, et al. 2018. Organophosphate flame-retardants alter adult mouse homeostasis and gene expression in a sex-dependent manner potentially through interactions with ERa. Toxicol Sci 162(1):212-224, PMID: 29112739, https://doi.org/10.1093/toxsci/ kfx238.

57. Vail GM, Walley SN, Yasrebi A, Maeng A, Degroat TJ, Conde KM, et al. 2022. Implications of peroxisome proliferator-activated receptor gamma (PPARY) with the intersection of organophosphate flame retardants and diet-induced obesity in adult mice. J Toxicol Environ Health A 85(9):381-396, PMID: 35000574, https://doi.org/10.1080/15287394.2021.2023716.

58. Vail GM, Walley SN, Yasrebi A, Maeng A, Degroat TJ, Conde KM, et al. 2022. Implications of estrogen receptor alpha (ERa) with the intersection of organo-phosphate flame retardants and diet-induced obesity in adult mice. J Toxicol Environ Health A 85(10):397-413, PMID: 35045790, https://doi.org/10.1080/ 15287394.2022.2026849.

59. Liu Y, Le Y, Xu M, Wang W, Chen H, Zhang Q, et al. 2022. Remodeling on adi-pocytic physiology of organophosphorus esters in mature adipocytes. Environ Pollut 305:119287, PMID: 35421551, https://doi.org/10.1016/j.envpol. 2022.119287.

60. Gu J, Su F, Hong P, Zhang Q, Zhao M. 2019. 1H NMR-based metabolomic anal-ysis of nine organophosphate flame retardants metabolic disturbance in Hep G2 cell line. Sci Total Environ 665:162-170, PMID: 30772545, https://doi.org/10. 1016/j.scitotenv.2019.02.055.

61. Cano-Sancho G, Smith A, La Merrill MA. 2017. Triphenyl phosphate enhances adipogenic differentiation, glucose uptake and lipolysis via endocrine and noradrenergic mechanisms. Toxicol In Vitro 40:280-288, PMID: 28163246, https://doi.org/10.1016/j.tiv.2017.01.021.

62. Hu W, Kang Q, Zhang C, Ma H, Xu C, Wan Y, et al. 2020. Triphenyl phosphate modulated saturation of phospholipids: induction of endoplasmic reticulum stress and inflammation. Environ Pollut 263(pt A):114474, PMID: 32259740, https://doi.org/10.1016/j.envpol.2020.114474.

63. Liu Y, Li Y, Dong S, Han L, Guo R, Fu Y, et al. 2021. The risk and impact of organophosphate esters on the development of female-specific cancers: comparative analysis of patients with benign and malignant tumors. J Hazard Mater 404(pt B):124020, PMID: 33049558, https://doi.org/10.1016/j.jhazmat.2020. 124020.

64. Ya M, Yu N, Zhang Y, Su H, Tang S, Su G. 2019. Biomonitoring of organophos-phate triesters and diesters in human blood in Jiangsu Province, eastern China: occurrences, associations, and suspect screening of novel metabo-lites. Environ Int 131:105056, PMID: 31369981, https://doi.org/10.1016/j.envint. 2019.105056.

65. Wang X, Wang L, Li F, Teng Y, Ji C, Wu H. 2022. Toxicity pathways of lipid meta-bolic disorders induced by typical replacement flame retardants via data-driven analysis, in silico and in vitro approaches. Chemosphere 287(pt 4):132419, PMID: 34600017, https://doi.org/10.1016/j.chemosphere.2021.132419.

66. Batista TM, Jayavelu AK, Wewer Albrechtsen NJ, Iovino S, Lebastchi J, Pan H, et al. 2020. A cell-autonomous signature of dysregulated protein phospho-rylation underlies muscle insulin resistance in type 2 diabetes. Cell Metab 32(5):844-859.e5, PMID: 32888406, https://doi.org/10.1016/j.cmet.2020.08.007.

67. Sacco F, Seelig A, Humphrey SJ, Krahmer N, Volta F, Reggio A, et al. 2019. Phosphoproteomics reveals the GSK3-PDX1 axis as a key pathogenic signal-ing node in diabetic islets. Cell Metab 29(6):1422-1432.e3, PMID: 30879985, https://doi.org/10.1016/j.cmet.2019.02.012.

68. Yu X, Yin H, Peng H, Lu G, Liu Z, Dang Z. 2019. OPFRs and BFRs induced A549 cell apoptosis by caspase-dependent mitochondrial pathway. Chemosphere 221:693-702, PMID: 30669111, https://doi.org/10.1016/j.chemosphere.2019.01.074.

69. Meng Y, Xu X, Niu D, Xu Y, Qiu Y, Zhu Z, et al. 2022. Organophosphate flame retardants induce oxidative stress and Chop/Caspase 3-related apoptosis via Sod1/p53/Map3k6/Fkbp5 in NCI-1975 cells. Sci Total Environ 819:153160, PMID: 35051466, https://doi.org/10.1016/j.scitotenv.2022.153160.

70. Deng Q, Jiang L, Mao L, Song XH, He CQ, Li XL, et al. 2020. The role of protein kinase C alpha in tri-ortho-cresyl phosphate-induced autophagy in human neuroblastoma SK-N-SH cells. J Appl Toxicol 40(11):1480-1490, PMID: 33020912, https://doi.org/10.1002/jat.3999.

71. Chen R, Hou R, Hong X, Yan S, Zha J. 2019. Organophosphate flame retardants (OPFRs) induce genotoxicity in vivo: a survey on apoptosis, DNA methylation, DNA oxidative damage, liver metabolites, and transcriptomics. Environ Int 130:104914, PMID: 31226563, https://doi.org/10.1016/j.envint.2019.104914.

72. de Guia RM, Rose AJ, Herzig S. 2014. Glucocorticoid hormones and energy homeostasis. Horm Mol Biol Clin Investig 19(2):117-128, PMID: 25390020, https://doi.org/10.1515/hmbci-2014-0021.

73. Kuo T, McQueen A, Chen TC, Wang JC. 2015. Regulation of glucose homeo-stasis by glucocorticoids. Adv Exp Med Biol 872:99-126, PMID: 26215992, https://doi.org/10.1007/978-1-4939-2895-8_5.

74. Kang S, Tsai LT, Zhou Y, Evertts A, Xu S, Griffin MJ, et al. 2015. Identification of nuclear hormone receptor pathways causing insulin resist-ance by transcriptional and epigenomic analysis. Nat Cell Biol 17(1):44-56, PMID: 25503565, https://doi.org/10.1038/ncb3080.

75. Umamaheswari S, Karthika P, Suvenitha K, Kadirvelu K, Ramesh M. 2021. Dose-dependent molecular responses of Labeo rohita to triphenyl phosphate. Chem Res Toxicol 34(12):2500-2511, PMID: 34847329, https://doi.org/10.1021/ acs.chemrestox.1c00281.

76. Ramesh M, Angitha S, Haritha S, Poopal RK, Ren Z, Umamaheswari S. 2020. Organophosphorus flame retardant induced hepatotoxicity and brain AChE inhibition on zebrafish (Danio rerio). Neurotoxicol Teratol 82:106919, PMID: 32853706, https://doi.org/10.1016/j.ntt.2020.106919.

77. Canbaz D, Logiantara A, van Ree R, van Rijt LS. 2017. Immunotoxicity of organophosphate flame retardants TPHP and TDCIPP on murine dendritic cells in vitro. Chemosphere 177:56-64, PMID: 28282624, https://doi.org/10.1016/ j.chemosphere.2017.02.149.

78. Van den Eede N, Cuykx M, Rodrigues RM, Laukens K, Neels H, Covaci A, et al. 2015. Metabolomics analysis of the toxicity pathways of triphenyl phosphate in HepaRG cells and comparison to oxidative stress mechanisms caused by acet-aminophen. Toxicol In Vitro 29(8):2045-2054, PMID: 26318275, https://doi.org/10. 1016/j.tiv.2015.08.012.

79. Mennillo E, Cappelli F, Arukwe A. 2019. Biotransformation and oxidative stress responses in rat hepatic cell-line (H4IIE) exposed to organophosphate esters (OPEs). Toxicol Appl Pharmacol 371:84-94, PMID: 30974155, https://doi.org/10. 1016/j.taap.2019.04.004.

80. Chen G, Zhang S, Jin Y, Wu Y, Liu L, Qian H, et al. 2015. TPP and TCEP induce oxidative stress and alter steroidogenesis in TM3 Leydig cells. Reprod Toxicol 57:100-110, PMID: 26049154, https://doi.org/10.1016/j.reprotox.2015.05.011.

81. Wang X, Li F, Liu J, Ji C, Wu H. 2020. Transcriptomic, proteomic and metabolo-mic profiling unravel the mechanisms of hepatotoxicity pathway induced by tri-phenyl phosphate (TPP). Ecotoxicol Environ Saf 205:111126, PMID: 32823070, https://doi.org/10.1016/j.ecoenv.2020.111126.

82. Schang G, Robaire B, Hales BF. 2016. Organophosphate flame retardants act as endocrine-disrupting chemicals in MA-10 mouse tumor Leydig cells. Toxicol Sci 150(2):499-509, PMID: 26794138, https://doi.org/10.1093/toxsci/kfw012.

83. An J, Hu J, Shang Y, Zhong Y, Zhang X, Yu Z. 2016. The cytotoxicity of organo-phosphate flame retardants on HepG2, A549 and Caco-2 cells. J Environ Sci Health A Tox Hazard Subst Environ Eng 51(11):980-988, PMID: 27336727, https://doi.org/10.1080/10934529.2016.1191819.

84. Chen G, Jin Y, Wu Y, Liu L, Fu Z. 2015. Exposure of male mice to two kinds of organophosphate flame retardants (OPFRs) induced oxidative stress and endocrine disruption. Environ Toxicol Pharmacol 40(1):310-318, PMID: 26183808, https://doi.org/10.1016/j.etap.2015.06.021.

85. Yuan S, Zhu K, Ma M, Zhu X, Rao K, Wang Z. 2020. In vitro oxidative stress, mitochondrial impairment and G1 phase cell cycle arrest induced by alkyl-phosphorus-containing flame retardants. Chemosphere 248:126026, PMID: 32006839, https://doi.org/10.1016/j.chemosphere.2020.126026.

86. Yao Y, Li M, Pan L, Duan Y, Duan X, Li Y, et al. 2021. Exposure to organophos-phate ester flame retardants and plasticizers during pregnancy: thyroid endo-crine disruption and mediation role of oxidative stress. Environ Int 146:106215, PMID: 33113466, https://doi.org/10.1016/j.envint.2020.106215.

87. Wang Y, Hong J, Shi M, Guo L, Liu L, Tang H, et al. 2021. Triphenyl phosphate disturbs the lipidome and induces endoplasmic reticulum stress and apopto-sis in JEG-3 cells. Chemosphere 275:129978, PMID: 33662732, https://doi.org/ 10.1016/j.chemosphere.2021.129978.

88. Hong J, Jiang M, Guo L, Lin J, Wang Y, Tang H, et al. 2022. Prenatal exposure to triphenyl phosphate activated PPAR? in placental trophoblasts and impaired pregnancy outcomes. Environ Pollut 301:119039, PMID: 35192884, https://doi.org/10.1016/j.envpol.2022.119039.

89. Wang X, Li F, Liu J, Li Q, Ji C, Wu H. 2021. New insights into the mechanism of hepatocyte apoptosis induced by typical organophosphate ester: an integrated in vitro and in silico approach. Ecotoxicol Environ Saf 219:112342, PMID: 34023725, https://doi.org/10.1016/j.ecoenv.2021.112342.

90. Bowen C, Childers G, Perry C, Martin N, McPherson CA, Lauten T, et al. 2020. Mitochondrial-related effects of pentabromophenol, tetrabromobisphenol A, and triphenyl phosphate on murine BV-2 microglia cells. Chemosphere 255:126919, PMID: 32402876, https://doi.org/10.1016/j.chemosphere.2020.126919.

91. Zhang Q, Wang J, Zhu J, Liu J, Zhao M. 2017. Potential glucocorticoid and mineralocorticoid effects of nine organophosphate flame retardants. Environ Sci Technol 51(10):5803-5810, PMID: 28430429, https://doi.org/10.1021/acs.est. 7b01237.

92. KojimaH,Takeuchi S,Van den EedeN,CovaciA.2016. Effectsofprimary metabo-lites of organophosphate flame retardants on transcriptional activity via human nuclear receptors. Toxicol Lett 245:31-39, PMID: 26778350, https://doi.org/10. 1016/j.toxlet.2016.01.004.

93. Kojima H, Takeuchi S, Itoh T, Iida M, Kobayashi S, Yoshida T. 2013. In vitro endocrine disruption potential of organophosphate flame retardants via human nuclear receptors. Toxicology 314(1):76-83, PMID: 24051214, https://doi.org/10. 1016/j.tox.2013.09.004.

94. Wang B, Chandrasekera PC, Pippin JJ. 2014. Leptin- and leptin receptor-deficient rodent models: relevance for human type 2 diabetes. Curr Diabetes Rev 10(2):131- 145, PMID: 24809394, https://doi.org/10.2174/1573399810666140508121012.

95. Bani D, Pini A, Yue SKS. 2012. Relaxin, insulin and diabetes: an intriguing con-nection. Curr Diabetes Rev 8(5):329-335, PMID: 22698078, https://doi.org/10. 2174/157339912802083487.

96. Yang YM, Kuen DS, Chung Y, Kurose H, Kim SG. 2020. Ga12/13 signaling in met-abolic diseases. Exp Mol Med 52(6):896-910, PMID: 32576930, https://doi.org/ 10.1038/s12276-020-0454-5.

97. Cartee GD. 2015. Mechanisms for greater insulin-stimulated glucose uptake in normal and insulin-resistant skeletal muscle after acute exercise. Am J Physiol Endocrinol Metab 309(12):E949-E959, PMID: 26487009, https://doi.org/ 10.1152/ajpendo.00416.2015.

98. Shieh JM, Wu HT, Cheng KC, Cheng JT. 2009. Melatonin ameliorates high fat diet-induced diabetes and stimulates glycogen synthesis via a PKC^-Akt-GSK3p pathway in hepatic cells. J Pineal Res 47(4):339-344, PMID: 19817973, https://doi.org/10.1111/j.1600-079X.2009.00720.x.

99. Saltiel AR, Kahn CR. 2001. Insulin signalling and the regulation of glucose and lipid metabolism. Nature 414(6865):799-806, PMID: 11742412, https://doi.org/10. 1038/414799a.

100. Huang X, Liu G, Guo J, Su Z. 2018. The PI3K/AKT pathway in obesity and type 2 diabetes. Int J Biol Sci 14(11):1483-1496, PMID: 30263000, https://doi.org/10. 7150/ijbs.27173.

101. Lehallier B, Gate D, Schaum N, Nanasi T, Lee SE, Yousef H, et al. 2019. Undulating changes in human plasma proteome profiles across the lifespan. Nat Med 25(12):1843-1850, PMID: 31806903, https://doi.org/10.1038/s41591-019-0673-2.