Circulating lncRNA- and miRNA-Associated ceRNA Network as a Potential Prognostic Biomarker for Non-Hodgkin Lymphoma: A Bioinformatics Analysis and a Pilot Study

Non-Hodgkin lymphoma (NHL) is characterized by a great variability in patient outcomes, resulting in the critical need for identifying new molecular prognostic biomarkers. This study aimed to identify novel circulating prognostic biomarkers based on an miRNA/lncRNA-associated ceRNA network for NHL. Using bioinformatic analysis, we identified the miRNA-lncRNA pairs, and using RT-qPCR, we analyzed their plasma levels in a cohort of 113 NHL patients to assess their prognostic value. Bioinformatic analysis identified SNHG16 and SNHG6 as hsa-miR-20a-5p and hsa-miR-181a-5p sponges, respectively. Plasma levels of hsa-miR-20a-5p/SNHG16 and hsa-miR-181a-5p/SNG6 were significantly associated with more aggressive disease and IPI/FLIPI scores. Moreover, we found that patients with risk expression profiles of hsa-miR-20a-5p/SNHG16 and hsa-miR-181a-5p/SNHG6 presented a higher risk of positive bone marrow involvement. Moreover, hsa-miR-20a-5p/SNHG16 and hsa-miR-181a-5p/SNHG6 pairs’ plasma levels were associated with overall survival and progression-free survival of NHL patients, being independent prognostic factors in a multivariate Cox analysis. The prediction models incorporating the ceRNA network expression analysis improved the predictive capacity compared to the model, which only considered the clinicopathological variables. There are still few studies on using the ceRNA network as a potential prognostic biomarker, particularly in NHL, which may permit the implementation of a more personalized management of these patients.


Introduction
Non-Hodgkin lymphomas (NHLs) are a heterogenous group of lymphoproliferative malignancies in which the majority of cases arise from B cells during different stages of normal B-cell differentiation [1]. The latest GLOBOCAN data from 2020 indicate that NHL represents the most common hematological malignancy worldwide, corresponding Moreover, additional studies are needed to further establish ncRNAs as potential new predictive and prognostic biomarkers for use in clinical practice in order to improve NHL patients' management. Therefore, in the present study, we analyzed the potential of miR-NAs and lncRNAs as prognostic biomarkers of NHL patients by investigating the plasma expression levels of has-miR-20a-5p ahashsa-miR-181a-5p and their respective lncRNA with which they form a ceRNA network. Hsa-miR-20a-5p is one of the components of the miR-17-92 cluster (comprising miR-17, miR-18a, miR-19a, miR-19b-1, miR-20a and miR-92-1), which was demonstrated to play a central role during the stages of B-cell development, and its deregulated expression was shown to have oncogenic potential [22][23][24]. In fact, B-cell-specific miR-17∼92 transgenic mice developed lymphomas with high penetrance which phenotypically resemble human lymphomas, including DLBCL, demonstrating the role of this cluster as a driver of B-cell lymphomagenesis [25]. Similarly, miR-181a-5p was also shown to be differentially expressed during the development stages of B cells; in particular, miR-181a-5p ectopic overexpression in common lymphoid progenitors results in increasing the total number of B cells, demonstrating its role in the fine-tuning of B-cell development [26][27][28].
The miRNA-targeted mRNAs, which were only validated by strong evidence methods, were retrieved using miRTarBase (https://bio.tools/mirtarbase (accessed on 3 June 2021)) [30]. The protein-protein interaction (PPI) networks were analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING) database. Construction and visualization of the protein interaction network of the selected target genes were realized using the STRINGapp of the Cytoscape software (v3.8.2, Cytoscape Consortium, San Diego, CA, USA)). STRING enrichment analysis tool was used to retrieve the functional enrichment analysis of Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathways. The enrichment results were filtered, and redundant terms were removed according to the Jaccard index. Cytoscape visualization software (https://cytoscape.org/, (accessed on 10 September 2021)) was used to construct the final interaction networks [31].

Study Population
The study included 113 patients diagnosed with B-cell NHL (high-grade lymphomas versus low-grade lymphomas), of Caucasian ethnicity, older than 18 years and without known familial cancer history, as described in a previous study [32,33]. Patients were admitted and treated at a Portuguese hospital between January 2016 and June 2020. Patients' clinical information is summarized in Table 1.

RNA Extractions and qPCR
Blood samples were obtained at baseline before starting therapy. Peripheral blood samples were centrifuged for preparation of platelet-free plasma (PFP) samples, as previously described [34]. Supernatant was aliquoted and stored at −80 • C until use. The GRS microRNA kit (Grisp ® , GRiSP Research Solutions, Porto, Portugal) was used to isolate the miRNA portion according to laboratory procedures, while Plasma/Serum RNA Purification Kit (Norgen ® , Biotek Corp. Thorold, ON, Canada)) was used for total RNA isolation, according to manufacturer's instructions. Using NanoDrop Lite spectrophotometer (Thermo Scientific ®, , Waltham, MA, USA), the RNA concentration and purity were determined. Taqman ® MicroRNA Reverse Transcription kit along with Taqman ® MicroRNA assay (Applied Biosystems ® , Waltham, MA, USA) were used to carry out the cDNA synthesis for the miRNAs. For the analysis of lncRNA expression, the cDNA synthesis was performed using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems ® , Waltham, MA, USA in accordance with manufacturer's instructions. Using qPCR, the miRNA and lncRNA expression levels were quantified using a StepOneTM qPCR Real-Time PCR machine, and the following reaction mix components: 1xTaqMan ® Gene Expression Master mix (Thermo Fisher Scientific ® ) and 1x probes TaqMan ® MicroRNA Assays (hsa-miR-20a-5p: 000580 and hsa-miR181a-5p: 000480) and TaqMan ® Noncoding RNA assays (SNHG16: Hs01598403_g1 and SNHG6: Hs00417251_m1) (Applied Biosystems ® , Waltham, MA, USA). MiRNAs expression levels were normalized to hsa-miR-16 (000391), and lncRNA expression levels were normalized to GAPDH (Hs99999905_m1) endogenous control. Each sample had two technical replicates. The amplification conditions were: a holding stage 95 • C for 20s, followed by 45 cycles of 95 • C for 1s and 60 • C for 20 s. Data analysis was performed using StepOneTM Sofware v2.2 (Applied Biosystems ® , Waltham, MA, USA), with the same baseline and threshold set for each plate, to generate threshold cycle (Ct) values for all the genes in each sample.

