RESEARCH ARTICLE

Open Access

<figure><img src="image_1.png" /></figure>

CrossMark

miRNome landscape analysis reveals a 30
miRNA core in retinoblastoma

Blanca Elena Castro-Magdonel¹,³, Manuela Orjuela², Javier Camacho³, Adda Jeanette García-Chéquer¹, Lourdes Cabrera-Muñoz⁴, Stanislaw Sadowinski-Pine⁴, Noé Durán-Figueroa⁵, María de Jesús Orozco-Romero⁵, Ana Claudia Velázquez-Wong⁶, Adriana Hernández-Ángeles¹, Claudia Hernández-Galván⁷, Citlali Lara-Molina⁸ and M. Verónica Ponce-Castañeda¹*

**Abstract**

**Background:** miRNAs exert their effect through a negative regulatory mechanism silencing expression upon hybridizing to their target mRNA, and have a prominent position in the control of many cellular processes including carcinogenesis. Previous miRNA studies on retinoblastoma (Rb) have been limited to specific miRNAs reported in other tumors or to medium density arrays. Here we report expression analysis of the whole miRNome on 12 retinoblastoma tumor samples using a high throughput microarray platform including 2578 mature miRNAs.

**Methods:** Twelve retinoblastoma tumor samples were analyzed using an Affymetrix platform including 2578 mature miRNAs. We applied RMA analysis to normalize raw data, obtained categorical data from detection call values, and also used signal intensity derived expression data. We used Diana-Tools-microT-CDS to find miRNA targets and ChromDraw to map miRNAs in chromosomes.

**Results:** We discovered a core-cluster of 30 miRNAs that were highly expressed in all the cases and a cluster of 993 miRNAs that were uniformly absent in all cases. Another 1022 miRNA were variably present in the samples reflecting heterogeneity between tumors. We explored mRNA targets, pathways and biological processes affected by some of these miRNAs. We propose that the core-cluster of 30 miRs represent miRNA machinery common to all Rb, and affecting most pathways considered hallmarks of cancer. In this core, we identified miR-3613 as a potential and critical down regulatory hub, because it is highly expressed in all the samples and its potential mRNA targets include at least 36 tumor suppressor genes, including RB1. In the variably expressed miRNA, 36 were differentially expressed between males and females. Some of the potential pathways targeted by these 36 miRNAs were associated with hormonal production.

**Conclusion:** These findings indicate that Rb tumor samples share a common miRNA expression profile regardless of tumor heterogeneity, and shed light on potential novel therapeutic targets such as mir-3613 This is the first work to delineate the miRNA landscape in retinoblastoma tumor samples using an unbiased approach.

**Keywords:** miRNome, Retinoblastoma, mir-3613, Tumor heterogeneity

<footnote>* Correspondence: vponce@ifc.unam.mx
¹Medical Research Unit in Infectious Diseases, Hospital de Pediatría, CMN SXXI, Instituto Mexicano del Seguro Social, Av. Cuauhtémoc 330, 06720 Mexico City, Mexico
Full list of author information is available at the end of the article</footnote>
## Background

MicroRNAs (miRNAs) are key biologic regulators, structurally they are small non-coding RNA sequences (22–25 nucleotides) that are complementary to specific 3'-UTR mRNAs. Their function once hybridized to their target is to silence the target mRNA through cleavage of the molecule or inhibiting translation. These negative regulatory mechanisms gives miRNAs a prominent position in the control of many cellular processes, including carcinogenesis; miRNAs are estimated to regulate posttranscriptional expression of approximately 70% of human genes. High throughput miRNAs expression profiles are more sensitive for discriminating different human tissues including different types of cancer than mRNA profiles. MiRNAs expression profiles also change with progressive stages of tumor development, and have been proposed as potential useful biomarkers for cancer. These changes in miRNA profiles have also been proposed as potential markers that can provide an opportunity and serve to dissect and improve our understanding of cellular functions and gene networks involved in cancer related clinical processes such as response to treatment [1–3].

Retinoblastoma (Rb) is an intraocular malignancy of early childhood and is a robust model of heritable predisposition to developing cancer [4]. Even though retinoblastoma is a rare tumor, its study has led to the understanding of critical molecular mechanisms in cancer development. Several miRNAs have been studied in Rb, including some studied initially in other malignancies. One such example is oncomir 1 also known as cluster-17-92, a potent oncogenic cluster previously studied in B-cell lymphomas, lung carcinoma, and prostate cancer, and found to be expressed in Rb [2, 5–7]. Studies in retinoblastoma cell lines examined miRNAs involved in specific process such as hypoxia in Rb, and proposed miR-181b, miR-30c-2, miR-125-3p, miR-497 and miR-491-3p as hypoxia regulated miRNAs [8]. MiRNA profiles obtained with low and medium density arrays have also been used to compare Rb tumor with normal retina in order to find differentially expressed miRNAs and to identify deregulated pathways between tissues [9–12].

While approximately 40 miRNAs have been described in Rb, reports have not examined whether reported miRNAs are uniformly expressed in all or most cases [13]. This is relevant because tumors are not homogeneous; they display both intra and inter tumor heterogeneity [14, 15]. While examination of tumor heterogeneity has been directed largely towards intra tumoral heterogeneity, inter tumoral heterogeneity has scarcely been addressed. Analysis of complete intra tumoral heterogeneity and subclonal architecture from human primary tumors is

particularly challenging in retinoblastoma, because of limited tumor volume and because diagnostic constraints preclude sampling of multiple tumor sites. Nonetheless expression profiles provide the opportunity to improve our understanding of miRNAs in Rb by permitting inter tumoral heterogeneity to be taken in to account. In this context, we can address unanswered questions such as to what extent are tumors similar and different between patients in terms of miRNA expression? or what is the miRNA machinery that defines Rb and is shared between tumors from different patients? And which miRNAs are similarly or differentially expressed in different patients?

