In the presented study, serum samples were collected in the Clinical Hospital in Bialystok by Biobank at the Medical University of Bialystok, with high standards of strict biobanking procedures for multi-omics studies (Niklinski et al., 2017; Michalska-Falkowska et al., 2023). The study group consisted of patients qualified for surgical treatment of adult-type diffuse gliomas at different stages of development (WHO, 2021 grade 2–4). All patients were grouped into 24 lower-grade gliomas CNS5 WHO grade 2 and 3 and 78 WHO grade 4 (Table 1). The serum was collected before any treatment of cancer patients was started. The serum samples were collected, centrifuged and stored at −80°C. A group of 102 individuals was qualified for further analysis. The exclusion criteria for cases were patients diagnosed with other types of cancer than an adult-type diffuse glioma, chemotherapy or radiotherapy before serum collection and no cancer history. All participants provided written informed consent and received detailed information on the study and associated risks before enrolment. This study was approved by the Bioethics Committee of the Medical University of Bialystok, Poland (approval numbers: R-I-002/357/2014, R-I-002/600/2019, and APK.002.171.2021) and was performed according to the principles of the Declaration of Helsinki.

Based on similar previous experiments and pilot data, we have calculated the minimal number of samples per experimental group (grade 4 or grade 2–3) to detect two-fold differences in relative expression levels between groups at the true positive detection powers of 80% and 90% (Hart et al., 2013). We have used the RNASeqPower R package to apply the statistics data covering obtained real counts and coefficients of variations per group. For the NanoString nCounter miRNA data, we have estimated that to obtain 80% power, we would need nine samples per group, whereas to obtain a high power of 90%, we would need 12 samples. Finally, our groups for NanoString analyses consisted of 78 (grade 4) and 24 (grade 2–3) samples, thus allowing for more than 90% power in any of the comparisons performed.

RNA isolation with miRNA fraction serum samples was performed using the miRNeasy Serum/Plasma Advanced Kit (Qiagen, Germany) according to the manufacturer’s instructions. A total of 102 samples were analysed using the nCounter®Analysis System (NanoString Technologies, WA, United States) and the nCounter Human v3 miRNA Panel. Briefly, as input material, 3 ng of isolated miRNA was used. Unique RNA tags were ligated onto the 3′ end of each mature miRNA, followed by an overnight hybridisation (65°C) to nCounter Reporter and Capture probes. After hybridisation, samples were placed into the nCounter Prep Station for sample purification and target/probe complexes immobilisation on the cartridge. For each assay, a high-density scan (555 fields of view) was performed on the nCounter Digital Analyzer (NanoString Technologies, United States) to count individual fluorescent barcodes and quantify target miRNA molecules present in each sample.

The seven DE miRNAs and two reference miRNAs were profiled using the miRCURY LNA SYBR Green PCR Kit (Qiagen, Germany). Reference miRNAs with stable expression across all samples, has-miR-103-3p and has-199b-5p, were selected based on the NanoString data using the NormFinder algorithm (Andersen et al., 2004). The miRCURY LNA RT Kit (Qiagen, Germany) was used for the reverse transcription reaction for miRNA assays. Assay IDs are assembled in Supplementary Table S1. Quantitative real-time PCR (RT-qPCR) was performed with miRCURY LNA SYBR^® Green PCR Kits (Qiagen, Germany) with specific commercial primers, optimized with LNA technology to enable sensitive and specific miRNA quantification. For STAT3, FOS and TLR4 mRNA analysis, we used the Transcription High Fidelity cDNA Synthesis Kit (Roche, Switzerland) and KiCqStart^® SYBR^® Green qPCR ReadyMix™ (Sigma-Aldrich, Germany). The following primers were used: STAT3 FW-CTTTGAGACCGAGGTGTATCACC; STAT3 RV-GGTCAGCATGTTGTACCACAGG; FOS FW-GCCTCTCTTACTACCACTCACC; FOS RV-AGATGGCAGTGACCGTGGGAAT; TLR4 FW-CCCTGAGGCATTTAGGCAGCTA; FOS RV-AGGTAGAGAGGTGGCTTAGGCT; GAPDH FW GTC​TCC​TCT​GAC​TTC​AAC​AGC​G; GAPDH RV-ACCACCCTGTTGCTGTAGCCAA. The temperature profile of the qPCT reaction was as follows: 2 min at 95°C and 45 cycles: 10 s at 95°C and 60 s at 56°C (for miRNA) and 60°C (for mRNA). Amplification was performed using LightCycler 480 (Roche, Switzerland). Subsequently, PCR threshold cycles (Ct) of the tested miRNA/mRNA and reference miRNA/mRNA were determined for the tested samples and the calibrator. The relative expression for each miRNA/mRNA was calculated with PCR efficiency correction (Pfaffl, 2001). Efficiency (E) was calculated from the slopes of the calibration curve according to the equation: E = 10 (−1/slope). Reactions with amplification efficiency below 1.6 were removed. The relative expression ratio of a target miRNA was computed based on its PCR efficiencies (E) and the Ct value difference (Δ) of unknown group samples (test) versus the control group (Δ Ct control-test). The relative calculation was based on the MEAN Ct of the experimental group.

