Breast cancer (BCa) is the most common malignancy in women, and the incidence and subtype distribution vary widely among women of different racial/ethnic backgrounds [1–3]. Historically, the prevalence of BCa was higher in Caucasian (CA) women compared to African American (AA) women, while mortality rates were significantly higher in AA women [2,4–8]. Updated statistics show that BCa incidence rates among AA women converged with those of CA women in 2015 [9,10], and these numbers have remained dismal over the past 10 years [11–14]. In addition, clinical outcomes for AA women with BCa remain less favorable than those observed in CA and Asian women in the United States, and the incidence rates among younger women (under 50 years of age) are also increasing. For greater than 50 years, it has been well documented by researchers, clinicians, epidemiologists, and funding agencies, etc., that cancer disparities exist, particularly in Black people including AAs, who have some of the lowest survival rates compared to any other racial or ethnic group for most types of cancer [15,16]. Additionally, to a greater degree, AAs experience the outcomes of socially constructed health inequities rooted in structural racism [14]. Although it is well known, biospecimens from AA women have traditionally been difficult to obtain for biomarker studies, leaving a significant knowledge gap. Thus, there remains a need for a more comprehensive understanding of these disparities, increased investment in research, strengthened prevention programs such as screening for earlier detection, intervention, equitable health care delivery, and sustained application of established cancer control strategies, including smoking prevention and cessation programs, as well as initiatives promoting healthier diets and increased physical activity. Collectively, these efforts are widely recognized as critical to improving cancer outcomes and overall health [11–13,17].

The pathological features of BCa follow a progressive sequence from the transformation of a normal cell to benign proliferation, hyperplasia, atypical ductal hyperplasia (ADH), ductal carcinoma in-situ (DCIS) to invasive ductal (or lobular) carcinoma (IDC) and metastatic diseases [18,19]. It has been reported that the disease phenotype is distinguishable in ADH and progresses along distinct pathways for each of the primary clinical subtypes [20,21]. The genetic signature for disease heterogeneity across subtypes is greater than the heterogeneity of progression from DCIS to IDC within a subtype, suggesting that the disease subtypes have distinct progression pathways [21] or mechanisms. Current gold standard methods of detection and initial disease assessment include physical examination and periodic mammography imaging. However, both methods have low sensitivity and specificity for predicting the stage and grade of disease and can miss smaller masses or generate false positive or false negative results [22–25]. Especially early-on, when masses are at the most effective treatable stage, and this can depend on multiple factors, like heterogeneity in breast tissue architecture based on density [26–28]. Definitive BCa diagnoses will likely always require tissue biopsy, allowing for well-established genetic and molecular evaluation; but this method is invasive, time-consuming and can risk overtreatment in cases where masses are found to be benign/non-malignant. Thus, there is a need for non-invasive screening methods that capture molecular and mechanistic changes defined by the presence of a tumor, such as new biomarkers [29], which can also potentially distinguish between less aggressive and aggressive forms of BCa that will progress and have the potential to become metastatic and a higher propensity for recurrence.

The cancer phenotype can be identified at an early stage because the genes altered in this stage result in further post-transcriptional and -translational modifications as the disease progresses [30]. In addition, our recent studies [31–38] and those of others [39–41] have shown an increased level of common metabolites in BCa tissue. Metabolomics can be used to detect these modifications and identify the contributing and critical factors [42]. Metabolic profiling is of special importance in cancer biology due to profound changes in central metabolism associated with many tumors as established early by biochemical studies [5,43]. Metabolite concentrations in cellular systems are highly dependent on the physiological, environmental and genetic status of an organism and are regarded as the ultimate result of cellular regulation, resulting in the visible phenotypes [44]. Therefore, the comprehensive analysis of metabolite levels and fluxes provides a powerful tool for assessing the degree of perturbation in biological systems [32,38]. We are especially focused on aggressive disease in the AA population, who in the past 10 years have had higher BCa incidence and the highest mortality rates of triple-negative (TN) disease compared to all other populations [8,12–14].

The goals of this exploratory pilot study were to identify novel candidate biomarkers and metabolic mechanisms that could potentially be leveraged for earlier detection and improved classification of BCa progression. As a first step toward this goal, we used breast tissue specimen from AA women with and without BCa and identified distinct metabolic profile differences associated with disease status compared with non-malignant tissues, as well as metabolic alterations across histological grades I-III. The library-matched metabolites were evaluated for associations with established BCa risk factors, including age, body mass index (BMI), smoking status, menopausal status, hormone receptor status, and other demographic and clinicopathological characteristics.

Study cohort participants and description of biospecimens

For this study, 48 breast tissue samples (Table 1) were selected from our biorepository at Howard University Cancer Center (HUCC) Tissue and Plasma Bank. Biospecimens were collected from AA women residing in the Washington DC Metropolitan area between 2007 and 2013. Participants included those screened for suspicious BCa lesions, those diagnosed with BCa, and/or undergoing reduction mammoplasty surgery (serving as control tissues) at the Howard University Hospital and Providence Hospital. In collaboration with the surgical teams and Pathology Departments at each institution, excess tissue biospecimens were obtained from the Pathologist, following histopathological conformation of tumor-dense regions. Immediately following collection, fresh tissues were maintained under ultra-cold storage at -80°C. Specimens were either placed directly on dry ice in the operating room or from transferred from the Pathology Department shortly after surgical procedures. Thus, according to the study protocols, no significant degradation of specimen integrity occurred between collection and long-term storage, and no freeze-thaw cycles were introduced that could potentially affect metabolomics analyses. The histopathological diagnosis of all control biospecimens was independently reviewed and confirmed by a board-certified pathologist. Specimens were classified as normal mammoplasty tissue, fibrocystic breast tissue, or fibroadenoma based on established histopathological criteria before inclusion in the metabolomics analysis cohort. In December 2015, biospecimens were shipped to RTI International, and metabolomics data acquisition was performed in 2016 by the Eastern Regional Comprehensive Metabolomics Recourse Core (RTI ERCMRC).

