Abstract

   (1) Background: Glioblastoma (GBM) is one of the most aggressive brain
   tumors with a poor prognosis. Therefore, new insights into GBM
   diagnosis and treatment are required. In addition to differentially
   expressed mRNAs, miRNAs may have the potential to be applied as
   diagnostic biomarkers. (2) Methods: In this study, profiling of human
   miRNAs in combination with mRNAs was performed on total RNA isolated
   from tissue samples of five control and five GBM patients, using a
   high-throughput RNA sequencing (RNA-Seq) approach. (3) Results: A total
   of 35 miRNAs and 365 mRNAs were upregulated, while 82 miRNAs and 1225
   mRNAs showed significant downregulation between tissue samples of GBM
   patients compared to the control samples using the iDEP tool to analyze
   RNA-Seq data. To validate our results, the expression of five miRNAs
   (hsa-miR-196a-5p, hsa-miR-21-3p, hsa-miR-10b-3p, hsa-miR-383-5p, and
   hsa-miR-490-3p) and fourteen mRNAs (E2F2, HOXD13, VEGFA, CDC45, AURKB,
   HOXC10, MYBL2, FABP6, PRLHR, NEUROD6, CBLN1, HRH3, HCN1, and RELN) was
   determined by RT-qPCR assay. The miRNet tool was used to build
   miRNA–target interaction. Furthermore, a protein–protein interaction
   (PPI) network was created from the miRNA targets by applying the
   NetworkAnalyst 3.0 tool. Based on the PPI network, a functional
   enrichment analysis of the target proteins was also carried out. (4)
   Conclusions: We identified an miRNA panel and several deregulated mRNAs
   that could play an important role in tumor development and distinguish
   GBM patients from healthy controls with high sensitivity and
   specificity using total RNA isolated from tissue samples.

   Keywords: glioblastoma, miRNAs, brain tissue, next-generation
   sequencing, biomarker