To address these questions, we determined the miRNA profiles of tumor tissues from 12 children with Rb using a high density microarray platform which included 2578 mature miRNAs. Using the detection call algorithm as an ‘expressed’ or ‘non-expressed’ dichotomized score for each miRNA in the array, we first generated a simplified miRNome landscape map of expressed and non-expressed miRNA in Rb. We identified 142 miRNAs present in all samples, and within these a central core of 30 miRNAs that were uniformly highly expressed. We also identified 993 miRNAs that were uniformly absent in all samples and 1022 miRNAs that were present in only some of the samples. We explored the targets and gene networks affected by miRNAs that were uniformly or variably expressed, using unsupervised clustering tools combined with clinical descriptors in order to examine the potential significance of the identified clusters. With this unbiased approach, we found miRNAs that had not previously been recognized to be involved in Rb. Our findings also reveal the magnitude and complexity of tumor heterogeneity in terms of miRNA expression and highlight some of the biological processes that are negatively regulated by these miRNAs.

## Methods

### Primary cultures

Tumor samples were obtained from 12 patients who were treated with enucleation prior to receiving any adjuvant radiation or chemotherapy. These were cultured in RPMI medium supplemented with 12% FBS for a week, using growth conditions reported for the retinoblastoma cell line Y79. All samples were collected after obtaining parental written informed consent, for participation in a larger IRB approved case-series study [16, 17] involving patients at the Hospital de Pediatría, Instituto Mexicano del Seguro Social (IMSS) and Hospital Infantil de México Federico Gómez Secretaría de Salud in Mexico City. The patient corresponding to T1 had a positive family history for Rb, while all other patients had sporadic Rb.
## RNA isolation

Total RNA was isolated from the primary cultures using TRI-zol® (Invitrogen, CA, USA) reagent and was then solubilized in nuclease free water. RNA was quantified using a Nanodrop (Thermo Fisher Scientific, USA) spectrometer and its quality was confirmed by electrophoresis using 1.5% agarose stained with SYBR® safe (Thermo Fisher Scientific, USA). We divided RNA in 200 ng aliquots and stored at 70 °C until use. We amplified miR-16 using stem-loop specific miRCURY LNA™ oligonucleotides in several randomly chosen samples to ensure that the RNA isolated was informative.

## Labeling and tailing samples for microarrays

Total RNA was tailed and biotinylated using Affymetrix Flash-tag biotin for miRNAs microarray (Affymetrix, USA) and spike-in control probes were added according to manufacturer instructions. Briefly, poly-A tailing on the 3' end was carried out at 37 °C in a 15 µl reaction volume containing 1× reaction buffer, MnCl₂ 25 mM, 1 µl of 1:50 ATP mix and 1 µl of phosphatidic acid phosphatase (PAP) polymerase. The biotin was incorporated to these poly-A tails at 25 °C after adding 4 µl of ligation mix containing biotin and 2 µl of T4 DNA ligase, yielding a total of 21 µl. An Enzyme Linked Oligosorbent assay (ELOSA) was performed to confirm biotin labeling.

## miRNAs microarrays

Hybridization cocktails were added to labeled samples containing 2X hybridization mix, 27% formamide, DMSO, 20× of hybridization control probes and nuclear free water. Each sample was injected in a GeneChip® miRNA 4.0 array and hybridized at 48 °C for 16–18 h. Every sample was washed in the GeneChip® Fluidics Station 450 following FS450–0002 protocol and array fluorescence was measured by the GeneChip® Scanner 3000 7G.

## Data analysis

Raw data from CEL files were analyzed using Affymetrix® Expression ConsoleTM software to normalize fluorescence signals and to obtain quality control data. To adjust background/signal a Robust Multichip Analysis (RMA) was carried out and once completed we observed that the values of spike-in control probes complied with the quality standards established by the manufacturer. ‘.chp’ files were created and then used for further expression analysis, the names of ‘.cel’ files (raw data) and ‘.chp’ files corresponding to each sample are in the Additional data Additional file 1. Using the ‘.chp’ files we obtained two types of data for each probe in the array: expression log2 intensity signal data, and detection call data categorizing each miRNA as ‘Present’ (P) or ‘Absent’ (A), along with

the corresponding *p* value estimates. The statistical algorithms used for the detection call metrics were developed by the supplier using an experimental design called the Latin Square [18, 19], on which naturally absent transcripts were spiked in a complex background at known concentrations. This analysis generates a detection *p*-value to determine the detection call, indicating whether a transcript is reliably detected (Present) or not detected (Absent) [20, 21].

## Bioinformatics analysis

