Regional and Temporal Differences in Gene Expression of LH_BETAT_AG Retinoblastoma Tumors

Abstract

Purpose.: The purpose of this study was to evaluate by microarray the hypothesis that LH_BETAT_AG retinoblastoma tumors exhibit regional and temporal variations in gene expression.

Methods.: LH_BETAT_AG mice aged 12, 16, and 20 weeks were euthanatized (n = 9). Specimens were taken from five tumor areas (apex, anterior lateral, center, base, and posterior lateral). Samples were hybridized to gene microarrays. The data were preprocessed and analyzed, and genes with a P < 0.01, according to the ANOVA models, and a log₂-fold change >2.5 were considered to be differentially expressed. Differentially expressed genes were analyzed for overlap with known networks by using pathway analysis tools.

Results.: There were significant temporal (P < 10⁻⁸) and regional differences in gene expression for LH_BETAT_AG retinoblastoma tumors. At P < 0.01 and log₂-fold change >2.5, there were significant changes in gene expression of 190 genes apically, 84 genes anterolaterally, 126 genes posteriorly, 56 genes centrally, and 134 genes at the base. Differentially expressed genes overlapped with known networks, with significant involvement in regulation of cellular proliferation and growth, response to oxygen levels and hypoxia, regulation of cellular processes, cellular signaling cascades, and angiogenesis.

Conclusions.: There are significant temporal and regional variations in the LH_BETAT_AG retinoblastoma model. Differentially expressed genes overlap with key pathways that may play pivotal roles in murine retinoblastoma development. These findings suggest the mechanisms involved in tumor growth and progression in murine retinoblastoma tumors and identify pathways for analysis at a functional level, to determine significance in human retinoblastoma. Microarray analysis of LH_BETAT_AG retinal tumors showed significant regional and temporal variations in gene expression, including dysregulation of genes involved in hypoxic responses and angiogenesis.

Retinoblastoma (RB) is the most common intraocular malignancy in children, affecting approximately 1 in 15,000, for an incidence of 250 to 300 new diagnoses a year in the United States. ^{1 –4} As treatment has progressed from external beam radiation therapy (EBRT) to chemoreduction combined with focal therapies, survival rates have climbed to 99% with a large percentage of children maintaining vision. ⁵ Despite the significant advancements in treatment and survival, current chemotherapy regimens and focal therapies may result in complications. Children are subjected to toxic chemotherapeutic drugs for multiple cycles, resulting in considerable risk for systemic toxicities. ^{6 –8} Focal consolidation therapies also contribute to morbidity, depending on the intraocular tumor size and location. ⁹ Finally, chemoreduction success varies depending on tumor classification, with more advanced eyes achieving tumor control in 47% to 83% of cases. ^10,11 As a result, a greater understanding of tumorigenesis is necessary to develop adjuvant therapies to potentially treat tumors that are unresponsive to current treatment protocols and to minimize local and systemic complications of treatment. 

The genetics of RB development have been studied, beginning with Knudson's “two hit” hypothesis. ¹² In RB development, mutations or epigenetic changes in both alleles of the RB1 gene lead to loss of retinoblastoma protein (pRB). pRB binds to E2F, which acts as a transcriptional regulator of the cell cycle. Loss of both RB1 alleles leads to susceptibility of retinal cells to formation of RB. It has been proposed that development of RB requires more than the two hits proposed by Knudson, with Corsen and Gallie ¹³ reviewing the literature for evidence of further genetic changes necessary for tumor development. 

The paradigm of cancer treatment and understanding has shifted from solely targeting hyperproliferative tumor cells and associated oncogenes/tumor suppressor genes to also targeting cancer stromal tissue, which consists of complex multicellular interactions, termed the tumor microenvironment. ^14,15 This environment consists of a plethora of cell types, including endothelial cells, fibroblasts, and inflammatory cells, ¹⁶ that, along with numerous growth factors and signaling molecules, contribute to tumorigenesis. Hypoxia has been strongly correlated with tumor growth, progression, resistance to therapy, and metastasis. ¹⁷ It has been shown that through O²-sensitive pathways, hypoxia alters tumor cell behavior, resulting in an integrated response of tumor cells to the tumor microenvironment, leading to altered gene expression and tumor adaptation and survival. Alterations to hypoxia include signaling through the mammalian target of rapamycin (mTOR), hypoxia inducible factor (HIF), and the unfolded protein response (UPR). These responses lead to altered cellular metabolism, angiogenesis, and other cell survival mechanisms. ¹⁸ The genetic changes associated with a tumor's adaptation to the microenvironment are prospective avenues for more specific and targeted therapy. In addition, an understanding of the timing of gene expression is fundamental in the optimal use of novel, multimodal adjuvant treatments. 

Hypoxia and angiogenesis have been shown to be significantly associated with tumor proliferation and metastasis. Our previous studies with the LH_BETAT_AG mouse model of RB have shown that 20% to 26% of tumoral areas are hypoxic and that there is a spatial difference, with hypoxia primarily in the central and basal regions of the tumor. ¹⁹ In addition, we have shown that there is a spatial distribution of blood vessel maturation, with mature blood vessels concentrated centrally, while immature neovessels radiate peripherally. As a result, we have proposed that RB tumors grow radially from the center, with the apex, anterior, and posterior margins serving as the leading edges. ²⁰  

We hypothesize that there are regional and temporal differences in gene expression of murine RB tumors corresponding to areas of tumor heterogeneity and variations in tumor microenvironment. These variations may provide valuable information regarding specific genes and pathways that facilitate tumor growth, progression, and resistance to treatment. The purpose of our study was to evaluate these regional and temporal differences by using mRNA microarray analysis in a transgenic model of RB. 

