A compendium of long non-coding RNAs transcriptional fingerprint in multiple myeloma

Multiple myeloma (MM) is a clonal proliferation of bone marrow plasma cells characterized by highly heterogeneous genetic background and clinical course, whose pathogenesis remains largely unknown. Long ncRNAs (lncRNAs) are a large class of non-protein-coding RNA, involved in many physiological cellular and genomic processes as well as in carcinogenesis and tumor evolution. Although still in its infancy, the role of lncRNAs in MM is progressively expanding. Besides studies on selected candidates, lncRNAs expression at genome-wide transcriptome level is confined to microarray technologies, thus investigating a limited collection of transcripts. In the present study investigating a cohort of 30 MM patients, a deep RNA-sequencing analysis overwhelmed previous array studies and allowed the most accurate definition of lncRNA transcripts structure and expression, ultimately providing a comprehensive catalogue of lncRNAs specifically associated with the main MM molecular subgroups and genetic alterations. Despite the small number of analyzed samples, the high accuracy of RNA-sequencing approach for complex transcriptome processing led to the identification of 391 deregulated lncRNAs, 67% of which were also detectable and validated by whole-transcript microarrays. In addition, we identified a list of lncRNAs, with potential relevance in MM, co-expressed and in close proximity to genes that might undergo a cis-regulatory relationship.

human lncRNAs and analyzed their genomic organization, modifications, cellular localizations and tissue expression profiles in different human cell line. LncRNAs contribute to several processes, e.g. maintenance of genomic integrity, X-chromosome inactivation, transcriptional regulation, genomic imprinting, cell differentiation and development 15,16 . Several lncRNAs have also been described to contribute to tumor formation and/or progression, as well as to metastatic processes, in many solid and hematologic tumors 17,18 , showing either oncogenic or tumor suppressive function.
Although the investigation of lncRNA in MM is still in its beginnings, our understanding of their role is progressively expanding. One of the most investigated lncRNA is MALAT1, deregulated in many solid tumors with a putatively oncogenic function 19,20 . MALAT1 is overexpressed in MM, where it has been shown to predict tumor progression 21 . Recent researches in MM have been focused on single lncRNAs already known as involved in different types of cancers, such as MEG3 functioning as tumor suppressor through both p53-dependent and p53-independent mechanisms 22 , or CRNDE overexpressed in association with poor clinical outcome 23 . Besides studies on selected candidates, lncRNA expression at genome-wide transcriptome level has been scarcely investigated in MM and the only two efforts reported so far are based on microarray data. In particular, Zhou et al. 24 investigated a repertoire of 2,330 lncRNAs in a publicly available clinically annotated cohort of 559 MM patients generating a four-lncRNA prognostic signature. In a previous study, our group analyzed the transcriptional patterns of 1,852 lncRNAs in 259 patients affected by the different forms of PC dyscrasia at onset, included in proprietary and publicly available datasets, identifying a series of deregulated lncRNAs associated either with disease progression or distinct molecular subgroups of MM 25 . However, these studies are limited to the detection of a relatively small number of sequences queried by the arrays, which were primarily designed to detect the coding transcriptome. Next-generation RNA sequencing (RNA-seq) addresses this shortcoming, but to date such an approach has not yet been pursued in MM.
In the present study, we investigated the lncRNA expression profiling in MM patients by RNA-seq, with the aim of providing a first exhaustive catalogue of lncRNAs specifically associated with the main molecular subgroups and genetic alterations in MM. Furthermore, we defined a repertoire of lncRNAs possibly involved in MM, as they meet the requirements of being both co-expressed and in close proximity to genes that have been described as relevant to this neoplasia, thus suggestive of a cis-regulatory relationship.
Overall, such a compendium and the free availability of RNA-seq data may provide the scientific community with valuable references for future research into the involvement of lncRNAs in MM.