The cohort study from which the biospecimens were collected was conducted in accordance with the Declaration of Helsinki. All participants provided written informed consent for the collection and research use of their biospecimens, including the tissue samples analyzed in this study. Collections were made under the protocols approved by the Institutional Review Boards (IRBs) of Howard University (IRB-12-CMED-39) and Providence Hospital (IRB-2014-0901). Biospecimens were collected along with corresponding demographic and clinical/pathological information. All study participants and samples were internally coded and subsequently re-coded by the ERCMRC conducting the metabolomics analyses to ensure patient confidentiality.

Metabolomics sample preparation and data acquisition

Metabolites were extracted from tissue samples and analyzed by broad spectrum (untargeted) metabolomics using proton (¹H) nuclear magnetic resonance (NMR) spectroscopy, as described previously [45–50]. Briefly, 48 frozen tissue samples (~100 mg per sample) were transferred to tubes containing homogenization beads. Sufficient extraction solvent (50:50 Acetonitrile:water) was added to each tissue sample to generate 1 mg/mL lysates using a MagNA Lyser (Roche Life Science, Pleasanton, CA, USA). Samples were centrifuged at 16,000 rpm at 4°C, and 500 µL of supernatant was transferred to clean tubes. An additional 125 µL aliquot from 23 randomly selected samples was pooled to generate analytical quality control (QC) samples. Samples were lyophilized to dryness and reconstituted in 700 µL of phosphate-buffered deuterium oxide (D₂O) Master mix containing the Chenomx internal standard (ISTD) [0.6 mM 4,4-dimethyl-4-silapentane-1-sulfonic acid (DSS-D6; Chemical Shift Indicator)] and 0.2 % NaN₃ with Imidazole. Tubes were vortexed, centrifuged, and a 600 µL aliquot of each sample was transferred to 5 mm NMR tubes (Bruker-BioSpin, Switzerland) for data acquisition.

NMR spectra were acquired for individual tissue samples and the pooled QC samples using a Bruker Avance III 700 MHz NMR spectrometer (Bruker-BioSpin, Rheinstetten, Germany), equipped with a cryogenically cooled 5 mm Automatic Tuning and Matching Accessory (ATMA) probe. Data were collected at 25°C and a 1D NOESY pre-saturation pulse sequence (noesypr1d) [45]. ¹H NMR spectra were preprocessed using ACD NMR Processor 12.0 software (Advanced Chemistry Development, Toronto, ON, Canada). All spectra were binned using intelligent bucketing with width of 0.04 ppm and a 50% looseness factor. Spectral regions from 0.70–9.00 ppm were included, while regions corresponding to water (4.77–4.90 ppm) and methanol (3.33–3.36 ppm) were excluded. Each bin integral was normalized to the total spectral integral. NMR signals in the spectra were matched to metabolites [32] (annotation or identification) using the Chenomx NMR Suite 8.1 Professional software (Chenomx, Edmonton, Alberta, Canada). The Chenomx software contains a pH sensitive internal metabolites library of standards run on a range of NMR frequencies and peak clusters modeled to chemical shift and J-coupling) [51]. In addition, annotations of lipid peaks were based on literature [52]. It should be noted that metabolites in NMR spectra may contain one or more peaks (chemical shift) and peak splitting (J-coupling) based on the chemical structure and chemical environment of the metabolite.

Multivariate analyses

Multivariate analyses, including principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA), were performed using the normalized binned data which were first Pareto-scaled and mean-centered. Analyses were performed using SIMCA 14.1 (Sartorius Stedim AB, Umeå, Sweden). For this exploratory pilot study, multivariate analysis as an a priori data analysis plan, was used to better understand the differences in phenotypic groups. The PCA scores plots were inspected to confirm that the phenotypic pooled samples clustered tightly within the center of their phenotypic groups and that the total study pool clustered centrally among all samples (Supplementary Figure S1). This QC approach is widely applied in metabolomic studies [53]. Loadings plots and variable influence on projection (VIP) plots were inspected and data lists were exported into Excel. The VIP statistic summarizes the contribution of each bin to group separation (e.g., BCa cases versus controls) in multivariate space [53]. All models were evaluated using 7-fold cross-validation to assess the predictive variation of the model (Q²).

Statistical analyses were performed using SAS 9.4 (SAS Institute Inc., Cary, NC) to calculate fold changes and nominal p-value for bins associated with differences in the VIP scores. For initial case-control comparisons of continuous variables, either t-test assuming unequal variances or the Wilcoxon Rank Sum Test were used because sample sizes were small. As this was an exploratory analysis, p-values were not adjusted for multiple testing [54], but the FDR p-values were calculated and reported when possible. Bins were considered important for distinguishing study groups if they met at least one of the following criteria: 1) VIP≥1.0 with a jack-knife confidence interval not including 0; 2) fold change magnitude ≥2; or 3) p-value <0.10 [55]. All peaks identified in bins that were important for differentiating the study groups were assigned to metabolites using Chenomx NMR Suite 8.1 Professional software (Edmonton, Alberta, Canada), which contained a reference library of approximately 350 compounds at the time of analysis. Subsequently, all 212 bins were evaluated using a consistent rank-based framework and FDR correction, as further described in Section 2.5. Positive signed fold change indicates a higher median in Group 1 than Group 2; where negative values indicate a lower median in Group 1. All raw and processed analytical data, along with associated de-identified metadata, have been deposited in the publicly accessible NIH Common Fund Metabolomics Data Repository:

http://www.metabolomicsworkbench.org/data/DRCCMetadata.php?Mode=Study& StudyID=ST000438&StudyType=NMR&ResultType=2.

Pathway analysis

The publicly available MetaboAnalyst 4.0 pathway enrichment platform (Xia Lab, McGill University and Wishart Research Group, University of Alberta, Canada) was used to identify perturbed canonical pathways and to assess enrichment of endogenous metabolic network differences between BCa cases and controls [56,57]. The threshold and p-value (p≤0.05) criteria were included based on the calculated fold-changes (threshold) and p-values (p≤0.10) derived from the binned data for metabolite-based analyses, unless otherwise specified in the table legend.

Statistical analysis