Methods

LH_BETAT_AG Mouse Model for RB

The study protocol was approved by the University of Miami Institutional Animal Care and Use Review Board Committee. The LH_BETAT_AG transgenic mouse model used in this study has been characterized previously. ²¹ This animal model develops bilateral multifocal retinal tumors that are stable and grow at a predictable rate (i.e., tumor at 4 weeks is grossly undetectable, at 8 weeks is small, at 12 weeks is medium, and at 16 weeks is large), with histopathologic, immunopathologic, and ultrastructural features that resemble human RB tumors. The LH_BETAT_AG mouse model has contributed to elucidating mechanisms in tumor development and progression and to providing a platform for the development of adjuvant therapies. 

Molecular Genomic Array Analysis in LH_BETAT_AG Retinal Tumor Growth

Transgenic mice with documented intraocular tumors were killed and the eyes enucleated at 12, 16, and 20 weeks of age (n = 9; three at each time point). These time points were chosen to represent early to advanced tumors and to correspond to time points studied in this model regarding angiogenesis, hypoxia, gelatinase expression, and tumor response to therapy. Five 3.37-mm³ sections were obtained from each tumor (apex, anterior lateral, center, base, and posterior lateral). Samples were meticulously dissected under a microscope by an experienced handler. With the time points chosen (12, 16, and 20 weeks), tumors had already grown to a macroscopic size, allowing an experienced handler to dissect them without obtaining normal retina, as the tumor has already expanded into the globe with boundaries distinct from normal retina. Although there is always a possibility of contamination, dissection protocols were meticulously used in ascertaining the extraction of sufficient and appropriate samples. Samples were obtained from the five areas based on the proposed mechanism of tumor growth and progression, with radial growth from the center. Leading edges have been shown to consist of more immature vasculature, whereas, central areas consisted of mature vasculature as well as a higher percentage of hypoxia. ^19,22 All specimens (n = 45) were placed in a lysis solution and stored at −20°C until analyzed. The samples were hybridized (16 hours) to a unique gene microarray chip that provides whole gene expression data (no 3′ bias) for over 28,000 genes (GeneChip Mouse Gene ST 1.0 arrays; Affymetrix, Santa Clara, CA). We used the Robust Multichip Average (RMA) Express method (http://rmaexpress.bmbolstad.com/ written by Ben Bolstad, University of California, Berkeley, and provided in the public domain) to measure differential gene and probe level expression measures (log₂), with a false-discovery rate (FDR) set at 5%. Quality control plots and summary measures were generated with R/Bioconductor 2.9.10. ^23,24 Gene level measures were analyzed by using analysis of variance (ANOVA) models for repeated measures, considering temporal or regional effects, using custom scripts (written for SAS ver. 9.2; SAS, Cary, NC). Genes with P < 0.01 from the ANOVA models and a log₂-fold change >2.5 were considered to be differentially expressed. Differentially expressed genes were analyzed for overlap with known networks, by using pathway analysis tools (GeneGo; St. Joseph, MI). Network significance was evaluated on the basis of the size of the intersection between our list of significantly differentially expressed genes and the set of genes/proteins corresponding to a network module/pathway. Each network was associated with a z-score that ranked the networks according to saturation with the objects from the experimental gene list. The z-score ranked the networks of the analyzed network algorithm with regard to their saturation with genes from the experiment. A high z-score means that the network is highly saturated with genes from the particular experiment. Each network was also associated with a g-score, which modifies the z-score on the basis of the number of canonical pathways used to build the network. If a network has a high g-score, it is saturated with expressed genes (from the z-score), and it contains many canonical pathways. A P value was determined by comparing the observed amount of intersection with the amount expected under the null hypothesis that the amount of overlap follows a hypergeometric distribution. 

Results

A total of 28,000 genes were assessed over the three different time points in the five regions. Significant temporal differences in gene expression were found between 12-, 16-, and 20-week LH_BETAT_AG RB tumors (P < 10⁻⁸, two-way analysis of variance, ANOVA). In addition, analysis identified genes with a greater than 2.5-fold difference in expression between the three time points that varied depending on region. There were significant differences in gene expression across time for 190 genes apically, 84 genes anterolaterally, 126 genes posteriorly, 56 genes centrally, and 134 genes at the base (Tables 1A, 1B). 

Table 1.

 

Differential Gene Expression by More Than 2.5-Fold in Five Tumor Regions with the Corresponding Ratio of Change

A. Upregulated Genes
Base		Center		Apex		Anterior-lateral		Posterior-lateral
Gene	Ratio of Change	Gene	Ratio of Change	Gene	Ratio of Change	Gene	Ratio of Change	Gene	Ratio of Change
KRT15	5.1448	GM9912	2.6366	CENPK	4.0012	KRT16	3.802	No gene name	3.6121
KRT4	5.0582			2810417H13R1	3.8653	No gene name	3.1226	GEN1	3.3636
CYP4A12B	4.6897			CASC5	3.7474	SPINK5	3.0778	CCNB1	3.2991
KRT13	4.5248			BUB1	3.7178	XIST	3.0497	HMMR	3.2884
KRT14	4.4163			RRM2	3.6037	KRT14	2.9956	CCNB1	3.2238
ANXA8	4.3149			CDH9	3.5187	S100A9	2.9409	EG665955	3.2117
KRT5	4.301			MKI67	3.5062	No gene name	2.9264	CCNB1	3.1662
TACSTD2	4.182			PBK	3.4692	KRT5	2.8168	NCAPG2	3.0262
KRT6B	4.1757			KRT4	3.466	BC100530	2.7341	LCN2	3.0002
CEACAM1	4.1426			MCM6	3.4537	SNORD116	2.7022	LPHN2	2.951
9930032O22RI	4.1396			KRT14	3.4365	No gene name	2.702	C79407	2.9377
UPK1B	4.0946			NUDT10	3.4343	DSG3	2.7003	CRYBA1	2.9014
DSG3	3.956			GEN1	3.3905	SNORD116	2.6892	PTPN3	2.8936
ENPP3	3.933			CCNB1	3.388	MMP3	2.6707	CDH9	2.8539
KRT6A	3.9257			FAM111A	3.3757	SNORD116	2.6558	SHCBP1	2.852
SLC6A14	3.8995			CCNB1	3.3655	S100A8	2.6524	BUB1	2.804
CBR2	3.8739			KIF11	3.3435	CRYBA1	2.628	PBK	2.7594
TMPRSS11B	3.7424			CCNB1	3.3391	SNORD116	2.5994	6720489N17RI	2.7205
DSC3	3.7068			KRT6B	3.3143	GM9912	2.5633	GINS1	2.7187
TRIM29	3.7054			KRT15	3.3032	KRT6B	2.5615	CASC5	2.6874
SCEL	3.6572			TOP2A	3.2889	SCG3	2.5175	TTK	2.6598
TMPRSS11A	3.6417			CENPH	3.2249	No gene name	2.5092	SERPINA3N	2.6336
ADH7	3.6356			ANLN	3.2177	SNORD116	2.5034	TPX2	2.629
ALDH3A1	3.6319			HIST1H2AB	3.1439			LPHN2	2.6035
AKR1B7	3.5832			BC100530	3.1268			SGOL2	2.5912
BC100530	3.5756			C330027C09RI	3.0901			KIF4	2.585
DSG1A	3.5619			SGOL2	3.074			C330027C09R1	2.5798
DSC2	3.5213			KRT5	3.0701			PPIL5	2.5518
PGLYRP1	3.4777			CKS2	3.0693			KIF14	2.5517
FABP6	3.4499			MNS1	3.0536			CENPK	2.5336
PLAC8	3.4191			KRT13	3.0486			RRM2	2.5299
DSP	3.3436			ACTC1	3.0294			CHEK1	2.528
S100A14	3.2942			KIF4	3.0204			SNORD116	2.5238
ADH1	3.2901			NCAPG	3.0026			KIF2C	2.5138
MUC4	3.2862			SPC24	3.0014			CLSPN	2.5114
RBP2	3.2579			HELLS	2.9981			NUF2	2.5102
CALML3	3.2514			NUF2	2.9795
CAPG	3.2201			KRT16	2.965			KIF14	2.5517
KRT19	3.1758			DSCC1	2.9644			CENPK	2.5336
LUM	3.1756			ARHGAP11A	2.9578			RRM2	2.5299
CMAH	3.1643			SMC2	2.9405			CHEK1	2.528
SERPINB5	3.1593			C79407	2.9212			SNORD116	2.5238
EMP1	3.1549			PRR11	2.8937			KIF2C	2.5138
PPL	3.1142			NDC80	2.8826			CLSPN	2.5114
PYHIN1	3.09			KRT6A	2.8734			NUF2	2.5102
IL1F9	3.0649			BIRC5	2.8706			KIF2C	2.5138
DEFB1	3.0552			LPHN2	2.8672			CLSPN	2.5114
SAMD9L	3			GINS1	2.8565			NUF2	2.5102
FGFBP1	2.9933			DSG3	2.8543			CLSPN	2.5114
ANXA2	2.9889			2810417H13RI	2.8429			NUF2	2.5102
9930023K05RI	2.9783			SHCBP1	2.8422			CLSPN	2.5114
IFITM1	2.9492			CKS2	2.8379			NUF2	2.5102
GSTA3	2.9068			CDC2A	2.8331			CLSPN	2.5114
GPR110	2.9005			LPHN2	2.8233			NUF2	2.5102
DCN	2.847			MASTL	2.8215
SPINK5	2.823			TPX2	2.8144
EHF	2.7891			KIF23	2.8003
KRT7	2.7817			GM9912	2.8001
1600029D21RI	2.769			TTK	2.7778
MALL	2.736			CKS2	2.7633
S100A9	2.7304			FBXO5	2.7586
OCM	2.7168			EG665955	2.7554
CLDN7	2.6877			PEG10	2.7433
				SPINK5	2.7412
GM9573	2.6769			NCAPH	2.7266
PERP	2.6515			KRT24	2.7084
ESRP1	2.6479			ECT2	2.695
GSTO1	2.6415			PPIL5	2.6805
GDA	2.6149			RBBP8	2.6706
PTGR1	2.6138			ANXA1	2.67
CYP2F2	2.6113			CENPF	2.6546
IFITM1	2.5923			CKAP2L	2.6376
C130090K23RI	2.5153			CYP4A12B	2.6161
GALNT3	2.5097			NCAPG2	2.5946
GM11428	2.5001			DLGAP5	2.5916
				LYZ2	2.588
				PRIM1	2.5849
				DBF4	2.5832
				LPHN2	2.5785
				ATAD2	2.5765
				TRIM59	2.5672
				MCM3	2.5618
				HMMR	2.5427
				NFIA	2.5368
				E2F8	2.534
				EXO1	2.5174
				2610039C10RI	2.5138

B. Downregulated Genes
Base		Center		Apex		Anterior-lateral		Posterior-lateral
Gene	Ratio of Change	Gene	Ratio of Change	Gene	Ratio of Change	Gene	Ratio of Change	Gene	Ratio of Change
CNGA1	−4.4701	GUCY2F	−4.1373	PDE6B	−5.0352	CNGA1	−4.6794	OPN1SW	−5.0585
RP1	−4.4608	PDE6B	−4.0316	CNGA1	−5.0191	RCVRN	−4.4971	GUCY2F	−4.6778
RCVRN	−4.3932	RCVRN	−3.8163	RCVRN	−4.8462	PDE6B	−4.4931	RCVRN	−4.5381
PDE6B	−4.3062	RHO	−3.7673	RPE65	−4.8052	RHO	−4.0053	CNGA1	−4.4204
GUCY2F	−3.9616	2610034M16R1	−3.7485	PDE6A	−4.7065	RHO	−3.9785	PDE6B	−4.3056
PDE6A	−3.9498	RP1	−3.5564	RP1	−4.6844	RP1	−3.9687	No gene name	−4.1752
RS1	−3.8567	PDE6A	−3.5066	RHO	−4.6441	PDE6A	−3.8735	RP1	−4.0818
RHO	−3.79	CNGA1	−3.4836	GUCY2F	−4.6367	GNAT1	−3.7975	2610034M16RI	−4.0801
TULP1	−3.6546	REEP6	−3.4382	RDH12	−4.5129	RS1	−3.6846	PDE6A	−4.0469
GRK1	−3.5238	RS1	−3.3294	GNAT1	−4.4516	SAG	−3.6444	RHO	−3.9724
RDH12	−3.4964	SLC24A1	−3.2645	VTN	−4.3343	GUCY2F	−3.638	GRK1	−3.8474
OPN1SW	−3.4265	RPE65	−3.2339	GUCA1A	−4.2548	VTN	−3.5997	RHO	−3.7115
GUCA1B	−3.4165	ARR3	−3.2278	SAG	−4.1852	SLC24A1	−3.449	CDS1	−3.6255
ADAMTS3	−3.4098	VTN	−3.1777	RHO	−4.1269	2610034M16RI	−3.4238	RS1	−3.5685
2610034M16RI	−3.3635	TULP1	−3.0543	TTR	−4.0968	IMPGI	−3.4029	RDH12	−3.5593
IMPG2	−3.3173	GUCA1B	−3.041	TULP1	−4.0863	REEP6	−3.4	GUCA1B	−3.4517
GM626	−3.3102	IMPG1	−3.0188	GUCA1B	−4.013	GRK1	−3.1713	IMPG1	−3.4482
PRPH2	−3.2885	PDE6G	−2.9703	PDC	−4.0107	GUCA1B	−3.1536	ARR3	−3.4129
SAG	−3.2643	OPNIMW	−2.9658	2610034M16RI	−4.0039	GUCA1A	−3.1167	FABP12	−3.3855
CALB1	−3.2542	CNGB1	−2.9629	SLC24AI	−3.9651	PDC	−3.0939	C030002C11RI	−3.3297
IMPG1	−3.2524	GRK1	−2.9333	IMPG1	−3.9048	TULP1	−3.0184	SLC24A1	−3.3201
PDE6G	−3.231	ME1	−2.9195	RGR	−3.8508	PRPH2	−3.0106	REEP6	−3.3172
WDR17	−3.2248	CALB1	−2.8971	REEP6	−3.8455	CNGB1	−2.9766	VTN	−3.3057
GNAT1	−3.2196	NRL	−2.8711	GM626	−3.8261	OPN1MW	−2.9231	SAG	−3.2897
RHO	−3.1948	SAG	−2.8696	RS1	−3.7915	PDE6G	−2.8083	ADAMTS3	−3.1062
SLC24A1	−3.188	ABCA4	−2.8646	OPN1SW	−3.7847	RPE65	−2.7849	WDR17	−3.0986
VTN	−3.1816	PRPH2	−2.8531	IMPG2	−3.7759	ABCA4	−2.7098	TULP1	−3.0714
FABP12	−3.1415	C030002C11RI	−2.8371	GRK1	−3.7481	NR2E3	−2.6847	UPK1B	−3.0451
ARR3	−3.0916	RDH12	−2.835	CDS1	−3.6674	FABP12	−2.6803	CNGB1	−3.0442
PDC	−3.0772	GLBIL2	−2.8234	PRPH2	−3.6533	PROM1	−2.652	OPN1MW	−3.0287
ROM1	−2.9229	PEX5L	−2.8231	ROM1	−3.6231	IMPG2	−2.6508	FAM161A	−3.0234
FAM161A	−2.9159	WDR78	−2.8094	FABP12	−3.6175	GM626	−2.6288	GLB1L2	−3.0131
A930003A15RI	−2.9133	NXNL1	−2.8054	PDE6G	−3.4549	GM10664	−2.61	PDC	−2.9972
WDR78	−2.9098	WDR17	−2.7996	ABCA4	−3.4477	NRL	−2.6025	TAC1	−2.9925
GM10664	−2.8845	NR2E3	−2.7562	PROM1	−3.3796	WDR17	−2.5942	ABCA4	−2.988
SAMD7	−2.8688	RGR	−2.7358	WDR78	−3.3315	GLB1L2	−2.5708	LRIT2	−2.9737
ADAMTS3	−2.8323	PROM1	−2.7204	CALB2	−3.2487	RDH12	−2.5568	PROM1	−2.9684
OPN1MW	−2.8125	SPATA1	−2.7178	GLB1L2	−3.2272	C530030P08RI	−2.5517	PDE6G	−2.9669
REEP6	−2.7962	RHO	−2.7129	NR2E3	−3.2055	WDR78	−2.5363	KRT12	−2.9552
LRIT2	−2.7883	PLA2R1	−2.7107	PPARGC1A	−3.2005	No gene name	−2.5165	PRPH2	−2.9497
GUCA1A	−2.7699	GNAT1	−2.701	CNGB1	−3.1234	LRIT2	−2.5045	WDR78	−2.9474
ABCA4	−2.7201	CDS1	−2.6796	MPP4	−3.104			GNAT1	−2.9465
C030002C11RI	−2.7019	FAM161A	−2.6649	WDR17	−3.1008			PSCA	−2.946
CNGB1	−2.6999	FABP12	−2.6637	NRN1	−3.0902			RPE65	−2.9378
LRRC2	−2.6992	LRIT2	−2.6099	NRL	−3.0787			PLA2R1	−2.9029
SPATA1	−2.6869	ME1	−2.5699	CALB1	−3.077			CDR2	−2.9013
GM626	−2.665	A930003A15RI	−2.5299	PEX5L	−3.0711			NRL	−2.9007
HIST2H3C2	−2.6462	3632451O06RI	−2.5294	GM11744	−3.0574			MPP4	−2.8886
A330023F24RI	−2.6434	BST1	−2.5082	LGI1	−2.993			GM626	−2.874
PEX5L	−2.6154	GM11744	−2.5007	FAM161A	−2.9842			GM11744	−2.8575
PDE6C	−2.6118			GM10664	−2.9635			RABGEF1	−2.8558
TAC1	−2.5092			CHRNB3	−2.9436			TDRD7	−2.8429
				SLC17A6	−2.9268			PEX5L	−2.8263
				OPN1MW	−2.9101			UNC13C	−2.7611
				LRIT2	−2.9017			NR2E3	−2.7492
				COX8B	−2.8988			APBH	−2.7348
				ENPP2	−2.8935			HERC3	−2.7337
				ARR3	−2.8684			PDE6C	−2.7175
				HERC3	−2.8628			GP2	−2.7166
				C030002C11RI	−2.8626			4833423E24RI	−2.7154
				GM626	−2.8602			DCN	−2.7131
				CD59A	−2.8517			CENPH	2.7061
				SPATA1	−2.809			RGR	−2.7039
				ALDOC	−2.8089			CALML3	−2.6967
				LDHB	−2.7821			STFA3	−2.684
				SAMD7	−2.7803			GUCA1A	−2.6827
				GM626	−2.7474			PPEF2	−2.6757
				UNC13B	−2.7459			NT5E	−2.667
				RLBP1	−2.726			NXNL1	−2.6568
				RGS9	−2.7235			IMPG2	−2.6363
				RP1L1	−2.7206			MLANA	−2.6293
				GNB1	−2.7097			PGLYRP1	−2.6228
				PLA2G5	−2.7038			C130021120RI	−2.619
				PLA2R1	−2.674			GM626	−2.6093
				ME1	−2.664			ATP8A2	−2.6083
				KRT18	−2.6561			S100A14	−2.6036
				TAC1	−2.6465			GM10664	−2.6028
				CLIC6	−2.6461			CALB2	−2.6025
				PPEF2	−2.644			CHRNB3	−2.5758
				ATP8A2	−2.6388			TYRP1	−2.5606
				HK2	−2.638			LRAT	−2.5599
				LRRC67	−2.6246			ADH1	−2.5317
				GRIA3	−2.6234			ROM1	−2.5257
				ZDHHC2	−2.5984			CSMD3	−2.5134
				LYNX1	−2.5834			PVALB	−2.5085
				PVALB	−2.5725
				GM626	−2.561
				NXNL1	−2.5568
				AIPL1	−2.5277
				LRAT	−2.5118
				ME1	−2.5094
				CDR2	−2.5021

To identify functional activity of the unique genes, we performed pathway analysis (GeneGo software). Tables 2 to 6 show the key networks involved in cellular growth and proliferation, hypoxia, cell signaling, and angiogenesis for the five regions with key objects, processes, pathways, and statistical significance. 

The basal regions (Table 2) of the RB tumors had significant differences (P = 1.54 × 10⁻¹¹ to 3.72 × 10⁻⁶¹) in regulation of Ras/G-protein signaling, regulation of cellular proliferation, mTOR/PKC signaling, and JAK/STAT signaling. The z-scores ranged from 13.56 to 48.21, and the g-scores ranged from 28.06 to 146.23. On analysis of networks, key elements mediating cellular proliferation included c-Myc, JAK/STAT, TGF-β, MDM2, RB protein, AMPK, let-7a microRNA, and p53. We also found cellular responses mediated by the PI3K/akt/mTOR pathway, as well as HIF, VEGF, NOTCH, and IGF-1 (Fig. 1). 

Table 2.

 

Key Networks with Associated Key Objects and Processes for Basal Tumor Regions

Top Networks	Key Network Objects	GO Processes	Total Nodes	Root Nodes	Pathways	P	z-Score	g-Score
3	STAT1, c-Raf-1, DNM1L (DRP1), NDPK B, 14-3-3 zeta/delta	Regulation of cell proliferation, Positive regulation of cellular process, regulation of apoptosis	50	18	23	3.46 × 10−22	22.44	51.19
4	mTOR, PKC, Tuberin, DLL1, NEURL1	Enzyme linked receptor protein signaling pathway, insulin-like growth factor receptor signaling pathway, cellular response to insulin stimulus	50	11	30	1.54 × 10−11	13.56	51.06

Figure 1.

 

Network in basal aspects of tumors involved in hypoxic responses (network 4; P < 10⁻¹¹, 30 pathways), with key elements including mammalian target of rapamycin (mTOR), phosphoinositide-3 kinase (PI3K), Akt or protein ki), hypoxia-inducible factor (HIF), vascular endothelial growth factor (VEGF), NOTCH, and insulin-like growth factor (IGF-1). Cyan lines: fragments of canonical pathways; red circles: upregulated genes; blue circles: downregulated genes.

