miRNA signature identification of retinoblastoma and the correlations between differentially expressed miRNAs during retinoblastoma progression.

PURPOSE: Retinoblastoma (RB) is a common pediatric cancer. The study aimed to uncover the mechanisms of RB progression and identify novel therapeutic biomarkers.
METHODS: The miRNA expression profile GSE7072, which includes three RB samples and three healthy retina samples, was used. After data normalization using the preprocessCore package, differentially expressed miRNAs (DE-miRs) were selected by the limma package. The targets of the DE-miRs were predicted based on two databases, followed by construction of the miRNA-target network. Pathway enrichment analysis was conducted for the targets of the DE-miRNAs using DAVID. The CTD database was used to predict RB-related genes, followed by clustering analysis using the pvclust package. The correlation network of DE-miRs was established. MiRNA expression was validated in another data set, GSE41321.
RESULTS: In total, 24 DE-miRs were identified whose targets were correlated with the cell cycle pathway. Among them, hsa-miR-373, hsa-miR-125b, and hsa-miR-181a were highlighted in the miRNA-target regulatory network; 14 DE-miRs, including hsa-miR-373, hsa-miR-125b, hsa-miR-18a, hsa-miR-25, hsa-miR-20a, and hsa-let-7 (a, b, c), were shown to distinguish RB from healthy tissue. In addition, hsa-miR-25, hsa-miR-18a, and hsa-miR-20a shared the common target BCL2L11; hsa-let-7b and hsa-miR-125b targeted the genes CDC25A, CDK6, and LIN28A. Expression of three miRNAs in GSE41321 was consistent with that in GSE7072.
CONCLUSIONS: Several critical miRNAs were identified in RB progression. Hsa-miR-373 might regulate RB invasion and metastasis, hsa-miR-181a might involve in the CDKN1B-mediated cell cycle pathway, and hsa-miR-125b and hsa-let-7b might serve as tumor suppressors by coregulating CDK6, CDC25A, and LIN28A. The miRNAs hsa-miR-25, hsa-miR-18a, and hsa-miR-20a might exert their function by coregulating BCL2L1.

Introduction

Retinoblastoma (RB) is the most common pediatric malignancy. The epidemiology of RB has been investigated in many population-based studies, and the incidence rate of RB is estimated at 40–60 cases per million live births in the world [1-4]. Although the mortality rate is low in patients with RB who receive aggressive multimodal therapy, in developed countries, nearly half of patients with advanced bilateral RB suffer from partial or full sight loss [5]. In developing countries, because the disease is often diagnosed at later stages, the survival rate is lower than in developed countries [6].

To reduce morbidity and preserve the sight of a child during the early stages of RB, numerous studies have explored more effective approaches, such as targeted therapy using bioinformatics methods. For instance, combining epigenetic analysis with gene expression analysis, Zhang et al. identified the important oncogene spleen tyrosine kinase (SYK; Gene ID: 6850; OMIM: 600085), which is elevated in RB and essential for RB tumor cell survival [7]. Another study also discovered 119 candidate genes, such as CDC25C (Gene ID: 995, OMIM: 157680), CDC6 (Gene ID: 990, OMIM: 602627), and TP53 (Gene ID: 7157 OMIM: 191170), for RB diagnosis [8]. MicroRNAs (miRNAs) are small noncoding RNAs that play significant roles in cellular functions and physiology. By regulating the expression of the target genes, miRNAs are confirmed to be involved in the development of various cancers, and thus have been suggested as tumor biomarkers [9,10]. Several miRNAs such as miR-30, miR-let-7e, miR-21, and miR-320 are dysregulated in RB samples and have been supposed to be diagnostic biomarkers for detecting RB [11,12]. Downregulated miR-204 is another indicator in RB prediction [13]. Martin et al., using a TaqMan Low Density Array, discovered a total of 41 differentially expressed miRNAs (DE-miRs) between 12 RB samples and three healthy retina samples in humans, including 13 previously identified miRNAs (miR-17, miR-20b, miR-22, miR-25, miR-34a, miR-34c-5p, miR-106a, miR-106b, miR-93, miR-129, miR-193–5p, miR-342–5p, miR-370) in human and mouse RB and many novel miRNAs, such as miR-18a, miR-138, miR-155, miR-382, and miR-504 [14]. Additionally, the miR-17–92 cluster has been demonstrated as an RB-collaborating gene that promotes RB development [15]. More recently, another 18 miRNAs have been newly implicated in RB and have great potential to serve as signatures in the detection of this disease [16]. However, the target genes of these miRNAs are rarely reported. Notably, using paired mRNA and miRNA expression profiles, Huang et al. identified several targets of miRNAs in RB samples and further verified CDC25A (Gene ID: 993 OMIM: 116947) and BCL7A (Gene ID: 605 OMIM: 601406) are the target genes of let-7b [17]. However, the researchers emphasized the roles of miRNA let-7b and did not mention other potential miRNAs or the correlations between them. In addition, the detailed regulation mechanisms of miRNAs to RB remain obscure.

