MicroRNA Profiling as a Predictive Indicator for Time to First Treatment in Chronic Lymphocytic Leukemia: Insights from the O-CLL1 Prospective Study.

1. Introduction

Chronic lymphocytic leukemia (CLL) is a B-cell disorder characterized by the monoclonal accumulation of CD5/CD23-positive lymphocytes at multiple sites [1,2,3,4] and a clinical heterogeneity [5,6] which has been correlated with different cytogenetic and molecular features of leukemic cells [7,8,9,10].

Many CLL patients do not require immediate therapeutic intervention, and for the majority, therapy can be initiated after a variable period of disease progression. To avoid useless treatment and to limit treatment only to the progressing cases, a” watch and wait” strategy is generally adopted in which only patients who progress, according to defined clinical criteria, are treated [2]. Due to the aforementioned variability of CLL, prognostic studies aimed at defining the risk of disease progression, warranting therapeutic intervention, have a relevant clinical impact [11,12].

While single prognostic factors, although relevant, may have a limited prognostic significance, the evaluation of multiple biomarkers with varying prognostic power may lead to more definite prognostic models [13]. To this end, prognostic scores have been developed, that incorporate different biomarkers into comprehensive models. The CLL-International Prognostic Index (CLL-IPI) has emerged as one of the well-established models, demonstrating its efficacy in predicting overall survival, as well as time to first treatment (TTFT) and progression-free survival prediction in the context of chemo-immunotherapy [14,15]. Furthermore, our group has documented the predictive value of CLL-IPI for TTFT in early-stage CLL [16]. In the same setting of Binet stage A CLL, Condoluci et al. introduced the International Prognostic Score for Early CLL patients (IPS-E) [17]. Further validated by Smolej et al. [18], the IPS-E and its alternative version (AIPS-E) effectively predicted TTFT and therefore may provide valuable guidance for clinical decision-making in early-stage CLL patients. In our prospective O-CLL1 study, evaluating newly diagnosed Binet stage A patients, both AIPS-E and IPS-E accurately predicted the need for therapy in early-stage CLL patients [19].

The advent of next-generation sequencing (NGS) has facilitated the discovery of an unprecedented number of gene mutations with potential prognostic and predictive value [10]. A recent study found that mutations were detected in roughly 35% of CLL patients, with frequencies ranging from 2.3% to 9.8% of cases, and that NOTCH1 mutations were the most common [20]. When the predictive impact of nine recurrently mutated genes was tested, eight of them (BIRC3, EGR2, NFKBIE, NOTCH1, POT1, SF3B1, TP53, and XPO1) were independently associated with shorter TTFT [20].

In recent years, the role of microRNAs (miRNAs) in CLL has garnered increasing attention [reviewed in [21,22,23]. MiRNAs are evolutionarily conserved, single-stranded non-coding RNA molecules, that can redundantly and simultaneously regulate the expression of multiple genes. Acting post-transcriptionally, miRNAs bind to the messenger RNA of target genes, affecting translation or inducing degradation [24,25].

Several miRNAs have been described to be associated with CLL prognosis [26,27,28,29,30,31,32] and pathogenesis [33,34]. In addition, miRNAs may contribute to the deregulation of apoptosis [33,34,35], BCR signaling [36,37,38], or metabolism [39] in CLL cells, and treatment of CLL patients with chemotherapy or BCR inhibitors were shown to affect the expression of miRNAs involved in these processes [40,41,42]. However, the pathways regulated by most miRNAs remain to be elucidated.

MiRNA expression profiling studies have disclosed correlations between certain miRNA signatures and cytogenetic features and/or IGHV gene mutational status [32,43,44], which are recognized prognostic markers in CLL. Finally, certain miRNA signatures are associated with disease progression and outcome [26,32,36,45,46] or with the onset of Richter’s transformation [47,48,49], a lethal condition characterized by the development of aggressive lymphoma in CLL patients [2,50].

Specifically, miR-15a and miR-16-1, which are located on 13q14 [33], act as tumor suppressors and were the first miRNAs used to predict outcomes in CLL [26]. Consistent with the notion that miR-15 and miR-16 regulate cell apoptosis and proliferation [35,51], our group showed that transfection of miR-15 and miR-16 mimics into del(13)(q14) CLL cells significantly hindered their growth in NOD/Shi-scid,γ(c)(null) (NSG) mice, resulting in substantial tumor regression [52]. Moreover, our previous research revealed that CLL has a miRNA expression profile that closely resembles that of antigen-experienced B cells; some miRNAs in this profile are likely to influence disease progression, as suggested by their correlation with the clinical outcomes [32,53].

In the context of our prospective O-CLL1 study, we herein explored the predictive role of 513 miRNAs on TTFT by: i) incorporating significant miRNAs into a basic prognostic model including IGHV mutational status, del11q and del17p, beta-2-microglobulin (B2M), NOTCH1 mutation, and Rai stage, which are known to predict TTFT in CLL [54,55,56,57] and ii) analyzing their potential association with 50 genes involved in TTFT, that we previously identified by an Artificial Intelligence (AI)-based model. Our approach identified 16 out of 513 miRNAs independently associated with TTFT. Moreover, enrichment analysis was performed to evaluate biological processes regulated by the selected miRNAs, while correlation analysis and multiple bioinformatics software were used to evaluate possible interactions between the 16 prioritized miRNAs and the 50 genes selected as a CLL TTFT predictive signature, defined by an AI model [57].