Statistical Analysis
Statistical analysis was performed using SPSS (version 26.0; IBM Company, Chicago, IL, USA) and GraphPad Prism software (version 7.0; San Diego, CA, USA). Student's t-test or Mann-Whitney U test was used to evaluate statistical differences in the normalized expression (−∆Ct) of the miRNAs and lncRNAs among the different groups. Additionally, the 2 −∆∆Ct method (Livak method) was used to calculate the relative changes in gene expression between the different groups.
Receiver operating characteristic (ROC) analysis was performed to assess the prognostic accuracy, and the AUC was calculated. Overall survival (OS) and progression-free survival (PFS) curves were generated using Kaplan-Meier method and compared using the log-rank test. OS time was determined from the date of diagnosis to the date of mortality or the last follow-up. PFS time was determined to extend from the date of diagnosis to the date of disease progression, recurrence, mortality or last follow-up. Cox regression was used to analyze the prognostic value of the miRNAs/lncRNAs expression levels on the progression-free and overall survival. p < 0.05 was considered statistically significant.

lncRNA-miRNA-mRNA Network Construction
LncRNAs have the ability to sponge a variety of miRNAs to inhibit their regulatory effect on target mRNAs, creating a ceRNA network [18]. The StarBase database was used to identify the lncRNAs that target hsa-miR-20a-5p and hsa-miR-181a-5p ( Figure 1). The analysis identified 40 lncRNAs targeting hsa-miR-20a-5p and 35 lncRNAs targeting hsa-miR-181a-5p. For the subsequent analysis, we filtered the results by gene type=processed transcript and gene name = SNHG to obtain the respective lncRNAs from the biotype small nucleolar RNA host gene (SNHG), an emergent class of lncRNAs that have been involved in the induction of proliferation, cell cycle progression, invasion and metastasis of cancer cells, making this class of transcripts a viable and attractive biomarker for cancer development and progression [33]. effect on target mRNAs, creating a ceRNA network [18]. The StarBase database was used to identify the lncRNAs that target hsa-miR-20a-5p and hsa-miR-181a-5p ( Figure 1). The analysis identified 40 lncRNAs targeting hsa-miR-20a-5p and 35 lncRNAs targeting hsa-miR-181a-5p. For the subsequent analysis, we filtered the results by gene type=processed transcript and gene name = SNHG to obtain the respective lncRNAs from the biotype small nucleolar RNA host gene (SNHG), an emergent class of lncRNAs that have been involved in the induction of proliferation, cell cycle progression, invasion and metastasis of cancer cells, making this class of transcripts a viable and attractive biomarker for cancer development and progression [33]. (c) (d) Figure 1. In silico analysis of the lncRNAs targeting hsa-miR-20a-5p and hsa-miR-181a-5p. (a) lncRNAs that target hsa-miR-20a-5p according to StarBase database analysis; (b) details about the binding site of hsa-miR-20a-5p on SNHG16, predicted by StarBase database; (c) lncRNAs that target hsa-miR-181a-5p according to StarBase database analysis; (d) details about the binding site of hsa-miR-181a-5p on SNHG6, predicted by StarBase database.
Subsequently, the mirTarBase database, the largest known online database of validated miRNA:mRNA interactions, was used to identify target mRNAs of hsa-miR-20a-5p and hsa-miR-181a-5p. According to miRTarBase analysis, we filtered the list of mRNAs by Homo Sapiens species and retrieved the interactions that were only validated with strong evidence methods (Western blot, qRT-PCR or luciferase assay), which are listed in Figure 2. In order to perform the functional annotation and enrichment analysis, we analyzed the validated targets with the STRINGapp Protein Query from Cytoscape software. The identified targets were filtered into a protein-protein interaction (PPI) network with 63 nodes and 270 edges for hsa-miR-20a-5p targets and 66 nodes and 256 edges for the hsa-miR-181a-5p targets, presenting a significant enrichment (p = 1 × 10 −16 and p = 1 × 10 −16 , respectively). We also applied Markov clustering (MCL), which resulted in the clustering of the proteins according to their STRING interaction score (Figure 2b,d). Finally, we integrated the data from hsa-miR-20a-5p and hsa-miR-181a-5p identified targets and constructed the complete lncRNA-miRNA-mRNA network ( Figure 3).
For the functional enrichment analysis, we used a false discovery rate (FDR) threshold of p < 0.01, and the redundant terms were eliminated using a redundancy cutoff of 0.5, Figure 1. In silico analysis of the lncRNAs targeting hsa-miR-20a-5p and hsa-miR-181a-5p. (a) lncRNAs that target hsa-miR-20a-5p according to StarBase database analysis; (b) details about the binding site of hsa-miR-20a-5p on SNHG16, predicted by StarBase database; (c) lncRNAs that target hsa-miR-181a-5p according to StarBase database analysis; (d) details about the binding site of hsa-miR-181a-5p on SNHG6, predicted by StarBase database.
Subsequently, the mirTarBase database, the largest known online database of validated miRNA:mRNA interactions, was used to identify target mRNAs of hsa-miR-20a-5p and hsa-miR-181a-5p. According to miRTarBase analysis, we filtered the list of mRNAs by Homo Sapiens species and retrieved the interactions that were only validated with strong evidence methods (Western blot, qRT-PCR or luciferase assay), which are listed in Figure 2. In order to perform the functional annotation and enrichment analysis, we analyzed the validated targets with the STRINGapp Protein Query from Cytoscape software. The identified targets were filtered into a protein-protein interaction (PPI) network with 63 nodes and 270 edges for hsa-miR-20a-5p targets and 66 nodes and 256 edges for the hsa-miR-181a-5p targets, presenting a significant enrichment (p = 1 × 10 −16 and p = 1 × 10 −16 , respectively). We also applied Markov clustering (MCL), which resulted in the clustering of the proteins according to their STRING interaction score (Figure 2b,d). Finally, we integrated the data from hsa-miR-20a-5p and hsa-miR-181a-5p identified targets and constructed the complete lncRNA-miRNA-mRNA network ( Figure 3).   For the functional enrichment analysis, we used a false discovery rate (FDR) threshold of p < 0.01, and the redundant terms were eliminated using a redundancy cutoff of 0.5, which resulted in a total of 236 enriched terms among the KEGG, Reactome and GO categories for hsa-miR-20a-5p and 218 enriched terms for hsa-miR-181a-5p. The enriched terms for each category after being filtered are represented in Figure 4, which shows only the top 20 enriched terms for the GO Biological Process and KEGG analysis.
Analyzing the functionally enriched terms for the hsa-miR-20a-5p targets, we could observe that the majority of the targets, including the CCND1, CCND2, SMAD7/4, E2F3/E2F1, STAT3, CDKN1A, MYC, PTEN and MAPK family, are involved in the regulation of central cell functions, such as regulation of cell proliferation, cell cycle and transcription, according to GO terms analysis, highlighting their involvement in hematopoietic or lymphoid organ development. Among the functionally enriched terms in the KEGG and Reactome pathways, we could find major cancer signaling pathways, such as FoxO, MAPK, JaK-STAT and p53, all of which are involved in NHL development. In the functional analysis of hsa-miR-181a-5p targets through GO terms analysis, we observed that several targets are involved in the cellular protein modification process, regulation of protein kinase activity, cell proliferation and apoptosis, with indications of their involvement in hematopoietic or lymphoid organ development. In the KEGG and Reactome pathways, we could find p53, JAK-STAT and MAPK signaling pathways and cytokine signaling in the immune system, which is highly associated with miRNA deregulation, demonstrating that all central signaling pathways are involved in the development and progression of NHL.

