Increased Plasma Levels of lncRNAs LINC01268, GAS5 and MALAT1 Correlate with Negative Prognostic Factors in Myelofibrosis

Simple Summary Myelofibrosis (MF) displays the worst prognosis among Philadelphia-negative chronic myeloproliferative neoplasms. There is no curative therapy for MF, except for bone marrow transplantation, which however has a consistent percentage of failure. There is thus an urgent need of novel biomarkers to complement current stratification models and to enable better management of patients. To address this issue, we herein measured the plasma levels of several long noncoding RNAs (lncRNAs). Circulating lncRNAs has been already largely described as potential non-invasive biomarkers in cancers. In our study we unveiled that LINC01268, MALAT1 (both p < 0.0001) and GAS5 (p = 0.0003) plasma levels are significantly higher in MF patients if compared with healthy donors, and their increased plasma levels correlate with several detrimental features in MF. Among them, LINC01268 is an independent variable for both OS (p = 0.0297) and LFS (p = 0.0479), thus representing a putative new biomarker suitable for integrate contemporary prognostic models. Abstract Long non-coding RNAs (lncRNAs) have been recently described as key mediators in the development of hematological malignancies. In the last years, circulating lncRNAs have been proposed as a new class of non-invasive biomarkers for cancer diagnosis and prognosis and to predict treatment response. The present study is aimed to investigate the potential of circulating lncRNAs as non-invasive prognostic biomarkers in myelofibrosis (MF), the most severe among Philadelphia-negative myeloproliferative neoplasms. We detected increased levels of seven circulating lncRNAs in plasma samples of MF patients (n = 143), compared to healthy controls (n = 65). Among these, high levels of LINC01268, MALAT1 or GAS5 correlate with detrimental clinical variables, such as high count of leukocytes and CD34+ cells, severe grade of bone marrow fibrosis and presence of splenomegaly. Strikingly, high plasma levels of LINC01268 (p = 0.0018), GAS5 (p = 0.0008) or MALAT1 (p = 0.0348) are also associated with a poor overall-survival while high levels of LINC01268 correlate with a shorter leukemia-free-survival. Finally, multivariate analysis demonstrated that the plasma level of LINC01268 is an independent prognostic variable, suggesting that, if confirmed in future in an independent patients’ cohort, it could be used for further studies to design an updated classification model for MF patients.


Introduction
The Philadelphia (Ph)-negative myeloproliferative neoplasms (MPNs) are clonal stemcell disorders characterized by an increased proliferation and abnormal differentiation of myeloid cells. They include polycythemia vera (PV), essential thrombocythemia (ET) and primary myelofibrosis (PMF), the latter identified as prefibrotic (pre-PMF) and overtly fibrotic (overt-PMF) stages according to the degree of bone marrow fibrosis [1]. In addition to PMF, which has a de novo onset, secondary myelofibrosis (SMF) might evolve from PV or ET, resulting in post-PV myelofibrosis (PPV-MF) or post-ET myelofibrosis (PET-PV), respectively, all generally referred to as myelofibrosis (MF) [2].
MPNs originate from somatic mutations occurring in hematopoietic stem cells. Mutations in JAK2, MPL and CALR genes have been identified as driver events of disease onset, although observed with different frequencies in the specific diseases [3]. In addition, "High Molecular Risk (HMR)" mutations (e.g., in ASXL1, IDH1/2, SRSF2 and EZH2 genes) have been associated to a worse prognosis and a more frequent leukemic transformation [4]. The serial acquisition of somatic mutations underlies the clonal evolution in MPNs and we have recently reconstructed by single cell analysis the sequence of mutational events associated to PMF progression [5]. MF is the most severe among the MPNs and is characterized by the occurrence of bone marrow fibrosis and a consequent abnormal increase in the number of circulating CD34+ hematopoietic progenitors [6]. The most frequent clinical manifestations of MF include splenomegaly, caused by extramedullary hematopoiesis, bleedings and constitutional symptoms. Besides infections, thrombosis and cardiovascular complications, the main cause of death is due to transformation to acute myeloid leukemia (AML), which occurs in 15-20% of cases and is not responsive to conventional therapeutic treatments [4].
Prediction of patients' survival, with the aim to tailor the best therapeutic approach, is based on current prognostic models, such as the International Prognostic Scoring System (IPSS) [7], dynamic-IPSS (DIPSS) [8] and the more recently developed Mutation-enhanced International Prognostic Score System (MIPSS70) [9] and MIPSS70+ version 2.0 [10,11]. Although allogenic stem cell transplantation (ASCT) is currently the only curative therapeutic approach for MF, it is still associated to frequent graft-related morbidities and deaths [12]. Thus, it is of great interest to establish new specific prognostic factors that may predict disease outcome to identify high risk patients eligible for ASCT. In this regard, we have recently proposed that gene expression profiles in MF granulocytes might integrate current prognostic models [13].
Long non-coding RNAs (lncRNAs) are RNA molecules longer than 200 nucleotides and lacking protein-coding potential. It has been demonstrated that lncRNAs are heavily involved in several cellular mechanisms, inter alia, the regulation of gene expression, either in cis or in trans, chromatin remodeling, RNA processing and organization of nuclear domains. They have been found embedded as key components in highly wired networks of gene expression regulation [14]. LncRNAs expression is finely regulated, displaying even higher tissue and cell specificity than protein-coding mRNAs and making them candidate biomarkers and potential targets for therapeutic approaches [15,16].
An exhaustive description of the roles of lncRNAs in normal and malignant hematopoiesis has yet to be provided. However, it has been demonstrated that they play a pivotal role in the differentiation processes within specific hematopoietic lineages (e.g., EGO in eosinophils commitment [17], HOTAIRM1 in RA-induced granulocytic differentiation [18]). In many cases, a dysregulation in a single lncRNAs expression has also been linked to blood malignancies development and progression [19,20].
Over the past decade, evidence has been accumulating about the detection of RNAs in biofluids and their putative application as biomarkers in the early diagnosis and prediction of prognosis and therapy response for cancer patients [21]. Circulating RNAs are secreted from cells through active processes or released after cell death [22]. Notably, lncRNAs might be extensively detected in body fluids, despite of their low levels of expression, since they form secondary structures that confer them high stability [23].
In the present study we took to explore the expression profile of circulating lncRNAs in plasma samples from MF patients, with the aim to assess their potential significance as prognostic indicators of disease progression in the light of performing personalized therapeutic approaches. We identified by quantitative reverse transcription PCR (qRT-PCR) increased levels of circulating LINC01268, MALAT1 and GAS5 in MF patients and we correlated their expression with the patients' clinical features and prognosis.