Therefore, we reanalyzed the miRNA expression profile GSE7072 [17] to obtain more relevant miRNAs using differential analysis. The targets of these miRNAs were also predicted using two experimental validated databases (miRecords and MirWalk). Relationships between these miRNAs were further explored to comprehensively uncover the underlying mechanisms of RB progression. We aimed to find novel miRNA biomarkers for the prevention and prognosis of RB development.

Methods

A flowchart of the analyses in the study is shown in Figure 1.

Figure 1

Flowchart of the analyses.

Microarray data

The miRNA expression profile data with the accession number GSE7072 [17], which is available in the public Gene Expression Omnibus (GEO) database, was employed in the present study. The data set comprised the total RNA information of a cohort of 160 human miRNAs from three RB samples and three replicates of a healthy retina, based on the platform of the GPL4879Human miRNA 2k custom array (Agilent Technologies, Palo Alto, CA). The annotation files on the platform were downloaded.

Data preprocessing and identification of DE-miRs

Based on the annotation information, the probe levels were converted into miRNA expression values. The probe that did not correspond to a specific miRNA was removed, and when more than one probe corresponded to a single miRNA, the average value at the probe level was calculated as the final expression value of this miRNA. Then the data were subjected to normalization using the median method in the preprocessCore package [18]. Afterwards, the DE-miRs between the RB and healthy retina samples were selected using the limma (Linear Models for Microarray Analysis) package of R [19]. The cut-off values for significant DE-miRs were p<0.05 and |log2 (fold change)| >0.58.

Construction of integrated miRNA–target network

Considering that a miRNA works through the regulation of the target in a spectrum of biologic processes, we further explored the potential target genes of these identified DE-miRs, by integrating the information in two experimentally validated databases, the miRecords [20] and MirWalk [21], in which miRNA–target interactions were experimentally validated. Only the miRNA–target interaction that existed in at least one of the two databases was screened out to establish the integrated miRNA–target network.

Pathway enrichment analysis of the predicted target genes

To further identify the altered biologic pathways of the target genes, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted with the Database for Annotation, Visualization and Integration Discovery (DAVID) software [22], based on the hypergeometric distribution. p<0.05 was set as the threshold for the selection of significant pathways.

Identification of transcription factors, tumor-associated genes, and tumor suppressor genes

The target genes were mapped into the TRASFAC database [23] to further screen the potential transcription factors (TFs). Meanwhile, by comparing these targets with the information from the tumor suppressor gene (TSG) [24] and tumor-associated gene (TAG) databases, all known oncogenes and TSGs among the target genes were extracted.

Prediction of RB-related genes among the target genes

