The Landscape of Circular RNA Expression in the Human Brain

Journal Pre-proof The Landscape Of Circular RNA Expression In The Human Brain Akira Gokool, Firoz Anwar, Irina Voineagu PII:

S0006-3223(19)31581-1

DOI:

https://doi.org/10.1016/j.biopsych.2019.07.029

Reference:

BPS 13944

To appear in:

Biological Psychiatry

Received Date: 24 January 2019 Revised Date:

19 June 2019

Accepted Date: 9 July 2019

Please cite this article as: Gokool A., Anwar F. & Voineagu I., The Landscape Of Circular RNA Expression In The Human Brain, Biological Psychiatry (2019), doi: https://doi.org/10.1016/ j.biopsych.2019.07.029. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc on behalf of Society of Biological Psychiatry.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

The Landscape Of Circular RNA Expression In The Human Brain Akira Gokool1, Firoz Anwar1, Irina Voineagu1,*

1. School of Biotechnology and Biomolecular Sciences, University of New South Wales, Sydney, Australia ∗ Corresponding author: Assoc. Prof. Irina Voineagu School of Biotechnology and Biomolecular Sciences University of New South Wales Kensington, Sydney NSW 2052 Australia Phone: +61 (02) 9385 2029 Email: [email protected] Short title: Landscape Of Circular RNA Expression In The Human Brain Keywords: Human brain; Circular RNA; Transcriptome; Autism spectrum disorder; Splicing; Gene regulation

1

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69

ABSTRACT Background Circular RNAs (circRNAs) are enriched in the mammalian brain and are upregulated in response to neuronal differentiation and depolarisation. These RNA molecules, formed by non-canonical backsplicing, have both regulatory and translational potential. Methods Here, we carried out an extensive characterisation of circRNA expression in the human brain, in nearly two hundred human brain samples, from both healthy individuals and autism cases. Results We identify hundreds of novel circRNAs and demonstrate that circRNAs are not expressed stochastically, but rather as major isoforms. We characterise inter-individual variability of circRNA expression in the human brain and show that inter-individual variability is less pronounced than variability between cerebral cortex and cerebellum. Finally, we identify a circRNA co-expression module upregulated in autism samples, thereby adding another layer of complexity to the transcriptome changes observed in autism brain. Conclusions These data provide a comprehensive catalogue of circRNAs as well as a deeper insight into their expression in the human brain, and are available as a free resource in browsable format at: http://www.voineagulab.unsw.edu.au/circ_rna INTRODUCTION Circular RNAs (circRNAs) are an emerging class of RNAs formed by the non-sequential back-splicing of pre-messenger RNAs (1). CircRNAs are expressed in a tissue-specific manner and in some cases are more efficiently generated than their linear cognate mRNAs (2, 3). Although circRNAs are considered to be primarily non-coding molecules, a subset of circRNAs can be translated (4). Due to their circular nature, circRNAs lack a 5’ cap structure and a poly-A tail, which in turn renders them resistant to enzymatic degradation. Consequently, circRNAs are exceptionally stable. Several studies have demonstrated that across a range of mammalian tissues, circRNA expression is most abundant in the brain, with an overall enrichment in the cerebellum (5-7). circRNAs also show dynamic expression during neuronal differentiation and depolarisation (6-9), and are highly concentrated in synaptosomes (6), indicating that these molecules likely play a functional role in neurons. The extent to which circRNAs influence the expression levels of their parental genes is yet to be elucidated, but at least two mechanisms are known to be at play: they can function as miRNA sponges (10, 11), and can increase the transcriptional efficiency of their parental genes (12). The former mechanism is important for pluripotency and differentiation (13), as well as astrocyte activation (14). Remarkably, the first study to investigate circRNA loss-of-function in vivo revealed that interactions between circRNAs and miRNAs are important for normal brain function (11). Cdr1as knock-out mice, which lack circRNA expression at the Cdr1as locus, display impaired sensory-motor gating and abnormal synaptic transmission. These phenotypes are associated with changes in miR-7 and miR-671 levels, and consequent transcriptional changes in immediate early gene expression (11). Despite accumulating evidence for an important role of circRNAs in the brain, our current understanding of their expression in the human brain is limited by relatively low sample sizes of most studies (5, 15-19). Therefore, a robust large-scale dataset of circRNA expression in the human brain across multiple brain regions is currently lacking. Given the limited data on circRNA expression in the human brain, several questions remain open: (a) How does circRNA expression vary across individuals and human brain regions? (b) Does circRNA expression change significantly in an age-dependent manner? (c) Given that the brain is characterised by extensive alternative splicing (AS), is the abundant circRNA expression observed in the brain primarily a reflection of the complexity of brain AS events? Here, we begin 2

70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120

to address these questions by investigating circRNA expression in nearly two hundred human brain samples, from three brain regions: frontal cortex, temporal cortex and cerebellum, from normal individuals and autism spectrum disorder (ASD) cases. METHODS RNA-seq data The RNA-seq data published by Parikshak et al. (20) (ribo-depleted, unstranded, 50 bp paired-end), was kindly provided by Prof. Daniel Geschwind before publication. Frontal cortex samples: Brodmann area (BA) 9; temporal cortex: BA 41/42 and 22, and cerebellum: vermis. Only samples with RIN> 5 were included (Supplemental Table S1). Since 3’ bias of the RNAsequencing data better reflects RNA degradation (21), median 3’ bias values (reported in Parikshak et al. (20)) were used as a covariate. The distribution of samples between the two datasets (DS1 and DS2) is listed in Supplemental Table S1, and was based on the order in which we obtained the data, which in turn was based on the order of RNA-seq data generation that randomised across groups (Figure 1a). The RNA-seq data generated for the benchmarking analysis is described in Supplemental Methods. The RNA-seq data from human brain organoids was downloaded from SRA, accession number GSE99951. RNA-seq data analysis is described in detail in Supplemental Methods. RT-PCR validation CircRNAs were randomly selected for validation among those present in more than a third of the brain tissue samples. 22 circRNAs for which primers could be designed effectively were included in the validation experiments. RT-PCR was performed with divergent primers, on three independent RNA samples. RNA was extracted using a Qiagen miRNeasy kit, with on-column DNase digest. PCR amplification was carried out using BioRad iTaq Polymerase for 35 cycles, on an ABI ViiA7 cycler, followed by Sanger sequencing. Primer sequences for RT-PCR experiments are listed in Supplemental Table S2. Good quality Sanger sequencing of the PCR product was obtained for 15 of the circRNAs tested, and 13 of these confirmed the back-spliced junction (Supplemental Figure S1). RESULTS Dataset Overview We investigated circRNA expression in a total of 197 human brain samples (20) from frontal cortex, temporal cortex and cerebellum (Methods; Supplemental Table S1) obtained from control (CTL) and ASD individuals (Figure 1a). These paired-end unstranded RNA-seq data were generated following ribosomal RNA-depletion (Methods; (20)). In order to assess the reproducibility of our findings, samples were divided into a discovery dataset (DS1, 144 samples) and a replication dataset (DS2, 53 samples). To allow adequate statistical power, the discovery dataset DS1 was assigned a larger sample size than the replication dataset DS2, where we would test for effects already identified in DS1. Given the well documented transcriptional similarity between frontal and temporal cortex (22), also observed for circRNAs (see below), we considered these regions as a single group (cerebral cortex; CTX). Within each dataset, there was no significant difference in age or gender ratios between CTX and cerebellum (CB), or between ASD and CTL (Supplemental Figure S2). To assess the quality of the RNA-seq data and of the circRNA quantification approach, we first carried out a benchmarking analysis (Supplemental Information). We compared (a) circRNA expression quantification in ribodepleted unstranded RNA-seq data (DS1, DS2) with polyAselected and ribodepleted stranded RNA-seq data generated from the same samples and (b) two circRNA quantification algorithms: CIRCexplorer (23) and DCC (Detect CircRNAs from Chimeric

3

121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172

reads; (24)). Benchmarking showed that the RNA-seq data were appropriate for circRNA detection and that DCC performed the best for circRNA quantification (Supplemental Information). CircRNA quantification in all 197 samples (including control and ASD samples) lead to the detection of a total of 43,872 circRNAs in DS1 and 28,251 circRNAs in DS2 (Supplemental Methods). Circular junction (i.e. back-splice) reads were normalised to total library size as counts per million (CPM). circRNA expression was also normalised to that of the parental transcript by calculating a circular-to-linear ratio (CLR), as well as a circularization index (CI; Figure 1b). CircRNAs were considered robustly expressed if they were detected by at least 2 circular junction reads per sample in a minimum of 5 distinct samples. The choice of filtering criteria was based on reproducibility between DS1 and DS2 (Figure 1c; Supplemental Information). The filtered datasets (DS1: 14,386 circRNAs, DS2: 9,440 circRNAs) were used for all downstream analyses. We assessed the circular nature of a subset of 22 circRNAs by RT-PCR with divergent primers (which would amplify on a circular but not a linear molecule), and found a 90.9% validation rate (Methods; Figure 2b). We also observed a strong correlation (Spearman rho=0.93) between the mean circRNA expression level (CLR) in DS1 and DS2, supporting the robustness of these data. The 14,386 DS1 circRNAs and 9,440 DS2 circRNAs (Supplemental Table S3) were expressed from a total of 4,555 and 3,650 genes, respectively. Gene ontology enrichment analysis of circRNA genes, with length correction (Supplemental Methods), showed an over-representation of genes functioning at the synapse, particularly those involved in synapse vesicle transport and localisation, but also genes involved in cell cycle G2/M phase transition, chromatin organisation and gene silencing (Figure 2c and Supplemental Table S4). Pathway enrichment analysis (Supplemental Methods; Supplemental Table S4), highlighted endocytosis, axon guidance and long-term potentiation as the top three enriched pathways. We further investigated whether circRNA genes are enriched among marker genes for specific neuronal subtypes. To this end, we used marker genes from human brain single-nucleus RNA-seq data; Lake et al. (25) (Supplemental Methods). CircRNA genes were eniched among markers of 5 inhibitory neuron subtypes, including Purkinje neurons, and 2 excitatory neuron subtypes, including granule cells (Supplemental Table S4). We also assessed whether circRNA expressing genes were overrepresented among genes associated with neurodevelopmental disorders: ASD, intellectual disability (ID) and schizophrenia (SCZ); Supplemental Methods. circRNA genes were overrepresented among all three disease categories (hypergeometric test, Bonferroni corrected p-values: 3.15 e-78 for ASD, 1.24 e-32 for ID, 3.64 e-06 for SCZ). The proportion of circRNAs for which the strict orthologous mouse coordinates were also detected as circRNAs in our human brain dataset was 8.6% in DS1 and 9.6% in DS2, consistent with previous observations (6). Identification of novel circRNAs We assessed whether our datasets contained novel circRNA compared to those curated from previous studies in circBase (26). We detected 1,548 novel circRNAs in DS1 and 692 in DS2, of which 83% were detected in both datasets (Supplemental Methods; Figure 2a). Hundreds of novel circRNAs were detected in more than 10 brain samples, and some were expressed in more than a hundred samples, demonstrating that they are frequently expressed in the human brain (Figure 2d). CircRNA annotation relative to genomic features showed that most circRNAs, including novel circRNAs, were formed between annotated exon-exon junctions (Supplemental Figure S3). The novel circRNAs identified included both novel isoforms from known circRNAproducing genes, as well several hundred circRNAs expressed from genes not previously reported to circularize. Throughout this manuscript, we use the term “circRNA isoforms” to refer to circRNAs produced by back-splicing of distinct exon-exon junctions of a gene. As an example, we outline circRNA expression from RIMS2, a gene that encodes a presynaptic protein involved in regulating synaptic membrane exocytosis (27). RIMS2 shows 4

