Circular RNA expression profile in peripheral whole blood of lung adenocarcinoma by high: Throughput sequencing.

BACKGROUND: Lung adenocarcinoma (LA) is a most common form of non-small cell lung cancer (NSCLC). To date, there are still no effective early diagnosis methods for patients to be cured in time. Noncoding RNA plays an important role in oncogenesis and tumor development. The expression profile of circular RNA (circRNA) in peripheral whole blood (PWB) of LA has not been systematically investigated. In this study, we identified the differentially expressed (DE) circRNAs in PWB of LA by high-throughput sequencing.
METHODS: Five paired LA and normal participants PWB samples were chosen to investigate the expression profile of circRNAs by high-throughput sequencing. Twenty LA and 10 normal controls PWB samples were subjected to reverse-transcription polymerase chain reaction (RT-PCR) for validation of circRNAs expression profile. Gene Ontology (GO) functional analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and circRNA-miRNA network analysis was also performed to predict the function of circRNAs in PWB.
RESULTS: A total of 10566 circRNAs were identified and annotated, most of the circRNAs were exonic (78.14%). Statistical analysis revealed 4390 DE circRNAs, in which were 3009 upregulated circRNAs and1381downregulated circRNAs in LA. RT-PCR results showed that circRNA expression in LA was higher than that in controls. GO functional analysis, KEGG pathway analysis, and circRNA-miRNA network analysis all showed that circRNAs correlated with tumor development and progression to a certain degree. The current study is the first to systematically characterize and annotate circRNA expression in PWB of LA. Some host genes of the DE circRNAs were involved in tumor signaling pathway and had complicated correlations with tumor related miRNAs, indicating that circRNAs might involve in development and progression of LA.
CONCLUSIONS: Our study revealed that circRNAs were abnormally expressed in PWB of LA, which might offer potential targets for the early diagnosis of the disease and new genetic insights into LA.

Introduction