2. Results

2.1. TTFT Prediction by the Basic Model

The estimated median TTFT for the 224 CLL cases from our prospective O-CLL1 study was 105.2 months (95% CI, 94.2-116.1 months), with 4-year and 8-year cumulative first-treatment free survival of 76.8% (95% CI, 71.1-82.5) and 57.4% (95% CI, 50.1-64.6), respectively (Figure 1). Only 6 out of 9 univariable predictors (i.e., Rai stage, B2M, IGVH and NOTCH1 mutational status, 11q and 17p deletions) remained significantly associated with TTFT when tested in a Cox multivariable model adopting a backward elimination strategy. Univariable and multivariable Cox analyses, focusing on established TTFT predictors are detailed in Table S1A-B. The six variables together gave a prognostic value for TTFT (i.e., a Harrell’C index) of 75% and an explained variation of 45.4%.

2.2. TTFT Prediction by miRNAs

Of the 513 available miRNAs tested, 73 were found to be significantly associated with TTFT in univariable Cox regression analyses (all P<0.05) (Table 1). These 73 miRNAs were then simultaneously included in the same multivariable Cox regression model, and those significantly associated with TTFT were identified by a backward elimination strategy.

Sixteen miRNAs retained an independent association with TTFT (Table 2). For 8 miRNAs (i.e., miR-582-3p, miR-33a-3p, miR-516a-5p, miR-99a-5p, and miR-296-3p, miR-502-5p, miR-625-5p, and miR-29c-3p), lower expression levels were associated with a higher likelihood of treatment (i.e., shorter TTFT). Conversely, for the remaining 8 miRNAs (i.e., miR-150-5p, miR-148a-3p, miR-28-5p, miR-144-5p, miR-671-5p, miR-1-3p, miR-193a-3p, and miR-124-3p), a higher expression was associated with shorter TTFT.

Using these 16 miRNAs, a risk prediction score (ranging from 0 to 100%) was derived (Figure 2).

Inclusion of the miRNA scores into the basic model (i.e., a Cox model including Rai stage, B2M, IGVH, and NOTCH1 mutational status, 11q and 17p deletions) significantly increased both the Harrell’C index (from 75.0% to 81.1%) and the explained variation in TTFT (from 45.4% to 63.3%). Remarkably, the inclusion of the miRNA scores into the basic model also yielded an IDI and an NRI of +14.9% and +44.2%, respectively (Table 3).

2.3. Correlation and Interaction Analysis between miRNAs and Genes Found to be Related to TTFT by an AI Model