173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223

conserved circRNA expression in human, mouse, and pig brain (6, 28), with 35 known isoforms of circRIMS2 in human brain (26). Here, we identify 15 RIMS2 circRNA isoforms detected in the human brain in both DS1 and DS2, of which 7 are novel circRNAs (Figure 2e). Notably, 3 of the novel circRNAs are highly expressed in both datasets (> 0.1 CPM in at least 2 distinct samples), and 2 of the novel isoforms were expressed in more than 50 brain samples. In addition to the high complexity of circRIMS2 isoform expression, we also find that RIMS2 shows circRNA isoform switching between cerebral cortex and cerebellum (Figure 2e). CircRNA expression is characterised by major isoform(s) To investigate the global properties of circRNA expression in the human brain, we first used the data from control samples (N=68 in DS1; N=29 in DS2). An important layer of regulation of mRNA expression consists of major isoform(s), which account for most of the transcriptional output of a gene in a given tissue or cell type. We thus asked whether circRNA isoforms are stochastically expressed, or show evidence of major isoform expression. The latter scenario would strongly indicate regulated circRNA expression, in a manner that is not confounded by the difference in circRNA vs. mRNA degradation rates (Figure 3). We assessed the major circRNA isoform relative CLR (i.e. the expression level of the most highly expressed circRNA relative to the total circRNA expression from a given gene, after correction for linear RNA expression levels; Figure 3b). We asked whether the observed major isoform relative CLR was higher than expected by chance. To this end, we randomly sampled back-splice junction reads and linear junction reads from each brain tissue sample 1,000 times, following a negative binomial distribution. Simulated back-splice and linear junction expression levels for each sample were calculated as an average of random sampling values, and used to calculate simulated CLRs. We classified genes by the number of circRNA isoforms expressed, and compared the observed and simulated distribution of mean major isoform relative CLR within each class. We found that the major circRNA isoform accounted for significantly more of the circRNA output than expected by chance, in both DS1 and DS2, demonstrating that circRNAs were not expressed stochastically (Figure 3c-d). We observed good agreement in major isoform calls between DS1 and DS2. 1,910 major isoforms were identified in DS1 (i.e. circRNAs called major isoform in at least 50% of the DS1 samples), of which 1,577 were expressed in DS2. Of these, 89% were called major isoforms in DS2 (Supplemental Figure S4). We then investigated whether major isoform expression is explained by sequence complementarity of flanking introns. Reverse-complementary sequence matches (RCM) of flanking introns were calculated for all circRNAs using autoBLAST ((29); Supplemental Methods). We found that circRNA major isoforms were only marginally more likely than expected by chance to have the highest sequence complementarity score (Supplemental Figure S5). Interplay between canonical- and back-splicing in the human brain We characterised alternative splicing events in the larger dataset (DS1, control samples) using rMATS (30), and contrasting CTX and CB (Supplemental Methods). Cassette exon (i.e. exon skipping) was the predominant AS event (72%), and thus we focused on the comparison of alternatively spliced cassette exons (referred to as “AS” from here on) and circRNA formation. The percentage of AS exons was significantly higher among circRNA-forming exons than among non-circ forming exons after correction for gene expression levels and intron length (Figure 4a). We then compared exon inclusion levels (i.e. the ratio between the inclusion and the skipping isoform of a given AS exon; Supplemental Methods) between circ-exons and non-circ-exons that undergo alternative splicing. Interestingly, circ-exons were formed primarily from exons with high inclusion rates (Figure 4b), indicating that circRNA formation is not primarily a by-product of exon skipping.

5

224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275

The number of circRNAs expressed per gene varied between one and sixty (DS1), and was highly correlated between DS1 and DS2 (Spearman rho = 0.82). Using circRNAs expressed in both DS1 and DS2, we found that 446 genes expressed more than 5 circRNAs. For these “circRNA hotspot genes” we found a significant correlation between the number of circRNAs expressed and the number of alternative splicing events detected per gene (rho=0.32, p< 2.2e-16). The correlation remained significant after correction for the total number of exons per gene (rho=0.14, p< 2.2e-16). Some of the major hotspot genes (Figure 4c), including RIMS1 (a schizophrenia associated gene (31)) and NRXN1 (involved in autism (32, 33)) were indeed characterised by extensive alternative splicing (Figure 4c). However, hotspots of circRNA expression were also observed for 18 genes showing no detectable AS event (Figure 4c). Inter-individual variability of circRNA expression is less pronounced than inter-region variability Gene expression data commonly clusters by broad brain region (such as CTX and CB) rather than by individual (22, 34). To begin to understand the inter-individual variability of circRNA expression in the human brain, we first compared the relationship between mean expression and variance for canonical splice junctions and circRNAs (i.e. back-splice junctions) using the DS1 control samples (N=68). We observed a higher coefficient of variance (CV) for circRNA backsplice junctions compared to canonical splice junctions (mean CV: 3.6 and 3.08 respectively; p < 2.2e-16, Wilcoxon rank sum test; Supplemental Figure S6a). This difference is explained by the fact that while most splice junctions are detected in nearly all samples, circRNAs are often expressed in a small proportion of samples (Supplemental Figure S6b). This observation is consistent with (a) the fact that the efficiency of canonical splicing is higher than that of backsplicing, and consequently back-splicing occurs less frequently (8) and (b) the well documented property of non-coding RNAs to be specifically rather than broadly expressed (35). In addition, over-dispersion of count-level data is also a likely contributing factor. Since circRNA formation is a more rare event than the transcription of the parental transcript, we then investigated whether its frequency was consistent across the two datasets. We found high agreement between the proportion of samples in which a circRNA was detected in DS1 and DS2, as well as high correlation between CLRs in the two datasets (Figure 5a-b). CircRNA expression data clustered by brain region, even after normalisation of circRNA levels to the parental transcript (i.e. CLR), suggesting that similarly to mRNA expression, circRNA expression variability between individuals is less pronounced than region-specific differences (CTX and CB, Figure 6a). Within CTX, frontal and temporal cortex samples clustered together, as commonly observed for gene expression data. ASD samples did not show distinct clustering based on CLR, rather followed the overall pattern of clustering based on brain CTX and CB region (Figure 6a). We also calculated Spearman correlation coefficients of circRNA expression either (i) between different individuals within the same brain region: CTX, CB, or (ii) between samples from different brain regions of the same individual (CTX vs. CB). The correlation values were significantly lower between tissues of the same individual than between individuals (Figure 6b). CircRNA expression differences between CTX and CB We assessed circRNA expression differences between control CTX and CB samples using a linear model that included as covariates age, sex, brain bank, median 3’bias and sequencing batch, as well as estimated proportion of neurons (Supplemental Methods). We identified 496 circRNAs for which expression normalised to the linear transcript (i.e. CI) was significantly different between CTX and CB (FDR < 0.05; Supplemental Methods). 205 of the 496 circRNAs were replicated as significant in DS2 (FDR < 0.05; Supplemental Methods; Supplemental Table S5a), with high agreement of directionality (Figure 6c). Interestingly, 97% of these region-specific circRNAs showed higher circularization rate in CB, suggesting the existence of trans-factors that favour circRNA formation in CB. Of 205 circRNAs differentially expressed between CTX and CB, only 10 involved exons differentially spliced between these regions. For comparison, we also report 6

276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326