miRNA and lncRNA Expression Levels in NHL Patients' Plasma Samples
We evaluated the plasma expression levels of hsa-miR-20a-5p and hsa-miR-181a-5p in 113 NHL patients using quantitative real-time PCR. The expression levels of hsa-miR-20a-5p were significantly higher (fold change = 2.83) in patients with high-grade lymphoma, whereas hsa-miR-181a-5p levels were significantly lower (fold change = 0.33) compared to low-grade lymphoma patients (p = 0.003 and 0.013, respectively) ( Figure 5).  Analyzing the functionally enriched terms for the hsa-miR-20a-5p targets, we could observe that the majority of the targets, including the CCND1, CCND2, SMAD7/4, E2F3/E2F1, STAT3, CDKN1A, MYC, PTEN and MAPK family, are involved in the regulation of central cell functions, such as regulation of cell proliferation, cell cycle and transcription, according to GO terms analysis, highlighting their involvement in hematopoietic or lymphoid organ development. Among the functionally enriched terms in the KEGG and Reactome pathways, we could find major cancer signaling pathways, such as FoxO, MAPK, JaK-STAT and p53, all of which are involved in NHL development. In the  . Expression levels of hsa-miR-20a-5p and hsa-miR-181a-5p in plasma samples of NHL patient groups according to lymphoma grade. Patients with high-grade lymphoma presented higher levels of hsa-miR-20a-5p and lower levels of hsa-miR-181a-5p compared to patients with low-grade lymphoma (p = 0.003 and 0.013, respectively). The graphics represent the −ΔCT of miRNA expression normalized to the endogenous control (mean ± SD).
Next, the plasma expression levels of SNHG16 and SNHG6 were evaluated, and we observed higher plasma levels of both lncRNAs in patients with high-grade lymphomas compared to low-grade lymphomas (SNHG16: fold change = 1.63; SNHG6: fold change = 2.08) (p = 0.029 and p < 0.001, respectively) ( Figure 6). . Expression levels of SNHG16 and SNHG6 in plasma samples of NHL patient groups according to lymphoma grade. SNHG16 and SNHG6 presented higher plasma levels in patients diagnosed with high-grade lymphoma compared to low-grade lymphoma (p = 0.029 and p < 0.001, respectively). The graphics represent the −ΔCT of miRNA expression normalized to the endogenous control (mean ± SD).

