The intestinal microbiota influences the microenvironment of metastatic colon cancer by targeting miRNAs.

This study aimed to investigate the molecular mechanisms through which the intestinal microbiota and microRNAs (miRNAs) participate in colon cancer metastasis. Intestinal flora data, and the GSE29621 (messenger RNA/long non-coding RNA [mRNA/lncRNA]) and GSE29622 (miRNA) datasets, were downloaded from The Cancer Gene Atlas and Gene Expression Omnibus databases, respectively. Immune-related cells in M1 vs. M0 samples were analyzed using the Wilcoxon test. Furthermore, an lncRNA-miRNA-mRNA (competing endogenous RNA [ceRNA]) network was constructed, and survival analysis of RNAs in the network was performed. A total of 16 miRNA-genus co-expression pairs containing eight microbial genera and 15 miRNAs were screened; notably, Porphyromonas and Bifidobacterium spp. were found to be associated with most miRNAs, and has-miR-3943 was targeted by most microbial genera. Furthermore, five immune cell types, including activated natural killer cells, M1 macrophages, resting mast cells, activated mast cells and neutrophils, were differentially accumulated between the M1 and M0 groups. Enrichment analysis suggested that mRNAs related to colon cancer metastasis were mainly involved in pathways related to bacterial and immune responses. Survival analysis revealed that TMEM176A and PALM3 in the ceRNA network were significantly associated with the prognosis of patients with colon cancer. In conclusion, this study revealed a potential mechanism by which the intestinal microbiota influences the colon cancer microenvironment by targeting miRNAs.

Introduction

Colon cancer, a common malignancy of the digestive tract, is currently the fourth most commonly diagnosed cancer type and is the second leading cause of cancer-related death among 36 cancer types worldwide (Bray et al. 2018). Colon cancer is a highly heterogeneous disease and a relevant public health issue in a growing number of countries (Siegel et al. 2021). Currently, the standards for clinical treatment and prognostic prediction of survival and recurrence of colon cancer rely mainly on the tumor-node-metastasis system and on the histopathological criteria established by the American Joint Committee on Cancer (Brierley). In general, early colonoscopy screening and treatment can improve patient outcomes; however, there are no obvious symptoms in the early stages of colorectal cancer, and approximately 15% to 25% of patients present with synchronous metastasis at diagnosis (Bonnot and Passot 2019). Moreover, the prognosis of patients with metastatic colon cancer appears to be much worse than that of patients in the early and intermediate stages (Tjandra and Chan 2007). Therefore, it is important to explore the potential molecular mechanisms underlying colon cancer development.

MicroRNAs (miRNAs), small single-stranded non-coding RNA molecules of 18–24 nucleotides in length, play key roles in the modulation of gene expression at the post-transcriptional level; moreover, miRNAs exert important biological effects in animals, participating in processes such as immune system development, immune response (Xiao and Rajewsky 2009) and metabolism (Vienberg et al. 2017). Accumulating evidence has shown that intestinal miRNAs are key factors for the maintenance of a healthy gastrointestinal environment. For example, it has been suggested that intestinal microorganisms and microRNAs may interact to regulate the expression of host genes (Dalmasso et al. 2011, Moein et al. 2019). Moreover, miRNAs exert a wide range of effects on the intestinal immune system and play an important role in the pathogenesis of intestinal diseases (Kalla et al. 2015). Additionally, an increasing number of studies have confirmed the functional role of miRNAs in mediating the communication between intestinal microorganisms and host intestinal epithelial cells (Takeda et al. 2011, Aguilar et al. 2019). However, although these reports have suggested a role of the intestinal microbiota and miRNAs in the development of colon cancer, it is not yet known whether the intestinal microbiota alters the tumor microenvironment by influencing miRNA expression or stability. Therefore, in order to explore the molecular mechanisms through which the intestinal microbiota and miRNAs influence colon cancer progression, this study examined intestinal microflora data from a large cohort derived from The Cancer Genome Atlas, as well as miRNA-seq and long non-coding RNA/messenger RNA (lncRNA/mRNA)-seq datasets. Subsequent correlation analysis, as well as the construction and examination of a lncRNA-miRNA-mRNA regulatory axis (Supplementary Fig. 1), revealed alterations in the microflora and miRNA expression during the formation of colon cancer metastasis.

