PDF

Theranostics 2022; 12(8):3794-3817. doi:10.7150/thno.68611 This issue Cite

Research Paper

Annexin-A1 deficiency attenuates stress-induced tumor growth via fatty acid metabolism in mice: an Integrated multiple omics analysis on the stress- microbiome-metabolite-epigenetic-oncology (SMMEO) axis

Jianzhou Cui^1,2,3, Karishma Sachaphibulkij^1,2,3, Wen Shiun Teo^1,2,3, Hong Meng Lim^1,2,3, Li Zou⁴, Choon Nam Ong^4,5, Rudi Alberts^1,6, Jinmiao Chen^1,7, Lina H. K. Lim^1,2,3

1. Immunology Translational Research Program, Yong Loo Lin School of Medicine, National University of Singapore (NUS), Singapore
2. Department of Physiology, Yong Loo Lin School of Medicine, NUS, Singapore
3. NUS Immunology Program, Life Sciences Institute, NUS, Singapore
4. Saw Swee Hock School of Public Health, NUS, Singapore
5. NUS Environmental Research Institute, NUS, Singapore
6. Department of Microbiology and Immunology, Yong Loo Lin School of Medicine, NUS, Singapore
7. Singapore Immunology Network (SIgN), Agency for Science, Technology and Research (A*STAR), Singapore

✉ Corresponding author: L.H.K.L. (linalimedu.sg) of the Immunology Translational Research Program #03-06J, 28 Medical Drive, Centre National University of Singapore, Singapore 117456.

Abstract

Background: High emotional or psychophysical stress levels have been correlated with an increased risk and progression of various diseases. How stress impacts the gut microbiota to influence metabolism and subsequent cancer progression is unclear.

Methods: Feces and serum samples from BALB/c ANXA1^+/+ and ANXA1^-/- mice with or without chronic restraint stress were used for 16S rRNA gene sequencing and GC-MS metabolomics analysis to investigate the effect of stress on microbiome and metabolomics during stress and breast tumorigenesis. Breast tumors samples from stressed and non-stressed mice were used to perform Whole-Genome Bisulfite Sequencing (WGBS) and RNAseq analysis to construct the potential network from candidate hub genes. Finally, machine learning and integrated analysis were used to map the axis from chronic restraint stress to breast cancer development.

Results: We report that chronic stress promotes breast tumor growth via a stress-microbiome-metabolite-epigenetic-oncology (SMMEO) axis. Chronic restraint stress in mice alters the microbiome composition and fatty acids metabolism and induces an epigenetic signature in tumors xenografted after stress. Subsequent machine learning and systemic modeling analyses identified a significant correlation among microbiome composition, metabolites, and differentially methylated regions in stressed tumors. Moreover, silencing Annexin-A1 inhibits the changes in the gut microbiome and fatty acid metabolism after stress as well as basal and stress-induced tumor growth.

Conclusions: These data support a physiological axis linking the microbiome and metabolites to cancer epigenetics and inflammation. The identification of this axis could propel the next phase of experimental discovery in further understanding the underlying molecular mechanism of tumorigenesis caused by physiological stress.

Keywords: 16S rRNA gene sequencing, gut microbiome, metabolic level, epigenetic signature, PICRUSt, differentially methylated regions, machine learning, restraint stress, serum metabolites, feces metabolites, WGBS, breast cancer, tumorigenesis

Introduction

Breast cancer, among women alone, accounts for nearly 30% of all cancer diagnoses, with an estimated 281,550 new cases and close to 43,600 deaths in 2020 [1]. Epidemiological and clinical studies have provided strong evidence for links between chronic stress and breast cancer progression [2]. Stress is becoming an increasingly inevitable part of people's lives in contemporary societies. An integrated definition describes stress as a constellation of events that involves a stimulus that precipitates a reaction in the brain and, in turn, activates physiological fight-or-flight responses in the body[3-5]. Stress has been proposed to modulate the metabolic, transcriptional, and epigenetic regulation of various diseases, including cancers [6-9], and affects cancer progression through the suppression of immunity and the exacerbation of chronic inflammation [10-12].

Stress factors are also linked to the disruption of the commensal microbiome homeostasis. Given that the commensal microbiota is a host-intrinsic regulator of systemic immune functions [13], it is not surprising that stress alters microbiota diversity, via a process known as dysbiosis, influencing the systemic immune environment [14], which is associated with a poorer outcome in multiple diseases [15-17]. Perturbations in the diversity of the microbiome can affect the abundance of specific bacterial species and increase the risk of disease [18-20].

Epigenetic changes in tumors are emerging as fundamental regulators of breast cancer development and progression [21]. During carcinogenesis, chromatin structure alterations mainly include DNA chemical modification such as CpG methylation and post-translational modification of DNA bound proteins, including histones [22]. Together with subclonal mutations and signals from the microenvironment, methylation modulates the cancer cell phenotype and affects the metastatic propensity of the tumor [23-25]. The direct connections between metabolism and chromatin dynamics now present important conceptual challenges to explain many aspects of tumorigenesis [26]. Changes to intracellular metabolism can alter the expression of specific histone methyltransferases and acetyltransferases, conferring widespread variations in epigenetic modification patterns [27, 28].

Direct connections between the microbiome, metabolism, and chromatin dynamics present important conceptual challenges to explain many aspects of tumorigenesis. Several components of the epigenetic machinery require intermediates of cellular metabolism for enzymatic function. Furthermore, changes to intracellular metabolism can alter the expression of specific histone methyltransferases and acetyltransferases, conferring widespread variations in epigenetic modification patterns.

Studies on one or two aspects of our analysis, such as microbiome, metabolism, or epigenetics in breast cancer, have been reported. However, systemic analysis of the effect of stress on breast cancer through multi-omic approaches has not been performed. The overall purpose of this study was to determine the links between stress and breast tumorigenesis via the gut microbiome, gut-serum metabolite levels, epigenetic signatures, and transcriptomic expression in tumors using integrated analysis and machine learning. Hereby, we hypothesize that stress hormones (e.g., cortisol) alert the gut microbiome composition via the brain-gut axis and causes changes in metabolites in the gut and blood. Stress will also further disturb epigenetic signatures and transcriptome profiling in mice tumors in a remote manner.

In particular, Annexin-A1 (ANXA1) is a glucocorticoid regulated gene which plays an essential role in breast cancer development both in vitro and in vivo. ANXA1 is implicated in multiple functions essential in cancer, including cell proliferation, apoptosis, chemosensitivity, metastasis, and invasion [55]. ANXA1 associates with, and regulates NF-κB, and increases c-Myc activity to promote breast cancer migration and metastasis [56, 57]. Importantly, ANXA1 deficient mice exhibit reduced tumor growth and enhanced survival [58]. In this study, we hypothesized that ANXA1 is involved in stress-related promotion of breast cancer development and determined its role in the regulation of microbiome, gut-serum metabolites.

To our knowledge, our study suggests a novel stress-microbiome-metabolite-epigenetic-oncology (SMMEO) axis that is induced by stress.

Results

Stress promotes breast tumor development

BALB/c mice were subjected to 10 days of restraint stress with one day of rest before orthotopic injection of murine 4T1-luc breast cancer cells into the mammary gland (Figure 1A). Stressed mice exhibited higher levels of luciferase activity seven weeks after stress (Figure 1B), and tumor size was markedly increased after stress from days 34 to 47 after injection (Figure 1C) (p < 0.05). Plasma and feces corticosterone levels were measured at day 0 (before stress, BS), 4 (during stress, S), and 14 (2 days after stress, AS) of the stress procedure (Figure 1A). Stress significantly increased (p < 0.05) the plasma and feces corticosterone levels on days 4 and 14 (Figure 1D). To further study the effect of stress on tumor development, multi-omics analyses were performed to investigate the underlying molecular mechanism in the process (Figure S1).

 Figure 1 

Chronic stress promotes breast tumor development in vivo (A) Schematic of experimental design and sample collection protocol for stress and non-stress; n ≥ 6 mice in each group. The treatment and sampling time were indicated with Before stress (BS), Stress (S), End of stress (ES), and After stress (AS) stages. (B) In vivo images of tumor samples by Xenogen IVIS imaging system in week seven after stress. The left mice were non-stressed; the right mice were stressed by restraining. 4T1 cells were injected after stress. The tumor area is shown in red. (C) Growth curve of tumors seven weeks after stress. ^#p = 0.0216. (D) Fecal and serum corticosterone levels at different stages of acute restraint stress when compared with the non-stressed group. Data are representative of three independent experiments. *p < 0.05, **p < 0.01, ***p < 0.001 in one-way ANOVA with Tukey's post hoc test.

Stress alters gut microbiome composition in mice

16S rRNA gene sequencing was used to investigate the restraint-stress-induced changes in the microbiome that could lead to increased tumor development in mice. Ace and chao1 analysis of the fecal microbiome of non-stressed (NS) vs. stressed (S) mice showed that the stressed group had higher microbial alpha-diversity, even though the difference was not significant (Figure 2A). Principal component analysis (PCA) demonstrated that the fecal microbiome in NS mice clustered separately from S mice (Figure S2A), suggesting both qualitative and quantitative differences in the two groups' gut microbiomes. Increased bacterial loads at the phylum, class, order level were observed in stress stages in stressed mice (Figure S2B-G), implying that stress regulated the gut microbiota response and altered the gut microbiome. Stress decreased the Firmicutes/Bacteroidetes (F/B) ratio (Figure 2B) in the S and AS stages. A Circos diagram shows the composition and abundance of the gut microbiome were different among the before stress, stress, and after stress stages (stress vs. non-stress) (Figure S2H). Compared with the even distribution of each OTU in the before stress group, OTU61 (o_Rhodospirillales), OTU59 (Clostridiales), and OTU62 (o_Clostridiales) were of highest abundance in the stress and after stress groups. The taxonomy of the differently distributed OTUs is listed in Table S1, while the microbial composition and disease link in stress is shown in Table S2. Three major taxonomic units, Clostridiales, Rhodospirillales, and Gastranaerophilales, are closely related to cancer.

 Figure 2 

Altered gut microbiome composition and metabolome in stress (A) Alpha-diversity (ace and chao1) between the NS and S groups at three stages (D0, D4 and D14). (B) The ratio of Firmicutes and Bacteroidetes in NS and S groups at three stages. (C) Linear discriminant analysis effect size (LDA) between the NS and S groups at D4 stage. Red, bacterial taxa statistically overrepresented in NS samples; green, bacterial taxa overrepresented in S samples. The length of the bar represents log10 transformed LDA score. The absolute values of the effect size provide an interpretation of the scale of the difference between two groups for a certain taxon. (D) Taxonomic cladogram obtained from LEfSe analysis showing differentially abundant bacteria taxa between NS and S groups at D4 stage. Green represents increased abundance in the S group; red is increased abundance in the NS group. (E) COGs of PICRUSt analysis between NS and S groups at three stages. In this figure and below, BS, S, and AS were used to indicate the three stages when the S group compared with the NS stages (D0, D4, and D14). n ≥ 6 mice in each group.

A linear discriminant analysis (LDA) effect size (LEfSe) was applied to compare the microbiomes, demonstrating fine-scale differences in the bacterial taxon abundances between the non-stress and stress groups. We found seven differentially abundant clades (α = 0.05), including g_Parabacteroides in the stress mice and one differentially abundant clade in non-stress mice (Figure 2C). The cladogram (Figure 2D) summarizes the LEfSe associations by representing the taxonomic relationships and abundances between the two groups. Five taxa (a, b, d, e, f), in particular, significantly represent two certain phylogenies in the stress group. Stress mice exhibited a high abundance of Lachnospiraceae UCG-001, Mucispirillum, and Rikenellaceae RC9 gut group. In contrast, the non-stress group showed high levels of Anaeroplasma and Ruminococcaceae UCG-014 (Figure S2I).

After sequencing, we predicted the functional genes in the fecal samples using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) analysis. The differences in the clusters of orthologous genes and pathways between the microbial communities were analyzed using the Clusters of Orthologous Genes (COG) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. After normalization, protein function classification demonstrated that stress stage samples had a lower abundance of all 24 COG gene families compared to the before stress microbial communities, but the abundance was restored in the after-stress stage. For example, COGs involved in chromatin structure and dynamics were more highly abundant in the regression period (after stress) (Figure 2E). Three COGs, (RNA processing and modification, cell motility, and the cytoskeleton), gradually decreased during stress from before stress to after stress.

