Breast cancer in postmenopausal women is associated with an altered gut metagenome

1The human gut harbors thousands of bacterial species, together making up a population as large as 1013–14 microbes, which encode 150-fold more genes than the human genome. The gut microbiota is composed of a large number of anaerobic microorganisms, predominantly Bacteroidetes and Firmicutes, which are affected by a multitude of factors including host genetics, lifestyle, and environment. The microbiota plays important roles in maintaining an intestinal mucosal barrier, antagonizing the colonization of pathogenic microorganisms, and contributing to metabolism and immune homeostasis.

2The gut microbiota has been linked to various diseases, such as inflammatory bowel disease, obesity, diabetes, rheumatoid arthritis, atopic manifestations, liver cirrhosis, cardiovascular diseases, mental diseases, and colorectal cancer. The gut microbiota exert an influence on both local and systemic metabolism and immunity, and alterations of gut microbiota have been associated with extra-intestinal cancers including hepatocellular carcinoma, to which they may contribute by triggering chronic inflammation and altering microenvironment and metabolism.

3Increasing evidence suggests that microbe-host interactions have the potential to influence or serve as a biomarker of breast cancer pathogenesis. A comparison of 11 breast cancer patients and 7 healthy controls revealed differences in the gut microbiota, with Clostridia, Enterobacterium, Lactobacilli, Bacteroides, and Escherichia coli enriched in patients. A comparison of 48 postmenopausal breast cancer case patients and 48 healthy controls revealed an altered, less diverse gut microbiota in patients: Clostridiales, Clostridiaceae, Faecalibacterium, and Ruminococcaceae were enriched in patients, while Dorea and Lachnospiraceae were relatively less abundant in patients. Among controls, microbiota diversity correlated with total estrogen levels, suggesting that gut microbiota may be implicated in breast cancer by responding to or affecting estrogen metabolism.

4Those previous studies have provided useful insights into the potential response of gut microbiota in breast cancer, but they have not been able to comprehensively catalog the taxonomies of the microbiota because they have relied on only biochemical analysis or 16S rRNA gene sequencing. In addition, previous studies did not explore the functional capacity of the gut microbiota in patients with breast cancer, which could provide more mechanistic insights into the role of the gut microbiota in this disease. As a result, how the gut microbiota and their biochemical and metabolic products change in breast cancer is unclear.

5To address these questions, we used shotgun metagenomic analysis to compare the gut microbial community and its functional capabilities between breast cancer patients and healthy controls.

6The study was approved by the Ethics Committee of the Affiliated Tumor Hospital of Guangxi Medical University (Nanning, China). Fecal samples were collected from 18 premenopausal breast cancer patients, 25 premenopausal healthy controls, 44 postmenopausal breast cancer patients, and 46 postmenopausal healthy controls (Table 1). All patients with breast cancer were diagnosed by pathological examination at the Affiliated Tumor Hospital, and healthy controls were recruited from the Medical Examination Center of the First Affiliated Hospital of Guangxi Medical University. Controls were free of breast cancer at medical examination. None of the study subjects had diarrhea, diabetes, ulcerative colitis, Crohn’s disease, or other infectious diseases. No subjects took antibiotics, steroid hormones, Chinese herbal medicine (including oral, intramuscular, or intravenous injection), or probiotics such as yogurt during the 3 months before fecal sample collection. Breast cancer patients did not receive chemotherapy, radiation, or surgery prior to fecal sample collection.

7Fecal samples were freshly collected from individuals and transported to the laboratory on ice. Samples were stored at − 80 °C until extraction. Bacterial DNA was extracted from fecal samples using the QIAampDNA Stool MiniKit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions.

8All samples were sequenced on the Illumina HiSeq × 10 platform. A paired-end library was constructed with 350-bp inserts for each sample. Low-quality reads and reads mapping to human DNA were removed from the raw data. For taxonomic assignments, the high-quality reads from each sample were aligned against the integrated reference catalog of the human gut microbiome (IGC) by bowtie2 using the criterion of “identity > 90%,” genes from the existing reference gene catalog IGC inherited their original taxonomic annotation, and the relative abundance of a taxon was calculated from the relative abundance of its genes. During KO profiling, genes from IGC inherited their original KO annotation, and KO abundance was calculated by summing the relative abundance of genes annotated to the same KO.

