Metabolomic and Transcriptomic Analysis of MCF-7 Cells Exposed to 23 Chemicals at Human-Relevant Levels: Estimation of Individual Chemical Contribution to Effects

Abstract

Background: Humans are constantly being exposed to various xenobiotics at relatively low concentrations. To date, limited evidence is available to ascertain whether a complex xenobiotic mixture at human-relevant levels causes any health effect. Moreover, there is no effective method to pinpoint the contribution of each chemical toward such an effect.

Objectives: This study aims to understand the responses of cells to a mixture containing 23 xenobiotics at human-relevant levels and develop a feasible method to decipher the chemical(s) that contribute significantly to the observed effect.

Methods: We characterized the metabolome and transcriptome of breast cancer cells (MCF-7) before and after exposure to the mixture at human-relevant levels; preexposure levels were derived from existing large-scale biomonitoring data. A high-throughput metabolomics-based “leave-one-out” method was proposed to understand the relative contribution of each component by comparing the metabolome with and without the particular chemical in the mixture.

Results: The metabolomic analysis suggested that the mixture altered metabolites associated with cell proliferation and oxidative stress. For the transcriptomes, gene ontology terms and pathways including “cell cycle,” “cell proliferation,” and “cell division” were significantly altered after mixture exposure. The mixture altered genes associated with pathways such as “genotoxicity” and “nuclear factor erythroid 2-related factor 2 (Nrf2).” Through joint pathways analysis, metabolites and genes were observed to be well-aligned in pyrimidine and purine metabolisms. The leave-one-out results showed that many chemicals made their contributions to specific metabolic pathways. The overall metabolome pattern of the absence of 2,4-dihyroxybenzophenone (DHB) or bisphenol A (BPA) showed great resemblance to controls, suggesting their higher relative contribution to the observed effect.

Discussion: The omics results showed that exposure to the mixture at human-relevant levels can induce significant in vitro cellular changes. Also, the leave one out method offers an effective approach for deconvoluting the effects of the mixture.

Introduction

Humans are ubiquitously exposed to many xenobiotics on a day-to-day basis. These chemicals include personal care products (PCPs), pharmaceuticals, plasticizers, and dietary additives, many of which are known to be toxic to humans and wildlife (Boxall et al. 2012; Guo and Kannan 2013). In a recent study, we listed over 300 chemicals found in indoor dust, along with their concentrations; the data indicated the widespread occurrences of environmental contaminants (Dong et al. 2019). Another study reported the co-occurrence of more than 1,000 small molecules in human blood (Rappaport et al. 2014), suggesting exposure to chemicals with vast complexity. Furthermore, most environmental contaminants have very low biological concentrations, such as sub-parts-per-billion (ppb), parts-per-trillion (ppt) or even lower (Rappaport et al. 2014). Therefore, the question of whether a “chemical mixture” of higher complexity at low human-relevant levels would trigger any health effects remains unaddressed, and it is an ongoing challenge for toxicological risk assessments.

Although single contaminants may induce notable effects under the circumstances of high exposure, most people are exposed to chemical mixtures at lower concentrations (Carpenter et al. 2002). The combined effects of a chemical mixture have been recognized as a concern, due to the possibility that chemicals have additive or even synergistic effects (Carpenter et al. 2002; Kortenkamp 2007). Several earlier studies on mixtures revealed the ability of chemicals to act together to induce some observable effect at low or environmentally relevant levels. For example, a combination of eight weak estrogenic chemicals at concentrations below the no observed effect concentration (NOEC) produced significant mixture effects in a yeast recombinant screen (Silva et al. 2002). In an in vivo study, five estrogenic chemicals at low-effect concentrations demonstrated the capacity to act in combination to affect the fitness and reproductive performance of fish (Brian et al. 2007). In another example, results showed that the metabolomes of four representative human gut microbiota were significantly perturbed by exposure to xenobiotic mixtures at human-relevant levels (Zhang et al. 2018). However, the majority of these studies usually focused on limited compounds with the same mode of action (MOA) and tested a single toxicological end point such as the estrogenic effect (Brian et al. 2007; Repouskou et al. 2019; Silva et al. 2002).

Research to date has not yet determined the possible effects of a complex mixture that involves more chemicals, particularly the ever-increasing organic contaminants with diverse MOA. Organic contaminants such as flame retardants, phthalates, and bisphenols contribute substantially to the human exposure, due to their widespread applications and occurrences in the environment (Asimakopoulos et al. 2014; Lioy et al. 2002; Liu et al. 2019; Meeker et al. 2013; Wang et al. 2015c; Yang et al. 2014); however, there have been very few studies on mixtures containing these contaminants. Additionally, a major challenge in the study of mixtures is evaluating the contribution of each component to the total observed effect (Carpenter et al. 2002). It is more challenging for mixtures at human-relevant levels considering that the effect triggered by an individual component may be too small to be distinguishable. The chemical interactions between the components, each with a distinct MOA, can be very complex, resulting in a heavy experimental workload. Therefore, it is not feasible to use the traditional bioassay with a specific end point or effect-directed analysis, which usually considers a couple of chemicals as the dominant drivers in the mixture effect to ascertain the causal compounds (Fang et al. 2014; Fang et al. 2015b). There is still a pressing need to develop a more effective proof-of-concept to elucidate the contribution of each component in a mixture.

Omics technologies provide a global view and quantitative information on changes in critical cellular components under the influence of environmental factors (Rappaport and Smith 2010) and are ideal to address these questions. In particular, metabolomics, a method of understanding metabolic regulations by studying naturally occurring, low molecular weight (LMW) organic metabolites at a specified time point, can be regarded as the end point of the “omics cascade” (Dettmer and Hammock 2004). Providing the ultimate response of an organism to environmental influences, metabolome changes enable us to understand both phenotypical fingerprints and affected biochemical pathways (Weckwerth 2003). Additionally, the ability to detect even slight changes in endogenous metabolites makes metabolomics suitable for complex mixtures with multiple outcomes or cases where growth phenotypes are lacking. The integration of omics technologies offers a new solution for examining the biological processes triggered by mixtures at the human-relevant levels and prioritizing the relative contribution of each component.

To address these questions, we used omics technologies (i.e., metabolomics and transcriptomics) to fully characterize the biological effects triggered by a xenobiotic mixture of 23 chemicals, at human-relevant levels, on breast cancer cells (MCF-7), a known sensitive hormone-responsive and suitable cell model to study endocrine-disrupting chemicals (EDCs). This mixture covered 23 organic contaminants with diverse MOAs and included bisphenols, phthalates, PCPs, flame retardants, and other xenobiotics. These chemicals were selected based on several criteria, including high-volume manufacturing, widespread occurrences in environmental and human samples, the exhibition of toxicity in vivo, and endocrine-disrupting behaviors. The dose concentrations of these xenobiotics were drawn from existing large-cohort biomonitoring data as detailed in Table 1. A high-throughput metabolomics-based leave one out method was also developed to evaluate the relative contribution of each component in the mixture. Provided the molecular change signal is observable during all experimental manipulations, this method may effectively characterize and compare the changes in the mixture effect before and after a single compound is removed. The main objectives of this study are to develop a platform and “proof-of-principle” using omics technologies for mixture effects characterization and evaluation of the relative contribution of each component to the observed effects.

Table 1 The GM, 95th percentile, and maximum (Max) urinary or blood concentrations at human-relevant levels for the 23 chemicals.

Table 1 has six main columns, namely, Number, Compounds Name (Abbreviation if available), Urine or Blood Concentrations (nanogram per milliliter) begin superscript lowercase a end superscript, Exposure Matrix, Sample Size (uppercase n), and References. Urine or Blood Concentrations (nanogram per milliliter) begin superscript lowercase a end superscript column has five sub columns, namely, Maximum, Percentage (percentage) in maximum mixture begin superscript lowercase b end superscript, Geometric Mean, and Ninety Fifth Percentile.

No.	Compound name (Abbreviation if available)	Urine or blood concentrations (ng/mL)a					Exposure Matrix	Sample Size (n)	References
		Maximum (Max)	Percentage (%) in Max mixtureb	GM	Percentage (%) in GM mixture	95th Percentile
1	2,4-dihydroxybenzophenone (DHB)c	1,911.7	36.9	21.9	10.7	1,191.8	Urine	15,593	CDC 2018
2	Bisphenol A (BPA)	12.9	0.2	1.9	0.9	11.2	Urine	15,593	CDC 2018
3	Bisphenol S (BPS)	4.7	0.1	0.4	0.2	3.6	Urine	2,682	CDC 2018
4	Triclocarban	24.5	0.3	0.2e	0.1	13.4	Urine	2,686	CDC 2018
5	Triclosan (TCS)	599.7	8.6	13.8	5.0	504.0	Urine	15,593	CDC 2018
6	Butyl paraben	18.28	0.4	1.2e	0.6	13.9	Urine	13,076	CDC 2018
7	Ethyl paraben	99.6	2.5	5.4e	3.4	71.5	Urine	13,076	CDC 2018
8	Methyl paraben	1,056.6	28.8	52.4	35.9	909.2	Urine	13,076	CDC 2018
9	Tris(1,3-dichloroisopropyl)phosphate (TDCPP)d	8.6	0.1	0.9	0.2	7.2	Urine	2,646	CDC 2018
10	Triphenyl phosphate (TPhP)d	16.6	0.2	2.7	0.9	NA	Blood	69	Li et al. 2015, 2017; Zhao et al. 2017
11	Di-n-butyl phthalate (DnBP)d	13.9	0.2	3.5	1.3	NA	Blood	184	Frederiksen et al. 2010; Högberg et al. 2008; Li et al. 2018a
12	Mono-n-butyl phthalate (MnBP)	110.2	2.1	16.8	7.9	97.1	Urine	2,1003	CDC 2018
13	Benzylbutyl phthalate (BzBP)d	232.8	3.1	16.7	5.6	168.8	Urine	18,462	CDC 2018
14	Mono-2-ethylhexyl phthalate (MEHP)	31.3	0.5	2.7	1.0	23.8	Urine	21,003	CDC 2018
15	Mono-benzyl phthalate (MBzP)	70.0	1.1	8.3	3.4	59.6	Urine	21,003	CDC 2018
16	Diisononyl phthalate (DiNP)d	13.8	0.1	2.0e	0.5	8.4	Urine	21,003	CDC 2018
17	Diisodecyl phthalate (DiDP)d	20.9	0.2	2.7	0.6	16.6	Urine	13,075	CDC 2018
18	Genistein	781.3	12.0	29.7	11.5	598.2	Urine	15,689	CDC 2018
19	Perfluorooctanoic acid (PFOA)	10.7	0.1	3.5	0.9	8.8	Blood	14,178	CDC 2018
20	Bisphenol A Diglycidyl ether (BADGE)	2.5	0.0	0.5	0.2	NA	Blood	223	Kim et al. 2015
21	Bisphenol A bis(2,3-dihydroxypropyl) ether (BADGE-2H2O )	9.5	0.1	2.3	0.6	NA	Blood	243	Kim et al. 2015; Wang et al. 2015a
22	Bisphenol-F-diglycidyl ether (BFDGE)	180	2.4	26.2	8.7	NA	Blood	20	Wang et al. 2015a
23	Tetrabromobisphenol A (TBBPA)	7.6	0.1	0.7	0.1	NA	Blood	265	Nagayama et al. 2000; Thomsen et al. 2002; Xiao et al. 2011

