MicroRNA expression profiles identify biomarkers for predicting the response to chemoradiotherapy in rectal cancer.

Neoadjuvant chemoradiotherapy (nCRT) following surgery significantly improves the survival rate of patients with rectal cancer. However, nCRT is associated with significant adverse symptoms and high medical costs. Therefore, it is important to investigate potential biomarkers for the prediction of the response to nCRT in patients with rectal cancer. The present study identified candidate biomarkers for predicting a complete response (CR) to nCRT in patients with rectal cancer and investigated the associated mechanisms. Microarray data (accession no. GSE29298) was downloaded from the Gene Expression Omnibus database. Differentially expressed microRNAs (miRNAs/miR) were screened between the pathological CR (pCR) group and no pCR (incomplete response) group. miRNA target genes were predicted using the miRWalk 2.0 online tool and subjected to Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis. Furthermore, a miRNA co‑regulatory network was constructed and disease‑associated genes were predicted. The results demonstrated that a total of 36 upregulated and 5 downregulated miRNAs were identified between the two groups. Among these differentially expressed miRNAs, miR‑548c‑5p, miR‑548d‑5p and miR‑663a were significantly associated with a CR to nCRT. The co‑regulatory network and pathway analysis indicated that miR‑548c‑5p and miR‑548d‑5p may function together through stem cell pluripotency and ubiquitin‑mediated proteolysis signaling pathways. Furthermore, the prediction of disease‑associated genes demonstrated that miR‑548c‑5p/miR‑548d‑5p and miR‑663a may regulate genes associated with rectal cancer, including mutated in colorectal cancers (MCC) and adenomatous polyposis coli (APC), and colorectal neoplasms, including interleukin‑6 signal transducer (IL6ST), cell cycle checkpoint kinase 2 (CHEK2), marker of proliferation Ki‑67 (MKI67), cadherin 7 (CDH7), calreticulin (CALR) and transforming growth factor β1 (TGFB1). Therefore, miR‑548c‑5p, miR‑548d‑5p and miR‑663a are promising candidate biomarkers for predicting a CR to nCRT. miR‑548c‑5p/miR‑548d‑5p may be associated with a CR by regulating IL6ST, CHEK2, MKI67 and MCC. In addition, it may function through the pluripotency of stem cells and ubiquitin‑mediated proteolysis signaling pathways. miR‑663a may be associated with a CR to nCRT by targeting CDH7, CALR, APC and TGFβ1. Thus, the miRNA biomarkers investigated in the present study may represent novel therapeutic targets for the prediction and eventual improvement of the response to nCRT in patients with rectal cancer.

Introduction

Colorectal cancer (CRC) is the third most common type of malignant tumor in the digestive tract and poses a substantial threat to human health worldwide (1). In addition, CRC is the third highest contributor to cancer-associated mortality worldwide, with ~26,580 and 24,790 mortalities annually for males and females, respectively (2). It is reported that ~20% of all patients with CRC are diagnosed with rectal cancer, while ~45% of patients with rectal cancer are further diagnosed with locally advanced rectal cancer (LARC) (3,4).

Local recurrence is common in patients with LARC. Neoadjuvant chemoradiotherapy (nCRT) may prevent local recurrence and improve overall survival, and has therefore emerged as the standard treatment for patients with LARC, followed by surgical resection (5). Overall, following nCRT, ~15–20% of patients with LARC demonstrate a complete response (CR), with complete tumor eradication (6–8) while 20% of patients exhibit a poor response to nCRT (9–11). Patients with a CR to nCRT subsequently undergo a simple watch-and-wait procedure or endoscopic resection (12). Patients with a poor response to nCRT are associated with substantial adverse effects and high medical costs (13). Thus, identifying valid predictive biomarkers is of high importance.

MicroRNAs (miRNAs/miR) have been proposed as potential biomarkers for predicting the response to treatment and prognosis in rectal cancers. The expression of miR-143 and miR-145 is reported to be significantly upregulated in rectal tumor tissues following treatment, and miR-145 expression demonstrated a close association with tumor regression in patients with rectal cancer (14). Additionally, the distinct expression signatures of miR-16, miR-590-5p and miR-153 in pretreated rectal cancer tissue may predict the response to treatment, miR-519c-3p and miR-561 accurately predict a good or poor response to nCRT in patients with rectal cancer (15). Furthermore, miR-21-5p was demonstrated to be overexpressed in patients with CR (16). However, the role of these biomarkers has not been fully clarified.