Methods

Data retrieval and processing

Intestinal flora data derived from colon cancer tissue samples at different M stages, including 89 M0-stage tissue samples and 13 M1-stage tissue samples, were downloaded from The Cancer Microbiome Atlas (Dohlman et al. 2021) (https://tcma.pratt.duke.edu/) database in The Cancer Genome Atlas. In addition, RNA-seq data (log2 (FPKM + 1)), miRNA-seq data (log2 (RPM + 1)), clinical information (including age, sex, tumor-node-metastasis stage, tumor stage and site of tumor occurrence) and survival information (overall survival [OS] and OS time) of the corresponding colon cancer samples were obtained from the UCSC-Xena platform (Goldman et al. 2018) (https://toil.xenahubs.net). After being intersected with the intestinal flora samples, RNA-seq data from a total of 100 samples, including 87 M0-stage tissue samples and 13 M1-stage tissues samples, were eventually included in subsequent analyses; among these, survival information (OS and OS time) was available for 94 tissue samples. Furthermore, 94 miRNA-seq samples were included in the follow-up study, consisting of 81 M0-stage tissues and 13 M1-stage tissues.

Additionally, the GSE29621 (mRNA/lncRNA) and GSE29622 (miRNA) datasets were obtained from NCBI Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) (Barrett et al. 2005) for external validation. The two datasets were obtained from the same samples, namely, 65 colon cancer samples with survival information (OS and OS time), of which 46 samples were in the M0 stage and 18 samples were in the M1 stage.