Note: GM, geometric mean; NA, not available.

^aThe urinary or blood concentrations of 23 chemicals or their metabolites were calculated based on the sample weighted maximum (Max)/GM concentrations from earlier biomonitoring studies; The blood here refers to whole blood, plasma or serum. In this study, 2GM and 5GM represent the 2-fold and 5-fold GM, and 5Max and 10Max represent the 5-fold and 10-fold maximum. Detailed information on the compound and concentration selection can be found in “Methods” section.

^bThe percentage is calculated based on the molar concentration of each compound.

^cUse metabolite as a replacement for the parent compound: 2,4-dihydroxybenzophenone (DHB) replaced its parent compound benzophenone-3 (BP-3).

^dUse parent compound as a replacement for the metabolite: tris(1,3-dichloroisopropyl)phosphate (TDCPP) replaced bis (1,3-dichloro-2-propyl) phosphate (BDCPP), triphenyl phosphate (TPhP) replaced diphenyl phosphate (DPhP), di-n-butyl phthalate (DnBP) replaced mono-3-hydroxybutyl phthalate (MHBP), benzylbutyl phthalate (BzBP) replaced mono-(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP), diisodecyl phthalate (DiDP) replaced mono-(carboxynonyl) phthalate (MCNP), and di-isononyl phthalate (DiNP) replaced mono-isononyl phthalate (MiNP).

^eThe geometric mean concentrations of triclocarban, butyl paraben, ethyl paraben and diisononyl phthalate (DiNP) were below detection limit and were replaced with the mean values of 75th percentile of samples.

Methods

Chemicals and Materials

Analytical standards of bisphenol A (BPA; ≥97% ), bisphenol S (BPS; ≥98% ), 10 mM Tris pH 8, and ammonium hydroxide were purchased from VWR. Bisphenol A diglycidyl ether (BADGE; ≥95% ); bisphenol A bis(2,3-dihyroxypropyl) ether (BADGE-2H2O ; ≥97% ); bisphenol F diglycidyl ether (BFDGE; ≥95% ); 2,4-dihyroxybenzophenone (DHB; ≥99% ); triphenyl phosphate (TPhP; ≥99% ); di-n-butyl phthalate (DnBP; ≥99% ); diisononyl phthalate (DiNP; ≥99% ), diisodecyl phthalate (DiDP; ≥99% ); genistein (≥98% ); perfluorooctanoic acid (PFOA; ≥95% ); tetrabromobisphenol A (TBBPA; ≥97% ); resazurin sodium salt (∼80% ); ammonium acetate (≥97% ); Triton X-100; and other analytical grade chemicals, including butyl paraben, ethyl paraben and methyl paraben, triclosan (TCS), triclocarban, tris(1,3-dichloro-2-propyl)phosphate (TDCPP), mono-n-butyl phthalate (MnBP), benzylbutyl phthalate (BzBP), and mono-benzyl phthalate (MBzP) were purchased from Sigma-Aldrich. Quant-iT PicoGreen dsDNA Reagent, EDTA pH 8.0, Tris-Cl pH 8.0, Pierce BCA protein assay kit, and dimethyl sulfoxide (DMSO; ≥99.5% ) were purchased from Thermo Fisher. Liquid chromatography–mass spectrometry (LC-MS) grade Acetonitrile were purchased from Aik Moh Paints & Chemicals. Electrospray ionization (ESI) tuning mix and reference mass solution kit were all purchased from Agilent Technologies. The 273 endogenous metabolite standards were purchased from Sigma-Aldrich (see Excel Table S1).

MCF-7 cells were obtained from the American Type Culture Collection (ATCC). Cells were cultured in humidified 37°C atmosphere containing 5% CO2 with Dulbecco’s Modified Eagle Medium (DMEM; 10% fetal bovine serum (FBS), L-glutamine, and sodium pyruvate; Gibco, Thermo Fisher Scientific).

Experimental Workflow

Mixtures consisting of 23 xenobiotic chemicals at several human-relevant levels from already available biomonitoring data were prepared to mimic actual human exposure in the United States (CDC 2018; Frederiksen et al. 2010; Högberg et al. 2008; Kim et al. 2015; Li et al. 2015, 2017, 2018a; Nagayama et al. 2000; Thomsen et al. 2002; Wang et al. 2015a; Xiao et al. 2011; Zhao et al. 2017). The MCF-7 cell was chosen as the experimental model because it is a sensitive hormone-responsive cell line (Levenson and Jordan 1997), and many selected chemicals have been previously reported to be xenoestrogens and EDCs. Figure 1 summarizes the experimental design and entire workflow of this study. First, omics technologies metabolomics and transcriptomics were employed to characterize the biological effect on cells of exposure to this mixture at several human-relevant levels. Second, the altered genes and metabolites were integrated to further interpret the biological changes triggered by the mixture exposure. Finally, to evaluate the relative contribution of each component to the observed effect, a high-throughput targeted metabolomics method (i.e., leave-one-out method; see below) was developed to screen the causal chemicals, by characterizing the metabolome changes of the mixtures with and without each tested compound.

[Image omitted - see PDF]

Selection of Xenobiotics

The primary aim of this study was to develop a platform for characterization of mixture effects and prioritization of the relative contribution of each component to the observed effect. Twenty-three compounds were selected and used at levels that mimic previously measured human exposure, though the individual exposure may be quite personalized and varied greatly. The xenobiotic mixture contained organic contaminants with diverse MOA, namely bisphenol plasticizers (n=5 ), phthalate plasticizers (n=7 ), PCPs (n=6 ), flame retardants (n=3 ), perfluorochemicals (n=1 ), and phytoestrogens (n=1 ). As mentioned above, the criteria for selecting these chemicals were as follows: high-volume manufacturing, endocrine-disrupting behaviors (e.g., bisphenols, parabens, flame retardants), showing toxicity, widespread occurrences in human or environmental samples, representing a wide-array of MOAs [e.g., endocrine activity (bisphenol A); lipid metabolism disturbance (triclosan); mitosis, cell cycle interference, and endocrine activity (genistein); peroxisome proliferator-activator receptors, thyroid hormone system, fatty acid homeostasis, and cell communication (perfluorooctanoic acid)], and existing large-cohort biomonitoring data. Here, the toxicity of parent or its metabolites is also another important criterion, and we preferred to select the one with higher or known toxicity. For example, mono-(2-ethylhexyl) phthalate (MEHP) is more potent than its parent bis(2-ethylhexyl) phthalate (DEHP) in many toxicological end points (Erkekoglu et al. 2010; Huber et al. 1996). In contrast, chemicals such as tris(1,3-dichloroisopropyl)phosphate (TDCPP) and triphenyl phosphate (TPhP) have been reported to be toxic in in vivo studies (McGee et al. 2012; Wang et al. 2015b), but there are very few reports on the toxicity of their metabolites. In these cases, we decided to choose metabolite MEHP and parent compounds TDCCP and TPhP for our test. For compounds like benzophenone-3 (BP-3; parent compound) with an extremely low half-life, we selected its major metabolite 2,4-dihydroxybezophenone (DHB). The half-life of BP-3 was predicted to be 1.8 h, based on the quantitative structure–activity relationship (QSAR) approach called Interactive Fragment Selection (IFS), which is a one-compartment pharmacokinetic mass balance model to estimate the primary biotransformation half-life (Arnot et al. 2014; Brown et al. 2012).

Selection of Human Exposure Source Data

For 16 of 23 compounds in this study, we referred to the urinary concentrations documented on the CDC Fourth National Report on Human National Report on Human Exposure to Environmental Chemicals (the Fourth Report) with cumulative biomonitoring data (1999–2016) with large sample size n>10,000 for most of chemicals (CDC 2018) (Table 1). The geometric mean (GM) values of these 16 compounds for urine samples (CDC 2018) were compared with the blood concentrations that were derived from existing biomonitoring data (calculated by sample weighted average if available) or predicted values using the steady-state model for those chemicals with no documented blood concentration (see Table S1). The calculation of the blood steady-state concentrations was based on their homologous series with similar structure by Equation 3 (Wadhwa and Cascella et al. 2020)

t1/2=0.693×VdCL [1]

Css=Dose×FID×CL. [2]

(1)+(2)

Css,a=DIa×t1/2,a×Vd,bDIb×t1/2,b×Vd,aCss,b [3]

DI=Cu×VuBW×Fue×MWdMWm [4]

where t1/2=Half-life was predicted based on the QSAR approach called IFS (Arnot et al. 2014; Brown et al. 2012); Vd=volume of distribution ; CI=clearance ; F=bioavailability ; ID=dosing interval ; DI=daily intake (Equation 4) (Guo et al. 2014); where Cu: Urinary concentration, Vu=2L (urine volume), BW: body weight, Fue : fractional urinary excretion, MWd: molecular weight of metabolites, MWm: molecular weight of parent compound; Vd was calculated based on the Risk Assessment Identification and Ranking-Indoor and Consumer Exposure (RAIDAR-ICE) model (Li et al. 2018b). The differences between urine and blood were all within the same order of magnitude range (0.27 to 21). The majority of urine-to-blood ratios were less than 5, which we judged to be within an acceptable range, considering the complexity of actual human exposure. Therefore, levels in both blood and urine were considered as human-relevant in this study. Another reason is the limited sample size and variation between each blood biomonitoring study. In general, there were less biomonitoring data (n<560 , except for genistein) for the blood concentrations, and the values reported in different studies also exhibited a large variation for many chemicals. For example, due to different dietary habits, the genistein blood mean concentrations of Japanese women (109.7 ng/mL ) were almost 70 times higher than those of Finnish women (1.67 ng/mL ) in the same study (Uehara et al. 2000). The genistein blood mean concentrations of people living in Asian countries such as Korea (211 ng/mL ) and Vietnam (289.5 ng/mL ) (Ko et al. 2018) were much higher than those of people living in United Kingdom (4.1 ng/mL ) (Grace et al. 2004). Last, the concentration of compounds in the blood or urine is just a snapshot of human exposure. The levels of the pollutants are expected to be dynamic and fluctuating, even for the same participant on the same day. There are also high uncertainties (e.g., different reference compounds, parent metabolite ratio, and different modeling method, etc.) of conversion using the steady state urine-blood conversion model. The simulated results often vary greatly with actual concentrations. After considering all these factors, the urinary concentrations were not converted to blood concentrations in this study.

Xenobiotic Mixture Selection