More robust analyses were performed using SPSS to assess any significant correlation between library-matched spectral bins, through multivariate analyses and library-matched metabolite annotations that distinguished breast tissue phenotypes in case versus control groups, with a range of clinical and pathological BCa risk factors. Associations were first evaluated individually with respect to age range (≤54 years old vs. ≥55 years), Body mass index (BMI) (<30 vs. ≥30), tumor Grade (I, II or III), tumor Stage (TNM stages I-IV), lymph node status, estrogen receptor (ER) status (positive + or negative -), progesterone (PR) status (+ or -), HER2 status (+ or -), triple negative breast cancer (TNBC) status (yes or no), family history of cancer (yes or no), smoking status (yes, no, or former), alcohol use (yes or no), and menopausal status (yes or no). Minimum detectable standardized mean differences at 80% power and two-sided alpha=0.05 were calculated to contextualize the limited precision of small subgroup analyses.

Following the univariate analysis, all 212 spectral bins were reanalyzed using two-sided Mann-Whitney U tests for two-group comparisons and Kruskal-Wallis tests for variables with more than two levels. A uniform rank-based framework was used because Shapiro-Wilk testing showed that 47%-72% of bins were non-normal in at least one group across the principal comparisons and because several subgroups were small. Benjamini-Hochberg FDR (BH-FDR) correction was applied separately across the 212 bins within each prespecified tissue comparison and each clinical/pathological variable. PCA/QC assessment repeated cross-validated PLS-DA sensitivity analyses, permutation testing, outlier sensitivity were also conducted. Both nominal p-values and adjusted q-values were evaluated and q<0.05 was considered FDR-supported. Rank-biserial correlations were calculated for two-group comparisons.

Cohort characteristics

Biospecimens utilized for this pilot study were obtained from a cohort of AA women, who provided informed consent to participate and were seen at the Howard University Hospital, the Howard University Cancer Center, or Providence Hospital in Washington DC. Tissue samples were previously collected from women aged 21–78 years and categorized into two age groups (≤54 years and ≥55 years) for analysis. Samples were also classified as one of three types for breast tissue controls with no BCa, as reduction mammoplasty-normal, mammoplasty-with fibrocystic changes or mammoplasty-with fibroadenoma and compared to BCa cases, which included clinically localized disease (atypical ductal hyperplasia [ADH] or ductal carcinoma in situ [DCIS]), or invasive BCa (invasive ductal carcinoma, infiltrating duct carcinoma, inflammatory carcinoma, and pleomorphic carcinoma). Cohort characteristics are detailed in Table 1, along with aggregated corresponding demographic, clinical, and pathological information. No significant differences were observed in the mean age at diagnosis between participants in the fibrocystic tissue group and those in the BCa tissue group (p=0.2628). Similarly, BMI didn’t differ significantly among the groups (p=0.5944). The breakdown of the numbers of samples in each phenotypic tissue group is detailed in Supplementary Table 1.

Metabolomic profiling distinguishes normal breast tissues from benign diseases

All tissue biospecimens were prepared as described in the Materials and Methods section for metabolomics data acquisition using ¹H-NMR spectroscopy. Following data acquisition, the raw spectral free induction decays (FIDs) were zero-filled to 64k, and an exponential line broadening factor of 0.5 Hz was applied prior to Fourier transformation. The spectra were then manually phased, baseline-corrected, and referenced to DSS as a part of the QC process. Subsequently, spectra were aligned and binned for further analyses. To determine tissue-level metabolomic differences among the three breast tissue control groups, we first performed pairwise multivariate analyses in SIMCA comparing mammoplasty-normal breast tissues to mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma. These groups represented benign and non-malignant breast lesions obtained from women undergoing reduction mammoplasty or breast mass biopsy procedures. Although all three tissue types were initially selected as controls to characterize the metabolic alterations expected to differentiate BCa tumor samples, the orthogonal partial least squares discriminant analysis (OPLS-DA) plot unexpectedly revealed differentiation (based on VIP≥1.0 with a jack-knife confidence interval not including 0; or 2) fold change magnitude ≥2; or 3) p-value <0.10) between the three mammoplasty sample groups (black dots vs orange dots), with the latter two clustering together as a single group in the preliminary multivariate analysis (Figure 1A). The mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma tissue groups were analyzed together due to sample numbers and pathobiological similarity compared to fully normal mammoplasty or BCa biospecimens. Even by unsupervised principal component analysis (PCA), differentiation can be seen (Supplementary Figure S1). Although several studies use mammoplasty-normal tissue alone, it is more common to see a mixed combination of all three being used as controls for comparison with BCa. Our results reflect that benign fibrocystic or fibroadenoma tissues differ from mammoplasty-normal tissues not only by histopathological features (i.e., underlying proliferative or structural alterations, microenvironmental inflammation, etc.) but as demonstrated here, also by metabolic profile. Though preliminary, because of unadjusted statistical analyses and small sample sizes, these metabolic differences identified 43 library-matched metabolites (Table 2) that may reflect group separation driven by morphological differences, as we also see distinction in corresponding pathway enrichment analysis (Figure 1B) and for perturbed endogenous metabolomic pathways (Figure 1C). The top five significantly enriched pathways (p<0.05), although for unadjusted measurements, included Ammonia Recycling, Glycine and Serine Metabolism, Aspartate Metabolism, the Urea Cycle and Arginine and Proline Metabolism.

Figure 1. Comparisons of control mammoplasty-normal versus mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma breast tissues. A. OPLS-DA multivariate plot demonstrating metabolic separation between control tissue samples from women who underwent reduction mammoplasty surgery versus women who underwent surgical removal of benign fibrocystic or non-malignant fibroadenoma breast lesions. Multivariate model statistics are R²X(cum) =0.599, R²Y(cum) =0.602, Q2 =0.629. B. Overview of metabolic enrichment analysis. C. Overview of significantly different endogenous pathways (p<0.05, unadjusted), where both B and C were analyzed by MetaboAnalyst. Corresponding PCA plot can be reviewed in Figure S2.

Metabolomic profiling distinguishes all controls from BCa biospecimens

Among the tissue cohort sample comparisons, 50 putatively identified, library-matched metabolites were identified as important for differentiating all biospecimens categorized as controls (mammoplasty-normal vs mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma) from women with BCa across all stages and grades analyzed (Table 3). These metabolites correspond to the OPLS-DA plot in Figure 2A, in which all breast tissue controls are represented by black dots and all BCa tissues are represented by pink dots. In this preliminary comparison, nine of the 25 top enriched biological pathways (Figure 2B) aligned with the same significantly (p<0.50, unadjusted) differential perturbed endogenous pathway perturbations (Figure 2C). These included alterations in pathways involved in Carnitine Synthesis, Phosphatidylethanolamine Biosynthesis, Methionine Metabolism, Ammonia Recycling, Phospholipid Biosynthesis, Arginine and Proline Metabolism, Betaine Metabolism, Phosphatidylcholine Biosynthesis and Glycine and Serine Metabolism.