Both categorical and expression data were analyzed. For this, we employed text format tables containing data from all hybridized samples. All data were analyzed using Multi-experiment-viewer (MEV) TM4 Microarray Software Suite [22]. Non-supervised and supervised analyses were performed such as hierarchical cluster analysis (HC) and significance analysis for microarray (SAM) which incorporates correction for multiple testing [23]. For the ‘one class SAM’ the response variable or outcome was a quantitative continuous variable (expression data), used to find the miRNAs that were most highly expressed across all the samples. For the two class SAM the outcome variable was a quantitative continuous variable (also expression data) for two groups, used to find miRNAs that were differentially expressed between male and female patients. Fig-tree was used to visualize dendrograms derived from HC analysis. DIANA-Tools-microT-CDS, miRPathv3, (at http://www.microrna.gr/), Target-Scan and Tarbase v.7 (at http://diana.imis.athena-innovation.gr/Diana-Tools/index.php?r=tarbase/index) were used to search for mRNA targets and signaling pathways affected by miRNAs of interest [24]. In addition, we used TargetScan (at http://www.targetscan.org/vert_71/) and miRDB (at http://www.mirdb.org/miRDB/) to search for mRNA targets [25, 26]. The TS gene database was used to search for tumor suppressor genes (at https://bioinfo.uth.edu/TSGene/) [27, 28]. In order to find the chromosomal location for miRNAs, we used annotation data from the miRBase [26] (at ftp://mirbase.org/pub/mirbase/CURRENT/genomes/ and at http://www.mirdb.org/miRDB/). ChromDraw was used to map all miRNAs classified as “Absent” in order to compare their location with respect to chromosomal regions that we had previously reported as regions of recurring loss in Rb [29].

## Validation by qRT-PCR

We used qRT-PCR with locked nucleic acid primers to validate our findings; miRNA-16 and 3613-3p were used to validate ‘present’ in all tumors, and 3613-5p and 4529 were used to validate ‘absent’ in all the samples. We prepared cDNA with specific primers for these miRNAs
and prepared qRT-PCR mixes using 1 µg total RNA tumor with Mercury LNA™ Universal RT microRNAs PCR Exiqon kit. Once cDNA was synthesized, each sample was run in duplicate using Illumina Eco Real-Time PCR System and results were plotted to compare cycle detection.

## Results

### Discretized data reveals miRNA landscape

To characterize the Rb miRNome we used the GeneChip miRNA 4.0 arrays from Affymetrix including probes for 2578 mature human miRNAs, and RNA from Rb primary cultures obtained from 12 tumors as described above. All children received enucleation as their first line of therapy. Corresponding clinical data is summarized in Table 1.

In order to obtain an initial panoramic view of how many miRNAs are expressed and how many are not expressed in these tumors, we used detection call metrics based on statistical criteria as an approximation for ‘expressed’ or ‘not-expressed’ categories for each miRNA (as detailed above).

With this discretized data we then used an unsupervised hierarchical clustering analysis, and produced a heat map using yellow for miRNAs classified as ‘expressed’ and black for those miRNAs classified as ‘not expressed’. The resulting heat map (Fig. 1) shows a miRNA landscape composed of three main clusters. The cluster at the bottom named “P” (for present) shows 561 miRNAs expressed in almost 90% of samples including a remarkable smaller cluster with 142 miRNAs that are expressed in all the samples. The second cluster at the center of the heat map named “A” (for absent) shows 995 miRNAs that are uniformly not expressed in all samples and a third cluster named “V” (for variably present) with 1022 miRNAs that are expressed in some

<table><caption>**Table 1** Clinical data from patients with Rb</caption><thead><tr><th>Patient</th><th>Age at diagnosis (months)</th><th>Gender</th><th>Laterality<sup>a</sup></th><th>Clinical stage<sup>c</sup></th></tr></thead><tbody><tr><td>1</td><td>26<sup>b</sup></td><td>F</td><td>B</td><td>II</td></tr><tr><td>2</td><td>60</td><td>M</td><td>U</td><td>II</td></tr><tr><td>3</td><td>37</td><td>F</td><td>U</td><td>II</td></tr><tr><td>4</td><td>52</td><td>F</td><td>U</td><td>II</td></tr><tr><td>5</td><td>9</td><td>M</td><td>B</td><td>II</td></tr><tr><td>6</td><td>29</td><td>M</td><td>B</td><td>II</td></tr><tr><td>7</td><td>19</td><td>M</td><td>U</td><td>II</td></tr><tr><td>8</td><td>12</td><td>M</td><td>B</td><td>IV</td></tr><tr><td>9</td><td>32</td><td>M</td><td>U</td><td>II</td></tr><tr><td>10</td><td>33</td><td>M</td><td>U</td><td>II</td></tr><tr><td>11</td><td>36</td><td>M</td><td>U</td><td>I</td></tr><tr><td>12</td><td>36</td><td>F</td><td>U</td><td>III</td></tr></tbody></table>

<figure><img src="image_2.png" /></figure>

but not all the samples. The list of miRNAs in these clusters is shown in the Additional file 2. We validated miR-16 as present and miR-4529-3p as absent using qRT-PCR as part of our initial approach to describe the general miRNOME landscape of miRNAs in the on/off state. Results of this validation are shown in Additional file 2: Figure S1.

### Functional core of 30 miRNA in Rb discovered

We subsequently focused on the 142 miRNAs present in all samples and asked if the level of expression among them could tell us which miRNAs are most abundant and thus more functionally relevant. To answer this we extracted the corresponding expression level data into a heat map (Fig. 2a) and also plotted the signal intensity median for each of the 142 miRNAs across the samples in a histogram (Fig. 2b). We applied to this group of miRNAs a Significant Analysis for Microarrays (SAM) [23] to determine the most highly expressed miRNAs among them. In this SAM analysis the response (outcome) variable was a quantitative continuous variable (expression data) in one group for one class SAM, corresponding to the analysis results shown in Fig. 2 which shows the most significant highly expressed miRNAs from the 142 miRNAs that were found to be present in all tumor samples. Using these criteria we found 30 out of 142 miRNAs to be very highly expressed (red bars in Fig. 2c), and 22 miRNAs out of 30 were top ranked (see

<footnote><sup>a</sup>(B bilateral, U unilateral); <sup>b</sup>Positive family history; <sup>c</sup>St Jude’s staging system</footnote>
<figure><img src="image_3.png" /></figure>

Fig. 2b). Some of the miRNAs plotted as top ranking were not significant by SAM since expression levels among samples were more variable despite similarly high medians.

We next investigated the potential mRNA targets, pathways and biological processes affected by this core of 30 miRNAs. Using the MicroT-CDS algorithm to predict miRNA targets from the Diana-Tools platform, we found 8120 potential target genes and according to miR-Pathv3 [24], 182 mRNAs in this group of predicted targets are related to cancer (Additional file 3). The function of miRNAs is to repress protein translation, thus we searched for tumor suppressor genes in the list of potential targets and found that 48 of these mRNA targets are tumor suppressors (TS) according to the TS gene database [27, 28]. The complete list can be found in Additional file 3. From the core of 30 miRNAs, we selected the five miRNAs with the greatest number of tumor suppressor gene targets: miR-3613-3p with 36 targets, miR-4668-5p with eight, miR-5787 with seven, miR-762 with two and miR-1273 g-3p with four targets. We noted that miR-3613-3p not only has 36 tumor suppressor targets, but also, using TarBase v.7 we were able to identify 400 mRNAs for which there is experimental evidence supporting them as targets for this miRNA. Using other algorithms for theoretical target predictions such as TargetScan and miRDB, we uncovered that miR-3613-3p has several thousand potential targets.

Once we had explored potential targets, we then focused on identifying pathways and biological processes

affected by this core of 30 miRNAs. We used DIANA-mirPATH v3.0 and found 59 pathways that are potentially regulated by this core of 30 miRNAs (*p* < 0.05). Many of these pathways are related to proliferation, angiogenesis and apoptosis, critical processes belonging to the so called-cancer hallmarks [30] (Fig. 3). From this list of 59 pathways, we found 12 pathways (*p* < 0.001) related to anti proliferation, regulation of pluripotency in stem cells, pathways in cancer, axon guidance and others (Additional file 3).

We confirmed miRNA-3613-3p expression by qRT-PCR (Additional file 2: Figure S2) and also searched in independent public cancer datasets generated with miRNA Affymetrix arrays v.3 and v.4 in samples from patients with hepatoblastoma [31], medulloepithelioma [32], Ewing's sarcoma (mostly pediatric malignant tumors) [33] and synovial sarcoma [34]. As shown in Additional file 2: Figure S3, median and mean expression levels of miR-3613-3p are consistently higher than 7, and are as high as 10 in medulloepithelioma than its counterpart miR-3613-5p in all data sets. Signal distribution in all data sets were similar and show the range of expression level across the miRNOME.

## Variable miRNAs are related to gender

To test if clusters 'P' and 'V' identified in the miRNome landscape were related to known clinical properties corresponding to these tumors we first focused on the 419 miRNAs from the cluster 'P', present in most but not all samples from the landscape map (Fig. 1). We extracted
<figure><img src="image_4.png" /></figure>

the corresponding expression data, performed unsupervised hierarchical clustering and superimposed clinical data (laterality, age at diagnosis, clinical stage, or gender) on the two branch tree we obtained. Surprisingly, cases appeared to group largely according to gender in the two branches with the exception of one case T5 (Fig. 4a and b). In order to identify miRNAs most likely to be related to gender differences, using SAM, we searched among the 419 in cluster P for differentially expressed miRNAs. The response (outcome) variable was again a quantitative continuous variable (expression data) on two groups, using two class SAM, corresponding to the differential analysis. Results shown in Fig. 4c show the genes that are most significantly differentially expressed when comparing males and females. This analysis yielded 36 miRNAs that are differentially expressed between males and females (Fig. 4c), the list of the miRNAs and the intensity plots can be found in Additional file 3.

We next searched for potential mRNA targets and pathways affected by these 36 miRNAs, and identified 2900 potential targets and 50 significant (*p* < 0.001) pathways. From these, eight correspond to hormonally and developmentally related processes (Fig. 4c).

A similar analysis was performed on the 1022 miRNAs contained in cluster V from the miRNome landscape. We obtained a three branched cluster using expression data and hierarchical clustering analysis. Subsequent search for mRNA targets led us to uncover pathways

related to pluripotency, stem cells, angiogenesis and migration. Although these pathways suggest a potential relationship to characteristics that might be found in more invasive tumors, we were unable to find correlations with microscopic features in the tumors or other clinical descriptors. The three branched figure and the list of pathways can be found in Additional file 3.

## Some miRNAs are consistently absent in all Rbs

Because many chromosomal regions have recurrently been described as lost in Rb, we next examined whether the 995 miRNAs in cluster 'A' that were classified as absent or undetectable in all samples, were located within loci or regions recurrently described as lost [35, 36]. Using annotated data from the miRBase we were able to map in the human genome 2573 out of 2578 miRNAs included in the Affymetrix chip. With a map that we had previously generated using NGS which demonstrated regions of recurrent losses in Rb [29], we identified in the 995 undetectable miRNAs, those located at loci recurrently lost in Rb by chromosome. In total, we found 144 miRNAs located in areas recurrently lost in Rb [37, 38] the complete list is shown in Additional file 3. Data plotted in Fig. 5a show a survey of miRNAs across the human genome, with three bars per chromosome representing the total number of miRNAs mapped to each chromosome, the corresponding undetectable miRs and the undetectable miRs mapped at recurrently lost regions per chromosome.
<figure><img src="image_5.png" /></figure>

Figure 5b shows an ideogram of chromosome 1's cyto- genetic map as an example, indicating the location of undetected miRNAs. Above the ideogram, data plots of gains and losses demonstrate that these miRNAs are in regions corresponding to areas of recurring chromo- somal loss in Rb.

### Analysis of cluster17–92 in Rb

Cluster 17–92 which is probably the best studied miRNA family in cancer, has been reported as highly expressed in Rb [29]. We explored the expression of cluster 17–92 within the more global perspective of inter-tumoral hetero- geneity presented in our approach. For this we first evalu- ated the detected/undetected score for the three human paralogs of this cluster [39] which are located in chromos- omes 7, 13 and X and applied unsupervised hierarchical clustering (Fig. 6a). To evaluate the expression levels of all members of this miRNA family in our 12 samples, we dis- aggregated each paralog member by chromosome cluster and plotted the median intensity signal for each of the samples (Fig. 6b). Four miRNAs from the paralog cluster

in chromosome 13 which includes mir-17-5p, mir-20a-3p, miR19b-3p and miR-92a-3p, have levels of expression greater than 6, while members of the paralog in chromo- some 7 are the most consistently expressed as a group.

### Discussion

With our initial step of using call detection algorithms, we obtained a simplified miRNOMIC landscape that is easier to interpret, and is composed of three categories: a central miRNA machinery shared by all the samples; an important number of miRNAs in the “off” state also shared by all the samples; and a third group of miRNAs that are variably present and that we believe correspond to biologic variability and tumor heterogeneity. We de- scribed a miRNOME landscape interpreting and situat- ing each miRNA present in the array, while not limiting this description to any biased grouping of miRNAs. Our findings are the product of data exploration with higher order bioinformatics analytical tools guided by clear and pre-defined questions. This approach allowed us to obtain a truly panoramic view of the miRNA landscape
<figure><img src="image_6.png" /></figure>

in Rb. These results are biological relevant despite the
small size of our sample. With this approach we were
able to “dissect” complex data using dichotomized data
to first determine how many and which miRNAs are
turned ‘on’ or ‘off’ in each tumor. Subsequent analysis of
expression patterns allowed us to interpret the potential
biologic significance of groups of miRNA with similar
levels of expression. We propose that this approach and
‘landscaping’ can provide a useful panoramic view of the
Rb miRNome, revealing a shared core machinery of 142
miRNAs expressed in all Rbs we studied, and also
uncovers the magnitude of inter tumoral heterogeneity
at the miRNA level. Our interpretation is that by using a
2578 miRNA array based tool, Rb may be defined by 142
miRNAs that are consistently detectable and 995
miRNAs that are consistently not detectable. Biologic

variability and heterogeneity might be defined by a set of
419 miRNAs that are detectable in most but not all the
cases and by 1022 miRNAs that are detectable only in
some cases.

Regarding the core of 30 miRNAs, and considering
that 1) arbitrary units of intensity in the plotted data run
from 0 to 14, 2) median signal intensity for all the mi-
croarrays is 2.01 and 3) that the average signal intensity
for the 30 miRNAs core is above 8, we interpret that
these highly expressed miRNAs constitute a functional
core of miRNAs shared by all tumors included in this
study. Notably most of these 30 miRNAs have been
described recently by NGS in cervical cancer [40]. The
majority of predicted and experimentally supported tar-
get genes identified for the complete 30 miR core are
targets for miR-3613-3p. Within these targets we found
<figure><img src="image_7.png" /></figure>