Next, we analyzed the relationship between the 16 miRNAs identified as predictive of TTFT and the previously 50 AI model prioritized genes (AI genes) that we found significantly associated with TTFT in the same patients of the O-CLL1 cohort [57]. Fifteen out of 16 miRNAs showed at least one significant negative correlation with the AI genes; in particular, 8 miRNAs (i.e., miR-29c-3p, miR-625-5p, miR-150-5p, miR-144-5p, miR-28-5p, and miR-516a-5p) showed a relationship with one of the top 10 genes selected by the neural network (i.e., CEACAM19, PIGP, FADD, FIBP, IGF1R, COL28A1, QTRT1, MKL1, GNE, SLC39A6) (Figure 3). Only hsa-miR-124-3p showed no correlation with any of the AI genes; in addition, miR-29c-3p had a significant negative correlation with the highest number of genes prioritized by the AI model.

Comparing the previous correlation analysis with the miRComb model, this analysis retained only 5 miRNA-mRNA annotated interactions, displaying a mild to moderate correlation index. miR-29c-3p showed again the highest number of interactions, 3 of which were previously predicted; specifically, miR-29c-3p interaction with PRICKLE1 was annotated in microCosm database, while interactions with ANKRD52 and ZBTB34 were annotated in TargetScan database. miR-29c-3p interaction with KDM5B was also predicted by TargetScan with experimental validation, although with less strong evidence. MiR-625-5p was found to be significantly anticorrelated with IGF1R, an interaction annotated by miRTarBase with a strong validation score. Interestingly, the MiRComb regulation score suggested an association between IGHV mutational status and the interactions between miR-29c-3p with the previously identified AI model selected genes KDM5B, ANKRD52, ZBTB34, PRICKLE1; this dysregulation was not evident when considering the interaction between miR-625-5p and IGF1R (Figure 4). Although the magnitude of the effect of IGHV mutational status on the former interaction is strong, this is largely due to an over-expression of miR-29c-3p in the IGHV mutated group compared to the unmutated one, and to a lesser extent to a slightly not significant downregulation of the previous 4 AI-selected genes in the IGHV mutated group compared to the IGHV unmutated one.

The aforementioned miRNA-mRNA interactions represented by miRComb modelwere found to be specific for CLL. Although the statistical power might be affected by the small sample size, miRNA and mRNA expression profiles of tonsil B cells from normal individuals did not reveal any significant interactions between any of the 16 miRNAs and 50 AI selected genes.

2.4. Pathways Regulation by the 16 miRNAs Linked to TTFT

To further investigate the biological significance of the 16 identified miRNAs, an enrichment analysis was carried out, showing that selected miRNAs are involved in cancer and cellular processes associated with disease progression and drug resistance.

Specifically, Kegg term overrepresentation analysis revealed a significant pathway (i.e., hsa05206) corresponding to the “MicroRNAs in cancer” pathway, associated with the miRNA genes MIR1-2, MIR28, MIR29C, MIR99A, MIR124-3, MIR150, and MIR625 encoding for miR-1-3p, miR-28-5p, miR-29c-3p, miR-99a-5p, miR-124-3p, and miR-150-5p, respectively (Table 4).

Overrepresentation analysis using Wikipath revealed six significant terms associated with the miRNA genes MIR29C, MIR33A, and MIR150 (Table 4). The terms WP1601, WP2249, and WP1545, signatures of “Fluoropyrimidine activity”, “Metastatic brain tumor”, and “miRNAs involved in DNA damage response”, respectively were enriched in the MIR29C gene product; the terms WP299 and WP430, signatures of “Nuclear receptors in lipid metabolism and toxicity” and “Statin inhibition of cholesterol production”, were enriched in the MIR33A gene product, respectively. Finally, the WP2023 term, the signature of the “Cell differentiation expanded index” was enriched in the MIR150 gene product.

Table 4. Summary of significantly enriched ontologies in the 16 analyzed miRNAs displaying ontology-associated miRNAs and the corresponding adjusted p values an d q values.

Enrichment	ID	Term specification	Associated miRNA genes	Adjusted p-value	q value
KEGG ORA summary	hsa05206	MicroRNAs in cancer	MIR1-2, MIR28, MIR29C, MIR99A, MIR124-3, MIR150, MIR625	1.41e-10	7.41e-10
WikiPathways ORA summary	WP299	Nuclear receptors in lipid metabolism and toxicity	MIR33A	0.03	0.009
	WP430	Statin inhibition of cholesterol production	MIR33A	0.03	0.009
	WP1545	miRNAs involved in DNA damage response	MIR29C	0.04	0.01
	WP1601	Fluoropyrimidine activity	MIR29C	0.03	0.009
	WP2023	Cell differentiation expanded index	MIR150	0.04	0.01
	WP2249	Metastatic brain tumor	MIR29C	0.03	0.009

