HFD modulate the expression of multiple miRNAs in murine prostate tumors
Previously, we reported that HFD significantly increased TRAMP-C1 tumor growth, oncogenes and oncomiRs expression in C57BL/6J male mice [4,7]. In this work, we investigated the impact of HFD on miRNA expression profile from these tumors using a high-throughput platform. GeneChiP miRNA 4.0 Affymetrix was hybridized with total RNA isolated from CD or HFD TRAMP-C1 tumors. After data normalization, we selected differentially expressed miRNAs with a FDR<0.05 and identified 18 down- and 6 up-regulated miRNAs in HFD tumors compared to CD group (Figure 1, Table 1). Particularly, mmu-miR-133a-3p, 133b and 1a-3p showed the most robust expression changes between the two groups. As miR-133a-3p and 1a-3p are expressed from the same miR-1/133a genomic clusters (located on chromosomes 2 and 18 in mouse and 18 and 20 in human and mmu-miR-133a-3p and 133b are part of the same miRNA gene family they were selected for further analyses.
Down-regulated miRNAs by HFD modulate cancer related pathways
To explore the functional role of the differentially expressed miRNAs and pathways, we used DIANA-miRPath v3 tool (http://snf-515788.vm.okeanos.grnet.gr/) and STRING database (https://string-db.org/). Since each miRNA can control several genes, we selected only the deregulated miRNAs that had orthologous in human to obtain a more robust result. We employed DIANA-TarBase v7 of the miRPath tool, which is a reference database of experimentally supported miRNA targets, and followed the Figure 2A workflow. Eight of the down-modulated mouse miRNAs (8 out of 18) and 10 of its orthologous human miRNAs (10 out of 18) showed validated target genes (Figure 2AI). Additionally, two (2 out of 6) of the up-regulated mouse miRNAs and 2 of its human orthologous presented validated target genes. Next, we excluded common target genes among down- and up-regulated miRNAs from the same species and selected common target genes between species (human and mice). This analysis revealed that down-regulated miRNAs showed 1,278 target genes common to human and mice, while up-regulated miRNAs did not present common target genes (Table S1). We continue working with down-modulated miRNAs and their target genes. For this group of genes, KEGG signaling pathways were identified with a FDR< 0.05 using STRING database (Figure 2B, Table S2). Target genes of the down-modulated miRNAs showed enriched processes mainly associated with cancer-related and signal transduction pathways (Figure 2B). To determine the relevance of down-regulated miRNAs-target genes in human samples, we performed a PCA of the 1,278 target genes using the normalized expression data from normal prostate (GTEX), NAT and prostate tumors samples (TCGA-PRAD). The 2-dimensional scatterplot of the first two principal components revealed marked differences in overall gene expression between normal prostate and PCA samples (Figure 2C).
In order to found relevant target genes for the down-modulated miRNAs, within the 1,278 genes, we performed a correlation matrix, using the expression data of the down-modulated miRNAs and their target genes of prostate tumors from patients available in TCGA-PRAD (UCSC Xena). From this analysis, we selected the target genes that showed a negative Spearman correlation coefficient rho <-0.2 and a p-Value <0.05 (Figure 2AII) with its regulatory miRNA. Next, to shorten the number of relevant miRNAs and target genes, we choose common target genes that negatively correlate with the expression of the previously selected miRNAs hsa-miR-133a-3p, 133b and 1-3p. We found nine common target genes (BCL7A, CENPF, CPSF6, ELAVL1, NRP1, PAICS, RALA, RBM15B and WHSC1) with a significant negative correlation with the mentioned miRNAs (Figure 3A, Table S3). Also, hsa-miR-133a-3p, 133b and 1-3p showed a strong positive correlation between them with a significant p-value (Figure 3A). Additionally, we performed a ssGSEA to explore the degree to which the 9 common target genes were coordinately up- or down-regulated within the prostate tumor samples in two different PCa datasets: TCGA-PRAD and Taylor (GSE21032). As shown in Figure 3B-C we found that the 9 target genes showed a positive ssGSEA-enrichment in most of PCa samples in both datasets. Moreover, a strong negative correlation between the ssGSEA-enrichment and the expression of miRs 133a-3p, miR-133b and miR-1-3p were found in the TCGA-PRAD dataset. A similar result was found for 133a-3p, miR-133b in Taylor dataset. Finally, we analyzed the contribution of the 9 target genes to specific gene sets annotated in MSigDB collections. As shown in Figure 3D, three of the nine target genes (WHSC1, CENPF and PAICS) were found to be significantly enriched in cancer-related gene sets.
Hsa-miR-133a-3p, miR-133b and miR-1-3p expression are decreased in PCa
In order to further understand the pattern expression of the three selected miRNAs in prostate tumors from patients, we performed a bioinformatic analysis using the TCGA-PRAD, and Taylor (GSE21032) datasets. First, we evaluated the expression of the selected miRNAs in prostate primary tumors in comparison to NAT and normal prostates. As shown in Figure 4A, we found hsa-miR-133a-3p, 133b and 1-3p significantly downregulated in PCa compared to NAT (TCGA-PRAD cohort). Additionally, hsa-miR-133a-3p, 133b and 1-3p expression were significantly decreased in PCa compared to normal prostate (Taylor cohort) (Figure 4B). These results proposed that hsa-miR-133a-3p, miR-133b and miR-1-3p might act as tumor suppressors in PCa.
Hsa-miR-133a-3p, miR-133b and miR-1-3p expression negatively correlated with Gleason grade
To correlate hsa-miR-133a-3p, 133b and 1-3p with PCa stage, we used TCGA-PRAD and Taylor (GSE21032) patient cohorts, which were divided into groups according to Gleason’s score on histopathology (6, 7 and 8-9). As shown in Figure 5, hsa-miR-133a-3p, 133b and 1-3p were found to be significantly down-regulated according to Gleason score increased in both datasets. Although the expression of the selected miRNAs did not affect overall survival in PCa patients (data not shown), a low expression of hsa-miR-133a-3p, 133b and 1-3p correlated with worst progression-free interval (Figure 5 A-C) in the TCGA-PRAD cohort. Thus, hsa-miR-133a-3p, 133b and 1-3p low expression correlated with high gleason score and worst progression free interval.
Hsa-miR-1-2/miR-133a-1 promoter hypermethylation is increased in PCa and negatively correlated with miRNA expression
The hsa-miR-1/133 family is located at three different loci (as clustered miRNAs) at chromosomes 18q11.2 (miR-1-2/miR-133a-1), 20q13.33 (miR-1-1/miR-133a-2), and 6p12.2 (miR-206/miR-133b) [17]. Here, we analyzed the methylation status at miR-1-2/miR-133a-1, miR-1-1/miR-133a-2 and miR-206/miR-133b promoters. Based upon the Ensembl and UCSC Xena resources, two, six and thirteen methylation probes were identified targeting the hsa-miR-1-2/miR-133a-1, miR-133b and miR-1-1/miR-133a-2 promoters, respectively (Figure 6, Table 2). We compared the methylation status at the three cluster promoters in 52 paired PCa and NAT samples. As shown in Figure 6, one probe (cg17106157) targeting miR-1-2/133a-1 (Figure 6A) and two probes (cg14487577 and cg22720139) targeting miR-133b promoters (Figure 6B) showed significantly higher beta values in PCa compared to NAT. Also, ten and one probes targeting miR-1-1/133a-2 promoter showed lower and higher beta values in PCa vs NAT respectively (Figure 6C). Then, a correlation analysis between promoter methylation and mature miRNA expression levels was performed in 497 PCa samples. As presented in Table 2, Spearman correlation analysis showed that hsa-miR-133a-3p and hsa-miR-1-3p expression negatively correlated with cg17106157 (chromosome 18) probe and hsa-miR-133a-3p and hsa-miR-1-3p expression negatively correlated with cg05898333 (chromosome 20) probe. Moreover, hsa-miR-133a-3p and hsa-miR-1-3p expression positively correlated with cg15580304, cg14523475, cg08148458 and cg22617703 probes targeting the miR-1-1/miR-133a-2 promoter. Therefore, the decreased expression observed in hsa-miR-1-3p/miR-133a-3p in prostate tumors might be due to a hsa-miR-1-2/miR-133a-1 cluster promoter hypermethylation.
CENPF and WHSC1 are increased in primary prostate tumors from patients and positively correlated with Gleason grade
We found that hsa-miR-133a-3p, 133b and 1-3p showed nine common target genes with negative correlation. Then, we analyzed the expression pattern of the nine target genes in prostate tumors and normal samples from TCGA-PRAD and Taylor (GSE21032) datasets. As shown in Figure 7A, we found that CENPF, WHSC1, CPSF6, RBM15B, PAICS and ELAVL1 expression was significantly increased in primary prostate tumors compared to NAT (TCGA-PRAD dataset). Moreover, CENPF, WHSC1 and PAICS expression were significantly increased in prostate tumors compared to normal prostates (Taylor dataset, Figure 7B). With these results we analyzed the correlation of CENPF, WHSC1 and PAICS with PCa stage. CENPF and WHSC1 were found to be significantly up-regulated according as Gleason score increases in TCGA-PRAD and Taylor datasets (Figure 8A-B). Moreover, a high expression of these genes correlates with a worst progression free interval (Figure 8A-B). These results collectively suggest that common target genes of hsa-miR-133a-3p, 133b and 1-3p act as oncogenes favoring PCa aggressiveness.