Altered expression pattern of circular RNAs in primary and metastatic sites of epithelial ovarian carcinoma

Abstract

Recently, a class of endogenous species of RNA called circular RNA (circRNA) has been shown to regulate gene expression in mammals and their role in cellular function is just beginning to be understood. To investigate the role of circRNAs in ovarian cancer, we performed paired-end RNA sequencing of primary sites, peritoneal and lymph node metastases from three patients with stage IIIC ovarian cancer. We developed an in-house computational pipeline to identify and characterize the circRNA expression from paired-end RNA-Seq libraries. This pipeline revealed thousands of circular isoforms in Epithelial Ovarian Carcinoma (EOC). These circRNAs are enriched for potentially effective miRNA seed matches. A significantly larger number of circRNAs are differentially expressed between tumor sites than mRNAs. Circular and linear expression exhibits an inverse trend for many cancer related pathways and signaling pathways like NFkB, PI3k/AKT and TGF-β typically activated for mRNA in metastases are inhibited for circRNA expression. Further, circRNAs show a more robust expression pattern across patients than mRNA forms indicating their suitability as biomarkers in highly heterogeneous cancer transcriptomes. The consistency of circular RNA expression may offer new candidates for cancer treatment and prognosis.

Keywords: circular RNA; lymph node metastasis; miRNA; ovarian cancer; peritoneal metastasis.

Figure 1

A. circRNA detection pipeline. circRNA candidates are inferred from alignment of paired-end RNA-seq data to a reference scrambled exome. Read-pair alignments that contain apparent backsplice junction are reported as either junctional or supportive evidence for the backspliced exon junctions. B. and C. Sensitivity and Precision of circRNA detection pipeline. Simulated reads at different read depths representing the expression level of randomly selected backsplice junctions mixed with a background of their canonical transcripts were run through the detection pipeline. Measurements of the correctly called circRNA junctions (true positives;TP), incorrectly called junctions (false positives, FP) and unidentified junctions (false negatives, FN) were used to define sensitivity and precision as TP/(TP+FN) and TP/(TP+FP), respectively.

Figure 2

A. Divergent primers become properly inward facing and identify backspliced junction for circRNAs. Two sets of outward facing primers (a & b) with respect to genomic sequence were designed to amplify the backsplice junction sequence for circRNA candidates from ARHGAP5, SPECC1 and NFATC3 genes and each produced an expected size band in the qPCR assay (Supplementary File S1). B. Sanger sequencing confirms head-to-tail splicing for NFATC3. RNA samples were reverse transcribed with thermo-stable reverse transcriptase and a primer complementary to the coding RNA strand. PCR on the cDNA produced with primers 1 and 2 resulted in large amounts of product that could be Sanger sequenced. Primers 1 to 4 were used in sequencing, with coverage from each primer indicated in green. Expansion of the sequence provided by primer 3 shows complete coverage of the non-canonical junction. C. circRNA candidates are enriched after the RNase R digestion. Total RNA from OVCAR3 and SKOV3 cell lines was treated with RNAse R to test the sensitivity of the identified circRNA candidates to exoribonuclease digestion. Boxplots represent FPKM (Fragments per kilobase of exons per million mapped fragments) fold change for RNAse R treated vs. untreated samples and shows a significant enrichment of both junctional and supported circRNAs (n=7903) compared to mRNA (n=2785) in the RNAse R digested samples. Only genes that produced at least one detectable circular isoform were considered for this analysis. D. qPCR confirms RNAse R enrichment for circular junctions. Total RNA from OVCAR3 cell line was first treated with DNase1 to avoid any DNA contamination and then part of it was subjected to the RNase R digestion followed by RT-PCR for each sample. The four putative circular RNAs show enrichment for the backsplice junction compared to the linear RNA from beta-actin used as a negative control. E. Candidate circRNAs are devoid of Poly-A tails. The four candidate circRNAs were also tested for the presence /absence of a ploy-A tail by using either oligo (dT) or random hexamers to amplify the total RNA followed by a qPCR assay with primers specific for backsplice junctions of candidate circRNAs and linear RNA of b-actin and GAPDH genes. While both oligo (dT) and random hexamers could amplify the RNA from beta-actin and GAPDH, only random hexamers showed detectable PCR products for the four tested candidate circRNAs. The red rectangles show absence of any PCR bands for the backsplice junctions.

