Analysis of exosomal circRNAs upon irradiation in pancreatic cancer cell repopulation

Background Pancreatic cancer is one of the most malignant tumors. However, radiotherapy can lead to tumor recurrence, which is caused by the residual surviving cells repopulation stimulated by some molecular released from dying cells. Exosomes may mediate cell-cell communication and transfer kinds of signals from the dying cells to the surviving cells for stimulating tumor repopulation. Circular RNAs (circRNAs) may be one vital kind of exosomal cargos involving in modulating cancer cell repopulation. Methods Next generation sequencing (NGS) and bioinformatics were performed to analyze and annotate the expression and function of exosome-derived circRNAs in pancreatic cancer cells after radiation. Four circRNAs were chosen for qRT-PCR analysis to validate the sequencing results. Results In this study, 3580 circRNAs were annotated in literatures and circBase among 12,572 identified circRNAs. There were 196 filtered differentially expressed circRNAs (the up-regulation and down-regulation respectively is 182 and 14, fold change > 2, p-value < 0.05). Regulation of metabolic process and lysine degradation were the main enriched biological processes and pathway according to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Conclusions The hsa_circ_0002130-hsa_miR_4482-3p-NBN interaction network suggested potential sponging miRNA and target mRNA. Our results provided potential functions of circRNAs to explore molecular mechanisms and therapeutic targets in pancreatic cancer cell repopulation upon irradiation.


Background
Pancreatic cancer is a highly fatal disease, which morbidity is almost the same as the mortality. Focusing on the global burden of cancer worldwide, the number of new cases and deaths of pancreatic cancer respectively was 458,918 and 432,242, taking up 2.5 and 4.5% of all the analyzed cancers in 185 countries [1]. The high mortality of pancreatic cancer might partly ascribe to the lack of effective strategies. Radiotherapy should be an optimal choice for those regional unresectable disease without detectable metastasis in most cancers. Although American Society for Radiation Oncology (ASTRO) recently issued a guideline for treatment of pancreatic cancer patients using radiotherapy, either conventional radiotherapy or SBRT, etc., the quality of evidence for recommendation using of radiotherapy was largely low [2]. Efforts to improve the efficacy of radiotherapy in pancreatic cancer is ongoing. New radiation dose, modality, fraction size, and sequencing, as well as in combination with other treatment strategies such as immunotherapy have been largely studied [3]. However, no significant progress has been achieved so far.
As one of the 5Rs of radiobiology, tumor accelerated repopulation after radiation was known as the main cause of treatment failure [4,5]. Researches by our team have revealed that cancer cell repopulation after radiotherapy were related to the dying cells caused by the treatment, which could inactivate or activate relevant signaling to stimulate the residual surviving cells to fast grow. The related molecules, such as PGE2, microRNAs and other molecules, played as a switch to promote the repopulation [6][7][8]. Exosomes, one type extracellular vesicle of about 30-150 nm, contain rich cargos such as DNA, mRNA and some non-coding RNAs [9]. Recently, it has been well reported that exosomes might transmit the key factors to the receptor cells, which promote tumorigenesis and progression [10,11]. Thus, further exploring the contents of pancreatic cancer cell-derived exosomes would help to reveal the molecular mechanism of repopulation upon irradiation.
CircRNAs are predominantly in the cytoplasm, which formed a covalently closed continuous loop without 5'cap or 3' polyadenylated tail [12]. Therefore, they are resistant towards exonucleases. Based on their biogenesis mechanism, different kinds of circRNAs could be characterized to mainly five parts, exonic circRNAs, intronic circRNAs, antisense circRNAs, sense overlapping cir-cRNAs and intergenic circRNAs [13]. Among them, the exonic circRNAs were the most detected [14]. As endogenous competitive RNAs, circRNAs can compete for miRNAs through MREs (miRNA recognition elements, MREs) [15,16]. Recently, circRNAs could be sorted into exosomes and participate in cancer progression [17]. Especially, exosome-derived circRNAs could promote invasive growth through miRNA sponge in pancreatic cancer [18]. Although sponging miRNAs has been partly demonstrated during radiation, the function of circRNAs remains largely unknown.
During last decades, next generation sequencing (NGS) has been widely used to identify the differentially expressed genes, annotate functional pathways at genomic, transcriptional and epigenetic level, which greatly assists clinicians in early diagnosis and screening of highrisk populations, and eventually allows the development of personalized therapy in pancreatic cancer [19,20]. In this study, we generated RNA sequencing data from 4 types of pancreatic cancer cells (PANC-1, SW1990, BxPC-3, MIA PaCa-2), which were treated with or without irradiation, and identified~12,570 circRNAs. The accuracy of sequencing was verified by quantitative real-time RT-PCR (qRT-PCR) of the differentially expressed circRNAs (DE-circRNAs). Bioinformatics analysis including GO and KEGG were then performed to annotate the selected DE-circRNA functions. A circRNA-miRNA-mRNA network was subsequently constructed to reveal the molecular regulatory networks. In brief, the whole work flow of our study was shown in Fig. 1.

