Comprehensive profile and contrastive analysis of circular RNA expression in cervical squamous carcinoma and adenocarcinoma

Specimens preparation

Three pairs of samples of stage IB1 cervical cancer tissue and adjacent normal tissue (ANT) were selected for RNA-seq from three HPV-16 positive patients with SCC, ADC or ASC (all negative to other HPVs). ADC was limited to endometrioid adenocarcinoma, the usual type of adenocarcinoma, and ASC was limited to the term with both squamous and glandular areas that were each clearly recognizable without the use of special stains, the usual type of adenosquamous to carcinomas, according to histopathology. Diagnosis was based on the International Federation of Gynecology and Obstetrics (FIGO) criteria. For surgery samples, tumor tissues were taken from the center of the tumor. The corresponding adjacent noncancerous cervical tissues were taken at least two cm away from the edge of the cancer. All tissues were placed in liquid nitrogen immediately after being dissected and stored the same way. Other six pairs of samples of HPV-16 mediated SCC and adjacent normal tissues, and six matched HPV-16 mediated ADC samples, were used for quantitative real time polymerase chain reaction (qRT-PCR) validation. All procedures were approved by the ethics committee of the Peking University People’s Hospital (NO.2019PHB212-01), and written informed consent was obtained from all patients.

RNA extraction and quality control

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The integrity of RNA was evaluated by standard denaturing agarose gel electrophoresis and using the Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Purity and concentration, as well as the preliminary quantification, were determined using a Nano Drop spectrophotometer (Nano Drop, Wilmington, DE, USA).

Library building for RNA-seq

After removing ribosomal RNA and building a special library, RNA-seq data collection was conducted using the Illumina PE150 platform by Novogene Bioinformatics Technology Co. Ltd. (Beijing, China). Clean reads for subsequent analyses were obtained through raw data filtering, sequencing error rate check and GC content distribution check. Clean reads were then compared and mapped to the reference genome. Find_circ and CIRI2 were used to detect and identify circRNAs (Gao et al., 2017; Memczak et al., 2013). Following identification of circRNAs, length distribution and sources of known or novel circRNAs were counted. Density statistics and circRNA locations on each chromosome were identified by Circos software for all circRNAs of each sample, and compared to all chromosomes (Krzywinski et al., 2009). Normalization of readcount by TPM was conducted before expression analysis (Zhou et al., 2010).

Raw sequencing data has been uploaded to NCBI’s Gene Expression Omnibus (GEO) and is accessible through GEO Series accession number GSE208089.

Bioinformatic analysis

Differentially expressed circRNAs between three paired samples (S_1 vs. S_2 (cervical squamous cell carcinomas tissue vs. paired adjacent normal tissue), A_1 vs. A_2 (cervical adenocarcinoma tissue vs. paired adjacent normal tissue), AS_1 vs. AS_2 (adenosquamous carcinoma of cervix vs. the paired adjacent normal tissue)) were determined by negative binomial distribution test using the DESeq R package (Anders & Huber, 2010), with criteria: |log2(fold change)| > 1 and q value < 0.05. Three other extrinsic matched combinations were explored CCT_S vs. ANT_S (cervical cancer containing squamous cell carcinoma tissues vs. corresponding adjacent normal tissues, S_1&AS_1 vs. S_2&AS_2), CCT_A vs. ANT_A (cervical cancer containing adenocarcinoma tissues vs. corresponding adjacent normal tissues, A_1&AS_1 vs. A_2&AS_2), CCT vs. ANT (cervical cancer containing three pathological types above vs. corresponding adjacent normal tissues,S_1&AS_1&A_1 vs. S_2&AS_2&A_2), with p adj < 0.05. P value < 0.05 was only adopted when the differentially expressed circRNA were too little.

Based on the threshold values, we performed a hierarchical clustering analysis and generated a series of volcano plots in order to filter circRNAs with modified expression levels. The IRESfinder software based on the logit model was used as a prediction tool to determine the coding potential of circRNAs (Zhao et al., 2018). Subsequently, we predicted the interactions between circRNAs and miRNAs using the miRanda databases, and the interactions between miRNAs and genes of the eight candidate circRNAs using TargetScan v7.2 (Agarwal et al., 2015; Enright et al., 2003).

