Ubiquitination-Related miRNA–mRNA Interaction Is a Potential Mechanism in the Progression of Retinoblastoma

Abstract

Purpose: Retinoblastoma (RB) is the most common primary malignant intraocular cancer. The etiology of RB is complex, and the mechanisms driving its progression remain unclear. Here, we used a series of bioinformatics approaches and experimental methods to investigate the potential regulatory mechanism involved in RB progression.

Methods: The common differentially expressed genes were obtained from the public dataset GSE97508. Protein–protein interaction (PPI) network, correlation, and functional enrichment analyses were carried out. The candidate genes were verified in different RB cell lines, and ARPE19 cells served as control. miRNA–mRNA interaction analysis was performed and confirmed by real-time PCR. The CCK-8 assay was conducted to detect cell viability, and the transwell assay was utilized for evaluating the abilities of cell migration and invasion.

Results: Overall, a total of 258 common differentially expressed genes associated with RB progression were screened out. The PPI network analysis further identified eight downregulated genes mainly enriched in the protein ubiquitination pathway. Moreover, we confirmed UBE2E1, SKP1, FBXO9, FBXO15, and RNF14 from among eight genes through experimental validation in vitro. Furthermore, miRNA–mRNA interaction and real-time PCR analysis of five hub genes revealed that ubiquitination-related miR-548k was involved in RB progression. Loss- and gain-of-function experiments demonstrated that miR-548k and its targets were essential for cell viability, migration, and invasion in the RB cells.

Conclusions: Our data indicate that the dysregulation of protein ubiquitination may play an important role in RB progression, and ubiquitination-related miR-548k may be a promising therapeutic target for RB.

Retinoblastoma (RB) is the most common primary intraocular malignancy of childhood, with an incidence rate of 1 in 15,000 to 20,000 live births, corresponding to about 9000 new cases each year.¹ Although the survival rate has been enhanced through early diagnosis and treatment over the past 20 years, RB still accounts for 3% of all pediatric malignancies.² The survival rate and retention of vision depend on tumor aggressiveness. The progression of RB has not been well defined so far; therefore, determining the evolution of RB would contribute to impeding tumor progression and increasing therapy efficiency. 

Recent advances in bioinformatics approaches have been important in analyzing, comparing, and interpreting data both on a global scale and at a macroscopic level.^3,4 Finding target genes or pathways using bioinformatics may save scientists literally years of work at the lab. Many studies have shown that using comprehensive bioinformatics analyses is a proper strategy to uncover the potential genes and underlying mechanisms in various types of tumors, such as breast, lung, and liver cancers.^5–7 Thus, we can use bioinformatics approaches to explore the gene expression characteristics and potential regulative mechanisms to gain a better understanding of the molecular basis of RB progression and identify novel genomic targets for therapeutic intervention. 

Considering the complexity and heterogeneity of RB, we adopted microarray technology and bioinformatics methods to systematically explore patterns of ubiquitination-related genes and potential signaling pathways involved in RB development. This study analyzed gene expression data for RB from the Gene Expression Omnibus (GEO) database, and the common differentially expressed genes (DEGs) were acquired from the intersection of two sub-datasets (noninvasive tumor vs. normal retina, invasive tumor vs. noninvasive tumor). Protein–protein interaction (PPI) and functional enrichment analyses were performed, and the results reflect the important link between ubiquitination and RB. The prediction of potential miRNAs and loss- or gain-of-function experiments further confirmed the molecules related to ubiquitination. This study could provide new insights useful for future investigations of the malignant transformation of RB. 

Materials and Methods

Data Acquisition

The microarray expression profiles were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).⁸ The GSE97508 dataset, already deposited in GPL15207 (Affymetrix Human Gene Expression Array), consisted of normal retina, noninvasive RB, and invasive RB from humans. 

Data Processing and DEG Screening

The raw data of GSE97508 were preprocessed using the R 4.0.0 affy package (version 1.66.0; R Foundation for Statistical Computing, Vienna, Austria) and the robust multichip averaging algorithm.^9,10 Gene probe identification, based on the biomaRt package (version 2.44.0) and tidyverse package (version 1.3.0), was matched to the corresponding gene symbol after background correction and data normalization. DEGs among various groups in GSE97508 were identified using the empirical Bayes approach in linear models form the limma package (version 3.43.11),¹¹ and a heatmap of overlapping DEGs was drawn using the pheatmap package (version 1.0.12; https://github.com/raivokolde/pheatmap) to visualize the differential expression of DEGs with the progression of RB. DEGs with adjusted P < 0.05 and log2FC > 1 were considered statistically significant. 

PPI Network Construction and Gene Cluster Identification

The online Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) is a primary source that can be used to describe and display the functional organization of the proteome.¹² We uploaded overlapping DEGs onto the STRING website and identified significant interactions with interaction scores ≥ 0.9 (high confidence). The STRING analysis results were imported into Cytoscape 3.8.0, and cluster analysis of differential genes was conducted via the Molecular Complex Detection (MCODE) plug-in.¹³ The nodes with degree > 1 were reserved in the PPI network. The top eight genes sorted by degree were considered hub genes. The MCODE parameter criteria were executed by default. 