The target genes were mapped into the publicly available Comparative Toxicogenomics Database (CTD) database, which provides curated chemical- and gene-disease interactions from published articles [25], to search the RB-related genes among these target genes. Then combined with the corresponding miRNAs, clustering analysis was performed between the RB and healthy retina samples with the pvclust package of R [26].

Correlation network construction of the identified DE-miRs

Among the DE-miRs, two miRNAs that share a common target gene were filtered out to build the correlation network.

Validation of miRNA expression with another data set

Another RB-related miRNA expression profile, GSE41321 [27], which also contained three RB samples and three healthy samples, was downloaded in the GEO database and used to check whether the DE-miRs identified in this profile were consistent with those in the GSE7072 data set. Samples in the GSE41321 profile were collected from pooled serum from children with advanced RB and healthy age-matched children. Likewise, after expression value conversion from probe level to gene level, normalization was performed for data using the median method in the preprocessCore package [18]. Then the miRNA expression value was attained. The CONOR package [28] was used to perform cross-platform normalization for the GSE7072 and GSE41321 data sets. DE-miRs between the RB and control samples were identified in the two data sets, respectively, using the limma package with the same threshold. Thereafter, we compared the miRNA expression in the two data sets, especially the 24 DE-miRs identified in GSE7072.

Results

The DE-miRs between the RB and healthy retina samples

As a result, a set of 24 DE-miRs were selected, including nine upregulated (e.g., hsa-miR-25, hsa-miR-373, and hsa-miR-20a) and 15 downregulated miRNAs (e.g., let-7b, let-7a, let-7c, hsa-miR-125b, and hsa-miR-181a) between the RB and healthy retina samples. The clustering analysis of the heat map is presented in Figure 2, revealing the expression pattern of the DE-miRs.

Figure 2

Heat map of the clustering analysis of differentially expressed miRNAs between two samples. The x-axis represents the samples (healthy 1–3 denote healthy retina samples, and Rb 1–3 denote retinoblastoma samples), and the y-axis represents miRNAs.

The integrated miRNA–target regulatory network

Based on the criterion, the integrated miRNA–target regulatory network was built, containing five upregulated miRNAs, including hsa-miR-373 (the most prominent one with 53 target genes), hsa-miR-20a (with nine target genes), hsa-miR-18a, hsa-miR-25, and hsa-miR-175p, and 12 downregulated miRNAs, such as hsa-miR-125b (the most remarkable one with 54 target genes), hsa-miR-7 (a, b, c), and hsa-miR-145 (Figure 3). Notably, hsa-miR-20a, hsa-miR-18a, and hsa-miR-25 shared the common target gene BCL2L11; hsa-miR-125b, hsa-miR-7a, and hsa-miR-7b cotargeted the gene LIN28A.

Figure 3

Integrated miRNA–target regulatory network. The blue circle represents genes, and the red rectangle represents upregulated miRNAs, while the green rectangle represents downregulated miRNAs. Arrows denote direct interactions between a miRNA and its target gene.

Dysregulated pathways of the target genes

After the target genes were mapped into the KEGG databases, the altered pathways were identified. As shown in Table 1, the target genes of the upregulated miRNAs were significantly enriched in numerous cancer-related pathways such as pathways in cancer (hsa05200), bladder cancer (hsa05219), non-small cell lung cancer (hsa05223), and pancreatic cancer (hsa05212). The over-represented pathways for the targets of the downregulated miRNAs were also significantly correlated with various cancers, as well as the ErbB signaling pathway (hsa04012) and the cell cycle pathway (hsa04110).

Table 1

Significantly enriched pathways among the target genes of the differentially expressed miRNAs.

Target genes	Pathway Term	Count	P value
Target genes of upregulated miRNAs	hsa05200:Pathways in cancer	6	3.62E-04
	hsa05219:Bladder cancer	3	3.49 E-03
	hsa05223:Non-small cell lung cancer	3	5.73 E-03
	hsa05212:Pancreatic cancer	3	1.00 E-02
	hsa05220:Chronic myeloid leukemia	3	1.08 E-02