9Genes in the catalog were aligned against proteins in the Virulence Factors of Pathogenic Bacteria database using BLAST (version 2.2.24) set to default parameters, except that -p blastx -e 1e-5 -F F -a 4 -m 8. We selected the matches with the highest-scoring annotated hit containing identity > 40% and high-scoring segment pair scoring > 60 bits. The relative abundance of a virulence factor was calculated by summing the relative abundance of genes annotated to a feature. Genes in the gene catalog were aligned against the proteins in the Pathogen-Host Interactions database using BLAST (version 2.2.24) set to default parameters except that -p blastx -e 1e-5 -F F -a 4 -m 8. We selected the matches with the highest-scoring annotated hits containing an identity > 40% and high-scoring segment pair scoring > 60 bits. The relative abundance of an interaction gene was calculated by summing the abundance of genes annotated to a feature.

10Based on the species profile, we calculated the within-sample (alpha) diversity to estimate gut microbiota richness and evenness based on the Shannon index and Chao1 index. High alpha diversity indicates high diversity of gut microbiota within a sample. Between-sample differences in microbial composition (beta diversity) were assessed in terms of the Jensen-Shannon divergence (JSD). The JSD was calculated by the following steps: (1) We first calculated JSD between each two samples within one group. (2) The mean of all JSD values between a sample and others within one group was computed (the mean value represented the similarity of the sample to others). We compared the mean JSD values to find if the beta diversity is different or not among the groups.

11Samples were clustered based on relative genus abundances using JSD distance and the “partitioning around medoids” (PAM) clustering algorithm. The Calinski-Harabasz (CH) index was used to calculate the optimal number of clusters. Principal component analysis was used to visualize the taxonomic drivers of clusters.

12Demographics were compared among groups using Student’s t test or the chi-square test in SPSS 16.0 software (IBM, Chicago, IL, USA). R software (version 3.3.2) was used to perform other analyses. The Wilcoxon rank sum test was used to identify significant differences in abundance of genes, genera, virulence factors, interaction genes, and KOs. Differentially enriched pathways and modules were identified according to their reporter score from the Z scores of significant KOs. A module with a reporter score of Z > 1.6 was defined as differentially enriched. P values were adjusted based on the false discovery rate (FDR) using the method of Benjamini and Hochberg. Permutational multivariate analysis of variance (PERMANOVA) using the “adonis” function in the R Vegan package was performed to assess effects of phenotype on gene/taxa profiles. The R package “ade4,” which involves instrumental principal component analysis, was used to visualize the taxonomic drivers of clusters during enterotyping. The “pheatmap” package (version 1.0.8) was used to generate heat maps, and the clustering method used in “pheatmap” function was “correlation.” Spearman’s rank correlation was used to find correlations of metagenomic features and clinical indices.

13A species-based classifier was trained using the random forest package in R. A tenfold cross-validation was performed on a random forest model using the relative species abundance profile. The minimum error was calculated using fivefold cross-validation with the “rfcv” function, and the minimum error plus the s.d. at that point was used as the cutoff. The optimal number of species was selected by cross-validation with one SE rule. The case probability was calculated using this set of species and a receiver operating characteristic (ROC) curve within the pROC package in R. The model was tested on the testing set, and the prediction error was determined. Differentially enriched genes were identified using the Wilcoxon rank test, and adjusted P values were estimated using the R package “q value” (version 2.2.2). All differentially enriched genes (q value < 0.05) were annotated to the butanoate metabolism pathways (using their original KO annotation which was inherited from the integrated reference catalog of the human gut microbiome database).

14A total of 133 stool samples were sequenced from premenopausal breast cancer patients (n = 18), premenopausal healthy controls (n = 25), postmenopausal breast cancer patients (n = 44), and postmenopausal healthy controls (n = 46). The premenopausal breast cancer patients and controls were similar for age, BMI, and ethnicity (P > 0.05, Table 1); the postmenopausal breast cancer patients and controls were similar for age, BMI, age at menopause, and ethnicity (P > 0.05, Table 1). A total of 965 million 150-bp paired-end reads were generated, with an average (s.d.) of 7.25 ± 1.13 million reads for each sample. After quality control, we obtained 902 million high-quality reads free of adaptor and human DNA contaminants, with an average (s.d.) of 6.78 ± 1.08 million reads per sample (Additional file 1: Table S1).

15To determine whether the sequencing adequately captured the gene diversity of the gut microbiota, rarefaction analysis was performed. The curves in all samples were near saturation, suggesting that the sequencing depth was sufficient to capture most gene diversity (Additional file 2: Figure S1).