Correlation Analysis

The correlation between each pair of hub genes in different groups was analyzed by Pearson's correlation and visualized using the corrplot package (version 0.84; https://github.com/taiyun/corrplot) in R. The criterion for the correlation was set as P < 0.05. 

Gene Ontology and Pathway Enrichment Analysis

Gene Ontology (GO) term enrichment is a technique used to evaluate the characteristics of sets of genes, and GO annotation includes biological process (BP), cellular component (CC), and molecular function (MF).¹⁴ The clusterProfiler package (version 3.16.0) was used to perform this analysis.¹⁵ Reactome pathway analysis and visualization were conducted with the ReactomePA package (version 3.7.1).¹⁶ The criterion for significant enrichment was set as adjusted P < 0.05. 

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) is a computational method that determines whether the preset gene set has statistically significant and concordant differences between two biological states.¹⁷ GSEA was performed on the gene expression matrix using the clusterProfiler package, and “c5.bp.v7.1.entrez.gmt” was selected as the reference gene set. In addition, we employed the enrichplot package (version 1.7.4; https://github.com/YuLab-SMU/enrichplot) to visualize the GSEA enrichment data. The criteria of P < 0.05, false discovery rate (FDR) < 0.25, and normalized enrichment score (NES) > 1 were set. 

Construction of the mRNA–miRNA Interaction Network

Hub genes were selected for assessment with miRWalk 2.0 (http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/), which incorporates 12 algorithms for target prediction.¹⁸ The target miRNAs identified by at least seven algorithms were selected for further analysis. The results obtained were further visualized with Cytoscape. miRNAs that targeted more than three genes were selected. The criteria were set as follows: P < 0.05, minimum seed sequence length of 7-mer, and 3′UTR target gene binding region. 

Cell Culture

A human normal retinal pigmented epithelium cell line, APRE-19, the low-invasive RB cell line WERI-RB1, and the high-invasive RB cell line Y79 were all purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). ARPE-19 cells were cultivated in Dulbecco's Modified Eagle's Medium (DMEM; Invitrogen, Carlsbad, CA, USA) containing 10% heat-inactivated fetal bovine serum (FBS). Y79 and WERI-RB1 cells were nurtured in modified Roswell Park Memorial Institute (RPMI)-1640 medium (Gibco, Carlsbad, CA, USA) including 10% FBS. All cells were cultured in a humidified atmosphere with a mixture of 1% O₂, 5% CO₂, and 94% N₂ at 37°C. 

Real-Time Quantitative Reverse Transcription PCR

Isolation of total RNAs from cells was executed with TRIzol Reagent (Invitrogen) and Dr. GenTLE Precipitation Carrier (Takara Biomedical Technology, Beijing, China). Real-time quantitative reverse transcription (qRT-PCR) assays were performed according to the manufacturer's protocol for the SYBR PrimeScript TM RT-PCR Kit or Mir-X miRNA First-Strand Synthesis Kit (Takara Biomedical Technology). Relative target gene expression was quantitated according to the comparative ΔCT method; that is, it was normalized to an endogenous control gene, β-actin or U6, and relative to a calibrator after calculating the efficiency coefficient: relative expression = 2–Δ∆CT, where Δ∆CT = CT (target gene) – CT (β-actin/U6). The fold differences of gene expression among different groups were calculated using the 2–ΔΔCt method. The primers of mRNA and miRNA are shown in Supplementary Table S1. 

Cell Transfection and Lentivirus Infection

The miR-548k mimics and miR-548k inhibitor were purchased from RiboBio Company (Guangzhou, China). The above factors were transfected into WERI-RB1 and Y79 cells using Lipofectamine 3000 (Invitrogen) for 48 hours. Small interfering RNA (siRNA) against SKP1 (si-SKP1), RNF14 (si-RNF14), and non-specific control siRNA (si-NC) were purchased from Ribo (Guangzhou, China). SKP1 and RNF14 were cloned into pLJM1 lentiviral vectors to infect Y79 cells, and empty vectors served as control. 

