Integrative multi-omics analysis identifies a prognostic miRNA signature and a targetable miR-21-3p/TSC2/mTOR axis in metastatic pheochromocytoma/paraganglioma

Rationale: Pheochromocytomas and paragangliomas (PPGLs) are rare neuroendocrine tumors that present variable outcomes. To date, no effective therapies or reliable prognostic markers are available for patients who develop metastatic PPGL (mPPGL). Our aim was to discover robust prognostic markers validated through in vitro models, and define specific therapeutic options according to tumor genomic features. Methods: We analyzed three PPGL miRNome datasets (n=443), validated candidate markers and assessed them in serum samples (n=36) to find a metastatic miRNA signature. An integrative study of miRNome, transcriptome and proteome was performed to find miRNA targets, which were further characterized in vitro. Results: A signature of six miRNAs (miR-21-3p, miR-183-5p, miR-182-5p, miR-96-5p, miR-551b-3p, and miR-202-5p) was associated with metastatic risk and time to progression. A higher expression of five of these miRNAs was also detected in PPGL patients' liquid biopsies compared with controls. The combined expression of miR-21-3p/miR-183-5p showed the best power to predict metastasis (AUC=0.804, P=4.67·10-18), and was found associated in vitro with pro-metastatic features, such as neuroendocrine-mesenchymal transition phenotype, and increased cell migration rate. A pan-cancer multi-omic integrative study correlated miR-21-3p levels with TSC2 expression, mTOR pathway activation, and a predictive signature for mTOR inhibitor-sensitivity in PPGLs and other cancers. Likewise, we demonstrated in vitro a TSC2 repression and an enhanced rapamycin sensitivity upon miR-21-3p expression. Conclusions: Our findings support the assessment of miR-21-3p/miR-183-5p, in tumors and liquid biopsies, as biomarkers for risk stratification to improve the PPGL patients' management. We propose miR-21-3p to select mPPGL patients who may benefit from mTOR inhibitors.

. The x-axis shows the log 2 fold change of miRNA expression, whereas the y-axis shows the -log 10 (FDR) or -log 10 (P) for each miRNA. MiRNAs with P<0.05 are colored to indicate whether they are up-regulated (red) or down-regulated (green) in tumors of patients with mPPGL. miRNAs selected for validation are indicated. (B) Expression of miRNAs selected for validation in the different series. Log 2 -normalized expression from the different series is displayed as a transformed z-score (centered at the mean of the non-metastatic group for each series). miRNA expression in primary tumors of patients with non-mPPGL is shown in green, in primary tumors of patients with mPPGL in orange, and in metastatic tissues in red (metastatic tissues were only available from 8 patients, whose paired primary tumors are included in the orange group and the expression was only evaluated for the validated miRNAs). Data are shown as mean ± SD. * FDR<0.05 for the discovery series obtained during the differential expression analysis, and P<0.05 from a nonparametric Mann-Whitney test in the validation series.      Plots show levels of the indicated miRNAs determined by ddPCR in conditioned medium preparations (cells incubated for 24h in serum-free media) equivalent to 60,000 cells. Levels of miR-202-5p and miR-551b-3p were undetectable. Error bars indicate standard error of the mean from three independent conditioned medium preparations. Unpaired t-test was applied to test for differences (*: P<0.05).

Table S1. Clinical data of patients included in the validation series and patients with primary-metastatic paired samples available
For the validation series, twenty-four cases were classified as non-metastatic, as patients were disease free at the time of the last clinical follow-up (≥ 4 years follow-up, median=8.5 years). The remaining twenty-five samples were classified as metastatic since the patients presented tumor cells at non-chromaffin sites. The proportion of SDHB cases per group (28% in the metastatic and 8% in the non-metastatic group) was similar to that observed in the discovery series (average in the three sub-series of 24.3% and 4.4% in metastatic and non-metastatic groups, respectively). WT: wild type; PGL: Paraganglioma; PCC: Pheochromocytoma.
Up Yes

miR-96-5p
Already reported to be upregulated in PPGLs [64][65][66][67] Oncogene Up Yes miR-183-5p miR-182-5p Information obtained by quick reviewing abstracts from references obtained after performing the following advanced search of the PubMed database: "(miRNA name[Title/Abstract] AND cancer[Title/Abstract])". Last literature review: 28/August/2018. Only miRNAs whose "potential role according to literature" was not in disagreement with "status in PPGL Discovery series" were selected for validation.

CCND2
-0.09 4.07E-01 A_24_P270235 -0. 45   RAI2 function has not been characterized, but its expression is induced by retinoic acid that plays a critical role in development, cell growth, and differentiation 3 Tumor suppressor reported in breast and colorectal cancer [71][72][73] Tumor suppressor YES

SMAD7, SMAD family member 7
Involved in regulation of the TGF-β pathway 310 Dual role in the regulation of cancer progression [74,75], mostly reported to be a negative regulator of the pathway involved in epithelial-to-mesenchymal transition and metastasis [76] Mostly tumor suppressor YES

PLCB4, phospholipase C beta 4
Catalyzes the formation of inositol trisphosphate and diacylglycerol from PIP2: involved in the transduction of signals in the retina 6 Gain-of function alterations reported in cancer [77][78][79][80][81] Oncogene NO

STK35, serine/threonine kinase 35
Regulator of actin stress fibers in non-muscle cells 3 Might have an oncogenic role in osteosarcoma [82] ? NO

CETN3, centrin 3
Calcium binding-protein; it plays a fundamental role in centrosome duplication and separation 5 Loss of CETN3 connected to elevated EGFR [83] ? NO

KLHDC3, kelch domain containing 3
Involved in V(D)J and meiotic recombination, specifically expressed in testis 0 -? NO