The GM, 95th percentiles, and maximum (Max) urinary or blood concentrations at human-relevant levels for the 23 chemicals are summarized in Table 1. Concentrations GM and Max were applied as the basis to reflect the common exposure scenario and the extreme exposure scenario, respectively. The concentrations of 17 chemicals or their metabolites (16 urine and 1 blood) were drawn from the Fourth Report. These biomonitoring data are derived from samples provided by participants of the National Health and Nutrition Examination Survey (NHANES) consisting of about 2,500 participants from United States with a wide age range (blood: 1 y and older; urine: 6 y and older) (CDC 2018). Chemical concentrations were calculated by the average of all detected Max/GM values from the year 1999 to 2016 (GM was replaced by 75th percentile concentration if below limit of quantification). The concentration of the remaining chemicals was calculated based on the sample size weighted Max/GM concentration from available biomonitoring studies (Frederiksen et al. 2010; Högberg et al. 2008; Kim et al. 2015; Li et al. 2015, 2017, 2018a; Nagayama et al. 2000; Thomsen et al. 2002; Wang et al. 2015a; Xiao et al. 2011; Zhao et al. 2017). The single-compound stock solution was prepared in DMSO initially, and the stocks of chemicals in the mixture were pooled to create a solution that was 1,000 times each chemical’s final concentration before dosing. Mixtures with varying concentrations [GM, 2GM (i.e., 2 times or 2-fold of GM levels), 5GM, Max, 5Max (i.e., 5 times or 5-fold of the Max levels) and 10Max] were independently mixed based on the biological concentrations of each chemical.

Cell Viability and DNA Quantification

Cell viability of MCF-7 cells after exposure to the mixture was evaluated by the resazurin assay (Fang et al. 2015c) with four replicates (n=4 ). In brief, cells were seeded in 96-well flat-bottomed plates (Thermo Fisher Scientific) with a seeding density of 30,000 cells/well in 100μL DMEM for 7 h. Subsequently, the 100μL medium dissolved with different levels of mixtures were added with final concentrations of GM, 2GM, 5GM, Max, 5Max, and 10Max in different wells. Nontreated negative control samples were prepared by adding DMSO with a final concentration of 0.1%. After exposure for 40 h, 20μL resazurin (0.50mg/mL ) was added into each well and further incubated for 8 h. The viability was measured by fluorescence using a 530 EX nm/590 EX nm filter setting. Cell viability was calculated by setting the viability of the negative DMSO control samples to 100%. DNA concentration of MCF-7 cells after exposure to the mixture was quantified by the PicoGreen Assay with three replicates (n=3 ). Cells were seeded with a density of 0.4×106 in 6-well plates. When about 60% confluency was attained, cells were exposed to varying levels of the mixture (GM, 2 GM, 5GM, Max, 5Max) for 40 h. After exposure, the cells were detached in a 1-mL lysis buffer (10 mM Tris pH 8.0, 1 mM EDTA and 0.2% (v/v) Triton X-100). Keeping on ice throughout, the samples were vortexed for 10 s every 5 min for half an hour. Samples were diluted by adding 30-μL samples to a 170-μL TE buffer (10 mM Tris-Cl pH 8.0, 1 mM EDTA pH 8.0). Next, 50μL of diluted samples and PicoGreen working reagent (diluted 200 times) were added to 96-well plates and further measured by fluorescence using the 485 EX nm/535 EX nm filter setting. DNA concentrations of the treatment groups were normalized by setting the concentration of negative control samples to 100%. One-way analysis of variance (ANOVA) and Duncan’s post hoc analysis were used to test whether one treatment group was statistically different from the nontreated control; an analysis result of p≤0.05 was considered statistically significant.

Mixture Exposure and Metabolite Extraction

The MCF-7 cells were seeded in 6-well plates with a density of 0.4×106 in each well. After reaching around 60% confluency, cells were exposed to xenobiotic mixtures (each level with four replicates, n=4 ) with varying levels (GM, 2GM, 5GM, Max, and 5Max) in the culture medium for 40 h. For single-chemical exposures, cells were exposed to Max level of DHB (1911.7 ng/mL ) or BPA (12.9 ng/mL ) in the culture medium for 40 h. Nontreated control samples were prepared by adding DMSO with a final concentration at 0.1%. The method for metabolite extraction was a modified version of that used in our previous studies (Beyer et al. 2018; Fang et al. 2015a). Cell metabolites were extracted with 1.6mL ice-cold methanol:acetonitrile:water (2:2:1, v/v/v) and then collected by a cell scraper. The intracellular metabolites were extracted by freeze–thaw steps conducted thrice using liquid nitrogen followed by sonication in an ice bath for 10 min. Afterward, the samples were placed at −40°C for 1 h, followed by 15-min centrifugation at a speed of 13,000r.p.m. at 4°C to precipitate proteins. The protein concentration of the cells was measured in the final pellet after centrifugation using Pierce BCA assay (n=3 ) according to kit protocols. A 100-μL diluted sample (100 times) was mixed with 100μL of working reagents in 96-well plates and incubated at 37°C for 4 h. The absorbance was measured at 562 nm on a plate reader. The resulting supernatant was removed and evaporated to dryness by a CentriVap centrifugal vacuum concentrator (Labconco). The dried extracts were then reconstituted in the appropriate volume of acetonitrile:water (1:1, v/v), normalized by the protein concentration (see Table S2) with the lowest concentration of 40μL , sonicated for 10 min, and centrifuged for 15 min at 13,000r.p.m. and 4°C to remove insoluble debris. The supernatants were transferred to high-performance liquid chromatography (HPLC) vials with inserts and stored at −40°C before analysis.

Metabolite Profiling and Quality Assurance/Quality Control (QA/QC)

Instrumental analysis was performed using a HPLC system (1200 series, Agilent Technologies) coupled to a 6550 iFunnel quadrupole time-of-flight (QTOF) mass spectrometer (Agilent Technologies) according to our previous studies (Beyer et al. 2018; Fang et al. 2015a). To have a better metabolite coverage, a Phenomenex Luna NH2 Column (particle size 3μm , i.d. 2×length 150mm , for polar metabolites) and Waters Atlantis T3 column (particle size 3μm , i.d. 2.1×length 100mm , for hydrophobic metabolites) were used for separation in the negative and positive modes, respectively. For NH2 column, the injection volume was 6μL , and the mobile phase was composed of A=20 mM ammonium acetate and 40 mM ammonium hydroxide in 95% water and 5% acetonitrile and B=95% acetonitrile and 5% water. The flow rate was 0.25mL/min and the column temperature was 30°C. The linear gradient was set as follows: 0 to 2 min : 100% B; 2 to 15 min : 100% B to 10% B; 15 to 17 min : 10% B to 0% B; 17 to 33 min : 0% B; 33 to 55 min : 0% B to 100% B. The T3 column method was slightly modified based on an earlier study (Wang et al. 2019) with an injection volume of 5μL . The mobile phase A was 0.1% formic acid in water, and the mobile phase B was 0.1% formic acid in acetonitrile. The linear gradient was set as follows: 0 to 1 min : 1% B; 1 to 8 min : 1% B to 100% B; 8 to 10 min : 100% B; 10 to 10.5 min : 100% B to 1% B; 10.5 to 12.5 min : 1% B. Flow rate was 0.5mL/min and column temperature was 25°C. Other parameters were set as follows: gas temperature: 225°C; gas flow: 14L/min ; sheath gas temperature: 275°C; sheath gas flow: 11L/min ; nozzle voltage: 1,000V ; nebulizer pressure: 35 psig; capillary voltage: 3,500V and −3,500V for positive (ESI+ , electrospray ionization) and negative mode (ESI− ); nozzle voltage: 1,000V ; fragmentor: 175V . The instrument was set to acquire data over the m/z range of 50 to∼1,200 , with the MS acquisition rate of 1 spectra/s. The system was tuned in extended dynamic range mode (2 GHz ). Mass calibration was enabled using reference masses of 112.9,856 and 1033.9,881 in ESI (− ) or 121.0508 and 922.0098 in ESI (+ ). For the MS/MS of selected precursors, the default isolation width was set as narrow (∼1.3 m/z ), with a MS acquisition rate at 2 spectra/s and MS/MS acquisition at 6.67 spectra/s to acquire over the m/z range of 50 to∼1,100 . Tandem mass spectrometry (MS/MS) data were acquired at the collision energy of 20V and 40V . The sample sequence was randomized to avoid systematic decreases in signals over sample sets. To correct the mass, retention time, and response drift, a mix of metabolites (QC sample) was prepared by pooling all treated and control cell samples. QC analysis was conducted once for every six injections of biological samples and a blank sample (acetonitrile:water, 1:1, v/v). The variation of four biological replicates for each treatment group was manually checked with relative standard error of most features being less than 15% to∼20% as detailed in Figure S1. Data-dependent acquisition (DDA) auto-MS/MS and targeted MS/MS of selected precursors were run with the QC sample to confirm the identity of metabolites. To compare the metabolic changes under different dosing levels, one-way ANOVA and Duncan’s post hoc analysis were used to test whether one treatment group was statistically different from nontreated controls; a finding of p≤0.05 was considered statistically significant.

Metabolite Identification and Metabolic Pathway Analysis

The metabolite profiling data were converted to mzXML files using Agilent MassHunter Acquisition Software 6.0 and were then uploaded to the cloud-based XCMS Online platform for data processing (Tautenhahn et al. 2012) including peak detection, retention time correction, profile alignment, and isotope annotation. Data were processed using both pairwise and multigroup comparison. The parameter settings were as follows: centWave for feature detection (Δm/z=15 ppm ; minimum peak width=10s ; and minimum peak width=120s ); OBI-Warp setting for retention time correction (profStep=1 ); and parameters for chromatogram alignment, including mzwid=0.015 , minfrac=0.5 , and bw=5 . The relative quantification of metabolite features was based on extracted ion chromatograph (EIC) areas. After alignment and annotation, the significant features were first screened with p≤0.05 , maximal intensity >10,000 and |fold change|>1.3 . Metabolite identification was carried out using a few filters, including MS/MS fragment match (auto/targeted MS/MS), accurate mass match, and in-house retention time match if available. Moreover, the metabolite standard verification by the identical HPLC method was applied to support the metabolite identification process if available. Metabolite features were searched against several open-source platforms, including METLIN (Guijas et al. 2018), KEGG (Kanehisa 2000), and HMDB (Wishart et al. 2007; Wishart et al. 2018). Pathway analysis was conducted using MetaboAnalyst (Chong et al. 2019). Graphs were made using GraphPad Prism (version 7). Statistical analysis was performed with SPSS statistics (version 22).

RNA Sequencing Preparation and Data Analysis