1. Introduction

   Glioblastoma (GBM; World Health Organization grade 4) is the most
   common and aggressive primary brain cancer in adults and is considered
   to be the most prevalent form of brain tumor leading to death [[42]1].
   This incurable malignant tumor has a median survival time of about 15
   months from diagnosis; the 5-year survival rate is only 10%. Typically,
   GBM appears in the sixth decade of life, and it is slightly predominant
   in males [[43]2]. The current standard of treatment (called “Stupp’s
   regimen”) includes surgical resection of the tumor followed by
   radiotherapy combined with adjuvant temozolomide (TMZ) chemotherapy
   [[44]3]. Additionally, personalized therapeutic agents against specific
   deregulated targets that could be responsible for the induction of
   tumor growth have already been tested in several clinical trials.
   However, almost all GBM patients undergo unavoidable tumor recurrence
   [[45]4]. The molecular mechanisms behind the development of GBM are
   still not completely understood despite the recent achievements in GBM
   research. Therefore, it is crucially important to identify other
   factors that could contribute to the onset and progression of GBM.
   According to that goal, epigenetics and epigenetic modulators have come
   into focus that are involved in the development of cancer [[46]5].

   The latest members of the epigenetic machinery are the noncoding RNAs
   (ncRNAs), which either have a limited ability to encode a protein or
   lack it. MicroRNAs (miRNAs) are small 21–25 nucleotide-long ncRNAs that
   function as major players in the post-transcriptional regulation of
   protein-coding genes via their sequence-specific binding to the 3′
   untranslated regions (UTR) of target mRNAs [[47]6]. The results of
   different research groups show that miRNAs as post-transcriptional
   regulators are involved in the regulation of several important
   processes, such as cell differentiation, cell division, apoptosis, cell
   metabolism, and patterning of the nervous system [[48]7]. Different
   mechanisms like defects in miRNA biogenesis, abnormalities in miRNA
   processing, and epigenetic alterations and mutations in the miRNA
   recognition sites of target genes can lead to changes in miRNA
   expression and, as a consequence, can initiate tumorigenesis
   [[49]8,[50]9]. In the case of various cancers, many miRNAs show
   tumor-specific expression patterns and significant deregulation. MiRNAs
   can function as either oncogenes (oncomiRs) or tumor suppressors;
   however, because of the large number of target genes, the same miRNA
   may play opposing roles in different tumor types [[51]10].
   Characterizing miRNA expression in GBM could be applied as a potential
   diagnostic or prognostic tool; furthermore, miRNAs and their targets
   may be helpful in selecting appropriate therapy [[52]11]. Additionally,
   the simultaneous analysis of mRNA expression can add more details to
   the analysis of regulatory networks of differently expressed (DE)
   miRNAs in GBM.

   The objective of this study was to profile miRNA expression in
   conjunction with mRNA expression in identical tissue samples of GBM
   patients (GPs), who were diagnosed and treated at the Department of
   Neurosurgery, Faculty of Medicine, University of Debrecen, Hungary, to
   detect those miRNAs and mRNAs whose significant deregulation correlates
   with tumorigenesis. Previously, it was shown that ethnicity seemed to
   be one of the potential sources of heterogeneity between studies
   [[53]12,[54]13]; therefore, it is inevitable to screen for
   significantly deregulated miRNAs that are specific for a patient cohort
   at a certain geographic location to create a miRNA panel that could be
   applied for diagnosis of GBM in that specific region. The limitation of
   the use of a single-miRNA biomarker is that its significant
   deregulation could be associated with various types of cancers, so it
   is sensible to apply a combination panel of miRNAs in cancer diagnosis
   [[55]14]. In our work, we applied a high-throughput next-generation RNA
   sequencing method and bioinformatics analysis to determine the
   significantly deregulated miRNAs, together with mRNAs in the same
   tissue samples of GPs. The identified miRNAs and their targets were
   used to build interaction networks, which were subjected to Gene
   Ontology (GO) and pathway enrichment analyses in order to identify the
   significant molecular pathways that could be affected by the
   deregulated miRNA pattern. The identification of miRNAs whose
   expression was significantly different and their respective targets in
   the tissue samples of GPs could provide further insights and facilitate
   our understanding of the pathogenesis of GBM.

2. Results

2.1. Identification of Differently Expressed (DE) miRNAs in Tissue Samples of
GBM Patients and Control Group

2.1.1. Next-Generation Sequencing (NGS)

   To analyze the alterations in miRNA expression patterns in tissue
   samples of GPs, next-generation sequencing (NGS) technology was
   employed using an Illumina NextSeq 500 instrument (Illumina, San Diego,
   CA, USA). Small RNA-Seq sequencing libraries were generated from total
   RNA prepared from the surgically removed tumor tissue of five GPs and
   the peripheral tumor region of five lower-grade (grade 1–2) glioma
   patients, serving as a control group.

2.1.2. Hierarchical Clustering with Heatmap and Principal Component (PCA)
Analysis

   To visualize the differences in the expression patterns of miRNAs
   between GPs and controls, a hierarchical cluster analysis was performed
   on the NGS dataset of five replicates for each of the two groups to
   identify miRNAs with similar expression patterns. MiRNAs were ranked
   based on their standard deviation across all samples, and the top 200
   genes were used in hierarchical clustering. Expression profiles of the
   top 50 most variable miRNAs in the samples of GPs and controls are
   depicted as a heatmap in [56]Figure 1.

Figure 1.

   [57]Figure 1
   [58]Open in a new tab

   Heatmap with hierarchical clustering dendrogram of miRNA expression in
   GBM and control samples based on the variation in miRNA expression. The
   expression profiles of the top 50 miRNAs are shown. Columns represent
   patient and control individuals (G—glioblastoma sample and C—control
   sample), and each row represents a single miRNA. The down- and
   upregulated miRNAs are labeled green and red, respectively.

   Principal component analysis (PCA) was also performed to increase
   interpretability and visualize the distribution of miRNA expression
   values. A PCA plot using the first and second principal components is
   presented in [59]Figure 2. GBM samples form a single cluster based on
   their expression pattern and are clearly separated from the control
   samples by the first principal component, which represents 42% of the
   variance. This data distribution suggests that GBM biogenesis induces a
   drastic change in the expression of several miRNAs.