Figure 2. Comparisons of all controls combined (mammoplasty-normal vs mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma) versus BCa tissues. A. OPLS-DA multivariate plots demonstrating metabolic separation between all tissue controls compared to tissues from women with BCa (stages I–IV and grades I–III), with multivariate model statistics of R²X(cum) =0.743, R²Y(cum) =0.605, Q2 =0.351. B. Overview of metabolic enrichment analysis. C. Summary of significantly altered endogenous pathways (p<0.05, unadjusted). Analyses presented in B and C were performed using MetaboAnalyst. Corresponding PCA plot containing all study samples can be reviewed in Figure S1, also containing study Pool QC samples denoted by grey dots and red arrows demonstrates low analytical variation.

Metabolites that distinguish histological groups of BCa correlate with BCa risk factors

Having the opportunity to work with biospecimens from an understudied population using metabolomics, we were interested in several other comparative differences and new hypotheses the analytical method might reveal in the exploratory dataset. Since the OPLS-DA scores plot showed a differentiation between the mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma samples and mammoplasty-normal BCa samples (Figure 1 and Table 2), we decided to use only the mammoplasty-normal samples in all subsequent analyses as the reference control. The level of biological heterogeneity among control group samples might have diminished detectable differences between comparisons with smaller sample size numbers in the subgroups, that we present in the remaining results.

In the next multivariate comparison between mammoplasty-normal control tissues versus grade (G) I-II BCa samples, the OPLS-DA plot in Figure 3A demonstrates clear separation for nearly all control samples from BCa cases. Interestingly, one control tissue clusters more closely with the cancer cases along the primary component (t[1]) of the OPLS-DA model, suggesting the possibility of a more subtle, progression-related metabolic shift and warrants further investigation. The importance of such future analysis will be discussed in the next section. Figure 3B shows less pronounced separation between control mammoplasty-normal tissues and grade three (GIII) BCa samples. In both panels of Figure 3, additional stratification is provided by coloring the GI–II BCa samples according to estrogen receptor (ER) status, where GI-II, ER+ samples are indicated by blue dots and GI-II, ER- samples are indicated by pink dots (Figure 3A). Likewise, GIII, ER+ samples are indicated by turquoise dots, while GIII, ER- samples are indicated by dark green spots (Figure 3B). In both pairwise comparisons, the mammoplasty-normal control samples are presented as black dots. Distinct metabolic signatures were again identified by multivariate analysis and best visualized in the Venn diagram (Figure 3C), where three metabolites separate early-grade BCa (GI–II), and ten metabolites uniquely separate later-grade disease (GIII) relative to the normal reference control. Supplementary Table S3 provides the full list of all cross-sectional grade-related metabolite differences. Metabolites that are unique to specific comparisons are highlighted in bold and shown in blue font in the table.

Figure 3. OPLS-DA plot comparisons of control mammoplasty-normal versus BCa tissues across clinical grade. A. Mammoplasty-normal versus GI-II BCa tissues in a pair-wise comparison, where the BCa tissues are different colors just to show how they naturally separate in multivariate space by estrogen receptor (ER) status. Multivariate model statistics are R²X(cum) = 0.604, R²Y(cum) = 0.611, Q2 = 0.041. B. Mammoplasty-normal versus GIII BCa tissues, again where BCa tissues represent one comparison group, but are colored differently to show whether the samples separate by ER status. Multivariate model statistics are R²X(cum) = 0.711, R²Y(cum) =0.699, Q2 =-0.477. Corresponding PCA plots can be viewed in Figure S3 for A. and Figure S4 for B. C. Venn diagram of the common and unique metabolites between these individual analyses. Three were GI-II–specific and ten were GIII–specific.

Several unique metabolites also distinguished ER+ from ER- tissues in both GI-II (Figure 4A) and GIII (Figure 4B) comparisons, where ER status (positivity or negativity) was used for pairwise stratification. Among the metabolites common to both comparisons, five exhibited opposite fold-change directions or magnitude when comparing GI-II to GIII. For example, the graph in Figure 4C shows that myo-inositol, lactate, and proline are elevated in ER- GI-II BCa samples and elevated at the later stage/grade GIII in ER+ tissues. In contrast, arginine and valine display a reciprocal pattern, being elevated in ER+ GI-II tissues, and elevated in ER- GIII BCa tissues. Thus, in ER+ tissues, myo-inositol, lactate, and proline are lower at the earlier GI-II stage of disease and increase in GIII tissues. Likewise, arginine and valine decreased from GI-II to GIII in ER+ BCa tissues. Conversely, for ER- tissues, myo-inositol, lactate, and proline decreased in comparisons between GI-II to GIII, while arginine and valine increased with advancing grade.

Figure 4. Comparison of BCa tissues based on estrogen receptor (ER) status. A. OPLS-DA plot of GI-II (ER^- versus ER⁺). Multivariate model statistics are R²X(cum) = 0.503, R²Y(cum) = 0.565, Q2 = 0.255. The corresponding PCA plot can be viewed in Figure S5. B. OPLS-DA plot of GIII (ER^- versus ER⁺) BCa tissues. Multivariate model statistics are R²X(cum) = 0.715, R²Y(cum) = 0.756, Q2 = 0.144. The corresponding PCA plot can be viewed in Figure S6. C. The graph illustrates differences among the most altered metabolites based on ER positivity, which have the potential to lead to the identification of subtype- or grade-specific biomarkers in AA women with BCa.