PEX12, peroxisomal biogenesis factor 12
Essential for the assembly of functional peroxisomes 1 -? NO

Sphingosine-1-Phosphate Lyase 1 -13
Sphingosine 1-phosphate (S1P) is a lipid with important roles in growth, survival and migration. It has been largely studied in cancer and SGPL1 is the phosphatase that irreversibly cleaves S1P and is the only exit point from sphingolipid pathways [84,85] Tumor suppressor YES

CCND2, cyclin D2
Forms a complex with CDK4/CDK6 and functions as a regulatory subunit of the complex, whose activity is required for cell cycle G1/S transition . By doing this search we aimed to identify those genes that have been largely reported to be involved in cancer. ¤ Information obtained by quick reviewing abstracts from references obtained in Ꙙ , which might be useful to establish a potential role of the gene in cancer.
Only those miRs whose "potential role according to literature" was in agreement with a potential involvement in the disease were selected for further study.

Discovery of differentially expressed miRNAs
Merging the three datasets was not possible since different platforms had been used, and therefore they were treated as three sub-series of a larger discovery series as specified bellow: Sub-series 1: Microarray image acquisition and analysis was done using a G2505C microarray scanner (Agilent). The text files with the data of the processed images were analyzed with limma R package (version 3.26.9). Background correction was done using the "normexp" method. Normalization between arrays was done applying the "quantile" method. Values for within-array replicate probes were replaced with their average by using the avereps function of limma.
After estimating the fold changes and standard errors by fitting a linear model for each probe and applying the empirical Bayes to smooth the standard errors, a list of probes was obtained that were differentially expressed between metastatic and non-metastatic samples. The miRNA names corresponding to the probes were converted to those of release 20 of miRBase [6] by using an in-house Perl script that processed the conversion files downloaded from miRBase Tracker [7].
Sub-series 2: The miRNA read files of the samples in FastQ format were aligned to release 20 of miRBase [6] using bowtie [8] (version 0.12.7) with seed length equal to 17; no mismatches were allowed in the seed and at most one alignment per read was allowed. Htseq-count [9] (version 0.5.3p9) was used to generate the read count data for each miRNA.
Only miRNAs having a minimum of 15 counts in at least 7.3 % of the samples in both the subset of metastatic samples and the subset of non-metastatic samples were selected for further the analysis.
The edgeR package (version 3.12.1) [10] was then used to normalize the matrix using the trimmed mean of M-values (TMM) [11] and extract the list of differentially expressed genes by using the exactTest function.
Only miRNAs with an average count >100 in any of the groups were considered further.
Sub-series 3: An in-house Perl script was used to a) sum all the read counts corresponding to all isoforms of each mature miRNA accession number; b) convert the mature miRNA sequence names into the names in release 20 of miRBase; and c) generate a matrix containing the counts of all miRNAs in all input TCGA samples.
Only miRNAs having a minimum of 15 counts in at least 5.4 % of the samples in both subsets of metastatic and non-metastatic samples were selected for further analysis.
The R package edgeR was then used to normalize the matrix using the TMM and obtain the genes differentially expressed between non-metastatic and metastatic samples by using the exactTest function.
Only miRNAs with an average count >100 in any of the groups were considered further, as they were deemed easily detectable markers.  [12].

Integration of miRNA and mRNA expression profiles of the discovery series
Step 2. Generation of normalized mRNA expression matrices including only step 1 genes (Targetome) for each sub-series was performed as specified below.
Sub-series 1: Microarray image acquisition and analysis was done with GenePix TM Pro (Axon Instruments Inc).
The produced .gpr files were read and analyzed with limma R package [13] (version 3.26.9). Background correction was done using the "normexp" method. Normalization within arrays was performed using the "loess" method. Normalization between arrays was done applying the "quantile" method.

Sub-series 2:
We downloaded the Affymetrix (GeneChip Human Genome U133 Plus 2.0) intensity (.CEL) files of 177 samples from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) experiment ID E-MTAB-733. We read and processed these .CEL files using the affy and limma packages of R (versions 1.48.0 and 3.26.9, respectively). Background correction was done using the "rma" method, and the normalization method used was "quantiles".
Sub-series 3: An in-house Perl script collapsed these files into a unique matrix containing the raw counts for all samples. R package edgeR (version 3.12.1) was then used to normalize the matrix using the TMM.
Generation of targetome matrices: For each series, miRNA and mRNA expression submatrices were obtained from the miRNA and mRNA expression matrices previously generated. These submatrices contained the expression values for the sample codes shared by the miRNA expression matrix and the mRNA expression matrix. The number of sample codes shared by miRNA expression and mRNA expression matrices in each series was 87 for sub-series 1, 168 for sub-series 2, and 179 for sub-series 3.
Next, a third submatrix was generated from these submatrices for each series, containing only those genes that were included in the lists generated in Step 1 (potential gene targets for either miR-21-3p or miR-183-5p), which we designated targetome of miR-21-3p and targetome of miR-183-5p.
qPCR Regarding the cell model, 100,000 cells of each cell line were seeded in duplicate in 6 well-plates. One of the duplicates was seeded using complete medium, and the other using complete medium plus 1µg/ml doxycycline (Sigma#D9891). RNA was extracted using TRI Reagent® (MRC#TR 118) following the vendor's instructions at different time points (0, 24, 48, 84, 144 and 168 h after plating). Medium was changed after 84h. In the validation series we used the RNA that had been extracted for miRNA quantification and validation in this series.
In all cases, cDNAs were prepared from 500ng of RNA using the qScript TM cDNA Synthesis Kit (#95047-100, Quanta Biosciences, Gaithersburg, MD) and mRNA levels were quantified by real-time PCR using the