miRNAs and LncRNAs Plasma Levels According to IPI and FLIPI Score
We next examined the association of the plasma expression levels of hsa-miR-20a-5p, hsa-miR-181a-5p and lncRNAs SNHG16 and SNHG6 with the IPI scores and FLIPI scores of the NHL patients. The analysis revealed that higher levels of hsa-miR-20a-5p, SNHG16 and SNHG6 were associated with higher IPI and FLIPI scores. On the contrary, we observed a negative association between hsa-miR-181a-5p levels and IPI/FLIPI scores. Specifically, lower levels of hsa-miR-181a-5p were associated with higher IPI and FLIPI scores ( Figure 7). . Expression levels of hsa-miR-20a-5p and hsa-miR-181a-5p in plasma samples of NHL patient groups according to lymphoma grade. Patients with high-grade lymphoma presented higher levels of hsa-miR-20a-5p and lower levels of hsa-miR-181a-5p compared to patients with low-grade lymphoma (p = 0.003 and 0.013, respectively). The graphics represent the −∆CT of miRNA expression normalized to the endogenous control (mean ± SD).
Next, the plasma expression levels of SNHG16 and SNHG6 were evaluated, and we observed higher plasma levels of both lncRNAs in patients with high-grade lymphomas compared to low-grade lymphomas (SNHG16: fold change = 1.63; SNHG6: fold change = 2.08) (p = 0.029 and p < 0.001, respectively) ( Figure 6).
Biomedicines 2022, 10, 1322 10 of 21 Figure 5. Expression levels of hsa-miR-20a-5p and hsa-miR-181a-5p in plasma samples of NHL patient groups according to lymphoma grade. Patients with high-grade lymphoma presented higher levels of hsa-miR-20a-5p and lower levels of hsa-miR-181a-5p compared to patients with low-grade lymphoma (p = 0.003 and 0.013, respectively). The graphics represent the −ΔCT of miRNA expression normalized to the endogenous control (mean ± SD).
Next, the plasma expression levels of SNHG16 and SNHG6 were evaluated, and we observed higher plasma levels of both lncRNAs in patients with high-grade lymphomas compared to low-grade lymphomas (SNHG16: fold change = 1.63; SNHG6: fold change = 2.08) (p = 0.029 and p < 0.001, respectively) ( Figure 6). . Expression levels of SNHG16 and SNHG6 in plasma samples of NHL patient groups according to lymphoma grade. SNHG16 and SNHG6 presented higher plasma levels in patients diagnosed with high-grade lymphoma compared to low-grade lymphoma (p = 0.029 and p < 0.001, respectively). The graphics represent the −ΔCT of miRNA expression normalized to the endogenous control (mean ± SD).

miRNAs and LncRNAs Plasma Levels According to IPI and FLIPI Score
We next examined the association of the plasma expression levels of hsa-miR-20a-5p, hsa-miR-181a-5p and lncRNAs SNHG16 and SNHG6 with the IPI scores and FLIPI scores of the NHL patients. The analysis revealed that higher levels of hsa-miR-20a-5p, SNHG16 and SNHG6 were associated with higher IPI and FLIPI scores. On the contrary, we observed a negative association between hsa-miR-181a-5p levels and IPI/FLIPI scores. Specifically, lower levels of hsa-miR-181a-5p were associated with higher IPI and FLIPI scores ( Figure 7). Figure 6. Expression levels of SNHG16 and SNHG6 in plasma samples of NHL patient groups according to lymphoma grade. SNHG16 and SNHG6 presented higher plasma levels in patients diagnosed with high-grade lymphoma compared to low-grade lymphoma (p = 0.029 and p < 0.001, respectively). The graphics represent the −∆CT of miRNA expression normalized to the endogenous control (mean ± SD).

miRNAs and LncRNAs Plasma Levels According to IPI and FLIPI Score
We next examined the association of the plasma expression levels of hsa-miR-20a-5p, hsa-miR-181a-5p and lncRNAs SNHG16 and SNHG6 with the IPI scores and FLIPI scores of the NHL patients. The analysis revealed that higher levels of hsa-miR-20a-5p, SNHG16 and SNHG6 were associated with higher IPI and FLIPI scores. On the contrary, we observed a negative association between hsa-miR-181a-5p levels and IPI/FLIPI scores. Specifically, lower levels of hsa-miR-181a-5p were associated with higher IPI and FLIPI scores (Figure 7).