The central regions (Table 3) of the tumors had significant differences (P = 4.44 × 10⁻⁷ to 3.22 × 10⁻⁶¹) in regulation of cellular proliferation, response to oxygen levels, and regulation of oxidoreductase activity, JAK/STAT signaling, and regulation of cellular metabolism. The z-scores ranged from 8.81 to 45.46, and the g-scores ranged from 33.60 to 133.81. Network analysis showed that cellular proliferation and regulation were mediated by TNF-α, ERK1/2, RB protein, MAPK, VEGF, ubiquitin, TGF-β, and MDM2. Tumor response to oxygen levels was mediated by MMP-2 and -9, TGF-β, E-cadherin, c-Myc, PAI, VEGF, and CDK2/CDK4. In addition, the cellular proliferation and signaling pathways were mediated by JAK/STAT, ERK1/2, cyclin D1, and 14-3-3. 

Table 3.

 

Key Networks with Associated Key Objects and Processes for Central Tumor Regions

Top Networks	Key Network Objects	GO Processes	Total Nodes	Root Nodes	Pathways	P	z-Score	g-Score
2	c-Raf-1, NCOA3, (pCIP/SRC3), H-Ras, NCOA2 (GRIP1/TIF2), TFF1	Regulation of cell proliferation, positive regulation of cellular process	51	11	85	1.05 × 10−10	12.29	118.54
5	SMAD2, NF-κB, SMURF2, SerRS, RPS25	Anatomic structure morphogenesis, response to oxygen levels	50	14	35	5.58 × 10−15	16.04	59.79
9	FOXO3A, 14-3-3 zeta/delta, Cyclin D2, GPIAP1, BTF	JAK-STAT cascade involved in growth hormone signaling pathway, regulation of cell proliferation	50	20	12	4.22 × 10−25	24.06	39.06

