Bioinformatics Analysis and Experimental Identification of Immune-Related Genes and Immune Cells in the Progression of Retinoblastoma

Abstract

Purpose: Retinoblastoma (RB) is the most common type of aggressive intraocular malignancy in children. The alteration of immunity during RB progression and invasion has not yet been well defined. This study investigated significantly altered immune-associated genes and cells related to RB invasion.

Methods: The differentially expressed immune-related genes (IRGs) in noninvasive RB and invasive RB were identified by analysis of two microarray datasets (GSE97508 and GSE110811). Hub IRGs were further identified by real time PCR. The single-sample gene set enrichment analysis algorithm and Pearson correlation analysis were used to define immune cell infiltration and the relationships between hub IRGs and immune cells. Cell viability and migration were evaluated by CCK-8 and Transwell assays. A xenograft mouse model was used to verify the relationship between Src homology 3 (SH3) domain GRB2‐like 2 (SH3GL2) expression and myeloid-derived suppressor cells (MDSCs).

Results: Eight upregulated genes and six downregulated IRGs were identified in invasive RB. Seven IRGs were confirmed by real-time PCR. Moreover, the proportions of MDSCs were higher in invasive RB tissues than in noninvasive RB tissues. Furthermore, correlation analysis of altered immune genes and cells suggested that SH3GL2, Langerhans cell protein 1 (LCP1) and transmembrane immune signaling adaptor TYROBP have strong connections with MDSCs. Specifically, decreased SH3GL2 expression promoted the migration of RB cells in vitro, increased the tumor size and weight, and increased the numbers of MDSCs in the tumor and spleen in vivo.

Conclusions: This study indicated that SH3GL2 and MDSCs play a critical role in RB progression and invasion and provide candidate targets for the treatment of RB.

Retinoblastoma (RB) is the most common pediatric ocular malignancy with a worldwide prevalence of approximately 1/16,000 in newborns.¹ The prognosis and treatment of RB have dramatically improved over the past ten years; the treatments include enucleation, cryotherapy, thermotherapy, brachytherapy, intra-arterial chemotherapy, intravitreal chemotherapy, and systemic chemotherapy. Although the growth of primary tumors has been successfully controlled, current therapies have limited effects on invasive RB.² There is still limited knowledge about the differences between invasive and noninvasive RB. Noninvasive RB means that tumor was localized to the eye, and the invasive RB refers to tumor with invasion to the choroid, optic nerve, and orbit/metastasis.³ Early diagnosis and treatment of invasive RB is critical for the visual quality and survival of patients. Thus distinguishing the characteristics of invasive RB from those of noninvasive RB is important for its clinical treatment. 

A growing amount of literature has indicated that immune genes and cells can infiltrate tumors and play critical roles in tumor progression.^4–6 For example, overexpression of PD-L1 in cancer cells and PD-L1+ immune cell infiltration in Epstein–Barr virus-associated gastric cancer showed a significant correlation with tumor invasion.⁷ Moreover, immune scoring based on the different types, densities and locations of immune cell infiltration in tumors can be applied for the prediction of tumor invasion and prognosis.⁸ Furthermore, immunotherapies for the prevention of tumor progression have evolved rapidly over the past decades, and some studies have indicated the diagnostic and prognostic value of dendritic cells, PD-1 and PD-L1 in RB.^2,9 Our previous study also reported that liver metastasis of RB could be accelerated by exosome-mediated immune cell inhibition.¹⁰ However, studies focusing on immune genes and cells in RB invasion are still quite limited. Considering the important roles of immune genes and cells in tumor invasion, a thorough and systematic investigation of the immune genes expressed and immune cells present in RB tumors with different degrees of invasion is required and may reveal potential underlying mechanisms and clinical implications for RB. 

Of course, the task of elucidating the effects of immune-related genes and cells in invasive progression of RB is daunting with regular experimental methods. In the past decade, bioinformatics analysis has been well utilized to highlight the potential mechanisms of and novel therapeutic targets in tumor progression and invasion.^11–13 Cao et al.¹⁴ identified the relationship between the subtype and progression of RB and demonstrated the role of cell cycle regulation in RB invasion by bioinformatics analysis of four public datasets. However, the immune genes and cells involved in RB invasion have not been investigated through bioinformatics analysis to date. Thus, in the present study, using microarray data analysis and the single-sample gene set enrichment analysis (ssGSEA) algorithm, we sought to determine the changes in immune genes and cells during RB invasion. Subsequently, using interference with gene expression, we investigated the correlations between immune genes and cells during RB progression both in vitro and in vivo. This study is expected to provide new insight into the progression of RB tumors. 