differential expression results between CTX and CB without correction for cell type proportions (Supplemental Table S5b). We next investigated how circRNA expression varies with age in control CB and CTX samples. We did not identify any circRNAs significantly changing with age, after correction for covariates and multiple testing (Supplemental Methods). The total number of circRNAs expressed showed a weak increasing trend with age in both brain regions, which did not reach statistical significance, using a linear model with correction for covariates (DS1: p=0.06; DS2: p=0.14; Supplemental Figure S7). We also assessed circRNA expression during cellular maturation in human brain organoids, using RNA-seq data from cells immunopanned with either neuronal or astrocyte markers at various time points of organoid development (0-495 days) (36); Supplemental Methods; Supplemental Table S6. Overall, 56% of the circRNAs identified in organoid cultures overlapped with those identified in post-mortem brain. The number of circRNAs expressed (>=0.1 CPM) was variable at early maturation time points (0-150 days), followed by a progressive decrease in mature neurons and astrocytes (>150 days; Supplemental Figure S8, r= -0.64, p=0.02). Integrative analysis of circRNA expression and miRNA expression in the human brain Since a number of circRNAs have been reported to function as miRNA sponges, we took advantage of published miRNA expression data (37) from a subset of 95 brain samples included in DS1 to investigate the interplay between circRNA and miRNA expression in the human brain. We first assessed whether we could observe previously reported functional interactions. It has recently been reported that in mouse brain, the circRNA CDR1as together with a lncRNA (CYRANO) and two miRNAs (mir-671 and mir-7) form a regulatory network (38). CYRANO sponges mir-7, while mir-7 enhances the mir-671-dependent degradation of CDR1as. If this regulatory network were at play in the human brain as well, one would expect a direct correlation between CYRANO and CDR1as. Furthermore, the model predicts that CDR1as should be anti-correlated with mir-7 and mir-671. Remarkably, the human brain data faithfully follows these predictions: CDR1as was significantly correlated with CYRANO (rho=0.63) and it was significantly anti-correlated with mir7-5p (rho=-0.63) and mir-671-5p (rho= -0.57); Figure 7. To begin to identify candidate miRNA-circRNA interactions, we next calculated Spearman correlations between all circRNAs expressed in at least half of the 95 samples and all miRNAs. We then used STARbase (39), which curates Ago CLIP-seq data, to select miRNA-circRNA pairs that show both evidence of direct interaction by CLIP-seq and significant expression correlation in the brain (absolute rho >0.5). We identified 101 such pairs, of which 64 were negatively correlated (as expected for a sponging effect), and 34 were positively correlated (Supplemental Table S7). Co-expression networks identify circRNA expression differences in cerebral cortex in ASD We carried out a co-expression network analysis of circRNA expression (DS1 data, CI) using the cortex samples, after regressing out all covariates except phenotype (Supplemental Methods). To address the problem of sparse data (i.e. most circRNAs being expressed in a small number of samples) we (a) used a similarity measure robust to outliers (biweight midcorrelation (40)), and (b) only included in the network analyses circRNAs expressed in at least a half of the 92 CTX samples (1,280 circRNAs). We identified 5 co-expression modules, of which one (M4) showed significant eigengene differences between ASD and controls, with higher levels in ASD (p=0.02, Wilcoxon rank-sum test, Bonferroni corrected, Figure 8a; Supplemental Methods). M4 contained 140 circRNAs (Supplemental Table S8), and was moderately preserved (41) in DS2 (z-score=2.85; Supplemental Figure S9, Supplemental Methods). Interestingly, one of the hubs of this coexpression module (defined as circRNAs with the strongest correlation with the module eigengene, i.e. highest kME), is a circRNA expressed from ZKSCAN1, which plays a role in cell proliferation and migration in hepatic cells but its role has not yet been investigated in the brain (Figure 8b). M4 also included circHIPK3, a circRNA implicated in cell proliferation and migration through a

7

327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369

miRNA sponge mechanism (Figure 8b) (42, 43). Consistent with previous data on gene expression (20), we did not identify any ASD-associated circRNA co-expression module in the cerebellum. M4 was enriched for markers of inhibitory and excitatory neurons, and for genes localised at the synapse (Supplemental Table S8). We assessed the enrichment of M4 genes for genes associated with ASD, ID and SCZ and found a significant enrichment for ASD and ID genes (p=4.4 e-05 and p=0.026 respectively; hypergeometric test) but not for SCZ genes (p=0.079; hypergeometric test). We also compared M4 with previously reported gene co-expression modules from ASD and controls (20). M4 significantly overlapped M8 from Parikshak et al. (20). However, M8 was not significantly associated with ASD, consistent with distinct regulation of circRNA and linear mRNA expression. DISCUSSION The data presented here represents a large-scale assessment of circRNA expression in the human brain, bringing further insight into an additional layer of brain transcriptome complexity. We focused on identifying circRNAs reproducibly expressed (in at least 5 samples in each dataset), rather than very rare back-splicing events. The analysis of this rich resource revealed several novel insights into circRNA expression in the human brain. We demonstrate that circRNA expression is characterised by major isoforms, adding a novel aspect to the notion that circRNAs are expressed in a regulated manner in the human brain. Although the expression of most circRNAs was restricted to a subset of samples, we observed a remarkably high correlation between DS1 and DS2 circRNA expression levels, suggesting that the circularization rate was an intrinsic property of a given back-splice junction. Consistent with previous data from human fibroblasts (44), we observed a significant association between AS and circRNA expression. The number of AS events per gene explained ~10% of the variance in the number of circRNAs produced, indicating that AS only partially accounted for circRNA formation, and suggested that the interplay between circRNA formation and AS was locus-specific. circRNA expression differences between CTX and CB were pronounced, overriding interindividual variability. 90% of differentially expressed circRNAs showed higher expression in CB, in agreement with a previous observation based on two frontal cortex and two cerebellum replicates (6). The initial observation was interpreted as a reflection of higher proportion of neurons in cerebellum than in cerebral cortex, a known property of this brain region. Here we demonstrate that the enrichment of circRNAs in the cerebellum is independent of cell-type composition. CircRNA-producing genes significantly overlapped genes associated with ASD, SCZ and ID. Since genes associated with neurodevelopmental disorders are known to be enriched for neuronal and synaptic genes (45, 46), and circRNA–producing genes showed a similar functional enrichment, the overlap likely reflects this common factor. Co-expression network analyses identified a circRNA module showing increased expression in ASD, the first identification of circRNA expression changes in this disorder. Elucidating the link between circRNA expression changes and mRNA expression changes in ASD brain will require further mechanistic studies on the gene regulatory function of circRNAs.

8

370 371 372 373 374 375 376 377 378 379 380 381 382 383

ACKNOWLEDGEMENTS The authors would like to thank Dr. Cristopher Pardy and Dr. Jim Fang for technical support in the initial stages of the project. This work was supported by an ARC Future fellowship to IV (FT170100359), a UNSW Scientia Fellowship to IV and a UNSW research training program scholarship to AG. This study has been posted on bioRxiv ahead of publication: doi: https://doi.org/10.1101/500991 CODE AVAILABILITY The data analysis code is provided in Supplemental Information Code. Data pre-processing code (STAR, DCC, rMATS) is briefly described in Supplemental Methods, and is available upon request. AUTHOR CONTRIBUTIONS IV conceived the study and supervised all aspects of the project. AG, FA, and IV analysed data, AG carried out experimental validations, AG and IV wrote the manuscript. COMPETING INTERESTS The authors declare no biomedical financial interests or potential conflicts of interest.

9

384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439

REFERENCES 1. Chen LL (2016): The biogenesis and emerging roles of circular RNAs. Nat Rev Mol Cell Biol. 17:205-211. 2. Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO (2013): Cell-type specific features of circular RNA expression. PLoS Genet. 9:e1003777. 3. Liang D, Tatomer DC, Luo Z, Wu H, Yang L, Chen LL, et al. (2017): The Output of Protein-Coding Genes Shifts to Circular RNAs When the Pre-mRNA Processing Machinery Is Limiting. Mol Cell. 68:940954 e943. 4. Pamudurti NR, Bartok O, Jens M, Ashwal-Fluss R, Stottmeister C, Ruhe L, et al. (2017): Translation of CircRNAs. Mol Cell. 66:9-21 e27. 5. Li L, Zheng YC, Kayani MUR, Xu W, Wang GQ, Sun P, et al. (2017): Comprehensive analysis of circRNA expression profiles in humans by RAISE. Int J Oncol. 51:1625-1638. 6. Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, et al. (2015): Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Molecular cell. 7. You X, Vlatkovic I, Babic A, Will T, Epstein I, Tushev G, et al. (2015): Neural circular RNAs are derived from synaptic genes and regulated by development and plasticity. Nature neuroscience. 18:603-610. 8. Zhang Y, Xue W, Li X, Zhang J, Chen S, Zhang J-L, et al. (2016): The biogenesis of nascent circular RNAs. Cell reports. 15:611-624. 9. Szabo L, Morey R, Palpant NJ, Wang PL, Afari N, Jiang C, et al. (2015): Statistically based splicing detection reveals neural enrichment and tissue-specific induction of circular RNA during human fetal development. Genome Biol. 16:126. 10. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, et al. (2013): Natural RNA circles function as efficient microRNA sponges. Nature. 495:384-388. 11. Piwecka M, Glazar P, Hernandez-Miranda LR, Memczak S, Wolf SA, Rybak-Wolf A, et al. (2017): Loss of a mammalian circular RNA locus causes miRNA deregulation and affects brain function. Science. 357. 12. Li Z, Huang C, Bao C, Chen L, Lin M, Wang X, et al. (2015): Exon-intron circular RNAs regulate transcription in the nucleus. Nature structural & molecular biology. 22:256-264. 13. Yu CY, Li TC, Wu YY, Yeh CH, Chiang W, Chuang CY, et al. (2017): The circular RNA circBIRC6 participates in the molecular circuitry controlling human pluripotency. Nat Commun. 8:1149. 14. Huang R, Zhang Y, Han B, Bai Y, Zhou R, Gan G, et al. (2017): Circular RNA HIPK2 regulates astrocyte activation via cooperation of autophagy and ER stress by targeting MIR124-2HG. Autophagy. 13:1722-1741. 15. Chen BJ, Mills JD, Takenaka K, Bliim N, Halliday GM, Janitz M (2016): Characterization of circular RNAs landscape in multiple system atrophy brain. J Neurochem. 139:485-496. 16. Song X, Zhang N, Han P, Moon BS, Lai RK, Wang K, et al. (2016): Circular RNA profile in gliomas revealed by identification tool UROBORUS. Nucleic Acids Res. 44:e87. 17. Maass PG, Glazar P, Memczak S, Dittmar G, Hollfinger I, Schreyer L, et al. (2017): A map of human circular RNAs in clinically relevant tissues. J Mol Med (Berl). 95:1179-1189. 18. Xia S, Feng J, Lei L, Hu J, Xia L, Wang J, et al. (2017): Comprehensive characterization of tissuespecific circular RNAs in the human and mouse genomes. Brief Bioinform. 18:984-992. 19. Sekar S, Cuyugan L, Adkins J, Geiger P, Liang WS (2018): Circular RNA expression and regulatory network prediction in posterior cingulate astrocytes in elderly subjects. BMC genomics. 19:340. 20. Parikshak NN, Swarup V, Belgard TG, Irimia M, Ramaswami G, Gandal MJ, et al. (2016): Genome-wide changes in lncRNA, splicing, and regional gene expression patterns in autism. Nature. 540:423-427. 21. Feng H, Zhang X, Zhang C (2015): mRIN for direct assessment of genome-wide and genespecific mRNA integrity from large-scale RNA-sequencing data. Nat Commun. 6:7816. 22. Oldham MC, Konopka G, Iwamoto K, Langfelder P, Kato T, Horvath S, et al. (2008): Functional organization of the transcriptome in human brain. Nat Neurosci. 11:1271-1282. 23. Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L (2014): Complementary sequencemediated exon circularization. Cell. 159:134-147. 24. Cheng J, Metge F, Dieterich C (2016): Specific identification and quantification of circular RNAs from sequencing data. Bioinformatics. 32:1094-1096. 10