Figure 2.

   [60]Figure 2
   [61]Open in a new tab

   PCA of miRNA expression based on their expression profile. A clear
   separation is visible between the GBM samples and the control samples
   along the first principal component.

2.1.3. Differentially Expressed Genes (DEGs)

   Using the DESeq2 algorithm of the iDEP.96 web tool and applying
   thresholds of false discovery rate (FDR) < 0.05 and fold-change > 2, we
   identified 117 miRNAs whose expression proved to be significantly
   different in GPs compared to the control samples. Among the 117
   differently expressed miRNAs, 35 showed upregulation (log[2]FC ≥ 1),
   and 82 showed downregulation (log[2]FC ≤ −1). The whole list of the 117
   deregulated miRNAs is presented in [62]Table S1. The four most strongly
   upregulated miRNAs were hsa-miR-10b-5p (log[2]FC = 6.9),
   hsa-miR-196a-5p (log[2]FC = 5.62), hsa-miR-10a-5p (log[2]FC = 5.48),
   and hsa-miR-21-3p (log[2]FC = 4.39). Conversely, the four most
   downregulated miRNAs were hsa-miR-383-5p (log[2]FC = −6.3),
   hsa-miR-129-5p (log[2]FC = −5.96), hsa-miR-129-2-3p (log[2]FC = −5.89),
   and hsa-miR-219a-2-3p (log[2]FC = −5.8).

   To assess the precise effect of GBM development on miRNA expression
   patterns compared to control samples, k-means clustering was used to
   cluster miRNAs into groups by their expression values ([63]Figure 3a).
   The heatmap and volcano plot of the expression values ([64]Figure 3b)
   show that GBM development results in a massive change in the miRNA
   transcriptome.

Figure 3.

   [65]Figure 3
   [66]Open in a new tab

   (a) K-means clustering of the 117 differentially expressed miRNA genes
   detected with the DESeq2 algorithm. The down- and upregulated miRNAs
   are labeled green and red, respectively. (b) The volcano plot shows
   that GBM development leads to a massive change in the miRNA
   transcriptome. The labeled miRNAs (hsa-miR-196a-5p, hsa-miR-21-3p,
   hsa-miR-10b-3p, hsa-miR-383-5p, and hsa-miR-490-3p) were included in
   the validation group by RT-qPCR.

2.1.4. miRNA Ranking by Network-Based Analysis

   miRNAs can function either as oncogenes or tumor suppressors; however,
   identifying their targets can facilitate the elucidation of their role
   in GBM development. Since a single miRNA regulates multiple genes, and
   a combination of miRNAs can co-modulate signaling pathways, the
   application of network-based approaches is inevitable to understand the
   contribution of an exact miRNA to GBM development. Our analysis was
   based on experimentally validated miRNA targets, applying the miRNet
   tool. A miRNA-centric network was constructed, including direct
   miRNA–target gene interactions and target gene-coded protein–protein
   interactions (PPI). In this heterogeneous network, hsa-miR-15a-5p has
   the highest degree value (260), followed by hsa-miR-424-5p (185) and
   hsa-miR-21-5p (131), reflecting their importance in the interactome
   network. Regarding the target genes, Nuclear FMR1 Interacting Protein 2
   (NUFIP2) was regulated by the highest number of miRNAs and had the most
   interacting protein partners (degree value 20), followed by Zing Finger
   Protein 460 (ZNF460; degree value 19), Cyclin D1 (CCND1; degree value
   15), and Cyclin-Dependent Kinase 6 (CDK6; degree value 15). The minimum
   network built from the miRNA–target gene and protein–protein
   interactions—as generated by the miRNet tool—is shown in [67]Figure 4.
   The red squares represent the upregulated miRNAs, while the green
   squares represent the downregulated miRNAs.