that 48 mRNA are tumor suppressors [27, 28], of which 36 are targets for miR-3613-3p. High expression of this miRNA in all samples, indicates that many more tumor suppressor genes beyond RB1 are likely down regulated or not expressed through miR-3613-3p's effect in Rb, suggesting that this miRNA may be a strategic down regulatory hub or node for this tumor.

The biological processes affected by this core of 30 miRNAs that are shared by all 12 cases in our study, are related to the malignant phenotype [41]. Therefore, we propose that this miRNA core belongs to a stable and likely critical genetic network that may characterize the identity profile of Rb.

Regarding validation miRNA expression values generated on different platforms cannot be directly compared

because unique labeling methods and probe sequences result in variable signal distributions for probes that hybridize to the same miRNAs [42], nonetheless, we were able to confirm miRNA-3613-3p expression using qRT-PCR. Results for miRNA 3613-5p and 4529 as non-detectable miRNAs in all the samples are in accordance with our microarray results, however it was challenging to verify by this technique, consistent high expression for miRNAs 3613-3p and 638 from the 30 miRNA core present in all samples by microarrays (Additional file 2: Figure S2). In amplification curves using Illumina Eco Real-Time PCR system it can be observed that not all the samples were consistently amplified, and signal detection occurred at later cycles of the run (Additional file 2: Figure S2a and b). In contrast qRT-PCR amplification
curves using Roche Light Cycler 480 system show more consistent amplification of miRNA 3613-3p in red and no amplification for miRNA 4529 in green (Additional file 2: Figure S2c). Quantitative agreement between results obtained with the two different amplification systems is thus moderate and concordance with microarray results is not robust. As an alternative for verification we searched in public datasets for experimental information that could be useful as an independent form of validation. We searched among datasets generated with the miRNA Affymetrix platform in primarily pediatric cancers [31–34], and found high expression of miRNA 3613-3p in all the datasets we consulted. In some datasets we found even higher expression levels of miR-3613-3p than in Rb, suggesting that this mir may have a relevant role not only in retinoblastoma. Although more robust methods for validation and experimental work are needed to confirm the relevance of mir-3613-3p in cancer cells, in our search for targets and pathways, exclusion of this miRNA resulted in a complete absence of pathways with the word ‘cancer’. Similarly, there was a complete absence of the 36 tumor suppressors that are targets linked to this miRNA, suggesting that presence of miR-3613-3p is critical to tumor suppression in Rb. High expression of mir-3613-3p could also potentially explain those rare Rb cases in which no oncogenic mutations have been found in RB1 [43]. From these observations we predict that by turning off mir-3613-3p, the inhibition of protein translation of any of the tumor suppressor genes silenced by this miRNA in Rb cells could be lifted, restoring function; this could potentially reverse the malignant phenotype or result in tumor cell death. The conspicuous presence of mir-3613-3p in all samples and its centrality shown by our analysis, suggests that we may have uncovered a so called ‘oncogene addiction’ state in these samples, in which there is a dependency on one or a few genes for maintenance of the malignant phenotype [44].

