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

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. https://doi.org/10.1289/EHP6641


Introduction
Humans are ubiquitously exposed to many xenobiotics on a dayto-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 . 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 . 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.

Experimental Workflow
Mixtures consisting of 23 xenobiotic chemicals at several humanrelevant 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. 2015Li et al. , 2017Li et al. , 2018aNagayama 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-oneout method; see below) was developed to screen the causal chemicals, by characterizing the metabolome changes of the mixtures with and without each tested compound.

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,3dichloroisopropyl)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)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(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) The 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. b The percentage is calculated based on the molar concentration of each compound.
c Use metabolite as a replacement for the parent compound: 2,4-dihydroxybenzophenone (DHB) replaced its parent compound benzophenone-3 (BP-3). The 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.
ð1Þ + ð2Þ where t 1=2 = Half-life was predicted based on the QSAR approach called IFS (Arnot et al. 2014;Brown et al. 2012); V d = volume of distribution; CI = clearance; F = bioavailability; I D = dosing interval; DI = daily intake (Equation 4) (Guo et al. 2014); where Cu: Urinary concentration, V u = 2L (urine volume), BW: body weight, F ue : 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 Figure 1. Workflow for metabolome and transcriptome analysis and assessment on the relative contribution of each component to the observed effect using a "leave-one-out" method. Metabolome and transcriptome analysis: Omics technologies metabolomics and transcriptomics were employed to characterize the biological effect on cells of exposure to the 23 chemical mixture at several human-relevant levels [GM (geometric mean), 2GM (2-fold geometric mean), 5GM, Max (maximum), and 5Max (5-fold maximum)]. The dysregulated genes and metabolites were also integrated to interpret the biological perturbation triggered by the mixture exposure. Leave-one-out method: A high-throughput targeted metabolomics method was developed to evaluate the relative contribution of each component to the observed effect, by characterizing the metabolome changes of the mixtures with and without each tested compound (Total 25 groups: including 23 treatment groups each with one component removed, one 23 chemical xenobiotic mixture labeled as "All Compounds" and one control group).
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. 2015Li et al. , 2017Li et al. , 2018aNagayama 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 flatbottomed plates (Thermo Fisher Scientific) with a seeding density of 30,000 cells/well in 100 lL DMEM for 7 h. Subsequently, the 100 lL 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 lL resazurin (0:50 mg=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 × 10 6 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-lL samples to a 170-lL TE buffer (10 mM Tris-Cl pH 8.0, 1 mM EDTA pH 8.0). Next, 50 lL of diluted samples and PicoGreen working reagent (diluted 200 times) were added to 96well 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 × 10 6 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:6 mL 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,000 r: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-lL diluted sample (100 times) was mixed with 100 lL of working reagents in 96well 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 lL, sonicated for 10 min, and centrifuged for 15 min at 13,000 r:p:m: and 4°C to remove insoluble debris. The supernatants were transferred to highperformance 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 NH 2 Column (particle size 3 lm, i.d. 2 × length 150 mm, for polar metabolites) and Waters Atlantis T3 column (particle size 3 lm, i.d. 2:1 × length 100 mm, for hydrophobic metabolites) were used for separation in the negative and positive modes, respectively. For NH 2 column, the injection volume was 6 lL, 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:25 mL=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 lL. 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:5 mL=min and column temperature was 25°C. Other parameters were set as follows: gas temperature: 225°C; gas flow: 14 L=min; sheath gas temperature: 275°C; sheath gas flow: 11 L=min; nozzle voltage: 1,000 V; nebulizer pressure: 35 psig; capillary voltage: 3,500 V and −3,500 V for positive (ESI+, electrospray ionization) and negative mode (ESI −); nozzle voltage: 1,000 V; fragmentor: 175 V. 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 20 V and 40 V. 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 (Dm=z = 15 ppm; minimum peak width = 10 s; and minimum peak width = 120 s); 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 jfold changej > 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 lg) was used to construct sequencing library. Subsequently, RNA-seq transcriptome library was prepared following TruSeqTM RNA sample preparation Kit from Illumina using 1 lg 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, doublestranded 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 j log 2 ðFold changeÞj ≥ 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. Twentythree 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 abovementioned 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 lm, i.d. 2:1 × length 100 mm) in ESI − with an injection volume of 10 lL, 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:5 mL=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.

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  Table S5). (B) Normalized DNA concentrations (percent of control samples) of cells after exposure to the mixture with varying levels using the PicoGreen assay (GM, 2GM, 5GM, Max, and 5Max; n = 3; mean ± SD; complete DNA data are available in Table S6); * and ** represent p ≤ 0:05 and p ≤ 0:01 when treatment group is statistically different from nontreated controls using one-way ANOVA and Duncan's post hoc analysis. (C) Global metabolomic profiling of cells after exposure to the mixture with varying levels (GM, 2GM, 5GM, Max, and 5Max; n = 4). Up-and down-regulated significant features detected by global profiling (significance criteria: fold change >1:3, p ≤ 0:05, and maximal intensity >10,000 ion counts; the percentage was calculated based on total detected feature numbers; complete profiling data are available in Table S7). (D-E) Classic Venn diagram for the number of shared and distinct features in ESI-negative mode and ESI-positive mode, respectively. 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.

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 jfold changej > 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.
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 downregulated, 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 upregulated. 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 0 -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.

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).

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, log 2 ð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,

