A map of human circular RNAs in clinically relevant tissues

1Cellular circular RNAs (circRNAs) represent a class of single-stranded, unusually stable RNAs originating from 5′-to-3′ transcription of coding gene exons or long non-coding RNAs (lncRNAs) that produce covalently closed head-to-tail (or back-spliced) circularized transcripts. Many circRNAs are tissue- and developmental stage-specifically expressed. The circRNA co-transcriptional splicing can compete with linear splicing events and can depend on the binding of the RNA-binding proteins, MBNL1 or QKI, in intronic sequences. Intronic complementary sequences, inter alia repetitive elements, and the RNA-editing protein, ADAR1, were linked to the circularization of exons. Earlier, the circRNAs CDR1as (ciRS-7) and circSRY were shown to exhibit important functions in sponging miRNAs and thereby functioning as post-transcriptional regulators. circRNAs are resistant to the exonuclease RNase R that solely digests linear transcripts. This feature can be used to validate circRNA candidates by comparing their abundance in the RNase R-treated and untreated samples. circRNAs found in clinical specimens, like blood, reveal that these abundant transcripts could serve as biomarkers. Here, we present a circRNA resource catalog to supplement existing databases with new circRNA transcripts in human cell types. By identifying and validating selected circRNAs in these human tissues relevant to clinical research, we provide multiple examples of abundant and highly tissue-specific circRNA expression in host genes that have been associated with pathogenesis of human disease.

2After approval by the ethics committee (Charité Medical Faculty Berlin and University Clinic Erlangen) and written, informed consent, we obtained human tissues. Mesenchymal stromal cells (MSCs) from a non-affected healthy female (23 years) donor were obtained, characterized, and differentiated as previously described. Fibroblasts from buttocks biopsies of a non-affected male donor (25 years) were cultivated until passage six in M-199, supplemented with 20% FCS. Single samples, each from one donor, were used for sequencing.

3Patients or patients’ parents signed informed consent on anonymous data collection for research studies conducted on biological samples of patients with primary immunodeficiencies (three Wiskott-Aldrich syndrome samples, four ADA-SCID samples) at San Raffaele Hospital (TIGET02), approved by the San Raffaele Scientific Institute’s Ethical Committee. Four T cell lines were generated from peripheral blood mononuclear cells purified by density gradient centrifugation on Ficoll-Hypaque (Nycomed Pharma, Oslo, Norway) and expanded.

4Adipose tissue was extracted during lipo-aspiration of MSCs from upper abdominal fat. The tissue was rinsed with PBS on Teflon fleece to wash out erythrocytes. Fat spheres were subsequently frozen in liquid nitrogen. Neutrophils were extracted from peripheral whole blood that was supplemented with 30% of dextran. Cells settled down in a syringe after 30 min. The upper phase was under-laid with histopaque 1083. After centrifugation at 4 °C for 30 min at 1200 rpm, speed was slowed down to 1050 rpm after 15 min for another 15 min. Pelleted cells were resuspended in 10 ml water for water lysis. For neutralization, 3.33 ml of 3.6% NaCl was added for 10 min. After 10 min centrifugation at 1050 rpm, pelleted neutrophils were resuspended in TRIzol® Reagent (Ambion). Plasma, serum, and platelets were prepared using the Vacutainer system. Whole blood in serum tubes was left undisturbed for clot formation. After 15–30 min, the clot was centrifuged at 1000×g for 10 min and the supernatant serum was immediately frozen at −80 °C. Plasma was prepared from whole blood in EDTA tubes. After a centrifugation at 2000×g for 10 min, the supernatant was frozen. Citrate tubes were used to obtain platelets. The whole blood was centrifuged for 15 min at 100×g without rotor break, preventing platelets’ activation. Endothelial progenitor cells (EPCs) were extracted from umbilical cord blood and expanded in vitro using the Lonza EGM™-2 kit. Human umbilical vein cells (HUVECs) were freshly prepared using standard techniques from umbilical cord and cultivated until passage four in EGM medium (Lonza). Adipose tissue, cortex, placenta, decidua, heart, vena cava, muscle, and umbilical cord were minced using pistils and homogenized with matrix beads in the MP FastPrep-24 Tissue and Cell Homogenizer.