Cell Viability Assay (CCK-8)

Forty-eight hours after transfection, the WERI-RB1 or Y79 cells were inoculated into a 96-well plate with every well containing 5 × 10³ cells. Cells were kept in the incubator with 5% CO₂ overnight at 37°C. Before the CCK-8 was added, the medium was changed. For every well, 100 µL of medium and 10 µL of CCK-8 were added. The dish was kept in the incubator with 5% CO₂ at 37°C for 4 hours. The optical density at 450 nm was measured via an enzyme immunoassay instrument (BioTek Instruments, Winooski, VT, USA). The optical density values were detected after culturing for 0, 24, 48, and 72 hours. 

Cell Migration and Invasion Assays

The cell migration assay was performed using a 24-well transwell chamber (Corning, Corning, NY, USA). Briefly, 24 hours after transfection, WERI-RB1 or Y79 cells (1 × 10⁵) were resuspended in 200 µL serum-free medium and seeded to the upper chamber with an 8-µm pore size insert precoated with Matrigel (BD Biosciences, Franklin Lakes, NJ, USA). The plate wells were filled with 600 µL RPMI-1640 containing 10% FBS. After incubation for 48 hours at 37°C, cells on the upper side of the membrane were removed with clean swabs, and cells adhering to the lower chambers were fixed with 4% paraformaldehyde and stained with 4′,6-diamidino-2-phenylindole (DAPI). Under an inverted fluorescence microscope, the image was magnified 200 times to count the number of invaded cells in each of the five randomly selected fields. The experiments were performed in triplicate. 

For the cell invasion assay, the operation was similar to the cell invasion assay, except that the 24-well transwell chambers were precoated with Matrigel. 

Luciferase Reporter Assay

Luciferase assays were measured at 48 hours following transfection using a Dual Luciferase Reporter Assay System (Promega, Madison, WI, USA). The Firefly luciferase activity was normalized to Renilla luciferase activity. 

Statistical Analysis

Data were expressed as mean ± standard deviation (SD). Statistical significance was estimated using R 4.0.0. Data were analyzed using Tukey's honestly significant difference (HSD) test for one-way analysis of variance. The Benjamini–Hochberg correction method is utilized for multiple hypothesis testing in gene set enrichment analysis with a significance level of 0.05. All of the tests were two sided, and P < 0.05 was considered to be statistically significant. 

Results

Screening and Identification of DEGs

First, we obtained 1164 upregulated and 1109 downregulated genes between noninvasive RB samples and normal subjects from the GSE97508 dataset. We identified 6531 DEGs, including 4146 upregulated and 2385 downregulated DEGs, when comparing the invasive RB group with the noninvasive RB group (Fig. 1A). The screening criteria for DEGs were defined as adjusted P < 0.05 and log2FC > 1. Moreover, we screened the 26 common upregulated and 232 common downregulated DEGs (Fig. 1B), and the results are shown in Supplementary Table S2. Based on the analysis, a heatmap plot was constructed to visualize the upregulated and downregulated genes related to the progression of RB (Figs. 1C, 1D). 

Figure 1.

 

Screening and identification of DEGs. (A) Bar charts show the number of both upregulated and downregulated DEGs from two sub-datasets. (B) Venn diagrams show the number of upregulated or downregulated DEGs; the overlapping part represents the common upregulated or downregulated DEGs in the two datasets. (C) Heatmap of common upregulated DEGs. (D) Heatmap of common downregulated DEGs. The color in the panel of each heatmap represents the relative expression level (log2-transformed).

Module Analysis and Gene Correlation Analysis Based on PPI Network

A PPI network based on DEGs (67 nodes, 90 edges) was obtained to explore the potential interaction among RB-related molecules; the network had an interaction score ≥ 0.9 according to the STRING online database. Based on the criterion of filtering of node degree > 8, the top eight hub genes were SKP1, WWP1, UBE2G2, RNF144B, UBE2E1, RNF14, FBXO15, and FBXO9, all of which were downregulated (Fig. 2A). Then we used the MCODE algorithm to analyze a subset of the network. When setting the node density cutoff to 1, node score cutoff to 0.2, k-core to 2, and maximum depth to 100, we found that genes in the first gene cluster with the highest score were the exact eight hub genes mentioned above. Notably, we found that these genes were all involved in ubiquitination by literature mining^19–23; two encoded the E2 ubiquitin-conjugating enzyme and the other six were related to E3 ubiquitin ligase (Fig. 2B). We also performed Pearson’s correlation test to investigate the relationship among the hub genes in different groups. As shown in Figure 2C, the relevance between each pair of correlations differed and there were even differences among the various stages of RB. Interestingly, FBXO15 was inversely related to SKP1 in different groups (red squares). The FBXO15 and UBE2G2 pair and the FBXO9 and UBE2E1 pair each had a strong positive relationship in both invasive and noninvasive tumor samples but inversely negative relevance in the control (green squares). Based on these findings, we conjectured that downregulation of eight ubiquitination-related genes might be involved in the progression of RB, and further experiments to confirm the interactions among these genes are necessary. 