Cell culture and treatment
Four pancreatic cancer cells were used for analysis, including PANC-1 (ATCC® CRL-1469), BxPC-3 (ATCC® CRL-1687), SW1990 (ATCC® CRL-2172) and MIA PaCa-2 (ATCC® CRL-1420). These cells were divided into two groups, group A as the control group and group B as the experimental group. Ten Gy of X-ray radiation was used in experimental group as we reported to induce tumor repopulation before [6,7]. Fatal bovine serum (FBS) centrifuged at 120,000 g for 18 h was used for cell culture to avoid the interference with experimental results [21].

Exosome isolation and identification
Differential centrifugation was used for exosome purification. Briefly, centrifugation was respectively performed at 300, 2000 and 10,000 g to remove the alive cells, dead cells and cell debris. The crude extract of exosomes was obtained by ultra-high speed centrifugation (> 100,000 g), and the operation was repeated twice to remove the contaminated protein for collecting the purified exosomes [22].

RNA extraction and purification
Total RNAs from two groups (group A and B) were extracted using TRIzol reagent (Life Technologies, USA), and then treated with DNase I (Takara, China) to reduce interference of genomic DNA, according to manufacturer's instructions. The quantity and purity of total RNAs was determined by NanoDrop ND-100 (Thermo Fisher Scientific, USA). Well qualified RNA was further studied in the following experiments.
circRNA library construction and sequencing rRNAs were removed by Ribo-Zero rRNA Removal Kits (Illumina, USA) from total RNAs. TruSeq Stranded Total RNA Library Prep Kit (Illumina, USA) was used for constructing RNA libraries following the manufacturer's instructions. The quality and quantification of libraries was detected by BioAnalyzer 2100 system (Agilent Technologies, USA). Ten pM libraries were denatured to singlestranded DNA molecules, which amplified in situ as clusters and sequenced for 150 cycles on Illumina HiSeq Sequencer according to the manufacturer's instructions.

Sequence mapping and circRNA annotation
Paired-end reads were obtained from Illumina HiSeq 4000 sequencer for quality control by Q30. After 3' adaptortrimming, Cutadapt software (v1.9.3) was used to remove low quality reads, and STAR software (v2.5.1b) was used to compare high quality reads with hg19 human reference genome. CircRNAs were detected and identified using DCC software (v0.4.4). EdgeR software (v3.16.5) was used for data normalization and DE-circRNA analysis. GO and KEGG analysis was performed to predict the function of associated genes of candidate circRNAs.