Regarding heterogeneity within the group of miRNAs present in most but not all the cases (419 miRNAs in cluster “P”), our analysis yielded a two branched tree containing two clusters that appeared closely associated with the gender of patients. The high intensity signals of these two groups of differentially expressed miRNAs, suggest that this represents an underlying relevant difference. Our finding of associated pathways related to development and hormone metabolism further supports a relationship with gender. These results are plausible and coherent considering that gender differences reflect ongoing and subtle biochemical processes present throughout childhood, involving and affecting all cells and therefore would be expected to also affect tumor cells [45]. Some miRNAs consistently absent in all Rbs might function as tumor suppressors, since we found

that 144 out of 995 miRNAs that are consistently “off” in all these tumors are indeed located in areas previously reported as recurrently deleted in Rb [35–37].

Analysis of the discretized data of cluster17–92 in Rb indicates that these miRNAs belong to cluster ‘P’ of the landscape miRNome map, which is detected in most but not all the samples under study. Our results shows heterogeneity in the expression of cluster 17–92, with different levels of expression for members encoded by the same primary transcript, consistent with other reports [46, 47]. Even though some paralog members of cluster 17–92 show high expression levels, including miR-19 which is able to recapitulate the oncogenic activity of the full cluster and is located in a locus previously reported as amplified in Rb [29], they are not expressed in all the samples and do not belong to the 30 miRNA core described previously.