Figure 4.

   [68]Figure 4
   [69]Open in a new tab

   The minimum network of the differentially expressed (DE) miRNAs and
   their experimentally validated target genes. The network was
   constructed using the miRNet tool. Color code: A red square indicates
   an upregulated miRNA, while a green square indicates a downregulated
   miRNA; the blue spots are proteins.

2.1.5. Gene Ontology (GO) and Pathway Enrichment Analysis of miRNA Targets

   Using the GO and Kyoto Encyclopedia of Genes and Genomes (KEGG)
   database options of the NetworkAnalyst 3.0 tool, functional enrichment
   and pathway analysis were performed based on the PPIs of experimentally
   validated miRNA targets. As shown in [70]Figure 5a, the protein
   products of the target genes of the upregulated miRNAs were enriched in
   biological processes such as positive regulation of the metabolic
   process (p = 2.11 × 10^−30), regulation of signal transduction (p =
   4.54 × 10^−26), and regulation of apoptotic process (p = 9.54 ×
   10^−26). In contrast, the proteins of the downregulated miRNAs were
   involved in the regulation of developmental processes (p = 1.37 ×
   10^−33), negative regulation of metabolic processes (p = 5.21 ×
   10^−33), and positive regulation of transcription (p = 1.4 × 10^−30)
   ([71]Figure 5b).

Figure 5.

   [72]Figure 5
   [73]Open in a new tab

   GO Biological Process-based functional enrichment annotation of (a)
   upregulated and (b) downregulated miRNA-related target genes using the
   NetworkAnalyst 3.0 tool. The significance of the detected biological
   processes is characterized by their FDR and −log[10] p-values. The size
   of the dots is proportional to the number of genes included in the
   given process.

   According to the KEGG database, we have found that the target genes of
   upregulated miRNAs are implicated in diverse cancer-related pathways,
   including the pathways in cancer (p = 1.22 × 10^−21), cell cycle (p =
   5.33 × 10^−189), and FoxO signaling pathway (p = 3.22 × 10^−16)
   ([74]Figure 6a), while the downregulated miRNAs are involved in the
   regulation of the AGE-RAGE signaling pathway (p = 1.3 × 10^−18) and
   proteoglycan in cancer (p = 2.25 × 10^−16) ([75]Figure 6b).

Figure 6.

   [76]Figure 6

   [77]Figure 6
   [78]Open in a new tab

   The KEGG pathway-based functional pathway enrichment analysis of target
   genes (a) of upregulated and (b) downregulated miRNAs using the
   NetworkAnalyst 3.0 tool. The significance of the detected KEGG pathways
   is specified by their FDR and −log10 p-values. The size of the dot
   reflects the number of genes included in the given pathway.

2.1.6. Validation of Differentially Expressed (DE) miRNAs by RT-qPCR in
Tissue Samples

   To validate the results obtained by NGS, three upregulated miRNAs
   (hsa-miR-196a-5p (log[2]FC = 5.6), hsa-miR-21-3p (log[2]FC = 4.39), and
   hsa-miR-10b-3p (log[2]FC = 3.66)) and two downregulated miRNAs
   (hsa-miR-383-5p (log[2]FC = −6.33) and hsa-miR-490-3p (log[2]FC =
   −5.61)) were selected for an RT-qPCR analysis. For validation, total
   RNA was isolated from tissue samples from 30 GPs and 28 control samples
   from the peripheral tumor region of patients with low-grade (grade 1–2)
   glioma, which served as a control group. The number of women and men
   was equal in both groups. RT-qPCR was used to measure the relative
   expression of these miRNAs using hsa-miR-103a-3p as a reference miRNA
   [[79]15]. All measurements were performed in triplicate. According to
   the results of the RT-qPCR reactions using the Mann–Whitney U test for
   calculation, we confirmed that the expression of hsa-miR-196a-5p,
   hsa-miR-21-3p, and hsa-miR-10b-3p was significantly upregulated, while
   hsa-miR-383-5p and hsa-miR-490-3p showed significant downregulation
   compared to their expression in the control samples ([80]Figure 7).