Figure 2.

 

The PPI networks and correlation analysis of the hub genes. (A) Construction of the PPI network based on the DEGs. The network nodes represent proteins, and the edges represent the protein–protein associations. The red ellipses represent upregulated genes, the blue ellipses represent downregulated genes, and the yellow ellipses represent hub genes. (B) The module with the highest scores in the PPI network; these genes were classified based on text mining in PubMed. (C) Correlation heatmap of eight hub genes. The size of the colored squares represents the strength of the correlation; red represents a positive correlation, and blue represents a negative correlation. The darker the color is, the stronger the correlation. P < 0.05 was considered statistically significant. The red frame indicates stable alterations of correlation results; the green frame indicates significant alterations of correlation results that were of interest.

Ubiquitin-Mediated Regulation May Be Involved in RB Progression

To reveal the potential function of hub genes, we performed functional and pathway enrichment analyses of the eight genes. As shown in Figure 3A, GO analysis results indicate that hub genes were mainly related to ubiquitin-associated pathways and regulation, such as protein polyubiquitination (BP, adjusted P = 3.00E-13), ubiquitin–protein transferase activity (MF, adjusted P = 2.43E-11), and ubiquitin ligase complex (CC, adjusted P = 8.88E-10). The top five of three sub-ontologies (in the CC category, there were only four items) were screened out according to the adjusted P < 0.05, as shown in Figure 3A and Supplementary Table S3. Moreover, the pathway analysis demonstrated that a total of 118 pathways were associated with at least one protein from hub genes. Forty-seven pathways were all significantly enriched (P < 0.05). Fourteen pathways remained statistically significant after Benjamini–Hochberg correction, which was used to control the FDR and decrease the number of false positives (q < 0.01) (Supplementary Table S4). Four of the top six enriched pathways were relevant to ubiquitination, including antigen processing: ubiquitination and proteasome degradation (adjusted P = 3.42E-11), synthesis of active ubiquitin and roles of E1 and E2 enzymes (adjusted P = 0.0053), neddylation (adjusted P = 0.0101), and protein ubiquitination (adjusted P = 0.0143) (Fig. 3B). It is known that DEG analysis may miss important effects on pathways; therefore, we performed the GSEA to further confirm the results above. The GSEA evaluates the significant pathways with genome-wide expression profiles, which could provide more reproducible, interpretable, and reliable results. As shown in Figure 3C, two gene sets (noninvasive RB vs. normal retina, invasive RB vs. noninvasive RB) were both significantly negatively enriched in the regulation of ubiquitin-dependent protein catabolic processes (Fig. 3C). These data further indicated that ubiquitin-mediated regulation is gradually suppressed with the progression of RB. 

Figure 3.

 

Functional annotation and enrichment analysis of the hub genes. (A) The general GO annotations for cellular component, molecular function, and biological processes; the horizontal axis represents the number of DEGs under the GO term. (B) The cnetplot of the hub gene pathway analysis. (C) GSEA was used to analyze the whole gene expression value of two sub-datasets. GSEA first filtered the gene set according to the number of genes contained in the gene set, with a minimum number of 15 genes and a maximum number of 500 genes set by default. Significant gene sets were cut off by P < 0.05, FDR < 0.25, and NES > 1.

Experimental Validation of Hub Genes

To verify whether these hub genes were credible, we assayed the expression pattern of these genes by qRT-PCR in the low-invasive RB cell line WERI-RB1 and high-invasive RB cell line Y79, and ARPE19 cells served as control. As shown in Figure 4, FBX015, SKP1, RNF14, UBE2E1, and FBXO9 levels were significantly reduced as the invasiveness of cells increased: FBXO15 (F_2,6 = 153.9, P < 0.001); SKP1 (F_2,6 = 19.16, P = 0.003); RNF14 (F_2,6 = 20.15, P = 0.002); RNF144B (F_2,6 = 6.591, P = 0.031); UBE2G2 (F_2,6 = 10.09, P = 0.012); UBE2E1 (F_2,6 = 216.1, P < 0.001); FBXO9 (F_2,6 = 46.29, P < 0.001); and WWP1 (F_2,6 = 5.509, P = 0.044). Thus, these results made our determination of five key genes more reliable. 