MCF-7 cells were exposed to the 23 chemical mixture at the concentration of Max for 40 h with three biological replicates (n=3 ) as described above (in the section on “Mixture Exposure and Metabolite Extraction”) and were further fixed by the RNAprotect reagent (Qiagen) immediately. The total RNA extraction, mRNA purification, library preparation, and library sequencing were carried out by a service company (Majorbio). For purification, the total RNA was extracted using TRIzol® Reagent according to the experimental protocol (Invitrogen) and genomic DNA was removed using DNase I (TaKara). The RNA quality was determined by 2100 Bioanalyzer (Agilent) and quantified using the ND-2000 (NanoDrop Technologies). Only high-quality RNA sample (OD260/280=1.8 to∼2.2 , OD260/230≥2.0 , RIN≥ 6.5 , 28S:18S≥1.0 , >2μg ) was used to construct sequencing library. Subsequently, RNA-seq transcriptome library was prepared following TruSeqTM RNA sample preparation Kit from Illumina using 1μg of total RNA. In brief, messenger RNA was isolated according to polyA selection method by oligo (dT) beads and then fragmented by fragmentation buffer. Subsequently, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen) with random hexamer primers (Illumina). Then the synthesized cDNA was subjected to end repair, phosphorylation, and “A” base addition according to Illumina’s library construction protocol. Libraries were size selected for cDNA target fragments of 200–300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 polymerase chain reaction (PCR) cycles. After quantified by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina NovaSeq 6000 (2×150 bp read length ). For read mapping, the raw paired-end reads were trimmed and quality controlled by SeqPrep (John 2011) and Sickle (Joshi and Fass 2011) with default parameters. Then clean reads were separately aligned to reference genome with orientation mode using TopHat software (Trapnell et al. 2009). The mapping criteria of bowtie was as follows: sequencing reads should be uniquely matched to the genome allowing up to 2 mismatches, without insertions or deletions. Then the region of gene was expanded following depths of sites and the operon was obtained. In addition, the whole genome was split into multiple 15 kb windows that share 5 kb . New transcribed regions were defined as more than 2 consecutive windows without overlapped region of gene, where at least 2 reads mapped per window in the same orientation.

Differential Expression Analysis, Functional Enrichment and Joint Pathway Analysis

To identify differential expression genes (DEGs) between two different groups, the expression level of each transcript was calculated according to the fragments per kilobase of exon per million mapped reads (FPKM) method. RSEM (Li and Dewey 2011) was used to quantify gene abundances. EdgeR (Empirical analysis of Digital Gene Expression in R) was used for differential expression analysis (Robinson et al. 2010). After gene annotation, filter criteria of | log2(Fold change)|≥0.68 (i.e., the fold change of 1.5 in the expression obtained by comparing treated samples with untreated control) and adjusted p≤0.05 were applied to find the significantly differentially expressed genes. Gene ontology (GO) analysis was conducted with the DAVID Bioinformatic Resources 6.8 (Dennis et al. 2003) with adjusted p≤0.05 , using the Benjamini-Hochberg method. Unsupervised pathways analysis was conducted using a web server, for gene functional annotation KOBAS 3.0 (Wu et al. 2006) with adjusted p≤0.05 using the Benjamini-Hochberg method, and enrichment with top three rankings for each database were listed for illustration. A supervised approach was also adopted to assess the gene enrichment of a few xenobiotic related pathways based on a previous study (Xia et al. 2017). The gene hit numbers were expressed as a percentage of the total gene numbers in the respective pathway. The p -value was calculated based on Fisher’s exact test (Fisher 1935) and corrected by the Holm-Bonferroni method (Holm 1979) using the built-in stats package version 3.6.2 of R (R Core Development Team) (hypergeometric and p .adjust function). The demo codes are provided in Text S1. The source code can be found at https://www.rdocumentation.org/packages/stats/versions/3.6.2. An adjusted p≤0.05 was considered significant. Joint pathway analysis was conducted using MetaboAnalyst by selecting 20 top −log  (p-value of the hypergeometric test) values and hit numbers less than 5 were excluded from the discussion (Chong et al. 2019).

Leave-One-Out Method

Considering that the effect triggered by an individual chemical at human-relevant level could be too low to be distinguishable from nontreated samples, a novel leave-one-out approach was used to evaluate the relative contribution of each chemical to the observed mixture effect. This method ensured stable observable molecular change by characterizing the metabolome of the mixtures before and after the removal of one component. Twenty-three different mixtures were prepared each one of which had one compound removed with the others remaining at Max level, for its comparably stable responses among the treatments (total 25 groups: including 23 treatment groups with one component removed, one 23 chemical mixture labeled as “All Compounds,” and one DMSO control group). Due to the large sample size (n>90 for each batch), a high-throughput metabolomics method was developed with the modification based on the above-mentioned global metabolomics method. The chemical exposure and extraction method were the same as mentioned above with three replicates (n=3 ). For instrumental analysis, an earlier method (Wang et al. 2019) was slightly modified using Waters ACQUITY UPLC® BEH amide column (particle size 1.7μm , i.d. 2.1×length 100mm ) in ESI− with an injection volume of 10μL , and the run time was reduced from 55 to 12 min per sample. The mobile phase A was 25 mM ammonium acetate and 25 mM ammonium hydroxide in water, and the mobile phase B was acetonitrile. The flow rate was 0.5mL/min , and the column temperature was 25°C. The linear gradient was set as follows: 0 to 0.5 min : 95% B; 0.5 to 7 min : 95% to 65% B; 7 to 8 min : 65% B to 40% B; 8 to 9 min : 40% B; 9 to 9.5 min : 40% B to 95% B; 9.5 to 12.5 min : 95% B. This new method can save several weeks in sample preparation and analysis. Within 24 significantly altered metabolites that were identified at Max level, 19 could be well detected with good peak shape using this high-throughput method. These metabolites were further used as targeted compounds, and their peaks were extracted and manually integrated using the Agilent QToF Quantitative analysis software. The data sets were corrected using a total abundance regression calibration method (Wang et al. 2013). Afterward, the peak intensities of all groups were analyzed by unsupervised hierarchical clustering analysis based on Ward’s clustering method with the Euclidean distance used as the similarity metric after standardization.

Results

Mixture Cell Viability and DNA Quantification

Among all tested concentrations, only 10 Max showed significantly lower cell viability (32% of control) when compared with the DMSO control (Figure 2A). The cell viability of 2GM, 5GM, and Max were 125%, 128%, and 119%, respectively, suggesting the possible proliferation of cells. No significant difference was observed between 5Max and the control samples. DNA concentration was quantified for cells exposed to different levels of the mixture and normalized by setting control samples as 100% (Figure 2B). The results were generally in agreement with the cell viability data. Although a lower DNA concentration was observed at 5Max level, DNA concentrations from GM to Max were comparably higher than control (GM: 124%; 2GM: 131%; 5GM: 134% and Max: 110%), further suggesting the proliferation after exposure to xenobiotic mixture.

[Image omitted - see PDF]

Metabolomic Profiling

Because it is difficult to identify cell metabolites completely, we first semiquantitatively evaluated the dose–response effect of the xenobiotic mixture (Figure 2C) by first screening the feature numbers with p≤0.05 , maximal intensity >10,000 ion counts , and |fold change|>1.3 , followed by checking all features manually to minimize false-positive hits. In general, higher mixture doses induced more responsive metabolite features. More up-regulated features (all more than 50% of the total significant features) were generated regardless of dose level and ESI mode. The distinct and shared significant features were identified between different dose levels in ESI negative (Figure 2D) and positive modes (Figure 2E). Overall, more distinct features were found at 5Max level, which accounted for 81% of its total features. In contrast, treatment at Max, 5GM, 2GM, and GM generated more overlapped features, which accounted for 58%, 60%, 64%, and 62% of their total features, respectively.

Dysregulated Metabolite Identification

In total, 54 significantly differentially detected metabolites were identified across all treatment conditions (see Table S3), including nucleoside phosphate compounds, lipids, carnitines, and amino acids. As seen in Figure 3A, 5Max treatment showed higher potency to result in differentially detected metabolites that were often down-regulated, whereas an increasing number of upregulated metabolites were observed from concentrations of GM to Max. After 5Max treatment, nearly all nucleoside phosphate compounds were down-regulated, whereas several lipids [e.g., PC (16:0/0:0), PC (O-15:0/2:0)] and fatty acids (e.g., stearic acid and palmitic acid) were up-regulated. For levels not higher than Max, the dysregulated metabolites were mainly nucleoside phosphate compounds, such as adenosine triphosphate (ATP), adenosine diphosphate (ADP), and cytidine-5′-triphosphate (CTP), as well as carnitines such as L-palmitoylcarnitine, octanoylcarnitine, and L-carnitine. Approximately 40% of detected metabolites exhibited a significant dose–response relationship from the levels of GM to Max. For example, the responses of uridine diphosphate (UDP), cytidine diphosphate (CDP), and CTP were significantly higher as the chemical concentration increased from GM to Max (Figure 3B). Cellular lactate concentrations were significantly higher at Max and 5Max (Figure 3B) when compared with control. Overall, many metabolites exhibited a significant dose–response relationship after the mixture exposure.

[Image omitted - see PDF]

Metabolomic Pathway Analysis

All the significantly differentially detected metabolites identified for Max (Figure 4) and 5Max (see Figure S2A) treatments were mapped onto metabolic pathways. Metabolic responses after the Max mixture exposure were primarily involved in “purine metabolism,” “pyrimidine metabolism,” and “ascorbate and aldarate metabolism.” Especially in pyrimidine (Figure 4B) and purine metabolisms (Figure 4C), the majority of metabolites [uridine triphosphate (UTP), guanosine diphosphate(GDP), guanosine triphosphate (GTP), CDP, CTP, ATP, and ADP] at GM, 2GM, and 5GM levels were consistent with Max, all showing an up-regulated trend, whereas 5Max often exhibited down-regulation or no responses. After 5Max treatment, the potential influenced pathways were “nitrogen metabolism,” “glycerophospholipid metabolism,” and “glutathione metabolism” (see Figure S2B and S2C).

[Image omitted - see PDF]

Transcriptomic Profiling