The apical regions (Table 4) of the tumors had significant differences (P = 1.91 × 10⁻⁷ to 1.29 × 10⁻⁸³) in PKA/G-protein signaling, regulation of cell proliferation, oxygen transport, and regulation of stress responses. The z-scores ranged from 9.40 to 57.83, and the g-scores ranged from 32.46 to 184.4. Tumor anatomic and morphogenic responses were found to be mediated by key elements, including VEGF, MMP-2, MMP-9, and TGF-β (Fig. 2, left). Cellular proliferation was found to be mediated by JAK/STAT, ERK1/2, c-Myc, NF-κB, IGF-1, TNF-α, and caspase-8. Finally, cellular responses to stress were found to be mediated by HIF-1, NF-κB, and c-Myc (Fig. 2, right). 

Table 4.

 

Key Networks with Associated Key Objects and Processes for Apical Tumor Regions

Top Networks	Key Network Objects	GO Processes	Total Nodes	Root Nodes	Pathways	P	z-Score	g-Score
1	TGF-β1, RUNX3, TGF-β receptor type II, GAS1, FALZ	Organ development, anatomic structure morphogenesis, system development	50	8	140	1.91 × 10−7	9.40	184.4
4	NF-κB, FCGRT, Cyclin D2, ISG20, SPT2	Positive regulation of cellular process, regulation of cell proliferation	50	14	27	1.21 × 10−15	17.05	50.80
9	NF-κB, COX17, Nectin-2, PSMD7, LOC283412	Regulation of oxidoreductase activity, positive regulation of defense response, positive regulation of oxidoreductase activity, regulation of response to stress	50	27	4	1.34 × 10−40	35.94	40.94