Figure 7.
Expression levels of hsa-miR-20a-5p, hsa-miR-181a-5p, SNHG16 and SNHG6 in plasma samples of NHL patient groups according to IPI and FLIPI scores. High levels of hsa-miR-20a-5p, SNH16 and SNHG6 were associated with higher IPI and FLIPI scores, while low levels of hsa-miR-181a-5p were associated with higher IPI and FLIPI scores. The figures represent the −ΔCT of miRNA expression normalized to the endogenous control (mean ± SD).
Given the ceRNAs expression analysis, three groups were established considering the combination of the plasma levels of each ceRNA pair, hsa-miR-20a-5p and SNHG16 plasma levels and hsa-miR-181a-5p and SNHG6 plasma levels, which allowed for the definition of high-, intermediate-and low-risk groups. For the analysis of the hsa-miR-20a-5p:SNHG16 pair, the low-risk group included patients with low levels of hsa-miR-20a-5p and SNHG16; the intermediate-risk group combined both patients with high hsa-miR-20a-5p and low SNHG16 levels and patients with low hsa-miR-20a-5p and high SNHG16 expression; and the high-risk group included patients with high levels of hsa-miR-20a-5p and high levels of SNHG16. Concerning the ceRNA pair hsa-miR-181a-5p:SNHG6, the low-risk group was composed of patients with high expression of hsa-miR-181a-5p and low expression of SNH6. The intermediate-risk group combined both patients with high hsa-miR-181a-5p and high SNHG6 expression and patients with low hsa-miR-181a-5p and low SNHG6 expression. The high-risk group combined patients with a low expression of hsa-miR-181a-5p and high expression of SNHG6 (Table 2). Table 2. Definition of high-, intermediate-and low-risk groups considering the combination of the plasma levels of each ceRNA pair, hsa-miR-20a-5p/SNHG16 and hsa-miR-181a-5p/SNHG6 plasma levels(upwards arrow represents transcript upregulation, downwards arrow represents transcript downregulation).
Given the ceRNAs expression analysis, three groups were established considering the combination of the plasma levels of each ceRNA pair, hsa-miR-20a-5p and SNHG16 plasma levels and hsa-miR-181a-5p and SNHG6 plasma levels, which allowed for the definition of high-, intermediate-and low-risk groups. For the analysis of the hsa-miR-20a-5p:SNHG16 pair, the low-risk group included patients with low levels of hsa-miR-20a-5p and SNHG16; the intermediate-risk group combined both patients with high hsa-miR-20a-5p and low SNHG16 levels and patients with low hsa-miR-20a-5p and high SNHG16 expression; and the high-risk group included patients with high levels of hsa-miR-20a-5p and high levels of SNHG16. Concerning the ceRNA pair hsa-miR-181a-5p:SNHG6, the low-risk group was composed of patients with high expression of hsa-miR-181a-5p and low expression of SNH6. The intermediate-risk group combined both patients with high hsa-miR-181a-5p and high SNHG6 expression and patients with low hsa-miR-181a-5p and low SNHG6 expression. The high-risk group combined patients with a low expression of hsa-miR-181a-5p and high expression of SNHG6 (Table 2). Table 2. Definition of high-, intermediate-and low-risk groups considering the combination of the plasma levels of each ceRNA pair, hsa-miR-20a-5p/SNHG16 and hsa-miR-181a-5p/SNHG6 plasma levels(upwards arrow represents transcript upregulation, downwards arrow represents transcript downregulation).

miRNA and lncRNA Expression Levels in NHL Patients' Plasma Samples according to Bone Marrow Involvement
Considering the risk groups previously defined based on the expression levels of each ceRNA pair, we could observe that high-risk profiles were associated with a higher risk of presenting BMI (Table 3). Patients with a low expression of hsa-miR-181a-5p with high expression levels of SNHG16 had a higher risk of having positive BMI (p = 0.010). Conversely, a higher risk of positive BMI was associated with high expression of hsa-miR-20a-5p and high expression of the SNHG6 profile (p = 0.023). Next, ROC analysis was performed to explore the clinical value of each ceRNA pair levels in the diagnosis of BMI in NHL patients (Figure 8). Results show an AUC of 0.704 for SNHG16/hsa-miR-20a-5p (95% CI 0.579-0.829, p = 0.003) and, finally, an AUC of 0.721 for the ceRNA pair SNHG6/hsa-miR-181a-5p (95% CI 0.579-0.803, p = 0.002).

miRNA and lncRNA Expression Levels in NHL Patients' Plasma Samples according to Bone Marrow Involvement
Considering the risk groups previously defined based on the expression levels of each ceRNA pair, we could observe that high-risk profiles were associated with a higher risk of presenting BMI (Table 3). Patients with a low expression of hsa-miR-181a-5p with high expression levels of SNHG16 had a higher risk of having positive BMI (p = 0.010). Conversely, a higher risk of positive BMI was associated with high expression of hsa-miR-20a-5p and high expression of the SNHG6 profile (p = 0.023). Next, ROC analysis was performed to explore the clinical value of each ceRNA pair levels in the diagnosis of BMI in NHL patients (Figure 8). Results show an AUC of 0.704 for SNHG16/hsa-miR-20a-5p (95% CI 0.579-0.829, p = 0.003) and, finally, an AUC of 0.721 for the ceRNA pair SNHG6/hsa-miR-181a-5p (95% CI 0.579-0.803, p = 0.002).

miRNAs' and LncRNAs' Impact on Overall Survival and Progression-Free Survival of NHL Patients
For the analysis of the OS and PFS, patients were divided in terciles according to the transcript levels using the -ΔCT values of each miRNA and lncRNA (high, intermediate and low levels), in order to analyze the association between the OS and PFS of NHL patients and the expression levels of hsa-miR-20a-5p, hsa-miR-181a-5p, SNHG16 and SNHG6. Kaplan-Meier analysis and log-rank tests revealed that lower plasma levels of hsa-miR-181a-5p were associated with lower OS and PFS (log-rank test: p = 0.017 and p = 0.033; HR: 0.200, p = 0.032; HR: 0.450, p = 0.048, respectively). On the other hand, higher plasma levels of hsa-miR-20a-5p, SHNG16 and SNHG6 were associated with worse OS Figure 8. Plasma levels of ceRNA pair as biomarkers of BM involvement in NHL patients. ROC curve analysis of plasma levels of (a) SNHG16/hsa-miR-20a-5p and (b) SNHG6/hsa-miR-181a-5p ceRNA pair as diagnostic biomarkers differentiating positive BM involvement patients from negative BM involvement patients.