Figure 3

A. Isoform diversity of candidate circRNA versus linear transcripts. Violin plots - a combination of boxplot and kernel density function compares the abundance levels of candidate circRNA isoforms and cufflinks assembled linear transcripts for 7121 genes for which at least one detectable backsplice junction was observed. The diversity of candidate circRNA isoforms is significantly higher (wilcoxon test; p-value 3.324e-07) than that of linear transcripts in ovarian cancer. B. Usage of backsplice capable exons as linear and circular forms. For the backsplice junction forming exons, ratio of read pairs supporting the scrambled form to the read pairs supporting the linear form defines the relative abundance of candidate circRNA expression. Most backsplice capable exons are relatively more abundant in the linear form, however a significant proportion (5-10%) of these are preferentially expressed as circRNA form. Colors indicate ratio of circular to linear expression. C. Heatplot of relative abundance levels. Each line in the plot is a circRNA candidate as defined by two backsplice capable exons with colors indicating relative abundance levels of circular over linear expression for the two exons. Hierarchically clustered heatmap of the relative abundance levels shows the clear trends of stronger intra-patient homogeneity of circRNA expression as indicated by a distinct dendrogram branch for each patient. Sub-branching for patients 1 & 2 further reveals a primary tumor and metastases specific expression trend for these groups.

Figure 4

A, B. Micro RNAs are predicted to interact with circRNAs. Scatterplot of minimum free energy (MFE) versus binding energy (ΔG binding) for predicted structural dimers of circRNA-miRNA pairs are analogous to that of 3′UTR-miRNA seed matches, indicating the thermodynamic feasibility of these interactions at least similar to 3′UTR-miRNA binding events. A high correlation between MFE and ΔG binding for a vast majority of circRNA-miRNA cofolded dimers indicates favorable interaction. Only a marginal fraction (2%) of predicted interactions have ΔG binding >=0 (blue line separated). C, D. circRNA-miRNA dimer cofolds are predicted to be more stable. Binding energy distribution curves are similar for 3′UTR-miRNA and circRNA-miRNA structural dimers. The left shift of circRNA-miRNA curve for minimum free energy of dimer cofolds is highly significant (Wilcoxon rank sum test; p-value < 0.00001) and indicates more stability for circRNA-miRNA predicted secondary structures compared to 3′UTR-miRNA dimers. E, F. circRNA candidates are enriched for miRNA seed matches. Boxplots and Empirical cumulative distribution function (ECDF) curves showing a higher proportion of circRNA sequences having greater density of miRNA seed matches in comparison to 3′UTR, 5′ UTR regions and CDS.

Figure 5. Clustering analysis indicates tumor stage specific expression trend for circRNAs

Heat plots of A. circRNA and B. mRNA expression in tumor lesions. Each line in the plots is an isoform with colors indicating average expression (Z-score of log2-cpm values) from three patient samples for ovary, peritoneum and lymph node. Compared to mRNA, a higher fraction of circRNAs show tumor stage specific expression trend.

Figure 6. Ingenuity® Pathway Activity Analysis for circRNA and mRNA expression trends

The figure represents a hierarchically clustered heatmap of the Pathway activity z-score which predicts the direction of change for a biological function or process. A positive z-score predicts that the biological process or disease is trending towards an increase and a negative score reflects a decrease. Signaling pathways that are typically over expressed for mRNA in metastatic tumors are downregulated for circRNAs.