Materials and Methods

Data Acquisition and Annotation

The gene expression profiles of RB were obtained from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) database. GSE97508, including three noninvasive and three invasive RB samples, was collected from the Gene Expression Omnibus. Another dataset was GSE110811, which consisted of 18 noninvasive and 10 invasive RB samples from humans. 

Differential Expression Analysis and Acquisition of Immune-Related Genes

R software (version 4.1.2) was used to perform data preprocessing. The affy (version 1.72.0) package in R language was used to conduct normalization and background correction.¹⁵ For multiple probes matching the same gene, the average expression value was chosen as the targeted gene expression value. The limma (version 3.50.0) package in R was used to obtain differentially expressed genes (DEGs) between the noninvasive and invasive RB samples in the expression data.¹⁶ Differences were considered statistically significant at P < 0.05, and a fold change (FC) > 2 or <−2. Immune-related genes (IRGs) were accessed from the ImmPort (https://immport.niaid.nih.gov) and AmiGO2 databases (http://amigo.geneontology.org/amigo). Differentially expressed IRGs were selected from the IRGs and DEGs. 

Immune Cell Infiltration

SsGSEA was used to compute the relative abundance of 28 immune cell types in the tumor microenvironment based on traditional microarray expression data.¹⁷ Feature gene panels for specific immune cell types were obtained.¹⁸ SsGSEA was applied to calculate a sample ssGSEA score with the GSVA package (version 1.42.0, https://github.com/rcastelo/GSVA) in R. 

Correlation Analysis

The correlation between immune genes and cells was conducted by Pearson analysis and visualized by the corrplot package (version 0.92, https://github.com/taiyun/corrplot) in R. A significance criterion alpha was set as 0.05. 

Cell Culture

Two kinds of RB cell lines (low-invasive WERI-RB1 and the high-invasive Y79) were both purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). These cell lines were both grown in modified Roswell Park Memorial Institute (RPMI)-1640 medium (Gibco; Thermo Fisher Scientific, Waltham, MA, USA) added with 10% fetal bovine serum and incubated at 37°C in an atmosphere of 5% carbon dioxide. The medium was changed every two days. 

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

TRIzol (Invitrogen, Carlsbad, CA, USA) was applied to extract total RNA from cells and tissues. QRT-PCR assays were performed according to the manufacturer's protocol. The 2-ΔΔCt method was used for the relative quantification of target gene expression levels. The primer sequences are shown in Supplementary Table S1. 

Migration Assay

Two kinds of cell lines WERI-RB1 or Y79 cells (1 × 10⁵) were transfected and then suspended in 200 µL RPMI-1640 medium without serum. Seeding them separately into the 8-µm pore size upper chamber (BD Biosciences, Franklin Lakes, NJ, USA) in the 24-well transwell plates (Corning, Corning, NY, USA). The chamber below contained 600 µL RPMI-1640 with 10% FBS. After 48 hours incubation, cells adhering to lower side of the membrane were washed with PBS and counted. Cells that migrated into the lower chamber were captured by imaging microscopy and counted. The experiments were performed in triplicate. 

Cell Viability Assay (CCK-8)

Forty-eight hours after transfection, different cells (1 × 10⁵) were inoculated into a 48-well plate with 5% CO₂ overnight at 37°C. For each well, 200 µL medium and 20 µL CCK-8 reagent were added. Two hours after incubation, the absorbance at 450 nm was detected using an enzyme immunoanalyzer (BioTek Instruments, Winooski, VT, USA). 

Western Blotting

Cell, tumor, and spleen tissues were lysed in RIPA buffer with PMSF. Equal amounts of protein were separated by the sodium dodecyl sulfate/polyacrylamide electrophoresis gel and then the protein was transferred to the nitrocellulose polyvinylidene fluoride membrane. Different primary antibodies were used: rabbit anti-SH3GL2 (1:1000, cat. no12435-1-AP; Proteintech, Chicago, IL, USA) and rabbit anti-TUBULIN (1:1,000; cat. no. sc-5274; Santa Cruz Biotechnology, Dallas, TX, USA). The proteins were then incubated with horseradish peroxidase–conjugated anti-rabbit IgG (1:10000, cat. no. 7074s; Cell Signaling Technology, Inc., Danvers, MA, USA) for one hour at 37°C, and the membranes were then visualized by an enhanced chemiluminescence detection system. 

Immunofluorescence Assay

Tumor and spleen tissue sections were fixed with ice-cold 4% paraformaldehyde for 15 minutes and treated with 0.2% Triton X-100 for 10 minutes at 37°C. Tissues were blocked with 5% bovine serum albumin for 30 minutes. Then the samples were incubated with different primary antibodies, including anti- matrix metalloproteinase 9 (MMP9; 1:200, cat. no., 10375-2-AP; Proteintech), anti-C-X-C chemokine receptor type 4 (CXCR4; 1:300, cat. no. ab181020; Abcam, Cambridge, MA, USA), anti-Gr-1 (1:100, cat. no. 31469; Cell Signaling Technology, Inc.) and anti-CD11b (1:100, cat. no. 17800, Cell Signaling Technology, Inc.). Then the samples were stained with secondary antibodies for one hour at 37°C (Alexa Fluor 555 anti-rabbit IgG, 1:500, cat. no. 4413, Cell Signaling Technology, Inc.; Alexa Fluor 488 goat anti-Rat IgG, 1:500, cat. no. A-11006; Invitrogen). MMP9 is a marker of tumor invasiveness, and CXCR4 (a chemokine receptor marker for migration) labels the tumor progression. MDSCs are labeled by Gr-1 and CD11b. 

H&E Staining

The Liver paraffin sections were put into xylene I for 20 minutes, xylene II for 20 minutes, anhydrous ethanol I for 10 minutes, anhydrous ethanol II for 10 minutes, 95% alcohol for five minutes, 90% alcohol for five minutes, 80% alcohol for five minutes, 70% alcohol for five minutes, and distilled water for five minutes. Then the sections were transferred into Harris hematoxylin for three to eight minutes, tap water, 1% hydrochloric acid alcohol for several seconds, tap water, 0.6% ammonia water, and tap water. Then the sections were stained with eosin staining solution for one to three minutes, transferred into 95% alcohol I for five minutes, 95% alcohol II for five minutes, absolute ethanol I for five minutes, absolute ethanol II for five minutes, xylene I for five minutes, xylene II for five minutes. Neutral gum was used to seal the liver sections. 

Animals

Ten female athymic nude mice four to six weeks of age were purchased from the Ophthalmic Animal Laboratory, Zhongshan Ophthalmic Center, Sun Yat-sen University (Guangzhou, China). All animal experiments adhered to the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and were approved and monitored by the Institutional Animal Care and Use Committee of Zhongshan Ophthalmic Center (approval no. SYXK [YUE] 2019-009). 

Mice were subcutaneously injected in the left subaxillary region with 300 µL of 5 × 10⁶ WERI-Rb1 cells mixed with Matrigel (vol/vol, 1:1). After one week, mice were successfully transplanted tumor and then injected with an equal volume of small interfering RNA (siRNA; in vivo–grade cholesterol- and methylated- conjugated siRNA, 2 nmol; RiboBio, Guangzhou, China) into the primary tumor site every three days for 12 days. Body weight and tumor size were observed and recorded. Tumor volumes were calculated every three days by the following formula: length × width²/2. When the tumor reached about 1000 mm³, animals were killed. The tumors were collected 12 days after injection and stored at −80°C. 

Si-RNA and Plasmid Transfection

In vitro and in vivo grade siRNA against SH3GL2 (si-SH3GL2) and nonspecific control siRNA (si-NC) were obtained from RiboBio (Guangzhou, China) and performed according to the manufacturer's protocol. The Y79 cells were transfected with plasmid mRuby2-N1 (empty vector, OE-NC, cat. no. 54614; Addgene, Watertown, MA, USA) or plasmid EndophilinA1-mRuby2 (overexpression SH3GL2, OE-SH3GL2, cat. no. 171945; Addgene) by Lipofectamine 3000 transfection reagent (Invitrogen) according to the manufacturer's protocol. 

Statistical Analysis

All analyses were conducted using R software (version 4.0.3). An unpaired Student t test was used to compare differences between two groups of normally distributed data. All statistical tests were two sided, and P < 0.05 was considered statistically significant. 

Results

Identification of Immune-Associated Genes Involved in RB Invasion

A flowchart of the analysis procedure for exploring the relationship between tumor invasion and immune regulation is shown in Figure 1 and described in detail in the Methods section. First, we identified 88 upregulated and 89 downregulated genes between 18 noninvasive and 10 invasive RB samples in the GSE110811 dataset. (Fig. 2A). Second, after analyzing the GSE97508 microarray, which included three noninvasive and three invasive RB samples, a total of 9011 DEGs, namely, 4944 upregulated and 4067 downregulated DEGs, were identified. The expression fold changes of these DEGs were displayed in a volcano plot (Fig. 2B). Subsequently, the expression levels of the DEGs described above were visualized in a heatmap (GSE110811: Fig. 2C, GSE97508: Fig. 2D). Moreover, 3255 genes involved in immune-related phenomena were obtained from the KEGG and AmiGO2 databases. After taking the intersection of the DEGs from the two datasets, 14 differentially expressed IRGs that may potentially regulate the aggressiveness of RB (Fig. 2E), namely, eight upregulated DEGs and six downregulated DEGs were identified. 

Figure 1.

 

Schematic representation of the study protocol. The workflow consisted of exploring altered immune-associated genes and immune infiltrating cells and analyzing their relationship during RB development.

Figure 2.

 

Screening and identification of differentially expressed IRGs. (A) Volcano plot of DEGs among invasive RB samples and noninvasive controls from the GSE110811 dataset; (B) Volcano plot of DEGs among invasive RB samples and noninvasive controls from the GSE97508 dataset; (C) Heatmap of DEGs from the GSE110811 dataset; (D) Heatmap of DEGs from the GSE97508 dataset; (E) Venn diagrams of overlapping DEGs from two datasets and IRGs derived from the AmiGO2 and KEGG databases.

Validation of Hub IRG Expression by qRT–PCR

To verify whether the 14 differentially expressed IRGs obtained from the previous analyses were valid, we verified the mRNA expression levels of these genes in two different RB cell lines, the low-invasive RB cell line WERI-RB1 and the high-invasive RB cell line Y79. As shown in Figure 3, four of the eight most upregulated genes and three of the six most downregulated genes were positively validated by qRT–PCR (P < 0.05), which showed significant differences in the expression of these genes as the invasiveness of the cells increased (TYROBP: **P < 0.01; low density lipoprotein receptor-related protein 1 (LRP1): *P < 0.05; LCP1: **P < 0.01; collectin subfamily member 12 (COLEC12): **P < 0.01; LPCAT1: *P < 0.05; SH3GL2: ***P < 0.001; RAR related orphan receptor B: **P < 0.01). Although there was a significant difference in the expression level of carbonic anhydrase 2 (CA2) between the two RB cell lines, the experimental results were contrary to our previous microarray analysis results. All mRNA expression levels were normalized to beta-actin mRNA expression levels. Thus half of the hub IRGs were further confirmed by experimental verification, suggesting the greater plausibility of their roles and making them more reliable targets in tumor invasion. 

Figure 3.

 

Validation of differentially expressed IRGs in the low-invasive RB cell line WERI-RB1 and the high-invasive RB cell line Y79. (A) Validation of upregulated IRGs in two different kinds of RB cell lines by qRT-PCR; (B) Validation of downregulated IRGs in two different RB cell lines by qRT-PCR.

Profile of Immune Cell Subtype Distribution Patterns

Based on previous studies^19–21 and our aforementioned enrichment analyses, we speculated that in addition to the IRGs, immune cell infiltration and interaction might contribute to the acceleration tumor invasion. We used the ssGSEA algorithm to quantify mRNA profiling data and calculate the score of 28 infiltrating immune cells in each RB sample from two different datasets (GSE110811 and GSE97508), in which immune cell types were identified by labeling of specific marker proteins. The enrichment score in ssGSEA represents the relative abundance of each type of immune cell. A heatmap was constructed to visualize the relative abundance of 28 infiltrating immune cell populations (Fig. 4A). We observed similarities in the immune landscape profiles in different datasets. For example, central memory CD4 T cells (highlighted in the blue box) were the main infiltrating cells with antitumor activity, and some cells exhibiting protumor functions (e.g., macrophages, neutrophils, and regulatory T cells) had relatively low abundances. 

Figure 4.

 

Results of immune cell infiltration analysis by the ssGSEA algorithm. (A) Heat map of the 28 immune cell proportions and compositions in invasive and noninvasive RB tissues. (B) Quantified comparison of the distribution of 28 infiltrating immune cells in invasive and noninvasive RB tissues.

Additionally, we further quantitatively determined whether there were significant differences in the proportions of immune cells between the noninvasive and invasive groups. Our data showed that the proportions of MDSCs and NKT cells were significantly higher in invasive RB tissues than in noninvasive tissues, and these results were consistent in the GSE110811 and GSE97508 datasets (MDSCs—GSE97508: P = 0.039, GSE110811: P = 0.015; NKT cells—GSE97508: P = 0.015, GSE110811: P = 0.018) (highlighted in the green boxes, Fig. 4B). These findings suggested that MDSCs and NKT cells, especially protumor MDSCs, could have potential clinical significance in RB and are worthy of further investigation because of their function.²² 

Correlation Analysis between Key IRGs and Infiltrating Immune Cells

Subsequently, we further explored the relationship between immune cells and 7 key IRGs by Pearson correlation analysis to determine whether the experimentally verified immune genes accurately reflected the state of the tumor immune environment. Correlation analysis showed that only three of the seven IRGs (SH3GL2, LCP1, and TYROBP) were significantly correlated with at least one overlapping cell type in the two datasets (Fig. 5). SH3GL2 was negatively correlated with the overlapping immune cells, and both LCP1 and TYROBP were positively correlated with the overlapping immune cells. Of particular interest to us, all three genes were highly correlated with MDSCs, and only TYROBP was also correlated with NKT cells. Taken together, these results suggest that MDSCs coud play an important role in RB development. 

Figure 5.

 

Correlation between hub IRGs and infiltrating immune cells from two datasets. The size of the dots represented the strength of the correlation. The color of the dots represented the P value, and a darker purple color indicated a lower P value. Red labels denote consistent and significant correlations; based on the red label, two key immune cells were labeled in gold. *P < 0.05 was considered statistically significant.

SH3GL2 Expression iIs Low in Invasive Retinoblastoma Cells and Affects Cell Viability and Migration In Vitro

Because of the low mRNA expression of LCP1 and TYROBP in both WERI-Rb1 and Y79 cells (data not shown), we focused on SH3GL2, which requires further elucidation. To verify the expression of SH3GL2, we evaluated the protein level of SH3GL2 in low-invasive WERI-Rb1 and high-invasive Y79 cells by western blot analysis. Consistent with the mRNA results described in Figure 4B, SH3GL2 protein level was significantly higher in WERI-Rb1 cells than in Y79 cells (Fig. 6A). Statistical data are represented as histograms (relative fold change: WERI-Rb1 cells, 1 ± 0.074; Y79 cells, 0.202 ± 0.313; *P < 0.05) (Fig. 6B). 

Figure 6.

 

SH3GL2 expression is low in invasive retinoblastoma cells and affects cell viability and migration in vitro. (A, B) The protein level of SH3GL2 in WERI-Rb1 or Y79 cells was determined by Western blot assay, and the quantitative data were shown as histograms (n = 6; *P < 0.05). (C, D) WERI-Rb1 cells were transfected with si-NC or si-SH3GL2, and Y79 cells were transfected with OE-NC or OE-SH3GL2. After transfection, SH3GL2 protein level was detected by western blot analysis, and the relative protein levels are represented as histograms. (E) Viability of RB cells after transfection (n = 4; *P < 0.05, **P < 0.01). (F–I) Transwell assays represented the relative migration levels of different RB cell lines after transfection (n = 3; *P < 0.05, ***P < 0.001).

To confirm that SH3GL2 plays an important role in RB progression, WERI-Rb1 cells were treated with SH3GL2 siRNA, and SH3GL2 was exogenously expressed in Y79 cells. As shown in Figures 6C and 6D, SH3GL2 expression in WERI-Rb1 cells was significantly decreased in the si-SH3GL2 group compared with the control group (relative fold change: si-NC, 1 ± 0.143; si-SH3GL2, 0.576 ± 0.073; *P < 0.05), and its expression in Y79 cells was significantly increased in the OE-SH3GL2 group compared with the control group (relative fold change: OE-NC, 1 ± 0.280; OE-SH3GL2, 1.508 ± 0.262; *P < 0.05). Moreover, downregulation of SH3GL2 significantly increased the viability of WERI-Rb1 cells, and SH3GL2 upregulation decreased Y79 cell viability (relative fold change: WERI-Rb1 cells: si-NC, 1 ± 0.467; si-SH3GL2, 1.484 ± 0.673; *P < 0.05; Y79 cells: OE-NC, 1 ± 0.235; OE-SH3GL2, 0.827 ± 0.215; *P < 0.01) (Fig. 6E). Furthermore, Transwell assays showed that more WERI-Rb1 (si-SH3GL2 group) cells and fewer Y79 (OE-SH3GL2 group) cells than the corresponding control cells migrated from the upper chamber into the lower chamber, indicating that there may be an inverse relationship between SH3GL2 expression and RB cell invasion (relative fold change: WERI-Rb1 cells: si-NC, 1 ± 0.460; si-SH3GL2, 3.982 ± 1.499; *P < 0.05; Y79 cells: OE-NC, 1 ± 0.364; OE-SH3GL2, 0.343 ± 0.145; *P < 0.05) (Figs. 6F–I). Thus these data uniformly show that SH3GL2 might be involved in the malignant behavior of RB cells. 

Downregulation of SH3GL2 in Retinoblastoma Cells Promotes Tumor Growth, Invasiveness and Immunosuppression In Vivo

To confirm the in vitro results and the correlation between SH3GL2 expression and MDSCs, we established a xenograft mouse model by subcutaneously injecting WERI-RB1 cells. Although BALB/c nude mice are athymic, it has been reported that MDSCs promote tumor growth in immunodeficient nude mice through a process independent of T cells.^23–25 Si-NC or si-SH3GL2 was intratumorally injected once every three days for 12 days after tumor cell inoculation. As shown in Figures 7A and 7B, the tumor volume in the si-SH3GL2 group was smaller than that in the control at the beginning of the experiment (tumor volume: day 1, si-NC: 306.56 ± 89.96 mm³; si-SH3GL2: 218.75 ± 39.57 mm³), whereas they were significantly increased at days 9 and 12 (tumor volume: day 9, si-NC: 426.23 ± 113.28 mm³; si-SH3GL2: 699.02 ± 115.89 mm³; day 12, si-NC: 440.24 ± 161.15 mm³; si-SH3GL2: 818.01 ± 84.85 mm³; *P < 0.05). Accordingly, the tumor weight was also significantly increased by si-SH3GL2 injection (tumor weight: day 12, si-NC: 0.146 ± 0.072 g; si-SH3GL2: 0.283 ± 0.039 g; *P < 0.05) (Fig. 7C). Moreover, the protein level of SH3GL2 was significantly decreased in tumor tissues by si-SH3GL2 injection (relative fold change: si-NC: 0.932 ± 0.168; si-SH3GL2: 0.451 ± 0.084; *P < 0.05) (Figs. 7D, 7E). 

Figure 7.

 

Downregulation of SH3GL2 in retinoblastoma promotes tumor growth in vivo. (A) Representative macroscopic images of tumors after intratumoral injection of si-NC or si-SH3GL2 for 12 days. (B) The line chart of the tumor volume from xenotransplantation model (n = 5; *P < 0.05). (C) The tumor weight of the xenotransplantation model (n = 5; *P < 0.05). (D) The protein levels of SH3GL2 in tumors from the si-NC group and si-SH3GL2 group were represented by western blot assay. (E). Relative quantification of SH3GL2 protein level in tumors (n = 5; *P < 0.05).

To detect the tumor invasiveness, representative liver images from si-NC and si-SH3GL2 groups are shown in Figure 8A (black arrowhead: metastatic liver nodules), and H&E staining of liver tissues revealed that downregulation of SH3GL2 led to more metastatic lesions. Relative metastatic area in liver was quantified in Figure 8B (si-NC: 0.082 ± 0.076; si-SH3GL2: 0.364 ± 0.148; *P < 0.05). In addition, we identified the MMP9 in tumor tissues by immunofluorescence staining. MMP9 is a significant invasive marker in retinoblastoma and other types of tumors^26–28 (Fig. 8C). As Figure 8D shows, downregulation of SH3GL2 significantly increased the MMP9 expression (si-NC: 0.099 ± 0.025; si-SH3GL2: 0.168 ± 0.028; *P < 0.05). Another gene CXCR4, which was related to tumor progression and malignancy,^29,30 was also increased in si-SH3GL2 group compared with the negative control group (si-NC: 0.110 ± 0.032; si-SH3GL2: 0.202 ± 0.046; *P < 0.05) (Supplementary Fig. S1). 

Figure 8.

 

Downregulation of SH3GL2 in retinoblastoma promotes tumor invasivenesss in vivo. (A) Representative macroscopic images and H&E staining of livers from si-NC or si-SH3GL2 groups. Bar: 400 µm. Magnification bar: 100 µm. (B) Quantification data showed relative metastatic area in livers from si-NC and si-SH3GL2 groups (n = 5; *P < 0.05). (C) Immunofluorescence analysis of MMP9 (red) and DAPI (blue). Bar: 100 µm. (D) Quantification data showed MMP9/DAPI ratio from si-NC and si-SH3GL2 groups (n = 5; *P < 0.05).

Furthermore, we verified the infiltration of MDSCs in the tumor and spleen after SH3GL2 was downregulated. The tissues were double stained for Gr-1 and CD11b, which are MDSC-specific markers. As shown in Figures 9A and 9B, the number of MDSCs (orange triangles) was significantly increased in the si-SH3GL2 group compared with the negative control group (si-NC: 13.14 ± 2.98; si-SH3GL2: 19.78 ± 3.54; *P < 0.05). In particular, consistent with our hypothesis, MDSCs were abundant in spleen tissue and were significantly increased by more than twofold in the si-SH3GL2 group compared with the control group (si-NC: 20.14 ± 10.83, si-SH3GL2: 52.98 ± 30.10; *P < 0.05) (Figs. 9C, 9D). This was consistent with a previous study,³¹ demonstrating that MDSCs can accumulate in the spleen (the largest immune organ) to curb antitumor immune responses, which are involved in immune escape. Based on these results, SH3GL2 plays an important role in the tumor invasion and immune response, and MDSCs might be the critical immune cells involved in RB progression. 

Figure 9.

 

Downregulation of SH3GL2 in retinoblastoma increases MDSCs in tumor and spleen tissues. (A) Immunofluorescence analysis of both Gr-1 (green)- and CD11b (red)-positive cells showed that MDSCs in tumor tissues were increased in the si-SH3GL2-injected group. Bar: 100 µm. Magnified images from the merged figures (yellow box); the yellow triangle indicated the MDSCs in the tumor tissues. Bar: 25 µm. (B) Quantification of the number of Gr-1- and CD11b-positive cells/field in tumors (n = 5; *P < 0.05). (C) Immunofluorescence analysis of both Gr-1 (green)- and CD11b (red)-positive cells showed that the MDSCs in the spleen tissues were increased in the si-SH3GL2-injected group. Bar: 100 µm. Magnified images from the merged figures (yellow box); the yellow triangle showed the MDSCs in the spleen tissues. Bar: 25 µm. (D) Quantification of the Gr-1- and CD11b-positive cells/field in spleens (n = 5; *P < 0.05).

SH3GL2 Shows Low Expression in the Human Invasive Retinoblastoma

To explore and further validate the SH3GL2 expression in human samples, we collected six noninvasive and six invasive human retinoblastoma samples. SH3GL2 immunoreactive area showed an apparent decline in the invasive tumors (Fig. 10A). The SH3GL2/DAPI area ratio was significantly decreased in the invasive group (noninvasive: 0.028 ± 0.014; invasive: 0.005 ± 0.003; **P < 0.01) (Fig. 10B). Magnification images of yellow box area from Figure 10A are shown in the Figure 10C. This result was consistent with our cell and animal experiments, indicating that SH3GL2 could serve as an important gene in retinoblastoma metastasis and progression. 

Figure 10.

 

SH3GL2 shows low expression in the human invasive retinoblastoma. (A) Immunofluorescence analysis of SH3GL2 (red) and DAPI (blue) in human tumor tissue. Bar: 100 µm (B) Quantification data showed SH3GL2/DAPI ratio from noninvasive and invasive groups (n = 6; **P < 0.01). (C) Magnified images from the merged figures (yellow box). Bar: 50 µm.

Discussion

In the present study, we used two public datasets to analyze the immune genes and cells that could indicate the risk of progression from noninvasive RB to invasive RB. Functional enrichment analysis identified eight upregulated genes and six downregulated genes. Seven genes were significantly differentially expressed between the noninvasive and invasive RB cell lines. Moreover, MDSCs were identified as important immune cells in RB invasion. Correlation analysis of immune genes and cells revealed that three of the seven experimentally verified genes involved in RB progression were associated with MDSCs. Specifically, we found that SH3GL2 played a critical role in tumor progression both in vitro and in vivo and was closely related to MDSCs. 

A strength of our study was the use of bioinformatics analysis. Because of the cost and time-consuming nature of obtaining sufficient patient samples and because the traditional method may require many years of effort, it was not considered feasible to acquire a variety of patient tissues. Bioinformatics analysis is a relatively recently developed tool and has progressed very rapidly in the last few years. It has been well used in various fields, including molecular modeling, molecular interaction analysis, phylogenetic analyses, molecular dynamics simulation, and drug design.³² On the one hand, bioinformatics tools can be used to analyze existing online data to save considerable time in collecting patient tissues. On the other hand, bioinformatics analysis has made it possible to quickly identify potential hub genes from sets of thousands of genes in diseases, reducing costs by orders of magnitude. Last, it permits virtual testing of hypotheses and, hence, allows us to make a thorough and informed decision before launching costly experiments. Therefore we used bioinformatics analysis to obtain a novel understanding of immune genes and cells involved in RB invasion and progression. 

Through bioinformatics tools, our integrated analyses first indicated that SH3GL2 expression showed a statistically significant reduction in invasive RB compared to noninvasive RB (Fig. 2). SH3GL2 was identified as an SH3 domain–containing protein.³³ Our experimental data confirmed that SH3GL2 is related to RB development. SH3GL2 exhibited higher expression in the RB cell line with lower invasive capability (WERI-Rb1), whereas downregulation of SH3GL2 by siRNA notably increased RB cell viability and invasiveness (Fig. 6). In addition, animal data verified that decreased expression of SH3GL2 significantly facilitated tumor growth and liver metastasis (Figs. 7, 8). Previous studies have reported that SH3GL2 is mainly expressed in the central nervous system and is considered to be a tumor suppressor gene in many cancers, such as breast cancer, urothelial carcinoma, head and neck squamous cell carcinoma, and glioma.^34–37 The data of human samples also validated that SH3GL2 showed lower expression in invasive groups than that of noninvasive groups (Fig. 10). The consistency of the bioinformatic and experimental data indicates that the bioinformatics analysis had relatively good sensitivity and specificity. 

Moreover, our bioinformatics analysis indicated that the number of MDSCs was significantly elevated in aggressive RB, suggesting that MDSCs can promote RB invasion (Fig. 4). MDSCs, which are important cells in tumor progression, are a heterogeneous population of immunosuppressive cells derived from bone marrow and include early bone marrow progenitor cells, immature granulocytes, macrophages and dendritic cells at different stages of differentiation.²² The prominent feature of these cells is their low level in the peripheral blood of healthy people, whereas they proliferate rapidly under conditions such as inflammation or infection, especially after tumors develop. Kumar et al.³⁸ reported that MDSCs mainly inhibit the proliferation and cytotoxicity of antitumor T cells and block the activation of NK cells to favor tumor progression. Additionally, MDSCs affect tumor biological behaviors, such as tumor cell invasion, metastasis, and other nonimmune mechanisms.³⁹ However, the relationship between SH3GL2 and MDSCs has not been reported to date. Our correlation analysis data showed that SH3GL2 was inversely associated with MDSCs (Fig. 5), a finding that was supported by the in vivo experiments. We discovered that SH3GL2 silencing led to RB progression and an increase in MDSCs in vivo (Fig. 9). This evidence indicated that SH3GL2 might be a tumor suppressor in RB by regulating MDSCs. However, further investigation is required to reveal the underlying mechanism of SH3GL2 and MDSCs in RB progression. In addition, we verified that LCP1 and TYROBP expression was increased in the invasive tumors and that these genes had a strong association with MDSCs. 

Furthermore, despite the availability of multiple databases and various bioinformatics tools, there are still important limitations, such as the incompleteness of data platforms and the lower accuracy of this method compared with experimental methods. For instance, we verified that CA2 expression was elevated in invasive RB cells compared with noninvasive RB cells, a result that was completely different from our bioinformatics analysis results. CA2 has been reported to be an invasion-associated gene in bladder cancer.⁴⁰ We speculate that this discrepancy is attributed to the differences between the cell growth environment in vitro and the tumor microenvironment in vivo or to the lower accuracy of bioinformatics analysis. Therefore, from this perspective, bioinformatics analysis absolutely demands experimental verification to eliminate false data, and there is an urgent requirement for better tools and software. 

In summary, our data suggest a potential important effect of immune genes and cells on RB invasion and progression. We successfully identified the roles of the novel immune gene SH3GL2 and MDSCs, revealing a promising biomarker and therapeutic target for RB invasion. 

Acknowledgments

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

Supported by grants from the National Natural Science Foundation of China (grant no. 81900850, 82101134). 

Disclosure: S. Chen, None; X. Chen, None; P. Zhang, None; S. Chen, None; X. Wang, None; Q. Luo, None; Z. Cui, None; Y. Huang, None; L. Wan, None; X. Hou, None; H. Yao, None; X. Liu, None; A. He, None; Z. Jiang, None; J. Qiu, None; Y. Li, None; K. Yu, None; J. Zhuang, None