Next, to assess whether BMI was a key differentiating factor in this small cohort, we performed multivariate analysis in SIMCA across all tissue samples and observed a distinct separation pattern in the OPLS-DA plot. However, the model statistics were not strong (R2X(cum)=0.771, R2Y(cum) =0.542, Q2=-0.0889), and the negative Q2 indicated that the model was unlikely to be reproducible in a different dataset so we did not conduct any additional analyses related to BMI. Finally, to determine whether age served as a significant metabolic differentiator independent of BCa in AA women, we defined age thresholds at ≤54 years and 55 years to approximate pre- to peri-menopausal states compared to the post-menopausal state. In a likewise manner, the OPLS-DA multivariate model statistics (R2X(cum)=0.551, R2Y(cum)=0.291, and Q2=-0.45) also had a negative Q2 value and we report no further findings for either multivariate analysis. Although the plots demonstrated separation, as seen in Supplementary Figure S7, we believe that the limitation of small sample numbers overall and an imbalance across some groups comparisons (Supplementary Table 1) are in-part responsible for the poor statistical results in these comparisons, which have been shown in larger datasets to render statistically sound metabolomics data.

More robust secondary analyses using SPSS were performed to evaluate all NMR bins (Supplementary Table 2) in previous comparisons by multivariate analysis and for grade, ER status, BMI, age, and other clinical/pathological variables. Exact sample sizes were: mammoplasty-normal versus Grade I-II, 8 versus 12; mammoplasty-normal versus Grade III, 8 versus 18; Grade I-II ER-positive versus ER-negative, 8 versus 4; Grade III ER-positive versus ER-negative, 5 versus 13; BMI <30 versus ≥30, 18 versus 18; and age ≤54 versus≥55, 28 versus 20. Only two bins (b_89 and b_172) survived BH-FDR correction (p=0.000116 and 0.000202; q=0.0214 for each). Ten bins survived BH-FDR correction in the mammoplasty-normal versus Grade III comparison (minimum q=0.0385). No bin survived FDR in the Grade I-II, ER-status, BMI, or age comparisons. None of the 13 clinical/pathological variables had an association that survived BH-FDR correction. Subgroup findings, reported by multivariate analysis are therefore preliminary and should not be interpreted as validated progression, subtype, or risk-factor signatures. They do, however, provide justification to conduct a follow-up study in a larger cohort, to evaluate the subgroup pairwise comparisons and full study analysis with sufficient power. Currently, this study is underway, with planned enrollment of 300 cases and 300 controls.

The primary finding from this investigation is our metabolomics analysis of control tissue biospecimens which represents the first report demonstrating metabolic differentiation among mammoplasty-normal compared to mammoplasty-with fibrocystic changes and mammoplasty-with fibroadenoma, by 1H-NMR spectroscopy in a cohort of AA women (Figure 1, Table 2). These observations, while preliminary, may point to biologically relevant metabolic changes that have the potential to predispose these benign/non-malignant lesions toward malignant changes. However, we are cautious with this new observation, because of the exploratory nature of this pilot study (n=8 versus n=9), and any such finding requires confirmation and validation in both larger and more diverse cohorts. Further investigation regarding the role differentiation of control reduction mammoplasty tissues has for interpreting data where samples with non-malignant changes are used as a reference for BCa remains to be seen.

A secondary observation, also based on multivariate analysis, is from the pairwise comparison between control mammoplasty-normal tissues versus GI-II breast tumor samples, where differentiation included the metabolite 2-phenylpropionate. Though our exploratory study was not adequately powered to evaluate significance, the finding is interesting and could be truly novel in BCa biology. Prior reports have only shown this metabolite as relevant in liver or urinary excreta involved in studies of Chinese natural product hepatotoxicity or long-term use of steroid medications (Figure 3, Table 3). Thus, the presence of 2-phenylpropionate being identified as important in distinguishing cases versus mammoplasty-normal controls also requires further validation studies to fully understand the role this metabolite may play in BC etiology or disease progression.

Beginning in the early 1990s, a significant shift occurred in the understanding of the relationship between benign breast conditions, such as fibrocystic disease and fibroadenoma disease, and subsequent risk of BCa. Earlier views, dating from the 1920s to the early 40s, suggested no association between these benign conditions and BCa risk compared with women without such diseases. This perspective evolved to recognize a modestly increased likelihood of carcinoma [58]; and later to evidence indicating a two- to fourfold elevated risk [59], particularly among women with a familial history of BCa in the preceding four decades. Notably, this shift in understanding was driven primarily by findings from pathological studies. Thus, conditions once considered entirely benign and unrelated to BCa development, were re-evaluated after accumulating experimental and epidemiological evidence alongside ongoing refinement in pathological classification, challenged this dogma [60–74]. These advances led to a more nuanced understanding, indicating no increased risk among women without a family history and/or in the absence of specific morphology-based pathological features.

In one of our previous investigations, we identified a significant epigenetic alteration characterized by elevated DNA methylation level in the P1 promoter region of ESRα (p=0.0001) in non-malignant fibroadenoma compared with ER-negative BCa tumors. In the same cohort of AA women, on which the metabolomics data presented herein are based, we also observed a two-fold increase in ESRα expression in both benign tissue types (fibrocystic and fibroadenoma) [75]. Even so, the literature remains inconsistent on this topic. It is therefore important to emphasize that the most common forms of fibrocystic breast changes typically observed in women between the ages of 30 and 50 years old, are not associated with an increased risk of developing BCa compared with women in the general population without such benign lesions, particularly in the absence of a family history of BCa or histological evidence of proliferative hyperplasia [76,77]. However, when uncertainty remains, the most appropriate course of action is to pursue diagnostic clarification, often through biopsy.

Only two prior studies reported findings using high-resolution liquid-chromatography coupled to mass-spectrometry (LC-MS). Hines et al, distinguished fibroadenoma from BCa in a total of six biospecimens [78], while Liu et al. more recently applied MS-based metabolomics to differentiate four TNBC samples from four fibroadenoma samples, in conjugation with microbiome profiles [79]. In a third study, Raman spectroscopy was used to define a biochemical fingerprint distinguishing breast phyllodes tumors from fibroadenomas [80]. Consistent with findings from the tissue- based MS studies, our data similarly suggest that glutathione level was higher in mammoplasty-normal tissue compared with the mammoplasty-with fibrocystic changes and mammoplasty-w fibroadenoma (Table 2).