Figure 2.

 

(A) Network involved in angiogenesis (network 1), in apical leading edges of tumors (P < 10⁻⁷, 140 pathways). Key elements include vascular endothelial growth factor (VEGF), matrix metalloproteinase-9 (MMP-9), transforming growth factor-β (TGF-β). (B) Network involved in hypoxia (network 9) in the apical leading edges of tumors (P < 10⁻⁴⁰, 4 pathways). Response to hypoxia mediated by hypoxia-inducible factor (HIF).

The anterior-lateral regions (Table 5) of the tumors had significant differences (P = 1.41 × 10⁻¹⁵ to 2.26 × 10⁻⁸⁶) in cellular proliferation, cytokine-mediated signaling pathways, leukocyte migration, and glycosaminoglycan biosynthetic processes. The z-scores ranged from 16.94 to 58.76, and the g-scores ranged from 34.63 to 58.76. Tumor signaling cascades and growth regulation were found to be mediated by TGF-β, MMP-2, MMP-9, fibronectin, c-Myc, and NF-κB. Cytokine-mediated signaling was found to involve TNF-α, NF-κB, ICAM, E-selectin, and CCL-2/CCL-3. 

Table 5.

 

Key Networks with Associated Key Objects and Processes for Anterolateral Tumor Regions

Top Networks	Key Network Objects	GO Processes	Total Nodes	Root Nodes	Pathways	P	z-Score	g-Score
5	SMURF2, CRELD2, eIF4A1, SARA, SAHH2	Response to organic cyclic substance, transforming growth factor beta receptor signaling pathway	50	14	24	1.41 × 10−15	16.94	46.94
8	CCL13, HMG1 (amphotercin), NDPK B, HRPT2, Symplekin	Translational elongation, leukocyte migration, cytokine-mediated signaling pathway, response to mechanical stimulus, translation	50	21	9	9.13 × 10−28	26.67	37.92