16A similar number of species was detected in premenopausal breast cancer patients and premenopausal healthy controls (P = 0.767, Wilcoxon rank sum test, Fig. 1a). Based on the species profile, we calculated the within-sample (alpha) diversity to estimate gut microbiota richness and evenness based on the Shannon index and Chao1 index. The mean Chao1 index was similar between premenopausal breast cancer patients and premenopausal controls (P = 0.777, Wilcoxon rank sum test, Fig. 1b). The mean Shannon index was higher for premenopausal patients than premenopausal controls (P = 0.027, Wilcoxon rank sum test, Fig. 1c).

17In contrast, the between-sample variability (beta diversity) of the gut microbiota community structure tended to be lower in premenopausal patients than in premenopausal controls (P = 0.056, Fig. 1d).

18The number of species was significantly higher in postmenopausal breast cancer patients than in postmenopausal controls (P = 0.003, Wilcoxon rank sum test, Fig. 1a). Consistently, the mean Chao1 index was higher in postmenopausal patients than in postmenopausal controls (P = 0.007, Wilcoxon rank sum test, Fig. 1b). However, mean Shannon index was similar between postmenopausal breast cancer patients and postmenopausal controls (P = 0.502, Wilcoxon rank sum test, Fig. 1c). Beta diversity was higher for postmenopausal patients than postmenopausal controls (P < 0.001, Fig. 1d).

19Previous studies have suggested that the human gut microbiome can be assigned to several robust enterotypes. To group the breast cancer patients and control samples into enterotype clusters, we applied the PAM method using JSD for the relative abundance of genera.

20The optimal number of enterotypes was 2 as indicated by the CH index (Additional file 2: Figure S2a). Principal component analysis was used to cluster the samples of the four groups into two enterotypes (Additional file 2: Figure S2b). Enterotype 1 had a relatively high level of Bacteroides, enterotype 2, a relatively high level of Prevotella (Additional file 2: Figure S2c). These two enterotypes have been observed in European and Chinese populations. However, we found no significant relationship between enterotype and breast cancer disease status, either when we compared premenopausal patients with premenopausal controls (P = 0.141) or when we compared postmenopausal patients with postmenopausal controls (P = 0.445; Fisher’s exact test in both cases; Additional file 1: Table S2; Additional file 2: Figure S2d).

21To further explore features of the gut microbial community in breast cancer patients, we compared the relative abundances of species between patients and controls. The taxonomic assignment for the metagenomic data was carried out using bowtie. The relative abundance of gut microbiota was calculated by summing the abundance of genes. Relative abundance of the gut microbiota in the four groups at the species level is shown in Fig. 2.

22There was no significant difference in gut microbiota species between premenopausal breast cancer patients and premenopausal healthy controls (q value > 0.05, Wilcoxon rank sum test; Additional file 1: Table S3). In contrast, 45 species differed significantly between postmenopausal patients and postmenopausal controls: 38 species were enriched in patients, including Escherichia coli, Klebsiella sp_1_1_55, and Prevotella amnii, while 7 species were reduced in patients, including Porphyromonas uenonis, Eubacterium eligens, and Lactobacillus vaginalis (q value < 0.05; Table 2; Additional file 1: Table S4, Fig. 3). PERMANOVA analysis showed that breast cancer status, menopause status, and age were significant factors for explaining the variation in the examined gut microbial samples (P < 0.05; Additional file 1: Table S5).

23To illustrate the potential diagnostic value of gut microbiota for breast cancer in postmenopausal women, we used a random forest classifier in an attempt to detect breast cancer samples from among a mixture of samples from postmenopausal patients and healthy controls. Tenfold cross-validation was repeated for five times with a training set consisting of 44 postmenopausal patients and 46 postmenopausal controls; 14 optimal species markers were selected, including Escherichia coli, Shigella sp_D9, Eubacterium eligens, Proteus mirabilis, and Fusobacterium varium (Additional file 2: Figure S6, Additional file 1: Table S6). ROC curves for the training set showed a remarkable performance in the training set when discriminating between postmenopausal breast cancer patients and postmenopausal healthy controls by specificity and sensitivity; the area under receiver operating curve (AUC) was 85.52; and 95% confidence interval (CI) is 77.57–93.47% (Fig. 4b). Next, we tested the same markers for their ability to detect breast cancer among 43 samples not used during training, comprising 18 premenopausal breast cancer patients and 25 premenopausal healthy controls. The AUC was 72% (95% CI 56.01–88.44%; Fig. 4d).

