High-Throughput Transcriptomics of Nontumorigenic Breast Cells Exposed to Environmentally Relevant Chemicals

Abstract

Background: There is a suite of chemicals, including metals, pesticides, and personal care product compounds, which are commonly detected at high levels in US Center for Disease Control's National Health and Nutrition Examination Survey (NHANES) chemical biomarker screens. Whether these chemicals influence development of breast cancer is not well understood. Objectives: The objectives were to perform an unbiased concentration-dependent assessment of these chemicals, to quantify differences in cancer-specific genes and pathways, to describe if these differences occur at human population-relevant concentrations, and to specifically test for differences in markers of stemness and cellular plasticity. Methods: We treated nontumorigenic mammary epithelial cells, MCF10A, with 21 chemicals at four concentrations (25 nM, 250 nM, 2.5 lM, and 25 lM) for 48 h. We conducted RNA-sequencing for these 408 samples, adapting the plexWell plate-based RNA-sequencing method to analyze differences in gene expression. We calculated gene and biological pathway-specific benchmark concentrations (BMCs) using BMDExpress3, identifying differentially expressed genes and generating the best fit benchmark concentration models for each chemical across all genes. We identified enriched biological processes and pathways for each chemical and tested whether chemical exposures change predicted cell type distributions. We contextual-ized benchmark concentrations relative to human population biomarker concentrations in NHANES. Results: We detected chemical concentration-dependent differences in gene expression for thousands of genes. Enrichment and cell type distri-bution analyses showed benchmark concentration responses correlated with differences in breast cancer-related pathways, including induction of basal-like characteristics for some chemicals, including arsenic, lead, copper, and methyl paraben. Comparison of benchmark data to NHANES chemical biomarker (urine or blood) concentrations indicated an overlap between exposure levels and levels sufficient to cause a gene expression response. Discussion: These analyses revealed that many of these 21 chemicals resulted in differences in genes and pathways involved in breast cancer in vitro at human exposure-relevant concentrations.

Headnote

Background: There is a suite of chemicals, including metals, pesticides, and personal care product compounds, which are commonly detected at high levels in US Center for Disease Control's National Health and Nutrition Examination Survey (NHANES) chemical biomarker screens. Whether these chemicals influence development of breast cancer is not well understood.

Objectives: The objectives were to perform an unbiased concentration-dependent assessment of these chemicals, to quantify differences in cancer-specific genes and pathways, to describe if these differences occur at human population-relevant concentrations, and to specifically test for differences in markers of stemness and cellular plasticity.

Methods: We treated nontumorigenic mammary epithelial cells, MCF10A, with 21 chemicals at four concentrations (25 nM, 250 nM, 2.5 lM, and 25 lM) for 48 h. We conducted RNA-sequencing for these 408 samples, adapting the plexWell plate-based RNA-sequencing method to analyze differences in gene expression. We calculated gene and biological pathway-specific benchmark concentrations (BMCs) using BMDExpress3, identifying differentially expressed genes and generating the best fit benchmark concentration models for each chemical across all genes. We identified enriched biological processes and pathways for each chemical and tested whether chemical exposures change predicted cell type distributions. We contextual-ized benchmark concentrations relative to human population biomarker concentrations in NHANES.

Results: We detected chemical concentration-dependent differences in gene expression for thousands of genes. Enrichment and cell type distri-bution analyses showed benchmark concentration responses correlated with differences in breast cancer-related pathways, including induction of basal-like characteristics for some chemicals, including arsenic, lead, copper, and methyl paraben. Comparison of benchmark data to NHANES chemical biomarker (urine or blood) concentrations indicated an overlap between exposure levels and levels sufficient to cause a gene expression response.

Discussion: These analyses revealed that many of these 21 chemicals resulted in differences in genes and pathways involved in breast cancer in vitro at human exposure-relevant concentrations. https://doi.org/10.1289/EHP12886

Introduction

There are an estimated 2.3 million new cases of breast cancer diagnosed globally and >685,000 people die from breast can-cer each year.1 Breast cancer is a heterogeneous disease- tumors can be grouped into subtypes based on their molecular signature and expression of estrogen receptor (ER), progester-one receptor (PR), and human epidermal growth factor receptor 2 (HER2).2 Prognosis and outcomes vary based on subtype. Triple-negative breast cancers (TNBCs), those which do not express ER, PR, or HER2, have the worst overall survival rates.3 Identifying the factors that promote aggressive breast can-cers, such as triple-negative breast cancers, is of high public health priority. In general, high penetrance genetic risk factors are only thought to explain <5-10% of all breast cancers.4 As such, there is likely a large impact of the environment on aggres-sive breast cancer risk. There is an urgent need to identify poten-tially modifiable risk factors, such as chemical exposures, that promote aggressive breast cancers so that we can develop new strategies for prevention.

As aggressive breast cancers form, they acquire the hallmarks of cancer, including sustained proliferative signaling, the activation of invasion pathways, immortality, the evasion of growth suppression, the reactivation of stem cell pathways, and the acquisition of cellular plasticity.5-7 As an example, triple-negative breast cancers are linked to an inflammatory microenvironment and enriched for expression of stem cell-associated pathways.8-12 Triple-negative breast cancers also often present with a basal-like phenotype.13 Experimental data surprisingly suggest that these basal-like breast cancers develop from dysregulated luminal cells, where these cells acquire cell state phenotypic plasticity and shift from a luminal-like to a basal-like expression pattern.14,15 Identifying chemicals that impact the development of aggressive breast cancers through the promotion of these various hallmarks, including altered stemness and cell state phenotypic plasticity pathways, would provide novel insights into the mechanisms by which environmental exposures impact breast cancer risk and outcomes.

There are hundreds or thousands of chemical exposures that could impact breast cancer risk, and a systematic and rapid method for evaluation is essential to identify those which may alter aggressive breast cancer hallmarks at human-relevant con-centrations.4 Here, we developed a high-throughput transcrip-tomic platform for the unbiased assessment of chemical effects in vitro. We quantified the concentration-dependent effects of 21 chemical exposures in the commonly used in vitro breast carcinogenesis model, nontumorigenic MCF10A cells. We pri-oritized the chemicals based on them being commonly detected in chemical biomarker screens, with documented exposure dispar-ities across demographic groups, and with putative links to breast cancer.4,16-22 We had four goals, as follows: first, to assess and report concentration-dependent gene expression differences in an unbiased fashion; second, to test the hypothesis that these gene expression differences will be enriched in established breast can-cer hallmark pathways; third, to contextualize the concentrations at which gene expression differences occur relative to biomarker concentrations in a representative sample of the US population; and, fourth, to test the hypothesis that treatment with these chemicals will lead to differences in cell state phenotypes consistent with a luminal-to-basal transition. Overall, we evaluated putative chemical contributors at human relevant concentrations on bio-logical pathways associated with breast cancers.

Methods

Cell Culture

The nontumorigenic mammary epithelial cell line, MCF10A, was obtained from American Type Culture Collection (ATCC). Cells were cultured and maintained in Dulbecco's modified Eagle's medium/Hams F12 50/50 mix (Corning) supplemented with 5% horse serum (Thermo Fisher), 2:5 mg=mL HEPES (Thermo Fisher), 5 lg=mL insulin (Gibco/Thermo Fisher), 96 lg=mL hydrocortisone (StemCell), 100 lg=mL cholera toxin (Sigma-Aldrich), and 20ng=mL epidermal growth factor (StemCell). The cells were maintained at 37oC in a humidified 5% CO2 condition.