Target genes of downregulated miRNAs	hsa05219:Bladder cancer	10	2.93E-10
	hsa05200:Pathways in cancer	19	2.71E-09
	hsa05220:Chronic myeloid leukemia	11	3.75E-09
	hsa05214:Glioma	10	1.31E-08
	hsa05215:Prostate cancer	11	2.06E-08
	hsa05218:Melanoma	10	3.87E-08
	hsa05223:Non-small cell lung cancer	9	6.76E-08
	hsa04012:ErbB signaling pathway	9	2.91E-06
	hsa04110:Cell cycle	10	5.14E-06
	hsa05212:Pancreatic cancer	8	8.85E-06

Count, the gene numbers enriched in a specific pathway.

Identified TFs, TAGs, and TSGs among target genes

As expected, the potential TFs, oncogenes, and TSGs were selected by comparing the targets with the information in relevant databases. Detailed gene information is presented in Table 2, in which target genes such as CBFB (Gene ID: 865; OMIM: 121360), CEBPG (Gene ID: 1054; OMIM: 138972), and HMGA2 (Gene ID: 8091; OMIM: 600698) were considered TFs; CCNA2 (Gene ID: 890; OMIM: 123835), CCND1 (Gene ID: 595; OMIM: 168461), and ERBB2 (Gene ID: 2064; OMIM: 164870) were oncogenes; and CDKN1B (Gene ID: 1027; OMIM: 600778), CDKN2A (Gene ID: 1029; OMIM: 600160), and E2F1 (Gene ID: 1869; OMIM: 189971) were TSGs.

Table 2

Predicted TFs, oncogenes and TSGs among the target genes of the differentially expressed miRNAs.

TF	Oncogene	TSG
CBFB,CEBPG,HMGA2,HOXA11,ID1,ID2,ID3,MEF2D,MITF,PLAGL1,PRDM1,RARG,RUNX1,SOX9,TFAP4,ZHX1	CCNA2,CCND1,ERBB2,ERBB3,ESR1,FGFR3,HMGA2,HRAS,KRAS,MYC,NRAS	CDKN1B,CDKN2A,E2F1,FOXO1,ITGB3,LATS2,MAP2K4,NF2,PLAGL1,PPP1CA,PRDM1,RECK,SMO,TP53,TP53INP1,TUSC2,UHRF2

TF: transcription factor; TSG: tumor suppressor gene.

RB-related miRNAs with their targets

A total of 183 target genes of the DE-miRs were identified. By mapping them into the CTD database, 87 were considered to correlate with RB. Additionally, these genes were targeted by 14 DE-miRs, including hsa-miR-373, hsa-miR-125b, hsa-miR-20a, hsa-miR-145, hsa-let-7 (a, b, c), hsa-miR-25, hsa-miR-18a, hsa-miR-182, hsa-miR-99a, hsa-miR-183, hsa-miR-451, and hsa-miR-181a, which completely distinguished the RB samples from the healthy samples (Table 3 and Figure 4). Therefore, these miRNAs might be used as the biomarkers of RB.

Table 3

Retinoblastoma-related miRNAs and their target genes.