The identification of DE miRNA target genes was performed using Ingenuity Pathway Analysis Software (Krämer et al., 2014) (IPA, Qiagen Inc. https://qiagenbioinformatics.com/products/ingenuity-pathway analysis), mirDB database (Chen and Wang, 2020) (https://mirdb.org), and TagetScanHuman 8.0 database (Agarwal et al., 2015) (https://www.targetscan.org/vert_80/). Gene Ontology Biological Process and KEGG Pathway analyses were conducted with the ClusterProfile R package (Yu et al., 2012). IPA was used to perform the core analysis to identify canonical pathways. Over-representation analyses used hypergeometric tests under α = 0.05 (p-values corrected with FDR). To construct a PPI network, we used STRING (https://string-db.org/; v.11.5). Genes with a confidence score ≥0.4 were chosen to build a network model visualised by Cytoscape v.3.9.1. Nine topological algorithms in plug-in cytoHubba (Chin et al., 2014), consisting of “MCC,” “MNC,” “Degree,” “Bottle Neck,” “EcCentricity,” “Closeness,” “Stress,” and “Radiality” were selected to identify the hub genes in PPI analysis. Analyses of functional interaction networks were based on the HumanNet v3 platform by applying the HumenNet-FN (functional gene network) mode (https://www.inetbio.org/humannet/) (Kim et al., 2022).

The discovery and the analysis of miRNA co-expression modules in the patient’s tumour samples based on NanoString data were performed with Weighted Gene Coexpression Network Analysis (WGCNA) (Langfelder and Horvath, 2008) by application of the WGCNA and CEMiTool R libraries (Russo et al., 2018). Outliers were removed based on standardised connectivity. Pearson correlation was used in the signed network construction with soft threshold selection (R2 > 0.8). Modules were merged based on a high eigengene similarity correlation threshold (0.85). Gene modules to trait relationships were evaluated with Spearman correlation. Within each WGCNA module, an interaction network was analysed using CEMiTool. Interaction networks integrate the miRNA co-expression results with miRNA-target gene annotations based on miRTarBase 8.0 (Huang et al., 2020). Database miRTarBase provides comprehensive information on experimentally validated miRNA-target interactions from various literature data and other databases that have been data mined and manually curated. The top ten most connected genes were identified as network hubs and are indicated with text labels and coloured based on the interaction type (miRNA co-expression pattern–blue; miRNA-target gene interaction–red; miRNA co-expression pattern and miRNA-target gene interaction–yellow). The size of each node is proportional to its degree, as indicated in the plot legend.

The normalised NanoString counts and RT-qPCR data were used for model development. The attributes were selected using Waikato Environment for Knowledge Analysis (WEKA) version 3.8.3. (c) 1999–2018 The University of Waikato, Hamilton, New Zealand. The InfoGainAttributeEval algorithm was used to select the best classifiers. Feature selection via information gain using InfoGainAttributeEval was based on the calculation of decreasing entropy by adding attributes and selection attributes that most strongly reduce entropy. This method was carried out on the testing set using the LOOCV (Wainer and Cawley, 2021; Rafało, 2022). A multivariate logistic regression model was built using repeated (n = 3) k-fold-cross-validation (k = 10) in R version 3.6.1 [R Core Team (2013) R: A language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org] (Kuhn, 2008). To validate the model, the ROC curve and the AUC were calculated by the pROC package (Robin et al., 2011). The confusion matrix, including information about TP (true positives), TN (true negatives), FP (false positives), and FN (false negatives) has been prepared. The evaluation of model classification has been based on the testing dataset. The choice of miRNAs was evaluated using the independent, external and publicly available dataset GSE112462 from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE112462), which contains miRNA profiles of 28 glioma patients (18 patients with grade 2–3 and 10 patients with grade 4) (Drusco et al., 2018). This external data was generated using NanoString human miRNA panel NS_H_MIR_V3A. A logistic regression model was developed for this data according to the above workflow.