440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494

25. Lake BB, Chen S, Sos BC, Fan J, Kaeser GE, Yung YC, et al. (2018): Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nat Biotechnol. 36:70-80. 26. Glazar P, Papavasileiou P, Rajewsky N (2014): circBase: a database for circular RNAs. RNA. 20:1666-1670. 27. Han Y, Kaeser PS, Sudhof TC, Schneggenburger R (2011): RIM determines Ca(2)+ channel density and vesicle docking at the presynaptic active zone. Neuron. 69:304-316. 28. Veno MT, Hansen TB, Veno ST, Clausen BH, Grebing M, Finsen B, et al. (2015): Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome Biol. 16:245. 29. Cortes-Lopez M, Gruner MR, Cooper DA, Gruner HN, Voda AI, van der Linden AM, et al. (2018): Global accumulation of circRNAs during aging in Caenorhabditis elegans. BMC Genomics. 19:8. 30. Shen S, Park JW, Lu ZX, Lin L, Henry MD, Wu YN, et al. (2014): rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci U S A. 111:E5593-5601. 31. Schizophrenia Working Group of the Psychiatric Genomics C (2014): Biological insights from 108 schizophrenia-associated genetic loci. Nature. 511:421-427. 32. Kim HG, Kishikawa S, Higgins AW, Seong IS, Donovan DJ, Shen Y, et al. (2008): Disruption of neurexin 1 associated with autism spectrum disorder. Am J Hum Genet. 82:199-207. 33. Glessner J, Wang K, Cai G, Korvatska O, Kim C, Wood S, et al. (2009): Autism genome-wide copy number variation reveals ubiquitin and neuronal genes. Nature. 459:569-573. 34. Voineagu I, Wang X, Johnston P, Lowe JK, Tian Y, Horvath S, et al. (2011): Transcriptomic analysis of autistic brain reveals convergent molecular pathology. Nature. 474:380-384. 35. Mercer TR, Dinger ME, Sunkin SM, Mehler MF, Mattick JS (2008): Specific expression of long noncoding RNAs in the mouse brain. Proceedings of the National Academy of Sciences. 105:716. 36. Sloan SA, Darmanis S, Huber N, Khan TA, Birey F, Caneda C, et al. (2017): Human Astrocyte Maturation Captured in 3D Cerebral Cortical Spheroids Derived from Pluripotent Stem Cells. Neuron. 95:779-790 e776. 37. Wu YE, Parikshak NN, Belgard TG, Geschwind DH (2016): Genome-wide, integrative analysis implicates microRNA dysregulation in autism spectrum disorder. Nat Neurosci. 19:1463-1476. 38. Kleaveland B, Shi CY, Stefano J, Bartel DP (2018): A Network of Noncoding Regulatory RNAs Acts in the Mammalian Brain. Cell. 174:350-362 e317. 39. Li JH, Liu S, Zhou H, Qu LH, Yang JH (2014): starBase v2.0: decoding miRNA-ceRNA, miRNAncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42:D92-97. 40. Langfelder P, Horvath S (2012): Fast R Functions for Robust Correlations and Hierarchical Clustering. J Stat Softw. 46. 41. Langfelder P, Luo R, Oldham MC, Horvath S (2011): Is my network module preserved and reproducible? PLoS Comput Biol. 7:e1001057. 42. Chen G, Shi Y, Liu M, Sun J (2018): circHIPK3 regulates cell proliferation and migration by sponging miR-124 and regulating AQP3 expression in hepatocellular carcinoma. Cell Death Dis. 9:175. 43. Yao Z, Luo J, Hu K, Lin J, Huang H, Wang Q, et al. (2017): ZKSCAN1 gene and its related circular RNA (circZKSCAN1) both inhibit hepatocellular carcinoma cell growth, migration, and invasion but through different signaling pathways. Mol Oncol. 11:422-437. 44. Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu J, et al. (2013): Circular RNAs are abundant, conserved, and associated with ALU repeats. Rna. 19:141-157. 45. Fromer M, Pocklington AJ, Kavanagh DH, Williams HJ, Dwyer S, Gormley P, et al. (2014): De novo mutations in schizophrenia implicate synaptic networks. Nature. 506:179-184. 46. Zoghbi HY, Bear MF (2012): Synaptic dysfunction in neurodevelopmental disorders associated with autism and intellectual disabilities. Cold Spring Harb Perspect Biol. 4.

FIGURE LEGENDS Figure 1. Dataset Overview. (A) Sample composition of dataset 1 (DS1) and dataset 2 (DS2). Left: schematic representation of brain regions included in the study. FC - frontal cortex; TC temporal cortex; CB - cerebellum. Right: barplot displaying the number of samples included in DS1 and DS2 categorised by brain region and phenotype (ASD and CTL). (B) Outline of data 11

495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546

analysis approach. Top: schematic representation of a gene expressing circRNA. Grey boxes exons; grey lines - introns; blue lines - linear junction reads; red lines - circular (i.e. backsplice) junction reads. ULJ - upstream linear junction; DLJ - downstream linear junction. CJ - circular junction. Bottom: data analysis pipeline. For details on circular-to-linear ratio calculation, see Supplemental Methods. (C) Overlap between circRNAs detected in DS1 and DS2 as a function of filtering parameters. y-axis: Overlap between DS1 and DS2 as a proportion of the smaller dataset (DS2); x-axis: number of independent samples in which detection (by either 2 back-splice junction reads or 0.1 CPM) is required. Based on these data, the chosen filtering parameter was 2 backsplice junction reads detected in a minimum of 5 independent samples (*). Figure 2. Dataset characterisation. (A) Venn diagram displaying the overlap between DS1, DS2 and circBase. The numbers between brackets show the total number of circRNAs in each dataset, after filtering for circRNAs expressed in a minimum of 5 samples in each dataset. (B) circRNA PCR validation Top: schematic display of divergent primers (black arrows) around a circRNA junction (dashed arc) which would amplify circRNAs, but not linear RNA molecules; negative control divergent primers anneal to exons that do not circularize, within the same gene. Bottom: Agarose gel electrophoresis of RT-PCR products using divergent primers around 22 selected circRNA junctions. Each RT-PCR is carried out on 3 brain samples (S1-S3). All circRNAs generate a PCR product at the expected size, except circRMST, and circDCUN1D4, for which the product is not of the expected size. Additional bands at higher length are likely rolling circle amplification products. For a subset of circRNAs, RT-PCR was carried out with primers around the circRNA junction, and negative control primers, on a distinct sample (S4; bottom panel). (C) Gene ontology enrichment of circRNA producing genes. FDR - false discovery rate. The most biologically relevant enriched terms are displayed. A full list of enriched terms is provided in Supplemental Table S4. (D) Histogram displaying novel circRNAs expressed in more than 10 samples. Left: DS1. Right: DS2. (E) Top: schematic display of seven highly expressed RIMS2 circRNAs. Top track: hg19 chromosomal coordinates and a representative Ensembl transcript annotation of the region (due to space limit only one of multiple RIMS2 transcripts annotated at this position is displayed). Each circRNA is displayed as a line spanning the interval between its start and end junctions. Red: Novel circRNAs. Black: circRNAs present in circBase. Bottom, left: Boxplots displaying circRNA expression differences between CB and CTX, in both DS1 and DS2. CI: circularization index. Only five of the seven circRNAs are displayed, which showed significant differences in CI levels between CB and CTX after correction for covariates and multiple testing (Supplemental Table S5). The horizontal line represents the median, boxes extend between the first and third quartiles, and whiskers extend from the box to 1.5x the inter-quartile range. Bottom, right: Barplot showing the number of samples in which each circRNA was detected. circRNA labels correspond to the labels from the top annotations track. Figure 3. Major circRNA isoform expression. (A) Histogram of the number of circRNA isoforms per gene (DS1). Left: all genes. Right: genes with more than 5 isoforms per gene (blowout of the data in the left panel). (B) Top: schematic representation of a hypothetical gene expressing three circRNAs (circRNA1-3), of which circRNA1 has the highest expression level. Circular junction (CJ) reads are shown in red; linear junction reads are shown in blue. Circular-to linear ratios (CLR) are calculated for each circRNA as shown in Figure 1. Bottom: Formulas for the major isoform relative expression and major isoform relative CLR for the hypothetical example shown in the top schematic. (C) Boxplots comparing the observed and random sampling values of mean major isoform relative CLR. Mean values were calculated across samples. Top: DS1, Bottom: DS2. Figure 4. Interplay between circRNA expression and alternative splicing. (A) The percent of alternatively spliced exons (Percent AS) is significantly higher among circ-forming exons (red line) than among non-circ forming exons (histogram displaying percent AS for 10,000 random samplings of non-circ forming exons with similar flanking intron length as circ-forming exons, sampled from genes with similar expression levels as circ-forming genes). P-values were calculated as the number of random sampling scores with more extreme values than the circ-exons score. Left: 12

547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576

