To obtain complete datasets of exosome miRNA (exo-miRNA) related gene expression in CRC, four public databases were retrieved using the keywords “colon cancer,” “rectal cancer,” “colorectal cancer,” “exosome,” “exosomes,” and “exosomal.” The databases Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) (26), ArrayExpress (https://www.ebi.ac.uk/arrayexpress), ExpressionAtlas (https://www.ebi.ac.uk/gxa), and Genomic Expression Archive (GEA, https://www.ddbj.nig.ac.jp/gea) (Supplementary Table 1). Afterward, the following criteria were used to filter the data: (i) the data objects were miRNA expression profiles in blood or cellular exosomes. (ii) The blood exosome samples contained clinical staging information. (iii) At least three sample replicates. (iv) They were completely free and downloadable. Finally, three CRC exosome datasets were obtained: one associated with serum exosomes (GSE39833), which included 20 patients with stage I, 20 patients with stage II, 20 patients with stage IIIA, 16 patients with stage IIIB, 12 patients with stage IV, and 11 healthy individuals; the other two were associated with CRC cell line exosomes (GSE39832 and GSE39814), which included four fetal bovine serum-treated controls, three normal colon cells (FHC) and 15 CRC cells (HCT116, HT29, RKO, SW48, SW480; three of each).

Background-corrected and normalized Agilent microarray expression matrixes from the datasets above were downloaded, negative values were removed, and interpolation was performed using the caret R package's preProcess function (27). All values were log2-transformed. The limma R package (28) was used to analyze differentially expressed miRNAs in serum exosomes from different stages of CRC. Differential significance was determined with p < 0.05 and |logFC| > 1 (Supplementary Figure 1), and differential intersection screening and definition were performed by Jvenn online tool (29). MiRBase 22.1 (http://www.mirbase.org/index.shtml) was used to transform the names of all miRNAs obtained. Heatmap data were normalized using the scale function, clustered using fully linked hierarchical clustering based on Euclidean distances (if necessary), and visualized by the pheatmap R package (30).

To determine exo-miRNA target gene expression and miRNA expression in tissues, the GDC-TCGA (31) colon adenocarcinoma (COAD) and rectal adenocarcinoma (READ) log2(count+1)-transformed mRNA expression matrix, log2(RPM+1)-transformed miRNA expression matrix and the corresponding clinical data were downloaded from UCSC-ZENA (https://xenabrowser.net/datapages/) (32). COAD included 41 normal samples and 471 samples from colon cancer patients and READ included 10 normal samples and 167 samples from rectal cancer patients. Moreover, dataset GSE115513 (33) was used to validate the miRNA expression in the tissue level, which included 649 normal intestinal mucosa samples and 731 CRC tumor tissue samples (including 231 in stage I, 181 in stage II, 219 in stage III, and 100 in stage IV). The normalized expression matrix was downloaded from GEO and log2 transformed. Key miRNAs were selected and compared among different stages, the Wilcoxon test was used to evaluate the differences between two stages and the Kruskal-Wallis test was applied to estimate the differences among all stages. Gene expression and drug sensitivity data of 22 CRC cell lines were downloaded from the CCLE database (34).

To obtain the targeted genes from differentially expressed exo-miRNAs, 11 databases were searched based on the multiMiR R package (35) and miRWalk 2.0 (36), including three experimental validation databases [miRecords 4.0 (37), mirTarBase 7.0 (38) and Tarbase 8.0 (39)] and eight prediction databases [DIANA-microT-CDS 5.0 (40), EIMMo 5.0 (41), MicroCosm 5.0 (http://www.ebi.ac.uk/enright-srv/microcosm/cgi-bin/targets/v5/download.pl), miRanda (42), miRDB 6.0 (43), PicTar 2.0 (44), PITA 6.0 (45) and TargetScan 7.2 (46)]. The top 20% of data with high confidence based on p-value and included in at least 4 (5,127 miRNA-gene pairs) of the 8 (133,738 miRNA-gene pairs) prediction databases were taken. miRanda, PITA 6.0, and TargetScan 7.2 were specifically chosen to include only predictions with conserved loci. Furthermore, data from at least 2 (2,825 miRNA-gene pairs) of the 3 (46,455 miRNA-gene pairs) experimental validation databases were chosen. Finally, 7,306 miRNA-gene pairs were obtained from the 11 databases. The accuracy of prediction and the richness of miRNA targeting prediction results were ensured by combining experimental validation databases and prediction databases.

The VarElect online tool was used to perform phenotypic correlation analysis based on the keywords “Colon Cancer” OR “Rectal Cancer” OR “Colorectal Cancer” to further screen for target genes associated with CRC (47). Phenotype scores greater than five, and only direct phenotypic associations were considered.

The Pipeline of Outlier miRNA Analysis (POMA) is an analytical process developed by Chen et al. (24) for identifying key miRNAs. A translational POMA (tPOMA) in exo-miRNA was created to identify key exo-miRNAs in CRC (Figure 2). tPOMA is constructed based on two key assumptions: (i) The number of targeted CRC-related genes/TFs can reflect the regulatory activity of CRC exo-miRNAs. (ii) If only one exo-miRNA target a CRC-related gene/TF, that exo-miRNA is more likely to show regulatory specificity to that gene/TF. Thus, the greater the number of CRC-related genes/TFs regulated by an exo-miRNA in a single line, the more important that miRNA is likely to be. Furthermore, we suppose many CRC exo-miRNAs regulate a CRC-related gene. In that case, the gene is more likely to be repressed (miRNAs typically repress gene expression), and such a gene was referred to as a hotspot gene regulated by exo-miRNAs.

Based on this assumption, a tPOMA-based screening was performed in four steps on all obtained SC exo-miRNAs: (i) First, a global targeting network of SC exo-miRNAs was built using the previously mentioned 11 databases. (ii) The miRNA single-line regulatory network was further constructed, in which genes only regulated by one exo-miRNA (other exo-miRNAs do not regulate the gene). (iii) To generate a hotspot gene network, single-line regulated genes were removed from the global miRNA-gene network (genes are regulated by at least two exo-miRNAs). (iv) Using the VarElect tool, CRC phenotype-related genes were selected to built the CRC-related single-line regulated gene network and hotspot gene network (47). (v) The TF was confirmed using the THTF database (http://humantfs.ccbr.utoronto.ca) (48). Finally, the network of exo-miRNAs targeting genes was mapped using Cytoscape (3.7.1) (49).

The online tool miRNet (50) developed by Chang et al. was used to construct a target network of SC exo-miRNAs (target molecules including mRNAs, lncRNAs, circRNAs, pseudogenes, and sncRNAs) and ranked the miRNAs by the degree to demonstrate the reliability of our tPOMA model. The top 10 miRNAs were compared to the top 10 miRNAs predicted by our tPOMA model (Supplementary Figure 2), and there is a high degree of overlap (6 out of 10 miRNAs were consistent).

The following three steps were used to build the SS exo-miRNA subnetworks: (i) The global network of SS exo-miRNAs and target genes was built using the results of the previously mentioned 11 databases. (ii) Genes only present at a specific stage were extracted by Jvenn. (iii) SS exo-miRNAs and target gene sub-networks were further constructed by Cytoscape.

To better understand the role of these targeted genes in the development of CRC, enrichment analysis of Gene Ontology (GO) (51) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (52) was performed by the clusterProfiler R package (53), and further ranking enrichment analysis was performed using the g: Profiler online tool (http://biit.cs.ut.ee/gprofiler/gost) (54). According to the importance of targets, the top 50 genes with the highest relevance to CRC were extracted (FDR < 0.05) to identify the more important biological entries and relevant pathways. The GOplot R package (55) was used to plot and visualize enrichment string maps.

To obtain exo-miRNAs with high AUC, the pROC R package (56) was used for ROC analysis (the ROC curves of all SS and SC exo-miRNAs in Supplementary Figure 3). For combined-multivariate ROC analysis, logistic regression models were first developed using the glmnet R package (57), and further assessed by pROC.

Furthermore, to determine whether these exo-miRNAs and their target hotspot genes affect CRC patients' overall survival (OS), the survival R package (58) was used to conduct a survival analysis of key SC and SS exo-miRNAs. The TCGA survival data included both COAD and READ data, with means calculated for high and low expression groups, and log-rank p < 0.05 was considered a significant survival impact. GEPIA (59) was used to assess targeted hotspot genes for expression and survival analysis [367 CRC samples and 667 normal samples from TCGA (31) and GTEx (60) were included, |logFC| > 0.5, p < 0.05, survival analysis grouped by median expression].

Twenty-two cell lines were grouped according to their median expression in the CCLE database (34). The sensitivity of two grouped cell lines to the EGFR-targeting drug lapatinib (as measured by IC50) was compared using a two-tailed t-test for each targeted EGFR tyrosine kinase inhibitor resistance pathway genes.

The following four steps were used to identify the key CRC exo-miRNA-associated ceRNA networks: (i) The lncRNAs positively associated with miRNA-targeted mRNA were extracted from cbioportal online tool (61) (640 COAD samples from TCGA; Spearman analysis, p < 0.05). (ii) The experimentally validated miRNA-lncRNA pairs were determined in the lncBase v2 database (62). (iii) The corresponding ceRNA was identified when lncRNAs and mRNAs linked to the same miRNA. (iv) The corresponding networks were built using Cytoscape (49). Targeted pathway analysis of 10 key miRNAs was accomplished using Cytoscape's ClueGO plugin (63), with relevant KEGG pathways selected at p < 0.05.

Differential expression analysis of serum exo-miRNAs in different stages of CRC patients yielded 89 SC exo-miRNAs and 94 SS exo-miRNAs (p < 0.05, |logFC| > 1) (Figure 3A). Using hierarchical clustering of heatmaps, 89 SC exo-miRNAs were discovered to effectively distinguish CRC patients from healthy individuals (Figure 3B). The expression of SS exo-miRNAs was also investigated (Figure 3C), with stage I and IIIB being poorly differentiated, but stages II, IIIA, and IV being better differentiated, with stage IV showing the most significant difference.

An exosomal miRNA-gene targeting network for the 89 SC exo-miRNAs was constructed (Figure 4A), which included 89 exo-miRNAs and their 7,306 predicted target genes from 11 databases (see Materials and methods section for detailed steps). Second, by establishing the tPOMA model (detailed steps are included in the Materials and methods section), a single-line regulatory network of exo-miRNAs was obtained in CRC (Figure 4B). Then, only genes with a high degree of relevance to CRC were selected from the overall miRNA-gene network and the miRNA single-line regulatory network (for detailed steps, see Materials and methods section). Concurrently, single-line regulated genes were excluded from the overall exosomal miRNA-gene network. Finally, a hotspot gene network was constructed including 56 exo-miRNAs and targeted 173 genes associated with CRC (at least two exo-miRNAs regulated each gene) (Figure 4C; Supplementary Table 2), as well as a miRNA single-line regulatory network constructed including 43 exo-miRNAs and targeted 186 genes associated with CRC (a single exo-miRNA regulated each gene) (Figure 4D; Supplementary Table 3).

To confirm that screened exo-miRNA-targeted genes are indeed related to the biological functions of exosomes and to compare the biological functions of exo-miRNA-regulated genes with varying targeting burdens, the corresponding enrichment analysis of biological functions was performed.

First, GO and KEGG enrichment analyses were performed for 173 hotspot genes and 186 single-line genes (Supplementary Figure 4A), respectively. It was found that these genes were enriched in signaling and regulation (e.g., signaling receptor binding, identical protein binding, enzyme binding, growth factor binding, and kinase activity); cell membrane complex structures (e.g., membrane-bounded organelle, membrane region, intracellular organelle, cell surface and receptor complex); regulation of substance secretion, transport and cellular response (e.g., positive regulation of metabolic process, positive regulation of cellular process, cellular response to chemical stimulus, positive regulation of biological process and cell population proliferation) and various biological pathways associated with cancer (e.g., colorectal cancer, gastric cancer, breast cancer, PI3K-AKT signaling pathway and proteoglycans in cancer). The findings show that these genes are important for secretory transport and exosome regulation.

Second, two classes of genes was ranked based on CRC-related scores calculated by VarElect, and their biological functions were further explored using ranking enrichment analysis. Some differences in GO entries between hotspot genes (Figure 5A; Supplementary Table 4) and single-line regulated genes (Figure 5B; Supplementary Table 4) were discovered, as the former is more related to the regulation of cell death and phosphorylation metabolism in biological processes, with molecular functions focused on the binding of enzymes, proteins, and transcription factors, with the majority of these molecules localized in organelle membranes. The latter is primarily involved in activating enzymes in biological processes, molecular functions in cell development regulation and are primarily localized on the cell surface and membrane-related components. Both types of molecules are enriched in several cancer-associated pathways, including colorectal cancer, EGFR tyrosine kinase inhibitor resistance, endocrine resistance, PI3K/AKT signaling pathway, and cancer proteoglycans. However, they were also enriched in various pathways that are important for cancer development, immune checkpoint, drug tolerance, and pluripotent differentiation of stem cells (e.g., P53 signaling pathway, HIF-1 signaling pathway, FOXO signaling pathway, specific to hotspot genes; PD-L1 expression and PD-1 checkpoint pathway in cancer, platinum drug resistance, signaling pathways regulating pluripotency of stem cells, specific to single-line regulated genes) (Figures 5C,D; Supplementary Table 4).

Finally, the expression of genes from cancer-associated pathways was investigated. The hotspot genes being more down-regulated in these pathways (Figure 5E), which may be related to their negative regulation by multiple exo-miRNAs. In contrast, genes regulated by a single exo-miRNA were almost identical in terms of up-and down-regulation in these pathways (Figure 5F).

Through traditional and further ranking enrichment analysis, distinct regulatory patterns for two distinct types of genes with varying targeting burdens were identified in colorectal carcinogenesis and development, paving the way for future research into more precise regulation of exo-miRNAs.

According to the tPOMA model construction hypothesis, miRNAs targeting more CRC-related genes or transcription factors (TF) in the single-line regulatory network of exo-miRNAs may be more important. Additionally, in the hotspot gene network, CRC-related genes/TF targeted by a greater number of exo-miRNAs may be more susceptible to repressed. The corresponding exo-miRNAs and hotspot genes were arranged in descending order according to NCGSM-TF and INDEGREE (Table 1). To narrow down the screening of key miRNAs/genes and identify molecules of greater importance. NCGSM-TF and NCGSM were combined to thoroughly examine the top 15 miRNAs and hotspot genes with INDEGREE ≥ 5 (26 genes). Five of these miRNAs (miR-23b-3p, miR-186-5p, miR-22-3p, miR-142-3p and miR-221-3p) had significantly lower expression (Figure 6A) and high AUC (AUC > 85%, Figure 6B) in serum exosomes of CRC patients compared with healthy individuals, while miR-146a-5p demonstrated higher expression.

Interestingly, most miRNAs in CRC patients' serum exosomes were down-regulated (data not shown), which is consistent with the general down-regulation of original miRNA levels in cancer patients (64) and suggests that most miRNAs in exosomes may act as cancer suppressors. Ten of 26 hotspot genes screened were upregulated in tissue expression in CRC patients, and seven were downregulated (|logFC| > 0.5, p < 0.05), and nine were below threshold (Supplementary Figure 5). Furthermore, five genes (CCNA2, MAPK1, PTPRD, MET, and CDKN1A) were significantly associated with CRC patients' better survival (Figure 6C). There is reason to believe that these genes with a high miRNA suppressive burden in CRC patients act as protective molecules in the compensatory increase after tumor occurrence and positively correlate with patient survival. Finally, the exo-miRNA regulatory network of five survival-related genes was constructed (Figure 6D). Three miRNAs were identified (miR-146a-5p, miR-22-3p, and miR-23b-3p) with the strong single-line regulatory ability and the ability to regulate hotspot genes in conjunction with other exo-miRNAs.

First, target genes and screened CRC-related genes were identified for all SS exo-miRNAs obtained based on 11 databases, and an SS miRNA-target network was constructed (Figure 7A). Subsequently, the CRC-related genes of each stage were identified by Venn diagram (Figure 7B), and finally, five SS CRC miRNA-target sub-networks were obtained (stage IV, Figure 7C; stage I, Figure 7D; stage II, Figure 7E; stage IIIA, Figure 7F; stage IIIB, Figure 7G). Additionally, GO and KEGG enrichment analyses were performed to explore the role of different SS subnetworks in CRC development.

According to the results of GO analysis, stage I focuses on cellular development and metabolic processes (e.g., cellular development process, cell differentiation, and positive regulation of metabolic process) and is found on the cell surface, endomembrane system, vesicle, and intracellular vesicle (Supplementary Figure 4B). Stages II and IIIA have increased signal transduction and response (e.g., regulation of cell communication, signaling, signal transduction, and chemical response). They are associated with transcriptional regulation and cellular junctions (e.g., transferase complex, chromosome, DNA polymerase complex, and cell-cell junction). The main focus of stage IV is on metabolic and phosphorylation-related processes (e.g., positive regulation of the cellular metabolic process, phosphorus metabolic process, and phosphorylation regulation), which are closely related to cell adhesion and transcriptional regulation (e.g., cell junction, anchoring junction, focal adhesion, and transcription regulator complex). All stages are linked to molecular binding and enzyme activation in molecular function. Furthermore, genes in different SS subnetworks showed some pathway similarities and differences.

The effects of exo-miRNAs-regulated genes on CRC progression and changes in biological processes were observed by exploring the biological functions of genes in different subnetworks, which aids in understanding the regulation of exo-miRNAs in different stages of CRC patients.

The expression and ROC curve diagnostic performance of all SS exo-miRNAs were evaluated. There were two I-specific miRNAs (miR-301a-3p, miR-548i), two IIIA-specific miRNAs (miR-23a-3p, miR-532-5p), and six IV-specific miRNAs (miR-1246, miR-194-3p, miR-33a-3p, miR-485-3p, miR-194-5p and miR-379-5p) present high Stage-AUC (Stage-AUC > 65%, it represent the performance of miRNA can distinguish one stage from other stages).

Their stage diagnostic performance (Stage-AUC) was 65–75%, and their overall diagnostic performance (the miRNA can distinguish CRC patients from healthy people) was 62–85% (Figure 8; Supplementary Table 5). Furthermore, seven miRNAs found in previously constructed SS subnetworks were selected as key SS exo-miRNAs, with two for stage I (miR-301a-3p and miR-548i), one for stage IIIA (miR-23a-3p), and four for stage IV (miR-194-3p, 33a-3p, miR-485-3p and miR-194-5p). Given the low stage differentiation properties of these SS exo-miRNAs, a combined diagnostic analysis of two I-specific exo-miRNAs and four IV-specific exo-miRNAs was performed, with the stage I combined diagnostic model achieving 72% AUC (Figure 9A) and the stage IV combined diagnostic model achieving 93% AUC (Figure 9B).

A survival analysis based on CRC data in TCGA was performed to further understand the impact of these SC/SS exo-miRNAs on the survival of CRC patients. It was discovered that low expression of mir-1246 was associated with worse overall survival (OS) of patients (Figure 9C), while low expression of mir-379 and mir-485 predicted worse survival after 3 years (Figure 9D). Additionally, the expression of key exo-miRNAs in colon cancer tissues and healthy tissues were explored in TCGA (Figure 9E), in which the miR-23b-3p, miR-146a-5p, miR-301a-3p, and miR-194-3p present the same expression changes compared with CRC serum exosomes. The miRNA expression was further validated in GSE115513, and most of the miRNA expression alterations were similar to that of the TCGA cohort, except for miR-194-5p (SS-IV) had the opposite expression (Supplementary Figure 6). Among them, both tissue and exosome expression of miR-194-3p (SS-IV) were significantly decreased in stage IV, suggesting its possible role as a potential tumor suppressor miRNA for exosomal transport. Finally, the expression of key exo-miRNAs in the exosomes of five different CRC cell lines was investigated (Supplementary Figure 7). The expression of exo-miRNAs in CRC cell line exosomes was largely consistent with that in serum exosomes (e.g., miR-22-3p, miR-23b-3p, miR-301a-3p, miR-194-3p, miR-485-3p). In summary, miR-23b-3p, miR-301a-3p and miR-194-3p were validated being the most stably expressed stage-associated miRNAs in CRC serum exosomes, cell exosomes and tissues (Table 2).

The EGFR tyrosine kinase inhibitor tolerance pathway was enriched in hotspot genes and single-line regulated genes. The response of exo-miRNAs regulatory genes in this pathway to the EGFR-targeted drug lapatinib was analyzed to further identify exo-miRNAs regulatory genes with tolerance response to EGFR tyrosine kinase inhibitors (of the three EGFR targeting drugs in the CCLE database, only lapatinib responded to these genes). Low PIK3R1 expression among the hotspot genes in the CRC Cell significantly increased resistance to lapatinib (Figure 9F), and the previously identified key SC exo-miRNAs: miR-23b-3p and miR-221-3p were able to regulate this gene directly. In contrast, only high EGFR expression among single-line regulated genes in the CRC cell increased resistance to lapatinib (Figure 9G), and EGFR was directly regulated by miR-146a-5p.

Finally, the potential pathways, ROC diagnostic capabilities, survival, and drug response of 10 key stage-associated exo-miRNAs were summarized to demonstrate the potential importance of each miRNA (Figure 9H). MiR-22-3p, miR-23a-3p and miR-23b-3p are the top3 important stage-associated exo-miRNAs, and they display close correlation with several vital cancer-associated pathways such as TGF-beta signaling, p53 signaling, and hippo signaling pathways.

In addition to targeting the corresponding mRNAs, inhibiting their expression, and influencing downstream functions, miRNAs can also interact with the corresponding lncRNAs, and a competitive relationship between lncRNAs and mRNAs that bind to the same miRNA, known as competing endogenous RNA (ceRNA) regulation, can form. To explore the potential ceRNA regulation of stage-associated exo-miRNAs, the ceRNA network was constructed (Figure 10A; Supplementary Table 6). NEAT1 and KCNQ1OT1 are lncRNAs associated with all three ceRNA regulatory modules in the SC exo-miRNA ceRNA network. NEAT1 and KCNQ1OT1 have been reported to be upregulated in a variety of cancers and promote tumorigenesis by altering the expression levels of sponged miRNAs (65, 66). Additionally, the corresponding ceRNA regulatory networks in stages I, IIIA, and IV were identified, which have corresponding roles in regulating the various stages. To better understand the impact of these miRNAs on downstream biological pathways, their miRNA targeting pathways were analyzed using the ClueGO plugin in Cytoscape (63). There were 59 significantly enriched pathways (Supplementary Table 7), among which Toll-like receptor signaling pathway is the most significant (adj. p = 2.60e−07, adjusted by Bonferroni step down) (Figure 10B). Additionally, AKT3 was identified as the gene that is most strongly associated with all pathways and is involved in a variety of cancer-related processes. AKT3 is a member of the protein encoding gene of AKT serine-Threonine protein kinase. AKT is usually overexpressed in CRC and serves as a key target of targeted therapy based on PI3K/AKT signaling pathway (67, 68).