Table 4. Summary of significantly enriched ontologies in the 16 analyzed miRNAs displaying ontology-associated miRNAs and the corresponding adjusted p values an d q values.

Enrichment	ID	Term specification	Associated miRNA genes	Adjusted p-value	q value
KEGG ORA summary	hsa05206	MicroRNAs in cancer	MIR1-2, MIR28, MIR29C, MIR99A, MIR124-3, MIR150, MIR625	1.41e-10	7.41e-10
WikiPathways ORA summary	WP299	Nuclear receptors in lipid metabolism and toxicity	MIR33A	0.03	0.009
	WP430	Statin inhibition of cholesterol production	MIR33A	0.03	0.009
	WP1545	miRNAs involved in DNA damage response	MIR29C	0.04	0.01
	WP1601	Fluoropyrimidine activity	MIR29C	0.03	0.009
	WP2023	Cell differentiation expanded index	MIR150	0.04	0.01
	WP2249	Metastatic brain tumor	MIR29C	0.03	0.009

3. Discussion

Predicting overall and progression-free survival in CLL is a dynamic field of research, that adapts to advances in biology and therapeutics, providing insights to personalized patient planning and guiding clinicians in tailored treatment decisions. Conversely, models predicting TTFT which is not influenced by the treatment choice, play a pivotal role in counseling, family planning, surveillance, and identification of high-risk candidates for potential early intervention in clinical trials [58,59]. Indeed, several predictive models have been published for this purpose [60,61,62,63,64]. Recently, the European Research Initiative on CLL (ERIC) introduced a score that categorizes patients based on IGHV mutation and somatic hypermutation status [20].

Nevertheless, most of these predictive models typically require a combination of clinically accessible factors and molecular factors, which are relatively straightforward to assess. In the O-CLL1 trial, multivariable analyses identified a foundational model comprising Rai stage, B2M, IGHV, and NOTCH1 mutational status, in addition to 17p and 11q deletions [54,55,56,57]. This model also served as the basis for validating the improved prognostic efficacy of lymphocyte doubling time (LDT) [56] and, more recently, for incorporating our previously discovered AI gene expression data into an enhanced model predicting TTFT [57].

The current study focused on analyzing the profile of miRNAs at the time of diagnosis among treatment-naive CLL cases enrolled in the prospective O-CLL1 trial. The primary aim was to identify specific miRNAs that could function as predictive markers to better identify patients in need of early treatment. The establishment of a risk prediction score based on 16 miRNAs identified as independently associated with TTFT, emphasizes the significant impact of miRNA data in predicting TTFT in CLL. Specifically, lower expression of 8 miRNAs (i.e., miR-582-3p, miR-33a-3p, miR-516a-5p, miR-99a-5p, and miR-296-3p, miR-502-5p, miR-625-5p, and miR-29c-3p) was associated with a shorter TTFT, while for the remaining 8 miRNAs (i.e., miR-150-5p, miR-148a-3p, miR-28-5p, miR-144-5p, miR-671-5p, miR-1-3p, miR-193a-3p, and miR-124-3p), the higher expression was associated with a shorter TTFT.

Importantly, the integration of these miRNAs into the basic model significantly improved predictive accuracy, as reflected by improved Harrell’s C index from 75.0% to 81.1% and explained variation in TTFT from 45.4% to 63.3%. Additionally, this integration yielded an IDI of +14.9% and an NRI of +44.2%. This suggests that miRNA expression data may provide meaningful insights, offering a promising opportunity for counseling, and a more reasonable follow-up scheduling through early identification of cases with a higher likelihood of requiring early therapy.