Finally, the posterior-lateral regions (Table 6) of the RB tumors had significant differences (P = 1.07 × 10⁻⁸ to 1.27 × 10⁻⁷²) in regulation of cellular proliferation, JAK/STAT signaling, and cytokine-mediated signaling. Network analysis showed anatomic morphogenesis and growth regulation to be mediated by TGF-β, MMP-2, MMP-9, ERK-1/2, MAPK, Akt, and VEGF. Cell proliferation and regulation of cellular processes were found to be mediated by TNF-α, ERK-1/2, VEGF, EGFR, MAPK, PKC, IGF-1, and COX-2. 

Table 6.

 

Key Networks with Associated Key Objects and Processes for Posterolateral Tumor Regions

Top Networks	Key Network Objects	GO Processes	Total Nodes	Root Nodes	Pathways	P	z-Score	g-Score
1	SMAD1, SMAD2, TGFb 1, SCAP, RPS27A	Anatomic structure morphogenesis, organ development, regulation of developmental process	50	11	67	2.42 × 10−11	13.26	97.01
3	SHP-2, Ubiquitin, IGF-2, NIX, GlyRS	Regulation of cell proliferation, regulation of cellular process	50	9	38	1.07 × 10−8	10.70	58.20

Discussion

In the present study, we identified genes that are differentially expressed in the LH_BETAT_AG murine model of RB in five tumor regions (base, center, apex, anterior-lateral margin, posterior-lateral margin), and at three time points in tumor development (12, 16, and 20 weeks). Overall, gene expression was shown to significantly differ temporally (P < 10⁻⁸), as well as regionally. Of the 28,000 gene probe sets analyzed, we found differential expression of 190 genes apically, 84 genes anterolaterally, 126 genes posteriorly, 56 genes centrally, and 134 genes at the base. Analysis showed that these dysregulated genes were associated with multiple networks and canonical pathways, including regulation of cellular proliferation, cellular signaling and stress responses, response to hypoxia, angiogenesis, as well as anatomic morphogenesis and growth regulation. 