Figure 7.

   [81]Figure 7
   [82]Open in a new tab

   Representation of normalized Ct values of significantly deregulated
   miRNAs. The p-values were calculated using the Mann–Whitney U test. ***
   p < 0.001. The p-value of all validated miRNAs is p < 0.01.

   We created ROC-AUC curves applying the normalized expression data
   resulting from the RT-qPCR measurements of hsa-miR-196a-5p,
   hsa-miR-21-3p, hsa-miR-10b-3p, hsa-miR-383-5p, and hsa-miR-490-3p in
   order to verify their diagnostic potential. The ROC-AUC values were
   0.96032, 0.97768, 0.99206, 0.9375, and 0.9648 in the case of
   hsa-miR-196a-5p, hsa-miR-21-3p, hsa-miR-10b-3p, hsa-miR-383-5p, and
   hsa-miR-490-3p, respectively. For the control group and GPs, the
   normalized Ct values were dichotomized by mapping the sensitivity
   values in relation to 1—specificity in the case of hsa-miR-196a-5p,
   hsa-miR-21-3p, hsa-miR-10b-3p, hsa-miR-383-5p, and hsa-miR-490-3p—to
   calculate optimal cut-off values. We observed the highest sensitivity
   value in the case of hsa-miR-383-5p (95%) and hsa-miR-490-3p (95%), a
   little bit lower value for hsa-miR-10b-3p (94%) and hsa-miR-21-3p
   (93.8%), while the lowest one is in the case of hsa-miR-196a-5p (88%).
   The sequence of specificity was different; the highest value belonged
   to hsa-miR-10b-3p (100%), followed by hsa-miR-21-3p (92.9%),
   hsa-miR-196a-5p (92%), hsa-miR-383-5p (95%), and hsa-miR-490-3p (85%)
   ([83]Figure 8).

Figure 8.

   [84]Figure 8
   [85]Open in a new tab

   (a) ROC (Receiver Operating Characteristics) curves with AUC (area
   under the curve) were created to detect and represent the sensitivity
   and specificity of hsa-miR-10b-3p, hsa-miR-196a-5p, hsa-miR-383-5p,
   hsa-miR-490-3p, and hsa-miR-21-3p selected for the validation of the
   NGS via RT-qPCR. (b) Calculated optimal cut-off point values of
   hsa-miR-10b-3p (normalized Ct value: 9.9), hsa-miR-196a-5p (normalized
   Ct value: 9.8), hsa-miR-383-5p (normalized Ct value: 9.8),
   hsa-miR-490-3p (normalized Ct value: 9.8), and hsa-miR-21-3p
   (normalized Ct value: 5.8). (c) Dot plot analysis of hsa-miR-10b-3p,
   hsa-miR-196a-5p, hsa-miR-383-5p, hsa-miR-490-3p, and hsa-miR-21-3p for
   GPs and control tissue samples.

2.2. Identification of Differently Expressed (DE) mRNAs in Tissue Samples of
GBM Patients and Control Group

2.2.1. Next-Generation Sequencing (NGS)

   To assess the dissimilarities in mRNA expression profiles between the
   tissue samples of GPs and control individuals, next-generation
   sequencing (NGS) was conducted using the Illumina NextSeq 500
   instrument (Illumina, San Diego, CA, USA). The RNA-Seq libraries were
   created from the same total RNA samples, which were isolated from
   resected tumor tissue of five GPs and the peripheral tumor region of
   five lower-grade (grade 1–2) glioma patients. These later samples were
   used as a control group.