5RNA was prepared using TRIzol® Reagent (Ambion) and phenol/chloroform precipitation. For Illumina sequencing, rRNA depletion was done with the RiboMinus™ eukaryote kit according to manufacturer’s recommendation (Life Technologies). Bioanalyzer measurement validated the successful rRNA depletion. The Illumina TruSeq sample preparation kit (v2) was used to generate the libraries for sequencing. For qRT-PCR, total RNA of the identical samples that were used for RNA-seq was digested with RNase R (3 U/μg RNA, Epicenter Technologies) and incubated for 15 min at 37 °C with following inactivation for 3 min at 95 °C. To reach similarly effective RNase R treatment, all samples were treated simultaneously in one approach. Then, the RNA was spiked with 10% of Caenorhabditis elegans total RNA. After phenol/chloroform precipitation, the RNA was reverse transcribed using RevertAid first strand cDNA synthesis kit (Fermentas) or Maxima RT kit (ThermoFisher Scientific), and SYBR-green quantification (Roche) was performed according to standard protocols on ABI 7500 or StepOnePlus (ThermoScientific). Oligonucleotides flanking the circRNA head-to-tail junctions were designed in Primer3 (v. 0.4.0). RNase R assays were normalized to C. elegans eif3d spiked-in RNA and to human GAPDH or Vinculin. For experiments on WAS and ADA-SCID samples, ΔCt was calculated compared to 28S rRNA. In general, expression was quantified applying the ΔΔCt method. qRT-PCR products were analyzed for amplicon size, specificity, and integrity on 3% agarose gels; concatemers were not taken into account. Sanger sequencing of qRT-PCR products was performed using the Big Dye® Terminator Cycle Sequencing on the 3130xl Genetic Analyzer (ABI) using Gene Mapper® Software Version 4.0. Kit v1.1 (ABI). SeqMan software (Lasergene Version 7.0; DNAStar) was used to analyze the traces.

6circRNAs were detected and annotated using the Memczak et al. (2013) pipeline. Human genome reference used for all analyses was hg19 (February 2009, GRCh37), downloaded from UCSC. Upon detection, candidates were annotated using RefSeq and GENCODE v17 gene models. Table S2 and circBase summarize the detected circRNAs across all cell types. Table S2 harbors genomic positions and annotated host transcripts, sense or antisense strand orientation, circBase IDs, genomic and spliced lengths, number of sequencing reads supporting a head-to-tail junction, as well as the number of either 5′ or 3′ linear reads for each circRNA candidate. We also calculated the circular-to-linear ratios and added the list of samples from other studies listed within circBase.

7RNA sequencing data have been deposited in the Gene Expression Omnibus (GEO) under the accession number GSE100242.

8The ratio of circular and linear isoforms (circular-to-linear ratio, CLR) was calculated as described in. For each circRNA candidate, we counted the number of reads overlapping the head-to-tail junction, and the number of reads spliced linearly over the 5′ and 3′ sites that gave rise to a circRNA. CLR was expressed as the number of reads spanning the head-to-tail junction divided by the number of linear reads mapped over the splice site (5′ or 3′) with the higher read count:

9To avoid division by zero when calculating CLR, a pseudocount of 1 was used where no linearly spliced reads were detected. To estimate the expression levels of circRNA host genes, we mapped RNA-seq libraries to the hg19 reference using STAR and counted the reads mapped to Ensembl (release 75) gene models using the htseq-count tool.

10circRNAs in Figs. 1c and 2a, c were sorted into three discrete expression classes: (i) “high”—top 10% (5% for platelets) expressed, measured by raw read counts, circRNAs in a particular sample; (ii) “detected”—all circRNAs that satisfied the minimum expression threshold of two unique reads overlapping head-to-tail junction; and (iii) “not detected”—circRNAs that were not detected in a particular sample. Only circRNAs that were assigned to “high” category in at least one of the compared samples were plotted.