As tumors grow, the proliferating cells experience an imbalance of oxygen metabolism, leading to a disorganized and irregular microvasculature network. As a result, there is reduced oxygen delivery, leading to a microenvironment with low oxygen partial pressure. ^{15,25 –27} In the harsh tumor microenvironment depleted of oxygen and nutrients, cells undergo a hypoxic and/or angiogenic switch to support further growth. Growth of immature neovessels is stimulated, along with cell adaptation to hypoxia, including increase in glycolysis for energy metabolism. ²⁶ These events represent a key transition in tumorigenesis as cells adapt, altering the gene expression necessary to drive further tumor proliferation. 

We have previously shown that the tumor vasculature is highly heterogenous, with higher concentrations of large, mature vessels toward the base/center and smaller neovessels radiating into the periphery. These observations in both human and the LH_BETAT_AG model suggest that RB tumors proliferate radially from the center. ^20,22 The present study is the first to show the regional variation in gene expression by microarray in the LH_BETAT_AG RB model. Gene expression was found to be dysregulated in a regional distribution, with the apex and base showing the most variation. We hypothesized that the apical areas correspond with the leading edges of the proliferating tumor consisting of highly metabolic, hyperproliferating cells, and neovessels dependent on growth factor support. Concurrently, the basal regions correspond to areas of hypoxic stress necessitating cellular adaptation. The current observations support our previous model of RB progression and blood vessel development, identifying potential canonical pathways that may mediate these responses in the transgenic model of RB tumors. ^{27 –29}  

In the apical regions, more advanced tumors had a significantly different expression of genes involved in cellular proliferation, hypoxia, angiogenesis, and regulation of stress responses. TGF-β, NF-κB, PTEN, and cyclin D2, as well as signaling through the JAK/STAT pathway were shown to be key mediators. TNF-α and caspase-8, which were found to be dysregulated in RB cell lines ³⁰ and human RBs, ³¹ respectively, were also found in the present study to be dysregulated in advanced LH_BETAT_AG tumors. Hypoxic responses, including angiogenesis, were also identified in the apical regions, areas of intense growth composed of hyperproliferating cells competing for oxygen and nutrients. These adaptations were mediated through hypoxia-inducible factor (HIF), vascular endothelial growth factor (VEGF), and matrix metalloproteinases (MMP-2 and -9), factors known to be significantly associated with advanced tumors with high degrees of hypoxia. 

The urokinase plasminogen activator (uPA) and the receptor (uPAR) may play an integral role in tumor proliferation and metastasis. The uPA system consists of serine proteases that lead to activation of plasmin, which in turn activates matrix metalloproteinases (MMPs). MMPs have been linked to tumor growth and metastasis for their role in the degradation of the extracellular matrix. ³² Recent studies have identified hypoxic elements within the genes that regulate the uPA system of several other tumors. ^33,34 The present study has identified dysregulated genes involved in pathways for the uPA system as well as angiogenic pathways that use MMPs. These findings support our prior work showing enhanced tumor control with reduced expression of MMPs using anecortave acetate. ³⁵ MMPs appear to play a key role in the tumor microenvironment, and studies are needed to further elucidate the mechanisms and effects of tumor treatment. 