Figure 4.

 

Validation of the eight selected genes using real-time qRT-PCR, which showed the relative mRNA levels of FBOX15, SKP1, WWP1, RNF14, UBE2G2, UBE2E1, RNF144B, and FBXO9 in different cell lines (n ≥ 3). ^*P < 0.05, ^**P < 0.01, ^***P < 0.001. n.s., no significance.

miRNA–mRNA Interaction Network Analysis

Based on the previous studies,^24–27 we speculated that miRNAs might play important roles in the inhibition of ubiquitination; therefore, we performed gene and miRNA analysis among these five key genes with miRWalk 2.0 software. The target miRNAs were selected and identified by at least seven of 12 algorithms. Cytoscape was used to create the interaction network, as shown in Figure 5. In addition, three candidate miRNAs (miR-548k, miR-7-1-3p, and miR-7-2-3p) that had a high number of cross-links (≥3) were selected for further experimental validation. 

Figure 5.

 

Interaction network of genes involved in ubiquitination and their targeting miRNAs. Genes are shown in blue, miRNAs in red, and miRNAs targeting more than three genes simultaneously in yellow.

miR-548k Regulated Cell Viability, Migration, and Invasion in RB Cells

To determine these candidate miRNAs, we performed qRT-PCR in the different RB cell lines. As shown in Figure 6A, only the level of miR-548k was significantly elevated as the invasiveness of cells increased: miR-548k (F_2,6 = 18.51, P < 0.001); miR-7-1-3p (F_2,6 = 31.53, P = 0.026); and miR-548k (F_2,6 = 35.53, P = 0.032). Afterward, further analyses and experiments were performed to explore the role of miR-548k, and 1098 targets were selected in the miR-548k–mRNA network by at least seven of 12 algorithms (Supplementary Table S5), and pathway enrichment analysis was performed. Seven pathways were significantly enriched, and one of the top three enriched pathways was relevant to ubiquitination (Supplementary Figure S1). As shown in Figure 6B, miR-548k was strikingly boosted by transfection of the miR-548k mimic in WERI-RB1 cells and significantly reduced by transfection of the miR-548k inhibitor in Y79 cells (WERI-RB1 miR-548k mimic vs. NC, P = 0.013); Y79 miR-548k inhibitor vs. NC, P = 0.004). We subsequently clarified the regulatory interrelation between miR-548k and its potential targets and confirmed that the mRNA expression of UBE2E1, SKP1, and RNF14 were reduced by upregulation of miR-548k in WERI-RB1 cells (UBE2E1, P < 0.001; SKP1, P < 0.001; RNF14, P = 0.039), whereas transfection of the miR-548k inhibitor significantly promoted the expression of these genes (UBE2E1, P = 0.039; SKP1, P = 0.009; RNF14, P = 0.026) (Figs. 6C, 6D). 

Figure 6.

 

miR-548k regulated cell viability, migration, and invasion in RB cells. (A) qRT-PCR showed the relative miRNA level of miR-548k, miR-7-1-3p, and miR-7-2-3p in different cell lines (n ≥ 3). ^*P < 0.05, ^***P < 0.001. (B) Relative expression of miR-548k in WERI-RB1 and Y79 cells after transfection with miR-548k mimic or inhibitor was detected by qRT-PCR (n ≥ 3). ^*P < 0.05, ^**P < 0.01. (C) Relative expression of UBE2E1, SKP1, RNF14 in WERI-RB1, and Y79 cells after transfection with miR-548k mimic or inhibitor was detected by qRT-PCR (n ≥ 3). ^*P < 0.05, ^**P < 0.01. (D) Cell viability in WERI-RB1 and Y79 cells after transfection was detected by CCK-8 assay after culturing for 0, 24, 48, and 72 hours (n ≥ 3). ^*P < 0.05, ^*P < 0.01, ^##P < 0.01. (E) Relative migration and invasion of WERI-RB1 and Y79 cells after transfection were detected by transwell assay (n ≥ 3). ^*P < 0.05, ^**P < 0.01.