Statistical analyses were performed using GraphPad Prism 9 (v.9.3.1) software. The Wilcoxon rank sum test was used to investigate the differences in BMI, age and size of tumours between the grade 4 patients’ group and grade 2–3. Raw miRNA data were analysed using nSolver Software v. 4.0 (NanoString Technologies, WA, United States). Code-set content normalisation was performed relative to the ligation controls for technical variations. Ratios were calculated by specifying the lower-grade samples as a baseline. According to Benjamin-Hochberg, correction for multiple testing was performed with a False Discovery Rate (FDR). The Pearson correlation coefficient (r) was used to estimate the correlation between the identified DE miRNAs and clinical parameters.

The expression levels of 798 miRNAs were quantified in the study group consisted of patients qualified for surgical treatment of adult-type diffuse gliomas of different stages of development (according of WHO, 2021 classification), including 78 grade 4 gliomas and 24 grade 2 and 3. Altogether, seven unique miRNAs showed significant differences in counts between grade 4 gliomas and grade 2–3 gliomas serum samples with fold change (FC) > |1.5| and false discovery rates (FDRs) < 0.05 (Figure 1A). Table 2 summarises FCs and FDRs for all seven differentially expressed (DE) miRNAs. Five miRNAs were upregulated: miR-630, miR-362-3p, miR-320e, miR-4454+miR-7975 (combined because the mature sequence of miR-7975 differs from miR-4454 by only one base), whereas two miRNAs were downregulated: miR-1253 and miR-6721-5p. To validate the diagnostic usefulness of identified miRNAs, we performed RT-qPCR validation. All miRNAs were detectable by quantitative PCR, and the expression profiles were similar to the ones generated using the NanoString platform (Figure 1B). Next, we used Pearson’s correlation to determine a statistically significant relationship between miRNA expression levels and clinical features such as grade, Karnofsky Performance Score (KPS) and size of tumours (Figure 1C). The analysis indicated a positive correlation between miR-362-3p expression and a negative correlation between miR-6721-5p expression and gliomas’ grade with a coefficient of 0.45 and −0.43, respectively (p-value = 1.87 × 10⁻⁶ and p-value = 6.79 × 10⁻⁶). It means that the expression of miR-362-3p increased in the case of higher-grade gliomas, while the change in miR-6721-5p expression level was opposite to the change in the grade of gliomas.

The receiver operating characteristic (ROC) curve analysis was conducted to evaluate the diagnostic value of the identified miRNAs as biomarkers for binary classification (glioma 4 vs. glioma 2–3) (Figure 2A). The area under the ROC curve (AUC) for all miRNAs was higher than 0.700 (p < 0.0001). The highest AUC values were found for miR-362-3p, miR-630, and miR-320e (0.782, 0.769, and 0.748, respectively). Corresponding AUC values, confidence intervals (CI), sensitivity and specificity for cut-off points were calculated and are shown in Supplementary Table S3.

To select the combination of miRNAs with the strongest classification power for discrimination between grade 4 of adult-type diffuse glioma patients and lower grades, we applied the Information Gain (InfoGainAttributeEval) selection method in WEKA software to analyse obtained differentially expressed miRNA data. The Information Gain method with the Ranker Search method is based on the calculation of decreasing entropy by adding attributes, and Correlation-based feature selection prioritises uncorrelated features (Wainer and Cawley, 2021; Rafało, 2022). The three attributes by InfoGainAttributeEval Algorithms with the Ranker Search Method were miR-362-3p, miR-6721-5p, and miR-320e.

To assess the diagnostic values of the selected miR-362-3p, miR-6721-5p, and miR-320e, multivariate logistic regression was applied to develop the diagnostic models of the miRNAs under combination conditions. K-fold cross-validation (K = 10) was used to develop stable models. Based on the normalised NanoString data, four models were derived: model 1 was built on the expression of miR-362-3p and miR-6721-5p as independent variables; model 2 included miR-362-3p and miR-320e; model 3 was based on the expression of miR-6721-5p and miR-320e, while model 4 included all three variables, miR-362-3p, miR-320e, and miR-6721-5p (Table 3). To assess the specificity and sensitivity of the four models, ROC curves were drawn, as shown in Fig. XX. The AUC value for model 1 was 0.894, with the diagnostic sensitivity at 94.8% and the specificity at 69.6% (Figure 2B). For model 2, AUC was 0.851, sensitivity 83.1% and specificity 63.6 (Supplementary Figure S2A). For model 3, AUC was 0.831, sensitivity 89.6% and specificity 24.2% (Supplementary Figure S2B); for model 4, AUC was 0.862, sensitivity 92.8% and specificity 76.1% (Supplementary Figure S2C). The data suggested the best classification accuracy was achieved by model 1 (Table 4). In addition, the combination model showed a higher sensitivity and better specificity than individual miRNA.