24To analyze proteins encoded by genes in gut microbiota, genes in the catalog were aligned against proteins in the Pathogen-Host Interactions (PHI) database and in the Virulence Factors of Pathogenic Bacteria database. The relative abundances of genes encoding virulence factors or pathogen-host interaction were calculated by summing the abundances of genes annotated to a feature.

25PHI genes were relatively more abundant in postmenopausal breast cancer patients than in postmenopausal controls (P = 0.021, Fig. 5a). The top 15 representation of Pathogen-Host Interactions genes implicated in several human diseases, such as urinary tract infections, other infections, and tuberculosis (q value < 0.05, Wilcoxon rank sum test, Fig. 5b; Additional file 1: Table S7).

26Virulence factors were relatively more abundant in postmenopausal patients than in controls (P = 0.016, Fig. 5c); these factors included LOS glycosyltransferase, peritrichous flagella, and type I fimbriae (q value < 0.05, Wilcoxon rank sum test, Fig. 5d; Additional file 1: Table S8).

27Next, we asked whether the 45 microbial species differing between postmenopausal breast cancer patients and postmenopausal controls correlated with well-established clinical indices of breast cancer based on Spearman correlation analysis (Fig. 6). Two gut microbiota species, Acinetobacter radioresistens (Spearman rho = 0.413, P = 0.015, q value > 0.05) and Vibrio cholerae (Spearman rho = 0.349, P = 0.043, q value > 0.05), positively but weakly correlated with C4 levels, whereas Yersinia enterocolitica (Spearman rho = − 0.345, P = 0.046, q value > 0.05) negatively but weakly correlated with C4 levels. Two species, Acinetobacter radioresistens (Spearman rho = 0.442, P = 0.009, q value > 0.05) and Enterococcus gallinarum (Spearman rho = 0.386, P = 0.024, q value > 0.05), positively but weakly correlated with levels of high-sensitivity C-reactive protein (hsCRP). Anaerococcus vaginalis (Spearman rho = 0.48, P = 0.002, q value > 0.05) and Porphyromonas uenonis (Spearman rho = 0.42, P = 0.009, q value > 0.05) positively but weakly correlated with CD19 levels, while Enterococcus gallinarum (Spearman rho = − 0.351, P = 0.031, q value > 0.05) negatively but weakly correlated with CD19 levels. Shewanella putrefaciens (Spearman rho = 0.379, P = 0.025, q value > 0.05) and Erwinia amylovora (Spearman rho = 0.351, P = 0.039, q value > 0.05) positively but weakly correlated with estradiol levels. Actinomyces sp. HPA0247 (Spearman rho = − 0.384, P = 0.017, q value > 0.05) negatively but weakly correlated with CD3+CD8+ T cell numbers. All these gut flora species were enriched in postmenopausal breast cancer patients relative to controls. Conversely, Eubacterium eligens was enriched in controls relative to patients, yet it negatively but weakly correlated with CD3+CD4+ T cell numbers (Spearman rho = − 0.349, P = 0.032, q value > 0.05) and IgA levels (Spearman rho = − 0.532, P = 0.001, q value > 0.05).

28We explored functional features of the gut microbiota across the four groups in our study by annotating the gene catalog based on the KEGG modules. Modules differing between breast cancer patients and healthy controls with a reporter score > 1.6 were identified.

29Among premenopausal women, 44 KEGG modules were significantly different between patients and controls (q value < 0.05, Wilcoxon rank sum test, Fig. 7a; Additional file 1: Table S9). Modules enriched in patients included the PTS system, secretion system, vitamin B12 transport system, amino acid transport system, and manganese/iron transport system. Modules enriched in controls included aminoacyl-tRNA biosynthesis, coenzyme A biosynthesis, nucleotide synthesis, and dicarboxylate-hydroxybutyrate cycle.

30Among postmenopausal women, 43 KEGG modules were enriched in patients, including lipopolysaccharide biosynthesis, iron complex transport system, vitamin B12 transport system, PTS system, secretion system, amino acid transport system, and beta-oxidation (q value < 0.05, Wilcoxon rank sum test, Fig. 7b; Additional file 1: Table S10).

31The genes annotated to butanoate metabolism pathways differentially enriched in gut microbiome of postmenopausal patients and controls. Fourteen butyrate-synthesis genes were found: 10 genes were enriched in controls and 4 genes were enriched in postmenopausal patients (q value < 0.05, Wilcoxon rank sum test; Additional file 1: Table S11).