Tumor heterogeneity can be considered to reflect biological variability, and heterogeneity can be thought of as a multilayered structure since there are so many aspects of tumor biology that show heterogeneity, including response to treatment or genetic and phenotypic traits. This work reveals consistent similarities and differences or heterogeneity in several Rb tumors at the miRNA level. There is a need to better understand tumor heterogeneity and we propose that tumor heterogeneity be considered as a ‘composite’ of similarities and differences among cancer cells. This definition can be useful and practical for interpretation of high throughput data, and with this work we explain at the miRNome level the nature of these similarities and differences. These results indicate that tumor cells share a fundamental unity despite of intra or inter tumoral heterogeneity and we interpret that the diversity found may account at least in part, to inter tumor heterogeneity in miRNA terms. Tumor cells must share a fundamental unity and such a ‘unity’ has yet to be defined for every type of tumor. This endeavor is particularly challenging for solid tumors which are composed of many cell types in addition to tumor cells. Nonetheless, defining those characteristics that are shared between tumors may help determine those critical to promoting tumor survival.

There are some methodological considerations to note. We cultured the tumors because the amount of tissue we can ethically obtain from an enucleated eye with Rb is limited and necrosis poses an additional constraint. Despite this, there are advantages to culturing cells for a week. First, because we discarded dead cells and propagated viable cells through a standardized procedure that allowed us to generate sufficient material, we were able to assume constant conditions for all tumors. Secondly, because Rb cells grow in suspension we were able to separate out any adherent non-tumor cells. The profiles we obtained, thus originate from a very highly enriched
collection of tumor cells with limited or no miRNAs' signals from non-tumor cells. Our approach assumes that we have detected the average miRNA signals that were contributed during RNA extraction from different clones or tumor cell populations that might have existed within the original tumor sample grown in culture. This is an important limitation since intra-tumoral heterogeneity requires sampling at multiple tumor sites, which we were unable to do.

Two additional and important limitations are inherent to our experimental design which does not include normal tissues. First normal retina from this age group is extremely difficult to obtain even from children that have died from accidental causes. Second, we were not interested in comparing normal versus tumor tissue, the most obvious comparison in the cancer field made since the invention of the two channel microarray platforms. Our goal was to address a knowledge gap in the field, specifically what is shared and what is not shared among tumors with the same diagnosis. The chief advantage of this design is that it allows uncovering what is commonly expressed among the samples and also unveils underlying variability. The disadvantage of our design is that it does not permit identifying what is common or different in comparison with the normal counterpart, whether retina or Rb's cell of origin.

## Conclusions

This miRNA landscape approach reveals the existence of a set of 142 miRNAs shared by all Rb and an additional set of miRNAs which are variably expressed. Our exploration of mRNA targets and pathways affected by a core group of the 30 most consistently highly expressed miRNAs within the 142 shared by all Rb, suggests that this core belongs to and impacts stable genetic networks that have structures shared across samples. Some of these pathways belong to fundamental biologic processes involved in cancer, while others are related to neural functions reflecting the tissue of origin. Still other miRNAs target basic processes such as metabolism, which may merely reflect the shared human eukaryotic origin and more fundamental genetic networks affected by this core.

Our results indicate that human Rbs, share a common and fundamental miRNA expression profile despite their heterogeneity. The therapeutic implication of targeting mir-3613-3p for Rb patients remains to be studied. Studying characteristics shared among tumor cells, given the selection processes occurring in malignant tumors, may allow finding critical targets or the long sought 'Achilles heel' of cancer cells. Integrating global mRNA and proteomic data from these tumors could also improve our understanding of the role of mir-3613-3p and the effect it exerts on its 36 potential targets that

are tumor suppressors. Furthermore data derived from the landscape presented here can also be used for modeling cancer genetic networks using a systems biology approach.

## Additional files

**Additional file 1:** Names of cel/chp files corresponding to each sample. (DOCX 14 kb)

**Additional file 2:** Validation data and list of miRNA clusters from the Rb landscape. (XLSX 726 kb)

**Additional file 3:** Analysis of target genes, pathways, cellular processes and miRNAs located in lost regions. (XLSX 148 kb)

## Abbreviations

HC: Hierarchical clustering analysis; MEV: Multi Experiment Viewer; miRNAs/ miRs: microRNAs; Rb: Retinoblastoma; RMA: Robust Multichip Analysis;

SAM: Significant Analysis for Microarrays

## Acknowledgements

We thank Josefina Romero-Rendón for excellent support for recruiting families for this study and Marco Ramírez, Daphne García and Yolanda Vázquez for access to patients.

## Funding