Gene Ontology (GO) analysis via GOseq software was carried out (Young et al., 2010). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed with the Kobas3.0 software to detect the involvement of circRNAs genes in different biological pathways (Kanehisa et al., 2008; Mao et al., 2005). We used the DIANA-miRPath v3.0 and the DIANA-microT-CDS algorithm to determine the functions of target genes involved in different biological pathways (Paraskevopoulou et al., 2013; Vlachos et al., 2015). The significance threshold was set at p value < 0.05 and FDR value < 0.01. The network connecting candidate circRNAs, target miRNAs and the mRNAs of host genes was generated and displayed via Cytoscape v3.6.1 software (https://cytoscape.org/).

Quantitative real-time PCR

To validate circRNAs data generated from RNA-seq, ten significantly dysregulated circRNAs in all six comparisons (|log2(fold change)| > 2) were randomly selected for qRT-PCR. Other eight candidate circRNAs extracted from functional enrichment analysis were also selected. We designed divergent primers across the circular junctions as listed in Table S1. Total RNA was reverse transcribed using the Invitrogen Superscript cDNA Synthesis kit (Invitrogen, Carlsbad, CA, USA). CircRNA expression was measured through qPCR (SYBR Green PCR Master Mix; Applied Biosystems, Foster City, CA, USA). Reactions were performed in triplicate according to the manufacturer’s protocol. Relative circRNAs expression levels were calculated via the 2^−ΔΔCt method and β-actin was used as the housekeeping gene. Results are expressed as mean ±  SD (Standard Deviation). Statistical analysis was conducted using the SPSS Statistics 19.0 software. Significant difference between comparison was defined as a minimum 2-fold change in normalized expressed level with p value < 0.05.

General expression profiles of circRNAs in different types of cervical cancer

Identity and abundance of circRNAs of three paired samples of cervical cancer and adjacent normal tissues as described in the Methods section were investigated by RNA-seq, as shown in Fig. 1. From sequencing we identified a total of 27,148 circRNAs, 11,685 of which had annotations in the circBASE database. Based on genomic features and locations, we observed that the lengths of spliced circRNAs were mostly below 1200nt, and the corresponding genes mainly consisted of exonic and intronic sequences, while a small subset was derived from intergenic sequences. CircRNAs were distributed across all autosomes, as well as chromosome X, as shown in circos plots (Fig. 2). We found Chr 1,2,3 contained more circRNAs than the others. A full list of the identified circRNAs is available in File S1, including annotations, chromosomal locations, strand orientations, and source genes. In view of the accumulating evidence suggesting circRNAs have coding potential, data was analyzed and the predicted results of the identified circRNAs are shown in File S2.

Differential circRNA clustering analysis was used to determine the clustering pattern of differentially expressed circRNAs in different comparisons. We used the TPM value of differential circRNA sets in each sample for hierarchical clustering analysis, and for K-means and SOM clustering analysis (File S3). Discrimination among different comparison groups is shown in the cluster heatmap (Fig. 3A).

Volcano plots were created to visualize differentially expressed circRNAs between the six comparison groups (Fig. 3B). When we raised the threshold to 4 fold-change, we identified 42 up-regulated and 98 down-regulated circRNAs in the three intrinsic paired samples (A_1 vs. A_2, AS_1 vs. AS_2, and S_1 vs. S_2). If statistics was done in combination with the other three groups of comparison (CCT_A vs. ANT_A, CCT_S vs. ANT_S, CCT vs. ANT), there were only 10 and eight circRNAs left for further analysis, respectively (Table 1).

There were 215 circRNAs up-regulated in SCC (S_1 vs. S_2), but down-regulated in ADC (A_1 vs. A_2), and 50 circRNAs down-regulated in SCC (S_1 vs. S_2), but up-regulated in ADC (A_1 vs. A_2),. All circRNAs were equally expressed (|log2(fold change)| < 1) in ADC and in the corresponding ANT (File S4).

Confirmation of results with RNA-seq for differentially expressed circRNAs

Ten overlapping significantly expressed circRNAs were randomly selected to confirm the reliability of our sequencing results: five up-regulated and five down-regulated level > 4 fold circRNAs in the six comparisons described above (Table 1). We validated the expression levels in six SCC and ADC samples versus their corresponding ANT samples. The results revealed that qRT-PCR-based fold-changes mirrored those observed in the sequencing data (Fig. 4, File S5).

Comprehensive analysis of function enrichment derived from different comparisons

In order to predict potential biological functions of differentially expressed circRNAs, we conducted GO and KEGG enrichment analysis of the relative genes (“source genes” or “host genes”) in three intrinsic comparison groups. We found the most significant terms in all GO categories (BP/CC/MF) were consistent among pathologic variants of cervical cancer, except the BP term “organelle organization” in ASC (Fig. 5). KEGG pathway analysis was conducted to determine the involvement of host genes in different biological pathways. We found the most enriched pathway in all the three pathologic variants of cervical cancer was “ubiquitin mediated proteolysis” (Fig. 6). In addition, KEGG analyses were performed on the reversely expressed circRNAs in SCC and ADC as listed in File S4, details are shown in File S6.

Relative expression level detection of key candidate circRNAs

In the three intrinsic comparisons considered, we found there were 69, 70, and 74 genes being involved in the “ubiquitin mediated proteolysis” pathway, 54 of which were common among all three (Fig. 7, File S7). Among these 54 circRNAs, we considered eight circRNAs that had a significant common dysregulation trend in the three intrinsic comparisons, or they had completely reverse expression trend between SCC and ADC but were equally-expressed (|log2(fold change)| < 1) in ASC (Table 2).

We detected the eight key candidate circRNAs by qRT-PCR, and identified that hsa_circ_0005325 and hsa_circ_0005728 were up-regulated in SCC and ADC compared with the ANT, hsa_circ_0035811 and hsa_circ_0059960 were down-regulated in SCC and ADC compared with the ANT. In addition, we found that the relative expression levels of the four circRNAs, hsa_circ_0000989, hsa_circ_0004503, hsa_circ_0018484, and hsa_circ_0004258, were up-regulated in SCC, but down-regulated in ADC compared to ANT, which is consistent to our sequencing results (Fig. 8, File S5).

Construction of a ceRNA network based on candidate circRNAs and their host genes

We predicted the target miRNAs of the eight key candidate circRNAs (File S8), and built a target pool containing 245 miRNAs using the miRanda databases. A total of 143/245 miRNAs were overlapped with the target miRNAs predicted by TargetScan v7.2. Two interaction networks were constructed to reflect the relationships of ceRNA with Cytoscape v3.6.1 (Figs. 9 & 10). With a particular focus on the 102 miRNAs that showed no relation with the “ubiquitin mediated proteolysis” biological pathway, we generated a visual chordal graph to demonstrate the functional enrichment of miRNA and target genes using the Novomagic online platform tool (https://magic.novogene.com) (Fig. 11). TGF-beta signaling pathway, the most significant pathway, involved 61/102 miRNAs and 66 target genes. The profile data derived from KEGG pathway analysis by DIANA-miRPath v3.0 and with DIANA-microT-CDS algorithm are shown in File S9.