2.2.2. Hierarchical Clustering with Heatmap and Principal Component (PCA)
Analysis

   In order to interpret the results of mRNA NGS, a hierarchical cluster
   analysis was performed on the NGS dataset of GPs and the control group.
   MRNAs were ranked based on their standard deviation across all samples,
   and the top 200 genes were used in hierarchical clustering. The
   expression profiles of the top 50 most variable mRNAs in samples of GPs
   and controls are depicted as a heatmap in [86]Figure 9.

Figure 9.

   [87]Figure 9
   [88]Open in a new tab

   Heatmap with hierarchical clustering dendrogram of mRNA expression in
   GBM and control samples based on the variation in mRNA expression. The
   expression profiles of the top 50 mRNAs are shown. Columns represent
   patient and control individuals (G—glioblastoma sample and C—control
   sample), and each row represents a single mRNA. The down- and
   upregulated mRNAs are labeled green and red, respectively.

   The distribution of mRNA expression values was visualized by Principal
   Component Analysis (PCA). A PCA plot using the first and second
   principal components is presented in [89]Figure 10. It is visible that
   based on their mRNA expression pattern, GBM samples form a single
   cluster and are clearly separated from the control samples by the first
   principal component, which represents 43%. This data distribution of
   the variance suggests that GBM biogenesis induces a drastic change in
   the expression of several mRNAs.

Figure 10.

   [90]Figure 10
   [91]Open in a new tab

   PCA of mRNA expression based on their expression profile. A clear
   separation is visible between the GBM samples and the control samples
   along the first principal component.

2.2.3. Differentially Expressed Genes (DEGs)

   According to the DESeq2 algorithm of the iDEP.96 web tool with an
   adjusted threshold of false discovery rate (FDR) < 0.05 and fold-change
   > 2, we detected 365 upregulated (log[2]FC ≥ 1) and 1225 downregulated
   (log[2]FC ≤ −1) mRNAs in the GPs compared to the control samples. The
   whole list of the deregulated mRNAs is presented in [92]Table S2. The
   five most strongly upregulated mRNAs were HOXD10 (log[2]FC = 7.06),
   SHOX2 (log[2]FC = 6.86), POSTN (log[2]FC = 6.09), TOP2A (log[2]FC =
   6.08), and HOXD11 (log[2]FC = 5.95). In contrast, the five most
   significantly downregulated mRNAs were RELN (log[2]FC = −8.56), GRIN1
   (log[2]FC = −8.29), UNC13C (log[2]FC = −8.04), OPALIN (log[2]FC =
   −7.94), and GABRA1 (log[2]FC = −7.89).

   To assess the influence of GBM on mRNA expression in comparison to
   control samples, a k-means cluster analysis was conducted to categorize
   mRNAs based on their expression values ([93]Figure 11a). The volcano
   plot of the expression values ([94]Figure 11b) demonstrates that GBM
   development results in a profound alteration of the mRNA transcriptome.

Figure 11.

   [95]Figure 11
   [96]Open in a new tab

   (a) The k-means clustering of differentially expressed (DE) mRNA genes
   was conducted using the DESeq2 algorithm. The down- and upregulated
   mRNAs are labeled green and red, respectively. (b) The volcano plot
   shows that GBM development leads to a massive change in the mRNA
   transcriptome. The down- and upregulated mRNAs are labeled blue and
   red, respectively.