32Here, we performed a comprehensive metagenomic comparison of gut microbiota in breast cancer patients and healthy controls. Microbiota were analyzed in terms of taxonomic profile, genetic functional capacity, and associations with clinical indices of breast cancer. Our results identify various compositional and functional features of the gut microbiota metagenome that differ between postmenopausal patients and healthy controls, suggesting that they may be associated with postmenopausal breast cancer.

33Significant taxonomic differences in gut microbiota were not detected between premenopausal breast cancer patients and controls. In contrast, several bacterial species were found to be enriched in postmenopausal patients relative to controls: Escherichia coli, Citrobacter koseri, Acinetobacter radioresistens, Enterococcus gallinarum, Shewanella putrefaciens, Erwinia amylovora, Actinomyces sp. HPA0247, Salmonella enterica, and Fusobacterium nucleatum. These results are consistent with previous research on several of these species, which has suggested associations with breast cancer. In addition, we found a weak positive correlation of Shewanella putrefaciens and Erwinia amylovora with estradiol (P < 0.05), although there is no statistical significance after correction for multiple testing (q value > 0.05), which may be related to the small sample, but the association can still be considered exploratory. These results are consistent with the idea that the gut microbiota can influence or be affected by estrogen metabolism and thereby provide an independent biomarker of breast cancer. Recent literature has demonstrated that the gut microbiota is the modulation of systemic estrogens. Elevated levels of circulating estrogens are associated with an increased risk of breast cancer.

34We found that Eubacterium eligens and Roseburia inulinivorans were less abundant in postmenopausal breast cancer patients than in postmenopausal controls. Roseburia inulinivorans produces butyrate. In order to explore a potential association between butyrate-producing bacteria and breast cancer, we had added a differentially enriched gene analysis to show the potential association, all differentially enriched genes were annotated to the butanoate metabolism pathways, and finally, 14 butyrate-synthesis genes were found: 10 genes were enriched in controls and 4 genes were enriched in postmenopausal patients. Notably, these butyrate-synthesis genes were reduced in postmenopausal patients, which may be related to the decrease in butyrate-producing bacteria.

35Butyrate acts as an anti-inflammatory agent, mainly by inhibiting the activation of nuclear factor κB (NF-κB) in intestinal epithelial cells. Butyrate can also act on immune cells via specific G-protein-coupled receptors expressed on immune cells. Reductions in colonic butyrate can promote inflammation. These findings provide further evidence for the idea that alterations in the gut microbial community are associated with breast cancer. For example, a decrease in Roseburia inulinivorans may render postmenopausal women more prone to inflammation and therefore at higher risk of breast cancer if the decrease in Roseburia had occurred at the time breast cancer was initiated.

36Dysbiosis was detected in the gut microbiomes of postmenopausal breast cancer patients, but it was not detected in premenopausal patients. Therefore, the dysbiosis observed in the gut microbiomes of breast cancer patients may depend on age and menopause status.

37Breast cancer-associated alterations in the gut microbial community likely translate into alterations in gut microbial functions. Among premenopausal women, breast cancer was associated in our study with enrichment in genes involved in the PTS system, secretion system, vitamin B12 transport system, and manganese/iron transport system. Among postmenopausal women, breast cancer was associated with enrichment in genes involved in lipopolysaccharide (LPS) biosynthesis, iron complex transport system, vitamin B12 transport system, PTS system, and secretion system. Iron enrichment affects the gut microbiome, increases pathogen abundance, and induces intestinal inflammation. The PTS and secretion systems are associated with diabetes, liver cirrhosis, and rheumatoid arthritis, while vitamin B12 status correlates positively with breast cancer risk in women. Lipopolysaccharide is a potent trigger of macrophage-mediated systemic inflammation, which has been suggested to play an important role in promoting the transformation of inflammation into tumorigenesis. Enrichment of the iron transport system and lipopolysaccharide biosynthesis in gut microbiota may cause systemic low-grade inflammation, thereby increasing the risk of breast cancer if the dysbiosis observed in the patient cohort was present in the same cohort prior to contracting breast cancer.

38In conclusion, we have found alterations of gut microbial community and functions in postmenopausal breast cancer patients. The gut microbiota may regulate or respond to host immunity and metabolic balance. In this way, our study suggests an association between gut microbiota and development of postmenopausal breast cancer. However, our data do not allow us to determine whether the altered gut metagenome is the consequence of the disease process or is somehow involved in its pathogenesis.