This work and manuscript has been funded by CONACYT FOSISS-I010/270/13 and NIH grants CA167833, CA192662, CA98180 (MO). AJ-GCH was supported by CONACYT postdoctoral fellowship.

## Availability of data and materials

The datasets generated during the current study are available in NCBI's Gene Expression Omnibus repository and are accessible through GEO Series accession number GSE84747 http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84747. The datasets are also available in the ArrayExpress repository (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-4977 https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-4977/.

## Authors' contributions

MVPC and BECM designed the experiments. BECM performed the majority of the experiments. MVPC, BECM interpreted the results and wrote the manuscript, MO was major contributor in writing the manuscript, JC, MDJO and NDF did critical readings of the manuscript, AJGC gave outstanding bioinformatics support, MDLC, SSP, CL and CHG facilitated tumor tissue, histopathological and clinical data, ACW critically supervised experiments, AHA cultured primary tissues and maintained tissue bank. All authors read and approved the final manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Consent for publication

Not applicable.

## Ethics approval and consent to participate

All samples were collected after informed written consent from parents of the children participating in a larger IRB (Comisión ética y científica IMSS R2012-785-039 and Comisión de investigación, ética y bioseguridad HIM/ 2012/054) approved case-series study [16, 17] and were treated at Hospital de Pediatria, Instituto Mexicano del Seguro Social (IMSS) and Hospital Infantil de México Federico Gómez in Mexico City.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author details

¹Medical Research Unit in Infectious Diseases, Hospital de Pediatría, CMN SXXI, Instituto Mexicano del Seguro Social, Av. Cuauhtémoc 330, 06720 Mexico City, Mexico. ²Epidemiology Department, Columbia University, New York, USA. ³Pharmacology Department, CINVESTAV, Mexico City, Mexico. ⁴Pathology Department, Hospital Infantil de México Federico Gómez, Secretaría de Salud, Mexico City, Mexico. ⁵Unidad Profesional Interdisciplinaria de Biotecnología, Instituto Politécnico Nacional, Mexico City, Mexico. ⁶Medical Research Unit in Human Genetics, Hospital de Pediatría, CMN SXXI, Instituto Mexicano del Seguro Social, Mexico City, Mexico. ⁷Ophthalmology Department, Hospital de Pediatría, CMN SXXI, Instituto Mexicano del Seguro Social, Mexico City, Mexico. ⁸Ophthalmology Department, Hospital Infantil de México Federico Gómez, Mexico City, Mexico.

Received: 10 September 2016 Accepted: 9 June 2017

Published online: 01 July 2017

References

1. Wei JS, Johansson P, Chen QR, et al. MicroRNA profiling identifies cancer-specific and prognostic signatures in pediatric malignancies. Clin Cancer Res. 2009;15(17):5560–8.

2. Volinia S, Calin GA, Liu C-G, et al. A microRNA expression signature of human solid tumors defines cancer gene targets. Proc Natl Acad Sci U S A. 2006;103(7):2257–61.

3. Lü J, Qian J, Chen F, et al. Differential expression of components of the microRNA machinery during mouse organogenesis. Biochem Biophys Res Commun. 2005;334(2):319–23.

4. Jenkinson H. Retinoblastoma: diagnosis and management—the UK perspective. Arch Dis Child. 2015;100(11):1070–5.

5. He L, Thomson JM, Hemann MT, et al. A microRNA polycistron as a potential human oncogene. Nature. 2005;435(7043):828–33.

6. Conkrite K, Sundby M, Mukai S, et al. Mir-17~92 cooperates with RB pathway mutations to promote retinoblastoma. Genes Dev. 2011;25(16):1734–45.

7. Olive V, Jiang I, He L. Mir-17-92, a cluster of miRNAs in the midst of the cancer network. Int J Biochem Cell Biol. 2010;42(8):1348–54.

8. Xu X, Jia R, Zhou Y, et al. Microarray-based analysis: Identification of hypoxia-regulated microRNAs in retinoblastoma cells. Int J Oncol. 2011;38(5):1385–93.

9. Beta M, Venkatesan N, Vasudevan M, et al. Identification and insilico analysis of retinoblastoma serum microRNA profile and gene targets towards prediction of novel serum biomarkers. Bioinform Biol Insights. 2013;7:21–34.

10. Martin J, Bryar P, Mets M, et al. Differentially expressed miRNAs in retinoblastoma. Gene. 2013;512(2):294–9.

11. Jo DH, Kim JH, Cho CS, et al. STAT3 inhibition suppresses proliferation of retinoblastoma through down-regulation of positive feedback loop of STAT3/miR-17-92 clusters. Oncotarget. 2014;5(22):11513–25.

12. Zhao JJ, Yang J, Lin J, et al. Identification of miRNAs associated with tumorigenesis of retinoblastoma by miRNA microarray analysis. Childs Nerv Syst. 2009;25(1):13–20.

13. Thériault BL, Dimaras H, Gallie BL, Corson TW. The genomic landscape of retinoblastoma: A review. Clin Exp Ophthalmol. 2014;42(1):33–52.

14. Burrell RA, McGranahan N, Bartek J, Swanton C. The causes and consequences of genetic heterogeneity in cancer evolution. Nature. 2013;501(7467):338–45.

15. Jamal-Hanjani M, Quezada SA, Larkin J, Swanton C. Translational Implications of Tumor Heterogeneity. Clin Cancer Res. 2015;21(6):1258–66.

16. Orjuela MA, Cabrera-Muñoz L, Paul L, et al. Risk of retinoblastoma is associated with a maternal polymorphism in dihydrofolatereductase (DHFR) and prenatal folic acid intake. Cancer. 2012;118(23):5912–9.