11Differential gene expression analysis was performed using the DESeq package. method=“blind” and sharingMode=“fit-only” options were used when running the estimateDispersions function, as suggested by the package documentation for experimental designs with no biological replicates.

12Peptides were generated using an automated setup. Briefly, samples were reduced with 1 mM tris(2-carboxyethyl) phosphine (TCEP) and free sulfhydryl groups carbamidomethylated using 5.5 mM chloroacetamide. Proteins were pre-digested with 0.5 μg sequencing grade endopeptidase LysC (Wako) for 3 h at room temperature and subsequently diluted with four volumes of 50 mM ammonium-bicarbonate (ABC). Tryptic digestion occurred for 10 h at room temperature using 1 μg sequencing grade trypsin (Promega). The reaction was stopped by adding trifluoroacetic acid (TFA) to a final concentration of 1% resulting in a final pH of 2. The peptides were purified by using C18 stage-tips (3 M).

13Peptides eluted from C18 stage-tips were run on an LC-MS setup. The fractionated and unfractionated samples were measured by LC-MS/MS on a Q Exactive orbitrap mass spectrometer (Thermo) connected to a Proxeon nano-LC system (Thermo) in data-dependent acquisition mode using the top 10 peaks for HCD fragmentation. Peptides were separated on an in-house prepared nano-LC column (0.074 mm × 250 mm, 3 μm Reprosil C18, Dr. Maisch GmbH). Five microlitres of the sample were injected and the peptides were eluted on a 3-h gradient of 4 to 76% ACN and 0.1% FA in water at flow rates of 0.25 μl/min. MS acquisition was performed at a resolution of 70,000 in the scan range from 300 to 1700 m/z, MS2 spectra were collected at a mass resolution of 17,000 with a fixed injection time of 120 ms. Dynamic exclusion was set to 30 s and the normalized collision energy was specified to 26. The eluent was directly sprayed into an Q Exactive mass spectrometer (Thermo Fisher Scientific) equipped with a nano electrospray ion source. The recorded spectra were analyzed using MaxQuant software package version 1.5.2.4, with an Andromeda search using the combined UniProt Homo sapiens and Oryctolagus cuniculus databases and a custom database for the circRNA-derived peptides with a false discovery rate of 1% (peptides and proteins). The fixed and variable modifications were set to carbamidomethylation of cysteines and methionine oxidation, respectively. For further data analysis, the R statistical software package was used (supplement: python script).

14We generated a circRNA resource catalog for various research interests by sequencing ribosomal RNA-depleted total RNA (Supplemental Table S1). The head-to-tail splice junction identification and the sequence analysis were done according to previously published protocols (Supplemental Fig. S1),. The circRNA catalog can be retrieved in Supplemental Table S2 or circBase, http://www.circbase.org/. We selected circRNA candidates for validation according to either of the following criteria: (i) the candidate originated from disease or developmental genes (Table 1); (ii) the linear host transcripts have been proposed as biomarkers; (iii) the circRNA was not present in circBase; (iv) the circRNA was encoded from a lncRNA; or (v) the circRNA showed extraordinary genomic length or expression determined via the read count. In total, we selected 112 candidates, of which we validated 71 circRNAs by RNase R assays (validation rate 63.4%, see “Methods”, Supplemental Fig. S2a–p). We normalized expression values to C. elegans eif3d spiked-in RNA and to human GAPDH or Vinculin. Concatamers in putative circRNA candidates were not taken into account.

15Of the 5225 circRNAs, 35.9% (1878 circRNAs) were new compared to circBase. circRNAs (3841) were unique for the investigated cell types (Fig. 1a). We detected 82.9% circRNAs in coding genes (exons, 5′ + 3′ UTRs), 2.2% antisense transcripts, 5.4% intron-derived circRNAs, 6% in non-coding genes, and 1.1% from intergenic regions (Fig. 1b).