CTX, DS1; right: CB, DS1. (B) Density plots of exon inclusion level for circ-exons (red) and noncirc exons (blue) in CTX (left) and CB (right). (C) Scatterplot of the number of alternative splicing (AS) events per gene (y-axis) and the number of circRNAs expressed per gene (x-axis) for genes expressing more than 10 circRNA isoforms. Each data point represents a gene. Figure 5. CircRNA expression variability. (A) Scatterplot showing the proportion of samples expressing a given circRNA at >= 0.1CPM in DS1 (x-axis) and DS2 (y-axis). Each data point represents a circRNA. Only circRNAs detected in both datasets are included. rho: Spearman correlation coefficient. (B) Scatterplot showing the mean CLR in DS1 (x-axis) and DS2 (y-axis). Each data point represents a circRNA. Only circRNAs detected in both datasets are included. rho: Spearman correlation coefficient. Figure 6. circRNA expression in CTX and CB. (A) Principal component plots of normalised circRNA expression (CLR) data. Left: DS1. Right: DS2. PC1, PC2 – first and second principal components. (B) Boxplot of pairwise Spearman correlation coefficients for brain samples of different individuals, within the same brain region (CTX, CB) or brain samples of the same individual but different brain regions (CTX-CB). *: p < 2.2 e-16. Wicoxon rank sum test. (C) Scatterplot of circRNA expression difference between CB and CTX in DS1 (x-axis) and DS2 (yaxis). Figure 7. CDR1as circRNA expression is correlated with CYRANO and anti-correlated with mir-671-5p and mir-7-5p. (A) Scatterplot of expression levels for CDR1as (x-axis) vs. CYRANO (y-axis) colored by mir-7-5p expression levels. (B) Scatterplot of expression levels for CDR1as (xaxis) vs. mir-671-5p (y-axis) colored by mir-7-5p expression levels. Figure 8. Co-expression networks identify circRNA expression changes in ASD. (A) Module eigengene values for M4 show significant differences between ASD and control CTX samples. p=0.006, Wilcoxon rank-sum test, Bonferroni corrected. Boxplots were generated using the boxplot function in R; the horizontal line represents the median, boxes extend between the first and third quartiles, and whiskers extend to 1.5 IQR (inter-quartile range) from the box. Notches mark +/-1.58 IQR/sqrt(n), where n represents the number of data points. (B) Network plot of M4, showing the top 20 circRNAs by kME (circles), and the top 50 connections between them as edges. Nodes are colored by M4 kME values. CircRNAs are named following the convention: Chromosome_StartCoordinate_EndCoordinate__Strand_GeneSymbol.

13

A. FC

Number of samples

60

TC

C.

Proportion of overlap between DS1 and DS2

CB

1.00 0.95 0.90 0.85 0.80 0.75 0.70

DS1 144 samples

30

40

20

20

10

0

0

TC FC CB

B.

DS2 53 samples

CTL ASD TC FC CB

2 reads 3 4 5 6 7 8 Number of samples

Gene-level expression

STAR Linear junction counts

Circular to linear ratio (CLR) CLR=CJ/max(ULJ, DLJ)

0.1 CPM 2

DLJ reads CJ reads

*

1

ULJ reads

9 10

Chimeric junction counts

DCC

Circularization index (CI) CI=CJ/(CJ + mean(ULJ, DLJ))

Circular RNAs

A.

DS1 (14,386)

627

1,081 4,587

8,091

B.

DS2 (9,440)

118 604

100 bp circCRKL ladder S1 S2 S3

127,508 circBase (140,790)

C.

Amplicon size (bp):

synaptic vesicle transport gene silencing synaptic vesicle cycle cellular metabolic process pos.reg. of dendrite development cell cycle G2/M phase synaptic vesicle localization synapse

Amplicon size (bp):

circFKBP8

S1

133

100 bp circRIMS1 ladder S1 S2 S3

300

0 20 40 60 80

200 100 0

20 40 60 80 100 120 140 Chromosome 8 (hg19) Number of samples

Amplicon size (bp):

E. RIMS2

100 bp ladder

10

20

30

40

S1

S2

148

S2

circHERC1

circHIPK3

103

297

S1

115

circRIMS2 S1

S3

circRMST

S3

S1

S2

143

S2

S3

S2

S3 S1

circSETD3 S1

155

circSMARCA5 S1

S3

S2

S2

S1

circHDAC9 S2

S3

S3

circMINDY3 circN4BP2L2 circPHF21A

S1

S2

S1

S1

S2

S3

S3

S3

S1

S2

S2

S3

S3

65

circFAM13B

circZNF292

S1

S1

S2

S3

S2

146

circDCUN1D4

S1

S2

88

91

circMAN1A2 S2

S3

circBRWD3

S3

159

S1

S2

110

circSLAIN1

S3 S1

86

280

264

186

247

70

ol ol ol ntr ntr ntr NA NA NA cR ve co ircR ve co ircR ve co cir c c HDAC9

HERC1

PHF21A

S1

S3

232

circPAIP2 S2

S3

121

S1: ANO2987_TC S2: ANO3632_FC S3: ANO4582_TC

S4: UMB5309_TC

Chromosome 8 (hg19)

50

Number of samples

104900000

circRNAs

circFUT8

S3

Divergent negative control primers

2

N novel circRNAs (DS2)

N novel circRNAs (DS1)

0.5 1 1.5 -log10(FDR)

S2

144

circSOGA2 100 bp ladder S1 S2 S3

0

105000000

105100000

ENST00000436393.2

circ3: chr8_105043908_105161076_+ circ4: chr8_105053554_105161076_+ circ5: chr8_105053557_105161076_+ circ6: chr8_105080740_105161076_+ circ7: chr8_105105699_105161076_+

&LUF

160

'6'6

'6'6 '6'6 '6'6

100 80

DS1

60

DS2

40 20 0

circ7

'6'6

120

circ6

Number of samples

 

CI

   $#
$59
$#
$59 $#
$59
$#
$59 $#
$59
$#
$59 $#
$59
$#
$59 $#
$59
$#
$59

140

circ5

&LUF

circ4

&LUF

circ3

&LUF



&LUF

circ2

circ1: chr8_104897526_104898451_+ circ2: chr8_104922549_104973361_+

circ1

D.

Divergent primers around a circRNA junction

100 0

1500 1000

5

10

20

30

40

50

N circRNA isoforms

60

0

500

N genes

2000

N genes

200

A.

0

10

20

30

40

N circRNA isoforms

B.

50

60

Major isoform CJ1 relative expression = (CJ1+CJ2+CJ3) CJ1

CJ2

CJ3

circRNA1

circRNA2

circRNA3

Major isoform CLR1 relative CLR = (CLR1+CLR2+CLR3)

C. 0HDQmDMRUiVRIRUPUHODWLYH&/5



●

●

 ●

● ● ● ●



●

5DQGRP6DPSOLQJ ● ●

●



● ● ● ●

● ●

●

●

● ●















0HDQmDMRUiVRIRUPUHODWLYH&/5

2EVHUYHG'6 









●

●

●

●

               

1XPEHURIFLUF51$LVRIRUPV

●

●

● ● ● ● ●

● ●

 ● ● ● ● ● ●

● ● ●

2EVHUYHG'6 



● ● ●

5DQGRP6DPSOLQJ ● ● ● ●



●

●

●

● ● ●























●









1XPEHURIFLUF51$LVRIRUPV

















1000

Frequency

1500

p < 1 e-5

500

CB p=1 e-4

0

0 0

5

10

15

20

30

25

Percent AS

0

5

10

15

20

25

30

Percent AS

CTX

CB

10 5

5

10

Density

15

15

B.

Density

1000 1500 2000 2500

CTX

500

Frequency

2000

A.

0.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

0.6

0.8

1.0

Exon inclusion level

Exon inclusion level

C.

$751/

●



●

●

5,06

37.

5,06

3,&$/0



●

●

'/*

● ●

●



● ● ● ● ● ● ●

● ●

●

● ● ● ● ● ●

●



●

●

● ● ●

●

● ●

● ● &+' ● &$'36 ●$5+*$3 ● 15;1 ● 5<5 ● $*73%3 ● 3&&$ ● $5+*$3

%,5&

1$/&1 ● ● ● ● ● ● ●

77&

●

●

●



1XPEHURI$6HYHQWVSHUJHQH





●

● ;32

● ● ● ● ●

●

%$,

● 7,$0

67.

●



 1XPEHURIFLUF51$VSHUJHQH UKR 





A.

5 0 -5 -15

-10

log2(Mean CLR) DS2

0.8 0.6 0.4 0.2 0.0

Proportion of DS2 samples

1.0

B.

0.0

0.2 0.4 0.6 0.8 Proportion of DS1 samples rho=0.79

1.0

-15

-10 -5 0 log2(Mean CLR) DS1 rho=0.93

5

A.

●

● ●

● ● ● ●● ● ●● ● ●● ● ● ●●● ● ●● ●● ● ●● ●● ● ● ●● ● ● ● ●● ● ● ●● ● ● ● ●

●

0.05

● ●

PC2: 7.51% -0.15 -0.05

0.10 0.00

●

-0.25

-0.10

PC2: 14.68%

CLR (DS2)

CLR (DS1) ● ● ● ●●● ● ● ●●● ● ● ● ●●● ● ● ● ● ● ●

-0.10

-0.05

0.00

0.05

0.10

● ●

● ●

●

●● ●

●

ASD Controls

● ● ● ● ● ● ● ●

●

●● ● ● ●

●

-0.15

-0.05 0.05 PC1: 20.38%

PC1: 23.97%

Mean CI difference in DS2

*

0.0

0.15

0.4

ASPHD1 RIMS2 VAMP3

RIMS2 FKBP3

ATP5C1 UBE3A KIAA0368 GRHPR TERF2 RAB23 ZBTB44 ADAM22 CLUAP1 MARK4 RIMS2

RIMS2 NEURL1

DPYSL2 TUSC3

TNRC6B

-0.2

Spearman rho 0.2 0.4 0.6

*

0.2

C.

0.0

B.

Green: Frontal Blue: Temporal Orange: Cerebellum

KCNH1

DGKI

RIMS2

KAT6B

CTX

CB

CTX-CB

-0.4

-0.2

0.0

0.2

0.4

Mean CI difference in DS1 rho = 0.91

0.6

A.

B. rho=0.63

●

!

2.5

● ● ● ●● ● ●

CYRANO

0

●

● ●

● ●

● ●

● ●● ● ●

● ●

●

● ●

−1

● ● ●

● ●● ●

●●

●

●

● ● ●

● ●

●

●● ● ●

● ●● ● ●

●

● ●● ● ● ●

●●

●

● ●● ● ● ●●● ● ● ● ● ●

●

●●

● ●

!!

●

●

●

! ●

mir-75p

● ●

2 ●

0 ● ● ●

−2

●

●

0.0

!

!

! ! ! ! ! !

!!

! ! !! ! ! ! ! ! !! ! !! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! !

mir-6715p

1

!

rho=-0.57

!

! !! ! ! !! !

! !

!

!2.5

!

mir-75p

!

! ! !!! !

!!

! !

!

! !

! !

!

!

!

2

! ! !

0

!

!

!2

! !

!

●

!

● ●

-5.0

!

! !

●

−4

−2

0

CDR1as

2

4

!4

!2

0

CDR1as

2

!

4

A.