miRNAs' and LncRNAs' Impact on Overall Survival and Progression-Free Survival of NHL Patients
For the analysis of the OS and PFS, patients were divided in terciles according to the transcript levels using the -∆CT values of each miRNA and lncRNA (high, intermediate and low levels), in order to analyze the association between the OS and PFS of NHL patients and the expression levels of hsa-miR-20a-5p, hsa-miR-181a-5p, SNHG16 and SNHG6. Kaplan-Meier analysis and log-rank tests revealed that lower plasma levels of hsa-miR-181a-5p were associated with lower OS and PFS (log-rank test: p = 0.017 and p = 0.033; HR: 0.200, p = 0.032; HR: 0.450, p = 0.048, respectively). On the other hand, higher plasma levels of hsa-miR-20a-5p, SHNG16 and SNHG6 were associated with worse OS and PFS (hsa-miR-20a-5p: p = 0.038 and p = 0.006, HR: 2.834 p = 0.037, HR: 3.898 p = 0.001; SNHG16: p = 0.004 and p = 0.022, HR: 4.481 p = 0.002, HR: 2.346 p = 0.029; SNHG6: p = 0.028 and p = 0.015, HR: 2.621 p = 0.043, HR: 2.325 p = 0.022, respectively) ( Figure 9 and Table 4).  (Figure 9 and Table 4).    When considering the previously defined risk groups, we observed that the significance of the prognostic value was improved when expression levels of the miRNAs were combined with the respective lncRNA pair (Figure 10), showing that patients from the high-risk group from each ceRNA pair presented a significant lower OS and PFS than patients from the intermediate-and low-risk groups. When considering the previously defined risk groups, we observed that the significance of the prognostic value was improved when expression levels of the miRNAs were combined with the respective lncRNA pair (Figure 10), showing that patients from the high-risk group from each ceRNA pair presented a significant lower OS and PFS than patients from the intermediate-and low-risk groups. Multivariate analysis was conducted with significant clinical parameters associated with patients' prognosis and revealed that high levels of both hsa-miR-20a-5p and SNHG16 were independent prognostic factors of poor patient prognosis regarding OS and PFS (OS: p = 0.034 and p = 0.007; PFS: p = 0.002 and p = 0.030) ( Tables 5 and 6). Moreover, low expression of hsa-miR-181a-5p and high expression of SNHG6 were independent prognostic factors associated with a poorer prognosis of NHL patients (OS: p = 0.035 and p = 0.047; PFS: p = 0.029 and p = 0.032) (Tables 5 and 6). In order to compare the predictive ability of death and progression of the different variables, the c index was calculated for each model. Thus, Model 1, which incorporates age at diagnosis, tumor stage, tumor grade, presence of B symptoms, LDH serum levels and ECOG status, presented a predictive capacity with a c index of 0.608 for OS and 0.645 for PFS. Interestingly, when miRNA and corresponding lncRNA plasma levels were added to the previous variables, we observed an increase in the predictive capacity for each model created for both OS and PFS.
Therefore, an increased risk of death and disease progression was found in patients with low hsa-miR-181a-5p levels/high SNHG6 levels and high hsa-miR-20a-5p levels/high SNHG16 levels. Multivariate analysis was conducted with significant clinical parameters associated with patients' prognosis and revealed that high levels of both hsa-miR-20a-5p and SNHG16 were independent prognostic factors of poor patient prognosis regarding OS and PFS (OS: p = 0.034 and p = 0.007; PFS: p = 0.002 and p = 0.030) (Tables 5 and 6). Moreover, low expression of hsa-miR-181a-5p and high expression of SNHG6 were independent prognostic factors associated with a poorer prognosis of NHL patients (OS: p = 0.035 and p = 0.047; PFS: p = 0.029 and p = 0.032) (Tables 5 and 6). In order to compare the predictive ability of death and progression of the different variables, the c index was calculated for each model. Thus, Model 1, which incorporates age at diagnosis, tumor stage, tumor grade, presence of B symptoms, LDH serum levels and ECOG status, presented a predictive capacity with a c index of 0.608 for OS and 0.645 for PFS. Interestingly, when miRNA and corresponding lncRNA plasma levels were added to the previous variables, we observed an increase in the predictive capacity for each model created for both OS and PFS. Therefore, an increased risk of death and disease progression was found in patients with low hsa-miR-181a-5p levels/high SNHG6 levels and high hsa-miR-20a-5p levels/high SNHG16 levels.