16First, we analyzed circRNA expression during MSC differentiation. MSCs were differentiated into proliferating chondrocytes, osteocytes, and vascular smooth muscle cells (VSMCs). Flow cytometry revealed CD105+, CD90+, CD73+, HLA-ABC+, CD31−, CD34−, CD45−, and HLA-DR− cells and the multi-lineage potential validated MSC properties. In MSCs, we detected 55 circRNAs, in contrast to 148 in MSC-derived chondrocytes, 104 in osteocytes, or 137 in VSMCs. In chondrocytes, we validated a circRNA deriving from PLOD2, a gene which forms collagen crosslinks and was differentially expressed in a model of osteoarthritis-related synovial fibrosis. During the abdominal fat aspiration to obtain MSCs, we additionally harvested adipose tissue that harbors the MSC niches. In the abdominal fat of the same MSC donor, we identified 507 circRNAs. We confirmed a circRNA in SORBS1, a gene inhibiting the induction of glucose transport by insulin, and two circRNAs in PLIN4, a gene stimulating lipolysis in adipocytes (Supplemental Fig. S2a–d). When comparing different MSC-derived tissues, we commonly observed different circRNA isoforms spliced from the same host genes in different MSC-derived cells. Solely eight circRNAs overlapped between adipose tissue, MSCs, and their derived cells (Supplemental Table S2). In eight different tissues of one healthy donor, we excluded interindividual differences and found tissue-restricted expression patterns (Fig. 1a). The absence of circRNA housekeeper and direct comparisons of circRNA expressions between different tissues are controversially discussed endeavors. Thus, we selected the top 10% circRNA candidates within MSCs and MSC-derived cells and clustered 31 circRNAs candidates based on their expression levels: “not detected,” “detected” with at least two reads, or “highly expressed” when detected within the top 10% of the candidates (Fig. 1c). We observed ubiquitous expression for some circRNAs and differential expression regarding the MSC-derived cells. We validated selected circRNA candidates within the cluster analysis showing RNase R resistance and confirmed the results of clustering (Fig. 1d and Supplemental Fig. 2q).

17Next, we compared physiologically neighboring tissues: umbilical cord, endothelial progenitor cells (EPCs), HUVECs, and HUAECs. Of the 211 circRNAs in HUAECs, only two overlapped with the other tissues. EFEMP1, an angiogenesis promoter, poor prognostic marker in ovarian cancer, and potential therapeutic target in glioblastoma treatment, harbored a circRNA in HUVECs (22 circRNAs). As previously described, we analyzed the top 10% of the circRNA candidates. A pool of 43 circRNAs showed differential expression and we validated two circRNAs (Fig. 2a, b and Supplemental Fig. 2q).

18In cerebral cortex (339 circRNAs), we validated circRNAs spliced from important cerebral genes. ERC2 is involved in neurotransmitter release and expresses a circRNA. Five circRNA isoforms were found in ATRNL1, a gene regulating the energy homeostasis by melanocortins in the hippocampus. Another circRNA is hosted by NTRK2, a gene that was associated with synaptic dysfunction in Huntington’s disease and neuronal differentiation and plasticity in hippocampus. We validated a circRNA in RTN4, a gene inhibiting axonal sprouting and modulating Alzheimer’s disease progression in a mouse model and one of two predicted circRNAs in HOMER1, a gene that is involved in synaptic activity and various neurological disorders. Three circRNAs were expressed in ATXN10 that maintains the survival of neurons, studied in the spinocerebellar ataxia type 10 (Supplemental Fig. S2g).

19In placenta, we detected 63 circRNAs, in comparison to 173 in decidua; 15 circRNAs overlapped. We confirmed a circRNA in PSG5, a gene encoding a pregnancy-specific glycoprotein. Low levels indicate pregnancy complications. Severe early onset preeclampsia, fetal growth restriction, and HELLP syndrome were associated with high expression of PAPPA2 that encoded two circRNAs. ALPP, less expressed and active in hyperglycemic and diabetic placentas of pregnant women infected with or without Trypanosoma cruzi, harbored a circRNA (Supplemental Fig. S2h).