In addition, in this pilot study we observed increased levels of glutathione, glutamine, carnitine, creatine, creatinine, choline, O-phosphocholine, O-acetylcholine, and sn-glycero-3-phosphocholine (Table 3), consistent with previous published findings from MS-studies in both tissue and serum comparing control and BCa samples. Our preliminary multivariate results demonstrate metabolic distinction between histologically normal breast tissues obtained from healthy women who underwent reduction mammoplasty surgery (n=8) and those from healthy women who had undergone surgery for benign fibrocystic disease or non-malignant fibroadenoma removal (n=9) (Figure 1A, Table 2), together with corresponding impacted endogenous biological pathways.

Although non-malignant, these lesions may be classified as non-proliferative or proliferative in nature and can exhibit obstructive features and morphological elements of noninvasive hyperplasia [59]. In a manner analogous to BCa, which is clinically characterized by six to seven genetically, molecularly, morphologically, and pathologically distinct subtypes, fibrocystic changes and fibroadenoma are also sub-divided into a spectrum of morphological and pathological variants. Accordingly, future technological and analytical advances may further refine our understanding of benign breast disease, potentially revealing a more complex elevated relationship with BCa etiology than is currently appreciated.

In our next analysis comparing phenotypic groups across all tissue types, (Figures 1A, 2A, 3A, and 3B), again distinct metabolic signatures emerged spanning mammoplasty-normal controls grouped with mammoplasty-w fibrocystic change and -w fibroadenoma controls, compared to BCa stratified by grade and stage (Supplementary Table S3) in discreet pairwise comparisons. These metabolites are most clearly illustrated in the Venn diagram (Figure 3C), which highlights several metabolites unique to each preliminary group comparison. For example, five metabolites were uniquely library-matched in Chenomx for the comparison differentiation between mammoplasty-normal controls and mammoplasty-w fibrocystic change and -w fibroadenoma, including 4-hydroxybenzoate, dimethylamine, glutathione, methionine, and UDP-N-acetylglucosamine. Because it is a pilot, it is fundamental to contextualize these findings within the literature, which substantiates that glutathione, a key determinant of cellular redox balance and detoxification capacity, may play a central role in regulating or reflecting the metabolic alterations observed in the other four metabolites, given its critical function in maintaining cellular homeostasis. It is also important to note that both 4-hydroxybenzoate and UDP-N-acetylglucosamine (UDP-GlcNAc) were likely elevated in mammoplasty-w fibrocystic changes and mammoplasty-w fibroadenoma tissues by approximately 2-to 3-fold (Table 2), whereas dimethylamine, glutathione and methionine were likely higher in the mammoplasty-normal tissues, with increases ranging from 1.6- to 2.4-fold compared with mammoplasty-w fibrocystic changes and mammoplasty-w fibroadenoma lesions.

Again, because our results are based on multivariate analysis findings in this small exploratory cohort, we did a deeper review of the literature to seek possible interpretation of the observed, putatively annotated spectral regions that matched to metabolites (Chenomx library) in the comparison between fully histologically normal breast tissue and those with benign breast diseases and found that 4-hydroxybenzoic acid esters belong to the class of compounds known as parabens. These compounds were established more than two decades ago to exhibit estrogenic activity and influence cell proliferation in BCa cellular models [81]. Despite improvements in environmental regulation of parabens to limit our vast exposure to these compounds because of their detrimental impacts on human health (i.e., genotoxicity and promoting cellular proliferation) [82], they are still used to some degree as a preservative in food products to increase shelf-life by blocking microbial growth, in cosmetics and pharmaceuticals, where benzyl-paraben compounds render some of the greatest toxic effects [83]. Downstream toxic metabolites from parabens and similar environmental exposures to compounds like polycyclic aromatic hydrocarbons (PAHs), which have also been linked to malignant initiation across several tissue types, all rely on glutathione and its associated transferases or GSTs for detoxification by the formation of glutathione conjugates for removal [83,84].

In a similar manner, the observed upregulation of UDP-GlcNAc in benign breast disease tissues, relative to normal reduction mammoplasty tissue, aligns with established evidence linking increased glucose flux and accumulation of glycolytic intermediates (i.e., UDP-sugars) to the metabolic reprogramming that supports high-energy demand states in malignant tissues, commonly described as the Warburg Effect [85]. Beyond glycolysis, UDP-GlcNAc may also contribute to increased synthesis of tumor-supportive extracellular matrix components, such as hyaluronan [86], thereby promoting a microenvironment that can facilitate pro-tumorigenic processes. Although not exclusive to any single comparison, both glutamine and glutamate emerged as important metabolites across all multivariate analyses and may be pointing to functioning in concert with glucose metabolism through glutaminyls and the hexosamine biosynthetic pathway. This activity has been shown to support synthesis through UDP-GlcNAc [87–90] as an alternative energy generating and biosynthetic route.

Finally, we propose and will test in future investigations, that the elevated glutathione observed in mammoplasty-normal breast tissue, identified as important in the discreet pairwise multivariate comparisons, could be reflecting a more robust detoxification capacity supporting cellular homeostasis. In contrast, this system comparatively diminished in mammoplasty-w fibrocystic changes and mammoplasty-w fibroadenoma tissues, which may suggest a shift away from the fully normal metabolic state observed in reduction mammoplasty-normal samples (Table 2). Though preliminary and needing further validation, this metabolic alteration was further supported by pathways analyses, where it emerged as a significantly perturbed endogenous pathway (unadjusted) in the non-discreet comparison grouping all control tissues versus all BCa tissues (Figures 2B and 2C).

In the 1990s, glutathione along with glutathione S-transferases (GSTs) and associated cytochrome P450 enzymes, acting in concert as a major antioxidant detoxification system, particularly in the liver, came to be understood in an increasingly complex light. Initially, considered as potential important determinants of breast tissue susceptibility to carcinogenic compounds (including chemotherapeutics) [91], they were subsequently linked to the regulation of key cancer hallmark, including the evasion of apoptosis and promotion of cellular proliferation [92]. Ultimately, these systems were recognized as major contributors to tumor-driven multi-drug resistance against chemotherapeutic agents and other xenobiotics [93]. From the 1990s to the present, a substantial body of evidence has accumulated detailing the roles of glutathione and its associated superfamily of enzymes across multiple malignancies, notably in BCa [32,94–100]. Collectively, this work nominally (until validation) underscores their broad involvement in virtually all major mechanistic hallmarks of cancer, further highlighting their central role in tumor biology and progression [85].

