Abstract
Long non-coding RNAs (lncRNAs) are post-transcriptional and epigenetic regulators, whose implication in neurodegenerative and autoimmune diseases remains poorly understood. We analyzed publicly available microarray data sets to identify dysregulated lncRNAs in multiple sclerosis (MS), a neuroinflammatory autoimmune disease. We found a consistent upregulation in MS of the lncRNA MALAT1 (2.7-fold increase; meta-analysis, P = 1.3 × 10−8; 190 cases, 182 controls), known to regulate alternative splicing (AS). We confirmed MALAT1 upregulation in two independent MS cohorts (1.5-fold increase; P < 0.01; 59 cases, 50 controls). We hence performed MALAT1 overexpression/knockdown in cell lines, demonstrating that its modulation impacts on endogenous expression of splicing factors (HNRNPF and HNRNPH1) and on AS of MS-associated genes (IL7R and SP140). Minigene-based splicing assays upon MALAT1 modulation recapitulated IL7R and SP140 isoform unbalances observed in patients. RNA-sequencing of MALAT1-knockdown Jurkat cells further highlighted MALAT1 role in splicing (approximately 1100 significantly-modulated AS events) and revealed its contribution to backsplicing (approximately 50 differentially expressed circular RNAs). Our study proposes a possible novel role for MALAT1 dysregulation and the consequent AS alteration in MS pathogenesis, based on anomalous splicing/backsplicing profiles of MS-relevant genes.
Introduction
Multiple sclerosis (MS) is an autoimmune disease in which peripherally activated T cells attack myelin sheets, leading to chronic inflammation, neuronal loss and axonal damage (1). The disease is classified in three clinical subtypes based on the different severity and progression of symptoms. The relapsing–remitting (RR) is the most common form (characterizing 80% of MS patients), and it is associated with attacks followed by complete or partial recovery periods; in most RR-MS patients (65%), chronicity and worsening of symptoms occur, leading to the secondary progressive course. The remaining 20% of patients experience the primary progressive (PP) form, with a chronic increase of symptoms and disability starting from the onset of the disease (2). MS is a complex disease in which the interplay between genetic and environmental factors is known to contribute to the pathogenesis (3); however, the molecular events responsible for MS onset are still poorly understood. Recently, alterations in the alternative splicing (AS) process are increasingly recognized as a possible pathogenic mechanism underlying MS (4–7).
AS is a post-transcriptional mechanism that amplifies the complexity of the transcriptome through the expression of a subset of different mRNAs from single genes, thus potentially generating distinct protein isoforms (8). AS regulation is mediated by regulatory cis-acting sequence elements (exonic or intronic splicing enhancers and silencers), recognized by trans-acting factors, such as serine–arginine (SR) proteins, heterogeneous nuclear ribonucleoproteins (hnRNPs) and other tissue-specific factors (9). AS alterations are caused by mutations in cis-acting elements, which impact on the splicing pattern of the involved genes, or by alterations of trans-acting factors, which affect the AS profile of multiple genes (10). Both these mechanisms have been demonstrated to play a role in MS. Indeed, functional studies showed that several MS-associated genetic variants alter the AS pattern of the relevant gene and eventually impact on the balance between distinct protein isoforms (11–15). One paramount example is the IL7R (interleukin 7 receptor) gene, whose MS-associated polymorphism affects the ratio between the membrane bound and the soluble form of the receptor and, consequently, the IL7R-mediated immune signaling pathway (12). Concerning trans-acting factors, a general alteration of proteins involved in RNA processing was found in RR-MS patients' blood by RNA-seq analysis (7). At the same time, our group evidenced a global alteration in the expression of splicing factors (including, among the others, CELF1, HNRNPH1 and SRSF1) in RR-MS patients' blood, by performing a meta-analysis of seven publicly available microarray studies (6). Interestingly, novel somatic variants of the RNA-binding protein (RBP) hnRNPA1 have been found in MS patients, and antibodies against this protein, detected in MS sera, were described to contribute to neurodegeneration (16,17).
Concerning AS regulation, the role of additional players, including chromatin structure and modifications, and long non-coding RNAs (lncRNAs), is clearly emerging (18). LncRNAs are a class of long (>200 nt) RNAs with no significant protein-coding capacity, being able to interact with other RNAs, proteins or DNA, thus regulating gene expression at many levels (19,20). LncRNAs mainly exert AS regulation by (i) RNA–RNA duplexing that masks cis-acting elements, as demonstrated for natural antisense transcripts (21–23) and (ii) modulating splicing-factor functional levels, through direct interaction, regulation of their phosphorylation status or their sequestration (23–26). Variations in the expression of lncRNAs have been already associated with several diseases, including autoimmune and neurodegenerative disorders (27,28). Considering MS, two explorative microarray studies specifically focused on lncRNAs in RR-MS patients: one evaluated the serum expression of lncRNAs already linked with immune response (29); the other globally detected lncRNAs and mRNAs in peripheral blood mononuclear cells (PBMCs) (30). In both studies, several lncRNAs were shown to be differentially expressed in MS patients compared with healthy controls.
In this work, we aimed at identifying differentially expressed lncRNAs in blood cells from RR-MS patients and controls, by a comparative analysis of publicly available transcriptome data sets. We highlighted a consistent upregulation in MS patients of the lncRNA MALAT1 (metastasis-associated lung adenocarcinoma transcript 1), already known to modulate AS (24). We hence dissected the MALAT1 role in the regulation of both splicing-factor expression and selected MS-associated AS events, providing evidence of its possible involvement in MS pathogenesis.
Results
The lncRNA MALAT1 is upregulated in MS patients
To identify differentially expressed lncRNAs in RR-MS patients, we performed an in silico analysis, using seven microarray data sets selected in our previous work (6) and containing blood cell expression data derived from a total of 190 patients and 182 healthy controls. For each data set, only the lncRNAs that survived multiple testing correction were considered. We found differentially expressed lncRNAs in 4 out of 7 data sets, for a total of 168 dysregulated genes (Supplementary Material, Table S1). In particular, 20 genes were shared by 2 data sets, and 2 genes were shared by 3 data sets, i.e. MALAT1 and HCG18 (human leukocyte antigen (HLA) complex group 18; Fig. 1A). A meta-analysis, performed by the NetworkAnalyst tool (31), revealed that MALAT1 is consistently upregulated in MS cases (mean fold change, 2.70; combined t-stat, 55.387; combined P-value, 1.33 × 10−8), whereas HCG18 did not survive the threshold of significance.
MALAT1 is upregulated in MS patients. (A) Venn diagram showing the lncRNAs shared at least by two data sets. The genes shared by three out of three data sets are marked in bold. The number of genes present only in one data set is also displayed. (B) Boxplots showing expression levels of the lncRNA MALAT1 measured by semi-quantitative real-time RT-PCR in PBMCs of an Italian and of a US MS case–control cohort. Boxes define the interquartile range; the thick line refers to the median. Results were normalized to expression levels of the HMBS housekeeping gene and were presented as rescaled values. The number of subjects belonging to each group is indicated (N). Significance levels of t-tests are shown. *P < 0.05; **P < 0.01.
We then evaluated MALAT1 expression levels in PBMCs of an Italian cohort of 23 RR-MS patients and 30 healthy controls. We found a 1.4-fold upregulation (P = 0.024) of the lncRNA in RR-MS patients (Fig. 1B), confirming the meta-analysis results. MALAT1 upregulation was further validated in an independent US case–control cohort, including 29 RR-MS, 7 PP-MS patients and 20 healthy controls; indeed, MALAT1 expression levels were increased in RR-MS (1.6-fold, P = 0.016) and even more in PP-MS patients (2.3-fold, P = 0.0064), compared with healthy controls (Fig. 1B).
MALAT1 modulation affects the expression of splicing regulatory genes in a cell-type specific manner
Due to the involvement of MALAT1 in AS regulation (24), we decided to investigate the effects of MALAT1 modulation on the expression of splicing regulatory genes in different cell lines. We first performed MALAT1 knockdown experiments in HeLa, HEK293 and SH-SY5Y cells, and we measured the expression levels of HNRNPF, HNRNPH1, CELF1 and SRSF1 genes by semi-quantitative real-time reverse-transcription polymerase chain reaction (RT-PCR) assays. In these experiments, we obtained MALAT1 silencing levels comprised between 60% and 80% (Supplementary Material, Fig. S1A). In HEK293 cells, we found a significant reduction of the expression of HNRNPF (36% decrease, P = 0.0017), CELF1 (37% decrease, P = 0.026) and SRSF1 (40% decrease, P = 0.013). MALAT1 knockdown in SH-SY5Y cells led instead to a significant increase of HNRNPF and CELF1 expression levels (1.3- and 1.2-fold, respectively; P < 0.05). Finally, HNRNPF expression in HeLa cells was significantly decreased (36% decrease, P = 0.045; Fig. 2A).
MALAT1 modulation affects expression levels of splicing regulatory genes. Expression levels of HNRNPF (heterogeneous nuclear ribonucleoprotein F), HNRNPH1 (heterogeneous nuclear ribonucleoprotein H1), CELF1 (CUGBP Elav-like family member 1) and SRSF1 (serine and arginine rich splicing factor 1) in HEK293, SH-SY5Y and HeLa cells transfected either with an LNA antisense oligonucleotide GapmeR specific for MALAT1 (A) or with a plasmid expressing MALAT1 (B). In all cases, expression levels were measured by semi-quantitative real-time RT–PCRs, setting 1 as the value of the mock sample (transfected with a negative control or with an empty vector in knockdown and overexpression experiments, respectively). For normalization, expression levels of HMBS, ACTB and GAPDH housekeeping genes were evaluated and the best combination of two normalizers (measured through the NormFinder software, (https://moma.dk/normfinder-software) was used for each cell line (for HEK293 and HeLa cells, GAPDH and ACTB; for SH-SY5Y, HMBS and ACTB). Bars represent means + standard error of the mean (SEM) of three independent experiments, each performed in triplicate. Significance levels of t-tests are shown. *P < 0.05; **P < 0.01.
To corroborate these results, we evaluated the expression of the same splicing regulatory genes upon MALAT1 overexpression (overexpression >2.4-fold in all cell lines; Supplementary Material, Fig. S1B). In HEK293, levels of HNRNPF, HNRNPH1 and CELF1 significantly increased (2-, 1.6- and 1.9-fold, respectively; P < 0.05), showing, as expected, an opposite trend compared with knockdown experiments. HNRNPF expression levels were also significantly increased upon MALAT1 overexpression in HeLa (1.5-fold, P = 0.0085). Conversely, SH-SY5Y cells showed a significant downregulation of HNRNPF (18% reduction, P = 0.0065) and HNRNPH1 (20% reduction, P = 0.029) genes, confirming that MALAT1 modulation exerts an opposite effect on the expression of splicing genes in SH-SY5Y compared with HEK293 and HeLa cells (Fig. 2).
MALAT1 modulates AS events associated with MS
Given the modulation of splicing regulatory genes mediated by MALAT1, we evaluated if the dysregulation of this lncRNA could also reflect on MS-associated AS events. We focused on IL7R and SP140, two genes characterized by the presence of a MS-associated polymorphism that was demonstrated to regulate a specific AS event and to impact on the repertoire of protein isoforms. Concerning IL7R, the disease-associated C allele of the rs6897932 C>T polymorphism, located in exon 6, decreases the inclusion of this exon in the mature transcript (12). Similarly, the rs28445040 C>T MS-associated variant, located in SP140 exon 7, has been associated with an increased skipping of exon 7 from the mature transcript (15).
IL7R exon 6 and SP140 exon 7 (containing either alleles of the relevant polymorphism), with their respective intronic flanking sequences, were cloned into the hybrid α-globin-fibronectin minigene plasmid (Fig. 3A). The minigene constructs were co-transfected either with the MALAT1-expressing plasmid or with an empty vector into HEK293 cells. Analysis of chimeric transcripts showed that MALAT1 overexpression significantly increases IL7R exon 6 skipping, both for the C allele (from 82–93% of skipping level, P = 0.0010) and for the T allele (from 70–83%, P = 0.018; Fig. 3B). MALAT1 overexpression was also shown to increase the skipping of SP140 exon 7 from the mature transcript for both the alleles (C allele: from 10–26% of skipping, P = 0.0020; T allele: from 25–49% of skipping, P = 0.00020; Fig. 3C).
MALAT1 overexpression affects MS-associated AS events in splicing minigene assays. (A) Schematic representation of the α-globin-fibronectin minigene plasmid (pBS-KS modified). α-globin (HBA) and fibronectin (FN1) exons are represented by boxes and are approximately drawn to scale; introns are represented by lines. IL7R exon 6 and SP140 exon 7 are depicted in black with their respectively flanking regions; the positions of the polymorphisms rs6897932 and rs28445040 are indicated by an arrow. The primer couple used in the fluorescent-competitive RT-PCR assay is also shown; the forward primer is labeled with the 6-FAM fluorophore. (B) and (CC) The left panel shows a GeneMapper window, representing an example of the fluorescent products obtained in the splicing assay performed after transfection of the IL7R (B) construct or the SP140 (C) construct. The filled peaks, shaded in gray, correspond to the RT-PCR products; empty peaks represent the size standard (ROX-500 HD; Thermo Fisher Scientific). The schematic representation of the obtained products is also shown. IL7R exon 6 and SP140 exon 7 are depicted in black. F, full length; ∆ is followed by the SE. In the right panel, histograms indicate the quantitative analysis on fluorescence peak areas obtained in the RT-PCR assays performed on RNA extracted by HEK293 cells transfected with an empty vector (mock) or overexpressing MALAT1, for both the alleles of the IL7R rs6897932 (B) or of the SP140 rs28445040 (C) polymorphism. Bars represent means + SEM of three independent experiments, each performed in triplicate. Significance levels of t-tests are shown. *P < 0.05; **P < 0.01; ***P < 0.001.
The results obtained using minigene constructs were substantiated by analyzing the endogenous levels of the IL7R transcript in HEK293 cells under MALAT1 overexpression or knockdown conditions by fluorescent-competitive RT-PCR assays (Fig. 4A). Differently from minigene construct overexpression experiments, the AS repertoire of endogenous IL7R transcripts was more complex, including, besides the exon 6-skipped isoform (here referred to as ∆6), also the already reported transcript skipping exons 5 and 6 (here named ∆5∆6) and the newly identified isoform skipping exon 5 (∆5). Upon MALAT1 overexpression, the levels of exon 6-skipping isoforms significantly increased (2-fold, P = 0.022), whereas they decreased upon MALAT1 knockdown (18% decrease, P = 0.026). The increase in IL7R exon 6 skipping was also nicely confirmed in PBMCs of RR-MS patients (which exhibit higher levels of MALAT1, see Fig. 1B), compared with healthy controls (P = 0.047; Fig. 4A).
MALAT1 overexpression affects endogenous IL7R and SP140 AS isoform levels. The upper schematic gene representation displays the IL7R genomic region between exons 4 and 7 (A) and SP140 genomic region between exons 6 and 8 (B); exons are approximately drawn to scale; alternative exons are depicted in light gray. The position of the primer couple used in the fluorescent-competitive RT-PCR assay is shown; the forward primer is labeled with the 6-FAM fluorophore. The left panel depicts the GeneMapper window, representing an example of the obtained fluorescent products. The peaks shaded in gray correspond to the RT-PCR products; those not shaded represent the size standard (ROX-500 HD). The schematic representation of the corresponding isoform product is displayed on the right. Names (F, full length; ∆N, exon skipped; *, 9-nucleotide-shorter isoform) and accession numbers (if present) of University of California, Santa Cruz Genome Browser-annotated isoforms are indicated. The central panel in (A) represents the levels of IL7R exon 6 skipping (calculated summing the percentage of ∆6 and ∆5∆6 isoforms respect to the total) on RNA extracted from HEK293 cells upon MALAT1 knockdown or overexpression. The percentages were rescaled setting 1 as the value of the mock sample (transfected with a negative control or with an empty vector in knockdown and overexpression experiments, respectively). Bars represent means + SEM of three independent experiments, each performed in triplicate. The boxplots on the right show the percentage of IL7R exon 6-skipping isoforms (A) and SP140 exon 7-skipping isoforms (B) measured in PBMCs of the Italian MS case–control cohort. Boxes define the interquartile range; the thick line refers to the median. The number of subjects belonging to each group is also indicated (N). Significance levels of t-tests are shown. *P < 0.05; **P < 0.01.
As for SP140, this transcript is not endogenously expressed by HEK293 cells; we hence evaluated the level of exon 7-skipping isoforms only in the case–control cohort. The fluorescent-competitive RT-PCR assay, besides the already reported exon 7-skipped isoform (∆7), showed two novel isoforms, both potentially in frame and characterized by the inclusion of exon 6* (a 9-nucleotides-shorter version of exon 6). We confirmed that exon 7-skipping is significantly increased in PBMCs isolated from RR-MS patients (1.7-fold, P = 0.0082; Fig. 4B).
MALAT1 modulation globally impacts on splicing and backsplicing processes
To explore the global effect of MALAT1 modulation on transcriptome, we performed high-coverage RNA-seq experiments of MALAT1 knockdown in Jurkat E6-1 T cells, a more relevant cellular model for MS. Three biological replicates for each condition were processed, obtaining an average of 80.5 million reads per sample and an average of 86.1% of uniquely mapped reads (Supplementary Material, Table S2). The principal component analysis showed the expected clustering of MALAT1 knockdown versus mock samples (Fig. 5A).
MALAT1 knockdown in Jurkat cells: whole-transcriptome differential gene expression analysis. (A) A multidimensional scale plot of all biological replicates of mock and MALAT1 knockdown samples is shown; kd, knockdown. (B) Heat map representing the 107 differentially expressed genes (FDR of <0.1) between mock and MALAT1 knockdown samples. The heatmap was built using the logCPM values of the top dysregulated genes, as calculated by the edgeR package. (C) Validation of RNA-seq differential gene expression results. Expression levels of NUP62CL (nucleoporin 62 C-terminal like) and PPP4R4 (protein phosphatase 4 regulatory subunit 4) were measured by semi-quantitative real-time RT–PCRs, using HMBS as housekeeping gene. Results are presented as normalized rescaled values, setting 1 as the value of the mock sample. Bars represent means + SEM of three independent experiments, each performed in triplicate. **P < 0.01. A table showing the log2FC of the same genes derived from RNA-seq analysis is also shown. (D) Top differentially expressed gene sets (FDR of <0.001) identified by GSEA, using the Reactome database, ordered by normalized enrichment score (NES), based on the gene set enrichment scores for all data set permutations. Positive and negative enrichment scores indicate upregulation and downregulation in MALAT1 knockdown cells, respectively. (E) GSEA results showing selected significantly enriched biological processes and signaling pathways. The green curve represents the enrichment score (vertical axis), calculated by GSEA, showing the measure to which the genes are overrepresented at the top or bottom of a ranked list of genes. The vertical black bars indicate the position in the ranked list of each gene, belonging to the shown gene set. Genes positioned in the red and blue sides are upregulated and downregulated in MALAT1 knockdown cells, respectively. For each pathway, NES and FDR are indicated.
Concerning gene expression analysis, we found 107 differentially expressed genes [false discovery rate (FDR) of <0.1; Supplementary Material, Table S3], of which 77 were downregulated (Fig. 5B). Two among the top dysregulated protein-coding genes (ordered on fold change, FDR of ≤0.05) were chosen to confirm RNA-seq data by semi-quantitative real-time RT-PCR assays; indeed, the NUP62CL and PPP4R4 genes showed the expected reduction to about half levels (Fig. 5C). We then performed a gene set enrichment analysis using GSEA, which identified 100 positively enriched gene sets and 1 negatively enriched one (FDR of <0.01; Supplementary Material, Table S4). Interestingly, the MS-related `respiratory electron transport’, `adaptive immune system’ and `metabolism of mRNA’ gene sets ranked among the top enriched ones, with 79, 443 and 213 genes involved, respectively (Fig. 5D and 5E).
To specifically dissect the role of MALAT1 in AS regulation, we performed a global analysis of AS events using the Multivariate Analysis of Transcript Splicing (rMATS) tool. We found 1114 AS events significantly modulated by MALAT1 (FDR of <0.05), the 62.7% of which is represented by skipped exons (SEs) (Fig. 6A; Supplementary Material, Table S3). Among ASs with at least five reads sustaining the splice junctions and an FDR of <0.01, we selected two skipping events for the validation step: IFNAR2 and EMC4. According to the sequencing data, IFNAR2 exon 8 and EMC4 exon 4 inclusions decreased upon MALAT1 knockdown (Fig. 6B), as also confirmed by semi-quantitative real-time RT-PCR assays (Fig. 6C). Next, we performed a gene set enrichment analysis focusing, among the different types of splicing events, on the SE category, being the most representative one. We found three positively enriched and eight negatively enriched gene sets (FDR of <0.05; Fig. 6D; Supplementary Material, Table S4), including the `metabolism of RNA’ gene set (160 genes) and several pathways related to the post-transcriptional and translational regulation. In addition, we performed a more stringent enrichment analysis, taking into account only the 290 AS events supported by at least 5 reads in 4 out of 6 samples (Fig. 6A; Supplementary Material, Fig. S2). In this case, besides pathways dealing with general biological processes (such as `gene expression’), we evidenced a significant enrichment for the `adaptive immune system’ category (adjusted P < 0.007; Fig. 6E). To identify a possible enrichment in RBP motifs in the exons modulated by MALAT1, we performed an analysis using rMAPS, which calculates by default possible enrichments in 115 different RBPs. Among significant enrichments, we found three RBP motifs as being particularly relevant for our study, i.e. those relative to proteins ANKHD1, SRSF1 and QKI (Fig. 6E). In fact, ANKHD1 emerges among the genes significantly modulated by MALAT1 in our gene expression analysis (Supplementary Material, Table S3), SRSF1 has been already demonstrated to interact with MALAT1 (24) and QKI is known to modulate the AS of myelin protein transcripts as well as to regulate the formation of circular RNAs (circRNAs) (32–34).
MALAT1 knockdown in Jurkat cells: landscape of AS. (A) Summary table of the significant splicing events (FDR of <0.05) identified by rMATS in mock and MALAT1 knockdown samples. The first column represents the different types of splicing events identified by rMATS; the second column (N) indicates the number of significant splicing events; the third column (N after filtering) specifies the number of significant events surviving a more stringent filtering step. In detail, only those events characterized by inclusion/exclusion junctions sustained by at least five counts in four out of six samples were retained. (B) Sashimi plots representing the skipping events detected in IFNAR2 (interferon alpha and beta receptor subunit 2) and EMC4 (ER membrane protein complex subunit 4) genes from RNA-seq data. The inclusion level of the alternative exons is shown in the upper part of the figure where the thickness of the connecting lines is proportional to the number of reads supporting the corresponding splicing event. The lower panel represents the isoforms deriving from the inclusion/skipping of the alternative exons in the mature transcripts. Exon numbering is based on the reference sequence indicated below the scheme. The difference in the inclusion level of the alternative exons (∆), together with the FDR as obtained by rMATS, is also shown. (C) Validation of RNA-seq AS analysis. Inclusion levels of the alternative exons detected in IFNAR2 and EMC4 genes were measured by semi-quantitative real-time RT–PCRs, using HMBS as housekeeping gene. Results are presented as normalized rescaled values, setting 1 as the value of the mock sample. Bars represent means + SEM of three independent experiments, each performed in triplicate. *P < 0.05. (D) Enrichment pathway analysis performed using a ranked list of genes characterized by the presence of at least one skipping event. The top differentially expressed gene sets (FDR of <0.05) identified by GSEA, using the Reactome database, and ordered by NES are listed. Positive and negative enrichment scores indicate upregulation and downregulation in MALAT1 knockdown cells, respectively. (E) Enrichment pathway analysis using all the genes characterized by the presence of at least one skipping event after filtering for read counts (i.e. those enumerated in the `N after filtering’ columns in table A). The top 10 enriched gene sets (FDR of <0.05), identified by Enrichr using the Reactome database, are listed. For each gene set, the corresponding log(FDR) (on the left) as well as the number of genes identified (on the right) are showed. (F) RNA-binding maps based on the results of rMATS analysis on SEs for three proteins whose recognition motives were found significantly enriched in upregulated or downregulated exons or their flanking regions. The red and blue solid lines represent the enrichment in RNA-binding motifs for upregulated and downregulated exons, respectively. The solid black line represents instead the enrichment score for the exons that are not significantly dysregulated (background). Red and blue dashed lines represent the significance of the enrichment for upregulated and downregulated exons, expressed as –log10(P-value). The scheme below each map depicts the alternative exon (in green), whereas numbers indicate the distance in base pairs from it.
Given the global impact of MALAT1 on splicing and the specific modulation of exons containing QKI-binding motifs, we decided to also evaluate the effect of MALAT1 modulation on backsplicing. This process in fact requires the canonical spliceosomal machinery involved in linear splicing and is tuned by the same splicing trans-acting factors (34–36); it leads to the production of circRNAs by joining downstream splice donor sites to upstream acceptor sites (37). Our circRNA analysis, using the DCC software, detected a total of 840 circRNAs (having at least two backspliced reads in three samples; Fig. 7A). Among them, 664 (79%) were already annotated in circBase (38). The identified circRNAs were uniformly distributed on all chromosomes, with chromosome 1 (the largest one) containing, as expected, the higher number of circRNAs (Supplementary Material, Fig. S3A). The expression analysis revealed 49 differentially expressed circRNAs (P < 0.05; Fig. 7B; Supplementary Material, Table S3). Among circRNAs showing a fold change of >2-fold and a P < 0.05, we selected two circRNAs, derived from RNF168 and SLC45A4 genes (both already annotated in circBase as hsa_circ_0123215 and hsa_circ_0001829, respectively) for validation. To this aim, RT-PCR assays were performed using specific divergent primer couples, and the amplification products were verified by direct sequencing (Supplementary Material, Fig. S3B). In addition, we confirmed the changes observed by the RNA-seq analysis by semi-quantitative real-time RT-PCR assays (Fig. 7C). Finally, a gene set enrichment analysis using GSEA evidenced a significant enrichment in the `mRNA splicing/processing’ pathway (nominal P-value of <0.05; Fig. 7D).
MALAT1 knockdown in Jurkat cells: landscape of circRNA expression. (A) Mapping results of linear and circRNA reads on human chromosomes. The outside circle represents the chromosomes, the six inner circles display the counts per million (CPMs), as calculated by edgeR, for each sample. Linear transcripts counts are represented in black and blue for mock and MALAT1-knockdown samples, respectively; those with CPMs higher than 1000 are represented in violet. CircRNAs counts are displayed in red for each sample, and the relative scale is represented at the beginning of each circle. (B) Volcano plot representing circRNAs differentially expressed in MALAT1 knockdown samples. Red and blue dots indicate significantly upregulated and downregulated circRNAs, respectively. The gray shades highlight significantly dysregulated circRNAs that are characterized by a log2FC of ≥1 or a log2FC of ≤−1. (C) Validation of RNA-seq circRNA differential expression results. Expression levels of two circRNAs derived from RNF168 (ring finger protein 168) and SLC45A4 (solute carrier family 45 member 4) genes were measured by semi-quantitative real-time RT–PCRs, using HMBS as housekeeping gene. Results are presented as normalized rescaled values, setting 1 as the value of the mock sample. Bars represent means + SEM of three independent experiments, each performed in triplicate. *P < 0.05, **P < 0.01. On the right, a table showing the log2FC of the same genes as obtained from RNA-seq analysis is also shown. (D) Enrichment pathway analysis performed using a ranked list of circRNAs. The top differentially expressed gene sets (cut-off P-value <0.1), identified by GSEA using the Reactome database and ordered by NES, are listed. Positive and negative enrichment scores indicate upregulation and downregulation in MALAT1 knockdown cells, respectively. Significant P-values (P < 0.05) correspond to NES values above 1.6 or below −1.6.
To confirm the results observed in Jurkat cells, we attempted to perform MALAT1 knockdown in human CD4+ primary T cells. Proof-of-principle experiments were carried out on cells from a healthy donor and resulted in ~60% of MALAT1 silencing level (Supplementary Material, Fig. S4). The measurement of the inclusion level of the IFNAR2 exon 8 and of the expression of the circRNA derived from RNF168 showed the same downregulation trend compared to the results obtained in Jurkat cells (Supplementary Material, Fig. S4).
Discussion
Although several studies showed a differential expression of lncRNAs in autoimmune and neurodegenerative disorders, their role in disease pathogenesis often remained unexplored. Among their functions, lncRNAs are known to affect gene expression by regulating epigenetic, transcriptional and post-transcriptional mechanisms, including the AS process (19,20). Here, we took advantage of publicly available microarray data to identify differentially expressed lncRNAs in MS patients' blood cells. Our analyses disclosed MALAT1 as the most consistently dysregulated lncRNA and led us to focus on the AS process, whose dysregulation has been recently proposed as a novel pathogenic mechanism in MS (4,6).
MALAT1, also known as NEAT2 (nuclear-enriched abundant transcript 2), is an abundant and highly conserved lncRNA in mammals (39,40), originally identified for its overexpression in metastatic early-stage non-small cell lung cancer (41). Its nascent transcript undergoes a post-transcriptional processing, generating a long mature transcript, and a shorter one, the latter characterized by a tRNA-like structure and a cytoplasmic localization (42). The MALAT1 long transcript is instead retained within nuclear speckles, subnuclear domains enriched in splicing factors and involved in the regulation of the pre-mRNA splicing machinery (43). In this respect, MALAT1 has been demonstrated to influence the AS of pre-mRNAs in HeLa and WI-38 cells and to interact and modulate the activity of SR proteins (24,44). Moreover, it has been shown to bind other splicing factors, including several hnRNPs, as well as to influence their expression (44,45). We hence evaluated if MALAT1 could impact on the expression of splicing factors found differentially expressed in MS patients' blood (CELF1, HNRNPH1 and SRSF1) (6) and/or already demonstrated to interact with MALAT1 (HNRNPF, HNRNPH1 and SRSF1) (24,45). The opposite results obtained in SH-SY5Y respect to HEK293 and HeLa cells highlighted a possible cell-specific modulation of splicing factors exerted by MALAT1, which was also confirmed by RNA-seq experiments of MALAT1 knockdown in Jurkat cells. In fact, unlike in the other analyzed cell lines, in Jurkat none of the four-selected splicing factors were differentially expressed, but the `metabolism of mRNA’ gene set still emerged among the top enriched ones, again pointing to a dysregulation of the splicing process.
In addition to several types of cancers (46), MALAT1 has been demonstrated to play a role in autoimmune diseases, such as systemic lupus erythematosus and myositis (47,48). Here, we added MS to the list of autoimmune disorders characterized by a dysregulation of this lncRNA. In particular, we functionally linked MALAT1 dysregulation with the IL7R and SP140 isoform unbalances seen in MS patients' PBMCs.
It remains difficult to establish a general model through which MALAT1 exerts its regulation on AS events. Nonetheless, the well-studied AS regulation of the IL7R transcript can give us some hints. Galarza-Muñoz et al. (49) identified the DDX39B protein as the trans-acting factor regulating IL7R exon 6 AS, also demonstrating its genetic association with MS and a functional epistasis with the IL7R-associated polymorphism. However, their proteomic screen, performed to identify trans-acting factors regulating this AS event (50), inevitably excluded lncRNAs from the analysis. Instead, our RNA-seq data evidenced the involvement of MALAT1 in the regulation of IL7R exon 6; MALAT1 knockdown increased IL7R exon 6 inclusion (FDR = 0.099; Supplementary Material, Fig. S5A), coherently with our data on minigene assays. Interestingly, DDX39B exon 4 was among the 699 skipped events significantly modulated by MALAT1 (FDR = 0.029; Supplementary Material, Fig. S5B), providing a possible link with IL7R exon 6 AS regulation. Indeed, the cascade linking MALAT1 overexpression with the dysregulation of the splicing factor DDX39B—and the consequent effect on IL7R exon 6 skipping—could represent a paradigmatic and not isolated example. One can speculate on the existence of similar regulatory pathways, involving other splicing regulators and impacting on other downstream MS-relevant targets. For instance, our motif analysis evidenced enrichment for the RBPs, QKI and SRSF1 in the exons modulated by MALAT1. Interestingly, QKI is already known to modulate the AS of myelin protein transcripts (32,33), whereas SRSF1 is the major player in the splicing regulation of the MS-associated CD6 transcript (51).
In this already complicated landscape, an additional layer of complexity is given by the observed enrichment for the `RNA metabolism’ gene set in the pathway analysis of transcripts having at least one differentially expressed skipping event upon MALAT1 knockdown (Fig. 6D). Hence, AS represents a mechanism to functionally regulate the RBPs, as documented by several examples in which AS generates isoforms lacking RNA-binding domains, thus potentially affecting the AS of target genes (52–54). Moreover, it is known that splicing factors undergo complex autoregulatory and crossregulatory circuits (55); therefore, even slight variations at the transcriptional or post-transcriptional level in few of them may impact on the activity of other RBPs, and, ultimately, on the AS of many transcripts.
A further mechanism potentially interfering with AS is represented by the backsplicing process, as it competes with the linear splicing mechanism generating circRNAs (36). We hence hypothesized that MALAT1 modulation could also result in a differential expression of circRNAs. Our analysis, highlighting several differentially expressed circRNAs, provides preliminary evidence of a new MALAT1 regulatory function that may also be potentially relevant to MS. Of note, our group and Iparraguirre et al. have already identified dysregulation of circRNAs in MS patient blood (56,57).
Even though MALAT1 represents a well-studied lncRNA, the molecular mechanisms responsible for its altered levels in pathologic settings were rarely elucidated. Several stress conditions were shown to upregulate MALAT1 expression levels, including hypoxia, LPS stimulation, flavivirus infection and H2O2 treatment (58–63). Importantly, the `respiratory electron transport’ gene set was the most significant positively enriched pathway in our differential gene expression analysis in Jurkat cells. Since defects in the mitochondrial electron transport chain complexes I and III are known to increase reactive oxygen species (ROS) production (64,65), we decided to treat Jurkat cells with the H2O2 ROS. We detected a 1.5-fold increase of MALAT1 expression levels (P = 0.016; Supplementary Material, Fig. S6A), perfectly comparable with the change seen in PBMCs of MS patients. ROS production by mitochondrial dysfunction and/or by enzymes expressed in activated microglia and macrophages, such as the NADPH (nicotinamide adenine dinucleotide phosphate) oxidase complex, is a known factor contributing to MS pathogenesis, mostly involving the neurodegeneration and progression processes (66–68). Defects in mitochondrial metabolism have also been found in T cells of MS patients, and the role of ROS in T cell activation, a key process responsible for autoimmune diseases, is becoming increasingly clear (69–71). Importantly, we observed a 1.6-fold increase of IL7R exon 6 skipping upon H2O2 treatment in Jurkat cells (P = 0.032; Supplementary Material, Fig. S6B) providing a mechanism through which an environmental factor could cause splicing isoform unbalances. AS dysregulation, which was found to be exacerbated in specific genes by the presence of MS-associated polymorphisms (12,14,15), may have direct pathogenic consequences, since the inclusion/exclusion of cassette exons in the transcript may potentially produce novel antigenic epitopes (72). Importantly, AS was demonstrated to occur at high rate in known autoantigen transcripts, and the extrathymic expression of alternative isoforms may contribute to the immunogenicity, as demonstrated for PLP1 (proteolipid protein 1), one of MS putative autoantigens. Indeed, only a PLP1 short isoform was shown to be expressed and tolerized by the thymus, while the expression of an alternative exon in the central nervous system was demonstrated to cause the autoimmune response in the experimental–autoimmune–encephalomyelitis mouse model of MS (73).
Of course, the involvement of ROS in inducing increased levels of MALAT1, thus affecting splicing/backsplicing, is just one of the possible mechanistic explanations linking together our experimental evidence. However, we can hypothesize alternative mechanisms, especially considering the many molecular functions that MALAT1 exerts in the nucleus. For instance, it is known that a large percentage of lncRNAs physically interacts with various chromatin regulatory proteins (i.e. readers, writers and erasers), thus activating or repressing gene expression via a chromatin recruitment mechanism (74,75). Among chromatin regulatory proteins, PRC2 (polycomb repressive complex 2) specifically methylates histone 3 at lysine 27, thus inducing transcription repression and influencing splicing regulation (74–76). Interestingly, it has been demonstrated a direct binding of MALAT1 to the PRC2 complex (77), so that it can be speculated a role for MALAT1 similar to that of the HTT transcript in Huntington disease (i.e. extended versions of the polyglutamine tract in HTT RNA enhances huntingtin function as a facilitator of the chromatin repressor PRC2) (74,75). In conclusion, we propose a possible novel pathogenic mechanism for MS in which MALAT1 upregulation modulates the AS of MS-associated genes through a complex regulation of splicing factors. Importantly, the dysregulation of the AS process has been increasingly recognized as a possible pathogenic mechanism underlying also other autoimmune diseases (4,78), so that the proposed mechanism could indeed not only be associated with MS but also with other autoimmune conditions.
Materials and Methods
LncRNA analysis of GEO data sets
Microarray data sets were retrieved from the National Center for Biotechnology Information GEO database according to the following criteria: (i) expression data obtained from blood-derived samples, (ii) case–control design and (iii) data available for at least 10 RR-MS cases and 10 healthy controls. Seven data sets (GSE21942, GSE41848, GSE41849, GSE41890, GSE17048, GSE43592 and GSE13732) matched the inclusion criteria; their characteristics are shown in Supplementary Material, Table S1. To identify differentially expressed genes between cases and controls, the GEO2R web application (79,80) was used, considering as significant a Benjamini–Hochberg-adjusted P-value of 0.05. The list of genes obtained for each data set was crossed with the list of lncRNAs obtained from GENCODE (using the comprehensive gene annotation of lncRNA genes on the reference chromosomes, version 25) (81). Finally, the lists of the lncRNAs identified in each data set—and surviving multiple testing—were compared to identify common genes.
The meta-analysis was performed by the NetworkAnalyst tool (31) using the data sets compatible with the software, i.e. GSE13732 (excluding the cases/controls resampled at a second-time point), GSE21942 and GSE43592. Data sets were uploaded, processed, annotated and checked for their integrity. The meta-analysis was performed using the `combining P-values’ option, based on the Fisher's method [−2*∑Log(P)], setting the threshold of significance at 0.05 (31). MALAT1 mean fold change was calculated using the logFC available for each data set.
Subjects
This study was approved by local ethical committees and conducted according to the Helsinki Declaration. All subjects signed an informed consent. All patients were affected by clinically definite MS according to the revised McDonald's criteria (82).
All patients of the Italian cohort (mean age, 48.0 ± 6.3; 87% females) were diagnosed with the RR-MS form and, at the time of blood withdrawal, were in the remission phase and had not received any immunomodulatory treatment at least for a month. The enrolled 30 healthy controls (mean age, 48.3 ± 16.0; 57% females) declared no familial history for autoimmune or neurodegenerative diseases.
The US population of patients included 29 RR-MS (being in the remission phase and free from treatment for at least a month at the time of blood collection) and 7 PP-MS patients. Their mean age was 53.9 ± 13.8, 67% were females and 92% were of Caucasian origin. The 20 healthy control subjects (mean age, 53.4 ± 15.7; 60% females and 85% Caucasians) declared no familial history for autoimmune or neurodegenerative diseases.
PBMCs were isolated by centrifugation on a Lympholyte Cell Separation Medium (Cedarlane Laboratories Ltd, Hornby, Ontario, Canada) gradient. Total RNA was isolated using the EuroGold Trifast kit (Euroclone, Wetherby, UK) for the Italian cohort and using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) for the US one. RNA concentrations were measured using the NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA).
Plasmid constructs
The plasmid expressing MALAT1 was a kind gift of Professor K.V. Prasanth (University of Illinois, Urbana, Illinois). MALAT1 most abundant isoform (whose transcription start site is 1283 nucleotides downstream that of the long form, GenBank accession number NR_002819.2) was cloned into the pEYFP-N1 vector (Clontech, Mountain View, CA) by replacing the EYFP cDNA with MALAT1, through the AgeI and NotI restriction enzymes.
For splicing assays, regions containing IL7R exon 6 or SP140 (SP140 nuclear body protein) exon 7, with their respective intronic flanking sequences, were amplified from the genomic DNA of a healthy control heterozygous for the MS-associated polymorphisms rs6897932 and rs28445040. Both allelic versions of each exon were cloned into the hybrid α-globin-fibronectin minigene plasmid (modified pBS-KS) (83), as previously described (84).
Cell cultures and transfections
HEK293 and HeLa cells were cultured in Dulbecco Modified Eagle's Medium (EuroClone), SH-SY5Y and Jurkat E6-1 cells in RPMI 1640 (EuroClone). In both media, 10% fetal bovine serum, 1% glutamine and antibiotics (100 U/ml penicillin and 100 μg/ml streptomycin; EuroClone) were added. Medium for Jurkat E6-1 cells was also added with 1% sodium pyruvate (1 mm; Sigma-Aldrich, Saint Louis, MO) and 1% non-essential amino acids. Cells were grown at 37°C in a humidified atmosphere of 5% CO2 and 95% air, according to standard procedures.
For MALAT1 overexpression, HEK293, SH-SY5Y and HeLa cells were transfected using 1 μg of the MALAT1-expressing plasmid or 1 μg of an empty vector. In each experiment, an equal number of cells (2.5 × 105 for HeLa; 3 × 105 for HEK293 and SH-SY5Y) were transfected with the Polyplus jetPRIME (EuroClone) in 6-well plates, as described by the manufacturer. Cells were collected to perform total RNA extraction using the EuroGold Trifast kit (Euroclone) 24 h after transfection.
For MALAT1 knockdown, a specific locked nucleic acid (LNA™) longRNA GapmeR (Exiqon, Vedbaek, Denmark) targeting MALAT1 was designed using a software available on the Exiqon website (sequence, 5′-GACAAGATTCATGAGT-3′). Exiqon GapmeR negative control A was used as negative control (sequence, 5′-AACACGTCTATACGC-3′). HEK293, SH-SY5Y and HeLa cells were transfected using 50 nm GapmeR targeting MALAT1 or 50 nm GapmeR negative control. To achieve a greater knockdown efficiency, for HEK293 and HeLa cells, the transfection was repeated after 24 h, using the same concentration of antisense oligonucleotides.
For splicing minigene assays, HEK293 were co-transfected as described above using (i) 1 μg of the IL7R/SP140 minigene and MALAT1-expressing plasmid in a 1:2 ratio or (ii) 1 μg of the IL7R/SP140 minigene and an empty vector in a 1:2 ratio.
For MALAT1 knockdown in Jurkat E6-1 cells, the transfection was achieved through a gymnotic delivery, as detailed by Fazil et al. (85). Briefly, 5 × 104 cells were incubated in a 24-well plate with 500 nm GapmeR targeting MALAT1 or 500 nm GapmeR negative control. After 24 h, cells were collected for total RNA extraction.
Semi-quantitative real-time RT-PCR and fluorescent-competitive RT-PCR
Random hexamers (Promega, Madison, WI) and the Superscript-III Reverse Transcriptase (Thermo Fisher Scientific) were used to perform first-strand cDNA synthesis, according to the manufacturer's instructions. A total of 1 μl (undiluted or properly diluted) of the RT reaction was used as template for amplifications.
Semi-quantitative real-time RT-PCRs were accomplished by using the SYBR Premix Ex Taq II (TaKaRa, Kusatsu, Japan) and a touchdown thermal protocol on a LightCycler 480 (Roche, Basel, Switzerland). HMBS (hydroxymethylbilane synthase), ACTB (beta actin) or GAPDH (glyceraldehyde 3-phosphate dehydrogenase) were used as housekeeping genes. Reactions were performed at least in triplicate, and expression data were analyzed using the geNorm software (86).
To quantitate the percentage of isoforms skipping IL7R exon 6 and SP140 exon 7, both for the minigene and endogenous assays, fluorescent-competitive RT-PCR assays were performed using, in each reaction, a 6-FAM-labeled primer; reactions were performed under standard conditions using the GoTaq DNA Polymerase (Promega) on a Mastercycler EP Gradient (Eppendorf, Hamburg, Germany). Amplified fragments were separated by capillary electrophoresis on an ABI 3500 Genetic Analyzer (Thermo Fisher Scientific) and the peak areas measured by the GeneMapper v4.0 software (Thermo Fisher Scientific). The percentage of each specific isoform was measured by calculating the ratio of the relevant peak over the sum of all the fluorescence peak areas (set as 100%).
The sequences of all primer couples will be provided upon request.
All t-test were performed using the R software (http://www.r-project.org/). P-values <0.05 were considered as statistically significant.
RNA sequencing
Total RNA was extracted from mock and MALAT1 knockdown Jurkat E6.1 samples using the Maxwell 16 LEV simplyRNA Cells Kit (Promega). RNA quality (RNA quality score of >8) was assessed by LabChip GX Touch (PerkinElmer, Waltham, MA). Libraries were prepared starting from 500 ng of total RNA, using the TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA) and following the manufacturer's instructions. Samples underwent a high-coverage paired-end 150 bp strand-specific sequencing using a NextSeq 500 platform (Illumina).
RNA-seq data analysis
For differential gene expression analysis, sequencing reads were mapped to the human genome (hg19) using STAR (version 2.5.2) (87). Transcript quantification from mapped reads was performed using HTSeq-count (version 0.6.1p1) (88) and the human transcripts annotations from Ensembl database (GRCh37 version) (89). Differential expression analysis was carried out using R and the edgeR package (90). An FDR of <0.1 was set as threshold.
For determining the splicing-isoform landscape, raw reads of each sample were first filtered, trimmed and properly paired again to retain only 149 bp long reads, using the NGSUtils suite (91). The obtained reads were mapped to hg19 using STAR, and the mapping results (in bam format) were submitted to rMATS (version 3.2.5) to perform the AS analysis (92) and identify five basic types of AS patterns: SEs, mutually exclusive exons, alternative 5′ and 3′ splice sites and retained introns. An FDR of <0.05 was set as threshold. Only the reads mapping to exon–exon junctions were used for all the analyses.
Subsequently, the output of the SE analysis was submitted to the rMAPS web server (93), to detect a possible enrichment of specific RBP motifs in the pools of upregulated or downregulated exons.
For circRNA expression analysis, we used the DCC software (94). In details, the read pairs from paired-end data were aligned both together and separately using STAR and the options suggested by the software pipeline. Filtering steps were the following: (i) only circRNAs with at least 2 reads in 3 samples were retained and (ii) circRNAs mapping on mitochondrial DNA or repeat regions were discarded. Once the list of circRNAs was obtained, the statistical significance was assessed using edgeR. More specifically, backsplicing read counts were added to the linear gene count list as new genes for the normalization purpose, as described by Dou et al. (95), and a standard differentially expression analysis was carried out. A difference in gene expression was considered significant if the unadjusted P-value was <0.05. For validation, semi-quantitative real-time RT-PCR assays were performed as described above, using specific divergent primer couples.
Gene set enrichment analysis
Enrichment analysis was performed using GSEA (v.3.0) (96). The canonical pathways Reactome database, containing 674 gene sets, was tested for enrichment on RNA-seq results. Rank scores for differential expression data were calculated as −log10(P-value) multiplied by the sign of the edgeR fold change for the gene expression results (97) and as log10(P-value) multiplied by the sign of the inclusion level difference for the SE analysis. The GSEAPreranked function was applied using 1000 permutations and the classic option of enrichment statistic.
Data availability
The RNA-seq data from this publication have been deposited to GEO database (accession number: GSE110525).
Conflict of Interest statement. None declared.
Funding
Harry Weaver Neuroscience Scholar Award of the National Multiple Sclerosis Society (JF 2144A2/1 to L.P.); Fondazione Italiana Sclerosi Multipla (FISM) fellowship (2012/B/1 to C.C.); NMSS fellowship (FG 2010-A1/2 to C.C.).