To date, miR-15a and miR-16-1, localized on chromosome 13q14 and functioning as bona fide tumor suppressors [33,35,51] have been utilized as pioneering prognostic indicators in CLL [26]. Subsequent investigations have highlighted the potential of miRNA profiling in enhancing the accuracy of CLL prognostication. Specifically, the expression levels of dysregulated miR-155, miR-181b, miR-29a/b/c, and miR-34a have been systematically correlated with established prognostic biomarkers, including IGHV and TP53 mutational status, as well as ZAP70 expression, thereby exerting a discernible impact on the clinical outcomes of CLL patients [27,30,43,44,46,65,66]. Indeed, our previous study unveiled several dysregulated miRNAs in CLL, indicating their potential role in the pathogenesis of the disease and their contribution to its progression, ultimately influencing the initiation of therapy [32]. In particular, 8 miRNAs (i.e., miR-146b-5p, miR-222-3p, miR-26a-5p, miR-29c-3p, miR-29c-5p, miR-502-3p, and miR-503-5p) out of 15 dysregulated miRNAs, demonstrated a significant role in predicting TTFT with miR-26a-5p and miR-532-3p remaining significantly associated with TTFT after data adjustment for confounders.

Notably, in our current study, only miR-29c-3p was found to be significantly associated with TTFT after applying the Cox multivariable regression model. However, it is essential to consider that, in contrast to our previous analyses and other published studies, our current investigation diverges from the conventional focus based on dysregulated miRNAs. Instead, we opted for a broader statistical strategy, encompassing the evaluation of all the assessable miRNAs to gain a more comprehensive understanding of the miRNA landscape in the context of TTFT prediction. In the present investigation, we have delineated a distinctive signature comprising 16 miRNAs that influence the risk of initiating therapy in CLL. Of note, elevated levels of eight of them correlated with higher treatment likelihood, suggesting their role in disease progression, while the remaining showed an inverse relationship, lowering treatment initiation likelihood, indicating their dual ability to act either as promoters of oncogenic processes or as regulators of tumor suppression.

We reasoned that the understanding of the regulatory networks involving the identified 16 miRNAs could be of interest to further elucidate the CLL pathogenesis; therefore, we adopted a comprehensive in silico approach employing correlation, interaction, and enrichment models to validate and refine miRNA-mRNA analyses in our O-CLL1 dataset. The miRNA-mRNA correlation analysis revealed at least one significant negative correlation between 15 of the identified miRNAs and the set of 50 AI model-based genes we previously identified, with miR-29c-3p being the one with the highest number of correlations as confirmed by multiple tools. MiR-625-5p was found significantly anti-correlated only with IGF1R; interestingly, IGF-1R was functionally validated as a direct target of miR-625-5p also in melanoma cancer cells, being positively modulated by the long non-coding RNA LINC01291 via miR-625-5p sponging [67]. Of particular interest, miR-29c was previously found to exhibit differential expression based on the IGHV status [26,32], which correlated with the expression of some of its specific target genes [66,68]. The present study expanded this observation to four AI-generated genes (ANKRD52, KDM5B, PRICKLE1, and ZBTB34) previously described to be strongly associated with TTFT in CLL [57]. Regulation of KDM5B by miR-29c was described in endometrial carcinoma, showing a correlation between elevated levels of KDM5B with tumor grade and paclitaxel resistance [69]. Finally, the enrichment analysis revealed the multifaceted involvement of miRNAs in cellular activities, from tumor development and resistance to lipid metabolism regulation.

Our previous investigation demonstrated that IL-23R expression on CLL cells independently estimated TTFT in the O-CLL1 cohort, supporting the notion of an essential IL-23 autocrine loop driving CLL expansion [70]. Furthermore, we found that miR-146b-5p modulates IL-12Rβ1 expression, which influences TTFT, with lower levels associated with shorter TTFT duration [53]. Notably, in the current study, miR-148a-3p emerged as an independent prognostic factor in the multivariable model, showcasing a striking complementarity of 18 out of 21 bases with a sequence within the interleukin 23 subunit alpha (IL-23A) mRNA (miRWalk: refseq ID NM_016584, http://mirwalk.umm.uni-heidelberg.de/human/gene/51561/), thus expanding the exploration of IL23A gene regulation.

In conclusion, the identification of specific miRNAs as predictors of TTFT in CLL represents a promising avenue for refining risk stratification and predicting therapeutic needs. The integration of miRNA data into predictive models holds the potential to improve the accuracy of clinical decision-making in CLL management. However, further validation studies and comprehensive functional analyses are needed to confirm the reliability and robustness of these findings.

4. Materials and Methods

4.1. Patient Population and Study Design

In the observational O-CLL1 study (clinicaltrials.gov identifier NCT00917540), 224 newly diagnosed Binet A CLL cases from 40 Italian institutions were prospectively enrolled for miRNAs analysis [32,71]. The data are deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/) and are accessible via GEO Series accession number GSE40533.