Discussion
Despite the remarkable improvements in outcomes, the intrinsic heterogeneity of NHL is reflected in the unpredictability of tumor behavior and, consequently, in the patient's outcome [4][5][6]. Therefore, there is an urgent need to improve current patient prognosis stratification schemes by identifying molecular biomarkers that reflect tumor heterogeneity and clinical behavior. Given the precision medicine era that are currently in, it is now recognized that tissue biopsy, such as BM biopsy and BM aspirate, does not reliably reflect, either temporally and spatially, the whole genomic profile of the tumor. Moreover, they are considered invasive, complex procedures for obtaining samples, and are difficult to reproduce. Circulating tumor-associated components, such as miRNAs and lncRNAs, which can be easily assessed, have potential as lymphoma biomarkers, allowing for a personalized patient follow-up. Among the different lncRNA functions, there is an increasing interest in studying lncRNAs as miRNA sponges, which function as an extra layer of the post-transcriptional regulatory machinery to modulate miRNAmediated gene expression. Moreover, integration of regulatory layers, such as the ceRNA network, represents not only an opportunity to dissect aberrant cellular functions behind the complex process of lymphomagenesis, but also an interesting approach to select more feasible biomarkers with functional relevance. The joint detection of lncRNAs and miRNAs can significantly improve the specificity and sensitivity of liquid-biopsy-based diagnosis and prognosis and increase our understanding of prognostic and predictive phenotypes, eventually leading to better patient follow-up. Therefore, in this study, we aimed to investigate for the first time the prognostic value of the ceRNA network by analyzing the lncRNA-miRNA dynamic pair's expression in plasma samples of NHL patients, not only as predictors of the overall clinical outcomes, but also as biomarkers of important clinical determinants of disease progression, such as BMI. Our results show that hsa-miR-20a-5p, hsa-miR-181a-5p, SNHG16 and SNGH6 were differentially expressed in patients with high-grade lymphomas when compared to low-grade lymphoma patients, and their levels were statistically associated with IPI/FLIPI score, demonstrating their involvement in NHL prognosis. Moreover, the defined risk groups based on ceRNA pair levels could help determine the presence of BMI. Specifically, we determined whether NHL patients with high-risk expression profiles of hsa-miR-20a-5p/SNHG16 and hsa-miR-181a-5p/SNHG6 presented a higher risk of presenting a positive BMI. Moreover, we observed that high levels of hsa-miR-20a-5p, SNHG16 and SNHG6, and low levels of hsa-miR-181a-5p were associated with shorter OS and PFS, which was supported by the multivariate analysis that demonstrated that each transcript is an independent prognosis predictor. These results are further reinforced by the C-index analysis, where the models incorporating the ceRNA network expression analysis clearly improved the predictive capacity compared to the model which only considered the clinicopathological variables.
According to our bioinformatics analysis, lncRNA SNHG16 was identified as one of the hsa-miR-20a-5p sponges. SNHG16 has been widely described as an oncogenic factor in a variety of cancers, including B-cell lymphoma [34]. In fact, Zhu et al. reported that SNHG16 was upregulated in diffuse large B-cell lymphoma tissues and cell lines [34]. Functionally, SNHG16 was shown to induce cell proliferation, cell cycle and invasion and to inhibit apoptosis in the majority of human cancers [35][36][37]. However, the biological function of SNHG16 and its underlying mechanism in NHL are still unknown. SNHG16 transcription has been shown to be regulated by several transcription factors, such as c-Myc, STAT3 and TFAP2A [37,38]. Li et al. demonstrated that c-Myc recruits histone acetyltransferase and induces RNA polymerase II clearance to upregulate SNHG16 transcription, resulting in enhanced cell proliferation, migration and invasion and inhibited cell apoptosis in tumor cells [39]. Furthermore, in another study by Christensen et al., SNHG16 expression was positively modulated by Wnt-regulated transcription factors, including c-Myc [37]. On the other hand, SNHG16 can act as a scaffold on interchromatin clusters regulating gene expression, for example, by interacting with EZH2. The expression of p21 was shown to be directly inhibited by SNHG16 via recruitment of EZH2, which induces cell cycle, cell proliferation and inhibits apoptosis [36]. Moreover, silencing of SNHG16 leads to p21 upregulation and cyclin D1 and cyclin B1 downregulation [40]. SNHG16 can also post-transcriptionally regulate gene expression by acting as a ceRNA sequestering miRNAs. In diffuse large B-cell lymphoma, SNHG16 was shown to sponge miR-497-5p, releasing PIM1 from miR-497-5p-mediated inhibition, which promotes proliferation, cell cycle and inhibits apoptosis of lymphoma cells [34]. Knockdown of SNHG16 in multiple myeloma cells suppressed cell proliferation, induced cell arrest and promoted the apoptosis via inducing cleaved-Caspase-3, cleaved-Caspase-9, Foxa3a and Bax expression, while inhibiting CCND1, Bcl-2, Cyclin D1, PI3K and p-AKT. The SNHG16 effect was shown to be due to sponging miR-342-3p [41]. Li et al. reported that SNHG16 competitively binds to miR-4500, upregulating STAT3 and leading to cell proliferation, migration, invasion and the epithelial-mesenchymal transition process as well as inhibiting cell apoptosis [42]. Zhang et al. reported that SNHG16 acts as a ceRNA, sponging miR-17-5p to upregulate p62, which culminates in the activation of the mTOR/PI3K/AKt pathway and NF-κB signaling to promote proliferation of tumor cells and to repress apoptosis [43]. Recently, it was shown that SNHG16 contains a binding site for hsa-miR-20a-5p, acting as a ceRNA to inhibit its inhibitory function.
MiR-20a-5p is a well-established NHL-associated miRNA belonging to the miR-17-92 cluster, which regulates different stages of B-cell development and central tolerance [22][23][24]. According to our results, similarly to its target lncRNA SNHG16, hsa-miR-20a-5p was also found to be upregulated in the plasma samples of NHL patients. In fact, the miR-17-92 cluster is frequently found to be overexpressed due to genomic amplification (q31.3) in several lymphomas, including DLBCL and FL [44]. Moreover, overexpression of miR-20a is associated with c-Myc expression, whose concomitant expression promoted the onset of tumors and increased their growth in a mouse model of B cell lymphoma [45]. c-Myc is not only involved in the transcription of SNHG16, but can also bind to the promoter region of miR-20a, inducing its expression [46]. Diverse targets of has-miR-20a-5p have been identified that could explain its oncogenic role in NHL. Most notably, miR-20a targets PTEN, whose inhibition results in the activation of the PI3K/AKT pathway, one of the central pathways in NHL development [47][48][49]. Additionally, miR-20a targets CDKN1A/p21, which is a cell cycle inhibitor enforcing cell cycle arrest in G1/S [50]. The oncogene function of miR-20a-5p was also shown to be associated with the inhibition of early growth response (EGR)2, thus promoting cell proliferation and cell cycle progression [51]. Therefore, we propose that the upregulation of both hsa-miR-20a-5p and the corresponding lncRNA pair SNHG16 is due to the preponderant action of the deregulated MYC in lymphomas, which acts as a transcription factor of both genes in order to promote lymphoma cells' proliferation and survival.
Concerning hsa-miR-181a-5p, SNHG6 was found to be one of the lncRNAs acting as a ceRNA. SNHG6 has been found to be significantly overexpressed in different tumors, and is also associated with poor clinical outcomes [52]. In a recent study analyzing lncRNAmediated ceRNA networks in Hodgkin lymphoma, SNHG6 was reported as upregulated, and highly associated with patients' relapse [53]. Therefore, SNHG6 has been characterized as an oncogenic lncRNA involved in the regulation of cell differentiation, proliferation, apoptosis and multidrug resistance [54,55]. Similar to SNHG16, SNHG6 regulates gene expression transcriptionally by recruiting EZH2 to promoter regions of different tumor suppressor genes, such as P27 and P21, and represses their expression through methylation of their promoters [55,56]. Moreover, SNHG6 was identified as a molecular sponge of miR-101, miR-214 and miR-4465, all of which target EZH2 [57][58][59]. Therefore, SNHG6 can modulate the function of EZH2 at multiple levels. By competitively sponging miR-101, SNHG6 also regulates the expression of ZEB1, promoting cell migration and EMT [55]. As determined by our bioinformatics analysis, SNHG6 was identified as targeting miR-181a. In fact, SNHG6 was also shown to act as a molecular decoy for all four members of the miR-181 family. Overexpression of SNHG6 represses miR-181a, which in turn induces JAK2 expression and promotes tumor cell proliferation [60]. Moreover, SNHG6-mediated inhibition of miR-181a was shown to induce proliferation, cell cycle progression and migration and invasion and inhibits apoptosis via upregulation of E2F5 [61]. Using a bioinformatics analysis, we identified an miR-181a-related PPI network linked to proteins such as E2F5, CDKN1B, the MAPK family, BCL-2/BCL2L11 and MET signaling (KRAS and STAT3). This PPI network was enriched in cytokine signaling in the immune system, positive regulation of cell proliferation, cell cycle and deregulation of miRNAs in cancer pathogenesis, as demonstrated in the KEGG and Reactome pathway analysis. Interestingly, when analyzing the Go terms, one of the most enriched was hematopoietic or lymphoid organ development, emphasizing the involvement of miR-181a in the hematologic system and its possible involvement in NHL pathogenesis. In fact, a study by Kozloski et al. reported that miR-181a is a negative regulator of the NF-kB signaling pathway in DLBCL cells, inhibiting tumor cell proliferation and viability [62]. Overall, SNHG6 and miR-181a were shown to directly interact with each, in which SNHG6 functions as an miR-181a decoy; furthermore, there was an overlap of downstream targets, specially involving the MAPK signaling pathway, whose aberrant activation has been previously described in NHL [63]. Therefore, we hypothesize a potential involvement of SNHG6-mediated suppression of miR-181a leading to promotion of proliferation signaling networks and inhibition of apoptosis assisting in NHL progression.