2.2.4. Protein–Protein Interaction (PPI) Network Analysis of Deregulated
mRNAs

   The construction of minimal networks based on the 50 most significantly
   upregulated and downregulated mRNAs was carried out using the
   NetworkAnalyst 3.0 tool ([97]Figure 12). Many of the major
   hubs—considered to be key nodes with significant biological
   relevance—are proteins already known to be involved in tumorigenesis.
   It is an intriguing finding that the protein with the highest number of
   interactions in both analyses (degree values of 25 and 13 for the up-
   and downregulated networks, respectively) is ubiquitin C (UBC).
   Ubiquitination, depending on the residues involved in conjugation, can
   be related to the degradation of proteins, DNA repair, kinase
   modification, endocytosis, and regulation of the cell cycle and cell
   signaling pathways [[98]16,[99]17]. As shown in [100]Figure 12a,
   following UBC, the small ubiquitin-like modifier 2 (SUMO2) and
   Insulin-Like Growth Factor 2 mRNA Binding Protein 3 (IGF2BP3) have
   relatively high degree values in the mRNA-based upregulated
   protein–protein interaction network (degree values of 12 and 9,
   respectively), while synaptotagmin-1 (SYT1, degree value 7) and Amyloid
   Beta Precursor Protein (APP, degree value 5) follow UBC in the
   downregulated mRNA-based protein–protein interaction network
   ([101]Figure 12b).

Figure 12.

   [102]Figure 12
   [103]Open in a new tab

   Topology of the (a) 50 most upregulated mRNA-based and (b) 50 most
   downregulated mRNA-based PPI minimal networks using the NetworkAnalyst
   3.0 tool. Nodes indicate proteins; the size of the nodes corresponds to
   their degree centrality.

2.2.5. Gene Ontology (GO) and Pathway Enrichment Analysis of mRNA Molecules

   Applying the GO and Kyoto Encyclopedia of Genes and Genomes (KEGG)
   database options of the NetworkAnalyst 3.0 tool, PPI functional
   enrichment and pathway analysis were performed using deregulated mRNAs.
   The GO Biological Process enrichment analysis demonstrates that the
   upregulated proteins were implicated in the positive regulation of
   nucleobase-containing compound metabolic processes (p = 1.27 × 10^−22),
   cell cycle (p = 2.28 × 10^−21), positive regulation of RNA metabolic
   processes (p = 2.66 × 10^−20), and chromosome organization (p = 3.35 ×
   10^−14) (see [104]Figure 13a). In contrast, the downregulated proteins
   were enriched in processes such as neuron development (p = 4.5 ×
   10^−23), neuron projection development (p = 5.48 × 10^−22), neuron
   differentiation (p = 1.3 × 10^−19), and synaptic transmission (p = 1.43
   × 10^−17) ([105]Figure 13b).

Figure 13.

   [106]Figure 13
   [107]Open in a new tab

   GO Biological Process-based functional enrichment annotation of (a) 50
   most upregulated and (b) 50 most downregulated mRNAs using the
   NetworkAnalyst 3.0 tool. The significance of the detected biological
   processes is indicated by their FDR and −log10 p-values. The size of
   the dots is proportional to the number of genes included in the given
   process.

   According to the KEGG database, we found that the upregulated mRNAs
   were involved in different cancer-related pathways, like cell cycle (p
   = 8.31 × 10^−22), cellular senescence (p = 1.69 × 10^−14), pathways in
   cancer (p = 8.33 × 10^−12), transcriptional misregulation in cancer (p
   = 3.06 × 10^−11), and the IL−17 signaling pathway (p = 5.38 × 10^−12)
   ([108]Figure 14a), while focal adhesion (p = 1.13 × 10^−12), pathways
   in cancer (p = 1.23 × 10^−11), long-term potentiation (p = 1.27 ×
   10^−11), and the ErbB signaling pathway (p = 1.85 × 10^−10) showed
   downregulation ([109]Figure 14b). It is also important to note that
   DEGs were highly enriched in pathways that are involved in the
   development of several cancer types ([110]Figure 14).

Figure 14.

   [111]Figure 14
   [112]Open in a new tab

   The KEGG pathway-based enrichment annotation of (a) 50 most upregulated
   and (b) 50 most downregulated mRNA-based analyses using the
   NetworkAnalyst 3.0 tool. The significance of the detected KEGG pathways
   is specified by their FDR and −log10 p-values. The size of the dot
   reflects the number of genes included in the given pathway.