All participants gave written informed consent, and the study was approved by the appropriate institutional review boards. Detailed inclusion and exclusion criteria have been previously outlined [57]. Specifically, the recruitment was limited to cases diagnosed within 12 months, age <70 years, and at Binet stage A [32,71].

4.2. Assessment of Biological Markers

The diagnosis was confirmed by flow cytometric analysis, which assessed the proportion of CD5/CD19/CD23 triple-positive B cells; monoclonal antibodies (mAbs) against CD19-FITC (BD Biosciences Pharmigen), CD23-PE (BD Biosciences, BD), and CD5-PC5 (Beckman Coulter Immunotech) were used for this purpose. CD38-positive leukemic cells were quantified through triple staining with CD19 FITC (BD Biosciences), CD38 PE (BD Biosciences), and CD5 PC5 (Beckman Coulter Immunotech) mABs. CD38 positivity was defined with a cutoff of ≥20%, following previous reports [71].

ZAP-70 was detected by flow cytometry using a ZAP-70-FITC (clone 2F3.2, Millipore) or an isotype control mAb (mouse IgG2a-FITC, BD Biosciences) as previously described [71]. Briefly, peripheral mononuclear cells purified from fresh heparinized CLL samples by Ficoll-Hypaque gradient were first incubated with CD3 PE-CY7, CD19 PE, and CD5 PC5 mAbs (Becton Dickinson & Co.), fixed, permeabilized with fix and perm reagents (Caltag Laboratories), and exposed to ZAP-70 or the isotype control mAb. A cutoff value of more than 30% for ZAP-70 positivity was used as previously reported [71] and calculated by receiver-operating characteristic analysis as the most suitable ZAP-70 cutoff value to discriminate IGHV-Unmutated (UM) from IGHV-Mutated (M) cases. All flow cytometric analyses were performed on a FACSCalibur flow cytometer (BD Biosciences).

Cytogenetic abnormalities including chromosome deletions 11q23 and 17p13, were examined by fluorescence in situ hybridization (FISH) in a purified CD19-positive population, according to established protocols [72]. The mutational status of the IGHV gene was assessed on cDNA samples [73,74,75] by aligning sequences to the IMGT directory and by analyzing them using IMGT/VQUEST software. The NOTCH1 mutation hotspot was determined by next-generation deep sequencing, as previously outlined [54].

4.3. miRNAs Analysis