For The Cancer Genome Atlas intestinal flora data, preprocessed data with deviation removal were directly downloaded from the database, and colon cancer tissue samples with M stage were extracted for subsequent analysis. For The Cancer Genome Atlas RNA-seq data, log2 (FPKM + 1) values were downloaded directly, and Ensemble gene IDs were converted to gene symbols according to the annotation information (hg38, gencode.v22.annotation.gene.probemap of the GENCODE database) (Harrow et al. 2012) (https://www.gencodegenes.org/). For GSE29621 (mRNA/lncRNA) and GSE29622 (miRNA) data, the processed and standardized probe expression matrices were downloaded directly, and the corresponding platform annotation file was downloaded to perform gene symbol transformation of the probes. For different probes corresponding to the same gene symbol, the average value was taken as the gene expression value for subsequent analysis.

Analysis of differentially accumulated microflora and miRNAs associated with colon cancer metastasis

First, an inter-group t-test in R 3.6.1 was used to compare the relative abundance of each microbial taxon between the M1 and M0 groups. Next, the ggplot package (version 3.2.1) in R 3.6.1 was used to draw a bar chart showing the relative abundance of each microbial taxon in the M1 and M0 groups. Moreover, differential miRNA expression analysis was performed for M1 vs. M0 samples using the classical Bayesian method provided in the limma package (version 3.34.7) in R 3.6.1 (Smyth 2013) (https://bioconductor.org/packages/release/bioc/html/limma.html).

As reported in a previous study (Bi et al. 2020), intestinal miRNAs from intestinal epithelial cells or externally derived through the diet interact with intestinal microorganisms and regulate the composition and distribution of intestinal microbial communities. Hence, to screen for closely related miRNAs and intestinal microorganisms that are present during the formation of colon cancer metastasis, the relative abundance and expression values of metastasis-related microbial taxa and miRNAs, respectively, were extracted from each sample. Moreover, the Pearson correlation coefficient R was calculated to obtain genus-miRNA co-expression pairs, with a threshold of |R| > 0.2 and P < 0.05.

Identification of metastasis-related immune cells and mRNAs/lncRNAs in the cancer microenvironment

The CIBERSORT deconvolution algorithm (Chen et al. 2018) was used to estimate the abundance of infiltrating immune cells in all samples, applying the LM22 dataset provided in the CIBERSORT website as the characteristic gene expression template. After calculating the abundance of 22 infiltrating cell types in each sample, differences in the abundance of immune-related cells in M1 vs. M0 samples were analyzed using the Wilcoxon test (P < 0.05). Violin plots were drawn using the R package vioplot (version 0.3.2).

To further identify the mRNAs and lncRNAs closely related to colon cancer metastasis, weighted gene co-expression network analysis (WGCNA) was conducted. Briefly, based on the top 3000 genes showing the greatest variation among samples, the R package WGCNA (Langfelder and Horvath 2008) (version 1.61, https://cran.r-project.org/web/packages/WGCNA/) was used to identify gene set modules with highly synergistic variation. Subsequently, an inter-group t-test was performed on the genes in the modules, and lncRNAs and mRNAs whose expression differed significantly between M1 and M0 samples were further selected as metastasis-related RNAs for subsequent studies (P < 0.05). Furthermore, Gene Ontology analysis of biological processes, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Enrichment Analysis of Reactome Gene Sets were conducted on metastasis-related mRNAs using Metascape (Zhou et al. 2019) (http://metascape.org) (parameters: minimum overlap = 3; P value cutoff = 0.01; minimum enrichment = 1.5).

Construction of the competing endogenous RNA network

Metastasis-related miRNAs and lncRNAs were selected to construct competing endogenous RNA (ceRNA) networks. miRNA-mRNA interactions were predicted using miRWalk 2.0 (Dweep and Gretz 2015) (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) and the miRWalk, miRANDA, miRDB, PITA, RNA22 and TargetScan databases. To better understand the potential functions and importance of the mRNAs in miRNA-mRNA pairs, KEGG pathway enrichment analysis was performed. Pearson correlation coefficients between miRNA-mRNA pairs were calculated, and pairs characterized by R values greater than 0 were regarded as interacting pairs. lncBaseV2 (Paraskevopoulou et al. 2016) (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex) was used to retrieve interacting lncRNA-miRNA pairs. Finally, Cytoscape software (version 3.4.0, http://chianti.ucsd.edu/cytoscape-3.4.0/) (Shannon et al. 2003) was used to build and visualize an interactive ceRNA network based on the predicted lncRNA-miRNA and miRNA-mRNA interactions. The R package ggalluvial (version 0.11.1) was used to draw a Sankey map.

Validation and Kaplan–Meier survival analysis of the selected RNAs

The expression of the molecules in the ceRNA network was validated in the Gene Expression Omnibus database and visualized using ggplot (version 3.2.1). To further assess whether the molecules in the ceRNA network were significantly associated with colon cancer prognosis, the optimal cutoff point was determined using the R package Survminer (version 0.4.3) according to the expression value, survival time and survival state of each RNA molecule. Finally, the R package survival (version 2.42–6) was used for survival analysis and to carry out a log-rank survival test, with a significance threshold of P < 0.05.

Results

Differentially accumulated intestinal microorganisms and miRNAs associated with colon cancer metastasis

Differences in the relative abundance of each microbial taxon between M1 and M0 tissue samples were analyzed using the inter-group t-test in R 3.6.1 (Table 1). The results demonstrated that the abundance of Coriobacteriia and Deltaproteobacteria at the class level; Desulfovibrionales, Bifidobacteriales and Coriobacteriales at the order level; Odoribacteraceae, Peptostreptococcaceae, Desulfovibrionaceae, Bifidobacteriaceae, Coriobacteriaceae and Porphyromonadaceae at the family level; and Odoribacter, Barnesiella, Tyzzerella, Bilophila, Klebsiella, Porphyromonas, Bifidobacterium, Collinsella, Dialister and Peptostreptococcus at the genus level differed significantly between the groups (P < 0.05).

Table 1.

Differential microbial community between M1 and M0 tissue samples at different levels.

Phylum	Class	Order	Family	Genus
Bacteroidetes	Bacteroidia	Bacteroidales	Odoribacteraceae*	Odoribacter*
			Barnesiellaceae	Barnesiella*
			Porphyromonadaceae*	Porphyromonas*
Proteobacteria	Deltaproteobacteria*	Desulfovibrionales*	Desulfovibrionaceae*	Bilophila*
	Gammaproteobacteria	Enterobacterales	Enterobacteriaceae	Klebsiella*
Actinobacteria	Actinobacteria	Bifidobacteriales*	Bifidobacteriaceae*	Bifidobacterium*
	Coriobacteriia*	Coriobacteriales*	Coriobacteriaceae*	Collinsella*
Firmicutes	Clostridia	Eubacteriales	Clostridiaceae	Tyzzerella*
			Peptostreptococcaceae*	Peptostreptococcus*
	Negativicutes	Veillonellales	Veillonellaceae	Dialister*

Note: * represents P < 0.05.

Additionally, 58 miRNAs differentially expressed between M1 and M0 samples were identified. The observation of a volcano plot confirmed that the expression levels of these miRNAs were significantly different between the two groups (Fig. 1A). Subsequently, correlation analysis of 10 distinct genera and 58 unique miRNAs was conducted. As illustrated in Fig. 1B, 16 miRNA-genus co-expression pairs containing eight microbial genera and 15 miRNAs were obtained. Within these co-expression pairs, Porphyromonas and Bifidobacterium spp. were associated with most miRNAs. Moreover, has-miR-3943 was targeted by both Porphyromonas and Dialister spp.

Figure 1.

Differential expression analysis of miRNAs. Volcano plot of differentially expressed miRNAs (A). Correlation network of genera and miRNAs (B). Gray circles represent bacteria with lower abundance at the genus level; hexagons represent miRNAs, upregulated in orange, downregulated in blue; red lines indicate positive correlations; and blue lines indicate negative correlations..

Immunoinfiltrating cells associated with colon cancer metastasis

Differences in the relative abundance of 22 infiltrating immune cell types between the M1 and M0 groups were analyzed using the Wilcoxon test. The findings revealed that five immune cell types, namely, activated natural killer cells, M1 macrophages, resting mast cells, activated mast cells and neutrophils, showed differential accumulation between the two groups (Fig. 2). Among these cell types, the relative abundance of infiltrating activated natural killer cells, M1 macrophages and resting mast cells in the M0 group was significantly higher than that in the M1 group, while the relative abundance of infiltrating activated mast cells and neutrophils was lower.

Figure 2.

Violin plot of the relative abundance of 22 infiltrating immune cell types. M0 and M1 samples are depicted in green and red, respectively.

Identification of mRNAs and lncRNAs associated with colon cancer metastasis

WGCNA analysis of the top 3000 genes yielded 10 distinct gene modules (Fig. 3A). Next, the correlation between the feature vector genes of each module and clinical phenotypes (including OS, OS time, tumor-node-metastasis stage, tumor stage and the abundance of the five immune cell populations that showed differential accumulation in the above screening) was estimated (Fig. 3B). The blue module (544 genes), which was significantly correlated with M stage expression, was selected as the key module, and was observed to be significantly correlated also with OS and the abundance of activated natural killer cells and M1 macrophages.

Figure 3.

mRNAs associated with colon cancer metastasis. WGCNA module clustering (A). Left panel, module clustering results. The vertical axis represents the coefficient of difference, and the blue line indicates a coefficient of difference of 0.3. Right panel, systematic clustering tree of genes and gene modules generated by dynamic cutting. Different colors represent different gene modules, while genes that could not be incorporated into any other module are included in the gray module. Correlation analysis between WGCNA modules and traits (B). Enrichment analysis of transfer-related mRNAs for functional pathways (C). Colors ranging from light to dark indicate decreasing P values, and strip length is inversely proportional to the Pvalue, with longer strips indicating more significantly enriched pathways.

To further identify lncRNAs and mRNAs significantly correlated to the M stage in this module, inter-group differential expression analysis of the genes in the module was conducted. The results showed that the expression of 122 mRNAs and 11 lncRNAs differed significantly between the M1 and M0 groups. Furthermore, enrichment analysis suggested that these mRNAs were mainly involved in 355 Gene Ontology biological processes, 67 KEGG pathways and 24 Reactome Gene Sets. More importantly, 30 functional clusters were identified based on their genetic similarity. As shown in Fig. 3C, pathways connected to bacterial and immune responses were enriched, suggesting that the bacterial community and immunity might play a key role in metastasis formation.

ceRNA network construction

Target gene prediction for the 15 miRNAs found to be closely related to the intestinal flora during the formation of colon cancer metastasis yielded a total of 63 miRNA-mRNA interacting pairs, including 12 miRNAs and 33 mRNAs (Fig. 4A). KEGG pathway enrichment analysis showed that only eight miRNA target genes were involved in KEGG pathways (Fig. 4B). Furthermore, co-expression analysis of the lncRNAs and mRNAs related to colon cancer metastasis revealed 445 lncRNA-mRNA co-expression pairs, including 11 lncRNAs and 122 mRNAs (Fig. 4C). Finally, based on the above lncRNA-miRNA, miRNA-mRNA and lncRNA-mRNA interacting pairs, a ceRNA network and a Sankey diagram were constructed, which contained three upregulated lncRNAs, six downregulated miRNAs and 17 upregulated mRNAs (Fig. 4D and E).

Figure 4.

ceRNA analysis. miRNA-target gene interaction network (A). Red circles represent upregulated mRNAs, green circles represent downregulated mRNAs, orange hexagons represent upregulated miRNAs, blue hexagons represent downregulated miRNAs, arrows indicate the target genes regulated by the miRNAs and node size represents the degree of connectivity. KEGG pathway enrichment analysis of miRNA target genes (B). The larger the nodes in the figure, the greater the ratio of the number of enriched genes to the total number of miRNA target genes, and colors ranging from blue to red indicate greater and greater P values. lncRNA-mRNA co-expression network (C). Red circles represent upregulated mRNAs, green circles represent downregulated mRNAs, light purple diamonds represent upregulated lncRNAs, dark purple diamonds represent downregulated lncRNAs and gray dotted lines represent co-expression relationships between mRNAs and lncRNAs. ceRNA network (D). Red circles represent upregulated mRNAs, light purple diamonds represent upregulated lncRNAs, blue hexagons represent downregulated miRNAs, green dotted lines represent the co-expression relationships between lncRNAs and mRNAs and gray arrows represent miRNA regulation. Blue T-shaped lines represent lncRNAs competing with miRNAs. Sankey diagram (E).

Survival analysis of the RNAs in the ceRNA network was performed. Analysis of the Kaplan–Meier curves revealed that TMEM176A and PALM3 were significantly associated with prognosis. Specifically, patients exhibiting high TMEM176A or PALM3 expression displayed shorter survival times than those of patients showing low expression of these genes (Fig. 5A).

Figure 5.

Validation analysis. Kaplan–Meier survival curve analysis of ceRNA molecules (A). Differential expression of the RNA molecules in the ceRNA network (B).

In addition, the expression profiles of the RNA molecules in the M0 and M1 groups were verified using Gene Expression Omnibus datasets. The results revealed that most of the mRNAs tended to be upregulated in M1 samples, consistent with the findings in the analyzed dataset. In particular, the levels of GGT7, PRDX5, PTP4A3 and TMEM176B in the M1 group were higher than those in the M0 group (P < 0.05; Fig. 5B). Moreover, the lncRNAs SATB2-AS1 and OSER1-AS1 were upregulated in the M1 group compared with the M0 group (P < 0.05).

Discussion

The survival of patients with tumor metastasis significantly affects the organization of individualized treatments. In colon cancer, patients with metastasis have worse survival outcomes than those without metastasis, with a 5-year survival rate of only 14.0% (Provenzale et al. 2018). Interactions between miRNAs from intestinal epithelial cells and intestinal microbes are critical for maintaining intestinal health and ameliorating gastrointestinal diseases such as colon cancer (Bi et al. 2020). Hence, the aim of this study was to investigate whether the gut microbiome alters the tumor microenvironment by influencing the expression of miRNAs that may promote cancer metastasis.

Our findings revealed a number of differentially accumulated genera and 58 differentially expressed miRNAs between M1 and M0 tissues, and a total of 16 miRNA-genus co-expression pairs containing eight microbial genera and 15 miRNAs were screened. Within these co-expression pairs, Porphyromonas and Bifidobacterium spp. were found to be associated with most miRNAs, and has-miR-3943 was targeted by both Porphyromonas and Dialister spp. A previous study suggested that oral administration of Porphyromonas gingivalis, altering the gut microbiome and the serum metabolome, is associated with impaired gut barrier function, resulting in endotoxemia and subsequent inflammation of the liver and adipose tissue (Kato et al. 2018). Moreover, Wang and colleagues (2021) have reported that P. gingivalis can promote colorectal carcinoma by activating the hematopoietic inflammasome. miR-3943 is involved in the development of resected gastric cancer and can be used as an independent prognostic biomarker (Woo et al. 2021); however, the role of miR-3943 in colon cancer has not yet been reported. Our findings revealed that has-miR-3943 is associated with colon cancer metastasis and targeted by Porphyromonas spp., suggesting that Porphyromonas spp. in the intestinal tract participate in colon cancer metastasis through the regulation of has-miR-3943 expression.

Antitumor immune memory is essential for long-term prevention of tumor recurrence and metastasis. Accumulating evidence suggests that the tumor immune microenvironment plays a crucial role in modulating antitumor immunity and is associated with tumor progression (Xia et al. 2021). Moreover, immune cells can respond rapidly to changes in the tumor microenvironment in different diseases (Wu et al. 2020, Liu and Li 2021). The data of this study indicated that five immune cell types, namely, activated natural killer cells, M1 macrophages, resting mast cells, activated mast cells and neutrophils, were differentially accumulated between the M1 and M0 groups. In addition, this study identified 122 mRNAs and 11 lncRNAs related to colon cancer metastasis, and enrichment analysis suggested that these mRNAs were mainly involved in pathways linked to the bacterial and immune responses, suggesting that the bacterial community and immunity might play a key role in the formation of metastasis.

In addition, based on the RNAs related to colon cancer metastasis, a ceRNA network containing three upregulated lncRNAs, six downregulated miRNAs and 17 upregulated mRNAs was constructed. Moreover, survival analysis of the RNAs in the ceRNA network showed that TMEM176A and PALM3 were significantly associated with the prognosis of patients with colon cancer. Transmembrane protein 176A (TMEM176A) is located on the human chromosome region 7q36.1, which often displays loss of heterozygosity (Kimmel et al. 2006). Notably, abnormal expression of TMEM176A has been reported to be related to cancer pathology; this also points at TMEM176A as a promising potential therapeutic target for the treatment of some cancer types (Cuajungco et al. 2012). In particular, Gao and colleagues (2017) have reported that TMEM176A methylation is involved in the progression of human colon cancer and can be used as an independent prognostic marker. Paralemmin 3 (PALM3), belonging to the PALM protein family, was first described in Xenopus laevis as Xlgv7/Xlcaax-1 (Cornish et al. 1992), and can act as an adaptor connecting intrinsic membrane proteins to each other, the cytoskeleton, or motor proteins (Hu et al. 2005, Chen et al. 2011). Previous studies have demonstrated that downregulation of PALM3 can improve the survival rate, ameliorate the severity of lung injury and inhibit the production of pro-inflammatory cytokines in rats (Chen et al. 2017); however, its role in colon cancer has not been investigated until now.

In conclusion, through the analysis in M1 vs. M0 samples, this study identified several intestinal microbial genera, such as Porphyromonas and Bifidobacterium, that are associated with miRNAs such as has-miR-3943 in metastatic colon cancer. Moreover, the mRNAs associated with colon cancer metastasis are mainly involved in pathways connected to bacterial and immune responses. Furthermore, TMEM176A and PALM3 in the ceRNA network were found to be significantly associated with the prognosis of patients with colon cancer. Taken together, the findings of this study reveal a potential mechanism by which the intestinal microbiota influences the microenvironment of colon cancers by targeting miRNAs.

Data availability

The data used to support the findings of this study are available from the corresponding author upon request.

Click here for additional data file.