Patients
Human CD34+ cells were purified by using magnetic beads separation system (CD34 MicroBead Kit UltraPure, Miltenyi Biotech, Auburn, CA, USA) from peripheral blood (PB) of 83 patients (MFs). As control samples, CD34+ cells were isolated from 14 PB samples and 12 bone marrow (BM) samples from healthy donors (HDs).
Plasma samples were separated from 65 HDs and 143 patients with a diagnosis of PMF or SMF, recruited from seven Italian research centers. PMF was diagnosed according to 2016 World Health Organization criteria [1], whereas International Working Group for Myeloproliferative neoplasms Research and Treatment criteria were exploited for the diagnosis of PET-MF and PPV-MF [2]. This study was conducted in accordance with the Declaration of Helsinki and after the approval by local ethics committees. All subjects provided informed written consent.

Plasma Isolation
Blood samples were collected via venipuncture in EDTA or sodium citrate containing tubes. Plasma was separated from blood cells by centrifugation at 1900× g at room temperature for 5 min. Cell debris were removed with a second centrifugation at 17,000× g at +4 • C for 5 min. Hemolyzed samples were identified by spectrophotometric analyses, measuring the absorbance of hemoglobin at 414 nm using the NanoDrop ND-1000 spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA). Samples with a hemoglobin absorbance ≥ 0.3 were considered hemolyzed and excluded from further analyses [24]. Plasma was then aliquoted and stored at −80 • C.

RNA Purification
Total RNA from CD34+ cells was extracted using the miRNeasy Micro Kit (Qiagen, Hilden, Germany) following manufacturer's instructions.
For plasma analysis, total RNA purification was obtained from 100 µL of plasma samples by using the MagMAX TM mirVana TM Total RNA Isolation Kit (Applied Biosystems, ThermoFisher Scientific, Waltham, MA, USA) according to the manufacturer's specifications for plasma samples. Samples were thawed at room temperature to avoid the formation of cryoprecipitates. Synthetic spike-in RNA (#A39179, Taqman Universal RNA spike In reverse Transcription control, Applied Biosystems, Waltham, MA, USA) was introduced in each sample during RNA extraction as an exogenous control. Total RNA was eluted in a final volume of 30 µL and then stored at −80 • C.
The recovery of total RNA obtained with this protocol was too low to evaluate RNA quantity and integrity by spectrophometric or fluorimetric measurements. For this reason, the following steps were performed starting from the same volume of RNA elution for each sample.