Total RNA was extracted from CD19-positive purified B-cell samples using TRizol reagent (Life Technologies). Subsequently, RNA quality was evaluated through the Agilent 2100 Bioanalyzer (Agilent Technologies). Any RNA sample with poor quality, as indicated by an RNA integrity number < 7, was excluded from microarray analyses. After sample collection, processing was performed according to the guidelines of the Agilent manual for the Human miRNA Microarray V2 platform (Agilent Technologies), which includes miRNAs from the Sanger miRBase (v10.1). Expression values of miRNAs were computed using Agilent Feature Extraction Software 10.1, followed by summarization and background subtraction. MiRNAs with low expression (all detection calls missing), non-human miRNAs, and miRNAs expired according to Sanger miRBase release 15 (April 2010) were excluded. The obtained raw data underwent quantile normalization, conversion to positive values with a minimum value of 1, and log2 transformation using the R-2.14 statistical environment (http://www.r-project.org/) [32]. Array old identifiers were converted to the last version with Mirbase (i.e., miR-33a*=miR-33a-3p, miR- miR-99a= 99a-5p, miR-625= miR-625-5p, miR-29c=miR-29c-3p, miR-150=miR-150-5p, miR-148a= miR-148a-3p, miR-144*=miR-144-5p, miR-1=miR-1-3p, and miR-124=miR-124-3p).

4.4. miRNA-mRNA Correlation, Interaction, and Enrichment Analyses

To elucidate potential relationships between the miRNAs significantly associated with TTFT and the previously identified 50 genes implicated in TTFT as detected by an AI-based model [57], a Spearman correlation analysis between miRNA and mRNA expression profiles was performed using the O-CLL1 database [71]. To retain miRNA-mRNA correlations implicating a miRNA-mRNA interaction, an additional Spearman correlation analysis between expression profiles was performed using the miRComb model [76], selecting significant correlations that were suggestive of an interaction annotated in at least one of the following prediction databases, namely microCosm (version 5.18) and TargetScan 6.2. The sensitivity of the previous model was improved by extending the annotation query to the miRTarBase database [77], which annotated experimentally proven miRNA-mRNA interactions accompanied by a qualitative validation score including the two values “Low evidence” and “Strong evidence” based on the experimental techniques involved in the validation protocol. At the same time, the specificity of the model was increased by selecting correlations with a corresponding P-value, adjusted for multiple comparisons inherent to the whole gene expression and miRNA profile analysis, below a threshold of 0.05. Additionally, the mirComb model was applied to the expression dataset to calculate a miRNA-mRNA regulation score, which allowed the assessment of the effect of IGHV mutational status on each miRNA-mRNA interaction retained; the regulation score was calculated as follows:

$$score=-2\left({\log }_{2}rati{o}_{mRNA}\mathrm{ }·\mathrm{ }{\log }_{2}rati{o}_{miRNA}\right)$$

where log₂ ratio_mRNA and log₂ ratio_miRNA equal to the log-transformed fold changes, of mRNA and miRNA expressions of a given interaction pair respectively, calculated by comparing the IGHV mutated group with the unmutated group. Specifically, the previous score measured the magnitude of interaction deregulation between a miRNA and an mRNA as a result of the effect mediated by a categorical variable: a positive score implied interaction deregulation, while a negative score indicated a preserved interaction. To ensure the specificity of the miRNA-mRNA interactions detected from the O-CLL1 cohort for CLL pathogenesis, a prior correlation analysis was also carried out on a control set of normal B-cell sample subpopulations collected from tonsils of healthy individuals. In particular, miRNA and mRNA expression profiles were collected for 3 samples of Naive B cells, 2 samples of Marginal Zone-like B cells, and 3 samples of memory cells (GSE51529)[43,78].

To elucidate the biological significance of the identified miRNAs that are correlated with clinical outcomes, we conducted an over-representation analysis (ORA) of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Wikipathways terms. This analysis was performed using the R package ClusterProfiler, a versatile tool designed for the interpretation of omics data [79]. ORA allows the identification of enriched biological pathways or processes associated with a given gene set, thereby providing insight into the potential functional roles of the miRNAs under investigation.

4.5. Statistical Analyses

TTFT was calculated during the watch-and-wait period, which lasted from the date of diagnosis to the start of therapy or the last follow-up. The prognostic impact of standard risk factors and miRNAs for TTFT was preliminarily investigated by univariable Cox regression analyses, with data presented as hazard ratios (HR) with 95% confidence intervals (CIs). To obtain parsimonious models, all univariate predictors of TTFT (i.e., variables with P<0.05) were tested in Cox analyses with a backward elimination strategy to identify relevant prognostic variables among standard risk factors or miRNAs. Harrell’s c-index, the explained variation in TTFT (an index combining calibration and discrimination), the integrated discrimination improvement (IDI), and the net reclassification index (NRI) were used to assess the accuracy of prognostic models and to measure the gain in prognostic accuracy attributable to relevant miRNAs. To generate a risk prediction rule based on miRNAs, a logistic regression model was fitted with TTFT as the dependent variable and the miRNAs which remained significantly associated with the study outcome by the backward elimination strategy as independent variables. This analysis resulted in a miRNA score (ranging from 0 to 100%), which was then used for further analyses. Statistical calculations were performed using SPSS for Windows v.21 (IBM, Chicago, Illinois, USA) and Stata 16 (StataCorp, Texas, USA).