Results
LncRNAs expression profile in multiple myeloma. The expression profile of lncRNAs has been investigated by RNA-seq in a cohort of 30 MM patients at diagnosis, whose molecular features were representative of those mainly characterizing the disease (Table 1).
We used a custom pipeline, based on the GENCODE encyclopedia that considered only those genes with unambiguously mapped transcripts, that allowed to annotate 14,202 lncRNAs; among them, we investigated the 9,540 lncRNAs detectable upon removal of those unexpressed across the whole dataset. Overall, lncRNAs are scarcely expressed. Indeed, for each lncRNA the sum of the read counts in the 30 samples spans a wide range of values (from 2 to 6,707,843; median: 57). However, 86% of the 9,540 lncRNAs have average read counts <30, whereas only 1% of lncRNAs show average values > 500. Notably, 12 lncRNAs are very highly expressed displaying an average read counts >5000, counting 64% of the reads assigned to lncRNAs (Supplementary Table S1); in particular, this group includes NEAT1, MALAT1, MIAT and TUG1 frequently deregulated in malignant B-cells 26 . Based on the rationale that a single cis-acting molecule might be able to target effectively a neighboring locus, thus suggesting that even low expressed lncRNAs may have a key regulatory role 27 , we considered all the 9,540 detectable lncRNAs for subsequent investigations.
To identify MM patient subgroups, we used an unsupervised-learning method based on expression data. This analysis showed clusters of common global lncRNAs transcriptional patterns that were associated with the major and prognostically relevant molecular features, namely t(11;14), t(4;14), MAF gene translocations or HD status. In fact, unsupervised analysis of the 500 lncRNAs with the highest variation coefficient clearly showed that MM molecular subtypes were mainly and significantly clustered together (Fig. 1a) Fig. 2). The five most significant differentially expressed lncRNAs in each comparison are reported in Table 2. As regards MM carrying t (11;14), of note, the three most significantly downregulated lncRNAs belonged to a cluster of 6 transcripts located in a region of about 332 kb at 19q12 ( Fig. 3a and Supplementary Table S2). In addition, this group showed also the downregulation of MIAT (Fig. 3b-d), a well-known lncRNA already reported as involved in different cancers.

Identification of lncRNA signatures associated with genetic lesions or somatic mutations.
Other genetic alterations occur at high frequency in MM and were associated by others and us to specific transcriptional profiles. Information was available for the 30 MM sample on numerical alterations and secondary events, i.e. somatic mutations (Table 1): therefore, we queried the RNA-seq dataset to evaluate the occurrence of differentially expressed lncRNAs in those MM genetic subtypes (Supplementary Table S3, with the exclusion of FAM46C and P53 due to the low number of samples). In MM patients with 1q gain, we found the significant modulation of 12 lncRNAs (4 down-and 8-upregulated), two of which located on chromosome 1q. A list of 109 lncRNAs (31 down-and 78 up-regulated) distinguished del(13)-positive from wild-type patients; notably, 7 of 31 downregulated lncRNAs (23%) are located on chromosome 13. Finally, only two lncRNAs have been found downregulated in del(17)-MM.
As regards patients harboring the mutations of BRAF, NRAS, KRAS, or DIS3 mutations, we identified the upregulation of 97 lncRNAs in DIS3 mutated samples, whereas 6 lncRNAs are upregulated in the samples grouped according to the presence of MAPK-pathway genes mutations.