miRNAs	Target genes
hsa-miR-373	GPSM2, KIF23, GLTP, LATS2, TFAP4, GBP3, PRC1, LMNB1, RNF149, INSIG2, TEX30, CYB5R4, RPIA, PBK, STK4, KDM1B, CENPF, TBC1D2, ARHGEF3, CD83
hsa-miR-125b	ID2, CASP6, ADAMTS1, NTRK3, ID3, JARID2, CBX7, UGT2B17, ELAVL1, BAK1, EIF4EBP1, TENM2, E2F3, MAN1A1, KRT19, TP53INP1, TP53, CDC25A, CDKN2A, CDK6, SMO, BMPR1B, CYP1A1, CASP7, BBC3, CBLN2, UBE2I, IL1RN, B3GALT4, CDC25C, ERBB3, ERBB2, PIGR
hsa-miR-20a	BCL2L11, CCND1, BMPR2, RUNX1, E2F1, MAP3K12
hsa-let-7b	CDC25A, CDK6, CCND1, HMGA2, KRT19
hsa-miR-145	FSCN1, SOX9, IGF1R, RTKN, CCNA2, MYC, PPP3CA
hsa-let-7a	NF2, SMOX, KRAS, ITGB3, CASP3, HMGA2, HRAS, NRAS
hsa-miR-182	MITF, IGF1R, RARG
hsa-miR-99a	MTOR, IGF1R, FGFR3
hsa-miR-25	MAP2K4, BCL2L11
hsa-miR-18a	CTGF, BCL2L11, ESR1
hsa-miR-181a	CDKN1B, ESR1
hsa-miR-183	FOXO1
hsa-let-7c	MYC, HMGA2
hsa-miR-451	ABCB1

Figure 4

Clustering analysis of the 14 differentially expressed miRNAs in two types of samples. Healthy 1–3 denote the healthy retina samples, and Rb 1–3 denote the retinoblastoma samples. The values on the edges of the clusters are p values (%). The red values are approximately unbiased (au) p values, and the green values are bootstrap probability (bp) values. Clusters with au larger than 95% are highlighted by rectangles, which are strongly supported by data.

The correlation network of the DE-miRs

The correlation network showing the common target genes of the DE-miRs is shown in Figure 5. As presented in this network, the predominant correlated miRNA interactions were hsa-let-7b and hsa-miR-125b (the common target genes were CDC25A, LIN28A, and CDK6); and hsa-miR-18a, hsa-miR-20a, and hsa-miR-25 (the common target gene was BCL2L11).

Figure 5

Correlation network of the differentially expressed miRNAs. The blue circle represents genes, and the red rectangle represents upregulated miRNAs, while the green rectangle represents downregulated miRNAs. The arrow denotes the direct interaction between miRNA and its target gene.

Validation of miRNA expression

As shown in Table 4, in comparison with the GSE7072 profile, ten miRNAs, including six identified DE-miRs (out of the 24 DE-miRs), were also differentially expressed in the RB samples in the GSE41321 profile. Notably, the expression of three DE-miRs, including upregulated hsa-miR-18a and downregulated hsa-let-7b and hsa-let-7c, was in accordance with our results.

Table 4

miRNA (10 DE-miRs) expression comparison between expression profiles of GSE41321 and GSE7072.

miRNAs	GSE41321		GSE7072
	Log FC	P value	Log FC	P value
hsa-miR-30d	0.699283069	8.34E-06	−1.4135091	1.28E-02
hsa-miR-18a	0.820197258	3.05E-05	1.486038867	3.38E-02
hsa-miR-181a	0.872595328	1.01E-06	−1.9342523	5.73E-04
hsa-miR-183	0.561178808	1.26E-04	−1.2667168	1.94E-02
hsa-miR-182	0.95936798	5.54E-07	−1.24761365	2.40E-02
hsa-miR-125b	0.457144179	4.40E-04	−2.520854685	3.56 E-04
hsa-miR-191	−0.483798552	1.04E-03	−1.80388745	7.45E-03
hsa-miR-373	−0.109119276	0.20	1.196955933	9.69E-03
hsa-let-7c	−0.838050781	6.21E-05	−2.213052267	3.67E-03
hsa-let-7b	−0.704607286	4.49E-05	−3.044672233	9.35E-05

DE-miRs: differentially expressed miRNAs (between retinoblastoma and control samples); FC: fold change; miRNAs with positive value of Log FC represent upregulated miRNAs in retinoblastoma samples; while with negative value indicate downregulated miRNAs in retinoblastoma samples.

Discussion