Identification of DE-circRNAs
DE-circRNAs in the two groups were identified using the edgeR software with quasi-likelihood F test (fold change > 2 and p-value < 0.05). The DE-circRNAs were log 2 transformed for standardization and further visualization. Visual graphics were generated using ggplot2 in R (https:// cran.r-project.org/web/packages/ggplot2/index.html).

Functional annotation of target miRNAs and prediction of interaction networks
In terms of functional annotation of candidate circRNAs, the source genes of DE-circRNAs were clustered in GO annotations (http://www.geneontology.org) and KEGG pathways (http://www.genome.jp/kegg/) by DAVID Bioinformatics Resources (http://david.ncifcrf.gov/home.jsp). P-value was calculated by fisher test and adjusted by Ben-jamini& Hochberg. The result of prediction of potential circRNA-miRNA binding sites was obtained through Tar-getScan (http://www.targetscan.org) and miRanda (http:// www.microrna.org/microrna/home.do). The mRNA-miRNA-circRNA interaction analysis was performed for top 3 expression quantity by using TargetScan, miRanda, Fig. 1 A brief workflow. In the chart, our work was mainly divided into three steps. Through these steps, potential circRNAs were harvested. CircRNA profiling indicated that dying cell derived-exosome carried circRNAs that could regulate modulators and influence the tumor repopulation in pancreatic cancer Circinteractome (https://circinteractome.nia.nih.gov/) and Circbank (http://www.circbank.cn/). The enrichment score was comprehensively considered to the score from TargetScan and the thermodynamic properties of the binding site from miRanda. The predicted target genes of the DE-circRNAs were further subjected to GO term and KEGG pathway analyses. The gene network analysis was performed using Cytoscape.

qRT-PCR
Total RNA was reverse transcribed using random primers with the PrimeScript RT reagent kit (RR037A, Takara), following the manufacturer's protocol. We randomly selected 4 DE-circRNAs, including 3 up-regulated and 1 downregulated circRNAs. Details of the primer sequences are summarized in Supplementary Table 4. Only primers whose PCR products showed a single peak in the melting curve were considered for qRT-PCR validation. Then, qRT-PCR was performed using SYBR Premix Ex Taq II (Takara, Dalian, China) on a QuantStudio 6 Flex (Life technologies, USA) according to the manufacturer's instructions. The expression level of candidate circRNAs was normalized to βactin and then calculated using the 2 -△△Ct method. All experiments were repeated three times and presented as the means ±SD.

RNase R treatment
DNase-treated total RNAs were incubated for 1 h at 37°C with or without RNase R (Epicentre). We used linear RNA FBXW7, GAPDH, β-actin, 18S rRNA (which are poly(A)-tailed and must be degraded by RNase R treatment) and circRNA FBXW7 as controls.

Patient-derived xenograft mouse model
All animal experiments were conducted according to protocols approved by Shanghai General Hospital's Institutional Animal Care and Use Committee (IACUC). Clinical surgical specimens from tumor patients were cut into pieces and subcutaneously implanted into pentobarbital-anaesthetized BALB/c nude mice (about 6-week-old; purchased from Shanghai SLAC Laboratory Animal Co.,Ltd) at right groin to build patient-derived xenograft (PDX) mouse model. When PDX tumor grew into about 500 mm 3 , mice were subjected to 10 Gy radiation or mock treatment (5 mice for each group and sample size was chosen based on routine of cancer investigations). One week after radiation, the mice were anaesthetized with pentobarbital, blood was collected through heart punctures and EDTA was used as the anticoagulation agent. Mice were then euthanized by cervical dislocation. Plasma was separated by centrifuge at 3000 rpm for 15 min. Exosomes were collected from plasma by ultracentrifugation.

Western blot
Protein was isolated by 10% SDS-polyacrylamide gel electrophoresis (SDS-page) and then transferred to a nitrocellulose membrane (Sigma-Aldrich, USA). After blocking with 5% defatted milk for 90 min, membranes were incubated with TSG101 antibody (1:1000; Abcam, UK) at 4°C overnight. Membranes were washed and incubated for 1.5 h with secondary antibodies (1:5000; Goat Anti-Rabbit IgG (H + L) Dylight 800, Odyssey, LI-COR, USA). The band intensity was measured using Empiria Studio. GAPDH was used as a loading control.

Statistical analysis
GraphPad Prism 7.0 was used for statistical calculation (version 7.0c; GraphPad Software, San Diego, USA). All the data were expressed as ± SD of the mean value for three independent measurements. t test was used to evaluate the differences between groups. p-value < 0.05 were considered statistically significant.

Exosome characterization and RNA-seq library construction
Exosomes secreted by 4 kinds of pancreatic cancer cells (PANC-1, SW1990, MIA Paca-2, BxPC-3) that were treated with (group B) or without (group A) irradiation were obtained by ultracentrifugation. Under transmission electron microscope (TEM), exosomes were irregular spheres ranging of 30-100 nm in diameter (Fig. 2a, b), which agreed well with the references [9][10][11]. The concentration of exosomes in irradiation group was higher than those in control one (Fig. 2c, d). Western blot was performed to analyze TSG101, a common biomarker of exosome [22]. Exosomes showed high expression of TSG101 (Fig. 2e), however, there was no expression inside the cells.
To systematically identify circRNAs from exosomes, 4 samples were sequenced using quality-controlled total RNA-seq (rRNA depleted), which generated approximately 380 million reads. A flowchart of sequencing process was shown in Fig. 2f Table 3). Whereafter, the data were normalized and analyzed to identify differentially expressed circRNA.
Expression profiling of circRNAs 12,572 cicRNAs were quantitated in whole samples, and 8992 of them were novel based on comparison with cir-cBase and published literatures. The majority of cir-cRNAs were exonic circRNAs in all samples, whereas the sense overlapping circRNAs just took up a small proportion. Differently, among differently expressed circRNAs (DE-circRNAs, fold change > 2, p-value < 0.05), the intergenic circRNAs took up the least (Fig. 3a). The details of the novel and known circRNAs in 5 isoforms were shown in Fig. 3b. The proportion rate among 5 isoforms of DE-circRNAs was further shown in Supplementary Fig. 1A. Among the exonic circRNAs, the mean length of known circRNAs were majorly concentrated on 200-300 nt (254 nt, Fig. 3c). Known circRNAs were distributed across all chromosomes, while DE-circRNAs were distributed across all chromosomes but chromosome 22 and chromosome Y. Significantly, the circRNAs across chromosome 1 and 2 were more enriched than other chromosomes (Fig. 3d). The dysregulated distribution of DE-circRNAs was also exhibited (Suppementary Fig. 1). Furthermore, the details of expression profiles of all cir-cRNAs quantitated in our study were presented in Supplementary Data 1 and 2. The data of chromosomal location of the DE-circRNAs were shown in Supplementary Data 3.  (Fig. 4a, fold change > 2, p-value < 0.05). Especially, all the down-regulated circRNAs were located on the chromosome M. We noted that 174 cir-cRNAs were 6 folds significantly up-regulated, 7 cir-cRNAs were 8 folds significantly up-regulated, and 1 circRNA (hsa_circ_0000284) were more than 8 folds significantly up-regulated. Moreover, 12 circRNAs were 6 folds significantly down-regulated, and 2 novel circRNAs (chrM:14131-15754-and chrM: 14131-115754+) were 8 folds significantly down-regulated (Fig. 4b). The volcano plot also revealed the dysregulated circRNAs clustered in red squares (Fig. 4c).
Five hundred eighty-six GO terms were enriched by analyzing the host genes of the up-regulated circRNAs, which contained biological process (BP), cellular component (CC)  (Fig. 4d). It was revealed that genes relating to the regulation of metabolic process, binding and intracellular were most enriched in GO analysis. Twenty six enriched KEGG pathways were associated with lysine degradation, phagosome, axon guidance, allograft rejection, HTLV-I infection, endocytosis, Type I diabetes mellitus and Graft-versus-host disease (Fig. 4e). Furthermore, the details of all terms of GO and KEGG were shown in Supplementary Data 4 and 5.

Experimental validation of circRNA-sequencing
Four circRNAs (hsa_circ_0000419, hsa_circ_0001523, hsa_circ_0000825, and chrM:14131-15754-) identified as DE-circRNAs by high-throughput sequencing were randomly selected for validation by qRT-PCR. The results were consistent with the sequencing results (Fig. 5a). We also selected the above mentioned circRNA, circFBXW7, and paired linear RNA FBXW7 [23] mRNA of common housekeeping genes (β-actin, GAPDH, and 18S rRNA) as controls for RNase R resistance experiment. The results revealed that the expression of all linear RNAs (FBXW7, β-actin, GAPDH and 18S rRNA) decreased after 1 h of RNase R digestion (** p-value < 0.001). This indicated that circRNAs were resistant to RNase R digestion, whereas linear RNAs were sensitive to R treatment (Fig. 5b).

Prediction of circRNA-miRNA interaction network
CircRNAs could act as "miRNA sponge" to regulate miRNA level [17,18]. To explore the function of cir-cRNAs, the binding sites between DE-circRNAs and miRNAs were conducted by online databases. Thus, we selected the top 5 DE-circRNAs (p-value < 0.0001) to build a circRNA-miRNA interaction network, whose sequence might paired with the seed sequence of miRNA (Fig. 6). The circRNA-miRNA interaction analysis for all DE-circRNAs was shown in Supplementary Data 6.

hsa_circ_0002130/hsa_miR_4482-3p/NBN interaction axis
According to the highest binding force of miRNA to mRNA, 3 top known DE-circRNAs (hsa_circ_0002130, hsa_circ_0000825 and hsa_circ_0005882, p-value < 0.001) were chosen to construct the mRNA-miRNA-circRNA interaction network axis (Fig. 7a), and the newly circRNAs were shown in Supplementary Fig. 2. For example, we displayed the hsa_circ_0002130 to further observe the influence of target genes. In previous studies, the irradiation-induced cellular damage included DNA strand breaks [24]. Nibrin (NBN) is the gene whose product is associated with DNA double-strand break repair and DNA damage-induced check point activation [25]. Therefore, we supposed NBN could be related with the repopulation in pancreatic cancer. Here, NBN might be regulated by hsa_miR_4482-3p, a hsa_ circ_0002130 binding miRNA, through two possible binding sites (Fig. 7b). To further validate the expression changes of the axis, we performed qPCR assay. The expression level of hsa_circ_0002130 in cancer cell-derived exosomes was upregulated after radiation ( Fig. 7c (Left)). RNase R experiments were shown in the Supplementary  Fig. 1C. To further validate the change of hsa_circ_ 0002130 in vivo, we performed qPCR assay in exosomes derived from plasma of the irradiated/unirradiated PDX tumor, which was generated in our previous work [26]. We identified that exosomal hsa_circ_0002130 was apparently higher in plasma from the irradiated mice compared with that in the untreated group (Fig. 7c (Right)). Furthermore, mRNA level of NBN was up-regulated in exosomes during irradiation (Fig. 7d). According to TCGA, the survival curve for NBN was shown in Fig. 7e (Log rank p-value = 0.0234), which showed the survival probability of patients with the higher level expression of NBN was shorter than the lower ones.

Discussion
In our study, 12,572 circRNAs were identified, and 3580 circRNAs were annotated as the known circRNAs. We obtained the reads of circRNAs from 5 parts depending on the location of the chromosome, including exonic, Results were represented as means ± standard deviation (SD). *p-value < 0.05. b CircRNAs were resistant to RNase R digestion, whereas the linear RNAs were sensitive to RNase R digestion. Four circRNAs were examined, which were qualified in expression level. And FBXW7, GAPDH, β-actin, and 18S rRNA were used as negative controls and known circRNA FBXW7 was used as a positive control. The results of quantitative Real-time PCR (qRT-PCR) were evaluated by 2 -△CT method. Results were represented as means ± standard deviation (SD) *p-value < 0.05, ** p-value < 0.001 intronic, sense overlapping, intergenic and antisense, which were similar to the publication [13]. Different cir-cRNAs were produced by selecting different splicing sites and existed mainly in exons. However, we also found most of the identified circRNAs (71.52%) were novel, which mainly concentrated on four parts except the exon (Fig. 3b). Besides, the mean length of these identified exonic circRNAs was focus on 200-300 nt (Fig. 3c). Interestingly, it was reported that circRNAs releasing could be influenced by their own length, which indicated that the length of circRNAs might affect its formation [27]. Moreover, the unique structure of cir-cRNA without a poly(A)-tail rendered themselves highly insensitive to ribonuclease. We showed that circRNAs were more stable than linear RNAs through RNase R digestion. Alternatively, Ribo-zero could also be used to circRNA validation [28].
The results in BP indicated that the highly expressed target genes might be involved in regulation of metabolic process, RNA biosynthetic and metabolic process, and nucleic acid-templated transcription. These GO terms might be associated with nucleotide metabolism, which was involved in radiotherapy [29]. Refer to MF, "Binding", "protein methyltransferase activity", "histone methyltransferase activity", "histone methyltransferase activity (H3-K4 specific)" and "histone-lysine N-methyltransferase activity" were highly enriched. As reported, VDR and p53 in irradiated HEK 293 T was related to various biological effects, Fig. 6 circRNA-miRNA network prediction. The top 15 significantly expressed circRNAs and top 5 predicted miRNAs were selected to generate a network map. The circRNA-miRNA network was built using bioinformatics online programs (TargetScan and miRBase). The red circle indicated the up-regulated circRNAs and the green ones indicated the down-regulated ones which mediated gene transcription [30]. Furthermore, aberrant methylation could lead to dysregulation of related transcription-factor genes, which might be closely bound up to irradiation [31,32]. KEGG analysis indicated that the lysine degradation was important, which was consistent with the results of GO analysis. Lysine degradation mainly participated in the control of gene expression and transcription. These findings indicated the cargos in exosomes could regulate different pathway and biological processes concerning cell proliferation and irradiation. However, further experimental verification is still needed, especially onto tumor repopulation owing to radiotherapy.
Refer to miRNA sponge, numerous miRNA binding sites were predicted by Circinteractome, Circbank, Miranda and TargetScan, which revealed circRNAs could interact with miRNAs. Based on bioinformatics, mRNA-miRNA-circRNA regulatory network was constructed to clarify the interactions between molecules. As shown (Fig. 7a), the hsa_circ_0002130/hsa_miR_4482-3p/NBN interaction axis was involved in the regulatory network. NBN gene is a member of the MRE11/RAD50 doublestrand break repair complex. Several studies confirmed that NBN was crucially involved in DNA double strand fracture repairing and DNA damage-induced checkpoint Fig. 7 Hsa_circ_0002130-hsa_miR_4482-3p-NBN interaction axis. a The circRNA-miRNA-mRNA interaction network was predicted using online bioinformatics programs (Circinteractome, Circbank, TargetScan, and miRBase). The red circle indicated the up-regulated circRNAs, the green indicated the down-regulated ones, the arrow and the hexagon indicated miRNAs and target genes respectively. b The predicted binding sites between hsa_circ_0002130 and hsa_miR_4482-3p were exhibited. c The expression of hsa_circ_0002130 in PANC-1 cells (Left) and PDX tumors (Right) were examined respectively by qRT-PCR. d The expression of NBN was determined by qRT-PCR in samples of sequencing, PANC-1 and SW1990. *p-value < 0.05 and data indicated the mean ± SD. e The survival data from TCGA (Log rank p-value = 0.0234) activation [33]. The expression level of hsa_circ_0002130 was also up-regulated during irradiation, consistent with the expression pattern of NBN. Given the potential binding sites among hsa_circ_0002130/hsa_miR_4482-3p/NBN, hsa_circ_0002130 was considered to promote the expression of NBN through the absorption of hsa-miR-4482-3p. As we recently found that irradiated dying tumor cells-derived exosomes promoted DNA damage repair and the survival of damaged tumor repopulating cells [26] the exosomal hsa_circ_0002130 might also mediate the vital role of the exosomes. Thus, further works could be directed to identify the role of dying tumor cell-derived exosomal hsa_circ_0002130 in potentiating tumor repopulation via promoting DNA damage repair, which was deduced to be mediated by the hsa_circ_0002130/hsa_ miR_4482-3p/NBN axis. Remarkably, some circRNAs were also shown to be translated into proteins [34,35]. Thus, the underlying mechanisms still needed further investigations.
Nevertheless, there are still several drawbacks inherent to our study. The small sample size was for RNAsequencing, and samples from the patients were hard for us to collect. Therefore, the influence from the patients should be further validated. On the other hand, even if we believe that they worked, we did not perform further experiments to confirm the regulatory relationship among circRNAs, miRNAs and mRNAs, as well as other functions of circRNAs.

Conclusions
Hereinabove, we first identified and annotated 12,572 circRNAs in exosomes of human pancreatic cancer cells upon irradiation using RNA-seq analysis. Among them, the DE-circRNAs were associated with methylation, which was in accord with the signaling pathways. In further prediction, hsa_circ_0002130 could target NBN by combing with the hsa_miR_4482-3p and high expression of NBN revealed a worse survival rate.