Compound Preparation

The 21 chemicals [sodium (meta) arsenite, cadmium chloride, cop-per (II) chloride, lead (II) acetate trihydrate, mercury (II) chloride, phenanthrene, p,p' -DDE (dichlorodiphenyldichloroethylene), thiram, 1,4-dichlorobenzene, 2,5-dichlorophenol, methylparaben, propyl-paraben, perfluorodecanoic acid (PFDA), perfluorononanoic acid (PFNA), perfluorooctanoic acid (PFOA), perfluorooctanesulfonic acid (PFOS), polychlorinated biphenyl 153 (PCB153), poly-chlorinated biphenyl 187 (PCB 187), bisphenol A, bisphenol F, bisphenol S] were purchased from Sigma-Aldrich and Cayman Chemicals (Table S1). The chemicals were weighed and dissolved in dimethyl sulfoxide (DMSO) (for most) and in sterile water (for some) at a concentration of 5 mg=mL and stored at -20°C for long-term storage. For dosing, we prepared an intermediate stock concentration of 5 mM for each of the chemicals, which was fur-ther diluted into final concentrations for dosing the MCF10A cells: 25 lM, 2:5 lM, 0:25 lM, and 0:025 lM. These final concentrations were prepared fresh in MCF10A media for each experiment. For chemicals dissolved in DMSO, the final DMSO concentration was 0.5%, as were the vehicle controls. For chemicals dissolved in water, an equivalent amount of water was included in the vehicle control.

Cell Plating

MCF10A cells were cultured and expanded at passage 105 for all rounds of cell culture, and different flasks were used for different plates. Cells were plated in black, clear flat bottom 384-well plates (Corning) using a 2:5- to 125-lL multichannel pip-ette. First, 10 lL of plain medium was added to prewarmed plates. Five hundred cells in 30 lL of medium were added to each well at a slow speed to make up a final volume of 40 lL. The cells were left undisturbed to allow the cells to attach and dis-tribute uniformly over the well. Each microplate was incubated at 37°C in a humidified incubator with 5% CO2 for 24 h prior to dosing. Cells were then exposed to test chemicals for 48 h with four bi-ological replicates per concentration, including four replicates for vehicle controls (DMSO or water). Experiments were run in two batches to include experimental replicates for all chemical-concentration combinations.

plexWell RNA-Sequencing Overview

To conduct high-throughput RNA-sequencing, we developed an adapted version of the plexWell (SeqWell) plate-based sequencing method to make the method applicable for high-throughput in vitro toxicological analysis. In the plexWell protocol, sequencing reac-tions were prepared in custom 96-well plates, where each well con-tained a unique oligo barcode. Plates were also assigned a unique oligo barcode. Cells were lysed and cDNA was directly prepared from the cellular lysate. A tagmentation reaction added the well-specific barcode to the cDNA from each well. The tagmented cDNAs were then pooled prior to library preparation, where plate-specific barcodes can also be incorporated for ultrahigh sample multiplexing prior to sequencing. For example, 384 wells of treated samples would be processed across four custom 96-well plates, where each of the samples received a well-specific barcode and each of the four plates received a plate-specific barcode. This led to each of the 384 samples receiving a unique combination of well and plate barcodes, which can be bioinformatically deconvoluted following sequencing. This protocol decreased the cost of library preparation by <90% relative to traditional RNA-sequencing, ena-bling unbiased high-throughput toxicological analyses.

Our key modification to the manufacturer protocol to stream-line the process was to avoid the step of RNA purification prior to cDNA preparation. The plexWell protocol calls for the input of 1 lL of RNA for cDNA preparation. To eliminate the step of RNA extraction, we devised a method of lysing the cells directly in the 384-well microplate. After 48 h of incubation with the test chemicals, the cells were washed twice with 40 lL of 1× Hanks' balanced salt solution (HBSS). After the second wash, 5 lL of 1× HBSS was left behind to avoid cell loss. We added 10 lL of plexWell lysis buffer (with RNase inhibitor) to each well using a 2:5- to 125-lL multichannel micropipette. The buffer was pipetted up and down to mix 10 times, and the cell lysates were collected in 1 mL Eppendorf tubes on ice to avoid RNA degradation. Cell lysis was confirmed via examination by microscopy.

plexWell cDNA Preparation

The plexWell manufacturer's protocol was followed to prepare cDNA for subsequent library preparation. In brief, 1 lL of cell lysate was used as an input for cDNA preparation and oligo dT annealing in a PCR 96-well plate (Dot Scientific). The cDNA synthesis reaction was run for 12 cycles on a Bio-Rad CFX96 Touch real-time PCR detection system. Equivalent amount of MAGwise paramagnetic beads was added to the amplified cDNA in each well of the 96-well plate. The cDNA was allowed to bind for 5 min. The plate was placed on a magnet to allow the beads to pellet on the inner walls of each well. Supernatant was discarded, and the beads were washed once using 80% ethanol. The cDNA was eluted using 20 lL of 10 mM Tris solution. Purified cDNA was stored at -20°C for short-term storage.

PicoGreen Assay

As per the plexWell protocol, we quantified the cDNA concentration of a subset of samples using the Quant-iT PicoGreen double stranded DNA (dsDNA) assay kit (Thermo Fisher). Standards ranging from 25 pg=mL to 25 ng=mL were prepared using the Lambda Standard DNA (100 lg=mL) provided in the kit. One microliter of prepared dsDNA was diluted in 99 lL of 1× TE buffer to get a 1:100 dilution. A 200-fold aqueous Quant-iT PicoGreen working solution was prepared in 1× TE buffer at a volume of 100 lL=well. The standards and samples were plated in a flat bottom Corning 96-well plate using a 125- to 1,250-lL multichannel pipette. Quant-iT PicoGreen working solution added at 100 lL=well using a 125- to 1,250-lL multichannel pipette. The samples were incubated at room temperature for 5 min and protected from light. The fluorescence from Quant-iT picoGreen working solution was read on SpectraMax M5e micro-plate reader (Molecular Devices). The preset protocol on SoftMaxPro software version 5.4 for PicoGreen assay for nucleic acid was used for analysis.

As per the plexWell protocol, six samples per 96-well plate were analyzed on an Agilent high-sensitivity DNA bioanalyzer (up to 15,000 base pairs) at the University of Michigan Advanced Genomics Core. The analyzed data showed electropherograms with a summary of fragment sizes. The electropherograms for submitted samples were compared to the example electropherograms provided in the plexWell protocol to check for fragment size discrepancies. A pool of 22 random samples was submitted for QuBit analysis at the Advanced Genomics Core per round of experiment to obtain an average concentration. According to the plexWell protocol, the reagents in the kit were formulated to tol-erate up to 10-fold difference in sample input across 96 samples. After fragment distribution and concentration checks, we calcu-lated a global dilution factor using the QuBit concentrations (one per library preparation) as described in the plexWell protocol. This global dilution factor gives the volume for dilution of cDNA prior to tagmentation. Four microliters of prepared cDNA was diluted in the calculated Tris-HCl volume.

Library Preparation

Post dilution, 4 lL of the diluted cDNA per sample was used for library preparation. Each sample was individually barcoded in a hard skirted Sample Barcode (SB) plate provided in the plexWell 384 Rapid Single Cell RNA library prep kit (SeqWell). Each sample was labeled with an i7 index, referred to as a Sample Barcode, using a tagmentation reaction. Post i7 tagging, the 96 samples were pooled to a final volume of 800-860 lL. An equivalent amount of MAGwise paramagnetic beads was added to the pooled SB reac-tions. The cDNA was allowed to bind to the beads for 5 min. The beads then formed a pellet on the inner wall of the tube as it was placed on a magnetic stand. The pellet was washed two times with 80% ethanol. Forty microliters of 10 mM Tris was used to elute out purified SB reaction pool. Next, the purified SB reaction pool was labeled with an i5 index, referred to as a Pool Barcode (PB), using a tagmentation reaction. Post i5 tagging, an equivalent amount of MAGwise paramagnetic beads was added to the PB reaction. The DNA was allowed to bind to the beads. The beads then formed a pellet on the magnetic stand and was washed two times using 80% ethanol. Purified PB product was amplified using Bio-Rad CFX96 Touch real-time PCR detection system for 12 cycles. Post amplifi-cation, the total volume for the reaction was measured, and the DNA was diluted to a total volume of 205 lL with 10mM Tris. Two hundred microliters of the diluted product was transferred to new 1:5-mL LoBind tubes for purification. Five microliters of the unpurified product was stored as control. 0.8 equivalents of MAGwise paramagnetic beads was added to the PB product, and DNA was allowed to bind to the beads. The beads were allowed to form a pellet on the inner wall of the tube. The pellet was washed two times using 80% ethanol. Purified library was eluted using 32 lL of 10 mM of Tris. Twenty-eight microliters of the purified product was transferred to a new 1:5-mL LoBind tube. Four libra-ries were prepared, each having a unique Pool Barcode. Purified libraries were stored at - 20°C for short-term storage. Library qual-ity control (QC) was done on the Agilent Bioanalyzer (High-Sensitivity DNA 5,000 kit) at the Advanced Genomics Core.

RNA-Sequencing and Data Processing

Libraries were sequenced on the Illumina NovaSeq 6000 at the University of Michigan Advanced Genomics Core. FASTQ reads were demultiplexed back into individual samples based on the i7 and i5 indices. Sequencing data were transferred to the University of Michigan Great Lakes high-performance computing cluster for analysis. Sequencing read quality was assessed via FastQC and MultiQC Reads were aligned to a splice junction aware build of the human genome (GRCh38) using STAR. Aligned reads were assigned to genes using featureCounts, where multimapping and multi-overlapping reads were not counted.

Differential Gene Expression Testing

Read count matrices from featureCounts were loaded into edgeR, and samples with fewer than 1,000,000 mapped reads were excluded from downstream analysis. Genes with low expression were excluded from analysis using the edgeR filterByExpr() function with default settings. Normalization factors and disper -sion were calculated prior to generating a log2-transformed nor-malized counts per million (cpm) matrix for downstream analysis. Differential gene expression between each treatment (average of four replicates) and the relevant controls (average of four replicates) was calculated using quasi-likelihood negative binomial generalized log-linear modeling in edgeR. Genes were considered differentially expressed between a treatment and control at a false discovery rate (FDR) adjustedp-value of <0:05.

Benchmark Concentration Analyses

Gene-specific benchmark concentration and best fit benchmark concentration model were identified using BMDExpress version 3.2, a free software package developed by the National Institute of Environmental Health Sciences (NIEHS) and available for download (https://github.com/auerbachs/BMDExpress-3/releases). Best practices for dose-response modeling for each chemical were con-ducted according BMDExpress online documentation (https:// bmdexpress-2.readthedocs.io/en/feature-readthedocs/basic_ workflow/). Normalized counts per million reads were imported into BMDExpress and prefiltered using one-way analysis of var-iance (ANOVA) for significance at ap-value of <0:05 to identify genes showing significant increasing or decreasing concentration responses. To determine concentration-response relationships, the filtered data were modeled in BMDExpress with hill, power, linear, polynomial (2°, 3°, 4°), and exponential (3 and 5) models and the best fit model was chosen. The benchmark response (BMR) was set to 1 standard deviation relative to control response, maximum itera-tions of the model was set to 250, a confidence level of 0.95 was used, and power was restricted to > 1 for the power model. Best fit models were chosen for each concentration-response relationship per gene using nested chi square to select for the best model followed by lowest Akaike information criterion (AIC). Hill models were flagged if its "k" parameter was smaller than one-third of the lowest positive concentration and flagged hill models were included in the data. Directionality in response (upregulated/downregulated) was derived from the "Best adverseDirection" of the winning model for a given gene. We conducted "Defined Category" analyses to additionally assess the concentration-response relationship to the Molecular Signature Database "hallmark" gene sets,23 embryonic stem cell genes,24 breast cancer molecular signatures,25 and breast stem cell genes.26 Benchmark concentrations with values above the highest tested concentration (25 lM) were excluded from individual gene and functional classification analyses. Defined category results were filtered at a significance at a Fisher's exact two-tailed p-value of <0:05 for BMDExpress testing, and results were visualized by heatmap using the R package pheatmap and heatmap.2 function from the R package gplots. All original expression data, ANOVA filtered results, BMC results, and subsequent BMDExpress analysis for all chemicals are included in a .bm2 file available in supplemen-tary information.

For hypergeometric enrichment testing of BMC data, we used the R package HypeR for hypergeometric enrichment testing for two user-defined gene sets (one Breast Cancer Screen or "BCScreen" gene set and gene set from Pal et al.).27-29 For hypergeometric enrichment testing, background gene expression included all human genes (23467). Hypergeometric enrichment results were filtered for significance at an FDR of <0:05.

Comparison ofNHANES Exposure Levels to Benchmark Concentrations

We compared bioactive concentrations of the assayed chemicals to human biomarker concentrations measured as part of the National Health and Nutrition Examination Survey (NHANES) following our previously established methodology.20 NHANES is designed to understand the health status of the US population, including a substantial number of chemical biomarker measure-ments in urine, blood, and serum. For this analysis, we used our curated NHANES dataset, which contains information on chemical biomarker concentrations of the 21 assessed chemicals in up to 57,786 women recruited from 1999-2018.30 Linkage between chemicals used in vitro and their corresponding exposure bio-markers in NHANES are shown in Table S2. Chemical biomarker concentrations were converted to molarity units and then compared to the benchmark concentrations identified through the benchmark concentration (BMC) analyses of the RNA-seq data using boxplots. These comparisons are designed to help us understand whether dysregulated gene expression due to toxicant exposure may be occurring at population-relevant levels. Overlapping NHANES biomarker concentrations and BMC concentrations indicate that chemical biomarkers measured in US women in NHANES occur within the estimated benchmark concentrations.

Secondary Analysis ofNormal Breast Single Cell Data for Cell Type Composition Prediction

To quantify whether treatments resulted in differences in the cel-lular composition of our samples, we performed bioinformatic deconvolution of our bulk RNA-seq data based on a single cell atlas of the normal human mammary gland,28 where 10 normal mammary gland samples were processed for single cell analysis using the 10 × Genomics Chromium platform. Sample-specific count matrices of the single cell RNA-seq data were downloaded from the Gene Expression Omnibus and processed via Seurat (version 4.1). Briefly, we excluded cells that expressed fewer than 200 or more than 5,000 genes, along with cells with >15% mitochondrial reads. Data were log-normalized and scaled and then clustered with 0.02 resolution, following Pal et al.'s meth-ods, to identify four cell clusters.28 Clusters were then assigned to "luminal progenitor" (CD49f+EpCAM+), "mature luminal" (CD49f"EpCAM+), and "myoepithelial" (CD49f+EpCAMlo/_), identities based on marker gene expression assigned by Pal et al.28 Gene expression signatures for each of the cell types were calcu-lated using the FindAllMarkers() function in Seurat.

We then used these single cell gene expression data to estimate cell type proportions in our data using the multisubject single-cell (MuSiC) deconvolution method version 1.0.0 in R.31 MuSiC pre-dicts cell type proportions in bulk RNA-seq data based on anno-tated single cell RNA expression data from multiple donors. MuSiC weights genes, which show cross-subject and cross-cell type consistency in the single cell RNA-seq data to then transfer cell type-specific gene expression information from single cell RNA-seq datasets to bulk RNA-seq datasets. Cell type proportions in bulk RNA-seq data are predicted using a nonnegative least squares regression-based approach with these weighted gene expression signatures, with the constraints that cell type propor-tions cannot be negative and must sum to 1. Using MuSiC, we esti-mated the cell type proportions in each bulk RNA-seq sample to test the hypothesis that chemical treatments may lead to cellular plasticity or differences in cell type marker expression. We com-pared two MuSiC estimates for our data: one where proportions of "mature luminal," "luminal progenitor," and "myoepithelial" cells were estimated and one where only "luminal progenitor" and "myoepithelial" cells were estimated. Based on model fit statistics (R2), the model prediction, which only included "luminal progeni-tor" and "myoepithelial" cells, was better (median R2 = 0:29 vs. median R2 = 0:22 when including the three cell populations). The decreased model performance when including "mature luminal" cells also aligns with biological expectations of the MFC10A cell line, which, as an immortalized line, is unlikely to contain a significant proportion of mature luminal cells when grown in two dimen-sions (2D). Thus, we moved forward with predicting proportions of "luminal progenitor" and "myoepithelial" proportions in our bulk RNA-seq data. To visualize the differences in predicted cell type proportions for each chemical by concentration were then visualized using ggboxplot() in R, with differences between each concentration and control for a given chemical tested using a t-test via the stat_compare_means() function from the ggpubr() package.

Data Availability

The sequencing data have been deposited at the Gene Expression Omnibus (accession number GSE220051).

Results

To understand the potential impacts of our chemicals of interest on dysregulated breast cancer biology, we treated nontumorigenic MCF10A cells with chemicals that represent exposures from metals, pesticides, personal care products, and other industrial and environmental sources, including per- and polyfluoroalkyl sub-stances (see Methods and Table S1 for a list of the chemicals used). Although the chemicals differ in cytotoxicity, we treated with the same four concentrations (25 nM, 250 nM, 2:5 lM, 25 lM) of each chemical for 48 h for consistency and to obtain a broad understanding of potential bioactivity. Genes were consid-ered differentially expressed between a treatment and control at a false discovery rate (FDR) adjusted p-value of <0:05. Different chemicals, particularly arsenic, cadmium, thiram, and DDE showed unique gene expression responses across concentrations, suggestive of unique transcriptional effects, even at high concen-trations (Figure S1 and S2). Fold change and significance of dif-ferential gene expression varied substantially by chemical and by concentration (Figure 1 and Table 1; Figure S2 and Table S3). Sequencing data can be found at the Gene Expression Omnibus (accession number GSE220051). For example, p,p0-DDE showed over 700 differentially expressed genes (DEGs) at all four con-centrations, while mercury only had more than 10 DEGs at the highest concentration of 25 lM. Arsenic and cadmium both had DEGs at all four concentrations, but had more than 8,000 and 4,000 DEGs at 25 lM, respectively. As an example of the differ-ences in effect of each chemical at a given concentration, volcano plots showing the results of differential expression analysis com-paring the 2:5 lM concentration to control are shown in Figure 1, with the comparisons for the other concentrations compared to con-trol included in Figure S2.

Genes related to breast cancer were commonly found to be differentially expressed by this panel of chemicals at various con-centrations. Table S4 lists genes differentially expressed across the 21 chemicals assessed, sorted by number of chemicals where there were expression differences in at least one of the four con-centrations tested. Of note, serum amyloid alpha 1 and 2 (SAA1 and SAA2) showed expression differences compared to controls result-ing from 14 and 15 chemical exposures, respectively. S100A8 and S100A9 cytokines were also differentially expressed after treatment with 16 and 11 chemicals, respectively. G protein alpha subunit (GNAS) expression differences were observed across 13 of the 21 chemicals. Additionally, aldo-keto reductases AKR1C1, AKR1C2, and AKR1C3 expression were different across 15, 12, and 9 chemicals, respectively. Keratins KRT14, KRT15, and KRT17 expression levels were different across 5, 6, and 11 chemicals. Serine protease inhibitor, clade E member 1 (SERPINE1) also showed differences in expression across nine chemicals.

To characterize the concentration-response relationship of each exposure on gene expression, we performed benchmark concentra-tion modeling to define the range of benchmark concentrations for the different chemicals (Figure 2). Best benchmark concentrations were estimated from a set of concentration-response models for each gene for each chemical (Figure 2A). The total number of genes filtered by ANOVA and fit to benchmark models was different for each chemical, with some chemicals (arsenic, cadmium, and thiram) inducing concentration-responsive differences in over 1,000 different genes. Concentration-response fit was assessed via examination of the best fit concentration-response curves for individual gene/ chemical combinations. Two example gene BMC curves are dis-played in Figure 2B,C, showing the gene expression of ARC follow-ing exposure to a range of arsenic concentrations in a linear best fit model (Figure 2B) and thiram's concentration-response regulation of gene HSPA6 in a hill best fit model (Figure 2C) across the range of treated concentrations. All benchmark concentration response data for each dysregulated gene for every chemical tested is included in Table S5 and in the supplementary bm2 file. Of the BMC modeled data, 14 genes were dysregulated by 10 or more chemicals, 1,239 were dysregulated by 5 or more chemicals, and many of the genes dysregulated by multiple chemicals often corre-sponded tothe same directionalityof response (Table S5).

We next wanted to contextualize the concentrations required in vitro to dysregulate gene expression with biomarker levels detectable in women in the United States. We compared bench-mark concentration results from our in vitro experiments to National Health and Nutrition Examination Survey (NHANES) biomarker concentration data for each chemical, plotting the range of blood or urine chemical biomarker concentrations in US women compared to the range of BMCs for each chemical from our in vitro data (Figure 3). For all chemicals, there were at least a few genes predicted to change within the range of human popu-lation exposures. For highly exposed individuals in NHANES, biomarker concentrations for some chemicals (arsenic, copper, lead, p,p'-DDE, thiram, 2,5-dichlorophenol, methyl paraben, and propyl paraben) exceeded the median predicted benchmark concentrations derived from the RNA-seq data.

We next wanted to evaluate the differences in biological proc-esses after chemical exposure and characterize the concentrations at which these processes may be dysregulated. BMC modeled genes were assessed using MSigDB's hallmark gene sets consist-ing of 50 biological processes. The supplemental bm2 file contains benchmark concentration gene data for every hallmark category for each chemical organized under "functional classifications." Figure 4A displays an accumulation plot depicting the median BMC at which hallmark gene sets were impacted (at least one gene in that gene set with effects significantly different from control). The accumulation plot demonstrates low concentration effects present for many of the hallmarks. Noticeably, this accumulation plot demonstrating pathway effects (Figure 4A) differs from the accumulation plot of gene expression results (Figure 2A), high-lighting potential patterns in gene expression results. Figure 4B shows a heatmap of significantly different hallmark processes-20 of the chemicals showed significant enrichment or downregulation for 1 or more of the hallmark gene sets (filtered for significant enrichment, whereas Figure 4A is not filtered for significance). Notably, many of these processes are related to carcinogenesis, including the epithelial mesenchymal transition and reactive oxy-gen species pathway. Compared to controls, genes related to epithelial to mesenchymal transition were upregulated by bisphenol A (BPA), bisphenol F (BPF), bisphenol S (BPS), copper, lead, methyl paraben, PFOA, and thiram at median BMCs of 2:56 lM, 3:21 lM, 5:23 lM, 0:32 lM, 0:38 lM, 0:64 lM, 1:54 lM, and 7:15 lM, respectively. Genes related to reactive oxygen species were upregulated by BPF, lead, and mercury at median BMCs of 11:19 lM, 0:24 lM, and 1:81 lM, respectively. Copper, p,p0-DDE, lead, PFDA, PFNA, and phenanthrene upregulated cell cycle- related genes, including E2F targets and G2M checkpoint (median BMCs of 2:11 lM, 0:02 lM, 1:12 lM, 0:28 lM, 1:54 lM, and 2:09 lM for E2F target genes, respectively, and 3:21 lM, 0:02 lM, 1:09 lM, 0:30 lM, 1:65 lM, and 2:25 lM for G2M checkpoint genes, respectively). E2F targets were downregulated by arsenic, cadmium, and thiram (median BMCs of 5:00 lM, 8:38 lM, and 7:16 lM, respectively), and G2M checkpoint genes were downregu-lated by cadmium and thiram at median BMCs of 8:22 lM and 8:30 lM, respectively, highlighting the potential for distinct mecha-nisms of promoting carcinogenesis depending on the chemical involved (BMC data for hallmark processes in Table S7).

To assess specific pathways related to breast cancer, additional enrichment analyses assessed additional breast cancer-associated bi-ological pathways in MSigDB's Chemical and Genetic Perturbations (CGP) gene sets (Figure 4C): gene sets from three studies related to breast cancer and stemness.25-27 Nineteen chemicals differentially regulated these gene sets at various concentrations. Genes up-regulated in embryonic stem cells were upregulated by 1,4-dichlorobenzene, p,p0-DDE, PFDA, PFNA, and PFOS (median BMCs of 7:75 lM, 0:02 lM, 1:43 lM, 1:64 lM, and 1:82 lM). Additionally, genes upregulated in basal subtypes of breast can-cer were found upregulated by BPF, copper, and 1,4-dichloro-benzene (median BMCs of 3:18 lM, 1:38 lM, and 2:50 lM, respectively). BPA, lead, and PCB187 also upregulated genes that are upregulated in normal-like breast cancers (median BMCs of 10:79 lM, 1:34 lM, and 8:35 lM, respectively). (BMC data for CGP processes inTable S7.)

To assess potential breast cancer-specific gene expression dif-ferences regulated by our chemicals of interest, we tested for gene enrichment with a publicly available breast carcinogenesis gene panel, BCScreen.29 The overlapping genes within different categories of breast carcinogenesis are shown in Table S8. The gene panel included 500 genes total divided into 14 categories, with 33 genes per category other than the "mammary" category, which includes 71 genes. Of these, cells exposed to nine chemi-cals (arsenic, cadmium, copper, p,p'-DDE, thiram, 2,5-dichloro-phenol, PFDA, PFNA, and BPF) showed differential expression of at least three genes in one or more categories, with cell cycle, genotoxicity, and mammary categories showing the largest num-ber of genes dysregulated.

We next wanted to test the hypothesis that exposure to the chemicals under assessment result in differences consistent with a luminal-to-basal transition, a process likely important in the devel-opment of basal breast cancers. We performed cell type deconvolu-tion based on a normal breast single cell RNA-seq reference28 to predict the proportions of myoepithelial and luminal progenitor cells in each of our exposures (Figure 5; Figures S3 and S4).28 The initial proportion of cell types varied between controls as the cell population of MCF10As is heterogeneous. Only cells exposed to PFNA showed significantly lower predicted myoepithelial cell proportions at the highest concentration of 25 lM (Figure S3). Cells exposed to arsenic and lead showed a concentration-dependent significantly higher proportion of myoepithelial cells and corresponding lower proportion of luminal progenitor cells across the three highest concentrations compared to control; those exposed to copper showed a significantly higher proportion of myoepithelial cells across all concentrations, and those exposed to methyl paraben showed a significantly higher myoepithelial cell proportion at all concentrations other than 2:5 lM (Figure 5). To assess for differences in proliferation profiles of chemically treated cells, expression of three cell cycle-associated genes (PCNA, MKI67, and MCM2) was assessed, and it was found that four chemical treatments (arsenic, copper, p,p'-DDE, and lead)

caused upregulation and one chemical (cadmium) caused down-regulation of PCNA, four chemicals (arsenic, copper, p,p0-DDE, and PFDA) caused upregulation and one chemical (cadmium) caused downregulation of MKI67, and one chemical (p,p0-DDE) caused upregulation and three chemicals (arsenic, PFOS, and thiram) caused downregulation of MCM2 (Table S10).

Discussion

Across all 21 chemicals representing exposures commonly affect-ing women in the United States,18 12,609 genes showed expression differences in exposed breast cells compared to controls.

While many of these differences were found at the highest con-centration, lower concentration effects were also seen, with BMCs falling within or near the range of NHANES biomarker levels for all chemicals. This suggests that the gene expression effects seen in this study may be biologically relevant to human exposure levels, and so our analysis may link chemical exposures to incidence and characteristics of breast cancer, although further study is needed to establish human relevance in addition to our model. Interestingly, of the chosen set of chemicals, only ten (arsenic, cadmium, lead, p,p0-DDE, 1,4-dichlorobenzene, 2,5-dichlorophenol, PFOA, PFOS, PCB187, and PCB153) have been classified by International Agency for

Research on Cancer (IARC) as possible, probable, or definite carcinogens to humans.32 This study adds new data on the concentration-dependent effects of this suite of common expo-sures on nontumorigenic breast cells and identifies potential relevance to breast carcinogenesis.

Specific gene differences included some commonalities across chemical treatments, and many of these were related tobreast cancer progression. This included SAA1 and SAA2, S100A8 and S100A9, GNAS, aldo-keto reductases, keratins, and SERPINE1. SAA1 has been found to be associated with the accumulation of suppressive neutrophils and SAA1 and SAA2 have been found to be associated with the progression of breast cancer as well as poor prognosis, par-ticularly in triple-negative breast cancer.33-35 Overexpression of proinflammatory cytokines S100A8 and S100A9 in MCF10As were associated with reprogramming of normal breast cells (MCF10As) into tumorigenic cells.36,37 GNAS has been implicated in breast cancer cell proliferation and epithelial to mesenchymal transition via modulation of the PI3K/AKT/Snail1/E-cadherin signaling axis.38 Aldo-keto reductases were shown to be involved in estrogen and progesterone production and thus involved in pro-proliferative signaling of hormone responsive tumors.39 Dysregulation of keratin expression levels has also been associated with breast cancer, with increased KRT14 expression associated with invasiveness, decreased KRT15 expression associated with poor prognosis, and increased KRT17 associated with the develop-ment of TNBC.40-42 SERPINE1 was found to be associated with re-sistance of TNBC to paclitaxel and knockdown of SERPINE1 reversed resistance of breast cancer to paclitaxel through downregu-lation of vascular endothelial growth factor A (VEGFA).43

In addition to common gene expression differences across chemical treatments, our enrichment analyses suggested that chemicals may alter breast cell biology by several distinct mecha-nisms. Treatment with metals including copper, lead, and mer-cury were found to result in enrichment of genes associated with epithelial to mesenchymal transition and reactive oxygen species pathways, consistent with previous literature on metal-induced carcinogenesis that places particular emphasis on oxidative stress- related carcinogenesis.44-46 In breast cancer, arsenic and cadmium have been associated with oxidative stress-related pathways of malignant transformation.47 While these metals as well as many of the other chemicals in our analysis have been shown to play an endocrine-disrupting role that could lead to breast tumorigenesis, these findings propose alternative routes of carcinogenesis.22,48,49 p,p0-DDE upregulated gene signatures related to cell cycle proc-esses in our study. Previous literature has found that a mixture of organochlorine pesticides including p,p0-DDE had potential to modulate proliferation of breast cancer cell lines partially through induction of cell cycle entry.50

Differences in stemness features were also seen across cells exposed to fifteen of the chemicals analyzed, important because the cancer stem cell hypothesis states that there is a small subset of cancer stem cells that result in cancer initiation and invasion, and so differences in stemness features can lead to increased tumorigenicity and metastasis.51 We found stemness signatures enhanced by 1,4-dichlorobenzene, p,p0-DDE, PFDA, PFNA, and PFOS, findings that could indicate novel mechanisms of chemi-cally induced carcinogenesis.

These findings may merit consideration in the context of breast cancer disparities. We previously identified that in the US, non-Hispanic black women were disproportionately exposed to many of the chemicals assayed here.18 For example, non-Hispanic black women, on average, had 140% higher levels of urinary methyl paraben, 75% higher levels of blood p,p0-DDE, 93% higher levels of urinary bisphenol S, 39% higher levels of blood mercury, and 17% higher levels of blood lead compared to non-Hispanic white women.18 Non-Hispanic black women were also three times more likely to be diagnosed with basal-like triple-negative breast cancers compared to white women.52,53 Triple-negative breast cancers of-ten have a basal-like phenotype and may derive from a luminal to basal cell state transition.13-15 Here, using enrichment analysis, we found upregulation of basal cancer subtype genes in cells exposed to copper, 1,4-dichlorobenzene, and BPF, and cell type proportion analysis found increased predicted myoepithelial proportions (through examination of gene expression patterns consistent with myoepithelial vs. luminal features) by arsenic, copper, lead, and methyl paraben anda lower basal cell proportion at the highest con-centration of PFNA. The findings for PFNA here also indicate an area for future study to understand potential divergence in predictions by different methods. We further assessed proliferation markers in order to determine whether an increase in myoepithelial proportion was indicative of a change in proliferation, and we found that for some chemicals that induced a higher proportion of myoepithelial markers, such as arsenic, there was an increase in proliferation markers, while for others such as methyl paraben, there was no change in proliferation markers (Table S10). Previous literature has shown that arsenic is capable of inducting transition from luminal to basal features, but this association has not been studied in detail for the other chemicals.54 These findings highlight an increased need to examine chemical exposures as a modifiable risk factor for aggressive breast cancers, particularly in the context of breast cancer disparities. One major limitation of the present study is that all exposures were done on an acute timescale (48 h) while actual human exposures to the chemicals analyzed here are likely chronic, over the course of years or a lifetime. Differences in the vehicle chemicals were dissolved in and compared (water vs. DMSO), which represents another limitation, as DMSO could be affecting gene expression. Another limitation is that only one cell type was used for this study, nontumorigenic immortalized breast cells. Ongoing work is considering treatment on a more chronic timescale as well as using different breast cell lines representing diverse donors.

Additionally, MCF10As are an imperfect model of human breast cancer, as a cell line does not capture the complex milieu of breast tissue. Future work could consider treatment ofdiverse breast cancer cell lines, ranging from more basal to more luminal tumor types to investigate differences in luminal or basal features related to chemical exposure. Future studies can also assess effects of a broader range of concentrations of each chemical to further under-stand effects across relevant ranges, including low concentration effects. Another limitation in our analyses integrating benchmark concentration levels with NHANES biomarker concentrations is that the NHANES chemical biomarkers are measured in urine and blood. These levels may or may not reflect mammary tissue levels. As additional research is conducted on biomonitoring of mammary tissues or physiologically based toxicokinetic models reflecting mammary gland concentrations are developed, these data and meth-ods will substantially inform our understanding of how chemical exposures impact breast cancer risk.

Overall, our transcriptomic analysis of chemically treated normal breast cells uncovers molecular mechanisms that may be associated with more aggressive forms of breast cancer. This includes both features related to known carcinogenic proper-ties of these chemicals as well as potentially novel mechanisms, opening the door to future avenues of investigation to elucidate how different exposures relate to breast carcinogenesis and breast cancer disparities.

Acknowledgments

This work was supported by grants from the National Institutes of Health (R01 ES028802, R01 AG072396, T32 ES 007062, P30 ES017885, P30 CA046592).

References

References

1. Arnold M, Morgan E, Rumgay H, Mafra A, Singh D, Laversanne M, et al. 2022. Current and future burden of breast cancer: global statistics for 2020 and 2040. Breast 66:15-23, PMID: 36084384, https://doi.org/10.1016/j.breast.2022.08.010.

2. Johnson KS, Conant EF, Soo MS. 2021. Molecular subtypes of breast cancer: a review for breast radiologists. J Breast Imaging 3:12-24.

3. Lord SJ, Bahlmann K, O'Connell DL, Kiely BE, Daniels B, Pearson SA, et al. 2022. De novo and recurrent metastatic breast cancer-a systematic review of population-level changes in survival since 1995. EClinicalMedicine 44:101282, PMID: 35128368, https://doi.org/10.1016/j.eclinm.2022.101282.

4. Rodgers KM, UdeskyJO, Rudel RA, BrodyJG.2018. Environmental chemicals and breast cancer: an updated review of epidemiological literature informed by bio-logical mechanisms. Environ Res 160:152-182, PMID: 28987728, https://doi.org/10. 1016/j.envres.2017.08.045.

5. Hanahan D, Weinberg RA. 2000. The hallmarks of cancer. Cell 100(1):57-70, PMID: 10647931, https://doi.org/10.1016/s0092-8674(00)81683-9.

6. Hanahan D, Weinberg RA. 2011. Hallmarks of cancer: the next generation. Cell 144(5):646-674, PMID: 21376230, https://doi.org/10.1016/j.cell.2011.02.013.

7. Hanahan D. 2022. Hallmarks of cancer: new dimensions. Cancer Discov 12(1):31 46, PMID: 35022204, https://doi.org/10.1158/2159-8290.CD-21-1059.

8. Colacino JA, Azizi E, Brooks MD, Harouaka R, Fouladdel S, McDermott SP, et al. 2018. Heterogeneity of human breast stem and progenitor cells as revealed by transcriptional profiling. Stem Cell Reports 10(5):1596-1609, PMID: 29606612, https://doi.org/10.1016/j.stemcr.2018.03.001.

9. Thong T, Wang Y, Brooks MD, Lee CT, Scott C, Balzano L, et al. 2020. Hybrid stem cell states: insights into the relationship between mammary development and breast cancer using single-cell transcriptomics. Front Cell Dev Biol 8:288, PMID: 32457901, https://doi.org/10.3389/fcell.2020.00288.

10. Spike BT, Engle DD, Lin JC, Cheung SK, La J, Wahl GM. 2012. A mammary stem cell population identified and characterized in late embryogenesis reveals similarities to human breast cancer. Cell Stem Cell 10(2):183-197, PMID: 22305568, https://doi.org/10.1016/j.stem.2011.12.018.

11. Giraddi RR, Chung C-Y, Heinz RE, Balcioglu O, Novotny M,Trejo CL, et al. 2018. Single-cell transcriptomes distinguish stem cell state changes and lineage specification programs in early mammary gland development. Cell Rep 24(6):1653-1666, PMID: 30089273, https://doi.org/10.1016/j.celrep.2018.07.025.

12. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. 2018. Machine learning identifies stemness features associated with onco-genic dedifferentiation. Cell 173(2):338-354, PMID: 29625051, https://doi.org/10. 1016/j.cell.2018.03.034.

13. Dent R, Trudeau M, Pritchard KI, Hanna WM, Kahn HK, Sawka CA, et al. 2007. Triple-negative breast cancer: clinical features and patterns of recurrence. Clin Cancer Res 13(15):4429-4434, PMID: 17671126, https://doi.org/10.1158/1078-0432.CCR-06-3045.

14. Lim E, Vaillant F, Wu D, Forrest NC, Pal B, Hart AH, et al. 2009. Aberrant luminal progenitors as the candidate target population for basal tumor development in BRCA1 mutation carriers. Nat Med 15(8):907-913, PMID: 19648928, https://doi.org/ 10.1038/nm.2000.

15. Rädler PD, Wehde BL, Triplett AA, Shrestha H, Shepherd JH, Pfefferle AD, et al. 2021. Highly metastatic claudin-low mammary cancers can originatefrom luminal epithelial cells. Nat Commun 12(1)3742, PMID: 34145248, https://doi.org/10.1038/ s41467-021-23957-5.

16. CDC (Centers for Disease Control and Prevention). 2022. National Report on Human Exposure to Environmental Chemicals. https://www.cdc.gov/ exposurereport/[accessed31 January2023].

17. Nguyen VK, Colacino J, Patel CJ, Sartor M, Jolliet O. 2022. Identification of occupations susceptible to high exposure and risk associated with multiple toxicants in an observational study: national health and nutrition examination survey 1999-2014. Exposome 2(1):osac004, PMID: 35832257, https://doi.org/10. 1093/exposome/osac004.

18. Nguyen VK, Kahana A, Heidt J, Polemi K, Kvasnicka J, Jolliet O, et al. 2020. A comprehensive analysis of racial disparities in chemical biomarker concentra-tions in United States women, 1999-2014. Environ Int 137:105496, PMID: 32113086, https://doi.org/10.1016/j.envint.2020.105496.

19. Nguyen VK, Colacino JA, Arnot JA, Kvasnicka J, Jolliet O. 2019. Characterization of age-based trends to identify chemical biomarkers of higher levels in children. Environ Int 122:117-129, PMID: 30528102, https://doi.org/10.1016/j.envint.2018.10. 042.

20. Polemi KM, Nguyen VK, Heidt J, Kahana A, Jolliet O, Colacino JA. 2021. Identifying the link between chemical exposures and breast cancer in African American women via integrated in vitro and exposure biomarker data. Toxicology 463:152964, PMID: 34600088, https://doi.org/10.1016/j.tox.2021.152964.

21. Zeinomar N, Oskar S, Kehm RD, Sahebzeda S, Terry MB. 2020. Environmental exposures and breast cancer risk in the context of underlying susceptibility: a systematic review of the epidemiological literature. Environ Res 187:109346, PMID: 32445942, https://doi.org/10.1016/j.envres.2020.109346.

22. Wan MLY, Co VA, El-Nezami H. 2022. Endocrine disrupting chemicals and breast cancer: a systematic review of epidemiological studies. Crit Rev Food Sci Nutr 62(24):6549-6576, PMID: 33819127, https://doi.org/10.1080/10408398. 2021.1903382.

23. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. 2015. The molecular signatures database (MsigDB) hallmark gene set collec-tion. Cell Syst 1(6):417-425, PMID: 26771021, https://doi.org/10.1016/j.cels.2015. 12.004.

24. Wong DJ, Liu H, Ridky TW, Cassarino D, Segal E, Chang HY. 2008. Module map of stem cell genes guides creation of epithelial cancer stem cells. Cell Stem Cell 2(4):333-344, PMID: 18397753, https://doi.org/10.1016/j.stem.2008.02.009.

25. Smid M, Wang Y, Zhang Y, Sieuwerts AM, Yu J, Klijn JGM, et al. 2008. Subtypes of breast cancer show preferential site of relapse. Cancer Res 68(9):3108-3114, PMID: 18451135, https://doi.org/10.1158/0008-5472.CAN-07-5644.

26. Lim E, Wu D, Pal B, Bouras T, Asselin-Labat M-L, Vaillant F, et al. 2010. Transcriptome analyses of mouse and human mammary cell subpopulations reveal multiple conserved genes and pathways. Breast Cancer Res 12(2):R21, PMID: 20346151, https://doi.org/10.1186/bcr2560.

27. Federico A, Monti S. 2020. hypeR: an R package for geneset enrichment work-flows. Bioinformatics 36(4):1307-1308, PMID: 31498385, https://doi.org/10.1093/ bioinformatics/btz700.

28. Pal B, Chen Y, Vaillant F, Capaldo BD, Joyce R, Song X, et al. 2021. A single-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast. Embo J 40(11):e107333, PMID: 33950524, https://doi.org/10. 15252/embj.2020107333.

29. Grashow RG, De La Rosa VY, Watford SM, Ackerman JM, Rudel RA. 2018. BCScreen: a gene panel to test for breast carcinogenesis in chemical safety screening. Comput Toxicol 5:16-24, PMID: 31218268, https://doi.org/10.1016/j. comtox.2017.11.003.

30. Nguyen VK, Middleton L, Huang L, Verly E, Kvasnicka J, Sagers L, et al. 2023. NHANES 1988-2018. Kaggle. https://doi.org/10.34740/KAGGLE/DSV/ 4856629 [accessed 19 December 2023].

31. Wang X, Park J, Susztak K, Zhang NR, Li M. 2019. Bulk tissue cell type decon-volution with multi-subject single-cell expression reference. Nat Commun 10(1):380, PMID: 30670690, https://doi.org/10.1038/s41467-018-08023-x.

32. WHO (World Health Organization). 2024. List of Classifications - IARC Monographs on the Identification of Carcinogenic Hazards to Humans. https:// monographs.iarc.who.int/list-of-classifications/ [accessed 30 November 2023].

33. Niu X, Yin L, Yang X, Yang Y, Gu Y, Sun Y, et al. 2022. Serum amyloid a 1 indu-ces suppressive neutrophils through the toll-like receptor 2-mediated signaling pathway to promote progression of breast cancer. Cancer Sci 113(4):1140- 1153, PMID: 35102665, https://doi.org/10.1111/cas.15287.

34. Yang M, Liu F, Higuchi K, Sawashita J, Fu X, Zhang L, et al. 2016. Serum amy-loid a expression in the breast cancer tissue is associated with poor progno-sis. Oncotarget 7(24):35843-35852, PMID: 27058895, https://doi.org/10.18632/ oncotarget.8561.

35. Ignacio RMC, Gibbs CR, Kim S, Lee ES, Adunyah SE, Son DS. 2019. Serum amy-loid a predisposes inflammatory tumor microenvironment in triple negative breast cancer. Oncotarget 10(4):511-526, PMID: 30728901, https://doi.org/10. 18632/oncotarget.26566.

36. Jo SH, Heo WH, Son H-Y, Quan M, Hong BS, Kim JH, et al. 2021. S100A8/A9 mediate the reprograming of normal mammary epithelial cells induced by dynamic cell-cell interactions with adjacent breast cancer cells. Sci Rep 11(1):1337, PMID: 33446797, https://doi.org/10.1038/s41598-020-80625-2.

37. Song R, Struhl K. 2021. S100A8/S100A9 cytokine acts as a transcriptional coac-tivator during breast cellular transformation. Science Advances 7(1):eabe5357, https://doi.org/10.1126/sciadv.abe5357.

38. Jin X, Zhu L, Cui Z, Tang J, Xie M, Ren G. 2019. Elevated expression of GNAS promotes breast cancer cell proliferation and migration via the PI3K/AKT/ Snail1/E-cadherin axis. Clin Transl Oncol 21(9):1207-1219, PMID: 30767161, https://doi.org/10.1007/s12094-019-02042-w.

39. Penning TM, Byrns MC. 2009. Steroid hormone transforming aldo-keto reduc-tases and cancer. Ann NY Acad Sci 1155:33-42, PMID: 19250190, https://doi.org/ 10.1111/j.1749-6632.2009.03700.x.

40. Zhong P, Shu R, Wu H, Liu Z, Shen X, Hu Y. 2021. Low KRT15 expression is associated with poor prognosis in patients with breast invasive carcinoma. Exp Ther Med 21(4):305, PMID: 33717248, https://doi.org/10.3892/etm.2021. 9736.

41. Zhang H, Zhang Y, Xia T, Lu L, Luo M, Chen Y, et al. 2022. The role of Keratin17 in human tumours. Front Cell Dev Biol 10:818416, PMID: 35281081, https://doi.org/ 10.3389/fcell.2022.818416.

42. Hanley CJ, Henriet E, Sirka OK, Thomas GJ, Ewald AJ. 2020. Tumor resident stromal cells promote breast cancer invasion through regulation of the basal phenotype. Mol Cancer Res 18(11):1615-1622, PMID: 32868298, https://doi.org/ 10.1158/1541-7786.MCR-20-0334.

43. Zhang Q, Lei L, Jing D. 2020. Knockdown of SERPINE1 reverses resistance of triple-negative breast cancer to paclitaxel via suppression of VEGFA. Oncol Rep 44(5):1875-1884, PMID: 33000256, https://doi.org/10.3892/or.2020.7770.

44. Harris GK, Shi X. 2003. Signaling by carcinogenic metals and metal-induced reac-tive oxygen species. Mutat Res 533(1-2):183-200, PMID: 14643420, https://doi.org/ 10.1016/j.mrfmmm.2003.08.025.

45. Eckstein M, Rea M, Fondufe-Mittendorf YN. 2017. Transient and permanent changes in DNA methylation patterns in inorganic arsenic-mediated epithelial-to-mesenchymal transition. Toxicol Appl Pharmacol 331:6-17, PMID: 28336213, https://doi.org/10.1016/j.taap.2017.03.017.

46. Chen QY, DesMarais T, Costa M. 2019. Metals and mechanisms of carcinogen-esis. Annu Rev Pharmacol Toxicol 59:537-554, PMID: 30625284, https://doi.org/ 10.1146/annurev-pharmtox-010818-021031.

47. Zimta A-A, Schitcu V, Gurzau E, Stavaru C, Manda G, Szedlacsek S, et al. 2019. Biological and molecular modifications induced by cadmium and arsenic

during breast and prostate cancer development. Environ Res 178:108700, PMID: 31520827, https://doi.org/10.1016/j.envres.2019.108700.

48. Davey JC, Bodwell JE, Gosse JA, Hamilton JW. 2007. Arsenic as an endocrine disruptor: effects of arsenic on estrogen receptor-mediated gene expression in vivo and in cell culture. Toxicol Sci off Sci 98(1):75-86, PMID: 17283378, https://doi.org/10.1093/toxsci/kfm013.

49. Bimonte VM, Besharat ZM, Antonioni A, Cella V, Lenzi A, Ferretti E, et al. 2021. The endocrine disruptor cadmium: a new player in the pathophysiology of metabolic diseases. J Endocrinol Invest 44(7):1363-1377, PMID: 33501614, https://doi.org/10.1007/s40618-021-01502-x.

50. Aubé M, Larochelle C, Ayotte P. 2011. Differential effects of a complex organo-chlorine mixture on the proliferation of breast cancer cell lines. Environ Res 111(3):337-347, PMID: 21295777, https://doi.org/10.1016/j.envres.2011.01.010.

51. Kakarala M, Wicha MS. 2008. Implications of the cancer stem-cell hypothesis for breast cancer prevention and therapy. J Clin Oncol 26(17):2813-2820, PMID: 18539959, https://doi.org/10.1200/JCO.2008.16.3931.

52. McCarthy AM, Friebel-Klingner T, Ehsan S, He W, Welch M, Chen J, et al. 2021. Relationship of established risk factors with breast cancer subtypes. Cancer Med 10(18):6456-6467, PMID: 34464510, https://doi.org/10.1002/cam4.4158.

53. Carey LA, Perou CM, Livasy CA, Dressler LG, Cowan D, Conway K, et al. 2006. Race, breast cancer subtypes, and survival in the Carolina breast cancer study. JAMA 295(21):2492-2502, PMID: 16757721, https://doi.org/10.1001/jama.295.21.2492.

54. Danes JM, de Abreu ALP, Kerketta R, Huang Y, Palma FR, Gantner BN, et al. 2020. Inorganic arsenic promotes luminal to basal transition and metastasis of breast cancer. FASEB J 34(12):16034-16048, PMID: 33047385, https://doi.org/10. 1096/fj.202001192R.