20We also obtained right atrial tissue and vena cava from two children with multiple cardiac defects. Atrium (340 circRNAs) and vena cava (702) had an intersection of 115 circRNAs; 51 overlapped with calf muscle. In atrium, we validated a circRNA in NPPA, a gene that was associated with development of type 2 diabetes. NPPA is reactivated in response to cardiovascular disorders and converted to its active form by CORIN, which harbored four circRNAs (validated in atrium and vena cava) and could be potential biomarkers for heart failure. Atrial fibrillation was linked to the dysfunction of RYR2, which also produces seven circRNA isoforms in atrium and ten in vena cava. Mutations in myosin heavy chains cause hypertrophic cardiomyopathies. In myosin MYH6, we validated one circRNA. QKI is involved in circRNA biogenesis, is responsible for cardiovascular development, and encodes two circRNAs in atrium and vena cava. Alterations in the regulation and expression of SLC8A1 (two circRNAs in atrium) contribute to various cardiovascular symptoms (Supplemental Fig. S2i–k).

21Platelets expressed 3324 circRNAs. Platelets derive from bone marrow megakaryocytes, lack nuclei and highly abundant mRNA reservoirs, although translational capabilities are intact. Previously, high circRNA expression and circRNA properties were described in platelets and this enrichment was associated with transcriptome degradation. We found in our data that platelets harbored much more abundantly expressed circRNAs than any other tissue. For example, circRNA expression in ACVR2A and SMARCA5 was extremely high compared to mRNA (Supplemental Fig. S3). We validated circRNAs in the phosphodiesterases PDE3A, PDE4D, and PDE5A. PDEs hydrolyze cAMP and cGMP to control blood vessel relaxation, cardiac contractility, and inhibition of platelet aggregation. PDE3A was previously associated with Mendelian hypertension. Moreover, the guanylate cyclase GUCY1B3 converting GTP into cGMP expressed a circRNA (Supplemental Fig. 2l–o).

22CircRNA expression in plasma (57), serum (39), suggested that circRNAs could be secreted, as it was shown earlier for micro- and other RNAs, and indicated by circRNAs identified in cell culture or serum exosomes.

23In neutrophils (274 circRNAs), TLR6 functions in the innate immune response and harbored a circRNA. Another key component in the immune system expressing a circRNA in neutrophils is MYO1F, a class I myosin regulating the host defense against infection. No overlap between plasma, serum, neutrophils, and platelets facilitates the idea of tissue-restricted circRNA expression. For clustering, we used all plasma and serum circRNAs, the top 10% of neutrophils and the top 5% of platelets and validated four candidates (Fig. 2c, d and Supplemental Fig. 2q).

24Due to the lack of nuclei and the highly abundant circRNAs in platelets, we hypothesized that circRNAs could serve as templates for translation as recently shown for few circRNA examples in human and fly. Thus, we used RNase R-treated whole platelet RNA to perform in vitro translation experiments followed by highly sensitive mass spectrometry. We derived putative open reading frames (ORF) that span head-to-tail junctions of our circRNA candidates. These predictions were compared to mass-spectrometrically detected peptides. Controls were reticulocyte lysate of the in vitro translation kit, non-RNase R-treated whole platelet RNA, and total protein of the same platelet-donor. Although we detected peptides in the RNase R-treated translated sample and in the cell lysate matching circRNAs in platelets, those candidates did not overlap with head-to-tail junctions (Supplemental Table S3).

25We also investigated circRNA isoforms, since we observed around 100 genes hosting more than five different circRNA isoforms. For example, we detected 18 circRNA isoforms derived either from PTPN12 or TTN in platelets, atrium, and vena cava (Supplemental Fig. S4). We compared circRNA expression directly to linear transcript expression, by counting linearly spliced and head-to-tail spliced reads. The number of reads overlapping with the head-to-tail splice junctions was divided by the number of linear splicing events with identical splice sites (Supplemental Table S2). The calculated value was plotted against the transcript copies per million transcripts (TPM) to describe circRNA expression as circular-to-linear ratio (Fig. 2e). Collectively, we detected high circular-to-linear expression ratios in tissues with abundant circRNA expression, e.g., a platelet circRNA in SMARCA5 had a circular-to-linear ratio of 151:1 (Supplemental Fig. S3b).