Stressed microbiota were involved in diverse pathways at all three KEGG orthology levels (levels 1, 2, and 3). As shown in Figure S2J, metabolism was significantly enriched during stress compared to before stress and after stress stages in the level 1 category. A heatmap obtained from the level 2 analysis (Figure S2K) shows that pathways involved in carbohydrate metabolism, environmental adaption, cell motility, membrane transport, and transcription were highly enriched in samples taken before stress, while pathways classified to the digestive system and excretory system were more enriched in stress and after stress stages, and pathways involved the endocrine system, cancers, immune system, cell growth and death, were highly enriched after stress. Level 3 category results indicated that most annotated pathways were evenly distributed among all samples, with some outliers indicated (Figure S2L). Ether lipid metabolism was significantly enriched before stress, while a significantly greater proportion of genes involved in Parkinson's disease; various types of N-glycan biosynthesis and apoptosis was observed in after stress. Collectively, these results suggest that stress not only altered the composition of the gut microbiome but also changed the intestinal microbiome response which may lead to enhanced tumor growth.

Stress alters fecal and serum metabolome

Next, we wished to understand the metabolite status in feces and serum after stress and to study the relationship between changes in the composition of the microbiome and metabolism. First, untargeted profiling GC-MS analysis of the metabolites from fecal and serum samples of the stress group was performed. A total of 45 and 44 different metabolites were identified from the fecal and serum samples, respectively (Table S3 and Table S4). PCA analysis demonstrates the obvious clustering and separation of the before stress, stress, and after stress stages according to the fecal and serum samples (Figure 3A-B).

The different fecal and serum metabolite levels in before stress, stress, and after stress, between the NS and S groups were compared using metabolite cluster analysis [29]. In the fecal samples (Figure 3C), the majority (35/45, 77.8%) of metabolites (Feces Cluster 1) increased during stress but reduced after stress, demonstrating a normal stress response. Approximately one-fourth (12/45, 26.6%) of the metabolites exhibited an adaptation response to stress (Feces Cluster 2), increasing from before stress to stress and remaining elevated after stress. In particular, some monounsaturated fatty acids, e.g., palmitoleic acid, linolenic acid, lactic acid, were up regulated in the stress and after stress fecal samples. Similarly, in the serum samples, the majority of metabolites showed the Cluster 1 pattern (33/44, 75.0%) (Figure 3D), while 11 (25.0%) metabolites followed the Cluster 2 pattern. Seven metabolites were common to the Cluster 2 patterns in both the serum and feces samples (Figure 3E), which suggested that these altered metabolites exhibited a conserved pattern between feces and serum samples during stress. KEGG pathway analysis revealed that these metabolites are enriched in biosynthesis of unsaturated fatty acids, galactose metabolism, fatty acid biosynthesis, and linoleic acid metabolism (Figure 3F).

 Figure 3 

Altered metabolome in stress (A, B) PCA analysis of fecal and serum samples at three stages between NS and S samples. (C) Cluster analysis of fecal metabolites of S vs. NS samples at BS, S and AS stages; p < 0.01 Numbers in Parenthesis stand for the numbers of metabolites in the given cluster. (D) Cluster analysis of serum metabolites of S vs. NS samples at BS, S and AS stages; p < 0.01 (E) Venn diagram showing the overlap between Cluster 2 metabolites from fecal and serum samples. (F) KEGG analysis of the overlapping metabolites of Cluster 2 from fecal and serum samples according to MetaboAnalyst 5.0 [30]. The pathway whose p-value was below 0.05 was chosen as a potential target pathway. n ≥ 6 mice in each group.

The gut microbiome and metabolites show a positive correlation during stress

A correlation matrix was generated to explore the functional correlation between the gut microbiome and fecal metabolite changes. Clear correlations were identified between the perturbed gut microbiome and altered metabolite profiles (p < 0.01). As shown in Figure S3A, the relative abundance of Cluster 2 fatty acids positively correlated with g_Candidatus Saccharimonas and g_Parabacteroides, (p < 0.05), while it negatively correlated with g_Ruminiclostridium 6, g_Anaeroplasma (p < 0.05) (Table S5). The high correlation indicates that these specific bacterial taxa may enhance fatty acid metabolism pathway during stress. Further, the metabolites in feces and serum are well conserved in terms of function and expression (Figure S3B). In particular, the fatty acid metabolites (Cluster 2 pattern) in feces were positively correlated with fatty acid metabolites in serum samples from stress mice (p < 0.01) (Figure 4A-B and Table S6).

 Figure 4 

Correlation analysis between the microbiome and metabolites (A-B) Correlation analysis between fecal fatty acids and serum metabolites. There was a high degree of statistically significant concordance between fecal fatty acids and serum fatty acids, with p < 0.01. (C-E) Variable importance of metabolites from S serum samples using three machine learning approaches. (C) Coefficient magnitudes generated by generalized linear modeling (GLM) with 3-fold cross-validation; training AUC = 1 and testing AUC = 0.96. (D) Scaled importance generated by gradient boosting machine (GBM) with 3-fold cross-validation; training AUC = 1 and testing AUC = 1. (E) Scaled importance generated by distributed random forests (DRF) with 3-fold cross-validation; training AUC = 1 and testing AUC = 0.96.

Three machine learning methods, namely generalized linear modeling (GLM), gradient boosting machine (GBM), and distributed random forests (DRF), were used to predict potential metabolite markers in serum samples from S and NS mice. All three methods reasonably predicted that L-methionine, sucrose, alanine 1/2, L-ornithine, threonine, and glucose are candidate biomarkers for the stress response (training area under the curve [AUC] = 1 and testing AUC ≥ 0.96) (Figure 4C-E). Therefore, stress exposure induced a significant taxonomic perturbation in the gut microbiome via the brain-gut axis, which substantially altered the metabolomic profile of the gut/feces metabolites and, thus, the serum metabolites.

Stress affects the epigenetic signature in mammary tumors

Genome-wide mCpG site analysis