HiSeq RNA sequencing was applied to profile the differentially expressed genes (DEGs) of MCF-7 cells exposed to the mixture at Max level. The heatmap showed the top 80 significantly differentially expressed genes (see Figure S3). In total, 1,619 significantly differentially expressed genes were identified, and 45% of them showed up-regulation (see Table S4). The volcano plot visualized these DEGs in response to the mixture exposure (see Figure S4). Among these genes, the maximally up-regulated ones include the imprinted maternally expressed transcript [H19, log2(FC)=2.39 and the early growth response 3 (EGR3, 2.32) as well as the down-regulated genes, the C-Type lectin domain family 3 member A (CLEC3A, −3.19 ), and the membrane metalloendopeptidase (MME, −3.06 ).

GO Analysis and Unsupervised Pathway Analysis

To better understand the functions and pathways of these DEGs, we performed unsupervised GO functional annotation and enrichment analysis (Figure 5A) as well as unsupervised pathway analysis (Figure 5B). Most DEGs were involved in the biological processes such as cell proliferation, cell division, and DNA replication. The top three enriched GO terms in cellular component were “cytoplasm,” “extracellular exosome,” and “membrane,” and two additional GO terms in molecular function were “protein binding” and “structural constituent of cytoskeleton.” By mapping these DEGs to the pathways from databases KEGG, PANTHER, and REACTOME (Figure 5B), the top enriched pathways across all databases were “axon guidance mediated by Slit/Robo,” “DNA replication,” and cell cycle.

[Image omitted - see PDF]

Supervised Pathway Analysis

A supervised approach was also adopted by mapping DEGs onto some cellular toxicity pathways (Xia et al. 2017) (Figure 5C). Overall, most of the responses were found at the end points of genotoxicity followed by adaptive stress response and endocrine disruption. Genotoxicity, nuclear factor erythroid 2-related factor 2 (Nrf2), tumor protein p53 (TP53), glucocorticoid receptor (GR) (with p≤0.01 ) and retinoic acid receptor-related orphan receptor (ROR), thyroid hormone receptor alpha (THRα ), and androgen receptor (AR) (with p≤0.05 ) were identified as significantly affected. For genotoxicity, 29 genes (25% of total genes in this pathway) were mapped onto this pathway. Some important replication proteins were all significantly up-regulated, such as minichromosome maintenance protein complex [MCM2: log2(FC)=0.81 ; MCM5: 0.7; MCM7: 0.76; MCM10: 0.98], and origin recognition complex (ORC6: 0.65; ORC1: 1.24). Regarding adaptive stress response, higher responses were found in Nrf2 (21 hits, 15%) and TP53 (5 hits, 29%). For Nrf2 pathway, genes of the solute carrier family (SLCA11: 1.17) and ATP-binding cassette transporter (ABC transporter: ABCC5, −0.73 ), functioning as drug-metabolizing enzymes, were significantly altered. Last, higher responses were found in GR (10 hits, 13%) among hormone receptors.

Joint Pathway Analysis

Joint pathway analysis was conducted to integrate the metabolite and gene data of cells in response to exposure to the mixture at Max. Overall, pathways such as pyrimidine metabolism, “arginine and proline metabolism” and purine metabolism were found with comparably higher hits (Figure 6A). Among these pathways, both gene and metabolite data could be well aligned and mapped onto pyrimidine (Figure 6B) and purine metabolism (Figure 6C), which is in line with the above-mentioned analysis.

[Image omitted - see PDF]

Relative Contribution of Each Chemical to the Mixture Effect

We performed the high-throughput metabolomics analysis for the cells using the leave-one-out method to prioritize the relative contribution of each chemical. The unsupervised hierarchical clustering analysis grouped similar objects among the targeted 19 metabolites across 25 treatment groups (Figure 7A). Overall, the metabolome patterns from removing DHB or BPA showed great resemblance with the nontreated control (DMSO), suggesting their relatively higher contribution among the chemicals in this mixture. The metabolome patterns of removing some chemicals such as bisphenol A bis(2,3-dihydroxypropyl) ether (BADGE-2H2O ) and bisphenol A diglycidyl ether (BADGE), showed great similarity with the “All Compounds” treatment, indicating their minimal contribution to the overall metabolome changes.

[Image omitted - see PDF]

The combined EIC showed the significantly differentially detected 19 metabolites after exposure to the xenobiotic mixture at Max level in ESI negative mode (see Figure S5). Several representative metabolites (Figures 7B, 7C, and 7D) were selected to illustrate whether chemicals in the same category exhibit similar behavior. For example, removing chemicals from several categories such as PCPs, bisphenols (BPA, BPS), and most of the phthalates (DnBP, MnBP, BzBP, MEHP) showed a significantly lower UDP level when compared with “All Compounds” treatment (Figure 7B). This finding suggests that these chemicals are all likely to have an important contribution to UDP production. After the removal of MnBP, BzBP, genistein, or BPA (Figure 7C), the cellular lactate level was comparably lower than “All Compounds” though not reaching the level of control yet, implying that these chemicals are likely to be responsible for excess lactate production. After removing most of the chemicals (Figure 7D), the cellular ATP levels were much lower, but there was no statistical evidence to suggest that any of them is primarily responsible for ATP dysregulation. We further characterized the changes of 19 metabolites for individual chemical DHB (1911.7 ng/mL ) and BPA (12.9 ng/mL ) at Max level (see Figure S6). In general, most metabolites did not exhibit significant differential change by exposure to single chemical, suggesting DHB or BPA alone was unlikely to initiate much metabolite dysregulation to the observed mixture effect. Especially, only the significant up-regulation of UDP-glucuronic acid was observed in the DHB treatment, though many metabolites were comparably higher than control; yet the trend was statistically insignificant. For BPA treatment, no significant change was observed for the 19 metabolites.

Discussion

In the present study, we employed metabolomic and transcriptomic profiling to reveal the impact of xenobiotic mixture exposure at a human-relevant level on the MCF-7 cellular biological processes and prioritized the relative contribution of each component to the observed effect. Considering the heavy workload of the omics characterization and the leave-one-out method (n>90 in a batch experiment for 23 compounds), we did not include more chemicals in this study. The selected chemicals with a wide array of MOAs in toxicity pose environmental concerns. The data suggested that the xenobiotic mixture significantly altered the cellular biological process of the MCF-7 cells. Many of the dysregulated metabolites or genes were strongly associated with the cell proliferation process and oxidative stress. A leave-one-out method was also developed using high-throughput metabolomics to evaluate the relative contribution of each component to the mixture effect. The results suggested that DHB and BPA had a relatively higher contribution to the observed mixture effect. These findings provided insights regarding cell perturbations when exposed to xenobiotic mixtures at a human-relevant level and offered a new proof of concept to decipher the low-dose mixture effect.

In metabolomics analysis, the differentially detected metabolites were associated primarily with cell proliferation and oxidative stress compared to controls, for doses not higher than Max. Taking Max as an example, some key metabolite patterns, such as generating ATP rapidly (log2(FC)=1.09 at Max ) and lactate production (0.79) (see Figure S7), are consistent with an increase in metabolism as is seen in proliferating cells (DeBerardinis et al. 2008). In addition, the trends in metabolite expression related to pyrimidine and purine metabolisms after exposure to the mixture in this study are similar to those seen in proliferating cells (Lee et al. 2017). In the pyrimidine pathway (Figure 4B), when compared with controls, downstream metabolites (UDP: 0.52; UTP: 1.8; CDP: 1.27 and CTP: 1.51) were up-regulated whereas the upstream metabolite uridine was down-regulated (−0.45 ). A similar and significant trend was observed for purine metabolism (Figure 4C). Pyrimidine and purine nucleotides are major energy carriers that play critical roles in DNA and RNA synthesis as well as membrane lipid biosynthesis and protein glycosylation (Moffatt and Ashihara 2002). Their metabolism is important for regulating cell cycle progression and maintaining the survival of the cells (Quéméneur et al. 2003; Wang et al. 2016). Regarding oxidative stress biomarkers, we observed the trend of higher levels of oxidized glutathione (GSSG, fold change (FC)>1.3 except for 5Max) and lower levels of glutathione (GSH) after mixture exposure at most doses. This alteration and comparably lower cellular GSH:GSSG ratio compared with controls suggested elevated oxidative stress. Additionally, some pathways related to differentially detected metabolites were also related to antioxidation activities such as ascorbate and aldarate metabolism (Max) and glutathione metabolism (5Max).

Transcriptome alterations induced by xenobiotic mixture were in line with metabolome findings suggesting cell proliferation after mixture exposure. Above all, the maximally up-regulated genes such as H19 and EGR3 were associated with the cell proliferation process (Berteaux et al. 2005; Xi and Kersh 2004). H19 transcription was up-regulated during the S-phase of growth-stimulated cells to promote cell cycle progression of breast cancer cells (Berteaux et al. 2005). ERG3-deficient mice had reduced thymocyte proliferation in response to pre–T cell receptor signals, revealing its role in promoting proliferation (Xi and Kersh 2004). Top biological process GO terms cell proliferation, cell division, and DNA replication agreed substantially with the metabolome results. Studies also reported the high affinity to protein binding sites of several selected chemicals, such as BADGE and phthalates (Petersen et al. 2008; Yue et al. 2014). These environmental xenobiotics can disrupt endocrine function by interfering with the binding of natural ligands to steroid receptors and binding proteins (Danzo 1997; Ishihara et al. 2003). Therefore, it is not surprising to find a high gene percentage (48%) in the molecular function term protein binding, raising concerns that these xenobiotics may disrupt or alter some vital function through moving or transporting to the site of action. For unsupervised pathway analysis, several studies have shown that axon guidance molecules (AGMs) and SLIT2/ROBO1 signaling is important in regulating cancer cell proliferation (Dallol et al. 2002; Harburg and Hinck 2011; Marlow et al. 2008; Strickland et al. 2006). SLIT2 overexpression or treating with SLIT2 conditioned medium reduced the proliferation of MCF-7 cells (Dallol et al. 2002; Marlow et al. 2008). In our study, the down-regulation of this gene (−1.03 ) suggested the proliferation of cells. Additionally, the E2F transcription factor family (E2F1:0.83; E2F2: 0.68; E2F6: 0.61; E2F7:0.83; E2F8: 0.65) showed a consistent up-regulated trend in cell cycle pathway. This family plays a crucial and well-established role in cell cycle progression that integrates with DNA repair and replication, chromatin assembly, condensation, and segregation (Dyson 1998; Nahle et al. 2002; Ren et al. 2002).

The supervised pathway analysis suggested that the enhancement of several biological activities after exposure to xenobiotic mixture are associated with genotoxicity and oxidative insults. For genotoxicity, earlier studies have reported that the minichromosome maintenance protein (MCM) family is implicated at the initiation step of DNA synthesis (Kearsey and Labib 1998; Tye 1999) and that the origin recognition complex (ORC) is also essential for the initiation of DNA replication (Bell and Stillman 1992; Bell et al. 1993). The differential expression of genes related to these families suggests that the xenobiotic mixture changed the DNA synthesis and replication process. The regulation of SLC and ABC transporter genes in the Nrf2 pathway, controls the expression of proteins critical in detoxification and elimination of reactive oxygen species (ROS) and electrophiles (Nguyen et al. 2009; Suzuki et al. 2013). The higher expression of genes in the Nrf2 pathway after treatment with the mixture suggests the generation of chemical or oxidative insults after mixture exposure, which is in line with the metabolomics analysis. Higher expression of members of Nrf2 pathway, such as elicited by the xenobiotic mixture, have also been associated with the promotion of tumorigenesis in fibroblasts (DeNicola et al. 2011). Last, because GRs regulate genes that control growth, metabolism, and immune response of an organism (Zhou and Cidlowski 2005), higher expression of these genes suggests that the mixture exposure may affect the anti-inflammatory and immunosuppressive actions of glucocorticoids and result in a failure to maintain homeostasis (Oakley and Cidlowski 2013).

Joint pathway analysis was consistent with omics results that mixture exposure at Max resulted in higher levels of genes and metabolites involved in pyrimidine and purine metabolism. The activation of pyrimidine may be initiated by the up-regulation of thymidine kinases (e.g., TK1: 0.62), which play an important role in pyrimidine nucleoside analogues activation (Al-Madhoun et al. 2004). The reason for proliferation after mixture exposure can be explained by the fact that MCF-7 is a hormone-responsive breast cancer cell line, and many chemicals in the mixture are xenoestrogens (e.g., genistein, parabens, phthalates, etc.). Previous studies have documented the ability of xenoestrogens to promote cell proliferation (Blom et al. 1998; Brown and Lamartiniere 1995). In addition, these xenobiotics can induce chemical insults/oxidative stress on cells through various mechanisms (Samet et al. 2018). Tumor cells can enhance glycolysis to provide adequate ATP for rapid proliferation, as a protection from oxidative stress (Brand and Hermfisse 1997). Therefore, we assume that the cells exposed to this mixture stimulated the production of nucleotides to accommodate an increasing need for RNA and DNA synthesis after chemical insults and oxidative stress; the possible mechanism for this is proposed and provided in Figure S7.

The relative contribution assessment reveals that DHB and BPA made relatively higher contributions to the overall metabolome changes, while each chemical may have a different effect on various metabolite regulations. The overall differential detection of certain metabolites related to purine and pyrimidine metabolism were generally quite similar among several chemicals. Taking UDP as an example, its concentrations were comparatively lower than “All Compounds” after removal of ethyl paraben, MnBP, BzBP, BPA, or DHB, indicating that these chemicals may have a substantial impact on the disruption of pyrimidine metabolism pathway. The explanation for these results could be that the doses of these chemicals were considerably higher (all >100 ng/mL except BPA: 12.9 ng/mL ) and studies have documented their ability to exert estrogenic effects or induce proliferation in human breast cancer cells (Arbuckle et al. 2016; Chen and Chien 2014; Engeli et al. 2017; Le Fol et al. 2017; Suzuki et al. 2005). However, it should be noted that some important chemicals such as DHB and BPA at maximal concentration alone failed to initiate the observed effects by themselves, implying the possibility of strong chemical interaction between those chemicals. Overall, the leave-one-out concept along with high-throughput metabolomics provided a useful and efficient tool to screen causal chemicals in a low-dose mixture. Compared with traditional methods (e.g., to examine the effect of each component first and then predict the mixture effect), this approach is more suitable for a human or environmentally relevant level mixture considering its observable molecular change during all experimental manipulations. In addition, taking the relative high cost of other omics tools such as sequencing and proteomics into account, the high-throughput analysis of metabolome has shown its superiority in the study of mixtures and their effects upon exposure at human-relevant levels.

Overall, we have demonstrated the cellular metabolome and transcriptome changes upon exposure to a mixture of 23 organic contaminants at human-relevant levels and prioritized the relative contribution of each chemical to this mixture effect. The omics results suggest that the cellular activity after xenobiotic mixture exposure was associated with cell proliferation and oxidative stress. In addition, with the high-throughput metabolomics-based leave-one-out method, different effects on metabolic pathway perturbation were triggered by different chemicals in this mixture, among which DHB and BPA showed comparably higher contributions. Considering that the current risk assessment usually focuses on individual chemicals, it is of great significance to evaluate influences triggered by exposure to mixtures with higher complexity at human-relevant levels. In addition, a proof of concept was constructed with the leave-one-out method to evaluate the relative contribution of each chemical, though the exact mechanism of interactions among the components remains unclear. The interaction between different metabolic pathways can be further investigated in the future. As a case study, only 23 chemicals were selected; thus, it is of great importance to further study a more complex mixture covering both inorganic and organic contaminants (e.g., permanent organic contaminants) that can mimic the actual human exposome. It should also be noted that the applied metabolomics method is not adequate to extract hydrophobic lipids, and more omics techniques such as lipidomics and proteomics can be incorporated in further studies. In this study, the use of a single cancer cell line is not adequate to reflect the physiological or pathological responses of the normal cells to environmental exposure. Therefore, more human-relevant models such as in vivo animals and ex vivo tissues (e.g., stem cell differentiated organs) can be considered for future mixture studies.

Acknowledgments

The authors would like to thank Editage ( www.editage.cn) for English language editing.

This work was funded by NTU-Harvard SusNano (M4082370.030) and Singapore Ministry of Health’s National Medical Research Council under its Clinician-Scientist Individual Research Grant (CS-IRG) (MOH-000,141) and Open Fund−Individual Research Grant (OFIRG/0076/2018).

References

Al-Madhoun AS, Tjarks W, Eriksson S . 2004. The role of thymidine kinases in the activation of pyrimidine nucleoside analogues. Mini Rev Med Chem 4(4):341–350, PMID:15134537,10.2174/1389557043403963 . Crossref, Medline,

Arbuckle TE, Fisher M, MacPherson S, Lang C, Provencher G, LeBlanc A, et al. 2016. Maternal and early life exposure to phthalates: the plastics and personal-care products use in pregnancy (P4) study. Sci Tot Environ 551-552:344–356,10.1016/j.scitotenv.2016.02.022 . Crossref, Medline,

Arnot JA, Brown TN, Wania F . 2014. Estimating screening-level organic chemical half-lives in humans. Environ Sci Technol 48(1):723–730, PMID:24298879,10.1021/es4029414 . Crossref, Medline,

Asimakopoulos AG, Thomaidis NS, Kannan K . 2014. Widespread occurrence of bisphenol a diglycidyl ethers, p-hydroxybenzoic acid esters (parabens), benzophenone type-UV filters, triclosan, and triclocarban in human urine from Athens, Greece. Sci Tot Environ 470-471:1243–1249,10.1016/j.scitotenv.2013.10.089 . Crossref, Medline,

Bell SP, Kobayashi R, Stillman B . 1993. Yeast origin recognition complex functions in transcription silencing and DNA replication. Science 262(5141):1844–1849, PMID:8266072,10.1126/science.8266072 . Crossref, Medline,

Bell SP, Stillman B . 1992. ATP-dependent recognition of eukaryotic origins of DNA replication by a multiprotein complex. Nature 357(6374):128–134, PMID:1579162,10.1038/357128a0 . Crossref, Medline,

Berteaux N, Lottin S, Monté D, Pinte S, Quatannens B, Coll J, et al. 2005. H19 mRNA-like noncoding RNA promotes breast cancer cell proliferation through positive control by E2F1. J Biol Chem 280(33):29625–29636, PMID:15985428,10.1074/jbc.M504033200 . Crossref, Medline,

Beyer BA, Fang M, Sadrian B, Montenegro-Burke JR, Plaisted WC, Kok BP, et al. 2018. Metabolomics-based discovery of a metabolite that enhances oligodendrocyte maturation. Nat Chem Biol 14(1):22–28, PMID:29131145,10.1038/nchembio.2517 . Crossref, Medline,

Blom A, Ekman E, Johannisson A, Norrgren L, Pesonen M . 1998. Effects of xenoestrogenic environmental pollutants on the proliferation of a human breast cancer cell line (MCF-7). Arch Environ Contam Toxicol 34(3):306–310, PMID:9504980,10.1007/s002449900322 . Crossref, Medline,

Boxall AB, Rudd MA, Brooks BW, Caldwell DJ, Choi K, Hickmann S, et al. 2012. Pharmaceuticals and personal care products in the environment: what are the big questions?Environ Health Perspect 120(9):1221–1229, PMID:22647657,10.1289/ehp.1104477 . Link,

Brand KA, Hermfisse U . 1997. Aerobic glycolysis by proliferating cells: a protective strategy against reactive oxygen species. FASEB J 11(5):388–395, PMID:9141507,10.1096/fasebj.11.5.9141507 . Crossref, Medline,

Brian JV, Harris CA, Scholze M, Kortenkamp A, Booy P, Lamoree M, et al. 2007. Evidence of estrogenic mixture effects on the reproductive performance of fish. Environ Sci Technol 41(1):337–344, PMID:17265968,10.1021/es0617439 . Crossref, Medline,

Brown TN, Arnot JA, Wania F . 2012. Iterative fragment selection: a group contribution approach to predicting fish biotransformation half-lives. Environ Sci Technol 46(15):8253–8260, PMID:22779755,10.1021/es301182a . Crossref, Medline,

Brown NM, Lamartiniere CA . 1995. Xenoestrogens alter mammary gland differentiation and cell proliferation in the rat. Environ Health Perspect 103(7–8):708–713, PMID:7588483,10.1289/ehp.95103708 . Link,

Carpenter DO, Arcaro K, Spink DC . 2002. Understanding the human health effects of chemical mixtures. Environ Health Perspect 110(suppl 1):25–42,10.1289/ehp.02110s125 . Link,

CDC (Centers for Disease Control). 2018. Fourth national report on human exposure to environmental chemicals. Atlanta, GA: Department of Health and Human Services.

Chen F-P, Chien M-H . 2014. Lower concentrations of phthalates induce proliferation in human breast cancer cells. Climacteric 17(4):377–384, PMID:24228746,10.3109/13697137.2013.865720 . Crossref, Medline,

Chong J, Wishart DS, Xia J . 2019. Using MetaboAnalyst 4.0 for comprehensive and integrative metabolomics data analysis. Curr Protoc Bioinformatics 68(1):e86, PMID:31756036,10.1002/cpbi.86 . Crossref, Medline,

Dallol A, Da Silva NF, Viacava P, Minna JD, Bieche I, Maher ER, et al. 2002. SLIT2, a human homologue of the Drosophila SLIT2 gene, has tumor suppressor activity and is frequently inactivated in lung and breast cancers. Cancer Res 62(20):5874–5880, PMID:12384551 , https://cancerres.aacrjournals.org/content/62/20/5874 [accessed30 August 2019 ]. Medline,

Danzo BJ . 1997. Environmental xenobiotics may disrupt normal endocrine function by interfering with the binding of physiological ligands to steroid receptors and binding proteins. Environ Health Perspect 105(3):294–301, PMID:9171990,10.1289/ehp.97105294 . Link,

DeBerardinis RJ, Lum JJ, Hatzivassiliou G, Thompson CB . 2008. The biology of cancer: metabolic reprogramming fuels cell growth and proliferation. Cell Metab 7(1):11–20, PMID:18177721,10.1016/j.cmet.2007.10.002 . Crossref, Medline,

DeNicola GM, Karreth FA, Humpton TJ, Gopinathan A, Wei C, Frese K, et al. 2011. Oncogene-induced Nrf2 transcription promotes ROS detoxification and tumorigenesis. Nature 475(7354):106–109, PMID:21734707,10.1038/nature10189 . Crossref, Medline,

Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. 2003. David: database for annotation, visualization, and integrated discovery. Genome Biol 4(9):R60,10.1186/gb-2003-4-9-r60 . Crossref,

Dettmer K, Hammock BD . 2004. Metabolomics–a new exciting field within the “omics” sciences. Environ Health Perspect 112(7):A396–A397, PMID:15159211,10.1289/ehp.112-1241997 . Link,

Dong T, Zhang Y, Jia S, Shang H, Fang W, Chen D, et al. 2019. Human indoor exposome of chemicals in dust and risk prioritization using EPA’s ToxCast database. Environ Sci Technol 53(12):7045–7054. PMID: 31081622,10.1021/acs.est.9b00280 . Crossref, Medline,

Dyson N . 1998. The regulation of E2F by pRB-family proteins. Genes Dev 12(15):2245–2262, PMID:9694791,10.1101/gad.12.15.2245 . Crossref, Medline,

Engeli R, Rohrer S, Vuorinen A, Herdlinger S, Kaserer T, Leugger S, et al. 2017. Interference of paraben compounds with estrogen metabolism by inhibition of 17β-hydroxysteroid dehydrogenases. Int J Mol Sci 18(9):2007,10.3390/ijms18092007 . Crossref,

Erkekoglu P, Rachidi W, Yuzugullu OG, Giray B, Favier A, Ozturk M, et al. 2010. Evaluation of cytotoxicity and oxidative DNA damaging effects of di (2-ethylhexyl)-phthalate (DEHP) and Mono (2-ethylhexyl)-phthalate (MEHP) on MA-10 Leydig cells and protection by selenium. Toxicol Appl Pharmacol 248(1):52–62, PMID:20659492,10.1016/j.taap.2010.07.016 . Crossref, Medline,

Fang M, Getzinger GJ, Cooper EM, Clark BW, Garner LV, Di Giulio RT, et al. 2014. Effect‐directed analysis of Elizabeth River porewater: developmental toxicity in zebrafish (Danio rerio). Environ Toxicol Chem 33(12):2767–2774, PMID:25196082,10.1002/etc.2738 . Crossref, Medline,

Fang M, Ivanisevic J, Benton HP, Johnson CH, Patti GJ, Hoang LT, et al. 2015a. Thermal degradation of small molecules: a global metabolomic investigation. Anal Chem 87(21):10935–10941, PMID:26434689,10.1021/acs.analchem.5b03003 . Crossref, Medline,

Fang M, Webster TF, Stapleton HM . 2015b. Effect-directed analysis of human peroxisome proliferator-activated nuclear receptors (PPARγ1) ligands in indoor dust. Environ Sci Technol 49(16):10065–10073, PMID:26172369,10.1021/acs.est.5b01524 . Crossref, Medline,

Fang M, Webster TF, Stapleton HM . 2015c. Activation of human peroxisome proliferator-activated nuclear receptors (PPARγ1) by semi-volatile compounds (SVOCs) and chemical mixtures in indoor dust. Environ Sci Technol 49(16):10057–10064,10.1021/acs.est.5b01523 . Crossref, Medline,

Fisher RA . 1935. The logic of inductive inference. J Royal Stat Society 98(1):39–82,10.2307/2342435 . Crossref,

Frederiksen H, Jørgensen N, Andersson AM . 2010. Correlations between phthalate metabolites in urine, serum, and seminal plasma from young Danish men determined by isotope dilution liquid chromatography tandem mass spectrometry. J Anal Toxicol 34(7):400–410, PMID:20822678,10.1093/jat/34.7.400 . Crossref, Medline,

Grace PB, Taylor JI, Low Y-L, Luben RN, Mulligan AA, Botting NP, et al. 2004. Phytoestrogen concentrations in serum and spot urine as biomarkers for dietary phytoestrogen intake and their relation to breast cancer risk in European Prospective Investigation of Cancer and Nutrition-Norfolk. Cancer Epidemiol Biomarkers Prev 13:698–708, PMID:15159299 , https://cebp.aacrjournals.org/content/13/5/698. Medline,

Guijas C, Montenegro-Burke JR, Domingo-Almenara X, Palermo A, Warth B, Hermann G, et al. 2018. METLIN: a technology platform for identifying knowns and unknowns. Anal Chem 90(5):3156–3164, PMID:29381867,10.1021/acs.analchem.7b04424 . Crossref, Medline,

Guo Y, Kannan K . 2013. A survey of phthalates and parabens in personal care products from the United States and its implications for human exposure. Environ Sci Technol 47(24):14442–14449, PMID:24261694,10.1021/es4042034 . Crossref, Medline,

Guo Y, Weck J, Sundaram R, Goldstone AE, Buck Louis G, Kannan K . 2014. Urinary concentrations of phthalates in couples planning pregnancy and its association with 8-hydroxy-2′-deoxyguanosine, a biomarker of oxidative stress: longitudinal investigation of fertility and the environment study. Environ Sci Technol 48(16):9804–9811, PMID:25068827,10.1021/es5024898 . Crossref, Medline,

Harburg GC, Hinck L . 2011. Navigating breast cancer: axon guidance molecules as breast cancer tumor suppressors and oncogenes. J Mammary Gland Biol Neoplasia 16(3):257–270, PMID:21818544,10.1007/s10911-011-9225-1 . Crossref, Medline,

Högberg J, Hanberg A, Berglund M, Skerfving S, Remberger M, Calafat AM, et al. 2008. Phthalate diesters and their metabolites in human breast milk, blood or serum, and urine as biomarkers of exposure in vulnerable populations. Environ Health Perspect 116(3):334–339, PMID:18335100,10.1289/ehp.10788 . Link,

Holm S . 1979. A simple sequentially rejective multiple test procedure. Scand J Stat: 65–70. http://www.jstor.org/stable/4615733 [accessed20 October 2020 ].

Huber WW, Grasl-Kraupp B, Schulte-Hermann R . 1996. Hepatocarcinogenic potential of di (2-ethylhexyl) phthalate in rodents and its implications on human risk. Crit Rev Toxicol 26(4):365–481, PMID:8817083,10.3109/10408449609048302 . Crossref, Medline,

Ishihara A, Sawatsubashi S, Yamauchi K . 2003. Endocrine disrupting chemicals: interference of thyroid hormone binding to transthyretins and to thyroid hormone receptors. Mol Cell Endocrinol 199(1–2):105–117, PMID:12581883,10.1016/s0303-7207(02)00302-7 . Crossref, Medline,

John JS . 2011. SeqPrep: tool for stripping adaptors and/or merging paired reads with overlap into single reads. https://githubcom/jstjohn/SeqPrep.

Joshi N, Fass J . 2011. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files (version 1.33) [software], https://github.com/najoshi/sickle.

Kanehisa M . 2000. Post-Genome Informatics. Oxford, UK: Oxford UP.

Kearsey SE, Labib K . 1998. MCM proteins: evolution, properties, and role in DNA replication. Biochim Biophysica Acta 1398(2):113–136,10.1016/S0167-4781(98)00033-5 . Crossref, Medline,

Kim SI, Yang YJ, Hong YP, Myung SC, Kim SC . 2015. Distribution of serum bisphenol a diglycidyl ether and its metabolite in Korean adult men and its association with reproductive hormone levels. Mol Cell Toxicol 11(1):71–78,10.1007/s13273-015-0009-3 . Crossref,

Ko K-P, Yeo Y, Yoon J-H, Kim C-S, Tokudome S, Koriyama C, et al. 2018. Plasma phytoestrogens concentration and risk of colorectal cancer in two different Asian populations. Clin Nutr 37(5):1675–1682, PMID:28778370,10.1016/j.clnu.2017.07.014 . Crossref, Medline,

Kortenkamp A . 2007. Ten years of mixing cocktails: a review of combination effects of endocrine-disrupting chemicals. Environ Health Perspect 115(suppl 1):98–105,10.1289/ehp.9357 . Link,

Le Fol V, Aït-Aïssa S, Sonavane M, Porcher J-M, Balaguer P, Cravedi J-P, et al. 2017. In vitro and in vivo estrogenic activity of BPA, BPF and BPS in zebrafish-specific assays. Ecotoxicol Environ Saf 142:150–156, PMID:28407500,10.1016/j.ecoenv.2017.04.009 . Crossref, Medline,

Lee H-J, Jedrychowski MP, Vinayagam A, Wu N, Shyh-Chang N, Hu Y, et al. 2017. Proteomic and metabolomic characterization of a mammalian cellular transition from quiescence to proliferation. Cell Rep 20(3):721–736, PMID:28723573,10.1016/j.celrep.2017.06.074 . Crossref, Medline,

Levenson AS, Jordan VC . 1997. MCF-7: the first hormone-responsive breast cancer cell line. Cancer Res 57(15):3071–3078, PMID:9242427 . Medline,

Li B, Dewey CN . 2011. RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics 12:323, PMID:21816040,10.1186/1471-2105-12-323 . Crossref, Medline,

Li P, Jin J, Wang Y, Hu J, Xu M, Sun Y, et al. 2017. Concentrations of organophosphorus, polybromobenzene, and polybrominated diphenyl ether flame retardants in human serum, and relationships between concentrations and donor ages. Chemosphere 171:654–660. PMID: 28056452,10.1016/j.chemosphere.2016.12.126 . Crossref, Medline,

Li X, Sun H, Yao Y, Zhao Z, Qin X, Duan Y, et al. 2018a. Distribution of phthalate metabolites between paired maternal–fetal samples. Environ Sci Technol 52(11):6626–6635, PMID:29754483,10.1021/acs.est.8b00838 . Crossref, Medline,

Li L, Westgate JN, Hughes L, Zhang X, Givehchi B, Toose L, et al. 2018b. A model for risk-based screening and prioritization of human exposure to chemicals from near-field sources. Environ Sci Technol 52:14235–14244, PMID:30407800,10.1021/acs.est.8b04059 . Crossref, Medline,

Lioy PJ, Freeman NC, Millette JR . 2002. Dust: a metric for use in residential and building exposure assessment and source characterization. Environ Health Perspect 110(10):969–983, PMID:12361921,10.1289/ehp.02110969 . Link,

Liu M, Jia S, Dong T, Han Y, Xue J, Wanjaya ER, et al. 2019. The occurrence of bisphenol plasticizers in paired dust and urine samples and its association with oxidative stress. Chemosphere 216:472–478, PMID:30388685,10.1016/j.chemosphere.2018.10.090 . Crossref, Medline,

Marlow R, Strickland P, Lee JS, Wu X, Pebenito M, Binnewies M, et al. 2008. SLITs suppress tumor growth in vivo by silencing Sdf1/Cxcr4 within breast epithelium. Cancer Res 68(19):7819–7827, PMID:18829537,10.1158/0008-5472.CAN-08-1357 . Crossref, Medline,

McGee SP, Cooper EM, Stapleton HM, Volz DC . 2012. Early zebrafish embryogenesis is susceptible to developmental TDCPP exposure. Environ Health Perspect 120(11):1585–1591, PMID:23017583,10.1289/ehp.1205316 . Link,

Meeker JD, Cooper EM, Stapleton HM, Hauser R . 2013. Urinary metabolites of organophosphate flame retardants: temporal variability and correlations with house dust concentrations. Environ Health Perspect 121(5):580–585,10.1289/ehp.1205907 . Link,

Moffatt BA, Ashihara H . 2002. Purine and pyrimidine nucleotide synthesis and metabolism. Arabidopsis Book 1:e0018. PMID: 22303196,10.1199/tab.0018 . Crossref, Medline,

Nagayama J, Tsuji H, Takasuga T . 2000. Comparison between brominated flame retardants and dioxins or organochlorine compounds in blood levels of Japanese adults. Organohalogen Compd 48:27–30.

Nahle Z, Polakoff J, Davuluri RV, McCurrach ME, Jacobson MD, Narita M, et al. 2002. Direct coupling of the cell cycle and cell death machinery by E2F. Nat Cell Biol 4(11):859–864, PMID:12389032,10.1038/ncb868 . Crossref, Medline,

Nguyen T, Nioi P, Pickett CB . 2009. The Nrf2-antioxidant response element signaling pathway and its activation by oxidative stress. J Biol Chem 284(20):13291–13295, PMID:19182219,10.1074/jbc.R900010200 . Crossref, Medline,

Oakley RH, Cidlowski JA . 2013. The biology of the glucocorticoid receptor: new signaling mechanisms in health and disease. J Allergy Clin Immunol 132(5):1033–1044, PMID:24084075,10.1016/j.jaci.2013.09.007 . Crossref, Medline,

Li P, Li Q-X, Ma Y-L, Jin J, Wang Y, Tian Y . 2015. Determination of organophosphate esters in human serum using gel permeation chromatograph and solid phase extraction coupled with gas chromatography-mass spectrometry. Chinese J Anal Chem 43(7):1033–1039,10.1016/S1872-2040(15)60840-4 . Crossref,

Petersen H, Biereichel A, Burseg K, Simat TJ, Steinhart H . 2008. Bisphenol a diglycidyl ether (BADGE) migrating from packaging material ‘disappears’ in food: reaction with food components. Food Addit Contam Part A Chem Anal Control Expo Risk Assess 25(7):911–920, PMID:18569010,10.1080/02652030701837399 . Crossref, Medline,

Quéméneur L, Gerland L-M, Flacher M, Ffrench M, Revillard J-P, Genestier L . 2003. Differential control of cell cycle, proliferation, and survival of primary T lymphocytes by purine and pyrimidine nucleotides. J Immunol 170(10):4986–4995, PMID:12734342,10.4049/jimmunol.170.10.4986 . Crossref, Medline,

Rappaport SM, Barupal DK, Wishart D, Vineis P, Scalbert A . 2014. The blood exposome and its role in discovering causes of disease. Environ Health Perspect 122(8):769–774, PMID:24659601,10.1289/ehp.1308015 . Link,

Rappaport SM, Smith MT . 2010. Environment and disease risks. Science 330(6003):460–461, PMID:20966241,10.1126/science.1192603 . Crossref, Medline,

Ren B, Cam H, Takahashi Y, Volkert T, Terragni J, Young RA . 2002. E2F integrates cell cycle progression with DNA repair, replication, and G2/M checkpoints. Genes Dev 16:245–256, PMID:11799067,10.1101/gad.949802 . Crossref, Medline,

Repouskou A, Panagiotidou E, Panagopoulou L, Bisting PL, Tuck AR, Sjödin MO, et al. 2019. Gestational exposure to an epidemiologically defined mixture of phthalates leads to gonadal dysfunction in mouse offspring of both sexes. Sci Rep 9(1):6424, PMID:31015488,10.1038/s41598-019-42377-6 . Crossref, Medline,

Robinson MD, McCarthy DJ, Smyth GK . 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139–140, PMID:19910308,10.1093/bioinformatics/btp616 . Crossref, Medline,

Samet JM, Wages PA . 2018. Oxidative stress from environmental exposures. Curr Opin Toxicol 7:60–66, PMID:30079382,10.1016/j.cotox.2017.10.008 . Crossref, Medline,

Silva E, Rajapakse N, Kortenkamp A . 2002. Something from “nothing”—eight weak estrogenic chemicals combined at concentrations below NOECs produce significant mixture effects. Environ Sci Technol 36(8):1751–1756, PMID:11993873,10.1021/es0101227 . Crossref, Medline,

Strickland P, Shin GC, Plump A, Tessier-Lavigne M, Hinck L . 2006. SLIT2 and netrin 1 act synergistically as adhesive cues to generate tubular bi-layers during ductal morphogenesis. Development 133(5):823–832, PMID:16439476,10.1242/dev.02261 . Crossref, Medline,

Suzuki T, Kitamura S, Khota R, Sugihara K, Fujimoto N, Ohta S . 2005. Estrogenic and antiandrogenic activities of 17 benzophenone derivatives used as UV stabilizers and sunscreens. Toxicol Appl Pharmacol 203(1):9–17, PMID:15694459,10.1016/j.taap.2004.07.005 . Crossref, Medline,

Suzuki T, Motohashi H, Yamamoto M . 2013. Toward clinical application of the Keap1–Nrf2 pathway. Trends Pharmacol Sci 34(6):340–346, PMID:23664668,10.1016/j.tips.2013.04.005 . Crossref, Medline,

Tautenhahn R, Patti GJ, Rinehart D, Siuzdak G . 2012. XCMS online: a web-based platform to process untargeted metabolomic data. Anal Chem 84(11):5035–5039, PMID:22533540,10.1021/ac300698c . Crossref, Medline,

Thomsen C, Lundanes E, Becher G . 2002. Brominated flame retardants in archived serum samples from Norway: a study on temporal trends and the role of age. Environ Sci Technol 36(7):1414–1418, PMID:11999045,10.1021/es0102282 . Crossref, Medline,

Trapnell C, Pachter L, Salzberg SL . 2009. TopHat: discovering splice junctions with RNA-seq. Bioinformatics 25(9):1105–1111, PMID:19289445,10.1093/bioinformatics/btp120 . Crossref, Medline,

Tye BK . 1999. MCM proteins in DNA replication. Annu Rev Biochem 68:649–686, PMID:10872463,10.1146/annurev.biochem.68.1.649 . Crossref, Medline,

Uehara M, Arai Y, Watanabe S, Adlercreutz H . 2000. Comparison of plasma and urinary phytoestrogens in Japanese and Finnish women by time‐resolved fluoroimmunoassay. Biofactors 12(1–4):217–225,10.1002/biof.5520120134 . Crossref, Medline,

Wadhwa RR, Marco C . Steady state concentration. StatPearls [Internet]. StatPearls Publishing, 2020, PMID:31985925 , https://www.ncbi.nlm.nih.gov/books/NBK553132/. Medline,

Wang S-Y, Kuo C-H, Tseng YJ . 2013. Batch normalizer: a fast total abundance regression calibration method to simultaneously adjust batch and injection order effects in liquid chromatography/time-of-flight mass spectrometry-based metabolomics data and comparison with current calibration methods. Anal Chem 85(2):1037–1046, PMID:23240878,10.1021/ac302877x . Crossref, Medline,

Wang L, Xue J, Kannan K . 2015a. Widespread occurrence and accumulation of bisphenol a diglycidyl ether (BADGE), bisphenol F diglycidyl ether (BFDGE) and their derivatives in human blood and adipose fat. Environ Sci Technol 49(5):3150–3157, PMID:25679329,10.1021/acs.est.5b00096 . Crossref, Medline,

Wang Q, Lam JC-W, Man Y-C, Lai NL-S, Kwok KY, Yong Guo Y, et al. 2015b. Bioconcentration, metabolism and neurotoxicity of the organophorous flame retardant 1, 3-dichloro 2-propyl phosphate (TDCPP) to zebrafish. Aquat Toxicol 158:108–115, PMID:25461749,10.1016/j.aquatox.2014.11.001 . Crossref, Medline,

Wang W, Abualnaja KO, Asimakopoulos AG, Covaci A, Gevao B, Johnson-Restrepo B, et al. 2015c. A comparative assessment of human exposure to tetrabromobisphenol A and eight bisphenols including bisphenol A via indoor dust ingestion in twelve countries. Environ Int 83:183–191. PMID: 26177148,10.1016/j.envint.2015.06.015 . Crossref, Medline,

Wang Y, Wang W, Xu L, Zhou X, Shokrollahi E, Felczak K, et al. 2016. Cross talk between nucleotide synthesis pathways with cellular immunity in constraining hepatitis E virus replication. Antimicrob Agents Chemother 60(5):2834–2848, PMID:26926637,10.1128/AAC.02700-15 . Crossref, Medline,

Wang Z, Cui B, Zhang F, Yang Y, Shen X, Li Z, et al. 2019. Development of a correlative strategy to discover colorectal tumor tissue derived metabolite biomarkers in plasma using untargeted metabolomics. Anal Chem 91(3):2401–2408, PMID:30580524,10.1021/acs.analchem.8b05177 . Crossref, Medline,

Weckwerth W . 2003. Metabolomics in systems biology. Annu Rev Plant Biol 54:669–689, PMID:14503007,10.1146/annurev.arplant.54.031902.135014 . Crossref, Medline,

Wishart DS, Feunang YD, Marcu A, Guo AC, Liang K, Vázquez-Fresno R, et al. 2018. HMDB 4.0: the human metabolome database for 2018. Nucleic Acids Res 46(D1):D608–D617, PMID:29140435,10.1093/nar/gkx1089 . Crossref, Medline,

Wishart DS, Tzur D, Knox C, Eisner R, Guo AC, Young N, et al. 2007. HMDB: the human metabolome database. Nucleic Acids Res 35(Database issue):D521–D526, PMID:17202168,10.1093/nar/gkl923 . Crossref, Medline,

Wu J, Mao X, Cai T, Luo J, Wei L . 2006. KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res 34(Web Server issue):W720–W724, PMID:16845106,10.1093/nar/gkl167 . Crossref, Medline,

Xi H, Kersh GJ . 2004. Early growth response gene 3 regulates thymocyte proliferation during the transition from CD4− CD8− to CD4 + CD8 + 1. J Immunol 172(2):964–971, PMID:14707069,10.4049/jimmunol.172.2.964 . Crossref, Medline,

Xia P, Zhang X, Zhang H, Wang P, Tian M, Yu H . 2017. Benchmarking water quality from wastewater to drinking waters using reduced transcriptome of human cells. Environ Sci Technol 51(16):9318–9326, PMID:28696678,10.1021/acs.est.7b02648 . Crossref, Medline,

Xiao Z, Feng J, Shi Z, Li J, Zhao Y, Wu Y . 2011. Determination of three brominated flame retardants in human serum using solid-phase extraction coupled with ultra-performance liquid chromatography-tandem mass spectrometry and gas chromatography-mass spectrometry. Chinese J Chromatog 29:1165–1172, PMID:22500441,10.3724/SP.J.1123.2011.01165 . Crossref, Medline,

Yang Y, Guan J, Yin J, Shao B, Li H . 2014. Urinary levels of bisphenol analogues in residents living near a manufacturing plant in south China. Chemosphere 112:481–486,10.1016/j.chemosphere.2014.05.004 . Crossref, Medline,

Yue Y, Liu J, Liu R, Sun Y, Li X, Fan J . 2014. The binding affinity of phthalate plasticizers-protein revealed by spectroscopic techniques and molecular modeling. Food Chem Toxicol 71:244–253, PMID:24994562,10.1016/j.fct.2014.06.022 . Crossref, Medline,

Zhang Y, Keerthisinghe TP, Han Y, Liu M, Wanjaya ER, Fang M . 2018. “Cocktail” of xenobiotics at human relevant levels reshapes the gut bacterial metabolome in a species-specific manner. Environ Sci Technol 52(19):11402–11410, PMID:30153011,10.1021/acs.est.8b02629 . Crossref, Medline,

Zhao F, Chen M, Gao F, Shen H, Hu J . 2017. Organophosphorus flame retardants in pregnant women and their transfer to chorionic villi. Environ Sci Technol 51(11):6489–6497, PMID:28516762,10.1021/acs.est.7b01122 . Crossref, Medline,

Zhou J, Cidlowski JA . 2005. The human glucocorticoid receptor: one gene, multiple proteins and diverse responses. Steroids 70(5–7):407–417, PMID:15862824,10.1016/j.steroids.2005.02.006 . Crossref, Medline,