The CCK-8 assay showed that cell viability was increased by the addition of miR-548k mimic at 3 days after treatment, and the addition of the miR-548k inhibitor caused a distinct decrease in cell viability in Y79 cells (WERI-RB1 miR-548k mimic vs. NC, 24-hour P = 0.027 and 72-hour P = 0.008; Y79 miR-548k inhibitor vs. NC, 72-hour P = 0.005) (Fig. 6D). Transwell assays determined that transfection of the miR-548k mimic promoted migration and invasion of WERI-RB1 cells, and inhibition of miR-548k significantly weakened migration and the invasion of Y79 cells: migration (WERI-RB1 miR-548k mimic vs. NC, P = 0.041; Y79 miR-548k inhibitor vs. NC, P = 0.009); invasion (WERI-RB1 miR-548k mimic vs. NC, P = 0.026; Y79 miR-548k inhibitor vs. NC, P = 0.017) (Fig. 6E). Taken together, ubiquitination-related miR-548k was involved in cell viability, migration, and invasion in RB cells. 

Both SKP1 and RNF14 Were Involved in RB Progression In Vitro

Considering that the E2 conjugation enzyme UBE2E1 has been frequently reported in a variety of human cancers^28,29 and E3 ligases that bind directly to substrate proteins commonly serve as important targets for anti-tumor therapy,³⁰ we further focused on the role of E3 ligases SKP1 and RNF14 in RB progression. Luciferase assays were performed to confirm the predicted interaction between miR-548k and the 3′-UTR of targets. The predicted interactions between miR-548k and the target site in the SKP1 3′-UTR or RNF14 3′-UTR are shown in Figure 7A. The results demonstrate that miR-548k significantly inhibited the luciferase activity of both SKP1 and RNF14 reporter (SKP1, F_2,6 = 83.25, P < 0.001; RNF14, F_2,6 = 10.94, P = 0.007) (Fig. 7B), indicating a direct interaction between miR-548k and its targets. 

Figure 7.

 

Both SKP1 and RNF14 were involved in RB progression in vitro. (A) Predicted binding sites of miR-548k in the 3′-UTR of SKP1 and RNF14. (B) Luciferase assay of Y79 cells co-transfected with miR-548k mimics and a luciferase reporter containing wild-type or mutant 3′-UTR of the target genes. (C) Relative expression of SKP1 and RNF14 in WERI-RB1 and Y79 cells after transfection with siRNA or the pLJM1 vector was detected by qRT-PCR (n ≥ 3). ^*P < 0.05, ^**P < 0.01. (D) Cell viability in WERI-RB1 and Y79 cells after transfection was detected by CCK-8 assay after culturing for 0, 24, 48, and 72 hours (n ≥ 3). ^*P < 0.05, ^**P < 0.01, ^#P < 0.05, ^##P < 0.01. (E, F) Relative migration and invasion of WERI-RB1 and Y79 cells after transfection were detected by transwell assay. Representative images and quantitative results were shown (n ≥ 3). ^*P < 0.05, ^**P < 0.01.

Gene disruption/overexpression experiments and subsequent functional experiments were performed to determine the effect of these two genes on RB malignancy. siRNAs and pLJM1 overexpression vectors were used for silencing and overexpressing the target genes SKP1 and RNF14, respectively. As presented in Figure 7C, the mRNA expression of target genes was significantly decreased in WERI-RB1 cells by siRNA transfection (SKP1, P = 0.027; RNF14, P = 0.015). In contrast, transfection of the pLJM1-SKP1 or pLJM1-RNF14 lentiviral vectors into Y79 cells increased respective mRNA levels as compared with empty vector–transfected cells (SKP1, P = 0.005; RNF14, P = 0.021). Transient siRNA silencing of both SKP1 and RNF14 promoted the proliferation, migration, and invasion of the low-invasive retinoblastoma cell line WERI-RB1: CCK-8 (si-SKP1 vs. NC, 48-hour P = 0.032 and 72-hour P = 0.026; si-RNF14 vs. NC, 72-hour P = 0.006); migration (si-SKP1 vs. NC, P = 0.021; si-RNF14 vs. NC, P = 0.008); invasion (si-SKP1 vs. NC, P = 0.011; si-RNF14 vs. NC, P = 0.019) (Figs. 7D–7F), whereas overexpression of these two genes weakened the deterioration of the high-invasive retinoblastoma cell line Y79: CCK-8 (SKP1 vs. vector, 48-hour P = 0.023 and 72-hour P = 0.028; RNF14 vs. vector, 48-hour P = 0.033 and 72-hour P = 0.016); migration (SKP1 vs. vector, P = 0.002; RNF14 vs. vector, P = 0.018); invasion (SKP1 vs. vector, P = 0.021; RNF14 vs. vector, P = 0.003) (Figs. 7D–7F). 

These data indicate that the downregulation of two ubiquitination-related proteins (SKP1 and RNF14) improves RB progression in vitro. 

Discussion