As the gut microbiome is distant from the site of breast tumors, we propose that a potential systemic influx of microbiome-regulated metabolites from the feces to the blood may affect the epigenetic signature in tumors during stress by chromatin dynamics, whose activity is directly dependent on metabolites [26, 31]. Bisulfite sequencing enabled the acquisition of the genome-wide DNA methylation landscape at the single-base resolution in tumors from NS and S mice. There were three contexts of C bases on the genome, CG, chlorhexidine gluconate (CHG), and CHH, where H represents any base of A, T, C. We first identified the percentage of methylated cytosines in each context. The average ratios of methylated mCG to total CG were 68.2% and 65.7% in the NS and S samples, respectively (Figure S4A). Dramatically low methylation statuses were found for the CHG and CHH sites. Among these detected mC sites, mCHH represented the highest proportion (~50%), mCpG (CG sites) made up a moderate proportion (~35%) and mCHG accounted for the smallest proportion (`15%) of methylation sites (Figure S4B). We observed significant differences (p < 0.05) between non-stressed and stressed samples for all three mC contexts.

Next, we calculated the methylation density and average methylation level in each gene segment or transcript in the promoter, 5′UTR, exon, intron, and 3′UTR areas. A higher density of mCHG and mCHH methylation was found in the five genomic areas in the S group, particularly in exons (Figure S4C). In comparison, a higher mCG methylation density was observed in almost all five genomic regions in the NS group. Compared with mCHG and mCHH, mCG showed the highest average methylation level in all five regions in both groups (Figure S4D). A significantly higher average CpG site methylation level was observed in the NS group compared to the S group (p < 0.0001) (Figure 5A).

Differentially methylated regions (DMR) and Gene (DMG) analysis

A total of 2,225 DMRs were identified under the CG context that included 254 significant DMGs (p < 0.01, |Log₂FC| ≥ 1) (Table S7) in NS vs S tumors. Among these DMGs, 106 had down-regulated methylation levels, and 148 had up-regulated methylation levels in the S tumors. The volcano plot in Figure 5B provides a general overview of the DMGs between the two groups of samples. We next investigated the number of mC sites in each given DMG. Lurap1l, Ick, Plxna2, Eva1a, Nos1, Nfatc1, Klf4, Syngap1, Inpp4b, and Gm37564 had the highest numbers of mC sites (Figure 5C). The distribution of DMRs in each genomic feature is presented in the Sankey diagrams (Figure S4E), which show the genomic region and chromosome distributions of the DMGs.

We divided these DMGs genes into those that were up or downregulated in the S group and performed separate GO category and KEGG pathway analyses. The genes of the up and down-regulated groups were enriched to very different functional categories. MAPK signaling pathway, pathways in cancer, and breast cancer were associated with upregulated DMGs (Figure 5D), while arginine and proline metabolism and arginine biosynthesis were associated with down-regulated DMGs in S tumors. In terms of GO enrichment analysis, histone methyltransferase complex was key in the upregulated group, and another methylation-related GO item, m7G (5') pppN diphosphatase activity, was associated with a down-regulated DMG set (NUDT10) from S tumors. Our data suggest that DNA methylation changes are involved in regulating metabolism, histone methyltransferase, and other critical breast cancer signaling pathways.

Stress induces genome-wide gene expression changes in tumors

We performed RNA-seq analysis on NS and S breast tumor samples to identify the functional consequences underlying the negatively and positively altered DMGs under stress. We obtained 930 DEGs, including 676 significantly up-regulated genes and 254 significantly down-regulated genes, in the S vs. NS tumors. A volcano plot is used to visualize the results, in which the top 13 significantly changed genes were plotted (Figure 5E). A heatmap is also used to compare the levels of DEGs between non-stress and stress groups (Figure S4F).

 Figure 5 

Stress affects the epigenetic signature and gene expression changes in mice breast tumors. (A) Distribution of the methylation rate of mCpG sites on the whole genome. (B) Summary of DMR-related gene (DMG) results. Volcano plot representation of DMGs in the S and NS tumor samples. The X-axis shows Log2fold-changes in methylation level, and the y-axis is the -Log 10 p-value of a DMG. (C) The top 10 DMGs for mC count. (D) GO and KEGG enrichment analysis of DMGs in the S tumor samples. Biological processes (BP), cellular components (CC), and molecular function (MF) (E) Venn diagram showing the co-expression and up/down-regulation of genes in NS and S tumor samples (top). Volcano plot of RNA-seq transcriptome data displaying significantly differentially expressed genes (DEGs) from tumor samples with or without stress (bottom). Significant DEGs (FDR-corrected P ≤ 0.05) are highlighted in blue, with the grey lines representing the boundary for identification of up- or down-regulated genes (Log2 FC>1). Selected top high-expression genes related to stress response are indicated. (F, G) KEGG pathway analysis of up- and down-regulated DEGs significantly enriched in functional categories (P ≤ 0.05). The gene ratio is the ratio of the DEG number and the number of all genes in a certain enrichment pathway. The dot size denotes the number of DEGs, while colors correspond to the adjusted p-value range. All transcriptome experiments were performed in biological triplicate. One-way analyses of variance (ANOVA) were used to determine the inter-group differences between two groups for one or two variables (***p < 0.001).

To further analyze the functional variation in the DEGs, we used GO and KEGG analyses. The most significantly enriched GO terms in the BP, CC, and MF categories in the up-regulated DEGs were regulation of membrane potential, cell-cell junction, and inorganic cation transmembrane transporter activity. Downregulated DEGs were significantly enriched in pathways related to inflammation (Figure S4G-H). KEGG pathway enrichment analysis showed that the up-regulated DEGs were mainly involved in fatty acid biosynthesis and metabolism and cancer (Figure 5F), while the down-regulated DEGs were again enriched in pathways and diseases related to inflammation (Figure 5G). In addition, gene set enrichment analysis (GSEA) revealed that the relatively high-ranking pentose and glucuronate interconversions, fat digestion and absorption, and steroid hormone biosynthesis gene sets were positively associated with stress; while three gene sets, cytokine-cytokine receptor interactions, glycosaminoglycan biosynthesis chondroitin sulfate dermatan sulfate, and primary immunodeficiency, were negatively associated with stress (Figure S5).

DNA methylome and transcriptome analyses reveal the correlated regulation of hub genes in stress-promoted tumorigenesis

DNA methylation at gene regulatory regions is usually considered to influence transcript expression levels [32]. We merged the DMGs and DEGs and identified 13 genes that showed significantly differential methylation and gene expression (Figure 6A). According to the methylation and gene expression patterns, the overlapping genes between the DMGs and DEGs were classified into three distinct groups (denoted G1, G2, and G3). Genes in G2 (Tbc1d9) and G3 (Cdh10, Lrrc4c) followed the canonical model (a negative correlation between promotor methylation and gene expression) (R = -1, p = 0.027) (Figure 6B) [33], with Tbc1d9 downregulated and Cdh10 upregulated in stress tumors, implying that methylation may play a direct role in the regulation of transcriptomic-level phenotypes in stress-induced tumor development. In contrast, genes in G1 showed a non-canonical and weak correlation pattern (R = 0.44, p = 0.2) (Figure S6A). KEGG pathway analysis demonstrated that G2 and G3 genes are enriched in cancer (Figure 6C). G1 genes were enriched in the MAPK signaling pathway and the oxytocin signaling pathway (Figure S6B). We validated the expression of Cdh10 and Tbc1d9 in breast cancer subsets using the TCGA dataset (http://tcgaportal.org/TCGA/Breast_TCGA_BRCA). Cdh10 expression was not different throughout the subsets, while a lower expression of Tbc1d9 (p < 0.001) was observed in the ER-negative basal-type and HER2+ breast cancers (p < 0.05) (Figure 6D-E). Using Kaplan-Meier survival curves (https://kmplot.com/analysis/) for distant metastasis free survival (DMFS) for breast cancer, higher expression of Cdh10 was associated with worse DMFS. In contrast, lower expression levels of Tbc1d9 were associated with prolonged DMFS (p < 1e^-16) (Figure 6D-E and Figure S7). Similarly, higher expression of Cdh10 was associated with significantly worse post-progression survival (PPS, p < 1e^-16) in gastric cancer datasets and progression free survival (PFS) in lung cancer patients, while lower expression levels of Tbc1d9 in gastric cancer and lung cancer patients resulted in better PPS and PFS, respectively (Figure 6D-E and Figure S7). This clinical data demonstrates a positive correlation to our mouse stress studies, linking stress to DNA methylated genes Cdh10 and Tbc1d9, related to cancer prognosis and survival.

Integrated analysis of multi-omics data reveals the regulation network in stress-enhanced tumorigenesis

We next performed intergraded network analysis combined with metabolites with the DMGs/DEGs data. The metabolite-gene-disease interaction network showed that the list of metabolites identified, such as ornithine, linoleic acid, l-alanine, are involved in Alzheimer disease and lung cancer, amongst others. 25 significantly differently expressed hub genes in the network indicated that this new emergent network is composed of metabolites, genes, and related disease pathways (Figure S8). Neuronal diseases and cancer pathways share some of the same metabolites and genes, showing a strong brain-cancer link. A more detailed gene-metabolite network analysis revealed the relationships between fatty acid metabolism and cancer pathways. Five metabolites and 18 hub genes were identified, and the correlation networks were constructed (Figure 6F). The up-regulated DEGs mainly contributed to the activation of palmitic acid and oleic acid metabolism or synthesis pathways. Two forked family genes (Foxl2 and Foxp2) contributed to L-alanine synthesis pathways, suggesting that the DMGs/DEGs are involved in stress-induced fatty acid biosynthesis during tumorigenesis.

These hub genes are enriched in the inflammatory response, neuroactive ligand-receptor interactions, and glycine, serine, and threonine metabolism pathways (KEGG, Figure S9A). Seventeen hub genes were significantly enriched in GO classes for terms such as positive regulation of metabolic processes, glutathione transmembrane transporter activity (GO, Figure S9B). Interestingly, all the metabolites mentioned above were also directly connected to brain disorders, such as Alzheimer's disease, as well as hyperglycinemia. These results further confirm the concept that the brain-gut axis is a leading player in stress-enhanced tumorigenesis. In addition to the fatty acid metabolic changes in tumorigenesis after stress, the downregulated DEGs were significantly enriched in the inflammation pathways. Based on our current results and our previous findings that ANXA1 is involved in breast cancer growth, we hypothesized that ANXA1 may also play a role in stress-induced tumor development.

 Figure 6 

Integrated gene-metabolite network analysis. (A) Venn diagram comparing DMGs and DEGs. A total of 13 genes were identified from up- and down-regulated DMGs and DEGs. There were three groups identified: G1 (n = 10) genes were hypermethylated in the promoter with high gene expression level in tumor; G2 (n = 1) genes were hypermethylated in the promoter with low expression in tumor; G3 (n = 2) genes were hypomethylated in the promoter with high expression in S tumors. (B) Pearson correlation analysis of G2 and G3 genes; R = -1, p = 0.027. Epi and RNAseq stand for the fold change of given genes in methylation levels and RNA expression in S samples. (C) KEGG pathway analysis of G2 and G3 genes; p > 0.05. (D-E). Kaplan-Meier survival curves show the correlation between marker gene expression and distant metastasis free survival (DMFS), post progression survival (PPS), and progression free survival (PFS) in breast, lung, and gastric cancer. Gene expression in different breast cancer types on the top. FPKM: fragments per kilobase of transcript per million fragments mapped. (F) Metabolite-gene interaction network analysis related to fatty acid metabolism and cancer pathways. The correlation network is composed of five metabolite compounds combined with 18 genes. Metabolites are represented by blue diamonds and genes by circles. Yellow: up-regulated DEGs. Grey: down-regulated DEGs. Red: hypermethylated DMGs. Green: hypomethylated DMGs. Metabolites with KEGG annotations from the merged data set were mapped to KEGG reference pathways, and interaction networks were generated in MetaboAnalyst5.0 (p < 0.005).

ANXA1 deficiency reshapes the murine gut microbiome

ANXA1 knockout (ANXA1^-/-) mice were subjected to restraint stress conditions, and fecal and serum samples were collected for microbiome and metabolism analyses (Figure 7A). ANXA1 deficiency significantly reduced tumor growth to non-existent levels in the NS mice and S mice (p = 0.0051). Even with stress, the ANXA1^-/- mice exhibited significantly smaller tumor volumes (Figure 7B) than wild-type (ANXA1^+/+) mice (Figure 1C). We performed identical analyses in ANXA1^-/- feces as described in Figure 2A. ANXA1^-/- exhibited significantly lower α diversity under NS conditions, suggesting that the gut microbiome in ANXA1^-/- were more sensitive to stress (Figure 7C). β-diversity composition of the fecal microbiome was significantly affected by ANXA1 deficiency (p < 0.001) (Figure 7D). The changes in the fecal microbiota after stress in ANXA1^-/- were explored at different taxon levels (Figure S10A-F). At the phylum level, the F/B ratio in BS ANXA1^-/- feces samples was significantly lower than ANXA1^+/+ BS feces and did not change upon stress (Figure 7E). The enrichment in gut microbiota functions of metabolism, cell growth and death; metabolic diseases; cancers; and immune system diseases in response to stress in the ANXA1^+/+ mice was abolished in the ANXA1^-/- mice (Figure 7F).

A cladogram representative of the structure of the fecal microbiota and their predominant bacteria in NS and S ANXA1^-/- feces is shown in Figure S11A-B, in which the greatest differences in taxa between the two communities are displayed using LEfSe. In addition, tracking individual OTU abundances revealed different dynamics within the dominant microbiota in ANXA1^-/- mice to those of ANXA1^+/+ samples (Figure S11C, Table S8). Microbiota species in the S ANXA1^-/- feces mainly belonged to Saccharimonadia, Lachnospiraceae, and Bacteroidales. These microbes, different from those in the ANXA1^+/+ sample, are lowered in many diseases, including colorectal cancer, obesity and type 2 diabetes (Table S9). Collectively, these changes in the fecal microbiota revealed that deleting ANXA1 expression attenuated the stress-induced gut microbiome dysbiosis.

ANXA1 deficiency influences metabolism and methylation

Using a clustering algorithm to classify the levels of metabolites in all stress stages (as described in Figure 3C-D), two major distinct clusters (C1 and C2) of patterns representing differentially regulated metabolites from fecal and serum samples of the ANXA1^-/- mice. Combined with the clustering data shown in Figure 3C-D, 11 metabolites from ANXA1^+/+ -C2 (enhanced at S and AS, Figure 3C) were suppressed to ANXA1^-/- -C2 (reduced at S and AS), while 15 metabolites from ANXA1^+/+ -C1 (enhanced at S only, Figure 3C) were increased to ANXA1^-/- -C2 in the ANXA1^-/- fecal samples (Figure 8A). The suppressed metabolites in ANXA1^-/- were enriched in the biosynthesis of unsaturated fatty acids, linoleic acid metabolism, and fatty acid biosynthesis, suggesting that ANXA1 alters fecal microbiome metabolism by enhancing fatty acid metabolism during stress (Figure 8B). Moreover, the increased metabolites in ANXA1^-/- included aminoacyl-tRNA biosynthesis, phenylalanine, tyrosine and tryptophan biosynthesis, and phenylalanine metabolism, suggesting that these pathways may be involved in tumor suppression in ANXA1^-/- (Figure 8C). Similarly, we observed clear differences in the clusters between the ANXA1^+/+ and ANXA1^-/- serum samples (Figure 8D). Four metabolites were suppressed in ANXA1^+/+ -C2 to ANXA1^-/- -C1 in serum samples, and 12 metabolites from ANXA1^+/+ -C1 to ANXA1^-/- -C2 were increased in ANXA1^-/- S serum. Similarly, the suppressed serum metabolites were enriched in fatty acid biosynthesis, while the increased serum metabolites were enriched in alanine, aspartate, and glutamate metabolism pathways (Figure 8E-F). When these differentially regulated metabolites in the feces and serum were compared, three pathways: biosynthesis of unsaturated fatty acids, fatty acid biosynthesis, and alpha-linolenic acid metabolism, were commonly suppressed. Three other pathways, aminoacyl tRNA biosynthesis, glycerolipid metabolism, and glycolysis/gluconeogenesis, were commonly increased). Moreover, correlation analysis further demonstrated that all the metabolites in the fecal samples were strongly connected to those in the serum samples from the ANXA1^-/- mice (Figure S12 A).

 Figure 7 

ANXA1 deficiency alters gut microbiome structure under stress (A) Schematic of experimental design and sample collection protocol for ANXA1-knockout mice. n ≥ 6 mice in each group. (B) Orthotopic tumor weight 36 days after injection of 4T1 cells. Results were obtained in three independent experiments. The data are shown as the mean ± SD; **p < 0.01. (C-D) Alpha-diversity and PCA analysis of bacteria varied across ANXA1^+/+ and ANXA1^-/- fecal samples; ***p < 0.001 (E) Ratio of Firmicutes and Bacteroidetes in ANXA1+/+ and ANXA1^-/- fecal samples. (F) KEGG analysis of ANXA1^+/+ and ANXA1^-/- feces samples under stress conditions. (WT, ANXA1^+/+; AKO, ANXA1^-/-).

 Figure 8 

ANXA1 deficiency changes feces and serum metabolites (A, D) Metabolite clustering results (time-series line) based on metabolite expression levels. The horizontal axis denotes stress procedure timepoints (BS, S and AS) at which fecal (A) and serum (D) samples of ANXA1^+/+ and ANXA1^-/- mice were collected. The vertical axis denotes the levels of the metabolites relative to the internal standard. Numbers in brackets indicate the number of metabolites in each cluster. (B) KEGG analysis of the 11 suppressed metabolites transferred from ANXA1^+/+ -C2 to ANXA1^-/- -C2 in fecal samples. (C) KEGG analysis of the 15 increased metabolites transferred from ANXA1^+/+ -C1 to ANXA1^-/- -C2 in fecal samples. (E) KEGG analysis of the 4 suppressed metabolites transferred from ANXA1^+/+ -C2 to ANXA1^-/- -C1 in serum samples. (F) KEGG analysis of the 12 increased metabolites transferred from ANXA1^+/+ -C1 to ANXA1^-/- -C2 in fecal samples. (G) Interaction network constructed with GeneMania for highly altered hub genes and Anxa1. The color of the lines connecting the genes depicts the type of interaction (purple for co-expression, yellow for predicted, blue for co-localization, and pink for physical interaction). The colors in the circles indicate the functions (see key).

Machine learning analysis showed consistent results among three approaches that L-5-Oxoproline, L-Methionine, and Mannose represent the high AUC. The prediction analysis demonstrated the S ANXA1^-/- mice would have a different metabolite composition compared to the ANXA1^+/+ (Figure 4C-E, Figure S12B-D). The machine learning prediction results further demonstrated that ANXA1 deficiency changes the levels of the potential metabolite markers in the S serum samples. Accordingly, ANXA1 deficiency may alert the gut microbiome and regulate fatty acid metabolism in the feces and serum, leading to the remote inhibition of tumorigenesis.

To further investigate the role of ANXA1 in fatty acid synthesis, the expression of fatty acid synthase (Fasn) and the upstream protein-ATP citrate lyase (Acly) were assessed. ANXA1 deficiency suppressed the protein expression of Fasn and Acly, but not mRNA (Figure S13A, B). Furthermore, the effect of ANXA1 deficiency on DNA methylation was assessed. DNA methyltransferase (cytosine-5) 1(Dnmt1) mRNA and protein expression were not significantly different between WT and ANXA1 -/- 4T1 cells (Figure S14A). However, mRNA and protein expression of lysine (K)-specific demethylase 1A (Lsd1) was significantly lower in 4T1 cells deficient in ANXA1 (Figure S14B).

Correlation of ANXA1 with hub genes to regulate tumor development

Finally, as no tumors developed in the ANXA1^-/- mice, we next determined if ANXA1 interacts with hub genes to enhance oncogenesis based on the RNA-seq and epigenetic analysis of the ANXA1^+/+ tumor samples. The correlation network results showed that Anxa1 is involved in vascular processes in the circulatory system, rhythmic processes, enzyme inhibitor activity, and fatty acid derivative transport via interactions with Ocln, Cbs, Ptgds, Pkia, Prok2, Gja1, Kcnma1, Oaz3, and Oxt (Figure 8G). In addition, two upregulated hub DEGs, Plin1 and Oxt, were positively correlated with Anxa1 expression in breast cancer patients (Figure S15A-B). Plin1 was found to function as a modulator of adipocyte lipid metabolism, and Oxt is involved in cognition, tolerance, adaptation, and stress responses [34, 35].

Therefore, we here propose a model that illustrates the entire process through which stress promotes tumorigenesis via modulation of the microbiome composition and regulation of metabolism and resulted in the alteration of epigenetic signatures in our mouse model of breast cancer. ANXA1 deficiency suppresses tumor growth by altering the gut microbiome, regulating serum metabolite levels, and interacting with hub genes. All these changes involved fatty acid metabolism, the inflammatory stress response, and changes to neuro-hormone levels (Graphical Abstract). We have defined this microbiome-mediated metabolite and epigenetic interactions that caused breast cancer cell proliferation and collective feedback loops as the SMMEO axis. Brain (stress)-gut-tumor crosstalk, especially in non-gastrointestinal cancers, remains a crucial area of discovery.

Discussion

Humans are increasingly exposed to environmental factors that influence cancer development as the intensities of their lifestyles increase. Many studies have investigated the impact of stress in the promotion of cancer. Moreover, researchers have also tried to illustrate the impact of the gut microbiota on treatment responses to cancer [36, 37]. To our knowledge, our study is the first to systematically evaluate the role of stress in breast tumorigenesis through a multi-omics approach. This study will assist our understanding in determining how relevant the gut microbiome, metabolism, and tumor epigenetic signatures are to human breast cancer and what interventional strategies could be employed to improve patient outcomes.

Our study has demonstrated that exposure to stress prior to the initiation of mammary cancer contributes to increased tumor growth, and that perturbations in the gut microbiome caused by stress can affect tumor growth at a distant site, which is in line with previous findings that suggest the gut microbiome has endocrine functions [18, 38, 39]. A number of metabolites associated with particular gut bacterial species correlated with DMG expression changes. Therefore, stress signals can disrupt the gut microbiome composition via the brain-gut axis and change gut and blood metabolites. Subsequently, these blood metabolites can mediate epigenetic and gene expression changes in the breast tumor microenvironment (Graphical Abstract).

The F/B ratio is regarded to be of significant relevance in human gut microbiota composition [40], where a lower F/B ratio was observed in breast cancer survivors [41]. Our results showed that stress decreased the F/B ratio, which is consistent with previous studies [42-44]. Interestingly, the profile and abundance of the microbiome may return to the baseline after stress, or the microbiome may adapt to a new state. This latter pattern may continuously affect downstream pathways to promote disease and tumorigenesis. Hence, the composition of the gut microbiome exists in a dynamic equilibrium that becomes more complex after stress [38].

An early study showed that the abundance of Clostridiales is closely related to steroidal estrogen production, contributing to breast carcinogenesis and stimulating tumor growth [45]. More detailed studies demonstrated that Clostridiales increases bile acid degradation in breast cancer patients [46]. A secondary bile acid, lithocholic acid (LCA), reduced breast tumor cell proliferation, aggressiveness, and metastatic potential of primary tumors through mesenchymal-to-epithelial transition and enhancing antitumor immune response [47]. Moreover, Faecalibacterium prausnitzii, a member of Clostridium cluster IV, suppressed the proliferation and invasion and promoted the apoptosis of breast cancer cells via inhibition of the secretion of IL-6 and the phosphorylation of JAK2/STAT3 in MCF-7 cells [48, 49]. The collagenase enzyme of the bacterium Clostridium histolyticum (CCH) was shown to reduce cell proliferation and wound healing in breast cancer MDA-MB-231 cells [50]. In addition, dietary fat affects Clostridiales level while smoking increases Rhodospirillales abundance and impacts the tumor microbiome in lung cancer [51]. Accordingly, different lifestyles, including stress may contribute to the gut microbiome composition and may thus be linked to cancer development.

Besides gut microbiota, breast microbiome can differ between healthy and cancer patients and contributes to the pathology of cancer via regulation of lipid signatures [52]. Accumulating evidence suggests that the microbial composition of breast tissue plays an important effect in cancer development [53-57]. It has been proposed that the endogenous and exogenous microbiome from the breast contributes to the maintenance of normal function of breast tissue either by stimulating resident immune cells or regulating metabolism [58]. Moreover, cholesterol (and lipid) metabolism has been implicated in breast cancer carcinogenesis [59-61]. In addition, cholesterol is a risk factor for breast cancer, through its impact on membrane fluidity and signaling pathways [62]. High plasma cholesterol levels can enhance the proliferation of tumor cells and tumor growth in mouse models [63]. Cholesterol and other lipid metabolites may act as biomarkers for breast cancer development and provide a novel target for cancer therapy.

Short-chain fatty acids (SCFAs) are known to play crucial roles in immune regulation and inflammatory diseases [64-66]. However, the underlying mechanisms of how long-chain fatty acids (LCFAs) influence brain physiology, gut microbiota, gut-brain communication, and tumor development have not been fully elucidated [67]. Cancer cells have an exuberant ability to synthesize fatty acids to maintain and promote cell growth and cancer progression [68]. In our study, stress significantly altered the levels of fatty acids in the feces and serum, which indicates that stress may promote breast tumor development via fatty acid metabolism and biosynthesis. A number of metabolites were directly shown to correlate with different microbiome taxonomic levels and DEGs/DMGs. Notably, the up-regulated DEGs mainly contribute, via the gut-metabolite-gene network, to the biosynthesis of LCFAs palmitic acid and oleic acid. Palmitic acid and oleic acid bridge the link between fatty acid biosynthesis, metabolism and cancer pathways, further supporting the notion that the brain-gut axis plays a vital role in tumorigenesis. Four hub genes, Npffr, Gpr132, Prok2, and Ptgds were commonly connected with palmitic acid and oleic acid, suggesting the functional conservation of these genes in fatty acid metabolism.

Bioactive metabolites such as short-chain fatty acids, amino acid metabolites, or secondary bile acids secreted from the gut microbiome play an essential role in regulating breast cancer development [69]. The fatty acids highlighted in Figure S3A have previously been shown to be involved in breast cancer development. Epidemiological studies have shown that elevated serum concentrations of oleic acid together with low levels of stearic acid are associated with an increased breast cancer risk [70, 71]. High levels of palmitoleic acid to palmitic acid concentrations are associated with an increased breast cancer risk and can enhance cancer development [72-74]. The influence of stearic acid on the inhibition of tumor cells in vitro and tumor development in vivo has been reported [75]. Stearate induces apoptosis preferentially in breast cancer cells this may be protein kinase C dependent [76]. The protective association provided by stearic acid has been previously reported in premenopausal women in a meta-analysis of fatty acids in biological samples and BC risk [77]. In a meta-analysis of 12 prospective studies, both linoleic acid intake and serum levels of linoleic acid were associated with decreased breast cancer risk. However, none of the associations was statistically significant [78]. In animal studies, α-linolenic acids have been shown to suppress the growth and proliferation of cancer cells and promotes breast cancer cell death [79-81]. α-linolenic acid reduces growth and inhibits the migration of both triple-negative and luminal breast cancer cells in high and low estrogen environments [82, 83]. Our results have confirmed that fatty acids play an essential and complex role in breast cancer development.

Similar to SCFAs [84], blood-borne microbial metabolites, including lithocholic acid (LCA) [85, 86], deconjugated estrogens [87], and amino acid degradation products can remotely execute their physiological function. Anaerobic bacteria, such as Clostridiales, are responsible for bile acid transformation from gut origin to breast [88, 89]. The synthesis of LCA from the microbiome was vastly reduced in the advanced stage of breast cancer, pointing towards an antineoplastic role of LCA. LCA inhibited epithelial-to-mesenchymal transition and metastasis formation in breast cancer cells [85], and a high concentration of LCA inhibits fatty acid biosynthesis which could promote cell death [86, 90]. Interestingly, a large-scale metabolomics analysis demonstrated that a high level of microbiome-derived deoxycholate was found in breast tumors which was inversely associated with the cell proliferation in breast tumors, suggesting that this bile acid accumulation in the breast could be used for breast cancer survival prediction. [91]. In terms of amino acid degradation, bacteria synthesize cadaverine from lysine using the enzymes LdcC and CadA[92]. Similar to LCA, the cadaverine expression in gut microbiome is decreased in the later stages of breast cancer[93]. Moreover, cadaverine inhibits cell proliferation, cell growth and invasion, and tumor infiltration to the surrounding tissues by changing metabolism or reducing the proportion of ALDH1+ cancer stem cells [93]. Cadaverine exerts its function through the trace amine-associated receptor-1, 2, 3, 5, 8, 9 (TAAR1, 2, 3, 5, 8, 9), of which TAAR1 has tumor suppressive roles [94]. These metabolites secreted from the microbiome are critical constituents of the tumor microenvironment involved in breast cancer development. The dysregulation of the same metabolic pathways in tumors and the breast tumor microbiome therefore suggest an interconnection between the tumor and the corresponding microbiome.

Previous reports have shown that microbiota can directly affect fatty acid metabolism levels. Multiple metagenomes analyses have demonstrated that g_Candidatus Saccharimonas can negatively regulate short-chain fatty acids metabolism [95]. Candidatus_Saccharimonas were significantly inhibited with high-fat diets (HFDs) treatment in C57BL/6J mice[96]. On the other hand, the abundance of Parabacteroides was elevated after HFD feeding 4-weeks, which suggested that fatty acid level regulated microbe abundance [97]. A strong correlation between gut microbiota and fatty acids has been reported in autistic rats [98]. In terms of specific fatty acid, mice injected with valproic acid showed higher abundance of Candidatus_Saccharimonas which suggests a direct link between gut microbiota and fecal metabolites [99]. In line with our results, a positive regulation between Parabacteroides and SCFAs has been demonstrated in human fecal microbiota identified by in vitro fermentation [100]. In addition, high levels of SCFAs promote Parabacteroides bacteria abundance [101].

The essential amino acid L-alanine, which is related to alanine, aspartate, and glutamate metabolism, is an ingredient of protein synthesis and inflammation [102]. Alterations can cause an imbalance in the energy metabolism and inflammatory responses involved in stress-induced tumorigenesis [103]. In addition to L-alanine, sucrose also showed strong correlations with DEGs and DMGs. A higher total sugar intake increases the risk of breast cancer development was associated with a poorer prognosis after breast cancer diagnosis [104, 105]. These findings suggest that carbohydrate and amino acid metabolism are also involved in stress-enhanced tumor growth.

Finally, our DMG/DEG combined analysis identified 2 genes, CDH10 and TBC1D9, which were upregulated, and downregulated in stressed tumors, respectively. CDH10 is the gene for cadherin 10, which is associated with autism[106]. CDH10 is shown to be highly mutated in colorectal cancer and associated with better survival[107]. However, our study showed that stress enhances CDH10, and a higher expression of CDH10 in breast, gastric and lung tumors is associated with poorer survival. TBC1D9 is a gene implicated as a long-term survival gene in breast cancer [108], its expression can differentiate TNBC (low) from non-TNBC (high) breast cancer samples, and overexpression leads to better prognosis [109], which is similar to our study. This shows that CDH10 and TBC1D9 are important genes in cancer and stress.

Immune-modulatory AnxA1 possesses multiple functions essential to cancer pathogenesis, including cell proliferation, apoptosis, metastasis, and invasion [110, 111]. AnxA1 -deficient mice exhibit reduced tumor growth and enhanced survival in vivo [112], and AnxA1 expression is correlated with the high metastatic ability[113], and ANXA1 modulates the immune response in cancer [114]. Recently, evidence has shown that ANXA1 could be used as a possible new therapeutic avenue via repairing blood-brain barrier damage in metabolic disease [115]. In addition, ANXA1 inhibits obesity, suggesting that ANXA1 improves metabolism in models of metabolically stressed animals [116]. Our current study further demonstrated that silencing ANXA1 slows tumor growth and reverses the effects of stress on the changes in the gut microbiome and fatty acid metabolism. This has not been previously reported and may explain the multiple functions of ANXA1 in cancer.

In conclusion, combining multi-omics analysis with machine learning has led us to define the SMMEO axis that results from the stress induction of microbiome-mediated metabolite and epigenetic interactions that eventually lead to the enhancement of tumorigenesis. One limitation of this study is we did not present the direct relationship between microbiota composition / metabolomic profile changes and tumor growth. Other host-related changes may also play an important role in stress. An additional work showing that microbiota depletion or supplementation orally with the metabolites of interest can induce similar changes in tumor growth would benefit this multi-omics study. The novel microbiome and epigenetics markers are potential therapeutic avenues to preventing cancer in women facing stress and thus have enormous potential for improving treatment outcomes.

Materials and Methods

Mice, Stress exposure and Bioluminescence imaging of 4T1 tumor

All animal work was approved by the NUS Institutional Animal Care and Use Committee (IACUC) and was in accordance with the National Advisory Committee for laboratory Animal Research (NACLAR) Guidelines (Guidelines on the Care and Use of Animals for Scientific Purposes). BALB/c ANXA1^+/+ and ANXA1^-/- mice were housed under a 12-h light/dark cycle with food and water provided ad libitum under pathogen-free conditions in the animal housing unit of the Comparative Medicine Department of the National University of Singapore.

Mice (ANXA1^+/+ and ANXA1^-/- ) were allocated to the NS control or tube-restraint stress groups (n ≥ 6 per group). Mice in the tube-restraint stress group were placed in a ventilated 50 mL polypropylene conical tube (Corning Inc.) and subjected to restraint stress for two h per day (9:00 AM to 11:00 AM), while the NS group mice were left undisturbed in their cages. Mice in the stress group were restrained for 10 consecutive days. The mice were then subcutaneously injected with 50 µL of stably transfected 4T1-luciferase cell suspension (10⁴ cells per mouse) into the mammary fat pad, and mice were monitored for up to 40 days. The size of the tumours and tissue metastasis was measured and analyzed via a bioluminescence imaging assay using the Xenogen IVIS Spectrum Imaging System, along with the machine's software (Comparative Medicine facility at NUS). Mice were injected with 150 µL of luciferin (150mg/kg) VivoGlo™ Luciferin, Promega) intraperitoneal. Tumour volumes were measured manually using a digital caliper and were calculated using the equation: V (mm³) = length (mm) × width (mm) × width (mm)/2). Mice were euthanized either at the end of the study or earlier if they displayed significant weight loss, signs of distress, or palpable tumours ≥2.0 cm in diameter.

16S ribosomal rDNA sequencing and microbiome analysis

The QIAamp DNA Stool Mini Kit was used for fecal DNA extraction in accordance with the manufacturer's protocols, and the DNA (20-30 ng) was used to generate amplicons libraries. To analyze the taxonomic composition of the bacterial community, amplicons containing the V3 and V4 regions of the prokaryotic 16S RNA gene obtained with the primers were selected for the subsequent pyrosequencing (16S Amplicon PCR Forward Primer: 5'-CCTACGGRRBGCASCAGKVRVGAAT-3'; 16S Amplicon PCR Reverse Primer: 5'- GGACTACNVGGGTWTCTAATCC-3'). Then, the library was purified with magnetic beads, the concentration was detected on a microplate reader, and fragment size was detected by agarose gel electrophoresis. The library was quantified to 10 nM, and PE250/FE300 paired-end sequencing was performed according to the Illumina MiSeq/Novaseq (Illumina, San Diego, CA, USA) manual. MiSeq Control Software (MCS)/Novaseq Control Software (NCS) was used to read sequence information.

After quality filtering, VSEARCH clustering (1.9.6) sequence (sequence similarity was set to 97%) was used for OTU clustering. Then the Ribosomal Database Program (RDP) classifier Bayesian algorithm was used to analyze the OTU species taxonomic representative sequences and the different species classification levels. Shannon and Chao1 analyses and PCA results were displayed based on the sample OTU abundances table.

Intergroup difference analysis at the phylum, class, order, family, genus, and species levels in each cluster were analyzed with the LEfSe method with default settings in Galaxy workflow framework (https://huttenhower.sph.harvard.edu/galaxy/root)[117]. LEfSe used the two-tailed nonparametric Kruskal-Wallis test to evaluate the significance of differences between OTUs in the non-stress and stress groups. A set of pairwise tests was performed using the unpaired Wilcoxon test. Finally, LDA was performed to estimate the effect size of each differentially abundant OTU. The results are expressed as the mean ± SEM. The gut microbiotas were considered significantly different if their differences had a p-value of < 0.05 and an LDA score of |log10| > 2[118].

The functions from the prokaryotic clades were first predicted using FAPROTAX [119] and visualized by ImageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/index.html). In addition, we used Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) software (v1.0.0) [120] to characterize the functional genes in the sample through a comparison of the bacterial composition information obtained from the 16S RNA gene sequencing data. The three following steps were performed in the analysis: (1) The selection of closed-reference OTUs from the obtained 16S rRNA gene sequences, comparison of the sequences with the Greengenes database, and the use of “nearest neighbor” in the database as the reference OTU. (2) Normalizing the OTU abundance matrix using the “nearest neighbor” rRNA gene sequence counts as reference. (3) Calculating and predicting the overall COG functions and pathways based on the “nearest neighbor” function profile derived from the KEGG/COG database. Next, we used the G-test and Fisher exact test to test the significance of the difference between two samples and t-test to compare between two groups.

Disbiome database was used to uncover the microbial composition changes in different kinds of diseases, managed by Ghent University [121]. The major taxonomic units in the S group were searched and returned the information related to the experiment (related disease/bacteria, abundance subject/control, control type, detection method, and related literature).

Metabolomics Analysis

Chemicals and reagents. Methanol (MeOH, MS grade), pyridine (anhydrous grade), N-(9-fluorenylmethoxycarbonyl)-glycine (FMOC-glycine), methoxyamine hydrochloride, and N-methyl-N-trimethyl-silyl-trifluoroacetamide (MSTFA) were purchased from Sigma-Aldrich (St. Louis, Missouri, USA). Deionized water was obtained from Milli-Q purification system (Bedford, MA, USA).

Sample preparation. The serum/feces sample preparation method was similar to one we used previously [122]. Briefly, 40 μL of serum was extracted with 280 μL of cold MeOH (FMOC-glycine as internal standard) to precipitate the proteins. Feces samples were extracted in an ice-cold methanol/water mixture (7:1, FMOC-glycine as internal standard) after by placing in a TissueLyser LT (QIAGEN, Germany) for 10 min and sonicating at 25 Hz for 10 min. The mixture was then centrifuged for 20 min at 14,000 rpm and 5 °C, and the resulting supernatant was filtered through a micro-centrifugal filter (Thermo Scientific 750-µL micro-centrifugal filter, PTFE membrane, 0.2-µm pore size, non-sterile). The sample was vacuumed-dried (CentriVap concentrator, Labconco, USA) and derivatized with 100 µL of methoxyamine-pyridine solution (5 mg/mL) for 2 h at 60 °C, and then with 100 µL of MSTFA (40 °C, 16 h) for GC-MS analysis. Pooled quality control (QC) samples were prepared by mixing a certain amount of each serum/feces sample. The QC samples were analyzed at the beginning, the end, and randomly throughout the whole assay to evaluate the stability and reproducibility of the GC-MS analytical system.

GC-MS analysis. GC-MS analysis was performed on an Agilent 7683B Series Injector (Agilent, Santa Clara, CA, USA) coupled with an Agilent 7890A Series gas chromatography system and a 5977B mass detector (Agilent). A fused silica capillary column HP-5MSI (60 m × 0.25 mm i.d., 0.25-μm film thickness) was used, and the injector was kept at 250 °C. A 1-μL aliquot of the sample was pulsed-split injected for each individual analysis. Helium was used as the carrier gas at a constant flow rate of 2 mL/min through the column. The GC oven temperature was maintained at 50 °C for 1 min, then increased to 250 °C at a rate of 8 °C/min, and further increased at 25 °C/min to 300 °C and held for 7 min. The transfer line temperature was kept at 280 °C. Detection was achieved using MS in electron impact mode (70 eV) and full-scan monitoring (m/z 50 to 650). The temperature of the ion source was set at 230 °C, and that of the quadrupole was set at 150 °C.

Metabolite Identification: The spectral data were exported as mzData files and pretreated with the online open-source XCMS (https://xcmsonline.scripps.edu/) for peak detection and peak alignment. Peak area normalization in each dataset was calculated by comparison with the internal standard. All identified metabolites were confirmed with standards or matched to the NIST library in GCMS (>80%, p < 0.05). The low relative standard deviation filtration was less than 30%, and the detection frequency was more than 100%.

Metabolite cluster and metabolomics data analysis: The identified metabolites were clustered by using Cluster Trend tools in Hiplot (https://hiplot.com.cn), a comprehensive web platform for scientific data visualization. The biomarker analyses, enrichment analysis, integrated metabolic pathway analysis, and network analysis were performed by MetaboAnalyst5.0 [123].

Whole-Genome Bisulfite Sequencing (WGBS)

Next-generation sequencing libraries were constructed following the manufacturer's protocol (Illumina, San Diego, CA, USA). For each tumor sample in the NS and S groups, 1 μg of genomic DNA was randomly fragmented to < 500 bp by sonication (Covaris S220). The fragments were treated with End Prep Enzyme Mix for end repairing, 5′ phosphorylation, and dA-tailing in a single reaction, followed by T-A ligation to add methylated adaptors to both ends. Size selection of the adaptor-ligated DNA was then performed using VAHTS DNA Clean Beads (Vazyme, China), and fragments of ~410 bp (with the insert of approximately 350 bp) were recovered. Then bisulfite conversion was performed using EZ DNA Methylation-Gold™ Kit (Zymo Research, CA, USA). Each sample was then amplified by PCR (ETC811 Thermal Cycler (EASTWIN, China)) for 10 cycles using the P5 (5' AAT GAT ACG GCG ACC ACC GA 3') and P7 (5' CAA GCA GAA GAC GGC ATA CGA GAT 3') primers, with both primers carrying sequences that can anneal with patterned flow cell technology (Illumina) to perform bridge PCR, and P7 primer carrying a six-base index that allows for multiplexing. The PCR products were cleaned up using VAHTS DNA Clean Beads (Vazyme, China), validated using an Qsep100 DNA analyzer (Bioptic, China), and quantified on a Qubit3.0 Fluorometer (Invitrogen, Carlsbad, USA). Then, libraries with different indices were multiplexed and loaded onto an Illumina instrument (NovaSeq 6000) according to the manufacturer's instructions (VAHTS Universal Pro DNA Library Prep Kit for Illumina, Catalog# ND608). Sequencing was carried out using a 2 × 150 paired-end configuration, and image analysis and base calling were conducted using the Illumina pipeline (image control software [HCS/MCS] + OLB + GAPipeline-1.6).

Data analysis with Cutadapt (V1.9.1) was performed to remove the sequences for adaptors, PCR primers, those containing more than 10% N bases, and bases of a quality lower than 20. Bismark (V0.7.12) was used to map clean data to a reference genome and determine the number of mC sites with the coverage2cytosine command. Differentially methylated cytosines between the S and NS groups were detected by methylKit (V0.9.5), and swDMR (V1.0.7) was used to reveal a set of DMRs. Enrichment for target genes was completed with the GO/KEGG database. We obtained 2.48-2.49 × 10⁹ uniquely mapped reads among all the samples to ensure concordant coverage. The average ratios of uniquely mapped reads for the NS and S tumor samples were 91.21% and 91.35%, respectively. The methylation level was calculated as an average of 100,572,417 and 113,164,957 methylated cytosines (mCs) in the NS and S samples, respectively.

RNA-sequencing

Library preparation for transcriptome sequencing. A total of 1 μg of RNA per tumor sample from the NS and S groups were used as input material for RNA sample preparation. Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit from Illumina, following the manufacturer's recommendations. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. RNA strand fragmentation was carried out using divalent cations under an elevated temperature (42 °C) in NEBNext First Strand Synthesis Reaction Buffer (5X) or using sonication with the Diagenode Bioruptor Pico. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, the dTTP in the dNTPs was replaced with dUTP. To preferentially select cDNA fragments of 250-300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μL of USER enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. The PCR was performed with 0.2 µl Phusion High-Fidelity DNA polymerase (Catalog # M0530S, NEB), universal oligo-dt primer [d(T)23VN]. Lastly, the products were purified with the AMPure XP system, and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Data analysis

After cluster generation, the prepared libraries were sequenced on an Illumina platform, and paired-end reads were generated. Raw data (raw reads) of FASTQ format were first processed using fastp. In this step, clean data (clean reads) were obtained by removing reads containing adapter and poly-N sequences and low-quality reads from the raw data. Paired-end clean reads were aligned to the reference genome using Spliced Transcripts Alignment to a Reference (STAR) software (https://github.com/alexdobin/STAR).

Differential expression and enrichment analysis. Differential expression analysis between the NS and S groups (three biological replicates per condition) was performed using the DESeq2 R package. DEGs were identified using the DESeq2 R package based on the read counts at the transcriptional level, which were identified using the absolute log2FC value ≥1 and adjusted p-value of < 0.05 as the statistical standards. DESeq2 provides statistical routines for determining differential expression from digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). Genes with an adjusted p-value of < 0.05 found by DESeq2 were assigned as differentially expressed.

GO enrichment analysis and biological pathway analysis were carried out using REACTOME (http://www.reactome.org/PathwayBrowser) and Enrichr (https://maayanlab.cloud/Enrichr/#) [124]. MetaboAnalyst 5.0 [123] was used for functional pathways, and the default setting was used for enrichment analysis of metabolites. The gene interaction network was created using highly altered hub genes with the GeneMania Prediction Server [125].

ELISA

Feces and blood samples collected from the NS and S groups were kept at -80 °C till use. The Corticosterone ELISA Kit (Cayman,501320, Ann Arbor, USA) was used to measure corticosterone levels according to the manufacturer's instructions. Samples were extracted before use in ELISA. Briefly, ELISA buffer (100 µl) and a corticosterone standard were added to a 96-wells plate, followed by 50 µl of sample per well at a minimum of two dilutions. After that, 50 µl each of AChE tracer and ELISA antiserum were added to the wells, except for the blank control. The plate was covered with plastic film and incubated overnight at 4 °C. To develop the color reaction, Ellman's reagent was added, followed by wash buffer. Lastly, the color in the plate was read at 412 nm wavelength.

Correlation and survival analysis

Correlation analysis between the microbiome and metabolites was performed using the OmicStudio tools at https://www.omicstudio.cn/tool/62. R version 3.6.1 (Pearson and spearman, 2019-07-05), ggplot2 (3.3.2), and heatmaply (http://talgalili.github.io/heatmaply/)[126]. Machine learning was performed by h2o-genmodel (https://docs.h2o.ai/h2o/latest-stable/h2o-docs/flow.html). The General Linear Model (GLM is just the sum of the coefficient times the value and then adjusts the threshold. Gradient boosting machine (GBM) and distributed random forests (DRF) takes the mean of 40 forests, then compare with threshold and decide 1 or 0. Correlations in the expression of selected DMGs and DEGs were revealed by the interactive web tool GEPIA [127]. Overall survival based on gene expression was analyzed with GEPIA (http://gepia.cancer-pku.cn/) and Kaplan Meier plotter (http://kmplot.com/analysis/index.php?p=background).

Data Analysis

Samples sizes were calculated using practical meta-analysis effect size calculator (https://www.campbellcollaboration.org/escalc/html/EffectSizeCalculator-Home.php) and G*Power version 3.1.9.7 (https://stats.idre.ucla.edu/other/gpower/#) (power > 0.8)[128]. Three independent experiments were performed, and a two-tailed Student's t-test was used to determine the statistical significance of pro-inflammatory gene expression between the control and treated mice. Pearson correlation coefficient for significant DMGs and DEGs was calculated based on their respective intensity values using the CorrelationCalculator v1.0.1, in MetScape 3.1[129]. One-way and two-way analyses of variance (ANOVA) were used to determine the inter-group differences between more than two groups for one or two variables. The level of statistical significance was taken to be p < 0.05 throughout the study.

Abbreviations

ANOVA: Analysis of variance; AUC: Area under the curve; SMMEO: Stress-Microbiome-Metabolite-Epigenetic-Oncology; CHG: Chlorhexidine gluconate; COG: Clusters of orthologous genes; DEG: Differentially expressed genes; DMFS: Distant metastasis-free survival; DMR: Differentially methylated regions; DRF: Distributed random forests; FC: Fold change; GBM: Gradient boosting machine; GC-MS: Gas chromatography-mass spectrometry; GLM: Generalized linear modeling; GO: Gene ontology; GSEA: Gene set enrichment analysis; KEGG: Kyoto encyclopedia of genes and genomes; LEfSe: Linear discriminant analysis effect Size; MAPK: Mitogen-activated protein kinase; OUT: Operational taxonomic unit; PCA: Principal component analysis; PFS: Progression free survival; PICRUSt: Phylogenetic investigation of communities by reconstruction of unobserved States; PPS: Post progression survival; TCGA: The cancer genome atlas; WGBS: Whole-genome bisulfite sequencing.

Supplementary Material

Supplementary figures and tables.

Acknowledgements

BALBc ANXA1 deficient mice were a kind gift from Prof. Roderick Flower, QMUL, UK. We would like to thank Professor Amaury Cazenave Gassiot for his advice for the fatty acid synthase assay.

Funding

Funding for this study was provided by grants from the Ministry of Education in Singapore (MOE2017-T2-2-044) awarded to L.H.K.L.

Author contributions

L.H.K.L. and J.Z.C. conceived the study. L.H.K.L. and J.Z.C. designed the experiments. J.Z.C., K.S., H.L., Z.L., and W.T. carried out sample isolation and experiments, J.Z.C., L.H.K.L. K.S., R.A. L.Z. C.O. and J.M.C. analyzed the data and prepared the figures. J.Z.C., J.M.C. and L.H.K.L. drafted and edited the manuscript. All authors discussed the results and implications and commented on the manuscript at all stages.

Data availability statement

16sRNS sequencing datasets (PRJNA794259), RNAseq datasets (PRJNA796197) and WGBS datasets (SAMN25128104, SAMN25128105, SAMN25128106, SAMN25128107, SAMN25128108, SAMN25128109) were deposited into NCBI.

Availability of data and materials: All data needed to evaluate the conclusions in the paper are presented in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Ethical Approval and Consent to participate-

All animal work was approved by the NUS Institutional Animal Care and Use Committee (IACUC) and was in accordance with the National Advisory Committee for laboratory Animal Research (NACLAR) Guidelines (Guidelines on the Care and Use of Animals for Scientific Purposes).

Competing Interests

The authors have declared that no competing interest exists.

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021;71:7-33

2. Moreno-Smith M, Lutgendorf SK, Sood AK. Impact of stress on cancer metastasis. Future Oncol. 2010;6:1863-81

3. Dhabhar FS. Effects of stress on immune function: the good, the bad, and the beautiful. Immunol Res. 2014;58:193-210

4. Dhabhar FS. The short-term stress response-Mother nature's mechanism for enhancing protection and performance under conditions of threat, challenge, and opportunity. Front Neuroendocrinol. 2018;49:175-92

5. Antoni MH, Dhabhar FS. The impact of psychosocial stress and stress management on immune responses in patients with cancer. Cancer. 2019;125:1417-31

6. Zhang L, Pan J, Chen W, Jiang J, Huang J. Chronic stress-induced immune dysregulation in cancer: implications for initiation, progression, metastasis, and treatment. Am J Cancer Res. 2020;10:1294

7. Lutgendorf SK, Sood AK, Antoni MH. Host factors and cancer progression: biobehavioral signaling pathways and interventions. J Clin Oncol. 2010;28:4094

8. Zhi X, Li B, Li Z, Zhang J, Yu J, Zhang L. et al. Adrenergic modulation of AMPK-dependent autophagy by chronic stress enhances cell proliferation and survival in gastric cancer. Int J Oncol. 2019;54:1625-38

9. Lutgendorf SK, Andersen BL. Biobehavioral approaches to cancer progression and survival: Mechanisms and interventions. Am Psychol. 2015;70:186

10. Dhabhar FS. Effects of stress on immune function: Implications for immunoprotection and immunopathology. 2011.

11. Powell N, Tarr A, Sheridan JF. Psychosocial stress and inflammation in cancer. Brain, Behav, Immun. 2013;30:S41-S7

12. Steptoe A, Hamer M, Chida Y. The effects of acute psychological stress on circulating inflammatory factors in humans: a review and meta-analysis. Brain, Behav, Immun. 2007;21:901-12

13. Schirmer M, Smeekens SP, Vlamakis H, Jaeger M, Oosting M, Franzosa EA. et al. Linking the human gut microbiome to inflammatory cytokine production capacity. Cell. 2016;167:1125-36 e8

14. Rosean CMB, Rutkowski MR. The influence of the commensal microbiota on distal tumor-promoting inflammation. Semin Immunol: Elsevier. 2017 p. 62-73

15. Borgi M, Collacchi B, Ortona E, Cirulli F. Stress and coping in women with breast cancer:unravelling the mechanisms to improve resilience. Neurosci Biobehav Rev. 2020;119:406-21

16. Jurjus A, Eid A, Al Kattar S, Zeenny MN, Gerges-Geagea A, Haydar H. et al. Inflammatory bowel disease, colorectal cancer and type 2 diabetes mellitus: The links. BBA clinical. 2016;5:16-24

17. Bhatt AP, Redinbo MR, Bultman SJ. The role of the microbiome in cancer development and therapy. CA Cancer J Clin. 2017;67:326-44

18. Mikó E, Kovács T, Sebő É, Tóth J, Csonka T, Ujlaki G. et al. Microbiome-microbial metabolome-cancer cell interactions in breast cancer-familiar, but unexplored. Cells. 2019;8:293

19. Goedert JJ, Jones G, Hua X, Xu X, Yu G, Flores R. et al. Investigation of the association between the fecal microbiota and breast cancer in postmenopausal women: a population-based case-control pilot study. JNCI: Journal of the National Cancer Institute. 2015 107

20. Buchta Rosean C, Bostic RR, Ferey JCM, Feng T-Y, Azar FN, Tung KS. et al. Preexisting Commensal Dysbiosis Is a Host-Intrinsic Regulator of Tissue Inflammation and Tumor Cell Dissemination in Hormone Receptor-Positive Breast Cancer. Cancer Res. 2019;79:3662-75

21. Hinshelwood RA, Clark SJ. Breast cancer epigenetics: normal human mammary epithelial cells as a model system. J Mol Med. 2008;86:1315-28

22. Pasculli B, Barbano R, Parrella P. Epigenetics of breast cancer: Biology and clinical implication in the era of precision medicine. Semin Cancer Biol. 2018;51:22-35

23. El Helou R, Wicinski J, Guille A, Adélaïde J, Finetti P, Bertucci F. et al. a distinct DNA methylation signature defines breast cancer stem cells and predicts cancer outcome. Stem Cells. 2014;32:3031-6

24. Pal B, Bouras T, Shi W, Vaillant F, Sheridan JM, Fu N. et al. Global changes in the mammary epigenome are induced by hormonal cues and coordinated by Ezh2. Cell Rep. 2013;3:411-26

25. Maruyama R, Choudhury S, Kowalczyk A, Bessarabova M, Beresford-Smith B, Conway T. et al. Epigenetic regulation of cell type-specific expression patterns in the human mammary epithelium. PLoS Genet. 2011;7:e1001369

26. Keating ST, El-Osta A. Epigenetics and Metabolism. Circul Res. 2015;116:715-36

27. Grunt TW. Interacting cancer machineries: cell signaling, lipid metabolism, and epigenetics. Trends Endocrinol Metab. 2018;29:86-98

28. Ham J, Lee S, Lee H, Jeong D, Park S, Kim SJ. Genome-Wide Methylation Analysis Identifies NOX4 and KDM5A as Key Regulators in Inhibiting Breast Cancer Cell Proliferation by Ginsenoside Rg3. The American Journal of Chinese Medicine. 2018;46:1333-55

29. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7:1-11

30. Xia J, Psychogios N, Young N, Wishart DS. MetaboAnalyst: a web server for metabolomic data analysis and interpretation. Nucleic Acids Res. 2009;37:W652-W60

31. Etchegaray J-P, Mostoslavsky R. Interplay between Metabolism and Epigenetics: A Nuclear Adaptation to Environmental Changes. Mol Cell. 2016;62:695-711

32. Cao W, Lee H, Wu W, Zaman A, McCorkle S, Yan M. et al. Multi-faceted epigenetic dysregulation of gene expression promotes esophageal squamous cell carcinoma. Nature Communications. 2020;11:3675

33. Hansen KD, Timp W, Bravo HC, Sabunciyan S, Langmead B, McDonald OG. et al. Increased methylation variation in epigenetic domains across cancer types. Nat Genet. 2011;43:768-75

34. Sohn JH, Lee YK, Han JS, Jeon YG, Kim JI, Choe SS. et al. Perilipin 1 (Plin1) deficiency promotes inflammatory responses in lean adipose tissue through lipid dysregulation. J Biol Chem. 2018;293:13974-88

35. Goetz L, Jarvers I, Schleicher D, Mikan K, Brunner R, Kandsperger S. The role of the endogenous oxytocin system under psychosocial stress conditions in adolescents suffering from anxiety disorder: study protocol for a parallel group controlled trial. BMC Psychology. 2021;9:61

36. Vétizou M, Pitt JM, Daillère R, Lepage P, Waldschmitt N, Flament C. et al. Anticancer immunotherapy by CTLA-4 blockade relies on the gut microbiota. Science. 2015;350:1079-84

37. Iida N, Dzutsev A, Stewart CA, Smith L, Bouladoux N, Weingarten RA. et al. Commensal bacteria control cancer response to therapy by modulating the tumor microenvironment. Science. 2013;342:967-70

38. Ingman WV. The Gut Microbiome: A New Player in Breast Cancer Metastasis. Cancer Res. 2019;79:3539-41

39. Fernández MF, Reina-Pérez I, Astorga JM, Rodríguez-Carrillo A, Plaza-Díaz J, Fontana L. Breast cancer and its relationship with the microbiota. Int J Env Res Public Health. 2018;15:1747

40. Ley RE, Turnbaugh PJ, Klein S, Gordon JI. Human gut microbes associated with obesity. Nature. 2006;444:1022-3

41. Parida S, Sharma D. The power of small changes: Comprehensive analyses of microbial dysbiosis in breast cancer. Biochimica et Biophysica Acta (BBA)-Reviews on Cancer. 2019;1871:392-405

42. Conterno L, Fava F, Viola R, Tuohy KM. Obesity and the gut microbiota: does up-regulating colonic fermentation protect against obesity and metabolic disease? Genes Nutr. 2011;6:241-60

43. Ruan J-W, Statt S, Huang C-T, Tsai Y-T, Kuo C-C, Chan H-L. et al. Dual-specificity phosphatase 6 deficiency regulates gut microbiome and transcriptome response against diet-induced obesity in mice. Nature Microbiology. 2016;2:16220

44. Zhang H, DiBaise JK, Zuccolo A, Kudrna D, Braidotti M, Yu Y. et al. Human gut microbiota in obesity and after gastric bypass. Proceedings of the National Academy of Sciences. 2009;106:2365-70

45. Hill MJ, Goddard P, Williams REO. GUT BACTERIA AND ÆTIOLOGY OF CANCER OF THE BREAST. The Lancet. 1971;298:472-3

46. Murray WR, Blackwood A, Calman KC, MacKay C. Faecal bile acids and clostridia in patients with breast cancer. Br J Cancer. 1980;42:856-60

47. Mikó E, Vida A, Kovács T, Ujlaki G, Trencsényi G, Márton J. et al. Lithocholic acid, a bacterial metabolite reduces breast cancer cell proliferation and aggressiveness. Biochim Biophys Acta Bioenerg. 2018;1859:958-74

48. Ma J, Sun L, Liu Y, Ren H, Shen Y, Bi F. et al. Alter between gut bacteria and blood metabolites and the anti-tumor effects of Faecalibacterium prausnitzii in breast cancer. BMC Microbiol. 2020;20:82

49. Lopez-Siles M, Duncan SH, Garcia-Gil LJ, Martinez-Medina M. Faecalibacterium prausnitzii: from microbiology to diagnostics and prognostics. The ISME Journal. 2017;11:841-52

50. Diehm YF, Marstaller K, Seckler A-M, Berger MR, Zepp M, Gaida MM. et al. The collagenase of the bacterium Clostridium histolyticum does not favor metastasis of breast cancer. Breast Cancer. 2022

51. Peters BA, Hayes RB, Goparaju C, Reid C, Pass HI, Ahn J. The microbiome in lung cancer tissue and recurrence-free survival. Cancer Epidemiology and Prevention Biomarkers. 2019;28:731-40

52. Giallourou N, Urbaniak C, Puebla-Barragan S, Vorkas PA, Swann JR, Reid G. Characterizing the breast cancer lipidome and its interaction with the tissue microbiota. Communications Biology. 2021;4:1229

53. Urbaniak C, Cummins J, Brackstone M, Macklaim JM, Gloor GB, Baban CK. et al. Microbiota of human breast tissue. Appl Environ Microbiol. 2014;80:3007-14

54. Urbaniak C, Gloor GB, Brackstone M, Scott L, Tangney M, Reid G. The microbiota of breast tissue and its association with breast cancer. Appl Environ Microbiol. 2016;82:5039-48

55. Hieken TJ, Chen J, Hoskin TL, Walther-Antonio M, Johnson S, Ramaker S. et al. The microbiome of aseptically collected human breast tissue in benign and malignant disease. Sci Rep. 2016;6:1-10

56. Chan AA, Bashir M, Rivas MN, Duvall K, Sieling PA, Pieber TR. et al. Characterization of the microbiome of nipple aspirate fluid of breast cancer survivors. Sci Rep. 2016;6:1-11

57. Costantini L, Magno S, Albanese D, Donati C, Molinari R, Filippone A. et al. Characterization of human breast tissue microbiota from core needle biopsies through the analysis of multi hypervariable 16S-rRNA gene regions. Sci Rep. 2018;8:1-9

58. Xuan C, Shamonki JM, Chung A, DiNome ML, Chung M, Sieling PA. et al. Microbial dysbiosis is associated with human breast cancer. PLoS One. 2014;9:e83744

59. Garcia-Estevez L, Moreno-Bueno G. Updating the role of obesity and cholesterol in breast cancer. Breast Cancer Res. 2019;21:35

60. Nazih H, Bard JM. Cholesterol, Oxysterols and LXRs in Breast Cancer Pathophysiology. Int J Mol Sci. 2020 21

61. Danilo C, Frank PG. Cholesterol and breast cancer development. Curr Opin Pharm. 2012;12:677-82

62. Nelson ER, Chang C-y, McDonnell DP. Cholesterol and breast cancer pathophysiology. Trends Endocrinol Metab. 2014;25:649-55

63. Llaverias G, Danilo C, Mercier I, Daumer K, Capozza F, Williams TM. et al. Role of cholesterol in the development and progression of breast cancer. The American journal of pathology. 2011;178:402-12

64. Silva YP, Bernardi A, Frozza RL. The Role of Short-Chain Fatty Acids From Gut Microbiota in Gut-Brain Communication. Front Endocrinol (Lausanne). 2020 11

65. Parada Venegas D, De la Fuente MK, Landskron G, González MJ, Quera R, Dijkstra G. et al. Short Chain Fatty Acids (SCFAs)-Mediated Gut Epithelial and Immune Regulation and Its Relevance for Inflammatory Bowel Diseases. Front Immunol. 2019 10

66. Nogal A, Valdes AM, Menni C. The role of short-chain fatty acids in the interplay between gut microbiota and diet in cardio-metabolic health. Gut Microbes. 2021;13:1-24

67. Zhao L, Huang Y, Lu L, Yang W, Huang T, Lin Z. et al. Saturated long-chain fatty acid-producing bacteria contribute to enhanced colonic motility in rats. Microbiome. 2018;6:107

68. Mentoor I, Engelbrecht AM, Nell T. Fatty acids: Adiposity and breast cancer chemotherapy, a bad synergy? Prostaglandins, Leukotrienes and Essential Fatty Acids. 2019;140:18-33

69. Mikó E, Kovács T, Sebő É, Tóth J, Csonka T, Ujlaki G. et al. Microbiome-Microbial Metabolome-Cancer Cell Interactions in Breast Cancer-Familiar, but Unexplored. Cells. 2019 8

70. Pala V, Krogh V, Muti P, Chajès V, Riboli E, Micheli A. et al. Erythrocyte membrane fatty acids and subsequent breast cancer: a prospective Italian study. J Natl Cancer Inst. 2001;93:1088-95

71. Chajès V, Hultén K, Van Kappel AL, Winkvist A, Kaaks R, Hallmans G. et al. Fatty-acid composition in serum phospholipids and risk of breast cancer: An incident case-control study in Sweden. Int J Cancer. 1999;83:585-90

72. Chajès V, Assi N, Biessy C, Ferrari P, Rinaldi S, Slimani N. et al. A prospective evaluation of plasma phospholipid fatty acids and breast cancer risk in the EPIC study. Ann Oncol. 2017;28:2836-42

73. Shannon J, King IB, Moshofsky R, Lampe JW, Li Gao D, Ray RM. et al. Erythrocyte fatty acids and breast cancer risk: a case-control study in Shanghai, China. The American journal of clinical nutrition. 2007;85:1090-7

74. Chajès V, Thiébaut AC, Rotival M, Gauthier E, Maillard V, Boutron-Ruault M-C. et al. Association between serum trans-monounsaturated fatty acids and breast cancer risk in the E3N-EPIC Study. Am J Epidemiol. 2008;167:1312-20

75. Habib N, Wood C, Apostolov K, Barker W, Hershman M, Aslam M. et al. Stearic acid and carcinogenesis. Br J Cancer. 1987;56:455-8

76. Evans LM, Cowey SL, Siegal GP, Hardy RW. Stearate preferentially induces apoptosis in human breast cancer cells. Nutr Cancer. 2009;61:746-53

77. Saadatian-Elahi M, Norat T, Goudable J, Riboli E. Biomarkers of dietary fatty acid intake and the risk of breast cancer: A meta-analysis. Int J Cancer. 2004;111:584-91

78. Zhou Y, Wang T, Zhai S, Li W, Meng Q. Linoleic acid and breast cancer risk: a meta-analysis. Public Health Nutr. 2016;19:1457-63

79. Calado A, Neves PM, Santos T, Ravasco P. The Effect of Flaxseed in Breast Cancer: A Literature Review. Front Nutr. 2018;5:4

80. Jiang WG, Bryce RP, Horrobin DF, Mansel RE. gamma-Linolenic acid blocks cell cycle progression by regulating phosphorylation of p27kip1 and p57kip2 and their interactions with other cycle regulators in cancer cells. Int J Oncol. 1998;13:611-7

81. Watkins G, Martin TA, Bryce R, Mansel RE, Jiang WG. Gamma-Linolenic acid regulates the expression and secretion of SPARC in human cancer cells. Prostaglandins Leukot Essent Fatty Acids. 2005;72:273-8

82. Wiggins AK, Kharotia S, Mason JK, Thompson LU. α-Linolenic Acid Reduces Growth of Both Triple Negative and Luminal Breast Cancer Cells in High and Low Estrogen Environments. Nutr Cancer. 2015;67:1001-9

83. Wang SC, Sun HL, Hsu YH, Liu SH, Lii CK, Tsai CH. et al. α-Linolenic acid inhibits the migration of human triple-negative breast cancer cells by attenuating Twist1 expression and suppressing Twist1-mediated epithelial-mesenchymal transition. Biochem Pharmacol. 2020;180:114152

84. Thirunavukkarasan M, Wang C, Rao A, Hind T, Teo YR, Siddiquee AA-M. et al. Short-chain fatty acid receptors inhibit invasive phenotypes in breast cancer cells. PLoS One. 2017;12:e0186334

85. Mikó E, Vida A, Kovács T, Ujlaki G, Trencsényi G, Márton J. et al. Lithocholic acid, a bacterial metabolite reduces breast cancer cell proliferation and aggressiveness. Biochimica et Biophysica Acta (BBA)-Bioenergetics. 2018;1859:958-74

86. Luu TH, Bard J-M, Carbonnelle D, Chaillou C, Huvelin J-M, Bobin-Dubigeon C. et al. Lithocholic bile acid inhibits lipogenesis and induces apoptosis in breast cancer cells. Cell Oncol. 2018;41:13-24

87. Fuhrman BJ, Feigelson HS, Flores R, Gail MH, Xu X, Ravel J. et al. Associations of the fecal microbiome with urinary estrogens and estrogen metabolites in postmenopausal women. The Journal of Clinical Endocrinology & Metabolism. 2014;99:4632-40

88. Ridlon JM, Kang D-J, Hylemon PB. Bile salt biotransformations by human intestinal bacteria. J Lipid Res. 2006;47:241-59

89. Javitt NB, Budai K, Raju U, Levitz M, Miller D, Cahan A. Breast-gut connection: origin of chenodeoxycholic acid in breast cyst fluid. The Lancet. 1994;343:633-5

90. Goldberg AA, Beach A, Davies GF, Harkness TA, LeBlanc A, Titorenko VI. Lithocholic bile acid selectively kills neuroblastoma cells, while sparing normal neuronal cells. Oncotarget. 2011;2:761

91. Tang W, Putluri V, Ambati CR, Dorsey TH, Putluri N, Ambs S. Liver- and Microbiome-derived Bile Acids Accumulate in Human Breast Tumors and Inhibit Growth and Improve Patient Survival. Clin Cancer Res. 2019;25:5972-83

92. Seiler N. Catabolism of polyamines. Amino Acids. 2004;26:217-33

93. Kovács T, Mikó E, Vida A, Sebő É, Toth J, Csonka T. et al. Cadaverine, a metabolite of the microbiome, reduces breast cancer aggressiveness through trace amino acid receptors. Sci Rep. 2019;9:1-14

94. Vattai A, Akyol E, Kuhn C, Hofmann S, Heidegger H, von Koch F. et al. Increased trace amine-associated receptor 1 (TAAR1) expression is associated with a positive survival rate in patients with breast cancer. J Cancer Res Clin Oncol. 2017;143:1637-47

95. Albertsen M, Hugenholtz P, Skarshewski A, Nielsen KL, Tyson GW, Nielsen PH. Genome sequences of rare, uncultured bacteria obtained by differential coverage binning of multiple metagenomes. Nat Biotechnol. 2013;31:533-8

96. Gao X, Chang S, Liu S, Peng L, Xie J, Dong W. et al. Correlations between &#x3b1;-Linolenic Acid-Improved Multitissue Homeostasis and Gut Microbiota in Mice Fed a High-Fat Diet. mSystems. 2020;5:e00391-20

97. Lin H, An Y, Hao F, Wang Y, Tang H. Correlations of Fecal Metabonomic and Microbiomic Changes Induced by High-fat Diet in the Pre-Obesity State. Sci Rep. 2016;6:21618 -

98. Wang J, Tian S, Yu H, Wang J, Zhu W. Response of Colonic Mucosa-Associated Microbiota Composition, Mucosal Immune Homeostasis, and Barrier Function to Early Life Galactooligosaccharides Intervention in Suckling Piglets. J Agric Food Chem. 2019;67:578-88

99. Gu Y-Y, Han Y, Liang J-J, Cui Y-N, Zhang B, Zhang Y. et al. Sex-specific differences in the gut microbiota and fecal metabolites in an adolescent valproic acid-induced rat autism model. Frontiers in bioscience (Landmark edition). 2021;26:1585-98

100. Wang S, Zhang X, Li H, Ren Y, Geng Y, Lu Z. et al. Similarities and differences of oligo/poly-saccharides' impact on human fecal microbiota identified by in vitro fermentation. Appl Microbiol Biotechnol. 2021;105:7475-86

101. Li X, Guo R, Wu X, Liu X, Ai L, Sheng Y. et al. Dynamic digestion of tamarind seed polysaccharide: Indigestibility in gastrointestinal simulations and gut microbiota changes in vitro. Carbohydr Polym. 2020;239:116194

102. Li P, Yin Y-L, Li D, Kim SW, Wu G. Amino acids and immune function. Br J Nutr. 2007;98:237-52

103. Cunningham GA, McClenaghan NH, Flatt PR, Newsholme P. L-Alanine induces changes in metabolic and signal transduction gene expression in a clonal rat pancreatic beta-cell line and protects from pro-inflammatory cytokine-induced apoptosis. Clin Sci (Lond). 2005;109:447-55

104. Farvid MS, Barnett JB, Spence ND, Rosner BA, Holmes MD. Types of carbohydrate intake and breast cancer survival. Eur J Nutr. 2021

105. Jiang Y, Pan Y, Rhea PR, Tan L, Gagea M, Cohen L. et al. A Sucrose-Enriched Diet Promotes Tumorigenesis in Mammary Gland in Part through the 12-Lipoxygenase Pathway. Cancer Res. 2016;76:24-9

106. Guo H, Peng Y, Hu Z, Li Y, Xun G, Ou J. et al. Genome-wide copy number variation analysis in a Chinese autism spectrum disorder cohort. Sci Rep. 2017;7:44155

107. Yu J, Wu WK, Li X, He J, Li XX, Ng SS. et al. Novel recurrently mutated genes and a prognostic mutation signature in colorectal cancer. Gut. 2015;64:636-45

108. Ning S, Li H, Qiao K, Wang Q, Shen M, Kang Y. et al. Identification of long-term survival-associated gene in breast cancer. Aging (Albany NY). 2020;12:20332-49

109. Kothari C, Clemenceau A, Ouellette G, Ennour-Idrissi K, Michaud A, C-Gaudreault R. et al. TBC1D9: An Important Modulator of Tumorigenesis in Breast Cancer. Cancers (Basel). 2021 13

110. Foo SL, Yap G, Cui J, Lim LHK. Annexin-A1 - A Blessing or a Curse in Cancer? Trends Mol Med. 2019;25:315-27

111. Ang EZ-F, Nguyen HT, Sim H-L, Putti TC, Lim LHK. Annexin-1 Regulates Growth Arrest Induced by High Levels of Estrogen in MCF-7 Breast Cancer Cells. Mol Cancer Res. 2009;7:266-74

112. Moraes LA, Kar S, Foo SL, Gu T, Toh YQ, Ampomah PB. et al. Annexin-A1 enhances breast cancer growth and migration by promoting alternative macrophage polarization in the tumour microenvironment. Sci Rep. 2017;7:17925

113. Bist P, Leow SC, Phua QH, Shu S, Zhuang Q, Loh WT. et al. Annexin-1 interacts with NEMO and RIP1 to constitutively activate IKK complex and NF-κB: implication in breast cancer metastasis. Oncogene. 2011;30:3174-85

114. Araújo TG, Mota STS, Ferreira HSV, Ribeiro MA, Goulart LR, Vecchi L. Annexin A1 as a Regulator of Immune Response in Cancer. Cells. 2021;10:2245

115. McArthur S, Loiola RA, Maggioli E, Errede M, Virgintino D, Solito E. The restorative role of annexin A1 at the blood-brain barrier. Fluids and Barriers of the CNS. 2016;13:17

116. Grewal T, Enrich C, Rentero C, Buechler C. Annexins in Adipose Tissue: Novel Players in Obesity. Int J Mol Sci. 2019;20:3449

117. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS. et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12:1-18

118. Saito K, Koido S, Odamaki T, Kajihara M, Kato K, Horiuchi S. et al. Metagenomic analyses of the gut microbiota associated with colorectal adenoma. PLoS One. 2019;14:e0212406

119. Louca S, Parfrey LW, Doebeli M. Decoupling function and taxonomy in the global ocean microbiome. Science. 2016;353:1272-7

120. Langille MGI, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA. et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol. 2013;31:814-21

121. Janssens Y, Nielandt J, Bronselaer A, Debunne N, Verbeke F, Wynendaele E. et al. Disbiome database: linking the microbiome to disease. BMC Microbiol. 2018;18:50

122. Lu Y, Zou L, Su J, Tai ES, Whitton C, Van Dam RM. et al. Meat and seafood consumption in relation to plasma metabolic profiles in a Chinese population: a combined untargeted and targeted metabolomics study. Nutrients. 2017;9:683

123. Pang Z, Chong J, Zhou G, de Lima Morais DA, Chang L, Barrette M. et al. MetaboAnalyst 5.0: narrowing the gap between raw spectra and functional insights. Nucleic Acids Res. 2021;49:W388-W96

124. Xie Z, Bailey A, Kuleshov MV, Clarke DJB, Evangelista JE, Jenkins SL. et al. Gene Set Knowledge Discovery with Enrichr. Current Protocols. 2021;1:e90

125. Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P. et al. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 2010;38:W214-W20

126. Galili T, O'Callaghan A, Sidi J, Sievert C. heatmaply: an R package for creating interactive cluster heatmaps for online publishing. Bioinformatics. 2017;34:1600-2

127. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017;45:W98-W102

128. Faul F, Erdfelder E, Lang A-G, Buchner A. G*Power 3: A flexible statistical power analysis program for the social, behavioral, and biomedical sciences. Behav Res Methods. 2007;39:175-91

129. Basu S, Duren W, Evans CR, Burant CF, Michailidis G, Karnovsky A. Sparse network modeling and metscape-based visualization methods for the analysis of large-scale metabolomics data. Bioinformatics. 2017;33:1545-53

Author contact

Corresponding author: L.H.K.L. (linalimedu.sg) of the Immunology Translational Research Program #03-06J, 28 Medical Drive, Centre National University of Singapore, Singapore 117456.

Received 2021-11-1
Accepted 2022-3-28
Published 2022-5-9

Citation styles

©2024 Ivyspring International Publisher. Terms of use