Next, RT-qPCR analysis was performed to further validate the value of circulating miR-362-3p and miR-6721-5p as prediction markers for different grated of gliomas. The expression levels of miR-362-3p and miR-6721-5p were normalised to the expression of reference miRNAs, miR-103-3p and miR-199-5p. A diagnostic classification model 5 was developed based on the normalised RT-qPCR data. As shown in Figure 2C; Table 5, model 5 provided very good discrimination between grade 4 glioma and lower grades, with the AUC value of 0.867, with a sensitivity of 82.9% and a specificity of 92.8%.

Next, further validation of the optimised miR-362-3p and miR-6721-5p signature was performed using an external dataset (GSE112462, NanoString human miRNA panel NS_H_MIR_V3A), which contains miRNA profiles of 28 glioma patients (18 patients with grade 2–3 and 10 patients with grade 4). We found that model 6, developed based on miR-362-3p and miR-6721-5p expression data in the GEO dataset, efficiently discriminated between grade 4 glioma patients and lower grade glioma individuals with AUC value of 0.842, sensitivity 66.7% and specificity 83.3% (Figure 2C; Table 6).

To reveal the biological function of the DE miRNAs, first, we identified 5,657 putative DE miRNAs target genes using the Ingenuity Pathway Analysis (IPA) (Krämer et al., 2014), 3,493 target genes via mirDB (Chen and Wang, 2020) and 11,944 target genes via TagetScanHuman 8.0 (Agarwal et al., 2015) (Figure 3A). For further functional studies, we used 2,070 DE miRNA target genes that overlapped between the above databases. Using Cytoscape v.3.9.1 (Shannon et al., 2003), we constructed the protein-protein interaction (PPI) network, including 2,067 nodes and 12,250 edges, to identify the critical hub genes among the above miRNA target genes. The publicly available NetworkAnalyst platform (Zhou et al., 2019) was used to conduct “Zero order” interaction network analysis in a layout format using a force atlas to visualise the network, which included 514 nodes and 876 edges (Supplementary Figure S1). Hub genes in the network were ranked using eight topological analysis methods, including both local- and global-based algorithms from Cytoscape software’s cytoHubba plugin Pole (Chin et al., 2014). We found that STAT3, FOS, TLR4, CPLX1, CPLX2, STX1A, STX1B, VAMP2, MAPK3, and VEGFA scores ranked in the top (Figure 3A and Supplementary Table S2). The highly connected hub gene STAT3 encodes signal transducer and activator of transcription. It is a multifunctional transcription factor involved in many biological functions ranging from initiating malignant transformation to tumour invasion, migration, metastasis and angiogenesis (Tolomeo and Cascio, 2021).

Next, we investigated whether and how these 2070 prioritised DE miRNA target genes could correlate with the stage of glioma. To this end, we used the HumanNet v3 platform (Kim et al., 2022). First, juxtaposed DE miRNAs target genes identified in our cohort with genes related to “adult glioblastoma” from the publicly available DisGeNET v7.0 database (Piñero et al., 2020). The 305 gene targets that overlapped between 2070 DE miRNAs target genes and 2,528 putative target genes related to “adult glioblastoma” were analysed in terms of known interactions using HumanNet v3 platform. Network analysis (Figure 3B) was performed based on HumanNet (functional gene network) which includes co-functional links (given by co-expression, co-essentiality, pathway database, protein domain profile associations, gene neighbourhood, and phylogenetic profile association) and protein-protein interactions. The top guide gene within this network was also STAT3 (score = 42), and the other guide genes within the top scores included TP53, MAPK3, SP1, PML, BCL2L1, E2F1, RELA, IGF1R, and CHEK1. Within the obtained network, we also identified downstream candidate genes that could be functionally connected to the 305 input guide genes. The top scores comprised STAT3 (score = 48.5), BAX, MDM2, EGFR, PRKACA, PRKACB, SRC, MAPK1, MAPK8, and CTNNB1.