RB is a primary pediatric cancer of the retina. In this study, we performed a series of analyses using powerful bioinformatics methods and finally identified a total of 24 DE-miRs between the RB and healthy retina samples. Interestingly, nine DE-miRs, including let-7c, hsa-miR-182, hsa-miR-99b, hsa-miR-125b, hsa-miR-191, hsa-miR-181a, hsa-miR-17–5p, hsa-miR-373, and let-7b, were consistent with work performed by Huang et al. [17], who identified 25 DE-miRs in RB samples. The potential reasons for the other different miRNAs might be due to different selection criteria. As Huang et al. aimed to discover the interacting miRNA–target relationships using paired miRNA and mRNA expression profiles, they might have applied a different criterion to choose the interacting miRNA target. In contrast, this study used their miRNA expression profile and predicted the corresponding targets based on the miRecords and MirWalk databases. However, the nine common DE-miRs suggest that our findings are also reliable.

Interestingly, the expression of three DE-miRs was consistent with that in another profile, GSE41321, including hsa-miR-18a, hsa-let-7b, and hsa-let-7c. Although the GSE41321 samples are collected from serum, Beta et al. confirmed that a set of 33 deregulated miRNAs, such as the upregulated hsa-miR-18a and the downregulated hsa-let-7c, are also found in RB tumor samples in comparison with the tumor RB miRNA profiles [27]. This finding indicates several miRNAs in serum might have the same expression patterns as in tumors and could be used as potential biomarkers for the prognosis and prediction of RB development. Enrichment analysis indicated that the target genes of the present identified DE-miRs were mainly involved in cancer-related pathways and the cell cycle pathway. Among these DE-miRs, hsa-miR-373, hsa-miR-125b, and hsa-miR-181a were highlighted with multiple target genes in the integrated miRNA–target regulatory network; meanwhile, 14 DE-miRs, such as hsa-miR-373, hsa-miR-125b, hsa-miR-20a, hsa-miR-25, hsa-miR-18a, and hsa-let-7 (a, b, c) were shown to distinguish RB from healthy samples. Notably, hsa-miR-25, hsa-miR-18a, and hsa-miR-20a shared the common target gene BCL2L11; hsa-let-7b and hsa-miR-125b targeted three genes, CDC25A, CDK6, and LIN28A.

MiR-373 is implicated in numerous cancer events. Previously, miR-373–3 was considered the signature of human embryonic stem cells [29]. In prostate cancer, in vitro experiments showed that miR-373 stimulates cell migration and invasion [30]. Additionally, the role of miR-373 as an activator in tumor invasion and metastasis has been verified in vivo and in vitro by inhibiting the expression of CD44 [31]. Furthermore, miR-373 is proproliferative in human epithelial ovarian cancer cells [32]. Several studies corroborated that miR-373 is highly expressed in RB tumors, in comparison with healthy retinas [12,33], which is consistent with our result. Notably, our enrichment analysis showed that the targets of hsa-miR-373 were tightly correlated with a cohort of cancer-related pathways, providing potent evidence that hsa-miR-373 might function in RB tumors as in other cancers. Given that RB has a tendency toward local invasion and metastasis [34], it might be speculated that miR-373 might also play remarkable roles in the regulation of RB invasion and metastasis.

