Advertisement

The circular RNA landscape in multiple sclerosis: Disease-specific associated variants and exon methylation shape circular RNA expression profile

Open AccessPublished:November 22, 2022DOI:https://doi.org/10.1016/j.msard.2022.104426

      Highlights

      • Circular RNA expression in multiple sclerosis patients is dysregulated.
      • Multiple-sclerosis-associated variants can influence circular RNA expression.
      • DNA methylation may play a role in defining circular RNA landscape.

      Abstract

      Background

      Circular RNAs (circRNAs) are a class of non-coding RNAs increasingly emerging as crucial actors in the pathogenesis of human diseases, including autoimmune and neurological disorders as multiple sclerosis (MS). Despite several efforts, the mechanisms regulating circRNAs expression are still largely unknown and the circRNA profile and regulation in MS-relevant cell models has not been completely investigated. In this work, we aimed at exploring the global landscape of circRNA expression in MS patients, also evaluating a possible correlation with their genetic and epigenetic background.

      Methods

      We performed RNA-seq experiments on circRNA-enriched samples, derived from peripheral blood mononuclear cells (PBMCs) of 10 MS patients and 10 matched controls and performed differential circRNA expression. The genetic background was evaluated using array genotyping, and an expression quantitative trait loci (eQTL) analysis was carried out.

      Results

      Expression analysis revealed 166 differentially expressed circRNAs in MS patients, 125 of which are downregulated. One of the top dysregulated circRNAs, hsa_circ_0007990, derives from the PGAP3 gene, encoding a protein relevant for the control of autoimmune responses. The downregulation of this circRNA was confirmed in two independent replication cohorts, suggesting its implementation as a possible RNA-based biomarker. The eQTL analysis evidenced a significant association between 89 MS-associated loci and the expression of at least one circRNA, suggesting that MS-associated variants could impact on disease pathogenesis by altering circRNA profiles. Finally, we found a significant correlation between exon methylation and circRNA expression levels, supporting the hypothesis that epigenetic features may play an important role in the definition of the cell circRNA pool.

      Conclusion

      We described the circRNA expression profile of PBMCs in MS patients, suggesting that MS-associated variants may tune the expression levels of circRNAs acting as “circ-QTLs”, and proposing a role for exon-based DNA methylation in regulating circRNA expression.

      Keywords

      1. Introduction

      Circular RNAs (circRNAs) are a new class of non-coding RNAs (ncRNAs) that are characterized by a covalently closed structure (
      • Lasda E.
      • Parker R.
      Circular RNAs: diversity of form and function.
      ). Thousands of circRNAs have been detected in the human transcriptome, especially thanks to the advancement in bioinformatic analyses of transcriptomic data (
      • Salzman J.
      • Gawad C.
      • Wang P.L.
      • Lacayo N.
      • Brown P.O.
      Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types.
      ;
      • Memczak S.
      • Jens M.
      • Elefsinioti A.
      • Torti F.
      • Krueger J.
      • Rybak A.
      • Maier L.
      • Mackowiak S.D.
      • Gregersen L.H.
      • Munschauer M.
      • Loewer A.
      • Ziebold U.
      • Landthaler M.
      • Kocks C.
      • le Noble F.
      • Rajewsky N.
      Circular RNAs are a large class of animal RNAs with regulatory potency.
      ;
      • Jeck W.R.
      • Sorrentino J.A.
      • Wang K.
      • Slevin M.K.
      • Burd C.E.
      • Liu J.
      • Marzluff W.F.
      • Sharpless N.E.
      Circular RNAs are abundant, conserved, and associated with ALU repeats.
      ). Several types of circRNAs have been described, such as exonic, intronic, exonic-intronic, and intergenic circRNAs (
      • Ebbesen K.K.
      • Kjems J.
      • Hansen T.B.
      Circular RNAs: identification, biogenesis and function.
      ). Exonic circRNAs represent the most common class, and originate from the back-splicing (BS) process, in which a downstream splice-donor site is ligated to an upstream splice-acceptor site (
      • Lasda E.
      • Parker R.
      Circular RNAs: diversity of form and function.
      ). This process is promoted by the presence of inverted repetitive sequences, by specific binding sites for splicing factors in the intronic flanking regions, or by the formation of a lariat deriving from exon skipping (
      • Jeck W.R.
      • Sorrentino J.A.
      • Wang K.
      • Slevin M.K.
      • Burd C.E.
      • Liu J.
      • Marzluff W.F.
      • Sharpless N.E.
      Circular RNAs are abundant, conserved, and associated with ALU repeats.
      ;
      • Ashwal-Fluss R.
      • Meyer M.
      • Pamudurti N.R.
      • Ivanov A.
      • Bartok O.
      • Hanan M.
      • Evantal N.
      • Memczak S.
      • Rajewsky N.
      • Kadener S.
      circRNA biogenesis competes with pre-mRNA splicing.
      ). Among their proposed functions, it has been demonstrated that circRNAs can act as molecular sponges for microRNAs (miRNAs), thus affecting the expression of their targets (
      • Memczak S.
      • Jens M.
      • Elefsinioti A.
      • Torti F.
      • Krueger J.
      • Rybak A.
      • Maier L.
      • Mackowiak S.D.
      • Gregersen L.H.
      • Munschauer M.
      • Loewer A.
      • Ziebold U.
      • Landthaler M.
      • Kocks C.
      • le Noble F.
      • Rajewsky N.
      Circular RNAs are a large class of animal RNAs with regulatory potency.
      ;
      • Hansen T.B.
      • Jensen T.I.
      • Clausen B.H.
      • Bramsen J.B.
      • Finsen B.
      • Damgaard C.K.
      • Kjems J.
      Natural RNA circles function as efficient microRNA sponges.
      ). Few circRNAs have been demonstrated to be translated through internal ribosome entry sites that allow a cap-independent translation (
      • Legnini I.
      • Di Timoteo G.
      • Rossi F.
      • Morlando M.
      • Briganti F.
      • Sthandier O.
      • Fatica A.
      • Santini T.
      • Andronache A.
      • Wade M.
      • Laneve P.
      • Rajewsky N.
      • Bozzoni I.
      Circ-ZNF609 is a circular RNA that can be translated and functions in myogenesis.
      ;
      • Pamudurti N.R.
      • Bartok O.
      • Jens M.
      • Ashwal-Fluss R.
      • Stottmeister C.
      • Ruhe L.
      • Hanan M.
      • Wyler E.
      • Perez-Hernandez D.
      • Ramberger E.
      • Shenzis S.
      • Samson M.
      • Dittmar G.
      • Landthaler M.
      • Chekulaeva M.
      • Rajewsky N.
      • Kadener S.
      Translation of CircRNAs.
      ). Moreover, the BS process directly competes with pre-mRNA splicing (
      • Ashwal-Fluss R.
      • Meyer M.
      • Pamudurti N.R.
      • Ivanov A.
      • Bartok O.
      • Hanan M.
      • Evantal N.
      • Memczak S.
      • Rajewsky N.
      • Kadener S.
      circRNA biogenesis competes with pre-mRNA splicing.
      ;
      • Conn S.J.
      • Pillman K.A.
      • Toubia J.
      • Conn V.M.
      • Salmanidis M.
      • Phillips C.A.
      • Roslan S.
      • Schreiber A.W.
      • Gregory P.A.
      • Goodall G.J.
      The RNA binding protein quaking regulates formation of circRNAs.
      ); therefore, the biogenesis of circRNAs interferes with the alternative splicing (AS) process, potentially inducing exon skipping from the remaining linear RNA or its degradation (
      • Chen L.L.
      • Yang L.
      Regulation of circRNA biogenesis.
      ). The tight interconnection between BS and AS is also supported by the following observations: (i) BS requires the canonical spliceosomal machinery involved in linear splicing (
      • Starke S.
      • Jost I.
      • Rossbach O.
      • Schneider T.
      • Schreiner S.
      • Hung L.H.
      • Bindereif A.
      Exon circularization requires canonical splice signals.
      ); (ii) BS is tuned by the same trans-acting splicing factors (
      • Ashwal-Fluss R.
      • Meyer M.
      • Pamudurti N.R.
      • Ivanov A.
      • Bartok O.
      • Hanan M.
      • Evantal N.
      • Memczak S.
      • Rajewsky N.
      • Kadener S.
      circRNA biogenesis competes with pre-mRNA splicing.
      ;
      • Conn S.J.
      • Pillman K.A.
      • Toubia J.
      • Conn V.M.
      • Salmanidis M.
      • Phillips C.A.
      • Roslan S.
      • Schreiber A.W.
      • Gregory P.A.
      • Goodall G.J.
      The RNA binding protein quaking regulates formation of circRNAs.
      ;
      • Kramer M.C.
      • Liang D.
      • Tatomer D.C.
      • Gold B.
      • March Z.M.
      • Cherry S.
      • Wilusz J.E.
      Combinatorial control of Drosophila circular RNA expression by intronic repeats, hnRNPs, and SR proteins.
      ); (iii) circRNAs have also been shown to bind and affect the function of splicing factors (
      • Ashwal-Fluss R.
      • Meyer M.
      • Pamudurti N.R.
      • Ivanov A.
      • Bartok O.
      • Hanan M.
      • Evantal N.
      • Memczak S.
      • Rajewsky N.
      • Kadener S.
      circRNA biogenesis competes with pre-mRNA splicing.
      ;
      • Abdelmohsen K.
      • Panda A.C.
      • Munk R.
      • Grammatikakis I.
      • Dudekula D.B.
      • De S.
      • Gorospe M.
      Identification of HuR target circular RNAs uncovers suppression of PABPN1 translation by CircPABPN1.
      ).
      Considering that circRNAs affect gene expression regulation at many levels, it is not surprising that they are involved in the pathogenesis and progression of human diseases, including cancer, cardiovascular, neurological, and autoimmune disorders (
      • Han B.
      • Chao J.
      • Yao H.
      Circular RNA and its mechanisms in disease: from the bench to the clinic.
      ). Due to their resistance to exonuclease-mediated degradation (
      • Cocquerelle C.
      • Mascrez B.
      • Hétuin D.
      • Bailleul B.
      Mis-splicing yields circular RNA molecules.
      ), circRNAs are highly stable in biological fluids, thus representing also potential biomarkers for several disorders (
      • Maass P.G.
      • Glažar P.
      • Memczak S.
      • Dittmar G.
      • Hollfinger I.
      • Schreyer L.
      • Sauer A.V
      • Toka O.
      • Aiuti A.
      • Luft F.C.
      • Rajewsky N.
      A map of human circular RNAs in clinically relevant tissues.
      ).
      Multiple sclerosis (MS) is a complex autoimmune disease of the central nervous system characterized by demyelination, chronic inflammation, neuronal loss, and axonal damage (
      • Thompson A.J.
      • Baranzini S.E.
      • Geurts J.
      • Hemmer B.
      • Ciccarelli O.
      Multiple sclerosis.
      ). MS is classified in different clinical courses, depending on the severity and progression of symptoms; the most common form, involving 80% of patients, is the relapsing-remitting (RR) one, characterized by attacks followed by complete or partial remissions (
      • Lublin F.D.
      • Reingold S.C.
      • Cohen J.A.
      • Cutter G.R.
      • Sørensen P.S.
      • Thompson A.J.
      • Wolinsky J.S.
      • Balcer L.J.
      • Banwell B.
      • Barkhof F.
      • Bebo Jr B.
      • Calabresi P.A.
      • Clanet M.
      • Comi G.
      • Fox R.J.
      • Freedman M.S.
      • Goodman A.D.
      • Inglese M.
      • Kappos L.
      • Kieseier B.C.
      • Lincoln J.A.
      • Lubetzki C.
      • Miller A.E.
      • Montalban X.
      • O'Connor P.W.
      • Petkau J.
      • Pozzilli C.
      • Rudick R.A.
      • Sormani M.P.
      • Stüve O.
      • Waubant E.
      • Polman C.H.
      Defining the clinical course of multiple sclerosis: the 2013 revisions.
      ). Although the pathogenic factors underlying MS remain largely unknown, in the last years several studies pointed to alterations in AS and RNA processing as new molecular mechanisms potentially involved in the disease (
      • Paraboschi E.M
      • Cardamone G.
      • Rimoldi V.
      • Gemmati D.
      • Spreafico M.
      • Duga S.
      • Soldà G.
      • Asselta R.
      Meta-analysis of multiple sclerosis microarray data reveals dysregulation in RNA splicing regulatory genes.
      ;
      • Spurlock 3rd C.F.
      • Tossberg J.T.
      • Guo Y.
      • Sriram S.
      • Crooke 3rd P.S.
      • Aune T.M.
      Defective structural RNA processing in relapsing-remitting multiple sclerosis.
      ;
      • Hecker M.
      • Rüge A.
      • Putscher E.
      • Boxberger N.
      • Rommer P.S.
      • Fitzner B.
      • Zettl U.K.
      Aberrant expression of alternative splicing variants in multiple sclerosis – a systematic review.
      ). Moreover, several AS events as well as non-coding RNAs have been proposed as novel RNA-based biomarkers for the disease (
      • Hecker M.
      • Rüge A.
      • Putscher E.
      • Boxberger N.
      • Rommer P.S.
      • Fitzner B.
      • Zettl U.K.
      Aberrant expression of alternative splicing variants in multiple sclerosis – a systematic review.
      ;
      • D'Ambrosio A.
      • Pontecorvo S.
      • Colasanti T.
      • Zamboni S.
      • Francia A.
      • Margutti P.
      Peripheral blood biomarkers in multiple sclerosis.
      ).
      Concerning circRNAs, our group identified, for the first time, one dysregulated circRNA in MS patients, the GSDMB hsa_circ_0106803 circRNA. It was shown to be upregulated in peripheral blood mononuclear cells (PBMCs) of 30 RR-MS patients compared to 30 healthy controls (
      • Cardamone G.
      • Paraboschi E.M.
      • Rimoldi V.
      • Duga S.
      • Soldà G.
      • Asselta R.
      The characterization of GSDMB splicing and backsplicing profiles identifies novel isoforms and a circular RNA that are dysregulated in multiple sclerosis.
      ). Subsequently, Iparraguirre and colleagues detected 406 differentially expressed circRNAs through a microarray analysis on peripheral blood leucocytes of four RR-MS patients and four healthy controls; two downregulated circRNAs, both deriving from the ANXA2 host gene, were validated in an independent cohort of 20 RR-MS patients and 18 healthy controls (
      • Iparraguirre L.
      • Muñoz-Culla M.
      • Prada-Luengo I.
      • Castillo-Triviño T.
      • Olascoaga J.
      • Otaegui D.
      Circular RNA profiling reveals that circular RNAs from ANXA2 can be used as new biomarkers for multiple sclerosis.
      ). Our group also demonstrated an enrichment of circRNAs at MS genome wide-associated loci, as well as a genotype-dependent regulation of a circRNA derived from the MS-associated STAT3 gene (
      • Paraboschi E.M.
      • Cardamone G.
      • Soldà G.
      • Duga S.
      • Asselta R.
      Interpreting non-coding genetic variation in multiple sclerosis genome-wide associated regions.
      ). In addition to genetic variants, further mechanisms have been recently proposed to control circRNA expression levels, including epigenetic factors and regulatory networks involving long non-coding RNAs (lncRNAs) (
      • Holdt L.M.
      • Stahringer A.
      • Sass K.
      • Pichler G.
      • Kulak N.A.
      • Wilfert W.
      • Kohlmaier A.
      • Herbst A.
      • Northoff B.H.
      • Nicolaou A.
      • Gäbel G.
      • Beutner F.
      • Scholz M.
      • Thiery J.
      • Musunuru K.
      • Krohn K.
      • Mann M.
      • Teupser D.
      Circular non-coding RNA ANRIL modulates ribosomal RNA maturation and atherosclerosis in humans.
      ;
      • Ferreira H.J.
      • Davalos V.
      • de Moura M.C.
      • Soler M.
      • Perez-Salvia M.
      • Bueno-Costa A.
      • Setien F.
      • Moran S.
      • Villanueva A.
      • Esteller M.
      Circular RNA CpG island hypermethylation-associated silencing in human cancer.
      ;
      • Kleaveland B.
      • Shi C.Y.
      • Stefano J.
      • Bartel D.P.
      A network of noncoding regulatory RNAs acts in the mammalian brain.
      ). Indeed, we recently demonstrated that the lncRNA MALAT1, upregulated in MS patients’ blood, was able to modulate the expression levels of several circRNAs (
      • Cardamone G.
      • Paraboschi E.M.
      • Soldà G.
      • Cantoni C.
      • Supino D.
      • Piccio L.
      • Duga S.
      • Asselta R.
      Not only cancer: the long non-coding RNA MALAT1 affects the repertoire of alternatively spliced transcripts and circular RNAs in multiple sclerosis.
      ). Despite this evidence suggesting a potential involvement of circRNAs in MS, the complete landscape of circRNA expression in MS patients is still missing.
      In this work, we aimed at exploring the global landscape of circRNA expression in MS patients, also evaluating a possible correlation with their genetic and epigenetic background.

      2. Materials and methods

      2.1 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 RR-MS and, at the time of blood collection, were in the remission phase. All the healthy control subjects declared no familial history for autoimmune or neurodegenerative diseases. The main characteristics of the subjects included in the discovery and the replication cohorts are listed in the Supplementary Table 1.

      2.2 DNA and RNA samples

      DNA samples were extracted from peripheral blood using the Maxwell® RSC Blood DNA Kit (Promega, Madison, WI, USA).
      PBMCs were isolated by centrifugation on a Lympholyte Cell separation medium (Cedarlane, Burlington, Ontario, Canada) gradient and total RNA was extracted using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) for the US cohorts and the EuroGold Trifast kit (Euroclone, Wetherby, UK) for the Italian one.
      PGAP3 hsa_circ_0007990 circRNA expression levels were measured by using RNAs both from a commercial panel of human tissues (Invitrogen, Carlsbad, California, USA) (for each tissue, RNAs are pooled from at least 3 donors), and from sorted Th1 and Th17 cells, isolated from pooled buffy coats of healthy donors by using a FACSAria Cell Sorter (BD Biosciences, San Jose, USA).
      DNA and RNA concentrations were measured using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The integrity of RNA samples (RNA quality score >8) was assessed by LabChip GX Touch (Perkin Elmer, Waltham, MA, USA).

      2.3 RNA sequencing

      One µg of RNA was treated with 3 units of RNase R (Epicentre, Madison, WI, USA) for 15 min at 37 °C. After phenol-chloroform purification, the entire elution volume of each RNA was used for the library preparation, using the TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) and following the manufacturer's instructions for degraded RNA samples. Samples were sequenced using the NextSeq 500 platform (Illumina) to obtain 150-nt-long paired-end reads.

      2.4 RNA-seq data analysis

      CircRNAs were detected using the DCC software (
      • Cheng J.
      • Metge F.
      • Dieterich C.
      Specific identification and quantification of circular RNAs from sequencing data.
      ). The read pairs from paired-end data were aligned both together and separately to the human genome (hg19) using STAR (version 2.5.2) (
      • Dobin A.
      • Davis C.A.
      • Schlesinger F.
      • Drenkow J.
      • Zaleski C.
      • Jha S.
      • Batut P.
      • Chaisson M.
      • Gingeras T.R.
      STAR: ultrafast universal RNA-seq aligner.
      ) and the options suggested by the software pipeline. Only circRNAs with at least 5 reads in 10 samples were retained; circRNAs mapping on mitochondrial DNA or repeat regions were discarded. Once the list of circRNAs was obtained, the statistical significance was assessed using the edgeR package (
      • Robinson M.D.
      • McCarthy D.J.
      • Smyth G.K.
      edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
      ). Back-splicing read counts were added to the linear gene count list as new genes for the normalization purpose (
      • Dou Y.
      • Cha D.J.
      • Franklin J.L.
      • Higginbotham J.N.
      • Jeppesen D.K.
      • Weaver A.M.
      • Prasad N.
      • Levy S.
      • Coffey R.J.
      • Patton J.G.
      • Zhang B.
      Circular RNAs are down-regulated in KRAS mutant colon cancer cells and can be transferred to exosomes.
      ), and a standard differential expression analysis was carried out. A difference in gene expression was considered significant if the unadjusted P value was <0.05.

      2.5 Gene ontology (GO) enrichment analysis

      Enrichment analysis for the host genes of the 166 differentially expressed circRNAs was performed using the “topGO” Bioconductor tool, the “biological processes”, the “molecular functions” database, and the elim algorithm. Results are presented as row p-values, because the parent-child relationships of GO terms are not truly independent, and the p-value of a GO term is therefore conditioned on the neighboring terms (

      A. Alexa, J. Rahnenfuhrer (2016), topGO: enrichment analysis for gene ontology. R package version 2.24.0.

      ).

      2.6 Semi-quantitative real-time RT-PCR

      Random hexamers (Promega) and the Superscript-IV Reverse Transcriptase (Invitrogen) were used to perform first-strand cDNA synthesis, according to the manufacturer's instructions. One µl of the RT reaction was used as template for the subsequent semi-quantitative real-time RT-PCR assays.
      Semi-quantitative real-time RT-PCRs were accomplished by using the LightCycler® 480 SYBR Green I Master mix (Roche, Basel, Switzerland) and a touchdown thermal protocol on a LightCycler 480 (Roche). For circRNA expression analysis, specific divergent primer couples were designed and HMBS (hydroxymethylbilane synthase) was used as housekeeping gene for normalization. Reactions were performed at least in triplicate, and expression data were analyzed using the 2^(-delta delta Ct) method. Normalized values were then rescaled using as a reference the sample characterized by the lowest expression of that specific gene/circRNA. The sequences of all primer couples will be provided upon request.
      Normality distribution of data was checked using Shapiro-Wilk test, and t-test or Mann-Whitney U test were accordingly performed using the R software (http://www.r-project.org/) (
      R Core Team
      R: A Language and Environment for Statistical Computing.
      ). P values <0.05 were considered as statistically significant.

      2.7 Array genotyping and analysis

      Genotyping was performed using the Infinium® HumanCore-24 v1.1 BeadChip (Illumina), containing 307,342 markers spread all over the genome. BeadChips were scanned using the iScan platform (Illumina); genotype calling was performed with GenomeStudioTM V2.0 software (Illumina).
      Imputation was performed remotely through the Michigan Imputation Server (https://imputationserver.sph.umich.edu), using the Human Reference Consortium (HRC) data r1.1 2016 as reference panel, the Eagle v.2.3 program for the phasing step, and Minimac3 as imputation software (
      • Das S.
      • Forer L.
      • Schönherr S.
      • Sidore C.
      • Locke A.E.
      • Kwong A.
      • Vrieze S.I.
      • Chew E.Y.
      • Levy S.
      • McGue M.
      • Schlessinger D.
      • Stambolian D.
      • Loh P.R.
      • Iacono W.G.
      • Swaroop A.
      • Scott L.J.
      • Cucca F.
      • Kronenberg F.
      • Boehnke M.
      • Abecasis G.R.
      • Fuchsberger C.
      Next-generation genotype imputation service and methods.
      ;
      • Loh P.R.
      • Danecek P.
      • Palamara P.F.
      • Fuchsberger C.
      • Reshef Y.A.
      • Finucane H.K.
      • Schoenherr S.
      • Forer L.
      • McCarthy S.
      • Abecasis G.R.
      • Durbin R.
      • Price A.L.
      Reference-based phasing using the haplotype reference consortium panel.
      ). The imputed dataset was then filtered to retain only those variants with a good imputation quality (r2>0.3). X chromosome was omitted from this analysis.
      For the eQTL analysis, MS-associated SNPs, excluding those mapping in the highly complex HLA region, were retrieved from the literature (
      International Multiple Sclerosis Genetics Consortium (IMSGC)
      Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility.
      ). The analysis was performed using the FastQTL tool (
      • Ongen H.
      • Buil A.
      • Brown A.A.
      • Dermitzakis E.T.
      • Delaneau O.
      Fast and efficient QTL mapper for thousands of molecular phenotypes.
      ), a mapping window of a 1 megabase, and the adaptive permutations mode.

      2.8 Sanger sequencing

      Direct sequencing was performed to genotype the rs7295386 polymorphism in the replication cohorts, and to confirm the structure of the circRNAs identified by bioinformatic analyses and tested in further experiments. For genotyping experiments, the genomic region containing the variant was amplified by standard PCR reactions performed on 20 ng of genomic DNA by using the GoTaq DNA polymerase (Promega). For circRNAs, the product of the amplification derived by divergent primer couples was verified for the presence of the back-spliced junction. The sequencing reactions were prepared with the BigDye Terminator Cycle Sequencing Ready Reaction Kit v1.1 (Thermo Fisher Scientific) and run on an ABI-3500 DX Genetic Analyzer (Thermo Fisher Scientific). Primer sequences will be provided upon request.

      2.9 DNA methylation analyses

      The correlation between methylation level and circRNA expression was evaluated in Jurkat cells. In detail, methylation data were obtained from UCSC Genome Browser website (https://genome.ucsc.edu), whereas the circRNA profile in this same cell type was derived from a previous work (
      • Paraboschi E.M.
      • Cardamone G.
      • Soldà G.
      • Duga S.
      • Asselta R.
      Interpreting non-coding genetic variation in multiple sclerosis genome-wide associated regions.
      ). The following method was applied:
      • 1
        only the exons with available methylation data were selected (methylation probes localized within the exon or at a maximum distance of 50 nt from it). We called this set “m_exons pool”;
      • 2
        among the “m_exons” pool, only those that were represented in the identified circRNAs were further retained;
      • 3
        each circRNA was assigned a specific methylation level, derived from the methylation signals of its exons. If a circRNA contained more than one “m_exons”, then the average of the methylation levels of the “m_exons” was calculated;
      • 4
        The same workflow was repeated for methylation at gene promoters (considering all those methylation probes localized within 500 bp upstream of the transcription start site).
      We hence retained only those circRNAs for which methylation data were available both at gene body and promoter level. The correlation between circRNA expression level (in terms of counts) and methylation level (in terms of beta values) was evaluated through a Spearman correlation test. The statistical analysis was performed using the R software.

      2.10 Availability of data and materials

      The RNA-seq data from this publication have been deposited in the GEO database (accession number: GSE161196).

      3. Results

      3.1 CircRNA expression profile in MS patients by RNA-seq

      To globally explore the landscape of circRNAs in MS, we performed total RNA-seq experiments on RNA extracted from PBMCs of an US cohort of 10 RR-MS patients and 10 healthy controls (discovery cohort); subjects were matched for age, gender, and geographical origin (Supplementary Table 1). Prior to RNA-seq library preparation, RNAs were treated with RNase R, with the aim of enriching the samples for resistant circular products. We obtained an average of 25.6 million reads per sample. We identified 5663 circRNAs (detected by at least 5 back-spliced reads in 10 samples) using the DCC software; 12 circRNAs were sustained by more than 500 reads. The circRNAs characterized by the highest number of back-spliced reads derived from the SMARCA5 and CSNK1G3 host genes (Fig. 1a ). The majority of the detected circRNAs (74%) were already annotated in circBase (
      • Glažar P.
      • Papavasileiou P.
      • Rajewsky N.
      circBase: a database for circular RNAs.
      ). The identified circRNAs were distributed on all chromosomes, with chromosome 1 (the largest one) containing the highest number of circRNAs (Fig. 1b). The majority of circRNAs map within exons, whereas only a small fraction of them is represented by intronic or intergenic circRNAs (Fig. 1c).
      Fig 1
      Fig. 1Characteristics of the circRNAs identified in PBMCs of MS patients/healthy controls.
      (a) The figure shows the number of the back-spliced reads supporting each circRNA, distributed along chromosomes. CircRNAs are grouped according to their position on chromosomes (X axis) and represented as colored dots. The average number of the different back-spliced reads is reported on the Y axis. The arrows indicate the two circRNAs characterized by the highest number of back-spliced reads; the name of the corresponding host gene is also shown.
      (b) Distribution of circRNAs according to their chromosomal location. The number of circRNAs is reported on the Y axis.
      (c) Classification of circRNAs according to the position of their start and end (in exonic, intronic, or intergenic regions), depicted in a logarithmic scale. Exon-intron and exon-intergenic categories also include circRNAs whose start and end map in intronic-exonic and in intergenic-exonic regions, respectively.
      Differential expression analysis revealed 166 dysregulated circRNAs in RR-MS patients (P<0.05), 125 of which are downregulated (Fig. 2a and b and Supplementary Table 2). For the validation step, we selected 3 circRNAs among the top dysregulated ones (having a P<0.005, an absolute logFC≥0.75, and supported on average by at least 8 back-spliced reads), deriving from PGAP3, MED13L, and SMARCC1 genes (circBase IDs: hsa_circ_0007990, hsa_circ_0003059, and hsa_circ_0065301, respectively). We designed RT-PCR assays using divergent primer couples to confirm the presence of the annotated back-spliced junction by direct sequencing (Supplementary Fig. 1a). Semi-quantitative real-time RT-PCR assays, performed on the same samples that underwent RNA-seq, confirmed the downregulation trend detected by RNA-seq analysis, although the statistical significance was reached only for the PGAP3 circRNA (P=0.040) (Fig. 2c). The hsa_circ_0007990 circRNA consists of PGAP3 (Post-GPI Attachment to Proteins 3) exons 3 and 2, joined by the back-spliced junction, as verified by Sanger sequencing (Supplementary Fig. 1a). The downregulation of this circRNA was then confirmed in a US replication cohort including 18 RR-MS patients and 9 healthy controls (log2FC: -0.44, P=0.049) (Supplementary Fig. 2a) and in an independent Italian case-control cohort, including 19 RR-MS patients and 20 healthy controls (log2FC: -0.40, P=0.0075) (Supplementary Fig. 2b). We also detected a broad expression of the PGAP3 circRNA in selected MS-relevant human tissues (Supplementary Fig. 2c).
      Fig 2
      Fig. 2CircRNA profile in MS patients.
      (a) The heat map, built using the read counts of the top dysregulated circRNAs, as calculated by the DCC software, represents the 166 differentially expressed circRNAs between MS cases and healthy controls.
      (b) Volcano plot representing circRNAs differentially expressed in MS patients respect to controls. Red and blue dots indicate significantly up- and downregulated circRNAs, respectively. The grey shades highlight significantly dysregulated circRNAs that are characterized by a log2FC≥1 or a log2FC≤-1. The arrows indicate the circRNAs chosen for the validation step by semi-quantitative real-time RT–PCR.
      (c) Violin plots showing the expression levels of three circRNAs derived from PGAP3 (post-GPI attachment to proteins 3), MED13L (mediator complex subunit 13 like) and SMARCC1 (SWI/SNF related, matrix associated, actin dependent regulator of chromatin subfamily c member 1) genes by semi-quantitative real-time RT–PCRs, using specific divergent primer couples and HMBS as housekeeping gene. Results are presented as normalized rescaled values. A table showing the log2FC as obtained from RNA-seq analysis and the real-time RT-PCR assays is also shown. After checking for normal distribution of data, t-tests were applied for PGAP3 and SMARCC1 circRNAs, whereas Mann-Whitney test was performed for MED13L circRNA. *: one-tailed P<0.05; ns: not significant.
      (d) Top 25 differentially expressed GO terms for dysregulated circRNAs, identified by the “topGO” Bioconductor tool, using the “biological processes” and the “molecular functions” databases, ordered by elimFisher P value. The width of the dots indicates the percentage of the enriched genes out of the total number of genes belonging to each term.
      A GO enrichment analysis was performed on the 166 differentially expressed circRNAs, querying the biological processes (BP) and molecular functions (MF) databases (Supplementary Table 3). Interestingly, among the BP top hits, we observed an enrichment in immunity-related terms (i.e., positive regulation of innate immune response, T cell costimulation, stimulatory C-type lectin receptor signaling pathway), that could underlie a crucial deregulation of genes involved in immunological responses. In addition, among the MF top hits, we observed terms related to gene transcription and chromatin regulation (Fig. 2d).

      3.2 MS-associated SNPs influence the expression of circRNAs

      Given the differential expression of circRNAs in MS patients, and the demonstration of a genotype-dependent regulation of the STAT3 circRNA (
      • Paraboschi E.M.
      • Cardamone G.
      • Soldà G.
      • Duga S.
      • Asselta R.
      Interpreting non-coding genetic variation in multiple sclerosis genome-wide associated regions.
      ), we decided to globally explore the possibility that circRNAs may be modulated by MS-associated genetic variants.
      We genotyped the same 10 RR-MS patients and 10 healthy controls, chosen for the RNA-seq experiment (discovery cohort, Supplementary Table 1), using the Infinium® HumanCore-24 v1.1 BeadChip, containing 307,342 markers spread all over the genome. A total of 305,441 markers passed the QC steps, with an average call rate of 99.57% and an average genotype missingness <10%. To increase the coverage of genotyping data, we performed an imputation step through the Michigan Imputation Server, obtaining a total of 8,801,472 high-quality autosomal markers. Among the imputed markers, the MS-associated variants derived from the most recent meta-analysis of MS GWAS (
      International Multiple Sclerosis Genetics Consortium (IMSGC)
      Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility.
      ) were selected to perform an expression quantitative trait loci (eQTL) analysis, aimed at correlating the genotypes of the variants with the expression levels of circRNAs. In details, we selected both genome-wide associated (P<10−8) and strongly suggestive (10−8<P<10−5) variants, for a total of 317 markers; 285 (89%) of these were genotyped or imputed with a good quality (r2>0.3) in our cohort. The genotype data of these 285 SNPs were hence integrated with RNA-seq data using the FastQTL tool. The analysis evidenced a total of 89 significant eQTLs (Supplementary Table 4); in details, 43 SNPs were associated with the expression of one circRNA and 20 SNPs with the expression of at least two different circRNAs (Fig. 3a ). The rs7295386 polymorphism is particularly interesting, since it affects the expression of four circRNAs, deriving from three different host genes (Fig. 3a).
      Fig 3
      Fig. 3MS-associated SNPs regulate the formation of circRNAs.
      (a) The figure shows the number of circRNAs regulated by MS-associated SNPs, that are distributed along the chromosomes, and represented as colored dots. Only significant correlations are shown. The green dot indicates the most significant hit.
      (b,c) Boxplots representing the expression levels of the ATF1 hsa_circ_0098746 in the US cohort (b), and in the Italian one (c) by semi-quantitative real-time RT–PCRs, using specific divergent primer couples and HMBS as housekeeping gene. Results are presented as normalized rescaled values. Subjects are grouped according to the rs7295386 genotype. *: P<0.05.
      To validate the in-silico data, we decided to measure the expression of the top hit derived from the eQTL analysis, the ATF1 hsa_circ_0098746 (whose back-spliced junction was verified by Sanger sequencing, Supplementary Fig. 1b), which was associated with the rs7295386 SNP. We performed semi-quantitative real-time RT-PCR in the US MS patients/healthy controls (belonging to the discovery and replication cohort), and grouped the subjects according to the genotype of the rs7295386 SNP. We found a significant correlation between the genotype of the SNP and the expression of the circRNA (ANOVA P=0.0077, Fig. 3b), and we further confirmed this correlation also in the independent Italian replication cohort (ANOVA P=0.0027, Fig. 3c).

      3.3 Methylation at gene body level is correlated with circRNA expression

      To verify whether the methylation profile could affect the expression level of circRNAs, we took advantage of a previous RNA-seq experiment, in which we analyzed the back-splicing profile of the Jurkat cell line (
      • Paraboschi E.M.
      • Cardamone G.
      • Soldà G.
      • Duga S.
      • Asselta R.
      Interpreting non-coding genetic variation in multiple sclerosis genome-wide associated regions.
      ). Methylation data for the same cell type, in terms of beta values, are available in the UCSC Genome Browser. We explored the possibility that the methylation status at gene promoter or gene body could be correlated with circRNA expression. The Spearman ranking analysis evidenced a significant positive correlation between circRNA levels and gene-body methylation status (P=0.0021, rho=0.16), whereas no correlation was observed with methylation at promoters (P=0.26, rho=0.06; Fig. 4a ). Interestingly, a negative correlation was highlighted between the level of methylation at gene body and promoters (P=0.0090, rho=-0.14), confirming the observation that expressed genes display an opposite methylation profile at promoters compared to gene bodies (
      • Maunakea A.K.
      • Chepelev I.
      • Cui K.
      • Zhao K.
      Intragenic DNA methylation modulates alternative splicing by recruiting MeCP2 to promote exon recognition.
      ;
      • Shayevitch R.
      • Askayo D.
      • Keydar I.
      • Ast G.
      The importance of DNA methylation of exons on alternative splicing.
      ).
      Fig 4
      Fig. 4Gene-body methylation is associated with circRNA levels.
      (a) Correlogram showing the correlation level between circRNA expression and methylation profile at gene promoters (beta promoter) or gene bodies (beta gene body). Circles are drawn if a significant correlation is present, their size represents the correlation level, whereas their color indicates the sign of the correlation. Rho values, as calculated by the Spearman rank test, are shown in each box.
      (b, c) Boxplots showing the expression of circRNAs, as log(counts), on the basis of gene-body (b) and promoter (c) methylation level, in terms of beta values. Outliers are presented as dots; *: p< 0.05; ns: not significant.
      To better understand how methylation impacts on the back-splicing profile, we clustered circRNAs, on the basis of the methylation beta value, in a “low” and a “high” methylation-level group. This analysis was performed both for promoter and gene body methylation, setting the median methylation level as threshold for the assignment to a specific group. Concerning gene body methylation, the Wilcoxon test evidenced a significant difference in circRNA expression between the “low” and “high” methylation groups (P= 0.018, Fig. 4b), with a higher methylation level associated with an increased circRNA expression (reads: median=9, IQR=10 vs median=7, IQR=6). No significant difference was instead observed in the case of methylation at promoters (P=0.29, reads: median=8, IQR=9 vs median=8, IQR=7, Fig. 4c).

      4. Discussion

      CircRNAs are increasingly involved in several diseases, including autoimmune and neurodegenerative disorders, although the molecular mechanisms regulating their expression are still largely unknown.
      Our previous demonstration of a dysregulated circRNA in MS patients (GSDMB circRNA) (
      • Cardamone G.
      • Paraboschi E.M.
      • Rimoldi V.
      • Duga S.
      • Soldà G.
      • Asselta R.
      The characterization of GSDMB splicing and backsplicing profiles identifies novel isoforms and a circular RNA that are dysregulated in multiple sclerosis.
      ) underlined the need to describe a complete map of circRNA expression in MS, a disease in which alteration of AS and RNA processing have been repeatedly reported (
      • Paraboschi E.M
      • Cardamone G.
      • Rimoldi V.
      • Gemmati D.
      • Spreafico M.
      • Duga S.
      • Soldà G.
      • Asselta R.
      Meta-analysis of multiple sclerosis microarray data reveals dysregulation in RNA splicing regulatory genes.
      ;
      • Spurlock 3rd C.F.
      • Tossberg J.T.
      • Guo Y.
      • Sriram S.
      • Crooke 3rd P.S.
      • Aune T.M.
      Defective structural RNA processing in relapsing-remitting multiple sclerosis.
      ;
      • Hecker M.
      • Rüge A.
      • Putscher E.
      • Boxberger N.
      • Rommer P.S.
      • Fitzner B.
      • Zettl U.K.
      Aberrant expression of alternative splicing variants in multiple sclerosis – a systematic review.
      ;
      • Cardamone G.
      • Paraboschi E.M.
      • Soldà G.
      • Cantoni C.
      • Supino D.
      • Piccio L.
      • Duga S.
      • Asselta R.
      Not only cancer: the long non-coding RNA MALAT1 affects the repertoire of alternatively spliced transcripts and circular RNAs in multiple sclerosis.
      ;
      • Evsyukova I.
      • Somarelli J.A.
      • Gregory S.G.
      • Garcia-Blanco M.A.
      Alternative splicing in multiple sclerosis and other autoimmune diseases.
      ). An attempt in this direction was indeed already performed by
      • Iparraguirre L.
      • Muñoz-Culla M.
      • Prada-Luengo I.
      • Castillo-Triviño T.
      • Olascoaga J.
      • Otaegui D.
      Circular RNA profiling reveals that circular RNAs from ANXA2 can be used as new biomarkers for multiple sclerosis.
      , but their work suffered from the limitation of using microarrays, which, by design, only include probes for already known circRNAs. More recently, the same research group also analyzed the circRNA profile in MS patients’ and controls’ leucocytes by RNA-seq analysis (
      • Iparraguirre L.
      • Alberro A.
      • Sepúlveda L.
      • Osorio-Querejeta I.
      • Moles L.
      • Castillo-Triviño T.
      • Hansen T.B.
      • Muñoz-Culla M.
      • Otaegui D.
      RNA-Seq profiling of leukocytes reveals a sex-dependent global circular RNA upregulation in multiple sclerosis and 6 candidate biomarkers.
      ). There was no overlap between the dysregulated circRNAs found in our work and the top hits listed in their manuscript, but different cell populations were analyzed (leucocytes, containing also neutrophils, and PBMCs). Additionally, other groups have previously analyzed the circRNA profile in PBMCs of RR-MS patients (either in the relapsing or remitting phase of the disease) and healthy controls (
      • Zurawska A.E.
      • Mycko M.P.
      • Selmaj I.
      • Raine C.S.
      • Selmaj K.W.
      Multiple sclerosis: circRNA profile defined reveals links to B-cell function.
      ;
      • Mycko M.P.
      • Zurawska A.E.
      • Selmaj I.
      • Selmaj K.W.
      Impact of diminished expression of circRNA on multiple sclerosis pathomechanisms.
      ) by microarray-based studies. Also in this case, however, we did not find overlaps between their top hits and the circRNAs described in this work, again possibly due to different technical approaches. Our exploratory RNA-seq experiment revealed 166 differentially expressed circRNAs, 75% of which are downregulated, in line with the results obtained by Iparraguirre et al. in the microarray analysis (
      • Iparraguirre L.
      • Muñoz-Culla M.
      • Prada-Luengo I.
      • Castillo-Triviño T.
      • Olascoaga J.
      • Otaegui D.
      Circular RNA profiling reveals that circular RNAs from ANXA2 can be used as new biomarkers for multiple sclerosis.
      ). This trend towards a prevalent circRNA downregulation was also described in PBMCs of patients affected by another autoimmune disease, systemic lupus erythematosus (
      • Liu C.X.
      • Li X.
      • Nan F.
      • Jiang S.
      • Gao X.
      • Guo S.K.
      • Xue W.
      • Cui Y.
      • Dong K.
      • Ding H.
      • Qu B.
      • Zhou Z.
      • Shen N.
      • Yang L.
      • Chen L.L.
      Structure and degradation of circular RNAs regulate PKR activation in innate immunity.
      ). Since circRNAs are known to control gene expression by miRNA sponging activity, the observed downregulation may cause an alteration in miRNA levels, thus leading to microenvironmental changes possibly modulating MS development, as suggested by
      • Mycko M.P.
      • Zurawska A.E.
      • Selmaj I.
      • Selmaj K.W.
      Impact of diminished expression of circRNA on multiple sclerosis pathomechanisms.
      .
      Among the top dysregulated circRNAs, the hsa_circ_0007990 seemed particularly interesting, as the PGAP3 host gene encodes a protein shown to be relevant for the control of autoimmune responses. Indeed, it encodes a protein responsible for the remodeling of fatty acids of glycosylphosphatidylinositol (GPI)-anchors proteins; more importantly, PGAP3 knockout mice showed an enhanced T cell response and exacerbate the experimental autoimmune encephalomyelitis phenotype, as well as several autoimmune-like symptoms (
      • Murakami H.
      • Wang Y.
      • Hasuwa H.
      • Maeda Y.
      • Kinoshita T.
      • Murakami Y.
      Enhanced response of T lymphocytes from Pgap3 knockout mouse: Insight into roles of fatty acid remodeling of GPI anchored proteins.
      ;
      • Wang Y.
      • Murakami Y.
      • Yasui T.
      • Wakana S.
      • Kikutani H.
      • Kinoshita T.
      • Maeda Y.
      Significance of glycosylphosphatidylinositol-anchored protein enrichment in lipid rafts for the control of autoimmunity.
      ). We confirmed the downregulation of the PGAP3 circRNA in the US and Italian replication cohorts, providing a possible novel RNA-based biomarker. Although we are aware that the number of patients considered in this analysis is small, the potential role as a diagnostic/prognostic factor, in terms of disease progression, should be considered, and this promising molecular marker needs to be further replicated in an enlarged cohort.
      Our previous observation that STAT3 hsa_circ_0043813 level was regulated by the genotype of an MS-associated SNP (
      • Paraboschi E.M.
      • Cardamone G.
      • Soldà G.
      • Duga S.
      • Asselta R.
      Interpreting non-coding genetic variation in multiple sclerosis genome-wide associated regions.
      ) prompted us to hypothesize that MS-associated variants may influence the biogenesis of circRNAs. A study performed by
      • Putscher E.
      • Hecker M.
      • Fitzner B.
      • Boxberger N.
      • Schwartz M.
      • Koczan D.
      • Lorenz P.
      • Zettl U.K.
      Genetic risk variants for multiple sclerosis are linked to differences in alternative pre-mRNA splicing.
      analyzed the effect on AS of SNPs located in genetic risk loci for MS, and identified a number of GWAS-associated SNPs that act as splice-QTLs in B lymphocytes, but only the effect on linear splicing was evaluated. By investigating the genetic background of all the patients analyzed for their circRNA profile, we indeed found 89 significant correlations between genotypes of MS-associated SNPs and the expression level of circRNAs, which might be defined as circ-QTLs. It has to be underlined that a total of 63 out of 285 tested SNPs were correlated with circRNA expression level, suggesting that about the 22% of genome-wide MS-associated SNPs could impact on the disease pathogenesis by altering circRNA profiles.
      Importantly, this analysis also revealed that the genotype of a SNP may influence the expression of multiple circRNAs, as in the case of the rs7295386 SNP regulating four circRNAs. In this context,
      • Ahmed I.
      • Karedath T.
      • Al-Dasim F.M.
      • Malek J.A.
      Identification of human genetic variants controlling circular RNA expression.
      , taking advantages of the RNA-seq data of 1000 Genomes Project on lymphoblastoid cell lines (
      • Lappalainen T.
      • Sammeth M.
      • Friedländer M.R.
      • 't Hoen P.A.
      • Monlong J.
      • Rivas M.A.
      • Gonzàlez-Porta M.
      • Kurbatova N.
      • Griebel T.
      • Ferreira P.G.
      • Barann M.
      • Wieland T.
      • Greger L.
      • van Iterson M.
      • Almlöf J.
      • Ribeca P.
      • Pulyakhina I.
      • Esser D.
      • Giger T.
      • Tikhonov A.
      • Sultan M.
      • Bertier G.
      • MacArthur D.G
      • Lek M.
      • Lizano E.
      • Buermans H.P.J
      • Padioleau I.
      • Schwarzmayr T.
      • Karlberg O.
      • Ongen H.
      • Kilpinen H.
      • Beltran S.
      • Gut M.
      • Kahlem K.
      • Amstislavskiy V.
      • Stegle O.
      • Pirinen M.
      • Montgomery S.B.
      • Donnelly P.
      • McCarthy M.I.
      • Flicek P.
      • Strom T.M.
      • Consortium Geuvadis
      • Lehrach H.
      • Schreiber S.
      • Sudbrak R.
      • Carracedo A.
      • Antonarakis S.E.
      • Häsler R.
      • Syvänen A.C.
      • van Ommen G.J.
      • Brazma A.
      • Meitinger T.
      • Rosenstiel P.
      • Guigó R.
      • Gut I.G.
      • Estivill X.
      • Dermitzakis E.T.
      Transcriptome and genome sequencing uncovers functional variation in humans.
      ), found 139,485 circ-QTLs for 2260 circRNAs, and also observed that most circ-QTLs are independent of eQTLs acting on gene expression and do not impact on the expression of the host genes.
      • Liu Z.
      • Ran Y.
      • Tao C.
      • Li S.
      • Chen J.
      • Yang E.
      Detection of circular RNA expression and related quantitative trait loci in the human dorsolateral prefrontal cortex.
      , instead, analyzed human dorsolateral prefrontal cortex samples from the CommonMind Consortium. Interestingly, the rs7295386 appeared among the 196,255 circ-QTLs, again being associated with the ATF1 hsa_circ_0098746, that was shown to be regulated by other 246 circ-QTLs. They also found that many of the most significant circ-QTLs (for each circRNA) are enriched within GWAS signals of several complex diseases, although their analysis did not include MS.
      Besides genetic variants, we hypothesized that epigenetic features, such as DNA methylation, may represent a further regulatory mechanism to tune circRNA expression. Indeed,
      • Ferreira H.J.
      • Davalos V.
      • de Moura M.C.
      • Soler M.
      • Perez-Salvia M.
      • Bueno-Costa A.
      • Setien F.
      • Moran S.
      • Villanueva A.
      • Esteller M.
      Circular RNA CpG island hypermethylation-associated silencing in human cancer.
      have already demonstrated that cancer-specific promoter CpG island hypermethylation is associated with a decreased expression of both circRNAs and their host genes, and intragenic methylation is known to modulate AS (
      • Maunakea A.K.
      • Chepelev I.
      • Cui K.
      • Zhao K.
      Intragenic DNA methylation modulates alternative splicing by recruiting MeCP2 to promote exon recognition.
      ;
      • Shayevitch R.
      • Askayo D.
      • Keydar I.
      • Ast G.
      The importance of DNA methylation of exons on alternative splicing.
      ). Moreover,
      • Xu T.
      • Wang L.
      • Jia P.
      • Song X.
      • Zhao Z.
      An integrative transcriptomic and methylation approach for identifying differentially expressed circular RNAs associated with DNA methylation change.
      showed that circRNAs with differentially methylated sites were expressed differently in tumors versus adjacent normal samples, while their parental genes were not, thus hypothesizing that some aberrant DNA methylation events might only affect the processing of pre-mRNA to generate circRNAs but not the process generating linear RNAs. In our work, we evaluated whether DNA methylation could impact on back-splicing modulation in the Jurkat cell line. Our analysis, although indirect, evidenced the existence of a positive correlation between circRNA expression and methylation at gene body level, thus suggesting that epigenetic features may play an important role in the definition of the cell circRNA pool. We also analyzed publicly available DNA methylation data derived from peripheral blood of 98 female MS patients and 104 sex-matched healthy controls (GSE106648) (
      • Kular L.
      • Liu Y.
      • Ruhrmann S.
      • Zheleznyakova G.
      • Marabita F.
      • Gomez-Cabrero D.
      • James T.
      • Ewing E.
      • Lindén M.
      • Górnikiewicz B.
      • Aeinehband S.
      • Stridh P.
      • Link J.
      • Andlauer T.F.M.
      • Gasperi C.
      • Wiendl H.
      • Zipp F.
      • Gold R.
      • Tackenberg B.
      • Weber F.
      • Hemmer B.
      • Strauch K.
      • Heilmann-Heimbach S.
      • Rawal R.
      • Schminke U.
      • Schmidt C.O.
      • Kacprowski T.
      • Franke A.
      • Laudes M.
      • Dilthey A.T.
      • Celius E.G.
      • Søndergaard H.B.
      • Tegnér J.
      • Harbo H.F.
      • Oturai A.B.
      • Olafsson S.
      • Eggertsson H.P.
      • Halldorsson B.V.
      • Hjaltason H.
      • Olafsson E.
      • Jonsdottir I.
      • Stefansson K.
      • Olsson T.
      • Piehl F.
      • Ekström T.J.
      • Kockum I.
      • Feinberg A.P.
      • Jagodic M.
      DNA methylation as a mediator of HLA-DRB1*15:01 and a protective variant in multiple sclerosis.
      ): 36 of the genes that were differentially methylated in cases were also evidenced in our RNA-seq analysis as having a dysregulated circRNA profile. Although we cannot directly integrate the results obtained from these two different cohorts, this analysis further supports the need to explore the correlation between circRNA and methylation profiles to gain more insights on the regulatory features modulating circRNA landscape, and ultimately on disease pathogenesis.
      In conclusion, this work establishes a comprehensive profile of circRNA expression by RNA-seq analysis in PBMCs of MS patients, suggesting that a significant proportion (at least 20%) of MS-associated variants may influence the expression levels of circRNAs, acting as circ-QTLs and providing support for the role of DNA methylation in regulating circRNA biogenesis.

      Funding

      This work was partially supported by FRRB (Fondazione Regionale per la Ricerca Biomedica), project number: 1734478.

      CRediT authorship contribution statement

      Giulia Cardamone: Conceptualization, Investigation, Writing – original draft. Elvezia Maria Paraboschi: Conceptualization, Formal analysis, Writing – review & editing. Giulia Soldà: Conceptualization, Writing – review & editing. Giuseppe Liberatore: Resources. Valeria Rimoldi: Investigation. Javier Cibella: Investigation. Federica Airi: Investigation. Veronica Tisato: Resources. Claudia Cantoni: Resources. Francesca Gallia: Resources. Donato Gemmati: Resources, Writing – review & editing. Laura Piccio: Resources, Writing – review & editing. Stefano Duga: Conceptualization, Writing – review & editing. Eduardo Nobile-Orazio: Conceptualization, Writing – review & editing. Rosanna Asselta: Conceptualization, Writing – original draft, Writing – review & editing.

      Declaration of Competing Interest

      The authors report no conflict of interest.

      Acknowledgments

      This manuscript is dedicated in memory of Prof. Stefano Duga, whose immense support and guidance had inspired us immensely.

      References

        • Lasda E.
        • Parker R.
        Circular RNAs: diversity of form and function.
        RNA. 2014; 20: 1829-1842https://doi.org/10.1261/rna.047126.114
        • Salzman J.
        • Gawad C.
        • Wang P.L.
        • Lacayo N.
        • Brown P.O.
        Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types.
        PLoS One. 2012; 7: e30733https://doi.org/10.1371/journal.pone.0030733
        • Memczak S.
        • Jens M.
        • Elefsinioti A.
        • Torti F.
        • Krueger J.
        • Rybak A.
        • Maier L.
        • Mackowiak S.D.
        • Gregersen L.H.
        • Munschauer M.
        • Loewer A.
        • Ziebold U.
        • Landthaler M.
        • Kocks C.
        • le Noble F.
        • Rajewsky N.
        Circular RNAs are a large class of animal RNAs with regulatory potency.
        Nature. 2013; 495: 333-338https://doi.org/10.1038/nature11928
        • Jeck W.R.
        • Sorrentino J.A.
        • Wang K.
        • Slevin M.K.
        • Burd C.E.
        • Liu J.
        • Marzluff W.F.
        • Sharpless N.E.
        Circular RNAs are abundant, conserved, and associated with ALU repeats.
        RNA. 2013; 19: 141-157https://doi.org/10.1261/rna.035667.112
        • Ebbesen K.K.
        • Kjems J.
        • Hansen T.B.
        Circular RNAs: identification, biogenesis and function.
        Biochim. Biophys. Acta. 2016; 1859: 163-168https://doi.org/10.1016/j.bbagrm.2015.07.007
        • Ashwal-Fluss R.
        • Meyer M.
        • Pamudurti N.R.
        • Ivanov A.
        • Bartok O.
        • Hanan M.
        • Evantal N.
        • Memczak S.
        • Rajewsky N.
        • Kadener S.
        circRNA biogenesis competes with pre-mRNA splicing.
        Mol. Cell. 2014; 56: 55-66https://doi.org/10.1016/j.molcel.2014.08.019
        • Hansen T.B.
        • Jensen T.I.
        • Clausen B.H.
        • Bramsen J.B.
        • Finsen B.
        • Damgaard C.K.
        • Kjems J.
        Natural RNA circles function as efficient microRNA sponges.
        Nature. 2013; 495: 384-388https://doi.org/10.1038/nature11993
        • Legnini I.
        • Di Timoteo G.
        • Rossi F.
        • Morlando M.
        • Briganti F.
        • Sthandier O.
        • Fatica A.
        • Santini T.
        • Andronache A.
        • Wade M.
        • Laneve P.
        • Rajewsky N.
        • Bozzoni I.
        Circ-ZNF609 is a circular RNA that can be translated and functions in myogenesis.
        Mol. Cell. 2017; 66: 22-37https://doi.org/10.1016/j.molcel.2017.02.017
        • Pamudurti N.R.
        • Bartok O.
        • Jens M.
        • Ashwal-Fluss R.
        • Stottmeister C.
        • Ruhe L.
        • Hanan M.
        • Wyler E.
        • Perez-Hernandez D.
        • Ramberger E.
        • Shenzis S.
        • Samson M.
        • Dittmar G.
        • Landthaler M.
        • Chekulaeva M.
        • Rajewsky N.
        • Kadener S.
        Translation of CircRNAs.
        Mol. Cell. 2017; 66: 9-21https://doi.org/10.1016/j.molcel.2017.02.021
        • Conn S.J.
        • Pillman K.A.
        • Toubia J.
        • Conn V.M.
        • Salmanidis M.
        • Phillips C.A.
        • Roslan S.
        • Schreiber A.W.
        • Gregory P.A.
        • Goodall G.J.
        The RNA binding protein quaking regulates formation of circRNAs.
        Cell. 2015; 160: 1125-1134https://doi.org/10.1016/j.cell.2015.02.014
        • Chen L.L.
        • Yang L.
        Regulation of circRNA biogenesis.
        RNA Biol. 2015; 12: 381-388https://doi.org/10.1080/15476286.2015.1020271
        • Starke S.
        • Jost I.
        • Rossbach O.
        • Schneider T.
        • Schreiner S.
        • Hung L.H.
        • Bindereif A.
        Exon circularization requires canonical splice signals.
        Cell Rep. 2015; 10: 103-111https://doi.org/10.1080/10.1016/j.celrep.2014.12.002
        • Kramer M.C.
        • Liang D.
        • Tatomer D.C.
        • Gold B.
        • March Z.M.
        • Cherry S.
        • Wilusz J.E.
        Combinatorial control of Drosophila circular RNA expression by intronic repeats, hnRNPs, and SR proteins.
        Genes Dev. 2015; 29: 2168-2182https://doi.org/10.1101/gad.270421.115
        • Abdelmohsen K.
        • Panda A.C.
        • Munk R.
        • Grammatikakis I.
        • Dudekula D.B.
        • De S.
        • Gorospe M.
        Identification of HuR target circular RNAs uncovers suppression of PABPN1 translation by CircPABPN1.
        RNA Biol. 2017; 14: 361-369https://doi.org/10.1080/15476286.2017.1279788
        • Han B.
        • Chao J.
        • Yao H.
        Circular RNA and its mechanisms in disease: from the bench to the clinic.
        Pharmacol. Ther. 2018; 187: 31-44https://doi.org/10.1016/j.pharmthera.2018.01.010
        • Cocquerelle C.
        • Mascrez B.
        • Hétuin D.
        • Bailleul B.
        Mis-splicing yields circular RNA molecules.
        FASEB J. 1993; 7: 155-160https://doi.org/10.1096/fasebj.7.1.7678559
        • Maass P.G.
        • Glažar P.
        • Memczak S.
        • Dittmar G.
        • Hollfinger I.
        • Schreyer L.
        • Sauer A.V
        • Toka O.
        • Aiuti A.
        • Luft F.C.
        • Rajewsky N.
        A map of human circular RNAs in clinically relevant tissues.
        J. Mol. Med. 2017; 95: 1179-1189https://doi.org/10.1007/s00109-017-1582-9
        • Thompson A.J.
        • Baranzini S.E.
        • Geurts J.
        • Hemmer B.
        • Ciccarelli O.
        Multiple sclerosis.
        Lancet. 2018; 391: 1622-1636https://doi.org/10.1016/S0140-6736(18)30481-1
        • Lublin F.D.
        • Reingold S.C.
        • Cohen J.A.
        • Cutter G.R.
        • Sørensen P.S.
        • Thompson A.J.
        • Wolinsky J.S.
        • Balcer L.J.
        • Banwell B.
        • Barkhof F.
        • Bebo Jr B.
        • Calabresi P.A.
        • Clanet M.
        • Comi G.
        • Fox R.J.
        • Freedman M.S.
        • Goodman A.D.
        • Inglese M.
        • Kappos L.
        • Kieseier B.C.
        • Lincoln J.A.
        • Lubetzki C.
        • Miller A.E.
        • Montalban X.
        • O'Connor P.W.
        • Petkau J.
        • Pozzilli C.
        • Rudick R.A.
        • Sormani M.P.
        • Stüve O.
        • Waubant E.
        • Polman C.H.
        Defining the clinical course of multiple sclerosis: the 2013 revisions.
        Neurology. 2014; 83: 278-286https://doi.org/10.1212/WNL.0000000000000560
        • Paraboschi E.M
        • Cardamone G.
        • Rimoldi V.
        • Gemmati D.
        • Spreafico M.
        • Duga S.
        • Soldà G.
        • Asselta R.
        Meta-analysis of multiple sclerosis microarray data reveals dysregulation in RNA splicing regulatory genes.
        Int. J. Mol. Sci. 2015; 16: 23463-23481https://doi.org/10.3390/ijms161023463
        • Spurlock 3rd C.F.
        • Tossberg J.T.
        • Guo Y.
        • Sriram S.
        • Crooke 3rd P.S.
        • Aune T.M.
        Defective structural RNA processing in relapsing-remitting multiple sclerosis.
        Genome Biol. 2015; 16: 58https://doi.org/10.1186/s13059-015-0629-x
        • Hecker M.
        • Rüge A.
        • Putscher E.
        • Boxberger N.
        • Rommer P.S.
        • Fitzner B.
        • Zettl U.K.
        Aberrant expression of alternative splicing variants in multiple sclerosis – a systematic review.
        Autoimmun. Rev. 2019; 18: 721-732https://doi.org/10.1016/j.autrev.2019.05.010
        • D'Ambrosio A.
        • Pontecorvo S.
        • Colasanti T.
        • Zamboni S.
        • Francia A.
        • Margutti P.
        Peripheral blood biomarkers in multiple sclerosis.
        Autoimmun. Rev. 2015; 14: 1097-1110https://doi.org/10.1016/j.autrev.2015.07.014
        • Cardamone G.
        • Paraboschi E.M.
        • Rimoldi V.
        • Duga S.
        • Soldà G.
        • Asselta R.
        The characterization of GSDMB splicing and backsplicing profiles identifies novel isoforms and a circular RNA that are dysregulated in multiple sclerosis.
        Int. J. Mol. Sci. 2017; 18: 576https://doi.org/10.3390/ijms18030576
        • Iparraguirre L.
        • Muñoz-Culla M.
        • Prada-Luengo I.
        • Castillo-Triviño T.
        • Olascoaga J.
        • Otaegui D.
        Circular RNA profiling reveals that circular RNAs from ANXA2 can be used as new biomarkers for multiple sclerosis.
        Hum. Mol. Genet. 2017; 26: 3564-3572https://doi.org/10.1093/hmg/ddx243
        • Paraboschi E.M.
        • Cardamone G.
        • Soldà G.
        • Duga S.
        • Asselta R.
        Interpreting non-coding genetic variation in multiple sclerosis genome-wide associated regions.
        Front Genet. 2018; 9: 647https://doi.org/10.3389/fgene.2018.00647
        • Holdt L.M.
        • Stahringer A.
        • Sass K.
        • Pichler G.
        • Kulak N.A.
        • Wilfert W.
        • Kohlmaier A.
        • Herbst A.
        • Northoff B.H.
        • Nicolaou A.
        • Gäbel G.
        • Beutner F.
        • Scholz M.
        • Thiery J.
        • Musunuru K.
        • Krohn K.
        • Mann M.
        • Teupser D.
        Circular non-coding RNA ANRIL modulates ribosomal RNA maturation and atherosclerosis in humans.
        Nat. Commun. 2016; 7: 12429https://doi.org/10.1038/ncomms12429
        • Ferreira H.J.
        • Davalos V.
        • de Moura M.C.
        • Soler M.
        • Perez-Salvia M.
        • Bueno-Costa A.
        • Setien F.
        • Moran S.
        • Villanueva A.
        • Esteller M.
        Circular RNA CpG island hypermethylation-associated silencing in human cancer.
        Oncotarget. 2018; 9: 29208-29219https://doi.org/10.18632/oncotarget.25673
        • Kleaveland B.
        • Shi C.Y.
        • Stefano J.
        • Bartel D.P.
        A network of noncoding regulatory RNAs acts in the mammalian brain.
        Cell. 2018; 174 (e17): 350-362https://doi.org/10.1016/j.cell.2018.05.022
        • Cardamone G.
        • Paraboschi E.M.
        • Soldà G.
        • Cantoni C.
        • Supino D.
        • Piccio L.
        • Duga S.
        • Asselta R.
        Not only cancer: the long non-coding RNA MALAT1 affects the repertoire of alternatively spliced transcripts and circular RNAs in multiple sclerosis.
        Hum. Mol. Genet. 2019; 28: 1414-1428https://doi.org/10.1093/hmg/ddy438
        • Cheng J.
        • Metge F.
        • Dieterich C.
        Specific identification and quantification of circular RNAs from sequencing data.
        Bioinformatics. 2016; 32: 1094-1096https://doi.org/10.1093/bioinformatics/btv656
        • Dobin A.
        • Davis C.A.
        • Schlesinger F.
        • Drenkow J.
        • Zaleski C.
        • Jha S.
        • Batut P.
        • Chaisson M.
        • Gingeras T.R.
        STAR: ultrafast universal RNA-seq aligner.
        Bioinformatics. 2013; 29: 15-21https://doi.org/10.1093/bioinformatics/bts635
        • Robinson M.D.
        • McCarthy D.J.
        • Smyth G.K.
        edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
        Bioinformatics. 2010; 26: 139-140https://doi.org/10.1093/bioinformatics/btp616
        • Dou Y.
        • Cha D.J.
        • Franklin J.L.
        • Higginbotham J.N.
        • Jeppesen D.K.
        • Weaver A.M.
        • Prasad N.
        • Levy S.
        • Coffey R.J.
        • Patton J.G.
        • Zhang B.
        Circular RNAs are down-regulated in KRAS mutant colon cancer cells and can be transferred to exosomes.
        Sci. Rep. 2016; 6: 37982https://doi.org/10.1038/srep37982
      1. A. Alexa, J. Rahnenfuhrer (2016), topGO: enrichment analysis for gene ontology. R package version 2.24.0.

        • R Core Team
        R: A Language and Environment for Statistical Computing.
        R Foundation for Statistical Computing, Vienna, Austria2018 (Available online at)
        • Das S.
        • Forer L.
        • Schönherr S.
        • Sidore C.
        • Locke A.E.
        • Kwong A.
        • Vrieze S.I.
        • Chew E.Y.
        • Levy S.
        • McGue M.
        • Schlessinger D.
        • Stambolian D.
        • Loh P.R.
        • Iacono W.G.
        • Swaroop A.
        • Scott L.J.
        • Cucca F.
        • Kronenberg F.
        • Boehnke M.
        • Abecasis G.R.
        • Fuchsberger C.
        Next-generation genotype imputation service and methods.
        Nat. Genet. 2016; 48: 1284-1287https://doi.org/10.1038/ng.3656
        • Loh P.R.
        • Danecek P.
        • Palamara P.F.
        • Fuchsberger C.
        • Reshef Y.A.
        • Finucane H.K.
        • Schoenherr S.
        • Forer L.
        • McCarthy S.
        • Abecasis G.R.
        • Durbin R.
        • Price A.L.
        Reference-based phasing using the haplotype reference consortium panel.
        Nat. Genet. 2016; 48: 1443-1448https://doi.org/10.1038/ng.3679
        • International Multiple Sclerosis Genetics Consortium (IMSGC)
        Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility.
        Science. 2019; 365: eaav7188https://doi.org/10.1126/science.aav7188
        • Ongen H.
        • Buil A.
        • Brown A.A.
        • Dermitzakis E.T.
        • Delaneau O.
        Fast and efficient QTL mapper for thousands of molecular phenotypes.
        Bioinformatics. 2016; 32: 1479-1485https://doi.org/10.1093/bioinformatics/btv722
        • Glažar P.
        • Papavasileiou P.
        • Rajewsky N.
        circBase: a database for circular RNAs.
        RNA. 2014; 20: 1666-1670https://doi.org/10.1261/rna.043687.113
        • Maunakea A.K.
        • Chepelev I.
        • Cui K.
        • Zhao K.
        Intragenic DNA methylation modulates alternative splicing by recruiting MeCP2 to promote exon recognition.
        Cell Res. 2013; 23: 1256-1269https://doi.org/10.1038/cr.2013.110
        • Shayevitch R.
        • Askayo D.
        • Keydar I.
        • Ast G.
        The importance of DNA methylation of exons on alternative splicing.
        RNA. 2018; 24: 1351-1362https://doi.org/10.1261/rna.064865.117
        • Evsyukova I.
        • Somarelli J.A.
        • Gregory S.G.
        • Garcia-Blanco M.A.
        Alternative splicing in multiple sclerosis and other autoimmune diseases.
        RNA Biol. 2010; 7: 462-473https://doi.org/10.4161/rna.7.4.12301
        • Iparraguirre L.
        • Alberro A.
        • Sepúlveda L.
        • Osorio-Querejeta I.
        • Moles L.
        • Castillo-Triviño T.
        • Hansen T.B.
        • Muñoz-Culla M.
        • Otaegui D.
        RNA-Seq profiling of leukocytes reveals a sex-dependent global circular RNA upregulation in multiple sclerosis and 6 candidate biomarkers.
        Hum. Mol. Genet. 2020; 18: 3361-3372https://doi.org/10.1093/hmg/ddaa219
        • Zurawska A.E.
        • Mycko M.P.
        • Selmaj I.
        • Raine C.S.
        • Selmaj K.W.
        Multiple sclerosis: circRNA profile defined reveals links to B-cell function.
        Neurol. Neuroimmunol. Neuroinflamm. 2021; 8: e1041https://doi.org/10.1212/NXI.0000000000001041
        • Mycko M.P.
        • Zurawska A.E.
        • Selmaj I.
        • Selmaj K.W.
        Impact of diminished expression of circRNA on multiple sclerosis pathomechanisms.
        Front Immunol. 2022; 13875994https://doi.org/10.3389/fimmu.2022.875994
        • Liu C.X.
        • Li X.
        • Nan F.
        • Jiang S.
        • Gao X.
        • Guo S.K.
        • Xue W.
        • Cui Y.
        • Dong K.
        • Ding H.
        • Qu B.
        • Zhou Z.
        • Shen N.
        • Yang L.
        • Chen L.L.
        Structure and degradation of circular RNAs regulate PKR activation in innate immunity.
        Cell. 2019; 177: 865-880https://doi.org/10.1016/j.cell.2019.03.046
        • Murakami H.
        • Wang Y.
        • Hasuwa H.
        • Maeda Y.
        • Kinoshita T.
        • Murakami Y.
        Enhanced response of T lymphocytes from Pgap3 knockout mouse: Insight into roles of fatty acid remodeling of GPI anchored proteins.
        Biochem. Biophys. Res. Commun. 2012; 427: 1235-1241https://doi.org/10.1016/j.bbrc.2011.12.116
        • Wang Y.
        • Murakami Y.
        • Yasui T.
        • Wakana S.
        • Kikutani H.
        • Kinoshita T.
        • Maeda Y.
        Significance of glycosylphosphatidylinositol-anchored protein enrichment in lipid rafts for the control of autoimmunity.
        J. Biol. Chem. 2013; 288: 25490-25499https://doi.org/10.1074/jbc.M113.492611
        • Putscher E.
        • Hecker M.
        • Fitzner B.
        • Boxberger N.
        • Schwartz M.
        • Koczan D.
        • Lorenz P.
        • Zettl U.K.
        Genetic risk variants for multiple sclerosis are linked to differences in alternative pre-mRNA splicing.
        Front. Immunol. 2022; 13931831https://doi.org/10.3389/fimmu.2022.931831
        • Ahmed I.
        • Karedath T.
        • Al-Dasim F.M.
        • Malek J.A.
        Identification of human genetic variants controlling circular RNA expression.
        RNA. 2019; 25: 1765-1778https://doi.org/10.1261/rna.071654.119
        • Lappalainen T.
        • Sammeth M.
        • Friedländer M.R.
        • 't Hoen P.A.
        • Monlong J.
        • Rivas M.A.
        • Gonzàlez-Porta M.
        • Kurbatova N.
        • Griebel T.
        • Ferreira P.G.
        • Barann M.
        • Wieland T.
        • Greger L.
        • van Iterson M.
        • Almlöf J.
        • Ribeca P.
        • Pulyakhina I.
        • Esser D.
        • Giger T.
        • Tikhonov A.
        • Sultan M.
        • Bertier G.
        • MacArthur D.G
        • Lek M.
        • Lizano E.
        • Buermans H.P.J
        • Padioleau I.
        • Schwarzmayr T.
        • Karlberg O.
        • Ongen H.
        • Kilpinen H.
        • Beltran S.
        • Gut M.
        • Kahlem K.
        • Amstislavskiy V.
        • Stegle O.
        • Pirinen M.
        • Montgomery S.B.
        • Donnelly P.
        • McCarthy M.I.
        • Flicek P.
        • Strom T.M.
        • Consortium Geuvadis
        • Lehrach H.
        • Schreiber S.
        • Sudbrak R.
        • Carracedo A.
        • Antonarakis S.E.
        • Häsler R.
        • Syvänen A.C.
        • van Ommen G.J.
        • Brazma A.
        • Meitinger T.
        • Rosenstiel P.
        • Guigó R.
        • Gut I.G.
        • Estivill X.
        • Dermitzakis E.T.
        Transcriptome and genome sequencing uncovers functional variation in humans.
        Nature. 2013; 501: 506-511https://doi.org/10.1038/nature12531
        • Liu Z.
        • Ran Y.
        • Tao C.
        • Li S.
        • Chen J.
        • Yang E.
        Detection of circular RNA expression and related quantitative trait loci in the human dorsolateral prefrontal cortex.
        Genome Biol. 2019; 20: 99https://doi.org/10.1186/s13059-019-1701-8
        • Xu T.
        • Wang L.
        • Jia P.
        • Song X.
        • Zhao Z.
        An integrative transcriptomic and methylation approach for identifying differentially expressed circular RNAs associated with DNA methylation change.
        Biomedicines. 2021; 9: 657https://doi.org/10.3390/biomedicines9060657
        • Kular L.
        • Liu Y.
        • Ruhrmann S.
        • Zheleznyakova G.
        • Marabita F.
        • Gomez-Cabrero D.
        • James T.
        • Ewing E.
        • Lindén M.
        • Górnikiewicz B.
        • Aeinehband S.
        • Stridh P.
        • Link J.
        • Andlauer T.F.M.
        • Gasperi C.
        • Wiendl H.
        • Zipp F.
        • Gold R.
        • Tackenberg B.
        • Weber F.
        • Hemmer B.
        • Strauch K.
        • Heilmann-Heimbach S.
        • Rawal R.
        • Schminke U.
        • Schmidt C.O.
        • Kacprowski T.
        • Franke A.
        • Laudes M.
        • Dilthey A.T.
        • Celius E.G.
        • Søndergaard H.B.
        • Tegnér J.
        • Harbo H.F.
        • Oturai A.B.
        • Olafsson S.
        • Eggertsson H.P.
        • Halldorsson B.V.
        • Hjaltason H.
        • Olafsson E.
        • Jonsdottir I.
        • Stefansson K.
        • Olsson T.
        • Piehl F.
        • Ekström T.J.
        • Kockum I.
        • Feinberg A.P.
        • Jagodic M.
        DNA methylation as a mediator of HLA-DRB1*15:01 and a protective variant in multiple sclerosis.
        Nat. Commun. 2018; 9: 2397https://doi.org/10.1038/s41467-018-04732-5