As neoplastic cells proliferate, they experience a highly anabolic state requiring altered function to provide sufficient energy and waste removal, thus preventing cell death signals. ^{36 –38} In response to increasing metabolic demands, as well as altered microenvironments inside and outside the cells, neoplastic cells alter cellular metabolism, adopting a metabolic phenotype through differential gene expression for key enzymes and regulators of cellular metabolism. ³⁹ Neoplastic cells preferentially use glycolysis, referred to as the Warburg effect, in normoxic and hypoxic conditions, rather than oxidative phosphorylation and its higher ATP yield. In basal regions, areas of RB tumors shown to consist of a high population of hypoxic cells, we found differential expression of genes involved in networks and pathways that are upstream regulators of tumor metabolism, including PI3K/Akt/mTOR, IGF-1 signaling, and AMPK. In addition, our study showed significant involvement of HIF and Akt, both of which have been shown to increase GLUT1, a glucose transport receptor. ³⁷ Our prior studies in the LH_BETAT_AG RB tumor model have demonstrated the efficacy of targeting cellular metabolism in enhancing tumor control. Using 2-DG, a glycolytic inhibitor, we have shown that tumor burden is significantly reduced when treated with both systemic and local subconjunctival delivery of 2-DG. We have also recently shown that local delivery of 2-DG, when combined with chemotherapy, further enhances tumor control over either treatment alone. ¹⁹  

In addition, similar to prior human RB microarray gene expression studies, the PI3K/Akt/mTOR pathway was found to be dysregulated, potentially implicating mTOR as a therapeutic target. ⁴⁰ We found significant differential expression of this network in basal regions of advanced tumors, regions shown to experience significant hypoxia compared with other regions. The mammalian target of rapamycin (mTOR) is a serine-threonine kinase that is composed of two multiprotein complexes. Activation of these complexes leads to phosphorylation of downstream effectors, leading to regulation of protein translation, cell growth, proliferation, and metabolism. ⁴¹ Hypoxia has been shown to be a negative regulator of mTOR complexes, thus potentially acting as an inhibitor of growth and progression. However, it has been proposed that hypoxia drives the mutations necessary to deregulate mTOR signaling, termed hypoxic tolerance. As a result, it may be important for tumor cells to retain control of mTOR signaling for continued growth and proliferation. ¹⁸ mTOR inhibitors have been investigated in other tumors, and further studies are needed to define the role of mTOR and inhibitors of this pathway in RB tumor control. Notably, early studies with focal delivery of rapamycin in the LH_BETAT_AG RB model show promise, as the mTOR inhibitor led to enhanced tumor control (Murray TG, et al. IOVS 2010;43:ARVO E-Abstract 2067). 

Our current and previous findings support an evolving and complex cellular metabolic phenotype that differs in early versus advanced tumors, as well as variations secondary to the local, heterogeneous tumor microenvironment. This study provides further evidence of dysregulation of pathways involved in cellular metabolism, including genes to assess at a functional level to determine the effect on tumor growth and development. 

In addition, the study is the first to show temporal variation in gene expression with advanced tumors compared to earlier tumors examined with microarray analysis. In transgenic RB tumors, differential gene expression and associated pathways suggest potential mechanisms of tumorigenesis. These specific genes and pathways require further functional analysis to determine whether direct targeting has an effect on tumor growth and development, or whether the changes seen are indirect effects of other tumorigenic processes and would not serve as useful targets for tumor control. As with other tumors, we propose that regions of RB tumors have differential gene expression and may respond differently to various treatments, depending on the tumor's age, tumor burden, location, hypoxia, vasculature, and cellular metabolism. Therefore, the present study provides evidence of regional and temporal tumor heterogeneity, emphasizing the potential importance of timing in gene expression in the development of optimally timed, multimodal treatments with agents that target tumor cells, as well as components of the tumor microenvironment, including angiogenesis, hypoxia, and cellular metabolism. We propose that future treatment of RB tumors must target not only proliferating tumor cells, but also the local tumor microenvironment, in a multimodal, optimally timed, local approach. 

As functional analyses are performed on specific genes and pathways, including studies on targeting key elements, it is important to identify alterations in gene expression after treatment as well as to identify key escape pathways that cells use during times of cellular stress. These escape pathways may prove to be important regarding resistance and adaptation of tumor cells to current and future targeted therapies. 

Limitations of the present study include the small sample size of three mice at each time point (12, 16, and 20-weeks), as well as for each different region at each of these time points. In addition, we investigated the effects of time and location in the LH_BETAT_AG murine model of RB, which has been shown to share many similarities with human RB, but the correlation in gene expression between human and mouse tumors has not been fully determined. As a result, before the current findings can be related to human RB, further functional studies are needed of transgenic RB tumors and human RB cell lines. 

In conclusion, the findings in our study have shown that gene expression in the LH_BETAT_AG model of RB has significant temporal variations. In addition, we have shown significant regional variations, with the apical and basal regions exhibiting dysregulation of genes in key pathways involved in the regulation of cellular proliferation and growth, response to oxygen levels and hypoxia, and regulation of cellular processes, cellular signaling cascades, and angiogenesis. Prior studies have defined the gene expression profiles of RB tumors ^40,42 ; our study of the LH_BETAT_AG murine model suggests that these gene expression profiles may be dynamic, varying in a temporal and regionally dependent fashion. We anticipate that future developments in ocular oncology will focus on locally delivered, optimally timed therapies that target tumor cells as well as the tumor microenvironment. 

Footnotes

 Supported by National Institutes of Health center Grants R01 EY013629, R01 EY12651, and P30 EY014801; the American Cancer Society, Sylvester Comprehensive Cancer Center; and an unrestricted grant to the University of Miami from Research to Prevent Blindness, Inc.

Footnotes

 Disclosure: S.K. Houston, None; Y. Pina, None; J. Clarke, None; T. Koru-Sengul, None; W.K. Scott, None; L. Nathanson, None; A.C. Schefler, None; T.G. Murray, None