2.2.6. Validation of Differentially Expressed (DE) mRNAs by RT-qPCR in Tissue
Samples

   The same total RNA samples from GPs and the control group were used to
   validate the mRNA NGS results. Seven upregulated mRNAs (E2F2 (log[2]FC
   = 3.59), HOXD13 (log[2]FC = 3.69), VEGFA (log[2]FC = 4.3), CDC45
   (log[2]FC = 4.31), AURKB (log[2]FC = 4.6), HOXC10 (log[2]FC = 4.9), and
   MYBL2 (log[2]FC = 5.73)) and seven downregulated mRNAs (FABP6 (log[2]FC
   = −2.3), PRLHR (log[2]FC = −4.37), NEUROD6 (log[2]FC = −5.72), CBLN1
   (log[2]FC = −6.16), HRH3 (log[2]FC = −6.39), HCN1 (log[2]FC = −7.36),
   and RELN (log[2]FC = −8.5)) were chosen to check their expression in
   GPs by RT-qPCR. GAPDH was used as an internal control for normalizing
   the results [[113]18]. In the case of the selection of deregulated
   mRNAs for the validation procedure, we primarily aimed to cover the
   entire significantly up- and downregulated expression spectrum. The
   list of primer sequences is presented in [114]Table S3. Each experiment
   was performed in triplicate. Based on the results of RT-qPCR validation
   using the Mann–Whitney U test, the upregulation of E2F2, HOXD13, VEGFA,
   CDC45, AURKB, HOXC10, and MYBL2 ([115]Figure 15a) and the
   downregulation of FABP6, PRLHR, NEUROD6, CBLN1, HRH3, HCN1, and RELN
   were confirmed ([116]Figure 15b).

Figure 15.

   [117]Figure 15
   [118]Open in a new tab

   Representation of normalized Ct values of (a) significantly upregulated
   and (b) significantly downregulated mRNAs. The p-values were calculated
   using the Mann–Whitney U test. * p < 0.05; ** p < 0.01; *** p < 0.001.
   In the case of the upregulated mRNAs, p < 0.001. Related to the
   downregulated mRNAs, p < 0.001 is characteristic for HRH3, CBLN1, RELN,
   and PRLHR, while p < 0.01 is calculated in case of HCN1 and NEUROD6 and
   p < 0.05 in case of FABP6.

2.3. The Correlation Between miRNA and mRNA Expression Determined by
Next-Generation Sequencing (NGS)

   Finally, we analyzed miRNA–mRNA interactions using our miRNA and mRNA
   expression datasets. miRNA targets were identified using the miRTarBase
   and miRTargetLink 2.0 databases. Note that only experimentally
   validated miRNA target gene interactions were included in our study.
   [119]Table 1 shows the inverse correlation between the expression of
   miRNAs and their experimentally validated target mRNAs in our NGS
   dataset. We hypothesize that miRNAs that show significant up- or
   downregulation in GPs may regulate the expression of genes that are
   involved in the cell cycle (AURKB, CDC45, and CDK6), cell proliferation
   (EGFR and VEGFA), and angiogenesis (VEGFA) that support tumor growth.
   In this context, upregulated transcription factors, such as E2F
   transcription factor 2 (E2F2) and MYB proto-oncogene-like 2 (MYBL2),
   regulate the transcription of genes involved in the cell cycle, cell
   differentiation, and cell proliferation. Other genes (AJAP1, MMP9,
   POSTN, and STC2) promote metastasis formation by regulating adhesion or
   migration. In addition, the upregulated Homeobox C10 (HOXC10) is
   involved in the transcription of genes that enhance migration capacity.
   The tumor microenvironment may also be influenced by the regulation of
   tumor-associated macrophages (through LTBP-1 and POSTN).

Table 1.

   Correlation between the deregulated miRNAs and their experimentally
   validated target genes found in our NGS dataset.
   Regulated Gene Biological Process Expression Status
   (Up/Down) miRNA Expression Status
   (Up/Down) References