Supervised Pathway Analysis
A supervised approach was also adopted by mapping DEGs onto some cellular toxicity pathways (Xia et al. 2017) ( Figure 5C).

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  Table S10. (B) Unsupervised pathway analysis for the DEGs using the KOBAS 3.0 with an adjusted p ≤ 0:05 (corrected by Benjamini-Hochberg method); enrichments with top three rankings for each database (PANTHER, KEGG, and REACTOME) were listed for illustration. Complete pathway data are available in Table S11. (C) Supervised in vitro bioassay pathway enrichment analysis; p-values were determined by Fisher's exact test and Bonferroni correction; ** represents the adjusted p ≤ 0:01. Complete pathway data are available in Table S12. Note: ACHE, acetylcholinesterase; AHR, aryl hydrogen receptor; AP-1, activator protein-1; AR, androgen receptor; CAR, constitutive androstane receptor; DEGs, differentially expressed genes; ER, estrogen receptor; FXR, farnesoid X receptor; GR, glucocorticoid receptor; HIF, hypoxia-inducible factor; HSE, heat-shock element; LXR, liver X receptor; Nrf2, nuclear factor erythroid 2-related factor; PPAR, peroxisome proliferator-activated receptor gamma; PR, progesterone receptor; PXR, pregnane X receptor; RAR, retinoic acid receptor; ROR, retinoid acid-related orphan receptor; RARE, retinoic acid response element; THRa, thyroid hormone receptor alpha; TP53, tumor protein 53.
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.

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-2H 2 O) and bisphenol A diglycidyl ether (BADGE), showed great similarity with the "All Compounds" treatment, indicating their minimal contribution to the overall metabolome changes. 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 upregulation 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 (log 2 ðFCÞ = 1:09 at Max) and lactate production (0.79) (see Figure S7), are consistent with an increase Figure 7. Assessment on the relative contribution of each component to the observed effect using the high-throughput metabolomics-based leave-one-out method (n = 3). Total 25 groups: including 23 treatment groups with one component removed and with the others remaining at Max level, one 23 chemical mixture at Max level labeled as "All Compounds" and one DMSO control. (A) Intuitive hierarchical clustering analysis for grouping similar objects among the targeted 19 metabolites across 25 groups; each colored cell on the map corresponds to a concentration value; clustering of metabolites was based on Ward's clustering method with the Euclidean distance. Bar graph of the fold change (mean ± SEM) of (B), UDP; (C) Lactate; (D) ATP after the removal of each component; y-axis represents the fold change by comparing treated samples with control samples; each color/pattern represents a chemical category; symbol a indicates one treatment group is statistically different from the DMSO controls by one-way ANOVA and Duncan's post hoc analysis, "a" denotes p ≤ 0:05 and "aa" denotes p ≤ 0:01; symbol b indicates one treatment group is statistically different from the "All Compounds" group; "b" denotes p ≤ 0:05, and "bb" denotes p ≤ 0:01. Complete data are available in Table S15. Note: ADP, adenosine diphosphate; ATP, adenosine triphosphate; BADGE, bisphenol A diglycidyl ether; BADGE-2H 2 O, bisphenol A bis(2,3-dihydroxypropyl) ether; BFDGE, bisphenol F-diglycidyl ether; b-G6P, beta-D-glucose-6-phosphate;BPA, bisphenol A; TDCPP, tris (1,3-dichloroisopropyl) phosphate; BPS, bisphenol S; TPhP, triphenyl phosphate; BzBP, benzylbutyl phthalate; CDP, cytidine 5 0 -diphosphate; CTP, cytidine triphosphate; DnBP, di-n-butyl phthalate; DHB: 2,4-dihydroxybenzophenone; DiDP, diisodecyl phthalate; DiNP, diisononyl phthalate; GDP, guanosine diphosphate; GDP-mannose, guanosine diphosphate mannose; GTP, guanosine triphosphate; MBzP, monobenzyl phthalate; MEHP, mono-ethylhexyl phthalate; MnBP, mono-n-butyl phthalate; PFOA, perfluorooctanoic acid; Phosphocolamine, O-Phosphorylethanlamine; TBBPA, tetrabromobisphenol A; UDP, uridine diphosphate; UDP-GlcA, uridine diphosphate glucuronic acid; UDP-GlcNAc, uridine diphosphate N-acetylglucosamine; UDP-glucose, uridine diphosphate glucose; UTP, uridine 5 0 triphosphate. 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 leaveone-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.