Selection of lncRNAs potentially relevant in MM.
After the annotation process, we established a set of criteria to recognize the lncRNAs potentially relevant in MM biology. In particular, we investigated the levels of expression of lncRNAs localized in proximity to genes associated with MM, based on the recurrent evidence that the transcription of mRNAs and lncRNAs appears to be closely regulated, leading to a cis-regulatory relationship between the two transcripts [28][29][30] . For this purpose, a list of 707 genes mapped to GRCh38 primary assembly and associated with MM (from now on, referred as "MM-genes") was downloaded from NCBI database (https://www.ncbi.nlm.nih.gov/gene, Supplementary Table S4). We analyzed the genomic context of the 707 MM-genes. In the boundaries of 409 of them, we found at least one of the 9,540 lncRNAs mapped within 4 Mb (Supplementary Table S4). Next, for each lncRNA/MM-gene couples we assessed the correlation of their  (Table 3).

Discussion
In the present study, we have provided an unprecedented view of the lncRNAs expression in MM. As it occurred for mRNAs, miRNAs, and snoRNAs 6,7,31,32 , the natural clustering of whole lncRNAs transcriptional configuration is significantly associated with the major molecular prognostic alterations in MM, namely 11q13, 4p16, 16q23/20q12 chromosomal translocations, or HD status. In details, for each MM subtype we defined a specific and exhaustive lncRNAs expression signature based on the 14,202 lncRNAs currently annotated on GENCODE database. In a previous study concerning about 1,800 lncRNAs detectable by microarrays, we had reported a number of differentially expressed lncRNAs among the same MM subgroups. However, with very few exceptions (MATN1-AS1 upregulated in MM with t (11;14), and CRYM-AS1 and LINC00158 upregulated in MAF translocated patients), RNA-seq data scarcely overlapped with our previous data. This discrepancy can be explained in all likelihood by two reason: first, the array annotation was based on a previous version of the LNCipedia repository (https://lncipedia.org) 33 that had included pseudogenes and miscellaneous RNA within lncRNA transcripts, which are conversely excluded in the current study focused on transcripts annotated as "pure" long non-coding RNA. Second, very little is still known about the processing and the prevalence of alternative transcripts for many lncRNAs, whose "splicing" products are often roughly defined and/or based on predictions. While RNA-seq allowed to evaluate the non-coding genes in their full extension (according to the provided annotations), microarrays evaluation is probe-position dependent and might be therefore affected by the number of transcripts in the queried region. This last aspect undoubtedly reinforce the highest accuracy of RNA-seq data for complex transcriptome processing. We are aware that the number of samples analyzed in this study does not allow drawing definitive conclusion, all the more true in that myeloma patients may share different primary/secondary molecular alterations. We kept this in close consideration when the cohort was selected, aimed at being representative of the major genetic lesions and avoiding as much as possible that confounding variables might affect data in differential analysis (graphical legend to the Fig. 1a). To further overcome these limitations, lncRNAs expression evaluated by RNA-seq technology was validated by high-density arrays, overall leading to the definition of a comprehensive background for future investigations of lncRNAs in plasma cell dyscrasias.
Considering the lncRNAs expression signatures, no information is currently available for the majority of the lncRNAs identified. Among the five most significant lncRNAs found differentially expressed in each comparison (Table 2), the well-known lncRNA MIAT resulted specifically downregulated in MM carrying t (11;14). Originally identified within a susceptible locus for myocardial infarction on chromosome 22q12.1, MIAT was then characterized as the RNA component of specific nuclear bodies where it may affect RNA splicing, ultimately regulating gene expression 34 . Recently Sattari et al. 35 found MIAT upregulation in leukemia/lymphoma lymphoid lineage with mature B-cell phenotype; interestingly, they demonstrated a higher incidence of MIAT upregulation in aggressive types of CLL and worst clinical outcome. In addition, this study described a positive feedback regulatory loop between MIAT and OCT4, acting on evading apoptotic cell death in malignant mature B cells. Overall, these findings suggest an involvement of MIAT in supporting proliferation of the malignant mature B-cells. In this perspective, lower MIAT expression in t(11;14)-positive patients might be associated with the better prognosis associated with this MM subtype 36 .
Among the most significant lncRNAs defining the signature of HD-MM, we found the downregulation of ASHL1-AS1 and KB-1471A8.1. Both lncRNAs resulted also from the analyses aimed at identifying lncRNAs that are located in proximity to, and concordantly expressed with genes important in the context of MM pathology (Table 3). In details, ASHL1-AS1 maps 1,89 Mb telomeric to ILF2, overexpressed in MM as a result of 1q21 amplification. ILF2 overexpression deregulates homologous recombination (HR) by stabilizing the mRNA splicing of critical HR effectors, which enables genomic instability, promotes adaptive mechanisms to genotoxic stress, and enhances cell survival, thereby promoting drug resistance and disease progression 37 . As regards KB-1471A8.1, it maps at 8q24 antisense to the 5′ region of DEPTOR, a crucial gene in the maintenance of the terminal differentiation of MM cells 38 . Since the overexpression of DEPTOR in MM has been associated with MAF translocations and the expression of CCND1 and CCND3 genes 39 , the downregulation of KB-1471A8.1 in HD-MM further suggest a cis-regulatory connection with DEPTOR.
Finally, among the most significant lncRNAs deregulated in MAF translocated MM, our data unraveled the MIR222HG sequence, from the maturation of which originate microRNAs 221 and 222 that were found accordingly downregulated in this MM subtype 31 . MIR222HG is located at Xp11 about 1.7 Mb telomeric to the TIMP1 gene encoding an inhibitor of metalloproteinases. As the balance between metalloproteinases and their inhibitors, including TIMP1, largely influences cell adhesion, proteolytic shedding, and cell signaling, it will be of great interest to clarify the putative regulation of TIMP1 expression by MIR222HG.
Overall, to our knowledge our study provides the first comprehensive catalogue by RNA-seq of lncRNAs in MM, which is highly beneficial as a valuable reference for future research on their involvement into the pathogenesis of the disease.

Methods
Samples. The molecular features of the 30 patients at diagnosis included in the study cohort are shown in Table 1. PCs purification has been previously described and led to >90% enrichment in all samples [Mattioli, Oncogene 2005]. According to already reported FISH procedure 40 , eight samples showed the t(11;14) translocation, with the consequent overexpression of either CCND1, and a non-hyperdiploid (HD) status; 8 MM were HD; seven patients showed high CCND2 levels and the presence of the t(4;14) translocation; and four expressed the highest levels of CCND2 in association with either the t(14;16) or t (14;20) translocations. Information on 17p13 and 13q14 deletions, and gain of 1q arm was also available. Mutation of BRAF, NRAS, KRAS, P53, FAM46C and DIS3 were investigated by next-generation sequencing [41][42][43][44]  RNA sequencing. Total RNA was extracted from purified PCs by using Trizol reagent. Quantitative assessment of the RNA was performed using Nanodrop ND-1000 Biophotometer (NanoDrop Technologies): the minimum OD 260/280 ratio to be considered acceptable is 1.98-2.10. Four-hundred ng of total RNA were used to prepare paired-end (PE) cDNA libraries using the TruSeq ® RNA Sample Preparation kit for total RNA (Illumina).
The libraries were sequenced to obtain strand-specific 100 bp PE reads on a HiSeq. 2000 (Illumina). Reads were aligned to the human genome using STAR under default conditions and Gencode v25 GTF file. STAR aligner was based on splice junctions from the Ensembl database version 87. Transcript abundance was estimated using featureCounts (default parameters). FPKM (Fragments Per Kilobase Million) quantification was performed on sorted BAM files using cufflinks default procedure. Differentially expressed genes were identified using DeSeq at FDR < 0.01, provided that expression across the whole dataset was not null. Quality Control (QC) analysis was performed using multiqc tool and the QC metrics were comparable for all samples. The annotation allowed to detect 14,202 lncRNAs, including the following Gencode biotypes: lincRNA, antisense, bidirectional promoter lncRNA, sense intronic, sense overlapping, 3′ overlapping ncRNA. The expression filter retained 9,540 lncRNAs in our dataset.

Gene Expression Profiling. Thirty MM samples and four normal controls (purchased from Voden, Medical
Instruments IT) were profiled onto GeneChip ® Human Gene 2.0 ST arrays (Affymetrix Inc., Santa Clara, CA).
Total RNA samples were processed according to manufacture's procedure. Normalized expression values were obtained using Robust Multi Array Average (RMA) procedure. A custom annotation pipeline was applied that combined GENCODE v25 (Ensembl v87) annotations with the CDF (Chip Definition File) version 21 for gene annotations freely available at http://brainarray.mbni.med.umich.edu/Brainarray/Database/CustomCDF/21.0.0/ genecodeg.asp, in order to withdraw probes that map to regions where ambiguous detection due to transcript overlap might occur. Therefrom, the expression levels of Ensembl genes specific for 10138 unique lncRNAs were obtained.
All the data have been deposited in the NCBI Gene Expression Omnibus database (GEO; http://www.ncbi. nlm.nih.gov/geo) and are accessible under accession #GSE109116.
Statistical analysis. Pearson's correlation as distance and centroid linkage were used in hierarchical agglomerative clustering analysis. Conventional statistical tests were applied as reported in the manuscript using standard packages for R software.