Defects in the gene RB1, which serves as a negative regulator of cell cycle and a tumor suppressor, are known to be the main cause for the genesis of RB in children. A cluster of miRNAs acts as regulators in cell proliferation, apoptosis, and senescence via interfering with the p53 pathway or RB1/E2F function [35,36]. MiR-125b has been confirmed to function as a tumor suppressor by inducing cellular senescence in cutaneous malignant melanoma [37]. Notably, miR-125b is closely related to the p53 pathway in colorectal cancer [38]. These collectively provide a clue that miR-125b might also be involved in the regulation of proliferation during RB development via the p53-mediated pathway. In human U251 glioma stem cells, the overexpression of miR-125b is highly correlated with the decreased expression of CDK6 and CDC25A, two cell cycle–mediated genes [39]. Moreover, CDK6 and CDC25A are mediated by miR-125b in other types of cancer, such as breast cancer [40] and bladder cancer [41]. The oncogene LIN28B has been demonstrated to be the downstream of miR-125b, and by binding to the 3′ untranslated region (UTR) region of LIN28B, miR-125b successfully inhibited the cell growth in hepatocellular carcinoma [42]. The regulatory relationships between hsa-let-7b and the target CDC25A have been well established in RB cells [17]. Additionally, let-7 is suggested to indirectly regulate the cell cycle by suppressing IMP1, an essential gene for CDC25A and CDK6 expression [43]. Furthermore, let-7 could also downregulate the mRNA of LIN28B during neural stem cell commitment [44]. Interestingly, our results showed that CDK6, CDC25A, and LIN28A, the paralog of LIN28B, were all targeted by miR-125b and hsa-let-7b, indicating these two miRNAs might play crucial roles in the development of RB, via the coregulation of the three cell cycle–related target genes.

Hsa-miR-181b and hsa-miR-181a are downregulated in human gliomas and have been implicated as tumor suppressors by inhibiting glioma proliferation [45]. In RB cells, elevated miR-181b has been verified to play a positive role in RB cell proliferation under hypoxic conditions [46]. Considering that RB is a cancer of the nervous system and miR-181b is dysregulated in the central nervous system disease schizophrenia, miR-181b has been suggested to be involved in the genesis of RB [46]. In addition, miR-181a, along with miR-183 and miR-124, presents the most frequent alterations among the miRNAs differentially expressed in RB [47]. The CDKN1B gene, an important regulator of the cell cycle, has been demonstrated as the target of miR-221 and miR-222 in the development of thyroid cancer [48]. Our results indicated that hsa-miR-181a was a crucial miRNA of RB that targeted multiple genes, including CDKN1B, which was predicted to be a TSG enriching in the cell cycle pathway. Therefore, it might be inferred that hsa-miR-181a might act as an important regulator in RB progression via the CDKN1B-mediated cell cycle pathway.

The BCL2L11 encoded protein belongs to the BCL-2 protein family, which act as important apoptotic activators in extensive cellular activities [49]. Transforming growth factor-β (TGF-β) is a cytokine that has a significant role in the apoptosis of gastrointestinal cells. miR-25 overexpression decreases TGFβ-induced apoptosis by interfering with the synthesis of BCL2L11 in gastric cancer [50]. Additionally, miR-20a has also been reported to target BCL2L1, while miR-18a could target Smad2 and Smad4 in the TGF-β signaling pathway, in which the activation of TGF-β is partly regulated by BCL2L11 [51]. In the present study, hsa-miR-25, hsa-miR-18a, and hsa-miR-20a simultaneously targeted BCL2L1 in the correlation network of the DE-miRs, suggesting that the three miRNAs might have vital roles in the progression of RB by coregulating BCL2L1. However, the potential interactions between these miRNAs must be further elucidated.

The present study has several limitations. First, the sample size was small because only three RB samples and three healthy samples were contained in the data sets. Moreover, although we used another data set to validate the miRNA expression, the expression of only three miRNAs was confirmed. Experimental validation of miRNA expression and miRNA–target regulation was lacking, which will be considered in follow-up studies.

In conclusion, a set of critical miRNA signatures was identified in RB progression. Among them, hsa-miR-373 might play significant roles in the regulation of RB invasion and metastasis, hsa-miR-125b and hsa-let-7b might exert their roles as tumor suppressors via the coregulation of CDK6, CDC25A, and LIN28A, which all mediated the cell cycle pathway, and hsa-miR-181a might involve in the CDKN1B-regulated cell cycle pathway. Additionally, hsa-miR-25, hsa-miR-18a, and hsa-miR-20a might affect the progression of RB by the coregulation of BCL2L1. However, additional experimental studies are needed to confirm these results.