Della Vittoria Scarpati et al (17), reported a set of 13 miRNAs with specific signatures that were associated with a pathological CR (pCR) to nCRT based on microarray data (accession no. GSE29298). However, the mechanisms of the biomarkers remain unclear. Therefore, the present study downloaded this data from the Gene Expression Omnibus (GEO) database. Based on the miRNA expression profile, differentially expressed miRNAs were identified between pCR and no pCR (incomplete response) groups, and their target genes were predicted, followed by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, miRNA co-regulatory network construction and disease-associated gene analysis. The aim of the current study was to investigate specific miRNA signatures as potential biomarkers for predicting a CR to nCRT in rectal cancer and to determine the potential mechanisms of the miRNAs.

Materials and methods

Data acquisition

The miRNA expression profile of GSE29298 (17) was downloaded from the GEO (http://www.ncbi.nlm.nih.gov/geo/), which was sequenced on an Agilent-021827 Human miRNA Microarray platform (Agilent Technologies, Inc., Santa Clara, CA, USA). The miRNA expression profile was analyzed by microarray using fresh frozen biopsies. A total of 38 patients with LARC (cT3-4/N+) were treated with capecitabine-oxaliplatin and pelvic conformal radiotherapy (45cGy) followed by surgery (after 6–8 weeks). Pathological responses were scored according to the tumor regression grade (TRG), as described by Mandard (9). The number of patients with TRG1, TRG2, TRG3 and TRG4 was 9, 16, 10 and 3, respectively. TRG1 indicates pCR, which represents complete tumor regression, while TRG >1 represents an incomplete response (no pCR) with residual tumor cells. Patients were divided into two groups: pCR group (TRG1, n=9) and no pCR group (TRG >1, n=29).

Data preprocessing

Microarray data were preprocessed using the Bioconductor package (version 3.5) Limma (version 3.28.21) (18), which involved background correction, normalization and expression calculation. The miRNA ID was transformed from the probe ID according to the probe annotation file. The most recent human miRNA annotation file was downloaded from the miRBase Database (19,20) to obtain the miRNA name from the miRNA ID.

Identification of differentially expressed miRNAs

Differentially expressed miRNAs in the pCR group compared with the no pCR group were identified by non-paired Student's t test analysis and linear models using the Bioconductor package Limma. P<0.05 and|log2 fold change (FC)|≥0.58 were set as the cutoff criteria for differentially expressed miRNAs.

Prediction of target genes for differentially expressed miRNAs

Target genes regulated by the differentially expressed miRNAs were predicted by ‘validated target module’ using the miRWalk 2.0 online tool (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2) (21).

KEGG pathway enrichment analysis

KEGG pathways enriched in miRNA target genes were analyzed using the clusterProfiler (version 3.3.1) of R package (version 3.3.2) (22). Pathways involving miRNAs were determined. The P-value was calculated using a hypergeometric distribution statistical approach based on the KEGG.db annotation package (version 3.2.3) (23). The cutoff value for a significant pathway was set at P<0.05.

miRNA co-regulatory network

According to the association between miRNAs and target genes, miRNA-miRNA pairs that regulate the same target genes were screened and an miRNA co-regulatory network was constructed by Cytoscape software (version 3.2.0) (24).

Disease-associated gene analysis

Marker or therapeutic genes associated with colorectal neoplasms were downloaded from the Comparative Toxicogenomics Database (25). Subsequently, the regulatory network between differentially expressed miRNAs and disease-associated genes was constructed by Cytoscape software (version 3.2.0) (24).

Results

Differentially expressed miRNAs

Following preprocessing, 835 miRNAs were identified and subjected to non-paired Student's t test analysis. According to P<0.05 and |log2FC| ≥0.58, a total of 41 miRNAs were differentially expressed in the pCR group compared with the no pCR group, including 36 upregulated and 5 downregulated miRNAs. The heat map demonstrated that the two groups were distinguished by the expression profiles of 41 differentially expressed miRNAs (Fig. 1).

Figure 1.

Heatmap for differentially expressed miRNAs. The pink bar represents the pCR group and the green bar represents the no pCR group. miR/miRNA, microRNA; nCRT, neoadjuvant chemoradiotherapy; pCR, complete response to nCRT; no pCR, incomplete response to nCRT.

Target genes regulated by differentially expressed miRNAs

As indicated in Table I, 3,989 target genes were predicted to be regulated by 41 differentially expressed miRNAs. The results revealed that a single miRNA regulates multiple genes. miR-548c-5p, miR-548d-5p and miR-663a were demonstrated to regulate the greatest number of genes, and were reported to regulate 274, 277 and 127 genes, respectively.

Table I.

Target genes regulated by microRNAs in pathological complete response group compared with the pathological incomplete response group.

A, Upregulated miRNAs
miRNA	Target gene count
hsa-miR-1183	108
hsa-miR-1224-5p	84
hsa-miR-1246	52
hsa-miR-125a-3p	292
hsa-miR-1268a	49
hsa-miR-134-5p	95
hsa-miR-1471	3
hsa-miR-188-5p	136
hsa-miR-483-5p	145
hsa-miR-548c-5p	274
hsa-miR-548d-5p	277
hsa-miR-575	141
hsa-miR-622	118
hsa-miR-630	83
hsa-miR-659-3p	68
hsa-miR-663a	127
hsa-miR-671-5p	164
hsa-miR-765	355
hsa-miR-1182	82
hsa-miR-1207-5p	223
hsa-miR-1225-5p	49
hsa-miR-1226-5p	78
hsa-miR-1275	155
hsa-miR-1299	84
hsa-miR-149-3p	733
hsa-miR-150-3p	44
hsa-miR-1909-5p	50
hsa-miR-1914-3p	117
hsa-miR-202-3p	253
hsa-miR-370-3p	111
hsa-miR-371a-5p	324
hsa-miR-494-3p	154
hsa-miR-557	207
hsa-miR-584-5p	79
hsa-miR-601	18
hsa-miR-623	262
B, Downregulated miRNAs
miRNA	Target gene count
hsa-let-7e-5p	711
hsa-miR-1260a	170
hsa-miR-192-3p	149
hsa-miR-26a-5p	541
hsa-miR-30b-5p	477

miR/miRNA, microRNA.

Significant KEGG pathways of differentially expressed miRNAs

Significant pathways involving the top 11 differentially expressed miRNAs are illustrated in Fig. 2. For miR-548c-5p and miR-548d-5p, the significantly enriched pathways included ‘signaling pathways regulating pluripotency of stem cells’ and ‘ubiquitin-mediated proteolysis’.

Figure 2.

Enriched pathways for the top 11 differentially expressed miRNAs. Different dot color and dot size represent the P-value and gene ratio, respectively. The number below each miRNA name represents the number of target genes involved in pathway analysis. miR/miRNA, microRNA; STAT, signal transducer and activator of transcription; HIF, hypoxia inducible factor; EGFR, epidermal growth factor receptor; HTLV, human T-lymphotropic virus; TGF, transforming growth factor.

miRNAs and their co-regulatory target genes

A co-regulatory network with 41 nodes and 596 edges was constructed, as demonstrated in Fig. 3. The results indicated that the miR-202-3p/let-7e-5p pairing is responsible for the co-regulation 156 genes, while the miR-548c-5p/miR-548d-5p complex was demonstrated to co-regulate 183 genes, including interleukin (IL)-6 signal transducer (IL6ST), cell cycle checkpoint kinase 2 (CHEK2), mutated in colorectal cancers (MCC), marker of proliferation Ki-67 (MKI6), damage-specific DNA binding protein 1 and actinin α4.

Figure 3.

Co-regulatory network among miRNAs. Grey filled ovals represent downregulated miRNA and non-filled ovals represent upregulated miRNA in the pCR group compared with the no pCR group. miR/miRNA, microRNA; nCRT, neoadjuvant chemoradiotherapy; pCR, complete response to nCRT; no pCR, incomplete response to nCRT.

Rectal cancer-associated genes

A total 33 differentially expressed miRNAs associated with the regulation of rectal cancer-associated genes were analyzed. The regulatory network between 33 miRNAs and 80 marker or therapeutic genes associated with colorectal neoplasms was constructed (Fig. 4). Four genes, including cyclin-dependent kinase inhibitor 1B (CDKN1B), adenomatous polyposis coli (APC), MCC and interferon β1 (IFNB1), were demonstrated to be associated with rectal cancer. The following seven miRNAs were identified as regulators of the four rectal cancer-associated genes: miR-548c-5p/miR-548d-5p (MCC), miR-663a (APC), miR-149-3p (CDKN1B), miR-1909-5p (CDKN1B), miR-494-3p (MCC) and miR-26a-5p (IFNB1). Specifically, the miR-548c-5p/miR-548d-5p complex was reported to co-regulate IL6ST, CHEK2, MCC and MKI67, and miR-663a was demonstrated to regulate cadherin 7 (CDH7), calreticulin (CALR), APC and transforming growth factor (TGF) β1 (TGFB1).

Figure 4.

Regulatory network between miRNAs and target genes associated with colorectal neoplasms. Diamonds represent target genes, with pink diamonds representing genes associated with rectal cancer and all other genes being associated with colorectal neoplasms. Grey filled ovals represent downregulated miRNA and non-filled ovals represent upregulated miRNA in the pCR group compared with the no pCR group. miR/miRNA, microRNA; nCRT, neoadjuvant chemoradiotherapy; pCR, complete response to nCRT; no pCR, incomplete response to nCRT; MCC, mutated in colorectal cancers; CDKN1B, cyclin-dependent kinase inhibitor 1B; IFNB1, interferon β1; APC, adenomatous polyposis coli; MKI67, marker of proliferation Ki-67; IL6ST, interleukin-6 signal transducer; CHEK2, cell cycle checkpoint kinase 2; CDH7, cadherin 7; CALR, calreticulin; TGFB1, transforming growth factor β1.

Discussion

It was previously reported that nCRT improves the outcome of patients with rectal cancer; however, not all patients benefit from this treatment (15). Thus, identifying predictive molecular biomarkers is of clinical importance (16). The present study aimed to identify specific miRNAs as predictive biomarkers based on their expression profiles. Furthermore, the potential underlying mechanisms of the biomarkers involved in the response to nCRT were investigated. In the current study, 36 upregulated and 5 downregulated miRNAs were detected between pCR and no pCR (incomplete response) groups.

KEGG pathway analysis demonstrated that the differentially expressed miR-134-5p and miR-630 were significantly associated with the ‘HIF-1 signaling pathway’. Hypoxia-inducible factor (HIF)-1 is a transcription factor that affects gene expression by regulating oxygen delivery and deprivation (26). In addition, the activation of the HIF pathway is reported to contribute to tumor growth. Clinical evidence demonstrated that HIF-1 was overexpressed in patients with a poor prognosis in various cancer types (27). In the present study, miR-134-5p and miR-630 were demonstrated to be significantly upregulated in patients that exhibited a CR to nCRT. The expression of miRNA is negatively associated with the expression of its target genes. Thus, the target genes involved in the HIF-1 signaling pathway were significantly downregulated in CR cases, which may inhibit the HIF-1 signaling pathway and contribute to tumor regression. These results indicated that the bioinformatic findings were promising in the present study. Additionally, three miRNAs among the top 11 miRNAs, including miR-548c-5p, miR-548d-5p and miR-663a, were demonstrated to regulate genes associated with colorectal neoplasms.

It has been previously reported that miR-548c-5p inhibited the proliferation, migration and invasion, and promoted the apoptosis, of liver cancer stem-like cells (28). In addition, miRNA-548d-5p was reported to have a complementary role in supporting oncogenicity in cervical cancer (29). However, few studies have reported the association between these two miRNAs and rectal cancer. Co-regulatory network analysis indicated that miR-548c-5p and miR-548d-5p may function as a complex to co-regulate four genes associated with colorectal neoplasms, including IL6ST, CHEK2, MKI67 and MCC. High IL-6 cytokine levels were previously reported to be associated with an advanced stage of disease and decreased survival of patients with CRC (30). A previous study demonstrated that IL6ST mediated tumor cell proliferation and apoptosis through glycoprotein 130 activation, with subsequent signaling through signal transducer and activator of transcription 3 in CRC (31). CHEK2 is a tumor suppressor that functions in the p53 pathway of the DNA damage response (32). The CHEK2 variants I157T and 1100delC were reported to increase the risk of CRC (33,34). MKI67 encodes a nuclear protein that is strictly associated with cellular proliferation, which means it is regarded as a proliferative marker (35). Antigen Ki-67, encoded by MKI67, is reported to be highly expressed in proliferative cells (36) and a previous study demonstrated that increased Ki-67 expression post-therapy was associated with improved disease-free survival in patients with rectal cancer (37). In addition, another study reported that the extent of the tumor response to nCRT was closely associated with the expression level of Ki-67 (38). Therefore, the expression of Ki-67 may be a potential biomarker for predicting the radiosensitivity of tumors for CRT in rectal cancers. A colorectal mutant cancer protein encoded by MCC, located on chromosome 5q21, was identified as a tumor suppressor gene and is considered to negatively regulate cell cycle progression (39). Allele loss in MCC is common in CRC (40). These findings indicate that miR-548c-5p and miR-548d-5p may be involved in the mechanism of the CR to nCRT by targeting genes associated with cell proliferation, DNA damage and the cell cycle.

Furthermore, KEGG pathway analysis revealed that ‘signaling pathways regulating pluripotency of stem cells’ and ‘ubiquitin-mediated proteolysis’ signaling pathways were significantly enriched for the miR-548c-5p/miR-548d-5p complex. Cancer stem cells have the capacity for tumor initiation by differentiating into heterogeneous lineages of cancer cells (41). The ubiquitin-dependent proteolysis system is involved in the regulation β-catenin turnover (42) and it is established that the accumulation of β-catenin leads to the transcription of pre-oncogenes (43). These findings indicate that the miR-548c-5p/miR-548d-5p complex may contribute to the CR to nCRT through the pluripotency of stem cells and ubiquitin-mediated proteolysis pathways.

It was previously reported that miR-663 exhibited an anti-cancer effect by targeting the TGFB1 transcript through the TGFβ signaling pathway in SW480 human colon cancer cells (44), which is consistent with the results of the present study. In the present study, three other genes associated with colorectal neoplasms were demonstrated to be regulated by miR-663a, including CDH7, CALR and APC. Furthermore, E-cadherin (CDH1), a member of the CDH family, determines cell-cell adhesion and increases the nuclear translocation of β-catenin (44,45). It is well-document that the APC protein, encoded by APC, prevents the accumulation of β-catenin, a key factor in cancer (46). CDH7 is also a member of the CDH family (47). Thus, it can be suggested that miR-663a may function by altering the levels of β-catenin through targeting CDH7 and APC. Calreticulin, encoded by CALR, was reported to be associated with the infiltration of T-cells, and low expression of CALR may represent a novel mechanism underlying immune escape in colon cancer (48). Thus, miR-663a may be closely associated with a CR to nCRT by targeting these four genes.

The aforementioned results support the predictive roles of miR-548c-5p, miR-548d-5p and miR-663a in the response of patients with rectal cancer to nCRT. The present study further analyzed the miRNA microarray data developed from patients treated with capecitabine-oxaliplatin and pelvic conformal radiotherapy (45 cGy) followed by surgery (after 6–8 weeks). Due to the lack of the clinical samples under the same conditions, experimental validations were not available, which is a limitation of the present study. Thus, evaluations in larger, prospective trials are required to validate these biomarkers.

In conclusion, miR-548c-5p, miR-548d-5p and miR-663a are potential biomarkers for predicting a CR to nCRT in rectal cancer. The miR-548c-5p/miR-548d-5p complex may function by targeting IL6ST, CHEK2, MKI67 and MCC. In addition, this complex may function through the pluripotency of stem cells and ubiquitin-mediated proteolysis signaling pathways. miR-663a may affect the tumor response to treatment by regulating CDH7, CALR, APC and TGFB1. Furthermore, targeting these miRNAs with oligonucleotides may represent a potential therapy for improving a poor response to nCRT in patients with rectal cancer.