M4 kME

B. FKUBBBB60$5&$

0.4

FKUBBBB/55&

0.62 FKUBBBB=.6&$1

FKUBBBB.,$$

0.1

0.2

FKUBBBB)$0$

FKUBBBB&61.*

-0.3 -0.2 -0.1 0.0

Module eigengene value

0.3

FKUBBBB'(.

FKUBBBB(5&

FKUBBBB$7;1

FKUBBBB'=,3

FKUBBBB70(0

ASD

FKUBBBB3,731%

CTL FKUBBBB5+'0

FKUBBBB8631

FKUBBBB&&'

FKUBBBB.&11

FKUBBBB903

0.79

The Landscape of Circular RNA Expression in the Human Brain Supplement 1

Table of Contents SUPPLEMENTAL METHODS ............................................................................................................... 1 RNA-seq DATA ANALYSIS ............................................................................................................... 1 Benchmarking analysis...................................................................................................................... 1 Human brain circRNA dataset .......................................................................................................... 2 Organoid dataset ................................................................................................................................ 9 SUPPLEMENTAL FIGURES................................................................................................................ 10 SUPPLEMENTAL TABLE LEGENDS ................................................................................................ 23 SUPPLEMENTAL INFORMATION CODE ........................................................................................ 24 SUPPLEMENTAL REFERENCES ....................................................................................................... 25

Gokool et al.

Supplement

SUPPLEMENTAL METHODS RNA-seq DATA ANALYSIS Benchmarking analysis RNA-seq data for benchmarking. CircRNA quantification was benchmarked on 5 DS1 tissue samples: 5115_ba9, 5278_ba9, 5297_ba41-42, 5308_ba41-42, 5309_ba41-42, for which we generated polyA+ and ribodepleted stranded RNA-seq data. Brain tissue from the same individuals and brain region had been obtained by our lab from the NICHD brain and tissue bank. Total RNA was extracted from ~ 100 mg of brain tissue using a Qiagen miRNeasy kit according to the manufacturer’s protocol, and treated with 1 μl DNase I (Thermo Fisher Scientific, #AM2238) per 10 μg of RNA. Each RNA sample was divided in two, for library preparation with either the TruSeq Stranded Total RNA Ribozero kit or the TruSeq Stranded mRNA kit (which includes polyA selection). Library preparation was carried out at the UNSW Ramaciotti Centre for Genomics, followed by sequencing on an Illumina NextSeq 500 sequencer to obtain 100 bp paired-end reads (Supplemental Table S1). PolyA+ and ribo-depleted stranded RNA-seq data, were mapped to the human genome (hg19) using STAR (1) with the same parameters as described below for the human brain dataset. The chimeric read alignments were used as input for either DCC (2) or CIRCexplorer (3), run with default parameters for stranded paired-end data. CircRNAs were filtered to include those detected by at least 2 back-spliced junction reads in a minimum of 2 distinct samples. CircRNA counts were normalised to library size to obtain counts-per-million (CPM). The false-positive detection rate (i.e. the percentage of circRNAs detected in DS1 data that were also detected in the polyA+ data from the same sample, normalized for sequencing depth) was < 1% for DCC and < 3% for CIRCexplorer (Supplemental Figure S10a). Notably, we observed low falsepositive rates despite the fact that the sequencing depth and the paired-end read length were higher for the polyA+ libraries than the original RNA-seq data (Supplemental Figure S10b). Since previously 1

Gokool et al.

Supplement

reported polyA+ based false positive rates were 2.7%-8% depending on algorithm (4), we concluded that both DCC and CIRCexplorer performed well on the brain dataset, with DCC showing a particularly low false-positive rate. We also assessed the effect of strand-specificity of the RNA-seq data on false-positive rate detection. We found that all false-positive circRNAs detected in the unstranded data were detected in the stranded data as well (Supplemental Figure S10c), demonstrating that lack of strand-specificity did not lead to false-positive circRNA detection. Across the five DS1 brain tissue samples, DCC and CIRCexplorer identified 5,706 and 5,426 circRNAs, respectively, with 91% of these being identified by both algorithms. Furthermore, the correlation between circRNA expression quantified by DCC and CIRCexplorer in the same sample was between 0.97 and 0.99 (Supplemental Figure S10d). This result indicated that circRNA quantification was robust to the choice of method, and due to its lower false-positive detection rate, DCC was chosen for downstream analyses.

Human brain circRNA dataset RNA sequencing reads from brain samples (DS1) were aligned to the human genome (hg19) using STAR (1) , with the following parameters, recommended for optimal circRNA detection (2) : --outSJfilterOverhangMin 15 15 15 15 --alignSJoverhangMin 15 --alignSJDBoverhangMin 15 -outFilterMultimapNmax

1

--outFilterScoreMin

outFilterMismatchNmax 2 --chimSegmentMin 15

1

--outFilterMatchNmin

1

--

--chimScoreMin 15 --chimScoreSeparation 10 -

-chimJunctionOverhangMin 15. Gene-level expression was assessed with featureCounts, as implemented in STAR. Samples with more than 3 standard deviations away from the mean of mean inter-sample correlations, as well as outliers on PCA analysis were eliminated from further analyses. For each dataset (DS1 and DS2), 2

Gokool et al.

Supplement

genes were filtered for expression at a minimum of 1 RPKM in at least 30% of the smallest sample group within each dataset (i.e. 6 samples in DS1 and 4 samples in DS2). CircRNAs were identified using STAR’s chimeric read alignments as input for DCC (2) (Supplemental Information Code), with default parameters for paired-end un-stranded data. Counts for circRNAs with identical genomic coordinates on opposite strands were summed, given the unstranded nature of the data. CircRNA counts were normalised to the total number of uniquely aligned reads, to obtain counts per million (CPM). For each circRNA, its corresponding linear junction counts were quantified using splice junction counts generated by STAR. The downstream and upstream linear junction counts were calculated as the sum of all linear junction reads spanning the start- and the end- circRNA coordinate respectively (schematic representation below). Linear junction reads were also normalised to library size (i.e. the total number of uniquely aligned reads). Upstream linear junction counts

Downstream linear junction counts

circRNA

 CircRNA annotation relative to genomic features Since the RNA-seq data was un-stranded, circRNAs were assigned a strand based on their overlap with exon junctions as follows: if either the start or the end overlapped an annotated exon junction, the circRNA was assigned the strand of the corresponding transcript. If neither start nor end overlapped an exon junction, or if they overlapped exon junctions on both strands, the strand was set as ambiguous. CircRNAs left with ambiguous strand annotation were next overlapped with gene intervals and assigned a strand in the same manner as above. CircRNA genomic annotation was done separately 3

Gokool et al.

Supplement

for the start and end coordinates with the following hierarchy: exon junction>exonic>intronic>genic (Supplemental Figure S3). CircRNAs were annotated relative to circBase downloaded 09.2017. All genomic annotations were performed relative to Ensembl GRCh37 transcript annotations downloaded 10.2014.  CircRNA filtering To select robustly expressed circRNAs we used the overlap between DS1 and DS2 as a guide for filtering criteria. We assessed the overlap between circRNAs detected in DS1 and DS2 when requiring circRNAs to be expressed at a minimum of either 2 read counts (a permissive criterion) or 0.1 CPM (stringent criterion) in a range of numbers of samples from 1 to 10 (Supplemental Figure S10e). We found that the overlap between DS1 and DS2 increased with the number of samples we required expression in, as expected, and plateaued at 5 samples. Notably, when requiring expression in a minimum of 5 samples, the overlap between DS1 and DS2 was above 90% using either the permissive or the stringent criteria, and thus we used this filtering parameter for inclusion of circRNAs in further analyses. However, we also include in the Supplemental information (Supplemental Table S3) the number of samples in which each circRNAs is expressed at a higher expression threshold (>= 0.1 CPM), to allow this resource to be easily used to identify highly expressed circRNAs. Enrichment analyses of circRNA-producing genes All enrichment analyses were carried out using the intersection of circRNA producing genes from DS1 and DS2, and the intersection of all genes expressed in DS1 and DS2 as background. Gene Ontology enrichment analysis was carried out using GOseq (5), with correction for gene length and Benjamini and Hochberg (BH) correction for multiple testing. KEGG pathway enrichment analysis was carried out using gProfiler (6) with KEGG pathways as source and Bonferroni multiple testing correction.

4

Gokool et al.

Supplement

The enrichment of circRNA genes among neuronal subtypes was carried out using marker genes from Lake et al. (7). We used as markers genes reported by Lake et al. to be significantly differentially expressed for each excitatory and inhibitory neuron cluster (Supplemental Table 3 of the Lake et al study) with the additional requirement of fold change > 2. Enrichment of circRNA genes among markers of each neuronal sub-type was assessed using a hypergeometric test as implemented in the phyper function in R, followed by Bonferroni correction. ASD-associated genes were obtained from the SFARI database (8), and we included all syndromic genes and genes with a score < 4. ID genes were obtained from Vissers et al (9). SCZ genes were obtained from SZDB (10), and we included genes supported by at least 2 lines of evidence. Enrichment of circRNA genes among ASD, ID and SCZ genes was assessed using a hypergeometric test as implemented in the phyper function in R, followed by Bonferroni correction. The threshold for significance for all enrichment analyses was set at p < 0.05 after multiple testing correction. Conservation of circRNA expression To assess circRNA expression conservation between human and mouse brain, we first converted hg19 human genomic coordinates to mm9 mouse coordinates using the liftOver tool from the UCSC genome browser. The converted coordinates were then interrogated using the mouse circRNA expression data from Rybak Wolf et al. (11) (Supplemental Table S1 of that study). The % conservation was calculated as the percentage of human circRNAs in DS1 and DS2 respectively, for which their converted genomic coordinates were found as circRNAs expressed in the mouse brain. CircRNA flanking-intron sequence complementarity analysis Introns flanking circRNA back-splice junctions were used as input for autoBLAST (12), which employs BLAST with the following settings: parameters: blastn, word size 7, output format 5, to determine sequence complementarity in circRNA flanking introns. Intron-pairing score for a given 5

Gokool et al.

Supplement

circRNA was defined as the number of reverse-complementary matches with a minimum bit score of 20. Alternative splicing analysis rMATS (13) was run with default parameters for 50 bp paired-end reads on DS1 samples contrasting CTX (Sample 1) and CB (Sample 2), and identified 145,995 cassette exons. Alternatively spliced exons were then filtered within each brain region, to include those supported by at least 2 inclusion junction reads and 2 skipping junction reads in two independent samples, leading to 33,207 AS exons in CTX and 24,380 AS exons in CB. These data are consistent with recent alternative splicing analyses in the human brain: Takata et al. 2017 identified 29,271 alternatively spliced exons using RNA-seq data from the Common Mind consortium (206 human dorsolateral-prefrontal cortex samples) (14). All comparisons between circRNA expression and AS were carried out in control samples. In-silico deconvolution One of the important properties of human brain regions is their cellular composition and layer structure, which in turn can affect the cellular composition of dissected brain samples. Therefore, we estimated in-silico the proportion of individual cell types in the brain tissue samples, using DeconRNAseq (15), and gene expression data from pure populations of immunopanned neurons, astrocytes, oligodendrocytes, microglia and endothelial cells (16) as reference transcriptomes. We found a significant difference in the proportion of neurons between CB and CTX (Supplemental Figure S11a), with lower and more homogeneous neuronal proportions in the CB samples. We also found that the first principle components of both gene expression and circRNA expression data strongly correlated with the estimated proportion of neurons (|rho| > 0.8 for gene expression PC1; |rho| > 0.6 for circRNA expression PC1; Supplemental Figure S11b-c), indicating that cellular composition is an important covariate when assessing gene and circRNA expression differences between brain regions. 6

Gokool et al.

Supplement

To estimate in-silico the proportion of individual cell types in brain tissue samples, we used DeconRNAseq, which is specifically designed for RNA-seq data (15), and two distinct reference transcriptome datasets: (a) CAGE data from cultured neurons and astrocytes from the FANTOM5 consortium (17), and (b) RNA-seq data from immunopanned human astrocytes and neurons from the Barres lab (Zhang et al., 2016 (16)). We found a very high correlation between the neuronal proportion estimates based on the two reference transcriptomes (Spearman rho > 0.9, Supplemental Figure S12a). We also found high correlation between neuronal proportion estimates obtained by FANTOM5 or Zhang et al., 2016 reference transcriptomes and the expression level of neuronal-specific genes (MAP2, RBFOX1, Supplemental Figure S12b). Given the highly similar results obtained with neuronal proportion estimates based on the two reference transcriptome datasets, we used the cell-type reference transcriptome data from Zhang et al., 2016 in all downstream analyses. We also assessed cellular composition using Neuroexpresso (18) using mouse cortex-derived marker genes for the CTX data, and cerebellum-derived marker genes for the CB data. We calculated marker-gene based proportion of neurons (MGP-neurons) for CB and CTX samples as the sum of cerebellum neuron subtypes and cortical neuron subtypes respectively. We reproduced the observation of lower neuronal estimates in CB compared to CTX in DS1 and DS2 using this additional approach (Supplemental Figure S12c). Differential expression with brain region and age To assess differential circRNA expression with brain region and age, we applied a linear model to circRNA expression levels normalised to their linear transcript (CI), using the control samples (DS1). Only circRNAs expressed in at least half of the samples were included in the analysis. Data was log2-transformed with an offset of 0.5. The following variables were included in the model: proportion of neurons, sex, median 3’bias, sequencing batch, brain bank, age, and brain region (CTX and CB). The same linear model was applied to gene-level expression data (RPKM). P-values for the linear 7

Gokool et al.

Supplement

model, obtained using the summary function in R, were corrected for multiple testing using a BH correction

as

implemented

in

the

Bioconductor

multtest

package

(https://www.bioconductor.org/packages/release/bioc/html/multtest.html). We did not identify any circRNAs differentially expressed with age (adjusted p < 0.05). We also assessed circRNA expression changes with age using spline regression, with knots at 0, 10, 20, and 40 years old, and degree=1, using the

basic

spline

(bs)

function

in

the

splines

R

package

(https://cran.r-

project.org/web/packages/splines/index.html), with the same result. The result remained the same, whether or not we included the proportion of neurons in the model. To assess whether this result is replicable in the smaller dataset (DS2), we applied the same data analysis approach as above in DS2. We included both ASD and control DS2 samples, in order to increase statistical power, and added phenotype to the list of co-variates. Co-expression network analyses Network analysis was carried out in the larger dataset (DS1) using circRNAs expressed in at least half of the CTX samples (DS1). Normalised circRNA expression values were first corrected for covariates using a linear model, and the residual values were used for network construction. The coexpression network was constructed using the blockwiseModules function in the WGCNA R package (19)

with

the

following

parameters:

power=10,

networkType="signed",

corFnc="bicor",

minModuleSize=10, mergeCutHeight=0.2. The beta power was chosen so that the network fulfilled scale-free topology (r2 > 0.5). CircRNAs were assigned to a module based on their correlation to the module eigengene value (kME > 0.1) and a significant BH-corrected p-value for this correlation (adjusted p < 0.05). kME values are listed in Supplemental Table S8. Gene ontology enrichment analysis of the M4 module was carried out using gProfiler (6). Cell type enrichment analysis was carried out using single-cell based cell type markers from PsychEncode (http://resource.psychencode.org/), and a hypergeometric test as implemented in the p.hyper function 8

Gokool et al.

Supplement

in R. The overlap with ASD, ID, and SCZ genes was assessed using a hypergeometric test, with genes expressed in DS1 as a background. Module preservation (20) between DS1 and DS2 was assessed based on circRNAs expressed in at least half of the samples in each dataset (Supplemental Figure S9), using the modulePreservation function from the WGCNA R package, with the DS1 network as reference. CircRNA-miRNA interaction analysis miRNA expression data was obtained from (21), which included 95 brain samples present in DS1 (Supplemental Table S7; 46 frontal cortex, 46 temporal cortex and 52 cerebellum samples). We selected circRNAs expressed in at least half of these brain samples. Residual values for both miRNA expression data and circRNA expression were obtained using a linear model to correct for covariates (age, sex, median 3’bias, sequencing batch, and brain bank). miRNA-circRNA Spearman correlation values were calculated using the residuals. CLIP-seq data on circRNA-miRNA interactions were obtained using the STARBASE API (starbase.sysu.edu.cn) (22) with the following parameters: assembly=hg19,

geneType=circRNA,

miRNA=all,

clipExpNum=5,

degraExpNum=0,

pancancerNum=0, programNum=1, program=PITA,RNA22, target=al, cellType=all. Organoid dataset Mapping and circRNA quantification for RNA-seq data from brain organoids were carried out as described above for human brain data. CircRNAs were included in the dataset if they were detected in at least 2 samples.

9

Gokool et al.

Supplement

SUPPLEMENTAL FIGURES Divergent primers

1

2

2

1 BSJ

circFAM13B

circMAN1A2

circSMARCA5

BSJ

BSJ

BSJ

circBRWD3

circSOGA2

circRIMS2

BSJ

BSJ

BSJ

circHERC1

circSLAIN1

circRMST BSJ

BSJ

BSJ

circFKBP8 Reference sequence Sanger sequence

circHIPK3 Reference sequence Sanger sequence BSJ

BSJ

circFUT8 Reference sequence Sanger sequence

circRIMS1 Reference sequence Sanger sequence

BSJ

BSJ

Sample: UMB5278_BA9

Supplemental Figure S1. PCR amplification of back-splice junction validated by Sanger sequencing. Top: Schematic display of divergent primers (purple arrows) that anneal to circRNA exons (grey boxes) and amplify across the circRNA back-splice junction (red line). Bottom: The chromatogram of a Sanger sequencing experiment confirms the presence of the back-splice junction (red dashed line) for 13 circRNAs. 10

Gokool et al.

Supplement

Supplemental Figure S2. Characterisation of DS1 and DS2 sample composition. (A) DS1. (B) DS2. p: Wilcoxon rank-sum test p-values for age, and Fisher test p-values for gender ratios. Boxplots were generated using the boxplot function in R; the horizontal line represents the median, boxes extend between the first and third quartiles, and whiskers extend to 1.5 IQR (inter-quartile range) from the box. Notches mark +/-1.58 IQR/sqrt(n), where n represents the number of data points.

11

Gokool et al.

Supplement

DS2 Novel

exonic

DS2 All

exonJ exonJ-exonJ DS1 Novel

intergenic intronic

DS1 All 0%

20%

40%

60%

80%

100%

Supplemental Figure S3. CircRNA annotation relative to genomic features. exonJ-exonJ: circRNAs for which both ends correspond to annotated exon junctions. The rest of circRNAs were annotated based on the overlap of at least one end with genomic features, using the following hierarchy: exonJ > exonic > intronic > intergenic.

12

Gokool et al.

Supplement

Supplemental Figure S4. Major isoform agreement across datasets. A. Barplot displaying the number of circRNAs called major isoforms in DS1, classified on whether they were also called a major isoform in DS2 (TRUE) or not (FALSE). B. Histogram displaying the proportion of samples in which DS1 major isoforms were identified as a major isoform in DS2. CircRNAs identified as a major isoform in at least 50% of samples in each dataset were called a major isoform in that dataset.

13

Gokool et al.

Supplement

Supplemental Figure S5. Intron-pairing rank of circRNA major isoforms. Genes are classified based on how many circRNAs they express, and for each class (X-axis), the intron pairing rank of the major isoform is plotted in % (Y-axis). Rank=1: highest intron pairing score; Rank=10: lowest intron pairing score. Only genes expressing up to 10 circRNAs are plotted.

14

Gokool et al.

Supplement

Supplemental Figure S6. CircRNA expression variability. (A) Scatterplots of mean expression (xaxis) and standard deviation (y-axis). Left: genes, centre: splice junctions, right: circRNAs. CV: mean coefficient of variation. (B) Histograms displaying the number of samples in which genes (left), splice junctions (centre), and circRNAs (right) were expressed.

15

Supplement

Number of circRNAs expressed

Gokool et al.

6000

BroadRegionID 4000




CB CTX

2000

0

20

40

60

Age

Supplemental Figure S7. CircRNA expression variation with age. Scatterplot displaying the total number of circRNAs expressed in each sample (y-axis; Supplemental Methods) versus age (x-axis) in CB and CTX. Each data point represents a sample. Regression lines were generated using a loess smoothing function (geom_smooth) in ggplot2.

16

Gokool et al.

Supplement

Number of circRNAs expressed

750

500

CellType astrocytes neurons

250

0

100

200

300

400

500

TimePoint (days)

Supplemental Figure S8. CircRNA expression changes during cellular maturation in neurons and astrocytes. The number of circRNAs expressed in astrocytes and neurons (y-axis) at >= 0.1 CPM vs. the maturation time point in organoid cultures (x-axis).

17

Gokool et al.

Supplement

Supplemental Figure S9. WGCNA module preservation across datasets. Scatterplot of module preservation Z-summary vs. module size, using the DS1 circRNA network as a reference. Z-summary is a composite z-score combining the z-score statistic for multiple preservation measures of DS1 circRNA network modules into the DS2 circRNA network modules.

18

Gokool et al.

Supplement

Supplemental Figure S10. Benchmarking of circRNA expression data. (A) Barplots of false positive rates defined as the % of circRNAs detected in each DS1 sample that were also detected in polyA+ data from the same sample, normalized for the sequencing depth. (B) Barplots of sequencing depth (millions of uniquely mapped reads) for DS1 libraries, as well as RNA-seq data generated in the present study (polyA+ and RD_stranded). RD: ribo-depleted. (C) Barplots of the number of circRNAs 19

Gokool et al.

Supplement

detected in polyA+ libraries (i.e. false-positives) for DS1 data, ribo-depleted stranded data, as well as those common between both types of data. RD: ribo-depleted. Left: DCC, Right: CIRCexplorer. (D) Scatterplots of normalized circRNA expression levels (CPM: counts per million) quantified by DCC and CIRCexplorer. r: Pearson correlation coefficient. (E) CircRNA detection using DCC across various expression thresholds. Number of circRNAs detected by using a threshold of either 2 backsplice junction reads (left), or 0.1 CPM (right), in a minimum of 1 to 10 independent samples.

20

Gokool et al.

Supplement

Supplemental Figure S11. Estimated proportion of neurons in DS1 and DS2. (A) Boxplots displaying the estimated proportion of neurons across brain regions and phenotypes, in DS1 (left) and DS2 (right). Boxplots were generated using the geom_boxplot and geom_violin functions in R; the horizontal line represents the median, boxes extend between the first and third quartiles, and whiskers extend to 1.5 IQR (inter-quartile range) from the box. Notches mark +/-1.58 IQR/sqrt(n), where n represents the number of data points. (B) Scatterplots of first principal component values of gene expression data (PC1, y-axis) vs. estimated proportion of neurons. Triangles: Controls; Squares: ASD. (C) Scatterplots of first principal component values of circRNA expression data (PC1, y-axis) vs. estimated proportion of neurons. All neuronal proportion estimates are based on reference transcriptome data from Zhang et al. 2016. Triangles: Controls; Squares: ASD. 21

Gokool et al.

Supplement

Supplemental Figure S12. Assessment of cellular composition estimates. (A) Scatterplot of estimated neuronal proportions based on the FANTOM5 and Zhang et al. reference transcriptome data in DS1 and DS2; rho: Spearman correlation coefficient. (B) Scatterplot of gene expression levels of two neuronal-specific genes (RBFOX1-top row, and MAP2-bottom row; y-axis) vs. estimated neuronal proportions (x-axis). rho: Spearman correlation coefficient. (C) Marker-gene based proportion estimates calculated using Neuroexpresso.

22

Gokool et al.

Supplement

SUPPLEMENTAL TABLE LEGENDS Supplemental Table S1. Sample description (A: DS1, B: DS2) and summary of mapping results for all RNA-seq data included in this study (C:DS1, D:DS2, E:Brain Organoid data, F: benchmarking RNA-seq data) Supplemental Table S2. Primer sequences for RT-PCR validation of circRNAs. Supplemental Table S3. CircRNA annotation (A: DS1, B: DS2) and normalised circRNA expression data in CPM (C: DS1, D: DS2). Supplemental Table S4. Functional characterisation of circRNA-producing genes. A: Gene ontology enrichment. B: Pathway enrichment. C. Enrichment for neuronal subtype markers from lake et al. 2018. D. Enrichment for Neurodevelopmental disease (NDD) genes. See Supplemental methods for details on all enrichment analyses. Supplemental Table S5. CircRNAs differentially expressed between CTX and CB. A: DE results with correction for proportion of neurons. B. DE results without correction for proportion of neurons. Supplemental Table S6. CircRNA expression during astrocyte and neuronal maturation in organoid cultures. Raw sequencing data was obtained from SRA: GSE99951. CircRNA expression values are in CPM. Sample names include SRA id, time point (days), and either “Hepa” for the astrocyte immunepanned samples, or “Thy-1” for neuronal immune-panned samples. Supplemental Table S7. miRNA-circRNA interactions. A predicted miRNA-circRNA interactions. The table includes all miRNA-circRNA pairs supported by Ago2 CLIP-seq data from STARBASE (Supplemental Methods), which showed an absolute Spearman correlation coefficient > 0.5 between circRNA and miRNA expression in the brain. B. List of samples for the circRNA-miRNA interaction analysis, for which both circRNA expression and miRNA expression data were available. Supplemental Table S8. CircRNA WGCNA results. A. kME data. kME: correlation between individual circRNA expression and module eigengene values. pvalBH: Student p-values for the kME correlation values, corrected for multiple testing using a Benjamini and Hochberg correction. B. Gene ontology enrichment of the M4 ASD-associated module. C. Cell type marker enrichment of the M4 ASD-associated module.

23

Gokool et al.

Supplement

SUPPLEMENTAL INFORMATION CODE Data analysis code is available as a Github repository: https://github.com/Voineagulab/circRNA_HumanBrain

24

Gokool et al.

Supplement

SUPPLEMENTAL REFERENCES 1.

Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. (2013): STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 29:15-21.

2.

Cheng J, Metge F, Dieterich C (2016): Specific identification and quantification of circular RNAs from sequencing data. Bioinformatics. 32:1094-1096.

3.

Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L (2014): Complementary sequence-mediated exon circularization. Cell. 159:134-147.

4.

Szabo L, Morey R, Palpant NJ, Wang PL, Afari N, Jiang C, et al. (2015): Statistically based splicing detection reveals neural enrichment and tissue-specific induction of circular RNA during human fetal development. Genome Biol. 16:126.

5.

Young MD, Wakefield MJ, Smyth GK, Oshlack A (2010): Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14.

6.

Reimand J, Isserlin R, Voisin V, Kucera M, Tannus-Lopes C, Rostamianfar A, et al. (2019): Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nat Protoc. 14:482-517.

7.

Lake BB, Chen S, Sos BC, Fan J, Kaeser GE, Yung YC, et al. (2018): Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nat Biotechnol. 36:70-80.

8.

Abrahams BS, Arking DE, Campbell DB, Mefford HC, Morrow EM, Weiss LA, et al. (2013): SFARI Gene 2.0: a community-driven knowledgebase for the autism spectrum disorders (ASDs). Mol Autism. 4:36.

9.

Vissers LE, Gilissen C, Veltman JA (2016): Genetic studies in intellectual disability and related disorders. Nat Rev Genet. 17:9-18.

10. Wu Y, Yao YG, Luo XJ (2017): SZDB: A Database for Schizophrenia Genetic Research. Schizophr Bull. 43:459-471. 11. Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, et al. (2015): Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Molecular cell. 12. Cortes-Lopez M, Gruner MR, Cooper DA, Gruner HN, Voda AI, van der Linden AM, et al. (2018): Global accumulation of circRNAs during aging in Caenorhabditis elegans. BMC Genomics. 19:8. 13. Shen S, Park JW, Lu ZX, Lin L, Henry MD, Wu YN, et al. (2014): rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci U S A. 111:E55935601. 14. Takata A, Matsumoto N, Kato T (2017): Genome-wide identification of splicing QTLs in the human brain and their enrichment among schizophrenia-associated loci. Nat Commun. 8:14519. 15. Gong T, Szustakowski JD (2013): DeconRNASeq: a statistical framework for deconvolution of heterogeneous tissue samples based on mRNA-Seq data. Bioinformatics. 29:1083-1085. 16. Zhang Y, Sloan SA, Clarke LE, Caneda C, Plaza CA, Blumenthal PD, et al. (2016): Purification and Characterization of Progenitor and Mature Human Astrocytes Reveals Transcriptional and Functional Differences with Mouse. Neuron. 89:37-53. 17. Consortium F, the RP, Clst, Forrest AR, Kawaji H, Rehli M, et al. (2014): A promoter-level mammalian expression atlas. Nature. 507:462-470. 18. Mancarci BO, Toker L, Tripathy SJ, Li B, Rocco B, Sibille E, et al. (2017): Cross-Laboratory Analysis of Brain Cell Type Transcriptomes with Applications to Interpretation of Bulk Tissue Data. eNeuro. 4.

25

Gokool et al.

Supplement

19. Langfelder P, Horvath S (2008): WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 9:559. 20. Langfelder P, Luo R, Oldham MC, Horvath S (2011): Is my network module preserved and reproducible? PLoS Comput Biol. 7:e1001057. 21. Wu YE, Parikshak NN, Belgard TG, Geschwind DH (2016): Genome-wide, integrative analysis implicates microRNA dysregulation in autism spectrum disorder. Nat Neurosci. 19:1463-1476. 22. Li JH, Liu S, Zhou H, Qu LH, Yang JH (2014): starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42:D92-97.

26

KEY RESOURCES TABLE Resource Type

Specific Reagent or Resource

Add additional rows as needed for each Include species and sex when applicable. resource type

Source or Reference

Identifiers

Additional Information

Include name of manufacturer, company, repository, individual, or research Include catalog numbers, stock numbers, database IDs or accession numbers, and/or RRIDs. RRIDs are highly lab. Include PMID or DOI for references; use “this paper” if new. Include any additional information or encouraged; search for RRIDs at https://scicrunch.org/resources. notes if necessary.

Deposited Data; Public Database

Human postmortem brain RNA-seq data

Deposited Data; Public Database

Human organoid brain RNA-seq data

https://www.ncbi.nlm.nih.gov/pubmed/28817799

Software; Algorithm

Data analysis code

this paper

https://www.ncbi.nlm.nih.gov/pubmed/27919067 https://github.com/Voineagulab/circRNA_HumanBrain