Lung cancer is one of the most important malignant tumors worldwide, and is the leading cause of death worldwide. The 5-year survival rate of lung cancer remains only 17.4% for NSCLC.[ Although the 10-year survival rate of stage Ia lung cancer could reach approximately 92% with optimum treatment, about 85% of patients with lung cancer are diagnosed at more advanced stages.[ Detection of early stage lung cancer is quite essential to improve the overall survival. Lung cancer is divided into SCLC and NSCLC. NSCLC accounts for 80%, including LA and lung squamous cell carcinoma (LSCC). LA remains the most common subtype of NSCLC, for which the mortality and morbidity have been increasing year by year.[ Therefore, it is important to identify a novel cancer specific biomarker for LA patients to help make early diagnosis and guide clinical treatment.

An increasing evidence has shown that noncoding RNAs, such as long noncoding RNAs (lncRNAs), microRNAs (miRNAs), and circular RNAs (circRNAs), play important roles in the development and progression of many tumors and could be used as therapeutic targets and prognostic factors for Lung cancer.[ CircRNAs are a class of noncoding RNA molecules that lack 5′-3′ ends and a poly A tail, covalently forming closed continuous loops.[ Generally, circRNAs are stable molecules, and some have an effective role of sponge gene regulation ability.[ CircRNAs, with their specific features, have superior potential to serve as a novel biomarker for human diseases.

Some studies have provided evidence that circRNAs are DE in tumors tissue of LA and play an important role in carcinogenesis because of participating in cancer related pathways.[ However, whether circRNAs in PWB of LA are sensitive and specific biomarkers remains largely unknown. In order to reveal the potential roles of circRNAs in LA patients, we performed circRNA expression profiling in PWB from LA patients and healthy controls. We identified a number of circRNAs that are upregulated or downregulated in LA patients. The results suggest that these circRNAs may be developed as novel noninvasive biomarkers for LA patients in the future.

Methods

This study has been cleared by the Ningbo Medical Center Lihuili Hospital Ethics Review Board for human studies. Written informed consent was obtained from all participants included in the study.

Patients

Five LA samples and 5 healthy controls were enrolled from the Ningbo LiHuili Hospital in December 2016. Twenty LA samples and 10 healthy controls for validation were enrolled from the Ningbo LiHuili Hospital from October to December 2018. Patients were in accordance with the following criteria:

a pathologic diagnosis of LA;

no previous cases of cancer;

HIV negative;

no receiving any preoperative treatment (chemotherapy and /or radiotherapy);

no other important organ system diseases.

Details are shown in Tables 1 and 2. Specimen were immediately frozen quickly in liquid nitrogen after extraction and then stored at −80°C until RNA extraction.

Table 1

Pathological characteristics of study population for sequencing (N = 10).

Table 2

Pathological characteristics of study population for validation (N = 30).

RNA extraction

Trizol (Invitrogen, USA) was used for extraction of total RNA from LA and control PWB according to manufacturer's instructions. The RNA concentration and purity were checked by OD A260/A280 (≥1.8) and A260/A230 (≥1.6), and the yield and quality were assessed using an Agilent2100 Bioanalyzer (Agilent Technologies, USA) and Ribo-Zero H/M/R Kit (Illumina, MRZH11124). The RNA integrity number (RIN) of extracted RNA was >7.0

Next-generation RNA sequencing

The cleaved RNA fragments were reverse-transcribed to create the cDNA, which were next used to synthesize U-labeled second-stranded DNAs with Escherichia coli DNA polymerase I, RNase H and dUTP. An A-base is then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters are ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme treatment of the U-labeled second-stranded DNAs, the ligated products are amplified with PCR by the following conditions: initial denaturation at 95°C for 3 minutes; 8 cycles of denaturation at 98°C for 15 seconds, annealing at 60°C for 15 seconds, and extension at 72°C for 30 seconds; and then final extension at 72°C for 5 minutes. The average insert size for the final cDNA library was 300 bp ± 50 bp. At last, we performed the 150 bp paired-end sequencing on an Illumina Hiseq 4000 (LC-Bio Hangzhou, China) following the vendor's recommended protocol.

Bioinformatics and data analysis

First, cutadapt was used to remove the reads that contained adaptor contamination, low quality bases and undetermined bases. Then sequence quality was verified using FastQC (LC-Bio, Hangzhou, China). Each base was measured with a corresponding mass value, which was used to measure the accuracy of sequencing. Q20 and Q30 indicated the percentage of base whose mass value was greater than or equal to 20 or 30. Q20 or30 was ≥95%. We used Bowtie2 and Tophat2 to map reads to the genome of species. Remaining reads (unmapped reads) were still mapped to genome using tophat-fusion. CIRCExplorer was used to denovo assemble the mapped reads to circularRNAs. Then, back splicing reads were identified in unmapped reads by tophat-fusion and CIRCExplorer. CircRNAs with fold change ≥2 and P < .05 were considered to be statistically significant.

MiRNA - binding sites on circRNAs and the assumed target genes of these miRNAs are all predicted by custom - written software based on Target-scan and Miranda software (LC-Bio, Hangzhou, China). Predicted tumor related miRNAs were picked. KEGG analysis was carried out for the differential expression circRNA associated target genes. CircRNAs including the miRNA - binding sites and the appropriate miRNAs were subjected to analysis with Cytoscape software (LC-Bio, Hangzhou, China) to construct miRNA-circRNA networks and display interactions.

Quantitative real-time reverse transcription PCR

Total RNAs were extracted from PWB samples with Trizol reagent (Invitrogen, Carlsbad, CA). Reverse transcription was performed using the Invitrogen Superscript cDNA Synthesis kit (Invitrogen, Carlsbad, CA). CircRNA expression was measured through RT-PCR using the SYBR Green PCR kit on an Analytikjena thermocycler (Analytikjena, QTOWER 2.2, GER). The PCR primers used in this study were as follows: circRNA5430, TCATTTCCCAACAGATTACCC (forward) and GCTTGCCAATGGAACACT (reverse); circRNA6783, TGTCCTGCAATTAGGTATCCGAAT (forward) and CTCTGGTTTATTTGTGGGGAAAGC (reverse). Samples were run in triplicate for analysis. Relative circRNA expression was calculated with the 2−ΔΔCt method.

Statistical analysis

The fold change in circRNA expression was calculated by comparing expression levels between cancers and controls in PWB. x -tests were applied to flag where the proportion of positive results showed a significant difference between the LA patients and healthy controls. Mann-Whitney U-test or student t test was used to evaluate the significance of the difference between the 2 groups. We used the filter criteria of fold change ≥2 and P < .05 to screen for DE circRNAs. Agilent feature extraction software (version 11.0.1.1, Agilent, Santa Clara, CA) was used to analyze the acquired images. R software (LC-Bio, Hangzhou, China) was used to execute quantile normalization and for GO and KEGG analysis.

Results

CircRNA expression profile in PWB

We first analyzed the expression profile of circRNAs in PWB of 5 LA and 5 controls specimens by high-throughput sequencing. A total of 10,566 circRNAs were identified from PWB of 5 LA and 5 controls and annotation was performed. DE circRNAs with statistical significance between the 2 groups were displayed through fold change and P value (fold change ≥2 and P < .05). 4390 circRNAs were identified to significantly express differentially between the 2 groups. The 3009 circRNAs were significantly upregulated, and 1831 circRNAs were remarkably downregulated more than 2-fold in LA samples group compared with controls group on volcano plots and differential expression histogram (Fig. 1). Volcano plots and histogram showed that the circRNA expression levels were clearly distinguished and clustered between PWB samples from LA and controls specimens.

Figure 1

A Histogram to identify differentially expressed circRNAs. Histogram was used to identify differentially expressed circRNAs in LA PWB vs normal PWB. The x-axis represents circRNAs between 2 groups (cancer and control group), while the y-axis represents the number of 2 groups. Red histogram represents upregulated circRNAs and green histogram represents downregulated circRNAs. B scatter plots and volcano plots to identify differentially expressed circRNAs. Scatter plots were used to identify differentially-expressed circRNAs in LA PWB vs normal PWB. The x-axis represents fold-change values (log2 scaled), while the y-axis represents P values (-log10 scaled). Red scaled dots represent differentially expressed circRNAs in 2 groups.

The top 20 upregulated and top 20 downregulated circRNAs for LA are listed in Table 3.

Table 3

Top 20 ranking significantly upregulated and downregulated circRNAs in lung adenocarcinoma (fold change ≧ 2, P ≦ .05).

According to the source, circRNAs were divided into 3 categories, including exonic circRNAs, intronic circRNAs, and intergenic circRNAs (Fig. 2).

Figure 2

The figure shows the 3 types of differentially expressed circRNAs, most of which originate from exonic circRNAs.

GO enrichment and KEGG analysis

GO enrichment analysis revealed the top 10 significantly enriched GO terms associated with the DE circRNAs, indicating that the DE genes might be most related to “protein binding and poly(A) RNA binding (Fig. 3C)”,“cellular response to cystrol (Fig. 3B)”, and “Biological process- viral process” (Fig. 3A).

Figure 3

CircRNA of GO function analysis. A, GO function analysis to identify the biological process of circRNA. B, GO function analysis to identify the cellular component of circRNA. C, GO function analysis to identify the molecular function of circRNA. y-axis, function of closely related to circRNAs; x-axis, enrichment score –log 10 (P value). circRNA = circular RNA, GO = gene ontology.

In order to confirm the pathways in which the DE circRNAs were involved, we analyze the genes that produced them by KEGG pathway analysis. KEGG enrichment demonstrated that the DE circRNAs associated with the process of ubiquitin mediated proteolysis (Fig. 4). A total of 43 genes were identified in this pathway. These results showed that DE genes might be related to tumor signaling pathways.

Figure 4

Kyoto encyclopedia of genes and genomes pathway analysis to identify the enriched circRNA. y-axis, signaling pathways closely related to circRNAs (TOP 20); x-axis, enrichment score –log10 (P value). circRNA = circular RNA.

Prediction of circRNA-miRNA interaction and network visualization

According to the magnitude of fold changes and P value of the DE circRNAs from 5 LA samples and the known functions of circRNAs related to tumor process. The top 20 upregulated and top 20 downregulated circRNAs predicted miRNA response elements (MREs) in Table 4.

Table 4

Top 20 upregulated and downregulated circRNAs predicted miRNA response elements (MREs).

In addition, we showed 5 upregulated and upregulated circRNAs listed for analyzing the interaction network between circRNAs and miRNAs. The analysis showed that all of 10 circRNAs contained their respective MREs. A total of 50 miRNAs regulated by the 10 circRNAs were displayed as a network generated by cytoscape software (Fig. 5).

Figure 5

CircRNA-miRNA network. Cytoscape was used to generate a circRNA-miRNA coexpression network. The network map consists of the previously identified 5 significantly upregulated circRNAs (represented by red nodes) and 5 significantly downregulated circRNAs (represented by green nodes) along with their 50 target miRNAs (represented by Green V line node). circRNA = circularRNA, miRNA = microRNA.

DE circRNA validation

CircRNA5430 and circRNA6783 were selected for validation because they had large expression difference (greater change of Log2fold and smaller P value). The result showed that circRNA5430 and circRNA6783 in 20 LA PWB were expressed higher than those in controls, accord with the high-throughput data (Fig. 6A and B).

Figure 6

A. DE circRNA validation. The relative levels of circRNA5430 by RT-PCR analysis in peripheral whole blood samples from the 20 lung adenocarcinoma (CA) compared with 10 normal controls (NC). The levels of circRNA5430 from RT-PCR were expressed as a ratio of average 2−ΔΔCt of circRNA5430 form CA to NC. Data are expressed as means ± SD (P < .05). circRNA, circular RNA; RT-PCR, reverse-transcription polymerase chain reaction. B. The relative levels of circRNA6783 by RT-PCR analysis in peripheral whole blood samples from the 20 lung adenocarcinoma (CA) compared with 10 normal controls (NC). The levels of circRNA6783 from RT-PCR were expressed as a ratio of average 2−ΔΔCt of circRNA6783form CA to NC. Data are expressed as means ± SD (P < .05). circRNA = circular RNA, RT-PCR = reverse-transcription polymerase chain reaction.

Discussion

CircRNAs are a novel class of extensive and stable endogenous RNAs that regulate gene expression in organisms.[ The covalently closed loop structures make circRNAs more stable than liner RNA and insensitive to RNA exonuclease or RNase.[ These features make circRNAs the potential ideal biomarkers for human diseases. Recent studies have showed that circRNAs can be used as diagnostic or predicted biomarkers for colon cancer,[ hepatocellular,[ gastric cancer,[ leukemia,[ and lung cancer.[ However, little is known about the role of circRNAs in LA of PWB. In this study, we performed a high-throughput sequencing of dysregulated circRNAs by comparing the expression profiles of PWB from patients with LA and those from healthy controls. These observations may assist future pathophysiology researches in LA and help determine if circRNAs in PWB could be used as novel, noninvasive biomarkers for LA diagnosis and treatment. At present, no studies on the role of circRNAs in PWB of LA have been reported; therefore, investigating the expression profile of circRNA in PWB of LA and the corresponding functional mechanism is particularly important.

In this study, we first analyzed the profiling of circRNAs in PWB of 5 LA and 5 controls by high-throughput sequencing. The 4390 circRNAs were identified to significantly express differentially between the 2 groups. DE circRNAs (3009 upregulated and 1381 downregulated) revealed an important role of circRNAs in LA.

Furthermore, we identified the major significantly changed GO terms, most correlated pathways and predicted circRNA- miRNA interactions with bioinformatics analysis.

GO and KEGG pathway analysis were used to confirm the functional parental genes of DE circRNAs. We found that these parental genes were functionally predicted to be related to cellular components, molecular function regulation, and related to biological process. KEGG pathway analysis showed that many DE circRNAs corresponded to tumor signaling pathway, such as Ubiquitin mediated proteolysis, Long-term potentiation, MAPK signaling pathway et al. In particular, the MAPK signaling pathway was involved in LA development and progression.[ We hypothesized that some circRNAs might change activation of the tumor signaling pathway to influence LA development and progression. Therefore, circRNAs in PWB could be used as LA biomarkers for early diagnosis and treatment.

The circRNA -miRNA network showed that each DE circRNAs related to the 5 tumor relevant miRNAs, affirming circRNA may involve in the development of tumor. Recent studies have shown that circRNAs negatively regulate miRNAs and dramatically expand the endogenous network in competition with endogenous RNA.[ CircRNAs act as miRNA sponges in LA, and their potential biological functions need further study.

The study explored circRNA expression profiles of LA by high-throughput sequencing and identified DE circRNA.

Our study has some limitations. First, the sample size was too small to make accurate conclusions. We would increase the number of specimens in the future investigations. Second, DE circRNAs discovered in this study need further validation. Third, the study on the role of circRNAs in human tumor was still in its infancy. In the future, a larger study should be performed to verify our findings and confirm whether these findings can serve as novel biomarkers for LA diagnosis treatment.

In conclusion, our study is the first to measure circRNA expression in PWB from LA patients and healthy controls. The findings may improve our understanding of the role of circRNAs in PWB of LA patients, which suggest that circRNAs might serve as novel biomarkers that may have promising functions and valuable clinical significance in LA.

Acknowledgments

We owe thanks to the patients and their family. We thank the staffs at Ningbo Medical Center Lihuili Hospital.

Author contributions

Conceptualization: Yinyu Mu, Qiaoping Wu.

Data curation: Chunnian Wang.

Formal analysis: Yinyu Mu, Yun Fei Huang.

Investigation: Dongdong Yang, Chunnian Wang.

Methodology: Fuyi Xie, Chunnian Wang.

Project administration: Yinyu Mu, Fuyi Xie.

Resources: Dongdong Yang, Guodong Xu.

Software: Guodong Xu.

Supervision: Guodong Xu.

Validation: Yun Fei Huang.

Writing – original draft: Yinyu Mu.

Writing – review & editing: Yinyu Mu.

YINYU MU orcid: 0000-0003-1239-2761.