17. Ramírez-Ortiz MA, Ponce-Castañeda MV, Cabrera-Muñoz ML, et al. Diagnostic delay and sociodemographic predictors of stage at diagnosis and mortality in unilateral and bilateral retinoblastoma. Cancer Epidemiol Biomark Prev. 2014;23(5):784–92.

18. Statistical Algorithms Reference Guide. http://media.affymetrix.com/support/technical/technotes/statistical_reference_guide.pdf. Accessed 23 June 2017.

19. Gene Chip Expression Analysis Data Analysis Fundamentals. http://media.affymetrix.com/support/downloads/manuals/data_analysis_fundamentals_manual.pdf. Accessed 23 June 2017.

20. Statistical Algorithms Description Document. http://www.affymetrix.com/support/technical/whitepapers/sadd_whitepaper.pdf. Accessed 23 June 2017.

21. Liu W, Mei R, Di X, et al. Analysis of high density expression microarrays with signed-rank call algorithms. Bioinformatics. 2002;18(12):1593–9.

22. Saeed AI, Sharov V, White J, et al. TM4: a free, open-source system for microarray data management and analysis. BioTechniques. 2003;34(2):374–8.

23. Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci U S A. 2001;98(9):5116–21.

24. Vlachos IS, Zagganas K, Paraskevopoulou MD, et al. DIANA-miRPath v3.0: deciphering microRNA function with experimental support. Nucleic Acids Res. 2015;43(W1):W460–6.

25. Agarwal V, Bell GW, Nam J-W, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. elife. 2015; doi:10.7554/eLife.05005.

26. Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36(Database issue):D154–8.

27. Zhao M, Sun J, Zhao Z. TSGene: A web resource for tumor suppressor genes. Nucleic Acids Res. 2013; doi:10.1093/nar/gks937.

28. Zhao M, Kim P, Mitra R, Zhao J, Zhao Z. TSGene 2.0: an updated literature-based knowledgebase for tumor suppressor genes. Nucleic Acids Res. 2016; 44(D1):D1023-31. doi: 10.1093/nar/gkv1268. Epub 2015 Nov 20. PubMed PMID: 26590405; PubMed Central PMCID: PMC4702895.

29. García-Chequer AJ, Méndez-Tenorio A, Olguín-Ruiz G, et al. Overview of recurrent chromosomal losses in retinoblastoma detected by low coverage next generation sequencing. Cancer Genet. 2016;209(3):57–69.

30. Hanahan D. The Hallmarks of Cancer. Cell. 2000;100(1):57–70.

31. Chatterjee A, Leichter AL, Fan V, et al. A cross comparison of technologies for the detection of microRNAs in clinical FFPE samples of hepatoblastoma patients. Sci Rep. 2015;5:10438.

32. Edward DP, Alkatan R, Hafiq Q, et al. MicroRNA Profiling in Intraocular Medulloepitheliomas. PLoS One. 2015;10(3):e0121706.

33. Tanaka K, Kawano M, Itonaga I, et al. Tumor suppressive microRNA-138 inhibits metastatic potential via the targeting of focal adhesion kinase in Ewing's sarcoma cells. Int J Oncol. 2016;48(3):1135–44.

34. Fricke A, Ullrich PV, Heinz J, et al. Identification of a blood-borne miRNA signature of synovial sarcoma. Mol Cancer. 2015;14(1):151.

35. Lillington DM, Kingston JE, Coen PG, et al. Comparative genomic hybridization of 49 primary retinoblastoma tumors identifies chromosomal regions associated with histopathology, progression, and patient outcome. Genes Chromosomes Cancer. 2003;36(2):121–8.

36. Ganguly A, Nichols KE, Grant G, et al. Molecular karyotype of sporadic unilateral retinoblastoma tumors. Retina. 2009;29(7):1002–12.

37. Mol BM, Massink MPG, van der Hout AH, et al. High resolution SNP array profiling identifies variability in retinoblastoma genome stability. Genes Chromosomes Cancer. 2014;53(1):1–14.

38. McEvoy J, Nagahawatte P, Finkelstein D, et al. RB1 gene inactivation by chromothripsis in human retinoblastoma. Oncotarget. 2014;5(2):438–50.

39. Mogilyansky E, Rigoutsos I. The miR-17/92 cluster: a comprehensive update on its genomics, genetics, functions and increasingly important and numerous roles in health and disease. Cell Death Differ. 2013;20(12):1603–14.

40. Witten D, Tibshirani R, Gu SG, et al. Ultra-high throughput sequencing-based small RNA discovery and discrete statistical biomarker analysis in a collection of cervical tumours and matched controls. BMC Biol. 2010;8:58.

41. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646-74. doi:10.1016/j.cell.2011.02.013. Review. PubMed PMID: 21376230.

42. Sato F, Tsuchiya S, Terasawa K, Tsujimoto G. Intra-Platform Repeatability and Inter-Platform Comparability of MicroRNA Microarray Technology. PLoS One. 2009;4(5):e5540.

43. Rushlow DE, Mol BM, Kennett JY, Yee S, Pajovic S, Thériault BL, et al. Characterisation of retinoblastomas without RB1 mutations: genomic, gene expression, and clinical studies. Lancet Oncol. 2013;14(4):327–34.

44. Weinstein IB, Joe A. Oncogene addiction. Cancer Res. 2008;68(9):3077–80. discussion 3080

45. Sun T, Plutynski A, Ward S, Rubin JB. An integrative view on sex differences in brain tumors. Cell Mol Life Sci. 2015; doi:10.1007/s00018-015-1930-2.

46. Concepcion CP, Bonetti C, Ventura A. The microRNA-17-92 family of microRNA clusters in development and disease. Cancer J. 2012;18(3):262–7.

47. Obernosterer G, Leuschner PJ, Alenius M, Martinez J. Post-transcriptional regulation of microRNA expression. RNA. 2006;12(7):1161–7.