Retinoblastoma, which originates from the progenitors of retinal sensory cells, is the most common primary malignant intraocular cancer in children.³¹ The development of RB is a complex process involving the interaction of multiple genes and signaling pathways in tumor cells and stromal cells. In this study, we sought to identify key genes in RB tumor progression and further explore the role of signaling pathways involved in RB. 

After analyzing two selected sub-datasets (noninvasive RB vs. normal retina, invasive RB vs. noninvasive), the overlapping DEGs were identified whose expression was gradually upregulated or downregulated with the progression of RB. Eight genes of interest (UBE2E1, UBE2G2, SKP1, FBXO9, FBXO15, RNF144B, RNF14, and WWP1) were identified by the PPI network and key module analyses (Fig. 2A). Further organizing the results of the literature mining revealed that the eight genes were all associated with ubiquitination, and these genes can be divided into two categories: E2 conjugation enzymes (UBE2E1 and UBE2G2) and E3 ligases (SKP1, FBXO9, FBXO15, RNF144B, RNF14, and WWP1) (Fig. 2B). GO analysis also demonstrated that these eight genes are involved in synthesis and degradation of ubiquitination proteasome (Fig. 3). RT-PCR assay further confirmed that five genes (FBXO15, SKP1, RNF14, UBE2E1, and FBXO9) were significantly decreased in RB tumors compared with normal samples (Fig. 4). Thus, the above illustration and results suggest that there are some inherent links between ubiquitination and RB progression. 

Our results are supported by previous studies. For example, a growing number of studies have demonstrated that ubiquitination mediates cellular homeostasis and contributes to various processes in pathophysiological states and diseases, ranging from cancer to infection and hereditary disorders.^32–34 The ubiquitin-proteasome system plays an important role in the occurrence and development of cancer.³⁵ UBE2E1 can modulate PTEN ubiquitination and transport and serves as a tumor suppressor.^28,29 Thus, it is reasonable to suggest that downregulation of UBE2E1 promotes RB progression. Moreover, the F-box proteins bound to SKP1 function as substrate recognition subunits.²⁰ Thompson et al.^36,37 reported that diminished SKP1 expression induced irregular increases in cyclin E1 levels, which led to DNA double-strand breaks and chromothripsis events and indirectly promoted tumor metastasis. Notably, we also found that downregulation of SKP1 improved RB progression in vitro (Figs. 7D–7F). It is worth mentioning that we established a xenotransplantation model that has been described in detail in a previous article,³⁸ and we found that the mRNA levels of SKP1 were significantly decreased when the deterioration of RB increased in vivo (Supplementary Fig. S2). The underlying mechanisms of SKP1-mediated tumor suppression require further investigation. Furthermore, SKP1 is a critical component of the Skp1–Cullin1–F-box (SCF) complex, which shares a common catalytic core consisting of SKP1, the scaffold protein Cullin1, and the RING finger protein Rbx1. Many studies have reported the role of the SCF complex in the development of breast, colon, prostate, lung, and gastric cancers, among others.³⁹ In the current study, we identified two F-box proteins, FBXO9 and FBXO15, that were involved in the progression of RB. A previous study reported that FBXO9 promoted survival in multiple myeloma through ubiquitination.⁴⁰ Loss of FBXO9 accelerated the progression of acute myeloid leukemia, serving as a tumor suppressor in acute myeloid leukemia.⁴¹ FBXO15 has been reported to be a critical gene in maintaining normal development and physiology in mice.¹⁹ The knockdown of FBXO15 increased cancer-associated drug resistance by P-glycoprotein accumulation.⁴² Above all, the relationship between the SCF complex and RB progression is concerning. 

Moreover, downregulation of genes is generally related to miRNA. Real-time PCR further verified that miR-548k was significantly upregulated as the invasiveness of cells increased. Importantly, the two genes described above (SKP1 and RNF14) interacted with and were regulated by miR-548k, as demonstrated by luciferase reporter assays and a series of gene disruption–overexpression experiments (Figs. 6C, 7A, 7B). miR-548k was involved in cell viability, migration, and invasion in RB cells (Figs. 6D, 6E). Bioinformatics analysis (Supplementary Fig. S1) and experimental results (Fig. 6) further confirmed that miR-548k impaired the ubiquitin-proteasome system and indirectly promoted the progression of retinoblastoma. These results were partially supported by previous studies. miR-548k has been reported to be an oncogene in esophageal squamous cell carcinoma (ESCC).^43–45 Elevated expression of miR-548k was detected in ESCC cell lines compared with normal esophageal cells, which was consistent with our results. Also, the upregulated miR-548k promoted advanced invasion depths and lymph node metastases in ESCC through lncRNA-LET.⁴⁴ Zhang et al.⁴³ demonstrated that miR-548k promoted ESCC metastasis by regulating the KLF10/EGFR axis. Moreover, dysregulated miR-548k may lead to immune deficiency by targeting CXCL13 and was found to be related to head and neck cancers.^46,47 