In analyses comparing all three types of control tissues with BCa tumor tissues (Figure 2A), a subset of control samples (black dots) clusters more closely with, and in some cases overlap, the BCa samples (pink dots) along the primary t[1] axis in the OPLS-DA multivariate plot. Beyond validation of this preliminary result in a larger cohort, future analyses will focus on delineating the metabolic differences among these more discreet sample subsets, with the possibility of identifying progression-specific metabolic phenotypes. Such insights may enable the development of more selective biomarkers that capture the transition from normalcy to benign disease and ultimately to malignancy. In addition, identification of these types of biomarkers could help refine early disease detection, complement imaging modalities, and support more precise characterization of disease onset prior to definitive diagnosis, staging, and grading through tissue biopsy. Additionally, three metabolites uniquely distinguished mammoplasty-normal from GI-II BCa tissues (Figure 3A), include the organic dicarboxylic acids fumarate and succinate. These metabolites are key intermediates of the citric acid cycle [101] and potentially point to the elevated energetic demands characteristic of tumor tissue relative to normal breast tissue. Notably, succinate has also been more recently implicated in anti-tumor immunity, where it has been shown to modulate mitochondrial function and influence CD8⁺ T-cell activity [102,103].

The third metabolite uniquely identified by library-matching in this group comparison was 2-phenylpropionate, for which no prior reports were found linking it to any type of breast tissue, either mammoplasty-normal or BCa tumors. However, two metabolomics-based studies have reported it in other biological contexts; where disruption of 2-phenylpropionate homeostasis was among four key metabolic perturbations associated with dose-dependent hepatotoxicity induced by the traditional Chinese medicine genipin [104]; and when examining the long-term effects of the glucocorticoid prednisone on gut mycobiome dysbiosis in a rat model [105], a positive correlation between levels of 2-phenylpropionate associated with fungal species. These reports are grounded in the observation that 2-phenylpropionate represents a key structural moiety within the compounds under investigation, such as genipin, or that it rises as an altered microbial metabolite associated with the metabolism of various fungal and bacterial species [105–108,]. This observation also supports our ongoing investigations into potential interactions between the gut microbiome and BCa etiology, in which dysbiosis may facilitate the translocation of microbe populations from the typically indolent gut microenvironment. Within this framework, such microbial shifts may plausibly acquire pathogenic potential and contribute to disease initiation and early progression [109–111].

In Figure 3B, mammoplasty-normal tissues were compared with GIII BCa tissues, yielding a more extensive list of nine unique metabolites that may reflect the hallmark of increased metabolic dysregulation with advancing disease. We want to highlight hypoxanthine and uridine for further discussion. The purine metabolite hypoxanthine was preliminarily observed to be elevated in GIII BCa samples, showing a 1.7-fold increase relative to mammoplasty-normal controls. In a case-control cohort of women from the United Arab Emirates, hypoxanthine was also identified as a significantly differential metabolite using MS-based metabolomics, and associated with disruptions in endocrine and immune-related metabolic pathways [112]. The same research group further reported its direct link in key hallmark mechanisms underlying BCa metastasis (i.e., upregulation of lipid metabolism in an ER+ cell line model) [113]. With respect to uridine in BCa progression [114], it has been associated previously with stemness in more advanced and aggressive TNBC cell lines [115], and has also been implicated in the promotion of a stem-like phenotype through stress-induced gut dysbiosis via LRP5/β-catenin signaling pathway [116].

Similar to findings from other metabolomics studies of BCa cell models, plasma, and tissue biospecimens, both the mammoplasty-normal versus GI-II (Figure 3A) and mammoplasty-normal versus GIII (Figure 3B) comparisons could be pointing to considerations that, although not exclusive to these groupings, low density lipoprotein (LDL) and very low-density lipoprotein (VLDL) levels were reduced in the BCa samples relative to controls. In the literature, lipoprotein changes have been suggested to be potentially protective against BCa risk [117] and reflect alterations in cancer related immunity processes, particularly in more aggressive GIII TNBC samples [118]. Collectively, these metabolites may therefore complement uniquely altered features as part of a broader biomarker panel for monitoring disease progression. Our preliminary results may help better understand what that biomarker panel will ultimately look like, after validation in larger and more diverse cohort investigations specifically designed to evaluate metabolic progression.

Pathway analysis across the comparison groups also revealed distinct differences, with Ammonia Recycling emerging as a significantly perturbed pathway (p= 0.02–0.008) in both the mammoplasty-normal versus mammoplasty-w fibrocystic change and mammoplasty-w fibroadenoma and the all controls grouped versus all BCa tissue comparisons. This pathway has previously been shown to facilitate tumor adaptation by enabling BCa cells and tissue to detoxify accumulated ammonia within the tumor microenvironment, while recycling it for amino acids synthesis. In this way, it supports enhanced cellular proliferation and increased biomass accumulation during tumor growth [119,120]. While preliminary and needing to be confirmed in larger cohort studies and with more robust statistical methods, the observed pathway changes typically reflect growth-supporting biosynthetic programs that are consistent with the increased biomass demands of proliferating cells. Notably, the metabolic process underlying the Warburg Effect, a fundamental hallmark of cancer metabolism, was also enriched in both comparisons, even though only the second analysis included BCa tumor tissues. This observation is interesting and clearly requires further investigation.