26Finally, we demonstrate differential circRNA expression in ADA-SCID and WAS, two primary immunodeficiencies which are caused by mutations in ADA or WAS, respectively. First, we compared linear transcripts from one patient compared to a non-affected control. We detected significantly differential expression (p ≤ 0.05) of 79 mRNAs in ADA-SCID and 19 mRNAs in WAS lymphoblastoid cells (LCLs) (Supplemental Table S4 and Supplemental Fig. S5a, b). The results were consistent with the molecular pathogenesis of both disease phenotypes. For example, upregulated BANK1 (p = 2.30 × 10−3, log2-fold change (lfc) = 3.8) or PBXIP1 (p = 3.53 × 10−2, lfc = 2.4) mRNAs in ADA-SCID were associated with impaired B cell receptor-induced calcium mobilization or early blocking of B cell development in the bone marrow (Supplemental Table S4).

27We next asked whether these differentially expressed linear transcripts harbor also circRNAs with differential expression between patients and controls. We found a circRNA in ROBO1, a gene upregulated in ADA-SCID (mRNA: p = 9.76 × 10−6, lfc = 5.3) and WAS (mRNA: p = 3.68 × 10−4, lfc = 8.31, Supplemental Table S4). Moreover, CDC42BPA expressed an upregulated circRNA in ADA-SCID (mRNA: p = 3.46 × 10−3, lfc = 4.7) and WAS (mRNA: p = 1.93 × 10−3, lfc = 7.4). Notably, ROBO1 and CDC42 are linked to the pathogenesis of WAS. Slit-2 and Robo-1 complexes have been described to inhibit the CXCR4/CXCL12-mediated chemotaxis of T cells. Moreover, ROBO1 and ROBO4 bind WAS to induce filopodia formation. Cdc42-dependent WAS activation was also reported. CDC42 is a major regulator of podosome formation and remodels actin during B cell signaling, whereas CDC42BPA is a downstream effector of CDC42. B cell signaling is impaired both in ADA-SCID and WAS. In ADA-SCID, we found a circRNA in TNFRSF11A (mRNA: p = 8.28 × 10−3, lfc = 3.3) TNF receptors participate in several pathways altered in ADA-SCID. We first validated the circRNAs in ROBO1, CDC42BPA, and TNFRSF11A (Fig. 2f and Supplemental S6a) and tested next their differential expression in three WAS and four ADA-SCID samples, compared to four non-affected LCL samples (Fig. 2g). circRNA expression of phenotypically relevant genes was higher in the disease samples (Fig. 2g).

28Collectively, we provide a circRNA catalog of human tissues relevant to various fields of clinical research. We provide evidence that circRNAs could serve as biomarkers and that circRNA expression profiles could be directly linked to clinically apparent phenotypes. We focused on detecting circRNAs in various single samples. For further analyzing the proposed circRNA candidates as suitable biomarkers, broader studies addressing tissue specificity vs. donor specificity are needed. Our data corroborate recent findings that circRNA expression is highly tissue-specific. We did not find evidence that platelet circRNAs were translated; however, our result does not provide conclusive evidence that circRNAs are not translated, as it highly depends on mass-spec sensitivity. As previously suggested, a resistance to RNA degradation can explain the high abundance of circRNAs in platelets. In the absence of transcription, the detected circRNAs could function independently of transcriptional regulation.

29As discussed previously, circRNAs uncovered in this study could contribute to regulatory networks governing coding gene expression by acting as miRNA target decoys, RNA-binding protein (RBP) sponges, scaffolding molecules, and transcriptional regulators. A circRNA function is further supported by the conserved nature of circRNA expression and the tissue-specific and regulated abundance. Although we can only speculate that currently disclosed circRNAs influence the functions of their linear counterparts, these new isoforms need to be considered when investigating disease-relevant genes. Since circRNA biogenesis can compete with pre-mRNA splicing, this opens up the possibility that the mRNA output from those, oftentimes well studied genes, is controlled by the hitherto unknown circRNA.