Based on our data and the previous studies mentioned above, we may illustrate the underlying mechanism of RB progression as shown in Figure 8. In RB cells, miR-548k inhibits the expression of UBE2E1, RNF14, and SKP1 and indirectly suppresses the ubiquitination of some molecules that are associated with the viability, migration, and invasion of RB cells. Thus, miR-548k exerts a positive effect on the invasiveness of RB cells. In summary, miR-548k may impair ubiquitination by targeting E2 conjugating enzymes UBE2E1 and E3 ligases (RNF14 and SKP1) and promote RB progression. Because we confirmed that both SKP1 and RNF14 were involved in RB progression (Figs. 7D–7F), it would be valuable to further explore which substrate proteins are activated due to dysregulation of the ubiquitin–proteasome system. 

Figure 8.

 

miR-548k impairs the ubiquitin–proteasome system and indirectly promotes RB progress by targeting the E2 ubiquitin-conjugating enzyme UBE2E1 and E3 ligases RNF14 and SKP1. (1) miR-548k downregulates the expression of the E2 ubiquitin-conjugating enzyme UBE2E1 in RB cells. (2) miR-548k inhibits the expression of the E3 ligases RNF14 and SKP1 in RB cells. (3) Downregulation of UBE2E1, SKP1, and RNF14 impedes the ubiquitination and degradation of some potential substrate proteins. (4) Potential substrate proteins accumulate in RB cells. (5) Potential substrate proteins result in tumor progression.

The expression level of three genes (UBE2G2, RNF144B, and WWP1) was not consistent in the bioinformatics analysis (Fig. 4). However, the ubiquitin-conjugating enzyme UBE2G2, known as the specific E2 of tumor autocrine motility factor receptor gp78, is involved in endoplasmic reticulum–associated degradation.⁴⁸ Zhao et al.⁴⁹ determined that UBE2G2 is positively associated with the overall survival of non–small cell lung cancer patients, and its expression is significantly decreased in tumor samples compared with normal lung tissues. RNF144B has been shown to play an oncogenic role in regulating chordoma,⁵⁰ and WWP1, a HECT-type E3 ubiquitin ligase, has been demonstrated to be overexpressed in various cancers.^23,51 Perhaps different mechanisms exist during the development of RB. In addition, the FBXO15–UBE2G2 and FBXO9–UBE2E1 gene pairs had strong positive relations in both noninvasive and invasive tumor samples, which were not related to the invasive ability of RB. Details regarding the interactions among genes remain unclear. 

Functional annotation and pathway analysis provided further evidence of the relationship between the eight hub genes and ubiquitination (Fig. 3). Notably, a fraction of the genes were enriched in pathways but it seemed to be unrelated to ubiquitination, including class I major histocompatibility complex (MHC)-mediated antigen processing and presentation and RUNX2 expression and activity regulation. Further correlations among these ubiquitination-related genes in RB remain to be determined. No previous study has reported the relationship between RNF14 and RB, but we found that RNF14 and SKP1 had similar roles in RB development (Figs. 7D–7F, Supplementary Fig. S2). However, we did not explore deeply the specific mechanisms behind the observed results in which RNF14 was involved. To address all of these uncertainties, further study is required in the future. 

In conclusion, by applying a series of bioinformatics approaches and experimental methods to gene expression profiles, we successfully identified five key genes (UBE2E1, SKP1, FBXO9, FBXO15, and RNF14) and ubiquitination-related miR-548k. Moreover, we propose that the dysregulation of protein ubiquitination may play an important role in the pathogenesis of RB, thus providing novel insights into the progression of RB. 

Acknowledgments

The authors thank Guangzhou Huaxiang Medical Biotechnology Co., Ltd., for bioinformatics support. 

Supported by the National Natural Science Foundation (81470626, 81670848, and 81900850) and Guangzhou Science Technology and Innovation Commission (201804010076). 

Disclosure: X. Chen, None; S. Chen, None; Z. Jiang, None; Q. Gong, None; D. Tang, None; Q. Luo, None; X. Liu, None; S. He, None; A. He, None; Y. Wu, None; J. Qiu, None; Y. Li, None; X. Wang, None; K. Yu, None; J. Zhuang, None