In addition, we evaluated the correlation between the expression of the top 3 hub genes, STAT3, FOS and TLR4, in the grade 4 glioma tissue and model 1 miRNAs, miR-362-3p and miR-6721-5p, in serum. Expression of all these genes increased in tumour samples (Figure 4A). The results, as shown in Figure 4B, demonstrate that the correlation between the expression of miR-362-3p and the expression of STAT3, FOS, and TRL4 was moderately positive (0.56, 0.47 and 0.38, respectively) with a significant statistical value of p < 0.0001. Conversely, a moderate negative correlation was observed between miR-6721-5p and STAT3, FOS and TRL4 (0.52, 0.36 and 0.37, respectively). The correlations were also statistically significant (p < 0.0001). It is worth emphasising that the highest value for the correlation coefficients was observed for both miRNAs and STAT3.

To explore the potential biological functions and mechanisms of the DE miRNAs, 2070 DE miRNAs target genes were analysed by KEGG, WikiPathways and GO-BP pathway enrichment analyses using the enrichR tool (Kuleshov et al., 2016), as well as the canonical pathway enrichment analyses by the IPA (Krämer et al., 2014) Figure 5. The KEGG enrichment analysis indicated that DE miRNA target genes were mainly enriched in the JAK/STAT signalling pathway (FDR-corrected p-value = 5.34 × 10⁻¹⁰), HIF-1 signalling pathway (FDR-corrected p-value = 8.21 × 10⁻¹⁰), signalling pathways of pluripotency of stem cells (FDR-corrected p-value = 2.89 × 10⁻⁹), glioma-related pathway (FDR-corrected p-value = 5.94 × 10⁻⁹), and MAPK signalling pathway (FDR-corrected p-value = 2.15 × 10⁻⁷) (Supplementary Figure S3). It is noteworthy that the top 3 dysregulated pathways identified by the WikiPathways tool included the brain-derived neurotrophic factor (BDNF) signalling pathway (FDR-corrected p-value = 3.78 × 10⁻⁹), glial cell differentiation (FDR-corrected p-value = 4.16 × 10⁻⁸), and glioblastoma signalling pathway (FDR-corrected p-value = 6.92 × 10⁻⁸). BDNF regulates cell growth, differentiation, migration and apoptosis in the nervous system and high-grade gliomas (Xiong et al., 2013). The IPA analysis showed the involvement of the following canonical pathways: molecular mechanisms of cancer (FDR-corrected p-value = 7.45 × 10⁻¹⁰), glioma signalling (FDR-corrected p-value = 3.53 × 10⁻⁸), STAT3 signalling (FDR-corrected p-value = 1.63 × 10⁻⁷), glioblastoma invasiveness signalling (FDR-corrected p-value = 4.47 × 10⁻⁷), and synaptogenesis signalling pathway (FDR-corrected p-value = 2.29 × 10⁻⁶). The results of GO enrichment analysis showed that DE miRNA target genes were mainly related to the regulation of transcription and peptidyl-serine phosphorylation.

To cluster DE miRNAs that are highly related to the clinical traits, we used the Weighted Gene Co-Expression Network (WGCNA) algorithm (Langfelder and Horvath, 2008) as a standard method to find cooperatively expressed mRNA or miRNA modules. The soft threshold value was set to 16, creating a scale-free system (Supplementary Figure S4). Modules were generated through dynamic tree cutting. After merging highly similar modules, we developed a total of 17 modules. Module M1 contained two DE miRNAs, miR-320e and miR-7975, whereas module M2 includes four DE miRNAs, miR-630, miR-362-3p, miR-6721-5p, and miR-4454. Additionally, we calculated the correlation between each module and clinical parameters. We identified that modules M1 and M2 were highly correlated with grade, KPS, and tumour diameter (Figure 6A). The modules M1 and M2 network graphs highlight significant miRNAs as network hubs, such as miR-320e in module M1, and miR-362-3p, miR-630, and miR-4454 in module M2 (Figure 6B). Gene set enrichment analysis revealed that modules M1 and M2 were positively associated with grade 4 gliomas (NES = 1.99 and NES = 2.16, respectively) and negatively associated with the lower grades. Finally, the results of the KEGG and canonical pathway enrichment analyses indicated that the miRNA target genes from modules M1 and M2 were mainly involved in the JAK-STAT signalling pathway (FDR-corrected p-value = 2.56 × 10⁻¹³), tumour microenvironment pathway (FDR-corrected p-value = 2.89 × 10⁻¹³), FoxO signalling pathway (FDR-corrected p-value = 1.08 × 10⁻⁹), glioblastoma multiforme signalling (FDR-corrected p-value = 3.11 × 10⁻¹³), and MAPK signalling pathway (FDR-corrected p-value = 4.26 × 10⁻⁸).