We acknowledge that as an exploratory pilot, our study has five major limitations. Firstly, although intentionally designed to focus on women who self-identified as AA, ancestry was not further characterized or genetically confirmed using ancestry-informative markers. While this may be considered a limitation for certain statistical analyses, we contextualized our findings by comparing them with studies of similar objectives that employed the same metabolomics platform (¹H NMR), in order to assess the generalizability of our results and determine whether observed metabolic patterns were specific to this cohort of women with disease. Secondly, participants were recruited primarily from the Washington, DC region in close proximity to Howard University, a community-based hospital serving a lower-resourced population within the metropolitan area. To help mitigate potential selection bias, we included participants recruited from a second clinical site, Providence Hospital, located approximately three miles northeast of Howard University, and about 15 minutes near the Washington DC-Maryland border. This additional recruitment site represented a distinct patient population with a different socioeconomic demographic of participants, thereby broadening the diversity of the cohort. Thirdly, our study was supported by pilot program funding, and the resulting sample size of tissue biospecimens, especially across subgroups, was limited, which constrains statistical power. In addition to the overall size of the cohort evaluated, there were also imbalances in the subgroups analyzed, which likely contributed to many metabolites falling out of significance (by unadjusted multivariate analyses) and resulting in no significant associations between specific metabolites and clinical/pathological parameters after BH-FDR correction was applied in secondary statistical analyses. Although larger cohorts will be necessary for full reproducibility and independent validation, metabolomics has been widely recognized as a highly sensitive analytical approach capable of detecting robust biochemical differences between distinct phenotypes, even in studies with relatively small sample size (e.g., n=5 per group) when evaluated using multivariate analytical framework [121–123]. Thus, we are confident in progressing toward validation of these findings in a larger cohort study and are currently conducting an observational trial at the Howard University Cancer Center. This study will include concordant biospecimens of plasma, urine, and tissue from an expanded cohort (n=300 cases and 300 controls). Fourthly, some of the corresponding metadata provided by the Tumor Registry were incomplete for certain variables used in statistical analyses. As a result, the total number of samples included in each model varied depending on the variability of complete metadata across comparative groups. Smoking status and alcohol use were the categories with the highest proportion of missing data. This is reflected in Table 1, where the totals for clinical and pathological variables do not always add up to the overall number of tissue samples (n=48). Lastly, results of this discovery study were based on binned NMR data. Although Chenomx software contains a pH sensitive experimental library of standards for library matching, the binning approach can result in more than one bin per metabolite or more than one metabolite in a bin. In addition, lipid annotations were based on the literature. Additional 2D and/ or multidimensional NMR experiments are needed to further confirm metabolite identity with confidence. In addition, findings of this study need to be confirmed by using a targeted metabolomics study using metabolite standards and calibration curves. Metabolite pathways found to be perturbed in this study can be interrogated by using metabolites in those pathways using a similar approach.

Our exploratory pilot study was designed to identify metabolites that distinguish between cases and controls and histological differences of BCa across stage and grade of disease. We identified several promising candidates that will serve as a starting point to identify associated proteins and enzymes related to metabolism that results in these metabolites after validation. Further, their identification allows for hypothesis development to discover potential targets that may serve as biomarkers for earlier detection of disease or to monitor disease progression in an understudied population, based on our exploratory cross-sectional grade findings. Our current results and ongoing research establish a path to achieving these goals. As presented findings derive from a pilot study, a targeted metabolomics approach could be used to confirm annotations, while further investigation in a larger cohort focused on AA women, as well as replication in more diverse populations, is required to determine whether the observed metabolite differences are specific to this group. This observational trial is underway. Currently, we are also conducting validation studies for enzymes related to the N-acetyl amino acid metabolites and testing them by immunohistochemistry in tissue microarrays representing over 200 BCa biospecimens from AA women with and without BCa. Ultimately, we hope to identify and develop clinically useful biomarkers for earlier detection of BCa, prior to confirmation by tissue biopsy, that can improve clinical care for all women with BCa; as well as identify ones that help explain metabolic contributions that result health disparate outcomes.

The following abbreviations are used in this manuscript:

BCa: Breast Cancer; TNBC: Triple Negative Breast Cancer; AA: African American; CA: Caucasian; 1H NMR: Proton Nuclear Magnetic Resonance Spectroscopy; ER: Estrogen Receptor; PR: Progesterone Receptor; HER2: HER2/neu, Human Epidermal Growth Factor 2 Receptor; HUCC: Howard University Cancer Center; QC: Quality Control; Mammo: Mammoplasty; GI-II: Grade one to two; GIII: Grade Three; HDL: High Density Lipoprotein; LDL: Low Density Lipoprotein; VLDL: Very Low-Density Lipoprotein; TMAO: Trimethylamine N-oxide; MS: Mass Spectrometry

The authors declare they have no conflicts of interest.

Tissue samples were retrieved from Dr. Kanaan’s biorepository, established and maintained with funding support from the Howard University College of Medicine (U400042) Bridge Funds and the Pilot Study Awards Program (BFPSAP), as well as the NIH/NCI Howard/Hopkins Partnership (2U54 CA091431-06). Metabolomics data acquisition and analysis were funded by the NIH Common Fund Metabolomics Phase 1 program (U24DK097193, Sumner PI), and NC State Funds to the Sumner Lab at UNC Chapel Hill.

Y.M.K., R.L.C. and T.J.N. conceptualized the study. D.A.S., Y.M.K., R.L.C., T.J.N., R.L.D., L.R.S., V.T.F. and S.J.S. designed the study. Y.M.K., T.J.N., R.L.D., V.T.F., K.B., E.C-S. and D.A.B. coordinated the clinical studies for the collection and repository of biospecimens from cohorts of African American women at Howard University Hospital, Howard University Cancer Center, and Providence Hospital. D.A.S., L.B., W.W.P. and S.L.M. prepared samples, acquired the data and performed the preliminary analyses. L.B., D.A.B. and E.C-S. organized patient sample metadata. S.L.M. and L.R.S. performed metabolomics data and statistical analyses used to generate figures and tables. D.A.S., L.B. and W.W.P. reviewed and summarized results to draft the preliminary manuscript. D.A.S., L.R.S. and Y.K.M. completed the manuscript. Y.M.K and S.J.S. secured necessary funding for this project. All authors

The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Boards at Howard University (IRB-12-CMED-39) and Providence Hospital (IRB-2014-0901).

Written informed consent was obtained from all participants enrolled in cohort studies, from whom fresh excess surgical breast tissue samples, as well as clinical and pathological data, were collected. The study protocols were approved by the Howard University IRB and Providence Hospital IRB.

All raw and normalized metabolomics data, analytical parameters, and associated metadata supporting the conclusions of this study are publicly available through the NIH Common Fund Data Repository and Coordinating Center Metabolomics Workbench (http://www.metabolomicsworkbench.org/; Project ID PR000338; http://dx.doi.org/10.21228/M8X01N), breast tissue dataset (Study ID: ST000438).

The authors thank Mr. Zachery Acuff, Mr. Andrew Novohatny, and Mrs. Tammy Cavallo for their assistance with manuscript preparation.