Conclusions
Although several studies have identified miRNAs and lncRNAs in NHL, the number of functional studies on miRNA-lncRNA interplay is still limited, which precludes the defined diagnostic and prognostic importance of the ceRNA network in NHL patients. This approach narrows the scope of research and enhances the prediction accuracy to identify candidate biomarkers with great potential for the diagnosis, prognosis and as therapeutic targets for NHL patients. Moreover, the majority of studies focus on the expression profile in tissue samples and cell lines, with very few studies analyzing their expression levels in the circulation of NHL patients and viewing them as novel less invasive complementary biomarkers in NHL prognosis. Plasma levels of the miR-20a/SNHG16 pair and miR-181a/SNHG6 could serve as new prognostic biomarkers for NHL. In our study, plasma levels of hsa-miR-20a-5/SNHG16 and hsa-miR-181a-5p/SNHG6 were associated with patient clinical outcome, where patients with high hsa-miR-20a-5p/high SNHG16 and low hsa-miR-181a-5p/high SNHG6 presented worse OS and PFS. Therefore, our results indicate that these ceRNA pairs could function as NHL prognostic biomarkers to better identify risk patients and consequently could help improve patients' management. Despite the new approach of our study to identify new circulating NHL biomarkers by introducing the analysis of ceRNA network components, it would be interesting to also analyze the expression of the target mRNAs involved in the identified ceRNA networks and, consequently, analyze their clinical value. Additionally, future studies should focus on validating these results on a larger patient cohort and clarify the biological function of the identified transcripts by performing in vitro studies that permit the modulation of the transcripts' expression and, consequently, investigate their influence on tumor cells' properties.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of Braga Hospital.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Data Availability Statement: No application.