Measuring lncRNAs Expression in CD34+ Cells
Reverse transcription reaction was performed starting from 100 ng of total RNA from CD34+ cells by using Superscript TM VILO TM cDNA Synthesis Kit (Applied Biosystems, ThermoFisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions.
Next, a preamplification step was run mixing 5 µL of cDNA with pooled Taqman assays and TaqMan TM PreAmp Master mix (Applied Biosystems, ThermoFisher Scientific, Waltham, MA, USA) to increase uniformly the amount of cDNA for each target lncRNAs. Pooled Taqman assays mix was prepared by combining equal volumes of assays for all the targets of interest (listed in Table S1) so that each assay is at a final concentration of 0.2×.
Amplification reactions were performed on Applied Biosystems TM QuantStudio TM 12K Flex Real-Time PCR System (ThermoFisher Scientific, Waltham, MA, USA). We exploited custom OpenArray plates (56-assays format) both with custom and predesigned Taqman assays. Each OpenArray plate allowed simultaneous detection of lncRNAs of 16 samples in triplicate. For each sample, pre-amplified cDNA was diluted 1:20, mixed with TaqMan TM OpenArray™ Real-Time PCR Master Mix and loaded onto OpenArray™ plate by means of QuantStudio™ OpenArray ® AccuFill™ System. Data analysis was performed by means of Relative Quantification App on Ther-moFisher Cloud (ThermoFisher Scientific, Waltham, MA, USA) applying the C RT method ("relative threshold" method). The C RT method, made necessary by the small volume loaded on each through-holes, differs from the common C T method since it defines a threshold for each reaction based on the efficiency of reaction. We included, for further analyses, all the reactions with a C T ≤ 30 and good amplification scores (Amp Score > 1.0, Cq confidence > 0.8). Expression data have been analyzed applying the ∆∆C T method [25], using the geometric mean of four endogenous controls (RPL19, IPO8, HPRT1 mRNAs and LINC01353 non-coding RNA) as normalization factor and the median ∆Ct score of the value of HDs as calibrator. Stability expression of candidate reference controls among all the analyzed samples was evaluated by means of GeNorm algorithm [26]. Fold change (FC) expression for each lncRNA was calculated by using the 2 −∆∆CT method comparing MF plasma samples with HDs.

Measuring lncRNAs Levels in Plasma
Sample workflow to detect lncRNAs in plasma followed the same steps above described, with the following modifications.
Reverse transcription reaction was performed starting from 10 µL of total RNA prepared from plasma samples. Amplification reactions were performed in duplicate for each sample, using the TaqMan TM Fast Advanced Master Mix and predesigned TaqMan probes specific for each target (listed in Table S2)   The evaluation of the previously added synthetic spike-in RNA was run simultaneously to assess the efficiency of purification and reverse transcription processes.
For data analysis we included all the reactions with a C T < 33, and maximum cycle value (C T = 40) was attributed to amplification reactions with C T ≥ 33. Expression data have been analyzed applying the ∆∆C T method [25], using H19 as endogenous control since resulted as the lncRNA more stably expressed between control and MF groups. In addition, it displayed a distribution in the two groups similar to miR-23a, a small RNA previously reported as reference for circulating RNAs ( Figure S1) [24]. On the contrary, GAPDH, HPRT1 and ACTB resulted significantly more enriched in MF plasma samples and, also being coding mRNAs, are thus excluded as endogenous controls.
The median ∆Ct value of HDs was exploited as calibrator. FC expression for each lncRNA was calculated by using the 2 −∆∆CT method comparing MF plasma samples with HDs.

Statistical Analysis
Comparisons between two groups of numerical variables were performed by using the non-parametric Mann-Whitney U Test. Fisher's exact test or the Chi-square (χ 2 ) test was performed to discriminate differences between categorical variables. Outliers were removed from further analyses by the ROUT method. Overall survival (OS) and leukemia-free survival (LFS) were calculated from the date of sample collection to the date of last follow-up or event occurrence (death or leukemic transformation, respectively). OS and LFS analyses were performed with the Kaplan-Meier method, and the log-rank test was used to compare two curves. Multivariate analyses were carried out by means of Cox proportional hazard regression for OS and LFS using R version 3.4.1. p-values (p) < 0.05 were considered as statistically significant. All graphs and statistical analyses were performed using GraphPad Prism version 8.4.0 for Mac (GraphPad Software, San Diego, CA, USA). PCA and hierarchical clustering were performed using Partek Genomics Suite (GS) software, version 7.0 (Partek Inc., St. Louis, MO, USA).

LncRNAs Expression Profile in CD34+ Cells from MF Patients
Eighty-three MF patients (n = 41 PMF, n = 26 PET-MF, n = 16 PPV-MF) and twenty-six HDs were analyzed for the expression of a list of 38 lncRNAs ( Figure 1 and Table S1) by OpenArray real-time PCR platform. Some of the targets have been addressed by using multiple assays to obtain the higher coverage of all their different isoforms. In order to display relationship between samples, we performed unsupervised Principal Component Analysis (PCA) of expression data that demonstrated that MF samples clustered together and clearly differed from healthy controls (Figure 1a). By ANOVA analysis, we identified 26 differentially expressed lncRNAs with a FC > |1.5| in MF samples compared to HD controls. According to hierarchical clustering analysis (FDR < 0.05), these differentially expressed lncRNAs clearly separated our data set into two main branches, one of them including control samples and the other containing all MF patients, respectively. Samples from PMF and SMF patients resulted interspersed and did not clearly separate (Figure 1b). Together, these results depicted the expression profile of lncRNAs in CD34+ cells and showed that some lncRNAs (such as MEG3, LINC01268, TCL6, LINC01296) are differentially expressed in samples isolated from MF if compared with their normal counterparts.
The plasma levels of the remaining nine circulating lncRNAs are shown in Figure 2, whereas the relative FC is reported in Table S2. Several circulating RNAs were assessed as putative reference endogenous controls. Among them, H19 was selected since recovered at similar levels in control and MF plasma samples. Moreover, its distribution mirrors the pattern displayed by RNAs (i.e., hsa-miR-23a) previously reported as endogenous control for circulating RNAs [24] Table S2. All targets detected in fewer than half of plasma samples of MF patients (TCL6, HOXB-AS3, MEG3 and MIR-155HG) were excluded from further analyses since considered not indicative (data not shown). The plasma levels of the remaining nine circulating lncRNAs are shown in Figure 2, whereas the relative FC is reported in Table S2. Several circulating RNAs were assessed as putative reference endogenous controls. Among them, H19 was selected since recovered at similar levels in control and MF plasma samples. Moreover, its distribution mirrors the pattern displayed by RNAs (i.e., hsa-miR-23a) previously reported as endogenous control for circulating RNAs [24]   Therefore, seven lncRNAs were selected as potential biomarkers in MF patients and tested for correlations with patients' clinical and molecular features. Therefore, seven lncRNAs were selected as potential biomarkers in MF patients and tested for correlations with patients' clinical and molecular features. We investigated the potential association between features considered detrimental for MF progression and the seven lncRNAs whose plasma levels has been found increased in MF samples.
In order to unveil potential association between a number of detrimental clinical features and high plasma levels of circulating lncRNAs, the patient cohort was split into two groups (low-or high-) according to the lncRNAs levels. For LINC01268, LINC00899 and CDKN2B-AS1, displaying a bimodal distribution within MF samples (Figures 2 and S2), a cutoff value was arbitrarily set between the two peaks. For the remaining target RNAs (MALAT1, GAS5, TUG1, NEAT1), displaying a normal distribution, the median value among the MF samples was used as cutoff according to the "median split" method ( Figures 2 and S2). A good correlation with detrimental clinical variables was observed for LINC01268, MALAT1, GAS5, LINC00899 and TUG1 (Table 1), whereas NEAT1 and CDKN2B-AS1 levels seems to be less related to MF prognosis (Table S3).
An increased plasma level of circulating MALAT1 or GAS5 was associated with a lower median hemoglobin level (p = 0.0139 and p = 0.0072, respectively) and a decrease in the median hematocrit value is displayed only by the group with high levels of MALAT1. Platelets median count was lower in patients with high levels of circulating LINC01268 (p = 0.0021) and GAS5 (p < 0.0001). The presence of constitutional symptoms was significantly associated with higher plasma levels of MALAT1 (p = 0.0302) or GAS5 (p = 0.0306), whereas splenomegaly was correlated to higher LINC01268 (p = 0.0012), MALAT1 (p < 0.0001) or GAS5 levels (p = 0.0006) ( Table 1). MALAT1 and GAS5 have been found associated to MF clinical subtypes (p = 0.0001 and p = 0.0204, respectively), with pre-PMF samples enriched within the group of patients with low levels of circulating MALAT1; on the contrary, high levels of GAS5 correlated with patient with overt-PMF (Table 1). Accordingly, MALAT1 (p = 0.0004) and GAS5 (p = 0.0023) levels positively correlate with a more severe grade of bone marrow fibrosis. Furthermore, high levels of LINC01268 (p = 0.0144) are also associated to a more marked grade of fibrosis ( Figure 4a,d,g). Higher plasma levels of LINC01268, MALAT1, GAS5, TUG1 and LINC00899 lncRNAs were associated with elevated white blood cells (WBCs), circulating CD34+ cells together with a marked lactate dehydrogenase (LDH) plasmatic activity (Figure 3a-o).  MALAT1 and GAS5 have been found associated to MF clinical subtypes (p = 0.0001 and p = 0.0204, respectively), with pre-PMF samples enriched within the group of patients with low levels of circulating MALAT1; on the contrary, high levels of GAS5 correlated with patient with overt-PMF (Table 1). Accordingly, MALAT1 (p = 0.0004) and GAS5 (p = 0.0023) levels positively correlate with a more severe grade of bone marrow fibrosis. Furthermore, high levels of LINC01268 (p = 0.0144) are also associated to a more marked grade of fibrosis (Figure 4a,d,g). . Correlation analysis of level of circulating lncRNAs with fibrosis and mutational status of MF patients. Histograms obtained from contingency tables computed to correlate lncRNAs plasma levels and prognostic estimates such as fibrosis grade (a,d,g,j,m), assignment to DIPSS categories (b,e,h,k,n) and HMR mutations (c,f,i,l,o). Samples with low or high levels of target RNA are represented in blue and red, respectively. p = p-value was computed by Fisher's exact test or the Chisquare (χ 2 ) test. Notably, an interesting correlation emerged between circulating lncRNAs and DIPSS prognostic score system. Patients belonging to DIPSS Intermediate-2 and High categories being enriched in high LINC01268 (p = 0.0007), MALAT1 (p = 0.0008) and GAS5 (p = 0.0001) groups, whereas the majority of patients classified as DIPSS Low are gathered in groups with low plasma levels of lncRNAs (Figure 4b,e,h).
Taken together, these data highlighted the association between clinical features commonly monitored in MF patients and circulating levels of several lncRNAs. In particular, we unveiled that LINC01268, MALAT1 and GAS5 are the lncRNAs more frequently associated to clinical detrimental features in MF patients and thus of great interest as putative markers of MF prognosis.

High Plasma Levels of LINC01268, MALAT1, GAS5, TUG1 and NEAT1 Are Associated with MF Patients' Molecular Features
We selected a cohort of patients displaying expected frequencies in driver mutations distribution: JAKV617F was the most common mutation, occurring in 61% of samples, MPL in 6%, CALR in 27%, whereas 6% of patients are triple negative. Mutation in the CALR gene was the sole driver mutation differentially distributed, being enriched in patients with high plasma levels of LINC01268 (p = 0.0184), TUG1 (p = 0.0311) and NEAT1 (p = 0.0059).
On the other hand, regarding HMR mutations, the frequency of gene variants in ≥1 HMR gene was increased in patients with high levels of LINC01268 (p = 0.0087) and MALAT1 (p = 0.0163) (Figure 4c,f). Furthermore, occurrence of ≥2 HMR mutations seemed to be associated with LINC01268 levels, being detected only in patients displaying high levels of circulating LINC01268.
In particular, ASXL1 gene variants were more frequent in patients with high levels of LINC01268 (p = 0.0025), MALAT1 (p = 0.0101) as well as GAS5 (p = 0.0183), while SRSF2 variants were correlated with low levels of NEAT1.

High Plasma Levels of LINC01268, GAS5 and MALAT1 Affect OS
Since patients with high levels of some circulating lncRNAs displayed several detrimental features, we evaluated the association between plasma levels of these seven lncR-NAs and OS. The Kaplan-Meier estimates demonstrated that patients with high LINC01268 had an inferior survival rate compared with patients with low LINC01268 (p = 0.0018, log-rank test; hazard ratio (HR) = 2.705) (Figure 5a). In particular, OS at six years of follow up was 32.1% and 64% in the group of high or low LINC01268, respectively (data not shown). Similarly, patients with high GAS5 plasma levels displayed an inferior survival compared to patients with low levels (p = 0.0008, log-rank test; HR = 2.704) (Figure 5c), with survival proportion of 25.4% and 62%, respectively (data not shown). A statistically significant difference in OS was also observed between groups dichotomized according to circulating levels of MALAT1 (p = 0.0348, log-rank test; HR = 1.803) (Figure 5e). No differences were instead observed in OS of MF patients according to plasmatic levels of TUG1, LINC00899, NEAT1 and CDKN2B-AS1 ( Figure S3a,c,e,g). In agreement, correlation analysis indeed revealed that, within groups displaying high levels of LINC01268 (p = 0.0035), GAS5 (p = 0.0115) and MALAT1 (p = 0.0436), deaths resulted more frequent. (Tables 1 and S3).
A similar approach was exploited to assess association between circulating lncRNA levels and LFS. Kaplan-Meier curves showed that high levels of LINC01268 associated with a shorter LFS compared to patients with lower levels (p = 0.0063, log-rank test; HR = 9.685) (Figure 5b). No correlation with LFS was instead observed for MALAT1 and GAS5 levels (p = 0.5716 and p = 0.1534, respectively, log-rank test) (Figure 5d,f) or for the remaining lncRNAs studied ( Figure S3b,d,f,h).
Together these results demonstrated that circulating lncRNAs LINC01268, GAS5 and MALAT1, all detected at higher levels in MF patients if compared to HDs and associated to several features detrimental for MF prognosis, correlated with a shorter OS. Patients' cohort was stratified into two groups (low and high) according to the plasma levels of target lncRNA, as described in the text. Differences between two survival curves was evaluated by Log-rank (Mantel-Cox) test. Blue and red curves represent patients with low or high levels of circulating target, respectively. HR = hazard ratio computed to determine the magnitude of differences between two curves; p-value was computed by log-rank test; 95% CI = 95% confidence interval.
A similar approach was exploited to assess association between circulating lncRNA levels and LFS. Kaplan-Meier curves showed that high levels of LINC01268 associated with a shorter LFS compared to patients with lower levels (p = 0.0063, log-rank test; HR = 9.685) (Figure 5b). No correlation with LFS was instead observed for MALAT1 and GAS5 levels (p = 0.5716 and p = 0.1534, respectively, log-rank test) (Figure 5d,f) or for the remaining lncRNAs studied ( Figure S3b,d,f,h).
Together these results demonstrated that circulating lncRNAs LINC01268, GAS5 and MALAT1, all detected at higher levels in MF patients if compared to HDs and associated to several features detrimental for MF prognosis, correlated with a shorter OS. Patients' cohort was stratified into two groups (low and high) according to the plasma levels of target lncRNA, as described in the text. Differences between two survival curves was evaluated by Log-rank (Mantel-Cox) test. Blue and red curves represent patients with low or high levels of circulating target, respectively. HR = hazard ratio computed to determine the magnitude of differences between two curves; p-value was computed by log-rank test; 95% CI = 95% confidence interval.

LINC01268 Plasma Level Is an Independent Prognostic Factor for OS and LFS
As mentioned before, we observed that plasma levels of LINC01268, MALAT1 and GAS5 correlated with DIPSS classification of MF patients (Figure 4). In order to evaluate the independent prognostic value of these lncRNAs we performed a DIPSS-adjusted multivariate analysis.
Kaplan-Meier curves demonstrated that OS was significantly reduced in patients displaying high plasma level of LINC01268 when considering DIPSS lowest-risk categories (Low and Intermediate-1) (p = 0.0069, log-rank test; HR = 10.11) (Figure 6a), but not DIPSS highest categories (Intermediate-2 and High) (Figure 6b). Moreover, multivariate analysis confirmed the independent prognostic relevance for OS of LINC01268 plasma level when considering DIPSS classification (HR = 2.104; confidence interval (CI) = 1.08-4.12; p = 0.0297) (Table 2). Similarly, LFS was reduced in the high LINC01268 group only in lowest DIPSS categories (p = 0.0360, log-rank test) (Figure 6c,d) and multivariate analysis confirmed the negative impact of higher LINC01268 on LFS independent from DIPSS prognostication (HR = 8.190; CI = 1.02-65.78; p = 0.0479) ( Table 2). Conversely, our analysis was not able to highlight an independent prognostic relevance for MALAT1 and GAS5 (Table 2).  Survival (a,b) and Leukemia-free Survival (c,d) in the settings of DIPSS categories according to LINC01268 plasma levels. Patients were dichotomized in two subgroups according to DIPSS classification: lowest categories (a,c) correspond to Low and Intermediate-1 classes while highest categories (b,d) include Intermediate-2 and High groups. Patients' cohort was then stratified according to low and high plasma levels of LINC01268, as previously described in the text. Differences between survival curves were evaluated by Log-rank (Mantel-Cox) test. Blue and red curves represent patients with low or high levels of circulating LINC01268, respectively. HR = hazard ratio computed to determine the magnitude of differences between two curves. † = No event occurred in the Low LINC01268 subgroup; p-value was computed by log-rank test; 95% CI = 95% confidence interval.  Survival (a,b) and Leukemia-free Survival (c,d) in the settings of DIPSS categories according to LINC01268 plasma levels. Patients were dichotomized in two subgroups according to DIPSS classification: lowest categories (a,c) correspond to Low and Intermediate-1 classes while highest categories (b,d) include Intermediate-2 and High groups. Patients' cohort was then stratified according to low and high plasma levels of LINC01268, as previously described in the text. Differences between survival curves were evaluated by Log-rank (Mantel-Cox) test. Blue and red curves represent patients with low or high levels of circulating LINC01268, respectively. HR = hazard ratio computed to determine the magnitude of differences between two curves. † = No event occurred in the Low LINC01268 subgroup; p-value was computed by log-rank test; 95% CI = 95% confidence interval. Table 2. Results of multivariate regression analysis for overall survival and risk of transformation into acute myeloid leukemia. Analysis of the prognostic impact of high plasma level of each lncRNA was performed after stratification of samples according to DIPSS risk category. Significant p-value (p < 0.05) are represented in bold.

Overall Survival
Transformation into Acute Leukemia Taken together, results of multivariate analysis demonstrated that LINC01268 plasma level in MF patients is an independent factor for both OS and LFS analysis.

Discussion
Myelofibrosis is a hematological neoplasm which originates from somatic mutations emerging in the hematopoietic stem cells compartment. Once occurring in the JAK2, MPL and CALR genes, these mutations might be considered as driver events. However, about 5% of patients display none of them and are thus classified as triple negative. A number of other mutations might be found in MF patients, such as those affecting ASXL1, IDH1/2, EZH2 and SRSF2 genes referred to as HMR, since these are associated with a diminished OS and/or LFS [35]. According to its onset, MF might be primary or secondary when progressed from PV and ET. In addition, a pre-fibrotic and an overt stage might be identified within PMF according to the degree of bone marrow fibrosis [1]. Although establishing different entities, they are managed following the same therapeutic strategies and ASCT is considered so far as the sole curative treatment, even if it is still associated to frequent graft-related deaths and comorbidities [12]. Thus, the identification of specific prognostic predictors, such as new biomarkers, might integrate current prognostic scoring systems [8][9][10] to enable a better management of MF patients.
In this light, expression profiling studies have been extensively performed in recent years in order to provide new signatures suitable to integrate or substitute standard prognostic and/or predictive factors for cancer patients. In particular, gene expression pattern has been exploited to discover new additional prognostic indicators and to identify the best therapy for breast cancer [36], ovarian cancer [37] and melanoma [38]. A 17-gene signature was also proposed to predict survival in AML [39]. Recently, we demonstrated that gene expression profiling of granulocytes provides complementary prognostic information to manage MF patients [13].
MF is characterized by an intricated interplay between malignant hematopoietic stem cells and stromal cells in the bone marrow microenvironment, as synthesized by the "bad seeds in bad soil" concept. Several types of molecules have been identified as involved in this crosstalk, such as cytokines and growth factors together with oxygen and calcium deregulation [40]. By contrast, an exhaustive characterization of the putative involvement of circulating nucleic acids (CNAs) in MF pathogenesis has not been provided so far.
The fibrotic conversion of bone marrow occurring in MF induces a release into circulation of CD34+ cells normally residing in bone marrow [41]. Circulating neoplastic cells derived from neoplastic clones might thus discharge an abnormal number of secreted molecules into plasma.
In order to assess the potential significance as prognostic indicators of disease progression, in this study we decided to measure the levels of circulating lncRNAs since their secondary structures confer a high stability in body fluids [23]. In fact, due to their high specific pattern of expression and functional diversity in a variety of solid and hematological disorders, lncRNAs have promising application in cancer diagnosis, prognosis and therapy.
The presence of RNA molecules in plasma was firstly described in 1999 [42,43], but several factors, such as their low abundance in liquid samples, have prevented their adoption as cancer biomarkers for a long time. Circulating small non-coding RNAs (i.e., microRNAs) are more frequently investigated in research works since interaction with AGO2 protein preserves them from RNase activity [21].
Notably, despite the increasing number of publications about the putative utility of circulating non-coding (ncRNAs) as clinical biomarkers, urinary lncRNA PCA3 is the only circulating RNA approved by the FDA for molecular testing and used to complement PSA for management of early prostate cancer [44].
This is the first report in which plasma levels of several lncRNAs have been assessed in MF patients. To this purpose, we selected a list of 13 lncRNAs deregulated in hematological malignancies both from literature and from preliminary gene-expression data in CD34+ cells. Among them, we identified seven lncRNAs, namely LINC01268, GAS5, MALAT1, LINC00899, TUG1, NEAT1, CDKN2B-AS1, as increased in plasma of MF patients if compared with healthy controls. Our results clearly demonstrated that patients' stratification based on plasma levels of these circulating lncRNAs are in agreement with the presence of detrimental prognostic factors. In particular, high levels of LINC01268, GAS5 and MALAT1, correlated with detrimental features of MF, such as high WBC count, increased number of circulating CD34+ cells, marked LDH activity, presence of splenomegaly and a more severe grade of bone marrow fibrosis. By contrast, low levels of these lncRNAs correlate with the absence of these detrimental features and with a better OS, as described in [45]. Moreover, the percentage of patients displaying ≥1 HMR mutation was increased in group of patients displaying high levels of LINC01268 or MALAT1 and, more specifically, patients harboring ASXL1 mutation are enriched in the group with high levels of LINC01268, MALAT1 or GAS5.
Strikingly, the group of patients with high GAS5, MALAT1 or LINC01268 were characterized by a shorter OS, the latter leading in addition to a worse LFS. Once stratified according to DIPSS scoring system, patients classified into the High category were more frequently characterized by high levels of LINC01268, GAS5 or MALAT1 while low levels of these lncRNAs were more frequent in DIPSS Low class. Therefore, we demonstrated that lncRNA plasma levels correlated with contemporary prognostic scoring system. In order to evaluate whether lncRNAs might provide additional independent prognostic information we performed multivariate analysis. Our results demonstrated that high plasma levels of LINC01268 was also able to predict a shorter OS and LFS independently from the DIPSS classification, in particular when considering DIPSS lowest-risk categories (Low and Intermediate-1). This might be important because only observation or palliative therapy in the presence of specific symptoms is suggested for patients belonging to DIPSS Low or Intermediate-1 classes. More effective treatments, such as ASCT, are to be considered only for patients belonging to Intermediate-2 or High risk categories [46]. Our results suggest that evaluation of LINC01268 plasma levels might be used to identify patients at higher risk within DIPSS lowest categories thus guiding more appropriate clinical decision.
From this study, it emerges that GAS5, MALAT1 and LINC01268 are differentially expressed in MF and correlated to detrimental clinical features in patients; therefore, we speculate that these lncRNAs could be involved in disease pathogenesis.
GAS5 (Growth Arrest Specific 5) has been commonly considered as a tumor suppressor, being downregulated in tissues from several solid tumors (as reviewed by [47]). Nonetheless, in line with our results, circulating levels of GAS5 were increased in plasma of patients affected by some tumors, like mesothelioma [48], and other pathological conditions, like osteoporosis [49].
The nuclear MALAT1 (Metastasis-Associated Lung Adenocarcinoma Transcript 1) has been described to display contradictory action both as lncRNA promoting and suppressing metastasis [50,51].
LINC01268 has been recently described as being upregulated in AML samples [52,53]. In a first work LINC01268 has been described as acting as competing endogenous RNA (ceRNA) via sponging the onco-suppressor miR-217 which in turns inhibits at the posttranscriptional level SOS1 mRNA. By modulating miR-217/SOS1 regulatory axis, LINC01268 promotes cell growth and inhibits apoptosis [52]. The downregulation of miR-217 in AML cells was described also by Xiao et al., being proposed acting as a tumor suppressor by directly targeting KRAS [54].
LINC01268 involvement in solid cancer is still poorly described and, to our knowledge, a description of circulating levels has not been assessed so far. Thus, we tried to design a network of gene expression regulation by LINC01268, which demonstrates that its overexpression is one of the mechanisms underlying leukemogenesis in AML and the induction of expression of inflammatory cytokines (Figure 7).
LINC01268 involvement in solid cancer is still poorly described and, to our knowledge, a description of circulating levels has not been assessed so far. Thus, we tried to design a network of gene expression regulation by LINC01268, which demonstrates that its overexpression is one of the mechanisms underlying leukemogenesis in AML and the induction of expression of inflammatory cytokines (Figure 7). LINC01268 has been recently described as being upregulated in AML samples [52,53]. In a first work LINC01268 has been described as acting as competing endogenous RNA (ceRNA) via sponging the onco-suppressor miR-217 which in turns inhibits at the posttranscriptional level SOS1 mRNA. By modulating miR-217/SOS1 regulatory axis, LINC01268 promotes cell growth and inhibits apoptosis [52]. The downregulation of miR-217 in AML cells was described also by Xiao et al., being proposed acting as a tumor suppressor by directly targeting KRAS [54].
An upregulation of LINC01268 (therein referred as LOC285758) was described also to induce invasion and viability in AML cells by directly sponging miR-204, which inhibited N-Cadherin and TWIST-1 and increased E-Cadherin [53]. A downregulation of miR-204 in AML was also described by the work of Butrym et al., where decreased serum level of miR-204 have been correlated with a poor prognosis of AML [55]. Once again, miR-204 has been described to exert an anti-apoptotic effect by inhibiting BIRC6 mRNA in AML cells [56]. Beyond its putative involvement in cancer pathogenesis, LINC01268 (also known as lnc-MARCKS or ROCKI) acts as a master regulator of inflammatory response in macrophages, inducing pro-inflammatory cytokines and chemokines by acting in cis on MARCKS promoter [57].

Conclusions
As a whole, this work demonstrated that increased levels of circulating LINC01268, GAS5 or MALAT1 are associated with a number of clinical and molecular detrimental features and correlate with an inferior OS in MF patients. Notably, multivariate analysis confirmed that LINC01268 plasma levels can be considered an independent variable. If the prognostic value of LINC01268 is confirmed in future in an independent cohort, it could be evaluated in a perspective clinical study and then possibly integrated in contemporary prognostic models. An upregulation of LINC01268 (therein referred as LOC285758) was described also to induce invasion and viability in AML cells by directly sponging miR-204, which inhibited N-Cadherin and TWIST-1 and increased E-Cadherin [53]. A downregulation of miR-204 in AML was also described by the work of Butrym et al., where decreased serum level of miR-204 have been correlated with a poor prognosis of AML [55]. Once again, miR-204 has been described to exert an anti-apoptotic effect by inhibiting BIRC6 mRNA in AML cells [56]. Beyond its putative involvement in cancer pathogenesis, LINC01268 (also known as lnc-MARCKS or ROCKI) acts as a master regulator of inflammatory response in macrophages, inducing pro-inflammatory cytokines and chemokines by acting in cis on MARCKS promoter [57].

Conclusions
As a whole, this work demonstrated that increased levels of circulating LINC01268, GAS5 or MALAT1 are associated with a number of clinical and molecular detrimental features and correlate with an inferior OS in MF patients. Notably, multivariate analysis confirmed that LINC01268 plasma levels can be considered an independent variable. If the prognostic value of LINC01268 is confirmed in future in an independent cohort, it could be evaluated in a perspective clinical study and then possibly integrated in contemporary prognostic models.
To our knowledge, this is the first study describing the profile of circulating lncRNAs in plasma of MF patients and focusing on their putative role as biomarkers in clinical practice. In addition, our work sets the basis for further studies regarding the mechanism(s) underlying the role of these lncRNA(s) in MF pathogenesis.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cancers13194744/s1, Table S1: List of lncRNAs whose expression has been evaluated in CD34+ cells, Table S2: List of lncRNAs evaluated in plasma samples, Table S3: Clinical and molecular features of MF patients included in our dataset, grouped according levels of NEAT1 and CDKN2B-AS1, Figure S1: Analysis of stability of putative endogenous controls, Figure S2: Distributions of circulating lncRNAs in MF patients, Figure S3